Vtgyjfy

I am looking for a Matlab code for the $2$ stage multistep method :
$$y_{n+2} - y_n = hleft[(1/3)f_{n+2} + (4/3)f_{n+1} + (1/3)f_nright]$$
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

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)=%15fn", Y(end,1), Y(end,2));

printf("exact values y1(1)=%.15f, y2(1)=%15fn", 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