1
$\begingroup$

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 enter image description here

$\endgroup$
2
$\begingroup$

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

| cite | improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ Yes, of course. Thank you. $\endgroup$ – Lutz Lehmann Dec 22 '19 at 14:51
  • $\begingroup$ 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 $\endgroup$ – kawrik Dec 23 '19 at 0:27
  • $\begingroup$ You need to change the sign of 2/11. A small term with big consequences. $\endgroup$ – Lutz Lehmann Dec 23 '19 at 8:00
  • $\begingroup$ Ohh yes you are right. Thanks!!! $\endgroup$ – kawrik Dec 23 '19 at 14:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.