MATLAB TUTORIAL for the First Course. Part III: Predictor-corrector 3

Prof. Vladimir A. Dobrushkin

This tutorial contains many matlab scripts.
You, as the user, are free to use all codes for your needs, and have the right to distribute this tutorial and refer to this tutorial as long as this tutorial is accredited appropriately. Any comments and/or contributions for this tutorial are welcome; you can send your remarks to

Email Vladimir Dobrushkin

Predictor-corrector 3

Consider the initial value problem \( y' = f(x,y) , \quad y(x_0 ) = y_0 . \) The third order multistep method is based on the following Adams--Bashforth-Moulton formulas:

\[ \begin{split} p_{n+1} &= y_{n} + \frac{h}{12} \left[ 23\, f_n -16\, f_{n-1} + 5\, f_{n-2} \right] , \\ y_{n+1} &= y_n + \frac{h}{12}\left[ 5\,f\left( x_{n+1} , p_{n+1} \right) + 8\,f\left( x_{n} , y_{n} \right) - f\left( x_{n-1} , y_{n-1} \right) \right] , \end{split} \]
where as usual \( f_n = f \left( x_{n} , y_{n} \right) . \) Recall that we work with a uniform mesh of points \( x_n = x_0 + n\,h , \quad n=0,1,2,\ldots ; \) with fixed step size h.

There is known another third order multistep method, called Hermite predictor-Milne corrector method:

\[ \begin{split} p_{n+1} &= 5\, y_{n-1} - 4\,y_{n} + h \left[ 4\, f_n -2\, f_{n-1} \right] , \\ y_{n+1} &= y_{n-1} + \frac{h}{3}\left[ f\left( x_{n+1} , p_{n+1} \right) + 4\,f\left( x_{n} , y_{n} \right) + f\left( x_{n-1} , y_{n-1} \right) \right] , \end{split} \]

Predictor-corrector of order 2:

ABM2[f_, {x_, x0_, xn_}, {y_, y0_, y1_}, steps_] :=
Block[{xold = x0, yold = y0, ynew = y1, sollist = {{x0, y0}}, h},
h = N[(xn - x0)/steps]; sollist = Append[sollist, {x0 + h, y1}];
Do[xnew = xold + h;
xnew2 = xold + h*2;
k1 = f /. {x -> xold, y -> yold};
k2 = f /. {x -> xnew, y -> ynew};
abnew = ynew + (3*h/2)*k2 - h*k1/2;
k3 = f /. {x -> xnew2, y -> abnew};
bm = ynew + h*k2/2 + h*k3/2;
sollist = Append[sollist, {xnew2, bm}];
xold = xnew;
yold = ynew; ynew = bm, {steps}];
Return[sollist[[steps + 1]]]]

To apply this method, use the following commands:

f[x_, y_] := x^2 + y^2;
ABM2[f[x, y], {x, 0, 1.0}, {y, 1, Exp[0.1*0.1]}, 10]

If you want to see the whole table of values:

ABM2[f_, {x_, x0_, xn_}, {y_, y0_, y1_}, steps_] :=
Block[{xold = x0, yold = y0, ynew = y1, sollist = {{x0, y0}}, h},
h = N[(xn - x0)/steps]; sollist = Append[sollist, {x0 + h, y1}];
Do[xnew = xold + h;
xnew2 = xold + h*2;
k1 = f /. {x -> xold, y -> yold};
k2 = f /. {x -> xnew, y -> ynew};
abnew = ynew + (3*h/2)*k2 - h*k1/2;
k3 = f /. {x -> xnew2, y -> abnew};
bm = ynew + h*k2/2 + h*k3/2;
sollist = Append[sollist, {xnew2, bm}];
xold = xnew;
yold = ynew; ynew = bm, {steps}];
Return[sollist]]

Example.

f[x_, y_] := 2*x*y;
ABM2[f[x, y], {x, 0, 1.0}, {y, 1, Exp[0.1*0.1]}, 10]
Out[3]= {1., 2.73005} {{0, 1}, x, y, {0.1, 1.01005}, {0.2, 1.04096}, {0.3, 1.09458}, {0.4, 1.1743}, {0.5, 1.2854}, {0.6, 1.43554}, {0.7, 1.63575}, {0.8, 1.9017}, {0.9, 2.25576}, {1., 2.73005}, {1.1, 3.37112}}

Predictor-corrector method of order 2:

function  Complete


Complete