简单的数值方法

发表于 2020-05-19 更新于 2023-10-19 分类于 Note 阅读次数： 46

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

常微分方程初值问题

考虑常微分方程初值问题，设

${\begin{aligned} \frac{d u}{d t} = f (t, u), 0 < t ⩽ T \\ u (0) = 0 \end{aligned}$

$f (t, u)$ .在区域 $G : 0 ⩽ t ⩽ T, | u | < \infty$ 上连续，求 $u = u (t)$ ：满足：

${\begin{aligned} \frac{d u}{d t} = f (t, u), 0 < t ⩽ T & u (0) = 0 \end{aligned}$

通常 $f$ 满足 Lipschitz 条件： $| f (t, u_{1}) - f (t, u_{2}) | ⩽ L | u_{1} - u_{2} |$

线性 ODE 例子

$$\left{ $\begin{aligned} \frac{d u}{d t} = t^{2} + t - u, t \in [0, 1] & u (0) = 0 \end{aligned}$ \right.$$

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

Euler 方法

$u_{n + 1} = u_{n} + h f (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;  clf
h=0.1;
x=0:h:1;                       % function interval
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','North')
set(gca,'fontsize',12)
xlabel('x','fontsize', 16), ylabel('u','fontsize',16,'Rotation',0)

% print -dpng -r600  Euler1.png

输出结果

后退 Euler 方法

$u_{n + 1} = u_{n} + h f (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;  clf
h=0.1;
x=0:h:1;                       % function interval
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','North')
set(gca,'fontsize',12)
xlabel('x','fontsize', 16), ylabel('u','fontsize',16,'Rotation',0)

% print -dpng -r600  Trapezoidal.png

输出结果

改进的 Euler 方法

预估校正

{\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}

MATLAB 程序

% EulerPro.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;  clf
h=0.1;
x=0:h:1;                     % function interval
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','North')
set(gca,'fontsize',12)
xlabel('x','fontsize', 16), ylabel('u','fontsize',16,'Rotation',0)

% print -dpng -r600  EulerPro.png

输出结果

参考书籍

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