Euler公式

数值计算

张瑞
中国科学技术大学数学科学学院

Euler公式

实际问题得到的数学模型往往都是微分方程。如： 物种的增长和蜕变， 物体的运动， 传染病的传播等。

  • 函数的自变量只有一个的微分方程称为常微分方程。
  • 函数的自变量为两个或两个以上的微分方程称为偏微分方程。

微分方程的解通常不唯一，因此，通常会加上一个或多个的附加条件。

  • 给定微分方程及其初始条件，称为初值问题。
  • 给定微分方程及其边界条件，称为边值问题。

下面讨论常微分方程的初值问题

\[\begin{cases} y'(x)=f(x,y) \\ y(a)=y_0 \end{cases} , x\in[a,b] \]

为保证解的存在唯一性及稳定性，要求$f(x,y)$满足对$y$的Lipschitz条件， 即存在$L>0$，使得对任意$y_1$, $y_2$，有

\[|f(x,y_1)-f(x,y_2)|<L|y_1-y_2| \]

只有一些特殊形式的$f(x,y)$才能得到方程的解析解，对于绝大多数常微分方程的初值问题， 只能计算它的数值解。

  • 微分方程的解是一个函数，而计算机没办法表达任意的函数
  • 一种有效的表示方法，就是给出函数在一系列点上的函数值。
    1. 这样，可以先将区间$[a,b]$作一个剖分，得到点列
      \[\{a=x_0, x_1, \cdots, x_n=b\} \]
    2. 利用微分方程，构造算法求解$y(x_i)$的值。
    3. 这种方法称为有限差分方法。数值计算中得到的通常是$y(x_i)$的近似值$y_i$，也称 $\{y_0, y_1, \cdots, y_n\}$为格点函数。
  • 这种将区间分割为小区间的方式，称为数值离散方法。它是求解微分方程数值解的基本手段。
  • 除了有限差分法外，还有有限体积法，有限元法都是求解微分方程的数值离散方法。

向前Euler公式

简单起见，将区间$[a,b]$做等距分割

\[\{x_i=a+i h, i=0,1,\cdots,n \}, \quad h=\frac{b-a}n \]

导数$y'(x_i)$可以用差商公式来近似

\[y'(x_i)=\frac{y(x_{i+1})-y(x_i)}h-\frac{h}2f''(\xi_i), \xi_i\in(x_i, x_{i+1}) \]

由微分方程$y'(x_i)=f(x_i, y(x_i))$，因而有

\[f(x_i, y(x_i))=\frac{y(x_{i+1})-y(x_i)}h-\frac{h}2f''(\xi_i) \]

整理后，有

\[y(x_{i+1})=y(x_i)+hf(x_i,y(x_i))+\frac{h^2}2f''(\xi_i) \]

得到计算$y(x_{i+1})$近似值$y_{i+1}$的公式

\[y_{i+1}=y_i+h f(x_i, y_i) \]

这个公式称为向前Euler公式，或向前差商公式。

利用初值，给出$y_0=y(x_0)=y(a)$后，代入上式即可以顺序得到$y_1$, $y_2$, $\cdots$, $y_n$。

定义 1.
关于格点函数的方程，称为差分方程。

例 1. 用向前Euler格式解方程

\[\begin{cases} y'(x)=-\sin(x) , \\ y(0)=1 \end{cases}, x\in[0,\pi] \]

解. 如图

num-ode-euler-forward

向前Euler公式的收敛性

由

\[y(x_{i+1})=y(x_i)+hf(x_i,y(x_i))+\frac{h^2}2f''(\xi_i) \]

到向前Euler公式

\[y_{i+1}=y_i+h f(x_i, y_i) \]

中丢掉的项$\frac{h^2}2f''(\xi_i)$是格式的局部截断误差。 它是假定前面的计算是精确的($y_i=y(x_i)$)，用格式计算下一个点上的函数值时产生的误差。

  • 实际上， 从计算$y_1$开始，每个$y_i$都有截断误差，
  • 而这个误差加上前面所有的点的误差， 会传播累积到$y_{i+1}$。
  • 这些局部截断误差传播、累积起来得到的就是整体截断误差。
\[\begin{aligned} e_{i+1}=&y(x_{i+1})-y_{i+1} \\ =&\bigg(y(x_i)+hf(x_i,y(x_i))-\frac{h^2}2f''(\xi_i)\bigg) -\bigg(y_i+h f(x_i, y_i)\bigg) \\ =& e_i+h\bigg(f(x_i,y(x_i))-f(x_i, y_i)\bigg)-\frac{h^2}2f''(\xi_i) \end{aligned} \]

记$T_{i+1}=-\frac{h^2}2f''(\xi_i)$是局部截断误差。 由$f(x,y)$满足对$y$的Lipschitz条件，则有

\[\begin{aligned} |e_{i+1}|\leq & |e_i|+h|f(x_i,y(x_i))-f(x_i, y_i)|+|T_{i+1}| \\ \leq & |e_i|+hL|y(x_i)-y_i|+|T_{i+1}| \\ \leq & |e_i|(1+hL)+|T_{i+1}| \end{aligned} \]

记$T=\max_{i}|T_i|=O(h^2)$，则

\[|e_{i+1}|\leq |e_i|(1+hL)+T \]

因而，

\[\begin{aligned} |e_{i+1}|\leq &(1+hL)|e_i|+T \\ \leq & (1+hL)((1+hL)|e_{i-1}|+T)+T\\ = & (1+hL)^2|e_{i-1}|+T(1+(1+hL)) \\ \leq & (1+hL)^2((1+hL)|e_{i-2}|+T)+T(1+(1+hL)) \\ \leq & (1+hL)^{3}|e_{i-2}|+T(1+(1+hL)+(1+hL)^3) \\ \leq & \cdots \\ \leq & (1+hL)^{i+1}|e_{0}|+T(1+(1+hL)+\cdots+(1+hL)^i) \\ =&(1+hL)^{i+1}|e_{0}|+T\frac{1-(1+hL)^{i+1}}{-hL} \\ \leq & (1+hL)^{i+1}\left(|e_{0}|+\frac{T}{hL}\right) \end{aligned} \]

由$(1+x)^n<e^{nx}$, $\forall x>0, n>0$，因而有

\[\begin{aligned} |e_{i+1}|\leq & (1+hL)^{i+1}\left(|e_{0}|+\frac{T}{hL}\right) \\ \leq & e^{(i+1)hL}\left(|e_{0}|+\frac{T}{hL}\right) \\ \leq & e^{nhL}\left(|e_{0}|+\frac{T}{hL}\right) \\ =&e^{(b-a)L}\left(|e_{0}|+\frac{T}{hL}\right) \end{aligned} \]

且$e^0=y(x_0)-y_0$是初值给定的，因而$e_0=0$，因此

向前公式的整体截断误差为

\[|e^{i}|\leq \frac{e^{(b-a)L}}{L}\frac{T}h=O(h) \]

可以看到：

  1. $h\to0$时，$e_i\to 0$，即向前Euler公式是收敛的。
  2. 公式的局部截断误差$T=O(h^2)$，而整体截断误差为$O(h)$。

定义 2.
若公式的局部截断误差是

\[T_{i}=O(h^{p+1}) \]

则称方法是p阶的，或具有p阶精度。

向后Euler公式

类似向前Euler公式，用向后差商来给出数值微分估计

\[y'(x_{i+1})=\frac{y(x_{i+1})-y(x_i)}{h}+\frac{h^2}2f''(\xi_i) \]

利用微分方程，有

\[f(x_{i+1},y(x_{i+1}))=\frac{y(x_{i+1})-y(x_i)}{h}+\frac{h^2}2f''(\xi_i) \]

整理得

\[y(x_{i+1})=y(x_i)+hf(x_{i+1}, y(x_{i+1}))-\frac{h^2}2f''(\xi_i) \]

得到格式

格式

\[y_{i+1}=y_i+f(x_{i+1}, y_{i+1}) \]

称为向后Euler公式。

该公式的局部截断误差为

\[T_{i+1}=-\frac{f''(\xi_i)}2 h^2=O(h^2) \]

因而，也是1阶收敛格式。

注意到向后Euler公式中

\[y_{i+1}=y_i+f(x_{i+1}, y_{i+1}) \]

$y_{i+1}$在等号的两边。在已知$y_i$时，要得到$y_{i+1}$可能需要解一个非线性方程(如$f(x,y)=\sin(y)$)。

  • 格式中$y_{i+1}$在等号两边，不能直接由前面的计算值得到$y_{i+1}$的格式，称为隐格式。
  • 相应地，可以直接由前面的计算值得到$y_{i+1}$的格式，称为显格式。 如：向前Euler公式。

在向后Euler公式中

\[y_{i+1}=y_i+f(x_{i+1}, y_{i+1}) \]

为了得到$y_{i+1}$，

  • 可以使用解非线性方程的Newton迭代，
  • 更常使用的方法就是直接使用格式构造迭代
    \[y_{i+1}^{(k+1)}=y_i+f(x_{i+1}, y_{i+1}^{(k)}) \]
  • 初始估计值可以使用显格式（向前Euler公式）得到。 由于显格式也是对$y_{i+1}$的估计，因此，这种迭代通常不超过3步即可。

在实际使用中，可以用向前Euler公式做初始值估计，用向后Euler格式做一步迭代，得到

\[\begin{cases} y_{i+1}^{(0)}=&y_i+h f(x_i, y_i) \\ y_{i+1} = &y_i+hf(x_{i+1}, y_{i+1}^{(0)}) \end{cases} \]

利用微分方程的初值$y_0=y(a)$，就可以依次得到$y_1$, $y_2$, $\cdots$, $y_n$。

定义 3.
这种用显格式做初始值估计，用隐格式做一步迭代的格式又称为预估-校正格式。

注. 为了保证局部截断误差的精度，预估-校正类格式中，显格式的精度最多只能比隐格式低1阶。

事实上，对预估-校正格式有

\[y_{i+1}=y_i+h f(x_{i+1}, y_i+hf(x_i,y_i)) \]

而

\[y(x_{i+1})=y(x_i)+hf(x_{i+1}, y(x_{i+1}))-\frac{h^2}2f''(\xi_i) \]

当$y_i=y(x_i)$时，

\[\begin{aligned} |y(x_{i+1})-y_{i+1}|=&h|f(x_{i+1}, y(x_{i+1}))-f(x_{i+1}, y_i+hf(x_i,y_i))|+T_{Im} \\ \leq &hL|y(x_{i+1})-(y_i+hf(x_i,y_i))|+T_{Im} \end{aligned} \]

其中$T_{Im}$是隐格式的局部截断误差。

而

\[|y(x_{i+1})-(y_i+hf(x_i,y_i))| \]

正好是显格式的局部截断误差，记为$T_{Ex}$。

则预估-校正格式的误差是

\[|e_{i+1}|=hL T_{Ex}+T_{Im} \]

也就是说，$T_{Ex}$的误差阶可以比$T_{Im}$的误差阶低1阶。

例 2. 用向后Euler格式解方程

\[\begin{cases} y'(x)=-\sin(x) , \\ y(0)=1 \end{cases}, x\in[0,\pi] \]

解. 如图

num-ode-euler-backward

中心Euler公式

用中心差商近似导数

\[y'(x_{i})=\frac{y_{i+1}-y_{i-1}}{2h}-\frac{h^2}6 f'''(\xi_i) \]

代入微分方程，可以得到中心差商公式

\[y_{i+1}=y_{i-1}+2hf(x_i, y_i) \]
  • 分析该格式的局部截断误差时，要假定前面的计算是精确的，即$y_{i-1}=y(x_{i-1})$, $y_i=y(x_i)$。
  • 易得，局部截断误差是$\frac{f'''(\xi_i)}3 h^3=O(h^3)$。 中心差商公式是2阶精度格式。

在中心差商格式中

\[y_{i+1}=y_{i-1}+2hf(x_i, y_i) \]
  • $y_{i+1}$除了跟$y_i$相关外，还跟$y_{i-1}$相关。这类格式也称为多步格式。
  • $y_{i+1}$只跟$y_i$相关的格式称为单步格式。如： 向前Euler公式、向后Euler公式。
  • 在多步格式中，只有初值$y_0$是不够的，还需要用其它的格式计算$y_1$等，才能使用格式开始计算。 这个过程称为多步格式的起步计算。
  • 为了保证整体截断误差的精度， 起步计算所用格式的精度最多只能比格式的精度低1阶。

例 3. 用中心差商格式解方程

\[\begin{cases} y'(x)=-\sin(x) , \\ y(0)=1 \end{cases}, x\in[0,\pi] \]

解. 如图

num-ode-euler-center

谢谢

vertical slide 2

目录

本节读完

例 4. 证明$(1+x)^n<e^{nx}$, $\forall x>-1$, $n\in\mathbb{Z}^+$

4. 令$f(x)=e^x-x-1$，则

\[f'(x)=e^x-1, \quad f''(x)=e^x \]

因而$f(0)$是最小值，即有

\[f(x)> f(0), \forall x\neq 0 \]

即

\[e^x>x+1, \forall x\neq 0 \]

由$1+x<e^x$, $\forall x\neq 0$，则

\[(1+x)^n<(e^x)^n=e^{nx}, \forall x\geq -1, n\in\mathbb{Z}^+ \]