计算数值方法收敛阶

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

误差阶

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

$$u_{n+1}=u_{n}+h \varphi(x_{n}, u_{n}, h).$$

它在 处的数值解为 ，而初值问题在 处的精确解为 ，记 称为整体截断误差。收敛性就是讨论当 固定且 时 的问题。

收敛性定义：若一种数值方法（如单步法对于固定的 ，当 时有 ，其中 是初值问题的准确解，则称该方法是收敛的。

收敛性定理：假设单步法具有 阶精度，且増量函数 关于 满足 Lipschitz 条件

$$|\varphi(x, u, h)-\varphi(x, \bar{u}, h)| \leqslant L_{\varphi}|u-\bar{u}|$$

又设初值 是准确的，即 ，则其整体截断误差

$$u(x_{n})-u_{n}=O(h^{p})$$

数值方法误差阶的计算

设数值方法的误差阶

$$\left\|u\left(x_{n}\right)-u_{n}\right\|_*=C h^{p}+O\left(h^{p+1}\right).$$

省略 的高阶无穷小，记误差

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

两边都取对数

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

为 关于 的斜率

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

如果函数区间为 ，等距剖分为 等份，，则

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

线性 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.$$

方程的真解：.

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