计算数值方法收敛阶

本文介绍如何计算数值方法的收敛阶。以 Euler 方法、梯形公式、改进 Euler 方法为例，并用 MATLAB 编程，最后展示了输出结果。

误差阶

数值解法的基本思想是，通过某种离散化手段将微分方程转化为差分方程，如单步法

$$u_{n+1}=u_n+h\phi (x_n,u_n,h).$$

它在 $x_n$ 处的数值解为 $u_n$，而初值问题在 $x_n$ 处的精确解为 $u(x_n)$，记 $e_n=u(x_n)-u_n$ 称为整体截断误差。收敛性就是讨论当 $x=x_n$ 固定且 $h=\frac{x_n-x_0}{n}\to 0$ 时 $e_n\to 0$ 的问题。

收敛性定义：若一种数值方法（如单步法对于固定的 $x=x_n+nh$，当 $h\to 0$ 时有 $u_n\to u(x_n)$，其中 $u(x)$ 是初值问题的准确解，则称该方法是收敛的。

收敛性定理：假设单步法具有 $p$ 阶精度，且増量函数 $\phi (x_n,u_n,h)$ 关于 $u$ 满足 Lipschitz 条件

$$|\phi (x,u,h)-\phi (x,\overline{u},h)|⩽L_{\phi }|u-\overline{u}|$$

又设初值 $u_0$ 是准确的，即 $u_0=u(x_0)$，则其整体截断误差

$$u(x_n)-u_n=O(h^p)$$

数值方法误差阶的计算

设数值方法的误差阶

$${\Vert u\left(x_n\right)-u_n\Vert }_{\ast }=C{h}^p+O\left(h^{p+1}\right).$$

省略 $h^p$ 的高阶无穷小，记误差

$$E=C{h}^p.$$

两边都取对数

$$\ln (E)=p\ast \ln (h)+C.$$

$p$ 为 $\ln (E)$ 关于 $\ln (h)$ 的斜率

$$p=\frac{\ln (E_1)-\ln (E_2)}{\ln (h_1)-\ln (h_2)}=\frac{\ln (E_1/E_2)}{\ln (h_1/h_2)}.$$

如果函数区间为 $[a,b]$，等距剖分为 $N$ 等份，$h=\frac{b-a}{N}$，则

$$p=-\frac{\ln (E_1/E_2)}{\ln (N_1/N_2)}.$$

线性 ODE 例子

$$\{\begin{array}{rl}& \frac{\mathrm{d}u}{\mathrm{d}t}=t^2+t-u,0<t⩽1,\\ & u(0)=0.\end{array}$$

方程的真解：$u(x)=-e^{-t}+t^2-t+1$.

Euler 方法

MATLAB 程序

% Euler1_error.m
% Euler method for the ODE model
% u'(x)=x^2+x-u, x in [0,1] 
% Initial condition: u(0)=0 ;
% Exact solution: u(x)=-exp(-x)+x^2-x+1.
clear all; close all;
Nvec=[10 50 100 500 1000];        % Number of partitions
Error=[];
fun=@(x,u) x.^2+x-u;              % RHS
for k=1:length(Nvec)
    N=Nvec(k);
    h=1/N;
    x=0:h:1;                      % interval partition
    u(1)=0;                       % initial value
    for n=1:N
        u(n+1)=u(n)+h.*fun(x(n),u(n));
    end
    ue=-exp(-x)+x.^2-x+1;         % exact solution
    error=max(abs(u-ue));
    Error=[Error,error];
end
plot(log10(Nvec),log10(Error),'ro-','MarkerFaceColor','w','LineWidth',1)
%loglog(Nvec,Error,'ro-','LineWidth',1.5)
hold on
%loglog(Nvec, Nvec.^(-1), '--')
plot(log10(Nvec), log10(Nvec.^(-1)), '--')
grid on
%title('Convergence of Euler 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 = -1 ','fontsize',14)

% computating 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 Euler1_error.png
% print -depsc2 Euler1_error.eps

输出结果

Error =  
    0.0459    0.0090    0.0045    0.0009    0.0004

order =  
    1.0123    1.0035    1.0012    1.0003

梯形公式

MATLAB 程序

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

% add annotation of slope
ax = [0.63 0.58];
ay = [0.70 0.65];
annotation('textarrow',ax,ay,'String','slope = -2 ','fontsize',14)

% computating 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 Trapezoidal_error.png
% print -depsc2 Trapezoidal_error.eps

输出结果

Error =
    1.0e-03 *
        0.3069    0.0123    0.0031    0.0001    0.0000

order =
    2.0006    2.0000    2.0000    2.0000

改进的 Euler 方法

MATLAB 程序

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

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

% computating 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 ModiEuler_error.png
% print -depsc2 ModiEuler_error.eps

输出结果

Error =
    0.0027    0.0001    0.0000    0.0000    0.0000

order =
    2.0218    2.0063    2.0022    2.0006