简单的数值方法

本文列举了简单的数值方法：Euler 方法、后退 Euler 方法、梯形公式、改进 Euler 方法，并用 MATLAB 求解线性 ODE，最后展示了输出结果。

常微分方程初值问题

考虑 常微分方程初值问题，设 $f(t,u)$ 在区域 $\rm G$ : $0\leqslant t \leqslant T$, $|u|< \infty$ 上连续，求 $u=u(t)$ 满足：

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

通常 $f$ 满足 Lipschitz 条件： $|f(t,u_1)-f(t,u_2)|\leqslant L|u_1-u_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.$$

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

Euler 方法

$$u_{n+1}=u_{n}+hf(t_n,u_n).$$

Euler 方法数值求解

MATLAB 程序

% Euler1.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;
h=0.1;
x=0:h:1;                       % interval partition
N=length(x)-1;
u(1)=0;                        % initial value
fun=@(t,u) t.^2+t-u;           % RHS
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
plot(x,ue,'b-',x,u,'r+','LineWidth',1)
legend('Exact','Numerical','location','northwest')
% title('Euler method','fontsize',12)
set(gca,'fontsize',12)
xlabel('x','fontsize',16), ylabel('u','fontsize',16)

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

输出结果

后退 Euler 方法

$$u_{n+1}=u_{n}+hf(t_{n+1},u_{n+1}).$$

梯形公式

$$u_{n+1}=u_{n}+\frac{h}{2} [f(t_n,u_n)+f(t_{n+1},u_{n+1})].$$

梯形公式数值求解

梯形公式与后退 Euler 方法类似，这里考虑梯形公式，对于线性 ODE 例子：

$$u_{n+1}=u_{n}+\frac{h}{2}[t_n^2+t_n-u_n+t_{n+1}^2+t_{n+1}-u_{n+1}],$$

可得

$$u_{n+1}=\frac{2-h}{2+h} u_n+\frac{h}{2+h}[t_n^2+t_n+t_{n+1}^2+t_{n+1}].$$

MATLAB 程序

% Trapezoidal.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; close all;
h=0.1;
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
plot(x,ue,'b-',x,u,'r+','LineWidth',1)
legend('Exact','Numerical','location','northwest')
%title('Trapezoidal rule','fontsize',12)
set(gca,'fontsize',12)
xlabel('x','fontsize',16), ylabel('u','fontsize',16)

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

输出结果

改进的 Euler 方法

预估校正

$$\left\{\begin{aligned} & \bar{y}_{n+1}=y_{n}+h f(x_{n}, y_{n}), \\ & y_{n+1}=y_{n}+\frac{h}{2}[f(x_{n}, y_{n})+f(x_{n+1}, \bar{y}_{n+1}] \end{aligned}\right.$$

MATLAB 程序

% ModiEuler.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;
h=0.1;
x=0:h:1;                     % interval partition
N=length(x)-1;
u(1)=0;                      % initial value
fun=@(x,u) x.^2+x-u;         % RHS
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
plot(x,ue,'b-',x,u,'r+','LineWidth',1)
legend('Exact','Numerical','location','northwest')
%title('Modified Euler Method','fontsize',12)
set(gca,'fontsize',12)
xlabel('x','fontsize',16), ylabel('u','fontsize',16)

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

输出结果

参考书籍

• 数值分析 第 5 版 （李庆扬等）
• MATLAB 微分方程高效解法