BDF2 方法

本文用 BDF2 方法（属于线性多步法）求解ODE，需要用到方程的牛顿迭代，然后用 MATLAB 编程，最后输出结果的误差阶为 2。

BDF 方法

考虑常微分方程初值问题

$$u^{\prime}=f(t, u), \quad u(0)=u_{0}.$$

向后微分公式 (Backward differentiation formula, BDF) 是数值求解常微分方程的隐式线性多步方法. 由于该方法是 A 稳定的, 通常用来求刚性微分方程.

BDF1 方法就是隐式 Euler 方法 (或后退 Euler 方法)：

$$u_{n+1} = u_n + h f(t_{n+1}, u_{n+1}),\quad n=0,1, \cdots, N-1.$$

BDF2 方法：

$$u_{n+2} - \tfrac{4}{3} u_{n+1} + \tfrac{1}{3} u_n = \tfrac{2}{3} h f(t_{n+2}, u_{n+2}),\quad n=0,1, \cdots, N-2.$$

即

$$u_{n+2} = \tfrac{4}{3} u_{n+1} - \tfrac{1}{3} u_n + \tfrac{2}{3} h f(t_{n+2}, u_{n+2}),\quad n=0,1, \cdots, N-2.$$

BDF3 方法：

$$u_{n+3} - \tfrac{18}{11} u_{n+2} + \tfrac{9}{11} u_{n+1} - \tfrac{2}{11} u_n = \tfrac6{11} h f(t_{n+3}, u_{n+3}),\quad n=0,1, \cdots, N-3.$$

下面使用 BDF2 方法求解常微分方程初值问题.

线性 ODE 例子

$$\left\{\begin{aligned} &\frac{\mathrm{d}u}{\mathrm{d}t}=t^2+t-u, \quad 0 < t \leqslant 1, \\ &u(0)=0. \end{aligned}\right.$$

方程的真解：.

BDF2 方法求解方程的步骤：

1. 使用梯形公式计算 ,

   $$u_1 = u_0 + \tfrac{h}{2}f(t_0,u_0) + \tfrac{h}{2}f(t_1,u_1).$$

   可知：

   $$u_1 = \tfrac{2}{2+h}u_0 + \tfrac{h}{2+h}(t_0^2+t_0 - u_0) + \tfrac{h}{2+h}(t_1^2+t_1).$$

2. 使用 BDF2 方法计算 , ,

   $$y_{n+2} = \tfrac{4}{3} y_{n+1} - \tfrac{1}{3} y_n + \tfrac{2}{3} h f(t_{n+2}, y_{n+2}),\quad n=0,1, \cdots, N-2.$$

   计算可得

   $$u_{n+2} = \tfrac{4}{3+2h} u_{n+1} - \tfrac{1}{3+2h} u_n + \tfrac{2h}{3+2h} (t_{n+2}^2 + t_{n+2}), \quad n=0,1, \cdots, N-2.$$

计算数值解

% BDF2.m
% BDF method for the ODE model
% u'(t)=t^2+t-u, t in [0,1] 
% Initial condition: u(0)=0 ;
% Exact solution: u(t)=-exp(-t)+t^2-t+1.
clear all; close all;
h=0.1;
t=0:h:1;                       % interval partition
N=length(t)-1;
un=zeros(1,N+1);
un(1)=0;                       % initial value
% Trapezoidal rule to compute un(2)
un(2)=2/(2+h)*un(1)+h/(2+h)*(t(1)^2+t(1)+un(1))+h/(2+h)*(t(2)^2+t(2));
% BDF2 method
for n=1:N-1
    un(n+2)=4/(3+2*h)*un(n+1)-1/(3+2*h)*un(n)+2*h/(3+2*h)*(t(n+2)^2+t(n+2));
end
ue=-exp(-t)+t.^2-t+1;          % exact solution
plot(t,ue,'b-',t,un,'r+','LineWidth',1)
legend('Exact','Numerical','location','northwest')
% title('BDF2 method','fontsize',12)
set(gca,'fontsize',12)
xlabel('t','fontsize',16), ylabel('u','fontsize',16)

% computing error
error=max(abs(un-ue))

% print -dpng -r600 BDF2.png
% print -depsc2 BDF2.eps

计算误差阶

% BDF2_error.m
% BDF method for the ODE model
% u'(t)=t^2+t-u, t in [0,1] 
% Initial condition: u(0)=0 ;
% Exact solution: u(t)=-exp(-t)+t^2-t+1.
clear all; close all;
Nvec=[10 50 100 500 1000];   % Number of partitions
Error=[];
for k=1:length(Nvec)
    N=Nvec(k);
    h=1/N;
    t=0:h:1;                 % interval partition
    un(1)=0;                 % initial value
    % Trapezoidal rule to compute un(2)
    un(2)=2/(2+h)*un(1)+h/(2+h)*(t(1)^2+t(1)+un(1))+h/(2+h)*(t(2)^2+t(2));
    % BDF2 method
    for n=1:N-1
        un(n+2)=4/(3+2*h)*un(n+1)-1/(3+2*h)*un(n)+2*h/(3+2*h)*(t(n+2)^2+t(n+2));
    end
    ue=-exp(-t)+t.^2-t+1;    % exact solution
    error=max(abs(un-ue));
    Error=[Error,error];
end
plot(log10(Nvec),log10(Error),'ro-','MarkerFaceColor','w','LineWidth',1)
hold on
plot(log10(Nvec),log10(Nvec.^(-2)),'--')
grid on
%title('Convergence of BDF2 method','fontsize',12)
set(gca,'fontsize',12)
xlabel('log_{10}N','fontsize',14), ylabel('log_{10}Error','fontsize',14)

% add annotation of slope
ax = [0.57 0.53];
ay = [0.68 0.63];
annotation('textarrow',ax,ay,'String','slope = -2','fontsize',14)

% computing convergence order
for n=1:length(Nvec)-1
    order(n)=-log(Error(n)/Error(n+1))/(log(Nvec(n)/Nvec(n+1)));
end
Error
order

% print -dpng -r600 BDF2_error.png
% print -depsc2 BDF2_error.eps

输出结果

Error =
   1.0e-06 *
    0.3778    0.0006    0.0000    0.0000    0.0000    0.0000

order =
    4.0004    4.0001    3.9990    3.9310    4.3726

非线性例子

$$\left\{\begin{aligned} &\frac{\mathrm{d} u}{\mathrm{d} t}=u-\frac{2t}{u}, \quad 0 < t \leqslant 1,\\ &u(0)=1. \end{aligned} \right.$$

方程的真解：.

BDF2 方法求解方程的步骤：

1. 使用梯形公式计算 ,

   $$u_1 = u_0 + \tfrac{h}{2}f(t_0,u_0) + \tfrac{h}{2}f(t_1,u_1).$$

   这里需要用牛顿迭代求解一个非线性方程.

2. 使用 BDF2 方法计算 , ,

   $$u_{n+2} = \tfrac{4}{3} u_{n+1} - \tfrac{1}{3} u_n + \tfrac{2}{3} h f(t_{n+2}, u_{n+2}),\quad n=0,1, \cdots, N-2.$$

   这里每一步都需要用牛顿迭代求解一个非线性方程.

% BDF2Non.m
% BDF method for the nonlinear ODE model
% u'(t)=u-2t/u, t in [0,1]
% Initial condition: u(0)=1;
% Exact solution: u(t)=sqrt(2*t+1)
clear all; close all;
h=0.1;
t=0:h:1;          % interval partition
N=length(t)-1;
un=zeros(1,N+1);
un(1)=1;          % initial value
NI(1,N)=0;        % Record the number of iterations
% Trapezoidal rule to compute un(2)
X=un(1);
Xprev=0;
while abs(X-Xprev) > abs(X)*1e-14
    Xprev=X;
    X=X-(X-h/2*X+h*t(2)/X-un(1)-h/2*(un(1)...
        -2*t(1)/un(1)))./(1-h/2-h*t(2)/(X^2));
    un(2)=X;
    NI(1)=NI(1)+1;
end
% BDF2 method
for n=1:N-1
    X=un(n+1);
    Xprev=0;
    while abs(X-Xprev) > abs(X)*1e-14
    Xprev=X;
    X=X-(X-2/3*h*X+4/3*h*t(n+2)/X-4/3*un(n+1)...
        +1/3*un(n))./(1-2/3*h-4/3*h*t(n+2)/(X^2));
    NI(n+1)=NI(n+1)+1;
    end
    un(n+2)=X;
end
ue=sqrt(2*t+1);    % exact solution
plot(t,ue,'b-',t,un,'r+','LineWidth',1)
legend('Exact','Numerical','location','northwest')
% title('BDF2 method','fontsize',12)
set(gca,'fontsize',12)
xlabel('t','fontsize',16), ylabel('u','fontsize',16)

% computing error
error=max(abs(un-ue))

% print -dpng -r600 BDF2Non.png
% print -depsc2 BDF2Non.eps

计算误差阶

% BDF2Non_error.m
% BDF method for the nonlinear ODE model
% u'(t)=u-2t/u, t in [0,1]
% Initial condition: u(0)=1;
% Exact solution: u(t)=sqrt(2*t+1)
clear all; close all;
Nvec=[100 500 1000 5000 10000];  % Number of partitions
Error=[];
for k=1:length(Nvec)
    N=Nvec(k);
    h=1/N;
    t=0:h:1;          % interval partition
    un(1)=1;          % initial value
    NI(1,N)=0;        % Record the number of iterations
    % Trapezoidal rule to compute un(2)
    X=un(1);
    Xprev=0;
    while abs(X-Xprev) > abs(X)*1e-12
        Xprev=X;
        X=X-(X-h/2*X+h*t(2)/X-un(1)-h/2*(un(1)...
            -2*t(1)/un(1)))./(1-h/2-h*t(2)/(X^2));
        un(2)=X;
        NI(1)=NI(1)+1;
    end
    % BDF2 method
    for n=1:N-1
        X=un(n+1);
        Xprev=0;
        while abs(X-Xprev) > abs(X)*1e-12
            Xprev=X;
            X=X-(X-2/3*h*X+4/3*h*t(n+2)/X-4/3*un(n+1)...
                +1/3*un(n))./(1-2/3*h-4/3*h*t(n+2)/(X^2));
            NI(n+1)=NI(n+1)+1;
        end
        un(n+2)=X;
    end
    ue=sqrt(2*t+1);    % exact solution
    error=max(abs(un-ue));
    Error=[Error,error];
end
plot(log10(Nvec),log10(Error),'ro-','MarkerFaceColor','w','LineWidth',1)
hold on
plot(log10(Nvec),log10(Nvec.^(-2)),'--')
grid on
% title('Convergence of BDF2 method','fontsize',12)
set(gca,'fontsize',12)
xlabel('log_{10}N','fontsize',14), ylabel('log_{10}Error','fontsize',14)

% add annotation of slope
ax = [0.57 0.53];
ay = [0.68 0.63];
annotation('textarrow',ax,ay,'String','slope = -2','fontsize',14)

% computing convergence order
for n=1:length(Nvec)-1
    order(n)=-log(Error(n)/Error(n+1))/(log(Nvec(n)/Nvec(n+1)));
end
Error
order

% print -dpng -r600 BDF2Non_error.png
% print -depsc2 BDF2Non_error.eps
````

**输出结果**

Error =
1.0e-04 *
0.8167 0.0334 0.0084 0.0003 0.0001

order =
1.9858 1.9958 1.9986 2.0000
```