3
$\begingroup$

I have a program for Adams-Bashforth method and I want to do some tries to see if it's working properly but I don't know how to do it.

My code for AdamsBashforth is :

AdamsBashforth = 
  Function[{ti, tf, f, h, yi}, Module[{l, t, y, k1, k2, k3, k4, fl, i},
    l = {{ti, yi}}; t = ti; y = yi;
    fl = {f[t, y]}; i = 4;
    While[t + h <= ti + 3*h,
      k1 = h*f[t, y];
      k2 = h*f[t + h/2, y + k1/2];
      k3 = h*f[t + h/2, y + k2/2];
      k4 = h*f[t + h, y + k3];
      t = t + h; y = y + (k1 + 2*k2 + 2*k3 + k4)/6;
      AppendTo[l, {t, y}]; AppendTo[fl, f[t, y]];]
     While[t + h <= tf,
      y = 
       y + h*(55/24*fl[[i]] - 59/24*fl[[i - 1]] + 37/24*fl[[i - 2]] - 
           9/24*fl[[i - 3]]);
      t = t + h; i = i + 1;
      AppendTo[l, {t, y}]; AppendTo[fl, f[t, y]]];
    l
    ]];

And I have this DOE:

system of three differential equations

$p=s(t)+i(t)+r(t)$.

Now I want to try my code with $p=92342$, $R0= 1+7/15$, $s(0)=p-5$ and $i(0)=6$ but I don't know how to do it...

$\endgroup$
2
  • 1
    $\begingroup$ Is there a special reason why you want to roll your own DE solver instead of using NDSolve? $\endgroup$ Commented Dec 9, 2020 at 12:44
  • $\begingroup$ @SjoerdSmit yes its for a homework. I need to do my own program. Btw this does Adams with order 4 $\endgroup$
    – John hall
    Commented Dec 9, 2020 at 14:17

1 Answer 1

7
$\begingroup$

There are typos in the code, and after minor correction we have

AdamsBashforth = 
  Function[{ti, tf, f, h, yi}, 
   Module[{l, t, y, k1, k2, k3, k4, fl, i}, l = {{ti, yi}}; t = ti; 
    y = yi;
    fl = {f[t, y]}; i = 4;
    While[t + h <= ti + 3*h, k1 = h*f[t, y];
     k2 = h*f[t + h/2, y + k1/2];
     k3 = h*f[t + h/2, y + k2/2];
     k4 = h*f[t + h, y + k3];
     t = t + h; y = y + (k1 + 2*k2 + 2*k3 + k4)/6;
     l = Join[l, {{t, y}}]; fl = Join[fl, {f[t, y]}];]; 
    While[t + h <= tf, 
     y = y + h*(55/24*fl[[i]] - 59/24*fl[[i - 1]] + 
          37/24*fl[[i - 2]] - 9/24*fl[[i - 3]]);
     t = t + h; i = i + 1;
     l = Join[l, {{t, y}}]; fl = Join[fl, {f[t, y]}]];
    l]];

First test we have compared with NDSolve. Tested code

f[t_, y_] := y (1 - y)/3 + Sin[t]
ab = AdamsBashforth[0, 10, f, 1/10., 0.];

Standard code

eq = y'[t] == y[t] (1 - y[t])/3 + Sin[t];
    ic = y[0] == 0;
sol = NDSolveValue[{eq, ic}, y, {t, 0, 10}]

Visualization

Show[Plot[sol[t], {t, 0, 10}, FrameLabel -> {"t", "y"}, 
  Frame -> True], ListPlot[ab, PlotStyle -> Red]]

Figure 1

Second test with ODE system as it has been required in this post. Code with NDSolve

p = 92342; R0 = 
 p - 5; eq1 = {s'[t] == -s[t] i[t], R0 i'[t] == (R0 s[t] - 1) i[t], 
  R0 r'[t] == i[t]}; ic1 = {i[0] == 6, s[0] == 15 (p - 6)/7, 
  r[0] == p - i[0] - s[0]};

sol1 = NDSolve[{eq1, ic1}, {i, s, r}, {t, 0, 10^-3}]

Code with AdamsBashforth

f1[t_, {s_, i_, r_}] := {-s i, i (s - 1/R0), i/R0};ic0 = {1385040/7, 6, 92342 - 6 - 1385040/7} // N

ab = AdamsBashforth[0, 20 10^-5, f1, 10^-6, ic0];

Visualization of two numerical solutions

{Show[Plot[s[t] /. sol1[[1]], {t, 0, 2 10^-4}, PlotRange -> All], 
  ListPlot[Transpose[{ab[[All, 1]], ab[[All, 2, 1]]}], 
   PlotStyle -> Red]], 
 Show[Plot[i[t] /. sol1[[1]], {t, 0, 2 10^-4}, PlotRange -> All], 
  ListPlot[Transpose[{ab[[All, 1]], ab[[All, 2, 2]]}], 
   PlotStyle -> Red]], 
 Show[Plot[r[t] /. sol1[[1]], {t, 0, 2 10^-4}, PlotRange -> All], 
  ListPlot[Transpose[{ab[[All, 1]], ab[[All, 2, 3]]}], 
   PlotStyle -> Red]]}

Figure 2

$\endgroup$
0

Your Answer

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

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