1.2 Euler 方法及其改进方法
1.2.1 Euler 方法
用 \(f(x_n, y_n)\) 代替式 \((1.2)\) 中的 \(\varphi_n\),得到差分方程初值问题:
\[\left\{
\begin{align*}
& y_{n+1} = y_{n} + h f(x_n, y_n) \\
& y_0 = y(a)
\end{align*} \quad n=0,1,\cdots
\right. \tag{1.2.1}
\]
以上式问题的解作为微分方程初值问题的数值解,即 \(y(x_n) = y_n\),称为 Euler 方法。容易看出式 \((1.3)\) 是一个单步显式格式。
1.2.2 向后 Euler 法
用 \(f(x_{n+1}, y_{n+1})\) 代替式 \((1.2)\) 中的 \(\varphi_n\),得到差分方程初值问题:
\[\left\{
\begin{align*}
& y_{n+1} = y_{n} + h f(x_{n+1}, y_{n+1}) \\
& y_0 = y(a)
\end{align*} \quad n=0,1,\cdots
\right. \tag{1.2.2}
\]
以上式问题的解作为微分方程初值问题的数值解,即 \(y(x_n) = y_n\),称为向后 Euler 方法。容易看出式 \((1.4)\) 是一个单步隐式格式,需要迭代法求解 \(y_{n+1}\)。使用简单迭代法可以构造迭代格式如下:
\[\left\{
\begin{align*}
& y_{n+1}^{(0)} = y_n + h f(x_n,y_n) \\
& y_{n+1}^{(k+1)} = y_n + h f(x_n,y_{n+1}^{(k)}) \quad k=0,1,\cdots
\end{align*}
\right. \tag{1.2.3}
\]
下面考虑迭代格式 \((1.5)\) 的收敛性,由于 \(f(x,y)\) 关于 \(y\) 满足 Lipschitz 条件,即
\[\left| f(x,y_1) - f(x,y_2) \right| \leqslant L\left| y_1 - y_2 \right|
\]
那么可得到迭代收敛条件:
\[\begin{align*}
\left| y_{n+1}^{(k+1)} - y_{n+1}^{(k)} \right|
& = h \left| f(x_{n+1},y_{n+1}^{(k)}) - f(x_{n+1},y_{n+1}^{(k-1)}) \right| \\
& \leqslant h L \left| y_{n+1}^{(k)} - y_{n+1}^{(k - 1)} \right| \leqslant \cdots \\
& \leqslant (hL)^{k} \left| y_{n+1}^{(1)} - y_{n+1}^{(0)} \right|
\end{align*} \tag{1.2.4}
\]
因此向后 Euler 方法迭代条件为:\(hL < 1\)。
1.2.3 梯形公式
考虑下列差分方程初值问题:
\[\left\{
\begin{align*}
& y_{n+1} = y_{n} + \frac{h}{2} \left[ f(x_n, y_n) + f(x_{n+1}, y_{n+1}) \right] \\
& y_0 = y(a)
\end{align*} \quad n=0,1,\cdots
\right. \tag{1.2.5}
\]
以上式问题的解作为微分方程初值问题的数值解,即 \(y(x_n) = y_n\),称为梯形公式。容易看出式 \((1.4)\) 是一个单步隐式格式,需要迭代法求解 \(y_{n+1}\)。使用简单迭代法可以构造迭代格式如下:
\[\left\{
\begin{align*}
& y_{n+1}^{(0)} = y_n + h f(x_n,y_n) \\
& y_{n+1}^{(k+1)} = y_n + \frac{h}{2} \left[ f(x_n, y_n) + f(x_{n+1}, y_{n+1}^{(k)}) \right] \quad k=0,1,\cdots
\end{align*}
\right. \tag{1.2.6}
\]
下面考虑迭代格式 \((1.8)\) 的收敛性。
\[\begin{align*}
\left| y_{n+1}^{(k+1)} - y_{n+1}^{(k)} \right|
& = \frac{h}{2} \left| f(x_{n+1},y_{n+1}^{(k)}) - f(x_{n+1},y_{n+1}^{(k-1)}) \right| \\
& \leqslant \frac{h L}{2} \left| y_{n+1}^{(k)} - y_{n+1}^{(k - 1)} \right| \leqslant \cdots \\
& \leqslant \left( \frac{hL}{2} \right)^{k} \left| y_{n+1}^{(1)} - y_{n+1}^{(0)} \right|
\end{align*} \tag{1.2.7}
\]
得梯形公式的收敛条件:\(hL/2 < 1\)。
1.2.4 改进 Euler 法
将 Euler 法和梯形公式做一个预测 - 校正系统。
\[\left\{
\begin{align*}
& \overline{y}_{n+1} = y_{n} + h f(x_n, y_n) \\
& y_{n+1}= y_{n} + \frac{h}{2} \left[ f(x_n, y_n) + f(x_{n+1}, \overline{y}_{n+1}) \right]
\end{align*}
\right. \tag{1.2.8}
\]
其中,式 \((1.2.8)\) 的第一个式子为预测,第二个式子为校正。
实际计算中,改写成下式:
\[\left\{
\begin{align*}
& y_p = y_{n} + h f(x_n, y_n) \\
& y_q = y_{n} + h f(x_{n+1}, y_p) \\
& y_{n+1} = (y_p + y_q) / 2
\end{align*} \quad n=0,1,\cdots,
\right. \tag{1.2.9}
\]
1.2.5 单步法的误差估计
单步法一般形式: \(y_{n+1} = y_n + h \varphi_n\)
给出局部截断误差的定义
定义 1.2.1 局部截断误差
利用 \(x_i \, (i \leqslant n)\),假定 \(y(x)\) 已知,利用公式递推一步计算 \(y_{n+1}\),与 \(y(x_{n+1})\) 的差的即为局部截断误差。
\(R_{n+1} = y(x_{n+1}) - \overline{y}_{n+1}\),其中 \(\overline{y}_{n+1} = y(x_n) + h \varphi_n\)
给出整体截断误差的定义
定义 1.2.2 整体截断误差
利用公式从最开始 \(x_0\) 一步一步递推计算 \(y_{n+1}\),与 \(y(x_{n+1})\) 的差的即为整体截断误差。
\(e_{n+1} = y(x_{n+1}) - y_{n+1}\)
给出 \(p\) 阶方法的定义:
定义 1.2.3 \(p\) 阶方法
若某种数值方法的局部截断误差 \(R_{n+1} = O(h^{p+1})\),则称该方法是 \(p\) 阶方法。若进一步满足下式
\[R_{n+1} = \psi(x_n, y(x_n)) h^{p+1} + O(h^{p+2})
\]
则称 \(\psi(x_n, y(x_n)) h^{p+1}\) 为局部截断误差的主项。
(1) Euler 方法的截断误差
按局部截断误差的定义计算 Euler 方法的局部截断误差:
\[\begin{align*}
R_{n+1}
& = y(x_{n+1}) - \overline{y}_{n+1} \\
& = y(x_{n+1}) - \left( y(x_n) + h f(x_n, y_n) \right) \\
& = y(x_{n+1}) - \left( y(x_n) + h y'(x_n) \right) \\
& = y(x_{n}) + h y'(x_n) + h^2 y''(\xi) / 2 + o(h^2) - \left( y(x_n) + h y'(x_n) \right) \\
& = h^2 y''(\xi) / 2 + o(h^2) \quad x_n < \xi < x_{n+1} \tag{1.2.10}
\end{align*}
\]
故 Euler 方法为一阶方法。
按整体截断误差的定义计算 Euler 方法的整体截断误差:
\[\begin{align*}
\left| e_{n+1} \right|
& = \left| y(x_{n+1}) - y_{n+1} \right|\\
& \leqslant \left| y(x_{n+1}) - \overline{y}_{n+1} \right| + \left| \overline{y}_{n+1} - y_{n+1} \right| \\
& \leqslant \left| R_{n+1} \right| + \left| y(x_n) - y_n \right| + h \left| f(x_n, y(x_n)) - f(x_n, y_n) \right| \\
& \leqslant \left| R_{n+1} \right| + \left| y(x_n) - y_n \right| + h L \left| y(x_n) - y_n \right|\\
& = \left| R_{n+1} \right| + (1 + h L)\left| e_n \right| \tag{1.2.11}
\end{align*}
\]
我们进一步假设 \(f(x,y)\) 充分光滑,有 \(y(x) \in C^{(2)}[a, b]\),则存在 \(M > 0\),使得 \(y''(x) \leqslant M\),此时有 \(\left| R_k \right| \leqslant M h^2 / 2\)。我们有 \(e_1 = R_1\),进一步计算整体截断误差:
\[\begin{align*}
\left| e_{n+1} \right|
& \leqslant M h^2 / 2 + (1 + h L)\left| e_n \right| \\
& \leqslant M h^2 / 2 + (1 + h L)\left[ M h^2 / 2 + (1 + h L)\left| e_{n-1} \right| \right] \\
& \leqslant M h^2 / 2 \left[1 + (1 + h L) + (1 + h L)^2 + \cdots + (1 + h L)^n \right] \\
& = \frac{hM}{2L} \left[ (1+hL)^{n+1} - 1 \right] \\
\end{align*}
\]
因为 \((n+1)h \leqslant b - a\),所以
\[(1 + hL)^{n+1} \leqslant (1 + hL)^{\frac{b-a}{h}} < e^{L(b-a)}
\]
因此有
\[\left| e_{n+1} \right| \leqslant \frac{hM}{2L} \left[ e^{L(b-a)} - 1 \right] = O(h)
\]
可见 Euler 方法的整体截断误差与 \(h\) 同阶,为一阶方法。且 \(h \to 0\) 时,整体截断误差 \(e \to 0\)。
(2) 向后 Euler 方法的截断误差
按局部截断误差的定义计算向后 Euler 方法的局部截断误差:
\[\begin{align*}
R_{n+1}
& = y(x_{n+1}) - \overline{y}_{n+1} \\
& = y(x_{n+1}) - \left( y(x_n) + h f(x_n, y_n) \right) \\
& = y(x_{n+1}) - \left( y(x_n) + h y'(x_n) \right) \\
& = y(x_{n}) + h y'(x_n) + h^2 y''(x_n) / 2 + o(h^2) - \left( y(x_n) + h y'(x_n) \right) \\
& = h^2 y''(x_n) / 2 + O(h^3) \quad x_n < \xi < x_{n+1} \tag{1.2.12}
\end{align*}
\]
(3) 梯形公式的截断误差
按局部截断误差的定义计算梯形公式的局部截断误差。
考虑 \(y(x_{n+1})\)
\[y(x_{n+1}) = y(x_n) + h y'(x_n) + \frac{h^2}{2} y''(x_n) + \frac{h^3}{6} y'''(x_n) + \frac{h^4}{24} y^{(4)}(\xi_1) \quad x_n < \xi_1 <x_{n+1}
\]
考虑 \(\overline{y}_{n+1}\)
\[\begin{align*}
\overline{y}_{n+1}
& = y(x_n) + \frac{h}{2} \left( f(x_n, y(x_n)) + f(x_{n+1}, y(x_{n+1})) \right)\\
& = y(x_n) + \frac{h}{2} ( y'(x_n) + y'(x_{n+1}) ) \\
& = y(x_n) + \frac{h}{2} ( y'(x_n) + y'(x_n) + h y''(x_n) + \frac{h^2}{2} y'''(x_n) + \frac{h^3}{6} y^{(4)}(\xi_2) ) \quad x_n < \xi_2 <x_{n+1} \\
& = y(x_n) + h y'(x_n) + \frac{h^2}{2} y''(x_n) + \frac{h^3}{4} y'''(x_n) + \frac{h^4}{12} y^{(4)}(\xi_2) \quad x_n < \xi_2 <x_{n+1} \\
\end{align*}
\]
可得梯形公式的局部截断误差:
\[\begin{align*}
R_{n+1}
= y(x_{n+1}) - \overline{y}_{n+1} = - \frac{h^3}{12} y'''(x_n) + O(h^4)
\end{align*}
\]