基于常微分方程中迭代数值解法的剖析——以Picard方法为例①

刘禹希， 王 密

(武汉大学a.数学与统计学院;b.遥感信息工程学院，湖北 武汉 430061)

1 Cauchy问题的一种数值解法

Cauchy问题的求解是常微分方程理论中十分重要的课题，也是一阶常微分方程初值问题中最基础的一环。关于这种初值问题解的存在性理论，有以下著名的Picard定理：

定理1(Picard定理)

设初值问题：

(1)

式(1)中f(x,y)在矩形区域

R:|x-x0|≤a,|y-y0|≤b

(2)

内连续，且对y满足Lipshitz条件，则(E)在区间I=[x0-h,x0+h]上有且仅有一个连续可微解，其中常数

本文仅研究形如(E)，满足条件和Lipshitz条件的常微分方程初值问题的数值解。

在该定理证明过程中将微分方程通过代数变形转化为积分方程，并构造了如下Picard序列：

(3)

式(3)中y(x0)=y0

事实上通过选定初值并对该序列进行迭代，经过多次操作便可得到该微分方程的一个数值解，其截断误差有估计:

(4)

将以上数值求解方法称为Picard方法。

定理 1表明：对一类特殊的Cauchy问题,其解一定可以唯一地定义在一个足够小的闭区间中。但是对于整体积分曲线的研究，仍需要研究其局部解的最大存在区间问题，于是便有如下的微分方程解的延拓定理。

以下不加证明地给出该定理：

定理2(解的延拓定理)

设p为区域G内一点，设T为微分方程经过p的任一积分曲线。则积分曲线T将在区域G内延伸到边界。

2 数值积分方法

解决简单数值积分的常见方法有：利用牛顿—柯斯特公式求解，利用复合求积法得到复合梯形公式和复合辛普森公式求解。前面介绍的求积方法可提高求积精度，如若精度仍不够，则可通过将步长逐次减半的方式来提高精度,将改进后的公式称之为龙贝格求积公式。

贯穿于数值积分中最经典的思想便是“分割求和”，于是有以下著名的复合梯形公式：

公式1(复合梯形公式)

将[a,b]分为n等份，其节点为xi=a+ih,i=0,1,2,…,n,h=b-a/n，则有：

(5)

如对复合梯形公式进行代数变形即可导出其递推公式(6):

(6)

式(6)中T2n表示在Tn基础上步长h减半后的积分值。

为了便于刻画复合梯形公式中的积分项，可以采用幂级数展开的方法逼近计算，于是便有著名的理查德森定理。

以下不加证明地给出该定理：

定理3(理查德森定理)

设f(x)∈C∞[a,b]，则有：

T(h)=I+α1h2+α2h4+…+αlh2l+…

(7)

其中T(h)=Tn,系数αl(l=1,2,…)与h无关。

从式(7)可以看出，用近似积分值I，误差阶为O(h4)，误差阶得到了提高。类似地，将计算I的近似值的误差阶由O(h2)提高到O(h4)的方法称为外推算法，也称为理查德森外推算法[3]。

以上是数值分析中一个重要的技巧，当真值与近似值的误差能够表示成h的幂级数时，都可以使用外推算法来提高精度。以下将充分利用这种方法，进行适当的改进，推导出高精度的龙贝格求积算法。

3 龙贝格求积算法

龙贝格算法是在积分区间逐次分半的过程中，对用复化梯形法产生的近似值进行加权平均，以获得准确度较高的一种方法。具有公式简练，使用方便，结果更可靠的特点。

事实上，对于I的计算可以转化为对于近似迭代序列{Tn}的计算，通过迭代序列各项的比对，经过代数变形整理即可。

龙贝格算法推导如下：

设m次外推加速为式(8)：

(8)

余项为式(9)：

Tm(h)=I+δ1h2(m+1)+δ2h2(m+2)+…

(9)

(10)

式(10)称为龙贝格求积算法[4]，算法运行过程如下：

3.求加速值。

4 算法应用实例

下面以一个此类型初值问题的数值求解来展示这种方法的实际应用和重要意义。

4.1 例子

(11)

这是一类标准的一阶线性的微分方程初值问题。由定理1，可知该初值问题存在唯一解；由定理2，可知该解可以延拓至整个区间[0,1]上成为一条积分曲线。

同时以上常系数线性微分方程可以通过求解其特征方程得到特解，并用定系数求得通解的方法得到其解析解为：y=ex+2x-x2.

4.2 Picard迭代算法

根据以上定理以及算法的铺垫，整理可得常微分方程数值求解中的新型Picard算法，流程如下：

(1)写出Picard迭代序列：

(3)将以上数值积分项代入Picard序列并计算出yn+1(x)；

(4)交替进行2和3直至迭代序列满足|yn+1(x)-yn(x)|≤0.01，停止迭代。

4.3 结果展示

以下将通过比对常微分方程数值求解中著名的梯形迭代格式算法来体现Picard迭代算法的优越性[5]。

利用MATLAB对以上Picard算法进行迭代，可得如下结果：

由图1可知：该序列仅仅迭代四次便可很好地收敛到解析解，收敛速度非常快，这体现了Picard迭代算法的高效性。

图1 迭代序列

以下是误差表格1。

表1 Picard迭代误差表

将三种迭代格式的结果用MATLAB作图呈现：

由图2可知，相比于常微分方程数值求解中的梯形迭代格式，本文构建的Picard迭代序列有着更好的收敛精度。

图2 算法比较

分别取x=0.2,0.4,0.6,0.8,1.0，将Picard序列和梯形迭代序列与标准解析解的误差做成误差表如下表2。

表2 不同数值求解方法误差表

由表2可知，当变量x增大时，Picard迭代误差及梯形迭代误差均增加，但前者增加的速度明显慢于后者。因此相比于传统数值迭代方法，新型Picard迭代方法的效率明显较高。

事实上，在面对一类特殊的常微分方程初值问题时，相比于传统迭代方法，往往用Picard迭代方法可以更快地收敛到解，这也是迭代不动点方法的重要作用之一。

5 结 语

对传统Picard迭代算法进行改进，引入复合龙贝格数值积分方法计算迭代序列，将求解精度进一步优化，并结合实例验证其迭代稳定性同时和传统数值方法进行对比，得到了很好的结果；同时，该算法也体现出迭代不动点数值方法极快的收敛速度，在面对更复杂的多项式型Cauchy问题时，该方法的求解所占内存相对较小，有着很好的应用前景。但是在面对一系列略病态的Cauchy问题时，该方法求龙贝格数值积分部分的运算复杂度仍然会很大，如何针对性地寻找更便捷的方法也是未来常微分方程数值求解方法应当着重研究的一个领域。