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 ?