# Bdf3 Method no convergence

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

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

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)$$

• Yes, of course. Thank you. – Lutz Lehmann Dec 22 '19 at 14:51
• Could you please tell me what exactly should I correct in my code because I understood what you said but I think that I already have done it. Thanks in advance – kawrik Dec 23 '19 at 0:27
• You need to change the sign of 2/11. A small term with big consequences. – Lutz Lehmann Dec 23 '19 at 8:00
• Ohh yes you are right. Thanks!!! – kawrik Dec 23 '19 at 14:28