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

Question

I am looking for a Matlab code for 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]$$ Note that $f_n = f(x_n,y_n)$. Now, since $f_{n+2} = f(x_{n+2},y_{n+2})$ and $y_{n+2}$ is also the output of the method, cleary this is an implicit method. That means that Newton-Raphson shall be used for the calculation of the final output. The point of that code is to use it in a paper-work I am assembling for some conclusions over methods. Theoritically, I have proven that this is a method of order $4$, meaning that : $|T_n| \leq ch^4$. I want to prove this experimentally by using a certain method, but this can't be done if I do not have the code for the given method.

Since I am not that good in programming, I am kindly requesting some help regarding the code. I would really appreciate any elaborated algorithm. Thanks in advance !

My attempt at a code using the rk4 function already defined for the first 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

Answer 1

Using fsolve to short-cut solving the implicit equation without explicitly coding a Newton method gives the following implementation

function [T, Y] = multistep(FunFcn,t0,tfinal,step,y0)
  % initialize arrays
  N = 1;
  T = [ t0 ];
  Y = [ y0 ];
  F = [ FunFcn(t0, y0) ];

  % perform an RK3 step to have initial error O(h^4)
  % to get the requisite two samples to start the multistep method,
  % or use RK4 to stay with the O(h^5) local error
  N = N+1;
  T(N) = T(N-1)+step;
  Y(N,:) = RK3step(FunFcn,T(N-1),step,Y(N-1,:));
  F(N,:) = FunFcn(T(N), Y(N,:));

  % set error tolerances for implicit solver, 
  % the implicit error (times step size) should 
  % contribute less than the discretization error
  options = optimset("TolX", 1e-3*step^4, "TolFun", 1e-6*step^4);

  % main loop
  while abs(T(N)+0.6*step - tfinal) > 0.501*step 
    N = N+1;
    T(N) = T(N-1)+step ;
    % perform the multi-step method
    % y2 -y0 = h/3*(f0+4*f1+f2)
    % y2 - h/3*f(t2,y2) = b = y0 + h/3*(f0+4*f1)
    b = Y(N-2,:) + step/3*(F(N-2,:)+4*F(N-1,:));
    % solve the implicit equation using the built-in fsolve
    % Coding your own solver would allow a better management 
    % of the Jacobian approximation
    correct = @(y) ( y-step/3*FunFcn(T(N), y) - b );
    % predict = Y(N-2,:)+2*step*F(N-1,:);
    % use a 4th order predictor
    predict = 5*Y(N-2,:)-4*Y(N-1,:)+2*step*(2*F(N-1,:)+F(N-2,:));
    Y(N,:) = fsolve(correct, predict, options );
    F(N,:) = FunFcn(T(N), Y(N,:));
  end

  % Now tfinal is between T(N) + 0.1*step and T(N)+1.1*step
  % compute the final value with an RK3 step to add no more than an O(h^4) error
  N = N+1;
  T(N) = tfinal;
  Y(N,:) = RK3step(FunFcn,T(N-1),T(N)-T(N-1), Y(N-1,:));
end    

Using RK3 to handle the first and last steps is the minimum required, it is probably better to use RK4 to have the error contributions negligible against the global error order 4.

The Runge-Kutta step functions implement the well-known methods:

function y1=RK3step(FunFcn, t0, step, y0)
  k1 = step*FunFcn(t0,y0);
  k2 = step*FunFcn(t0+0.5*step, y0+0.5*k1);
  k3 = step*FunFcn(t0+    step, y0+2*k2-k1);
  y1 = y0 + (k1+4*k2+k3)/6;
end

function y1=RK4step(FunFcn, t0, step, y0)
  k1 = step*FunFcn(t0         , y0       );
  k2 = step*FunFcn(t0+0.5*step, y0+0.5*k1);
  k3 = step*FunFcn(t0+0.5*step, y0+0.5*k2);
  k4 = step*FunFcn(t0+    step, y0+    k3);
  y1 = y0 + (k1+2*k2+2*k3+k4)/6;
end

Now test the error order on a problem with known exact solution

%% test problem y''+sin(y) = r = p''+sin(p), y(0)=p(0), y'(0)=p'(0)
%% p=sin(t), p(0)=0, p'(0)=1, r(t)=-sin(t)+sin(sin(t))

TestODE = @(t,u) [ u(2), -sin(u(1))-sin(t)+sin(sin(t)) ];

t0 = 0; tf = 2;
[ T, Y ] = multistep(TestODE, t0,tf,0.2, [ 0, 1 ]);

printf("values found y1(1)=%.15f, y2(1)=%15f\n", Y(end,1), Y(end,2));
printf("exact values y1(1)=%.15f, y2(1)=%15f\n", sin(tf)  , cos(tf)  );
figure(1);
plot(Y(:,1), Y(:,2));

And to confirm the order, logarithmically plot the error against the step size

H = 0.1.^[1:0.5:3];
E1=[]; E2=[];
for h = H
  [ T, Y ] = multistep(@(t,u)TestODE(t,u),t0,tf,h,[ 0, 1 ]);
  E1(end+1) = max(abs( Y(:,1)'-sin(T) ));
  E2(end+1) = max(abs( Y(:,2)'-cos(T) ));
end
figure(2)
loglog(H,E1,H,E2)
end