数学 - 微分方程数值解 - 第 1 章一阶常微分方程初值问题 - 1.2 Euler 方法及其改进方法 - Black_x

数学 - 微分方程数值解 - 第 1 章一阶常微分方程初值问题 - 1.2 Euler 方法及其改进方法

1.2 Euler 方法及其改进方法

1.2.1 Euler 方法

用 $f(x_n, y_n)$ 代替式 $(1.2)$ 中的 $\varphi_n$ ，得到差分方程初值问题：

\begin{matrix} (1.2.1) & {\begin{aligned} y_{n + 1} = y_{n} + h f (x_{n}, y_{n}) \\ y_{0} = y (a) \end{aligned} n = 0, 1, \dots \end{matrix}

$\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$ ，得到差分方程初值问题：

\begin{matrix} (1.2.2) & {\begin{aligned} y_{n + 1} = y_{n} + h f (x_{n + 1}, y_{n + 1}) \\ y_{0} = y (a) \end{aligned} n = 0, 1, \dots \end{matrix}

$\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}$ 。使用简单迭代法可以构造迭代格式如下：

\begin{matrix} (1.2.3) & {\begin{aligned} 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)}) k = 0, 1, \dots \end{aligned} \end{matrix}

$\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 条件，即

| f (x, y_{1}) - f (x, y_{2}) | ⩽ L | y_{1} - y_{2} |

$\left| f(x,y_1) - f(x,y_2) \right| \leqslant L\left| y_1 - y_2 \right|$

那么可得到迭代收敛条件：

\begin{matrix} (1.2.4) & \begin{aligned} | y_{n + 1}^{(k + 1)} - y_{n + 1}^{(k)} | & = h | f (x_{n + 1}, y_{n + 1}^{(k)}) - f (x_{n + 1}, y_{n + 1}^{(k - 1)}) | \\ ⩽ h L | y_{n + 1}^{(k)} - y_{n + 1}^{(k - 1)} | ⩽ \dots \\ ⩽ (h L)^{k} | y_{n + 1}^{(1)} - y_{n + 1}^{(0)} | \end{aligned} \end{matrix}

$\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 梯形公式

考虑下列差分方程初值问题：

\begin{matrix} (1.2.5) & {\begin{aligned} y_{n + 1} = y_{n} + \frac{h}{2} [f (x_{n}, y_{n}) + f (x_{n + 1}, y_{n + 1})] \\ y_{0} = y (a) \end{aligned} n = 0, 1, \dots \end{matrix}

$\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}$ 。使用简单迭代法可以构造迭代格式如下：

\begin{matrix} (1.2.6) & {\begin{aligned} y_{n + 1}^{(0)} = y_{n} + h f (x_{n}, y_{n}) \\ y_{n + 1}^{(k + 1)} = y_{n} + \frac{h}{2} [f (x_{n}, y_{n}) + f (x_{n + 1}, y_{n + 1}^{(k)})] k = 0, 1, \dots \end{aligned} \end{matrix}

$\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{matrix} (1.2.7) & \begin{aligned} | y_{n + 1}^{(k + 1)} - y_{n + 1}^{(k)} | & = \frac{h}{2} | f (x_{n + 1}, y_{n + 1}^{(k)}) - f (x_{n + 1}, y_{n + 1}^{(k - 1)}) | \\ ⩽ \frac{h L}{2} | y_{n + 1}^{(k)} - y_{n + 1}^{(k - 1)} | ⩽ \dots \\ ⩽ {(\frac{h L}{2})}^{k} | y_{n + 1}^{(1)} - y_{n + 1}^{(0)} | \end{aligned} \end{matrix}

$\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 法和梯形公式做一个预测 - 校正系统。

\begin{matrix} (1.2.8) & {\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} \end{matrix}

$\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)$ 的第一个式子为预测，第二个式子为校正。

实际计算中，改写成下式：

\begin{matrix} (1.2.9) & {\begin{aligned} 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{aligned} n = 0, 1, \dots, \end{matrix}

$\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{aligned} R_{n + 1} & = y (x_{n + 1}) - {\bar{y}}_{n + 1} \\ = y (x_{n + 1}) - (y (x_{n}) + h f (x_{n}, y_{n})) \\ = y (x_{n + 1}) - (y (x_{n}) + h y^{'} (x_{n})) \\ = y (x_{n}) + h y^{'} (x_{n}) + h^{2} y^{″} (ξ) / 2 + o (h^{2}) - (y (x_{n}) + h y^{'} (x_{n})) \\ (1.2.10) & = h^{2} y^{″} (ξ) / 2 + o (h^{2}) x_{n} < ξ < x_{n + 1} \end{aligned}

$\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{aligned} | e_{n + 1} | & = | y (x_{n + 1}) - y_{n + 1} | \\ ⩽ | y (x_{n + 1}) - {\bar{y}}_{n + 1} | + | {\bar{y}}_{n + 1} - y_{n + 1} | \\ ⩽ | R_{n + 1} | + | y (x_{n}) - y_{n} | + h | f (x_{n}, y (x_{n})) - f (x_{n}, y_{n}) | \\ ⩽ | R_{n + 1} | + | y (x_{n}) - y_{n} | + h L | y (x_{n}) - y_{n} | \\ (1.2.11) & = | R_{n + 1} | + (1 + h L) | e_{n} | \end{aligned}

$\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{aligned} | e_{n + 1} | & ⩽ M h^{2} / 2 + (1 + h L) | e_{n} | \\ ⩽ M h^{2} / 2 + (1 + h L) [M h^{2} / 2 + (1 + h L) | e_{n - 1} |] \\ ⩽ M h^{2} / 2 [1 + (1 + h L) + (1 + h L)^{2} + \dots + (1 + h L)^{n}] \\ = \frac{h M}{2 L} [(1 + h L)^{n + 1} - 1] \end{aligned}

$\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 + h L)^{n + 1} ⩽ (1 + h L)^{\frac{b - a}{h}} < e^{L (b - a)}

$(1 + hL)^{n+1} \leqslant (1 + hL)^{\frac{b-a}{h}} < e^{L(b-a)}$

因此有

| e_{n + 1} | ⩽ \frac{h M}{2 L} [e^{L (b - a)} - 1] = O (h)

$\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{aligned} R_{n + 1} & = y (x_{n + 1}) - {\bar{y}}_{n + 1} \\ = y (x_{n + 1}) - (y (x_{n}) + h f (x_{n}, y_{n})) \\ = y (x_{n + 1}) - (y (x_{n}) + h y^{'} (x_{n})) \\ = y (x_{n}) + h y^{'} (x_{n}) + h^{2} y^{″} (x_{n}) / 2 + o (h^{2}) - (y (x_{n}) + h y^{'} (x_{n})) \\ (1.2.12) & = h^{2} y^{″} (x_{n}) / 2 + O (h^{3}) x_{n} < ξ < x_{n + 1} \end{aligned}

$\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)} (ξ_{1}) x_{n} < ξ_{1} < 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{aligned} {\bar{y}}_{n + 1} & = y (x_{n}) + \frac{h}{2} (f (x_{n}, y (x_{n})) + f (x_{n + 1}, y (x_{n + 1}))) \\ = 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)} (ξ_{2})) x_{n} < ξ_{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)} (ξ_{2}) x_{n} < ξ_{2} < x_{n + 1} \end{aligned}

$\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{aligned} R_{n + 1} = y (x_{n + 1}) - {\bar{y}}_{n + 1} = - \frac{h^{3}}{12} y^{‴} (x_{n}) + O (h^{4}) \end{aligned}

$\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*}$

posted on 2022-03-01 00:10 Black_x 阅读(1397) 评论(0) 编辑收藏举报

刷新页面返回顶部

登录后才能查看或发表评论，立即登录或者逛逛博客园首页

阅读排行：
· 分享一个免费、快速、无限量使用的满血 DeepSeek R1 模型，支持深度思考和联网搜索！
· 基于 Docker 搭建 FRP 内网穿透开源项目（很简单哒）
· 25岁的心里话
· ollama系列01：轻松3步本地部署deepseek，普通电脑可用
· 按钮权限的设计及实现

Black_x