Bdf3 Method no convergence

Question

I try to implement the Backward differentiation formula 3 (BDF3) on Octave and the plot shows no convergence. Here's my attemp. For the first 4 steps I use the rk4 method and since bdf3 is an indirect method I use Newton-Raphson to solve the non-lineal problem $y_{n+1}=G(y_{n})$ Here is the BDF3 Method: $y_{n+3}=\frac{18}{11}y_{n+2}-\frac{9}{11}y_{n+1}-\frac{2}{11}y_{n}+\frac{6}{11}f(x_{n+3},y_{n+3})$ I already know that the implementation of rk4 is right

function [tout, yout] = rk4(FunFcn,t0,tfinal,step,y0)
% Costant Stepsize 4 steop order 4 rk4.m
% Initialization
ceta = [1/2  1/2  1]';
alpha  = [ [ 1/2      0      0     0    ]
          [   0      1/2      0    0    ]
          [   0       0       1   0    ]]';
beta =  [1/6  1/3 1/3 1/6 ]';
stages=4;
t = t0; y = y0(:);f = y*zeros(1,stages);tout = t;yout = y.';
% The main loop
while abs(t- tfinal)> 1e-6 
   if t + step > tfinal, step = tfinal - t; end
  % Compute the slopes     
    temp = feval(FunFcn,t,y);
    f(:,1) = temp(:);
    for j = 1:stages-1
      temp = feval(FunFcn, t+ceta(j)*step, y+step*f*alpha(:,j));
      f(:,j+1) = temp(:);
    end
    t = t + step;
    y = y + step*f*beta(:,1);
    tout = [tout; t];
    yout = [yout; y.'];
end;

Here is the BDF3 Method

function [tout, yout] = bdf3(FunFcn, t0, tfinal, step, y0)
% calls rk4 for 4 steps
 tolbdf3 = 1e-6;
 tolnr = 1e-9;
 maxiter = 50;
 diffdelta = 1e-6;
 stages = 4;
 [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).';
 while abs(t - tfinal)> tolbdf3
 if t + step > tfinal, step = tfinal - t; end
 t = t + step;
 yp0 = y;
 ypf = yp0;
 yp = inf;
 iter = 0;
 while (abs(yp - ypf)>= tolnr) && (iter < maxiter)
 df = 1/diffdelta * (feval(FunFcn, t, yp0+diffdelta) - feval(FunFcn, t, yp0));
 yp = yp0 - 1/(6/11*step*df - 1) * (18/11*yout(end) -9/11*yout(end-1) -2/11*yout(end-2)+ 6/11*step*feval(FunFcn, t, yp0) - yp0);
 ypf = yp0;
 yp0 = yp;
 iter = iter + 1;
 end
 y = yp;
 tout = [tout; t];
 yout = [yout; y.'];
 end
end

Here are the function and the real function

function yout=gefunc(t,y)
yout=2*t-y;
end

function [ytrue]=getrue(t)
    ytrue=exp(-t)+2*t-2;
    end

And the run script

t0=0;
tfinal=5;
h=0.04;
y0=-1;
[tout3, yout3] = bdf3('gefunc',t0,tfinal,h,y0);
plot(tout3,yout3,'m',tout3,getrue(tout3),'g');

The result is

Answer 1

Your BDF coefficients are wrong. The $y$-coefficients on the right side have to add to 1, or if you put all $y$ terms on the left side, the sum of coefficients has to be zero.

Correct to $$ y_{n+3}=\frac{18}{11}y_{n+2}-\frac{9}{11}y_{n+1}\color{blue}{+\frac{2}{11}y_{n}}+\frac{6}{11}f(x_{n+3},y_{n+3})\Delta x $$

See

• How do you derive the backward differentiation formula of 3rd order using interpolating polynomials?

• wikipedia

This will correct the artificial growing oscillations connected to the roots $0.900 \pm 0.550i=1.05e^{0.175\pi i}$ of the wrong characteristic polynomial to the correct ones $0.318\pm 0.284i= 0.426e^{0.232\pi i}$ of a decaying oscillation.

On the method

The coefficients for BDF methods are easy to derive. The differentiation operator $D$ is to be expressed in the translation/shift/propagation operator $e^{-hD}$ and its powers modulo $O(h^{p})$. Now $$ D=-\frac1h\ln(1-(1-e^{-hD})=\frac1h\left((1-e^{-hD})+\frac{(1-e^{-hD})^2}2+\frac{(1-e^{-hD})^3}3+...\right) $$ Truncation using $1-e^{-hD}=O(h)$ and expansion gives $$ D=\frac1h\left((1-e^{-hD})+\frac{1-2e^{-hD}+e^{-2hD}}2+\frac{1-3e^{-hD}+3e^{-2hD}-e^{-3hD}}3\right)+O(h^3) \\ =\frac1{6h}(11-18e^{-hD}+9e^{-2hD}-2e^{-3hD})+O(h^3) $$