MATLAB TUTORIAL for the First Course. Part III: Milne--Simpson Method

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

Milne--Simpson Method

Another popular predictor-corrector scheme is known as the Milne or Milne--Simpson method. See
Milne, W. E., Numerical Solutions of Differential Equations, Wiley, New York, 1953.
Its predictor is based on integration of the slope function f(t, y(t)) over the interval \( \left[ x_{n-3} , x_{n+1} \right] \) and then applying the Simpson rule:

\[ y(x_{n+1}) = y(x_{n-3}) + \int_{x_{n-3}}^{x_{n+1}} f(t, y(t)) \,{\text d}t . \]
The predictor uses the Lagrange polynomial approximation for f(t, y(t)) based on four mesh points \( (x_{n-3} , f_{n-3} ), \ (x_{n-2} , f_{n-2} ), \ (x_{n-1} , f_{n-1} ), \ (x_{n} , f_{n} ). \) It is integrated over the interval \( \left[ x_{n-3} , x_{n+1} \right] . \) This produces the Milne predictor:
\[ p_{n+1} = y_{n-3} + \frac{4h}{3} \left( 2\,f_{n-2} - f_{n-1} + 2\,f_n \right) , \qquad n=3,4,\ldots . \]
The above formula was discovered by William Edmund Milne in 1953. This formula has local truncation error of order O(h5). The Milne corrector is developed similarly.
\[ y_{n+1} = y_{n-1} + \frac{h}{3} \left( f_{n-1} +4\, f_{n} + f_{n+1} \right) , \qquad n=3,4,\ldots ; \]
where \( f_{n+1} = f \left( x_{n+1} , p_{n+1} \right) \) is based on the predictor value pn+1. The Milne--Simpson multistep method is sometimes unstable. The Milne method converges when
\[ h < \frac{3}{\left\vert f_y (x_n , y_n ) \right\vert} , \qquad n=0,1,2,\ldots . \]
The error terms for the numerical integration formulas used to obtain both the predictor and corrector are of the order O(h5). The local truncation errors for both, predictor and corector are
\[ \begin{split} \phi (x_{n+1}) - p_{n+1} &= \frac{28}{90} \, \phi^{(5)} (\xi_{n+1})\,h^5 , \\ \phi (x_{n+1}) - y_{n+1} &= \frac{-1}{90} \, \phi^{(5)} (\eta_{n+1})\,h^5 , \end{split} \]
where \( y = \phi (x) \) is teh actual solution to the initial value problem \( y' = f(x,y) , \quad y(x_0 ) = y_0 \) and slope function f(x,y) is sufficiently smooth function on the interval of interest.

Example. Let us start with the Riccati equation \( y' = x^2 + y^2 \) subject to the initial condition \( y(0) =-1 . \) Its solution is expressed through Bessel functions:

\[ d(x)\, y(x) = x \left[ J_{-3/4} \left( \frac{x^2}{2} \right) \left( \Gamma^2 \left( \frac{3}{4} \right) + \pi \right) - Y_{-3/4} \left( \frac{x^2}{2} \right) \Gamma^2 \left( \frac{3}{4} \right) \right] . \]
where the denominator
\[ d(x) = Y_{1/4} \left( \frac{x^2}{2} \right) \Gamma^2 \left( \frac{3}{4} \right) - \left[ J_{1/4} \left( \frac{x^2}{2} \right) \left( \Gamma^2 \left( \frac{3}{4} \right) +\pi \right) \right] \]
has the first positive null at x = 2.223378383 as the figure shows
d[x_] = BesselY[1/4, x^2 /2]*Gamma[3/4]^2 - BesselJ[1/4, x^2 /2]*(Gamma[3/4]^2 + Pi)
Plot[d[x], {x, 0, 5.5}, PlotStyle -> Thick]
The sollution blows up near the point 2.223378.

 

function  M=milne(f,t,y)
% Input:    f is the function entered as a string 'f'
%           t is the vector of abscissas
%           y is the vector of ordinates
% Remark:   The first four coordinates of t and y must have 
%           starting values obtained with RK4
% Output:   M=[t',y'] where t is the vector of abscissas
%           y is the vector of ordinates
n=length(t);
if n < 5, break,end;
F=zeros(1,4); 
F=feval(f,t(1:4),y(1:4));
h=t(2)-t(1);
pold=0;
yold=0;
for k=4:n-1
    % Predictor
    pnew=y(k-3)+(4*h/3)*(F(2:4)*[2 -1 2]'); 
    % Modifier
    pmod=pnew+28*(yold-pold)/29;
    t(k+1)=t(1)+h*k;
    F=[F(2) F(3) F(4) feval(f,t(k+1),pmod)];
    % Corrector
    y(k+1)=y(k-1)+(h/3)*(F(2:4)*[1 4 1]')); 
    pold=pnew;
    yold=y(k+1);
    F(4)=feval(f,t(k+1),y(k+1));
end
M=[t',y'];

Example. Milne's solution to \( y' = 30-5y, \quad y(0) =1 \) over interval [0,10] with N = 93 steps produces oscillations. It is stabilized when N = 110. 


Complete