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$
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$ – user2799609 Jul 23 '18 at 15:04
  • $\begingroup$ What you need to print and why this part of the code For[y = 0, y < 30, y++, ? $\endgroup$ – Alex Trounev Jul 23 '18 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$ – user2799609 Jul 23 '18 at 17:47
  • $\begingroup$ The ode system runs for 120 times steps or {t,0,120}? $\endgroup$ – Alex Trounev Jul 23 '18 at 18:42
  • $\begingroup$ Sorry, {t,0,120} $\endgroup$ – user2799609 Jul 23 '18 at 19:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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