常微分方程欧拉折线法的Picard迭代改进

Improvement of Euler Method for Ordinary Differential Equations by Picard Iteration

摘要: 考虑常微分方程数值解法，对欧拉折线法加以改进。主要方法是，对每个节点处局部线性近似解进行一次Picard迭代修正，用局部修正解估算下一节点处解值。证明了算法局部截断误差比欧拉折线法高一阶，分析了算法的收敛性与稳定性。通过方程实例数值计算进一步验证了本文算法所得近似解比欧拉折线法近似解更接近精确解。

Abstract: In this paper we study numerical algorithms for ordinary differential equations, and improve the Euler method. We refine local linear approximate solution at each node by Picard iteration, and then estimate the numerical solution at next node by the refined local approximate solution. The local truncation error of our algorithm is smaller than that of Euler method by one order. The con-vergence and stability are proven. Finally, we show that our algorithm gives approximate solution closer to the exact solution than Euler method by an example.

1. 引言

一阶常微分方程初值问题形如

$\{\begin{array}{l}\frac{\text{d}y}{\text{d}x}=f\left(x,y\right),\\ {y|}_{x=x_0}=y_0.\end{array}$ (1.1)

为保证解的存在唯一性， $f\left(x,y\right)$ 在含点 $\left(x_0,y_0\right)$ 的某区域连续，且关于y局部Lipschitz连续。当 $f\left(x,y\right)$ 表达式较复杂时，难以用分离变量积分等办法求出解析解，只能通过数值方法计算近似解了。

数值方法主要思想是，用差商近似代替微商，将方程差分化，离散化。经典方法主要有Euler方法和Runge-Kutta方法。Euler方法的简单体现形式为显式Euler公式，即Euler折线法：

$y_{n+1}=y_n+hf\left(x_n,y_n\right)$ , $n=0,1,\cdots$ .

其局部截断误差为 $O\left(h^2\right)$ ，其中h为离散化步长。

若用向后差商 $\frac{y\left(x_{n+1}\right)-y\left(x_n\right)}{h}$ 近似代替导数 $y^{\prime }\left(x_{n+1}\right)$ ，还可得隐式或后退Euler公式：

$y_{n+1}=y_n+hf\left(x_{n+1},y_{n+1}\right)$ .

其局部截断误差也是 $O\left(h^2\right)$ 。具体计算时要先用显式Euler公式给出初步近似解，再由隐式公式进行迭代修正。

Euler方法的高级体现形式主要为梯形公式和预估–校正公式。梯形公式可表述为

$y_{n+1}=y_n+\frac{h}{2}\left[f\left(x_n,y_n\right)+f\left(x_{n+1},y_{n+1}\right)\right]$ ,

其局部截断误差为 $O\left(h^3\right)$ 。由于是隐式，具体计算时也需先用显式Euler公式给出初步近似解，再由隐式公式进行迭代修正。

预估–校正公式即对显式Euler公式所给初步近似解进行一次迭代修正，可表述为

${\bar{y}}_{n+1}=y_n+hf\left(x_n,y_n\right)$ ,

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

其局部截断误差为 $O\left(h^3\right)$ 。

Runge-Kutta算法思想是，用几个点处f预估值加权平均作为当前节点处 $y^{\prime }$ 估计值，以提高计算精度。

为了进一步提高计算精度，还可利用当前节点 $x_{n+1}$ 前面多个节点 $x_n,x_{n-1},\cdots ,x_{n-r}$ 处的数值结果 $y_n,y_{n-1},\cdots ,y_{n-r}$ 估算 $y_{n+1}$ ，即线性多步法。

据我们所查阅数值分析教材和参考书所述( [1] - [6] )，一阶常微分方程初值问题经典数值解法主要就是这些。

期刊文献上也发表了一些探索性工作。如2011年，钟巍提出的一种数值解法，先将方程变形为

$\frac{\text{d}y}{y}=\frac{f\left(x,y\right)}{y}\text{d}x$ ,

再根据右端在当前节点值的正负给出下一节点数值解。文章举例说明了算法的近似效果，但未给出收敛性证明和误差估计( [7] )。

现代人工神经网络理论和应用的发展也丰富了微分方程数值解法，通过训练带权重参数的试验解(trial solution)，使其逐渐适应方程，这便是微分方程的神经网络数值解法( [8] )。应该说，常微分方程人工神经网络数值解法具有一定开创性。

2022年，孙波、陈慧雄对一阶常微分方程初值问题解的存在唯一性Picard迭代证明法加以改进，提高了迭代收敛速度 [9] 。微分方程初值问题(1.1)可转化为等价积分方程

$y=y_0+{\displaystyle {\int }_{x_0}^{x}f\left(x,y\right)\text{d}x}$ .(1.2)

通过构造迭代函数序列

${\phi }_{n+1}\left(x\right)=y_0+{\displaystyle {\int }_{x_0}^{x}f\left(x,{\phi }_n\left(x\right)\right)\text{d}x}$ (1.3)

可推导解的存在性。一般教科书拿 ${\phi }_0\left(x\right)\equiv y_0$ 起始迭代，只要 $f\left(x,y\right)$ 在 $\left(x_0,y_0\right)$ 的某矩形邻域内连续，且关于y李普希兹连续，则 $\left\{{\phi }_n\left(x\right)\right\}$ 在 $x_0$ 的某邻域内一致收敛，其极限函数 $\phi \left(x\right)$ 便是积分方程(1.2)即初值问题(1.1)的解。Picard迭代的作用是改进局部近似解，每迭代一步使近似解更接近精确解。孙波、陈慧雄将迭代起始函数改为

${\phi }_0\left(x\right)\equiv y_0+f\left(x_0,y_0\right)\left(x-x_0\right)$ ,

并证明了所得迭代函数列也一致收敛，且比以 ${\phi }_0\left(x\right)\equiv y_0$ 起始迭代生成的函数列更快收敛到精确解。

考虑到Picard迭代具有修正局部近似解作用，而欧拉折线法实质上就是分段局部线性近似求解。若将每个节点处局部线性近似解用Picard迭代修正一次，应更接近精确解。基于这一考虑，我们用Picard迭代对欧拉折线法每个节点处局部线性近似解修正一次。第二节陈述算法，证明其收敛性与稳定性，并给出误差估计。第三节通过计算实例说明算法效果，并和欧拉折线法加以比较。

2. 基于Picard迭代修正的数值解法

本文算法可表述如下：

第一步：将初始点 $x_0$ 处一阶近似解

$y=y_0+f\left(x_0,y_0\right)\left(x-x_0\right)$

代入(1.2)或(1.3)右边，得修正的局部解

$y=y_0+{\displaystyle {\int }_{x_0}^{x}f\left(x,y_0+f\left(x_0,y_0\right)\left(x-x_0\right)\right)\text{d}x}$ .(2.1)

然后用(2.1)估算节点 $x_1=x_0+h$ 处的解值：

$y_1=y_0+{\displaystyle {\int }_{x_0}^{x_1}f\left(x,y_0+f\left(x_0,y_0\right)\left(x-x_0\right)\right)\text{d}x}$ .

第二步：对 $x_1$ 处线性近似：

$y=y_1+f\left(x_1,y_1\right)\left(x-x_1\right)$

进行Picard修正：

$y=y_1+{\displaystyle {\int }_{x_1}^{x}f\left(x,y_1+f\left(x_1,y_1\right)\left(x-x_1\right)\right)\text{d}x}$ ,

并用它估算 $x=x_2$ 处的解值。

对后续节点 $x=x_2,x_3,\cdots$ 重复第二步，节点 $x_{n+1}$ 处的数值解可表示为

$y_{n+1}=y_n+{\displaystyle {\int }_{x_n}^{x_{n+1}}f}\left(x,y_n+f\left(x_n,y_n\right)\left(x-x_n\right)\right)\text{d}x$ , $n=0,1,2,\cdots$ . (2.2)

上述算法的局部截断误差定义为

${\bar{e}}_{n+1}=y\left(x_{n+1}\right)-{\bar{y}}_{n+1}$ ,

其中 $y\left(x_{n+1}\right)$ 为精确解在 $x=x_{n+1}$ 处的值， ${\bar{y}}_{n+1}$ 表示取 $y_n=y\left(x_n\right)$ 时 $x=x_{n+1}$ 处的数值解，即

${\bar{y}}_{n+1}=y\left(x_n\right)+{\displaystyle {\int }_{x_n}^{x_{n+1}}f}\left(x,y\left(x_n\right)+f\left(x_n,y\left(x_n\right)\right)\left(x-x_n\right)\right)\text{d}x$ .(2.3)

于是

$\begin{array}{c}{\bar{e}}_{n+1}=y\left(x_{n+1}\right)-y\left(x_n\right)-{\displaystyle {\int }_{x_n}^{x_{n+1}}f}\left(x,y\left(x_n\right)+f\left(x_n,y\left(x_n\right)\right)\left(x-x_n\right)\right)\text{d}x\\ ={\displaystyle {\int }_{x_n}^{x_{n+1}}{y}^{\prime }\left(x\right)\text{d}x}-{\displaystyle {\int }_{x_n}^{x_{n+1}}f}\left(x,y\left(x_n\right)+f\left(x_n,y\left(x_n\right)\right)\left(x-x_n\right)\right)\text{d}x\\ ={\displaystyle {\int }_{x_n}^{x_{n+1}}f}\left(x,y\left(x\right)\right)\text{d}x-{\displaystyle {\int }_{x_n}^{x_{n+1}}f}\left(x,y\left(x_n\right)+f\left(x_n,y\left(x_n\right)\right)\left(x-x_n\right)\right)\text{d}x,\end{array}$ (2.4)

对(2.4)右端积分式中被积函数在 $x_n$ 进行一阶泰勒展开：

$f\left(x,y\left(x\right)\right)=f\left(x_n,y\left(x_n\right)\right)+\left[f_x\left(x_n,y\left(x_n\right)\right)+f_y\left(x_n,y\left(x_n\right)\right)f\left(x_n,y\left(x_n\right)\right)\right]\left(x-x_n\right)+o\left(x-x_n\right)$ ,

$\begin{array}{l}f\left(x,y\left(x_n\right)+f\left(x_n,y\left(x_n\right)\right)\left(x-x_n\right)\right)\\ =f\left(x_n,y\left(x_n\right)\right)+\left[f_x\left(x_n,y\left(x_n\right)\right)+f_y\left(x_n,y\left(x_n\right)\right)f\left(x_n,y\left(x_n\right)\right)\right]\left(x-x_n\right)+o\left(x-x_n\right),\end{array}$

其中 $o\left(x-x_n\right)$ 表示 $h\to 0$ 时 $x-x_n$ 的高阶无穷小。

将它们代入(2.4)式，整理得

${\bar{e}}_{n+1}={\displaystyle {\int }_{x_n}^{x_{n+1}}o\left(x-x_n\right)\text{d}x}=o\left(h^2\right)$ . (2.5)

欧拉折线法的局部截断误差为 $O\left(h^2\right)$ ，我们的改进算法比较欧拉折线法至少提高了一阶精度。

由局部截断误差可进一步估计整体截断误差

$e_{n+1}=y\left(x_{n+1}\right)-y_{n+1}$ .

由(2.2)和(2.3)可得

$\left|{\bar{y}}_{n+1}-y_{n+1}\right|\le \left|y\left(x_n\right)-y_n\right|+{\displaystyle {\int }_{x_n}^{x_{n+1}}\left|f\left(x,y\left(x_n\right)+f\left(x_n,y\left(x_n\right)\right)\left(x-x_n\right)\right)-f\left(x,y_n+f\left(x_n,y_n\right)\left(x-x_n\right)\right)\right|\text{d}x}$ . (2.6)

设 $f\left(x,y\right)$ 关于y的李普希兹常数为L，则由(2.6)进一步有

$\left|{\bar{y}}_{n+1}-y_{n+1}\right|\le \left(1+Lh+\frac{1}{2}{L}^2{h}^2\right)\left|y\left(x_n\right)-y_n\right|$ .

而局部截断误差结论(2.5)还可表述为

$\left|y\left(x_{n+1}\right)-{\bar{y}}_{n+1}\right|\le C{h}^3$

C为某正常数。

于是，

$\left|y\left(x_{n+1}\right)-y_{n+1}\right|\le \left|y\left(x_{n+1}\right)-{\bar{y}}_{n+1}\right|+\left|{\bar{y}}_{n+1}-y_{n+1}\right|\le \left(1+Lh+\frac{1}{2}{L}^2{h}^2\right)\left|y\left(x_n\right)-y_n\right|+C{h}^3$ .

即整体截断误差满足递推关系

$e_{n+1}\le \left(1+Lh+\frac{1}{2}{L}^2{h}^2\right)e_n+C{h}^3$ .

向前递推可得

$e_n\le {\left(1+Lh+\frac{1}{2}{L}^2{h}^2\right)}^n{e}_0+\frac{{\left(1+Lh+\frac{1}{2}{L}^2{h}^2\right)}^n-1}{L+\frac{1}{2}{L}^{2}h}C{h}^2$ .

因 $e_0=0$ ，进一步有

$e_n\le \frac{{\left(1+Lh+\frac{1}{2}{L}^2{h}^2\right)}^n-1}{L+\frac{1}{2}{L}^{2}h}C{h}^2$ .

由此可推知，本文算法整体截断误差为 $o\left(h\right)$ ，算法收敛性不言而喻。

从局部截断误差阶数看，本文算法达到了梯形公式和预估–校正公式的计算精度。但梯形公式和预估–校正公式为隐式，而我们的算法公式为显式。

3. 数值模拟

例3.1 考虑初值问题

$y^{\prime }=2xy$ ， $y\left(0\right)=1$ . (3.1)

用分离变量法可得其解析解 $y=\exp \left(x^2\right)$ 。为检验本文算法，在区间 $0\le x\le 1$ 内分别取离散化步长 $h=0.2$ 和0.1，分别用欧拉折线法和Picard迭代修正法进行数值计算，比较二者的计算误差。

用欧拉折线法近似求解初值问题(3.1)的公式为

$y_0=1$ , $y_{n+1}=y_n+2x_n{y}_{n}h$ ,

Picard迭代修正法近似公式为

$y_0=1$ ,

$\begin{array}{c}{y}_{n+1}=y_n+{\displaystyle {\int }_{x_n}^{x_{n+1}}2x}\left[y_n+2x_n{y}_n\left(x-x_n\right)\right]\text{d}x\\ =y_n+y_n\left(x_{n+1}^2-x_n^2\right)+\frac{4}{3}{x}_n{y}_n\left(x_{n+1}^3-x_n^3\right)-2x_n^2\left(x_{n+1}^2-x_n^2\right)\\ =y_n\left\{1+\left[\left(x_n+x_{n+1}\right)\left(1-2x_n^2\right)+\frac{4}{3}{x}_n\left(x_n^2+x_n{x}_{n+1}+x_{n+1}^2\right)\right]h\right\}.\end{array}$

先取 $h=0.2$ ，分别用欧拉折线法和Picard迭代修正法近似求解(3.1)，精确解、欧拉折线法近似解和Picard迭代近似解在各节点处的数值比较如表1。

由表中数值计算结果可以看出，Picard迭代修正法所得数值解比欧拉折线法所得数值解更接近精确解。欧拉折线法在五个节点处数值解与精确解最大相差为0.1526，Picard迭代修正法在五个节点处数值解与精确解最大相差为0.0586，约为前者三分之一。

同一坐标框内解析解曲线、欧拉折线法和Picard迭代修正近似解折线图如图1所示，图中标识符“Euler”、“precision”和“Picard”所指曲线分别为Euler折线数值解图象、精确解图象和Picard迭代修正解图象。

Table 1. Values of the exact solution, approximate solutions of (3.1) by Euler method and Picard iteration respectively, h = 0.2

表1. 初值问题(3.1)精确解、欧拉折线近似解和Picard迭代修正近似解在节点处的值， $h=0.2$

Figure 1. Graphs of the exact solution, approximate solution of (3.1) by Euler method and Picard iteration respectively, $h=0.2$

图1. 初值问题(3.1)精确解曲线，欧拉法及Picard迭代修正解折线图比较， $h=0.2$

然后取 $h=0.1$ ，分别用欧拉折线法和Picard迭代修正法近似求解(3.1)，精确解、欧拉折线法近似解和Picard迭代近似解在各节点处的数值比较如表2。

Table 2. Values of the exact solution, approximate solution of (3.1) by Euler method and Picard iteration respectively, h = 0.1

表2. 初值问题(3.1)精确解、欧拉折线法近似解和Picard迭代修正法近似解在节点处的值， $h=0.1$

取 $h=0.1$ 所得欧拉折线法近似解与精确解在各节点处相差最大值为0.086725，Picard迭代修正法近似解与精确解在各节点相差最大值为0.016316，不到前者五分之一。欧拉折线数值解最大误差比 $h=0.2$ 时缩小将近一半，Picard迭代修正数值解最大误差不到 $h=0.2$ 时的三分之一，体现出Picard迭代修正数值解比欧拉折线数值解收敛性也更好。

Figure 2. Graphs of the exact solution, approximate solution of (3.1) by Euler method and Picard iteration respectively, $h=0.1$

图2. 初值问题(3.1)精确解曲线、欧拉折线法及Picard迭代修正解折线图比较， $h=0.1$

同一坐标框内解析解曲线、欧拉折线数值解和Picard迭代修正解折线图如图2。

参考文献