Wrong results for $2$ stage multistep method $y_{n+2} - y_n = h\left[(1/3)f_{n+2} + (4/3)f_{n+1} + (1/3)f_n\right]$

Question

I need to fix a code to utilise the $2$ stage multistep method :

$$y_{n+2} - y_n = h\left[(1/3)f_{n+2} + (4/3)f_{n+1} + (1/3)f_n\right]$$

Since this is an implicit method, I used a Newton-Raphson approach for the final determination of $y_{n+2}$.

Below, is my attempt at a code implementing the already given rk4 method (runge-kutta 4) for the first approximation steps :

function [tout, yout] = newcorrect(FunFcn,t0,tfinal,step,y0)
maxiter = 1000;
tolnr   = 1e-9;
diffdelta = 1e-6;
stages=2;
[tout,yout]=rk4(FunFcn,t0,t0+(stages-1)*step,step,y0);
tout=tout(1:stages);
yout=yout(1:stages);
t = tout(stages);
y = yout(stages).';
% The main loop
 while abs(t- tfinal)> 1e-6 
    if t + step > tfinal, step = tfinal - t; end
    t = t + step;
    yn0 = y;
    ynf = yn0;
    yn = inf;
    iter = 0;
    while (abs(yn - ynf)>= tolnr) && (iter < maxiter)
        df = 1/diffdelta * (feval(FunFcn,t, yn0+diffdelta) - feval(FunFcn, t, yn0));
        yn = yn0 - 1/(1/3*step*df - 1) * (4/3*step*feval(FunFcn,tout(end),yout(end)) + 1/3*step*feval(FunFcn,tout(end-1),yout(end-1)) + 1/3*step*feval(FunFcn, t, yn0) -yn0 + yout(end-1));
        ynf = yn0;
        yn0 = yn;
        iter = iter + 1;
    end
    y = yn;
    tout = [tout; t];
    yout = [yout; y.'];
 end
end

This is to be used for showing experimentally that when you divide the step by two, the fraction of the maximum absolute errors per consecutive different steps $h$ is approximately equal to $2^{-p}$ where $p$ is the order of the method. The method is proven to be of order $4$. The exercise whishes the initial starting point to be $k0=2$ thus I created the following script :

clear all;
t0 = 1;
tfinal = 3;
y1 = 2;
tout  =t0:0.01:tfinal;
k0 = 2;
kf = input('enter final k:')
for k = k0:kf
  h(k-1) = 2^(-k);
  [tout,yout] = newcorrect('f0', t0, tfinal, h(k-1), y1);
  outputs{k-1} = [tout,yout];
end
for i = 1:(kf-1)
   maxabserror(i) = max(abs(outputs{1,i}(:,2)-f0true(outputs{1,i}(:,1))));
end
 for i = 1:(kf-2)
   consmax(i) = maxabserror(i+1)/maxabserror(i);
end

The function f0 and f0true are :

function yprime = f0(t,y)
yprime = (t.^2 + y.^2)/(2.*t.*y);
end

function y = f0true(t);
y = sqrt(t.*(t+3));
end

When I run the method, though, the problem is that despite providing very good approximations for the solution values, the experimental conclusion of the order cannot be carried out (the fraction vary and do not converge to $2^{-p}$.

**After running the script for $kf = 10$, I get the consecutive absolute max errors fraction vector :

$$\textbf{consmax} = \big[0.0704, \; 0.0657, \; 0.0639, \; 0.0632, \; 0.0628, \; 0.0627, \; 0.0656, \; 0.2400\big] $$

As it is easy to see, every term of it tends to $1/16 = 2^{-4}$ which indeed is what we want as $4$ is the order of the method. BUT what is it going on with that last value ?

Where is my mistake ? Is this something that can be explained ?

Answer 1

The local error in each integration step is composed of 3 parts:

• the discretization error of the method
• the floating point error of actually performing the method in precision-limited data types, and
• specifically for implicit methods, the error in solving the implicit equation.

The first error is of size $O(h^{p+1})$, the second of size $\sim \mu$, $\mu$ the machine constant of the number type. The last error depends on the implementation.

Usually one uses a predictor-corrector scheme with an explicit method step as predictor, often of the same order $p$ so that ensuing corrector steps towards the solution of the implicit equation only influence the coefficient of the leading power in $h$ of the error. Here one easy non-trivial predictor is the midpoint step $y_{n+1}^{[0]}=y_{n-1}+2hf(t_n,y_n)$ with error $O(h^3)$. .

Every Newton-like corrector implementation (like the simplified Newton method) will reduce that error by a factor of $\approx Ch^2$ so that in general after the first corrector iteration the order of the method is reached. If $Ch^2\ll1$, the update of the corrector step may be taken as proxy for the remaining distance to the exact solution, so instead of a fixed number of 1 to 3 corrector iterations one can also stop the corrector loop if the update is smaller than for example $0.1\,h^4\|y_{n+1}^{[0]}-y_n\|$ (it is $O(h^5)$ and scales with the scaling of the problem).

Obviously, setting this last tolerance independent of $h$ will introduce a lower bound for the local and thus global error. You used the previous value as predictor, $y_{n+1}^{[0]}=y_n$. This has error $O(h)$. One Newton step is always performed, so that the first error that is really compared to tolnr = 1e-9; has magnitude $O(h^3)$. A second step down to error $O(h^5)$ and thus the order of the method is not performed if $h^3\lesssim 10^{-9}\implies h\lesssim 10^{-3}\sim 2^{-10}$, which is about what you observed.