0
$\begingroup$

I am trying to run an ode system for a set number of time steps, take the size at the final time step for select equations and use these to calculate initial conditions for the next set of time steps. I am able to do this fine using a For loop, but now I would like to also vary one or more parameters. I have tried using both Table and Map to no avail. Any suggestions would be much appreciated. Here is a shortened version of my code:

lst = Table[X0 = 600;  
  Y0 = 300;  \[Mu]x = 0.006; \[Mu]y = 0.006; \[Mu]2 = 
   0.00001; \[Lambda] = 1; \[Sigma]x = 1;  \[Sigma]y = 1;  \[Gamma]x =
   0.1; \[Gamma]y = 0.2; \[Epsilon]y = 10; \[Phi] = 10; \[Alpha] = 
   0.008;
  gx[t_, \[Phi]_, \[Lambda]_, \[Epsilon]x_] := 
   If[t < \[Epsilon]x, 0, 
    PDF[GammaDistribution[\[Phi], 1/\[Lambda]], t - \[Epsilon]x]];
  gy[t_, \[Phi]_, \[Lambda]_, \[Epsilon]y_] := 
   If[t < \[Epsilon]y, 0, 
    PDF[GammaDistribution[\[Phi], 1/\[Lambda]], t - \[Epsilon]y]];
  soln[\[Epsilon]x_] =
   sol := NDSolve[{
      \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(x1[t]\)\) == 
       x3[t] gx[t, \[Phi], \[Lambda], \[Epsilon]x] - 
        x1[t] (\[Gamma]x + \[Mu]x),
      \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(x2[
         t]\)\) == \[Gamma]x x1[t] - \[Mu]2 x2[t],
      \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(y1[t]\)\) == 
       y3[t] gy[t, \[Phi], \[Lambda], \[Epsilon]y] - 
        y1[t] (\[Gamma]y + \[Mu]y),
      \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(y2[
         t]\)\) == \[Gamma]y y1[t] - \[Mu]2 y2[t],
      \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(x3[t]\)\) == 0,
      \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(y3[t]\)\) == 0,
      x1[0] == 0, x2[0] == 0, y1[0] == 0, y2[0] == 0, x3[0] == X0, 
      y3[0] == Y0}, 
     {x1, x2, x3, y1, y2, y3}, {t, 0, 120}];
  Xtotal = {X0}; Ytotal = {Y0};
  For[y = 0, y < 30, y++,
   loop = sol;
   X0 = \[Sigma]x y2[120] /. loop[[1, 4]];
   Y0 = \[Sigma]y x2[120]/(1 + \[Alpha] x2[120]) /. loop[[1, 3]];
   Xtotal = Append[Xtotal, X0];
   Ytotal = Append[Ytotal, Y0];
   ];, {\[Epsilon]x, 0, 10, 5}]
ListPlot[Xtotal]
ListPlot[Ytotal]
$\endgroup$

1 Answer 1

0
$\begingroup$

Taking into account our discussion, I propose the second variant depending on the parameters:

    tp = 120; \[Mu]x = 0.006; \[Mu]y = 0.006; \[Mu]2 = 
    0.00001; \[Lambda] = 1; \[Sigma]x = 1;  \[Sigma]y = 1;  \[Gamma]x =
    0.1; \[Gamma]y = 0.2; \[Epsilon]y = 10; \[Phi] = 10; \[Alpha] = 
    0.008;
  gx[t_, \[Phi]_, \[Lambda]_, \[Epsilon]x_] := 
     If[t < \[Epsilon]x, 0, 
       PDF[GammaDistribution[\[Phi], 1/\[Lambda]], t - \[Epsilon]x]];
  gy[t_, \[Phi]_, \[Lambda]_, \[Epsilon]y_] := 
     If[t < \[Epsilon]y, 0, 
       PDF[GammaDistribution[\[Phi], 1/\[Lambda]], t - \[Epsilon]y]];
     sol = ParametricNDSolve[{
          \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(x1[t]\)\) == 
            x3[t] gx[t, \[Phi], \[Lambda], \[Epsilon]x] - 
              x1[t] (\[Gamma]x + \[Mu]x),
          \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(x2[
       t]\)\) == \[Gamma]x x1[t] - \[Mu]2 x2[t],
          \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(y1[t]\)\) == 
            y3[t] gy[t, \[Phi], \[Lambda], \[Epsilon]y] - 
              y1[t] (\[Gamma]y + \[Mu]y),
          \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(y2[
       t]\)\) == \[Gamma]y y1[t] - \[Mu]2 y2[t],
          \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(x3[t]\)\) == 0,
          \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(y3[t]\)\) == 0,
          x1[0] == 0, x2[0] == 0, y1[0] == 0, y2[0] == 0, x3[0] == X0, 
          y3[0] == Y0}, 
        {x1, x2, x3, y1, y2, y3}, {t, 0, tp}, {\[Epsilon]x, X0, Y0}];
Table[{X0 = 600;  
     Y0 = 300;   Xtotal [\[Epsilon]x] = {X0}; 
   Ytotal [\[Epsilon]x] = {Y0};
     For[y = 0, y < 30, y++,
          {X0 = \[Sigma]x y2[\[Epsilon]x, X0, Y0][tp] /. sol;
        Y0 = \[Sigma]y x2[\[Epsilon]x, X0, Y0][
          tp]/(1 + \[Alpha] x2[\[Epsilon]x, X0, Y0][tp]) /. sol;
        Xtotal [\[Epsilon]x] = Append[Xtotal[\[Epsilon]x], X0];
        Ytotal [\[Epsilon]x] = Append[Ytotal[\[Epsilon]x], Y0];}
       ];}, {\[Epsilon]x, 5, 10, 5}];
Table[{ListPlot[Xtotal[\[Epsilon]x], PlotRange -> All, 
   PlotStyle -> Red, PlotLabel -> Row[{"\[Epsilon]x=", \[Epsilon]x}]],
   ListPlot[Ytotal[\[Epsilon]x], PlotRange -> All, 
   PlotStyle -> Orange]}, {\[Epsilon]x, 5, 10, 5}]

fig1

$\endgroup$
6
  • $\begingroup$ My apologies for the horrid description. What I would like is a set of plots for each value of [Epsilon]x $\endgroup$ Jul 23, 2018 at 15:04
  • $\begingroup$ What you need to print and why this part of the code For[y = 0, y < 30, y++, ? $\endgroup$ Jul 23, 2018 at 16:43
  • $\begingroup$ The ode system runs for 120 times steps, the For loop is there to repeatedly take the values of x2 and y2 after each 120 time steps. I just want the output of Xtotal and Ytotal (what the system is outputting after each 120 time steps) $\endgroup$ Jul 23, 2018 at 17:47
  • $\begingroup$ The ode system runs for 120 times steps or {t,0,120}? $\endgroup$ Jul 23, 2018 at 18:42
  • $\begingroup$ Sorry, {t,0,120} $\endgroup$ Jul 23, 2018 at 19:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.