常微分方程初值问题数值解法

常微分方程初值问题数值解法9.1引言9.2简单的数值方法与基本概念9.3龙格－库塔方法9.4单步法的收敛性与稳定性9.5线性多步法9.6方程组和高阶方程9.1

引言

本章讨论一阶常微分方程的初值问题：只要函数适当光滑—如满足利普希茨条件：理论上就能保证初值问题的解存在并且唯一。所谓数值解法，就是寻求解在一系列离散点上的近似值，相邻两个点间的距离称为步长。一般情况下我们取为常数，这是节点为：9.1

引言

初值问题的求解有一个基本特点，它们都是采取“步进式”求解的，即，一步一步地求函数的值。求解的主要方法是：先对方程进行离散化，建立求数值解的递推公式，一类在计算时只用到前面一步的，称为单步法。另一类在计算时除了用到理论上就能保证初值问题的解存在并且唯一。

所谓数值解法，就是寻求解在一系列离散点上的近似值，相邻两个点间的距离称为步长。一般情况下我们取为常数，这是节点为：

初值问题的求解有一个基本特点，它们都是采取“步进式”求解的，即，一步一步地求函数的值。求解的主要方法是：先对方程进行离散化，建立求数值解的递推公式，一类在计算时只用到前面一步的，称为单步法。另一类在计算时除了用到前利用前面一步的，还要用到前面的等，这种方称为多步法。其次，要研究公式的局部截断误差和阶。数值解和精确解的误差估计和收敛性，还有递推迭代公式的数值稳定性问题。显式欧拉法在平面上，微分方程初值问题的解称作方程的积分曲线。积分曲线上每一点的斜率等于该点的函数值。如果按函数在平面上建立一个方向场，那么，积分曲线上每一点的切线方向均与方向场在该点的方向一致。方称为多步法。其次，要研究公式的局部截断误差和阶。数值解和精确解的误差估计和收敛性，还有递推迭代公式的数值稳定性问题。9.2

简单的数值方法与基本概念显式欧拉法在平面上，微分方程初值问题的解称作方程的积分曲线。积分曲线上每一点的斜率等于该点的函数值。如果按函数在平面上建立一个方向场，那么，积分曲线上每一点的切线方向均与方向场在该点的方向一致。基于上面的几何解释，我们从初始点出发，先依方向场在该点的方向推进到上一点，然后再从依方向场的方向推进到上一点，循此前进做出一条折线。9.2

简单的数值方法与基本概念

一般地，设已做出该折线上的顶点，过依方向场的方向再推进到，显然两个顶点的坐标有如下关系即

基于上面的几何解释，我们从初始点出发，先依方向场在该点的方向推进到上一点，然后再从依方向场的方向推进到上一点，循此前进做出一条折线。

一般地，设已做出该折线上的顶点，过依方向场的方向再推进到，显然两个顶点的坐标有如下关系即这就是著名的欧拉公式。若初值已知，则依该公式可逐步算出：

例1求解初值问题：(其解为)解：根据欧拉方法，得到：这就是著名的欧拉公式。若初值已知，则依该公式可逐步算出：

例1求解初值问题：(其解为)解：根据欧拉方法，得到：取步长，计算得到如下结果：0.11.10001.09540.61.50901.48320.21.19181.18320.71.58031.54920.31.27741.26490.81.64981.61250.41.35821.34160.91.71781.67330.51.43151.41421.01.78481.7321隐式欧拉法在前面的讨论中，近似计算公式可以看成是由在区间上积分得到，而右边的积分是利用左矩形公式近似，再以代替得到，现在右端的积分用右矩形公式，则得到：取步长，计算得到如下结果：0.11.10001.09540.61.50901.48320.21.19181.18320.71.58031.54920.31.27741.26490.81.64981.61250.41.35821.34160.91.71781.67330.51.43151.41421.01.78481.7321隐式欧拉法在前面的讨论中，近似计算公式可以看成是由在区间上积分得到，而右边的积分是利用左矩形公式近似，再以代替得到，现在右端的积分用右矩形公式，则得到：此式的右端包含，是一种隐式的单步法，称为隐式欧拉方法。利用此法，每一步都要把该式作为的方程来求解。从数值积分的误差来分析，很难期望隐式欧拉方法比显式欧拉方法精确。为了得到更精确的方法，我们用梯形公式来近似前面的积分，得到梯形方法：它也是一种隐式单步法。改进的欧拉法从梯形法的推导，可望它比欧拉法更精确，但它计算量较大，在实际计算中，可取欧拉法的结果为迭代计算的初值，然后再用梯形公式计算一次，得到：欧拉方法。利用此法，每一步都要把该式作为的方程来求解。从数值积分的误差来分析，很难期望隐式欧拉方法比显式欧拉方法精确。为了得到更精确的方法，我们用梯形公式来近似前面的积分，得到梯形方法：它也是一种隐式单步法。改进的欧拉法从梯形法的推导，可望它比欧拉法更精确，但它计算量较大，在实际计算中，可取欧拉法的结果为迭代计算的初值，然后再用梯形公式计算一次，得到：或写成：改进的欧拉法是一种显式单步法，有时为了计算方便，也用下面的公式：其中：单步法的截断误差与阶初值问题的单步法可用一般形式表示为：其中多元函数与有关，当含有时，方法是隐式的，不含时则为显式的，所以显式单步法可表示为：称为增量函数，对于欧拉法，或写成：改进的欧拉法是一种显式单步法，有时为了计算方便，也用下面的公式：其中：单步法的截断误差与阶初值问题的单步法可用一般形式表示为：其中多元函数与有关，当含有时，方法是隐式的，不含时则为显式的，所以显式单步法可表示为：称为增量函数，对于欧拉法，为了分析其截断误差，我们将在处作泰勒展开，得到：假设是精确的，于是可得欧拉法计算公式的误差为：

定义1设是初值问题的准确解，称为显式单步法的局部截断误差。之所以称为局部的，是假设在前各步没有误差，当时，计算一步，则有开，得到：假设是精确的，于是可得欧拉法计算公式的误差为：

定义1设是初值问题的准确解，称为显式单步法的局部截断误差。之所以称为局部的，是假设在前各步没有误差，当时，计算一步，则有所以，局部截断误差可以理解为用计算的方法计算一步的误差，这样，欧拉法的局部截断误差为：这里称为局部截断误差主项，显然一般情况下的定义如下：

定义2设是初值问题的准确解，若存在最大整数使显式单步法的局部误差满足则称该方法具有阶精度。若将上式展开成则称为局部截断误差主项。所以，局部截断误差可以理解为用计算的方法计算一步的误差，这样，欧拉法的局部截断误差为：这里称为局部截断误差主项，显然一般情况下的定义如下：

定义2设是初值问题的准确解，若存在最大整数使显式单步法的局部误差满足则称该方法具有阶精度。若将上式展开成则称为局部截断误差主项。以上定义对隐式单步法也适用，如对后退欧拉法，其局部截断误差为：这里，是1阶方法，局部截断误差主项为同样，对梯形方法有：所以梯形方法是2阶的，其局部误差主项为：

以上定义对隐式单步法也适用，如对后退欧拉法，其局部截断误差为：这里，是1阶方法，局部截断误差主项为

同样，对梯形方法有：所以梯形方法是2阶的，其局部误差主项为：改进欧拉法的例梯形法精度得到提高，但算法复杂，计算量大。为了控制计算量，通常只迭代一两次就转入下一步的计算，这样就可以简化计算。具体做法：先用欧拉公式求一个初值，称为预测值，再用梯形公式校正一次，这样建立的公式通常称为改进的欧拉公式：预测：校正：或表示为：所以梯形方法是2阶的，其局部误差主项为：改进欧拉法的例梯形法精度得到提高，但算法复杂，计算量大。为了控制计算量，通常只迭代一两次就转入下一步的计算，这样就可以简化计算。具体做法：先用欧拉公式求一个初值，称为预测值，再用梯形公式校正一次，这样建立的公式通常称为改进的欧拉公式：预测：校正：或表示为：例：用改进的欧拉方法求解初值问题：解：仍取，计算结果如表0.11.09591.09540.61.48601.48320.21.18411.18320.71.55251.54920.31.26621.26490.81.61531.61250.41.34341.34160.91.67821.67330.51.41641.41421.01.73791.7321显式龙格－库塔法的一般形式欧拉法的局部截断误差为，是1阶的方法对于改进的欧拉法，它可表示为：此时增量函数例：用改进的欧拉方法求解初值问题：解：仍取，计算结果如表0.11.09591.09540.61.48601.48320.21.18411.18320.71.55251.54920.31.26621.26490.81.61531.61250.41.34341.34160.91.67821.67330.51.41641.41421.01.73791.73219.3

龙格－库塔方法显式龙格－库塔法的一般形式欧拉法的局部截断误差为，是1阶的方法对于改进的欧拉法，它可表示为：此时增量函数它比欧拉法的多计算了一个函数值，可望，若要使得到的公式阶数更大，就必须包含更多的值，实际上从方程等价的积分形式，即

(*)9.3

龙格－库塔方法要使公式的阶数提高，就必须使右端积分的数值求积公式精度提高，它必然要增加求积节点，为此可将(*)的右端用求积公式表示为：一般来说，点数越多，精度越高，上式右端相当于增量函数，为得到便于计算的显式方法，可它比欧拉法的多计算了一个函数值，可望，若要使得到的公式阶数更大，就必须包含更多的值，实际上从方程等价的积分形式，即

(*)要使公式的阶数提高，就必须使右端积分的数值求积公式精度提高，它必然要增加求积节点，为此可将(*)的右端用求积公式表示为：一般来说，点数越多，精度越高，上式右端相当于增量函数，为得到便于计算的显式方法，可类似于改进欧拉法，将公式表示为：其中这里均为常数，方法称为级显式R-K方法

当时，就是欧拉法，此时当时，改进欧拉法就是其中的一种，下面证明，要使公式有更高的阶，就要增加点数。下面我们仅就推导R-K方法，并给出时的常用公式，其推导方法与时类似，只是计算较为复杂。增量函数，为得到便于计算的显式方法，可类似于改进欧拉法，将公式表示为：其中这里均为常数，方法称为级显式R-K方法

当时，就是欧拉法，此时当时，改进欧拉法就是其中的一种，下面证明，要使公式有更高的阶，就要增加点数。下面我们仅就推导R-K方法，并给出时的常用公式，其推导方法与时类似，只是计算较为复杂。二阶显式R-K方法对于的R-K方法，可以得到如下的计算公式：这里均为待定常数，我们希望适当选取这些系数，使公式阶数尽量高。根据局部截断误差的定义，上式的局部截断误差为：

这里为得到的阶，要将上式各项在处做泰勒展开，由于是二元函数，故要用到二元泰勒展开，各项展开为：二阶显式R-K方法对于的R-K方法，可以得到如下的计算公式：这里均为待定常数，我们希望适当选取这些系数，使公式阶数尽量高。根据局部截断误差的定义，上式的局部截断误差为：这里为得到的阶，要将上式各项在处做泰勒展开，由于是二元函数，故要用到二元泰勒展开，各项展开为：其中将以上结果代入：得到：其中将以上结果代入：得到：要使公式成为2阶的，必须有：即：上式的解不唯一，可以令，则得这样得到的公式称为二阶R-K方法，如取，则，这就是改进的欧拉法。若取，则，得到：要使公式成为2阶的，必须有：即：上式的解不唯一，可以令，则得这样得到的公式称为二阶R-K方法，如取，则，这就是改进的欧拉法。若取，则，得到：称为中点公式，相当于数值积分的中矩形公式，也即

对的R-K公式能否使局部误差提高到？为此需要把多展开一项，通过分析可以得到不可能的结论。因此时只能得到2阶的公式。三阶与四阶显式R-K方法要得到三阶的显式R-K方法，必须，此时计算若取，则，得到：称为中点公式，相当于数值积分的中矩形公式，也即

对的R-K公式能否使局部误差提高到？为此需要把多展开一项，通过分析可以得到不可能的结论。因此时只能得到2阶的公式。三阶与四阶显式R-K方法要得到三阶的显式R-K方法，必须，此时计算公式为：其中及均为待定参数，误差为只要将按二元函数的泰勒公式展开，使可得待定参数需满足的方程为：

要得到三阶的显式R-K方法，必须，此时计算公式为：其中及均为待定参数，误差为只要将按二元函数的泰勒公式展开，使可得待定参数需满足的方程为：这里为8个未知数，6个方程，解也是不唯一的，可以得到很多公式，下面是其中称为库塔的公式这里为8个未知数，6个方程，解也是不唯一的，可以得到很多公式，下面是其中称为库塔的公式继续上述过程，经过较复杂的演算，可以导出各种四阶龙格－库塔公式，下面是常用的经典公式：例3：用四阶龙格－库塔方法求解定解问题：取步长，从到。解：写出具体的四阶龙格－库塔的具体格式为：继续上述过程，经过较复杂的演算，可以导出各种四阶龙格－库塔公式，下面是常用的经典公式：例3：用四阶龙格－库塔方法求解定解问题：取步长，从到。解：写出具体的四阶龙格－库塔的具体格式为：

列出计算结果如下表：0.21.18321.18320.41.34171.34160.61.48331.48320.81.61251.61251.01.73211.7321列出计算结果如下表：可见计算精度较以前大为提高。变步长的龙格－库塔方法单从每一步看，步长越小，截断误差就越小，但随着步长的缩小，在求解区间内的步数就增加了，这在增加计算量的同时也增加了计算的舍入误差，因此也有一个选择步长的问题。0.21.18321.18320.41.34171.34160.61.48331.48320.81.61251.61251.01.73211.7321选择步长需要考虑的两个问题：怎样衡量和检验计算结果的精度？如何依据所获得的精度处理步长？我们考虑经典的四阶龙格－库塔公式，从节点出发，先以为步长求出一个近似值，记为，由于公式的局部截断误差为，故变步长的龙格－库塔方法单从每一步看，步长越小，截断误差就越小，但随着步长的缩小，在求解区间内的步数就增加了，这在增加计算量的同时也增加了计算的舍入误差，因此也有一个选择步长的问题。选择步长需要考虑的两个问题：怎样衡量和检验计算结果的精度？如何依据所获得的精度处理步长？我们考虑经典的四阶龙格－库塔公式，从节点出发，先以为步长求出一个近似值，记为，由于公式的局部截断误差为，故然后将步长折半，即取为步长，从跨两步到再来求得一个近似值，每一步的截断误差是因此有比较前后两式，我们得到，步长折半后，误差大约减少到，即有由此易得下列事后估计式：这样，我们可以通过检查步长，折半前后两次计算结果的偏差来判定所选的步长是否合适，公式的局部截断误差为，故然后将步长折半，即取为步长，从跨两步到再来求得一个近似值，每一步的截断误差是因此有

比较前后两式，我们得到，步长折半后，误差大约减少到，即有由此易得下列事后估计式：这样，我们可以通过检查步长，折半前后两次计算结果的偏差来判定所选的步长是否合适，分下列两种情况处理，称为变步长方法：

1)对于给定的精度，如果，我们反复将步长折半进行计算，直到为止，这时取最终得到的作为结果；

2)如果，我们将反复将步长加倍，直到为止，这时再将步长折半一次，得到结果。9.4

单步法的收敛性与稳定性收敛性与相容性数值解法的基本思想是通过离散将微分方程转化为差分方程，如单步法

(1)它在处的解为，而初值问题的解为，记分下列两种情况处理，称为变步长方法：

1)对于给定的精度，如果，我们反复将步长折半进行计算，直到为止，这时取最终得到的作为结果；

2)如果，我们将反复将步长加倍，直到为止，这时再将步长折半一次，得到结果。9.4

单步法的收敛性与稳定性收敛性与相容性数值解法的基本思想是通过离散将微分方程转化为差分方程，如单步法

(1)它在处的解为，而初值问题的解为，记，称为整体截断误差，收敛性就是讨论当固定，且时的问题。定义3若一种数值方法对于固定的，当时有，则称该方法是收敛的。显然，数值方法收敛是指，对单步法(1)有下述收敛性定理：9.4

单步法的收敛性与稳定性

定理1假设单步法(1)具有阶精度，且增量函数关于满足利普希茨条件

(2)又设初值准确，即，则其整体截断误差为

(3)证明：设以表示取用公式(1)求得的结果,

固定，且时的问题。定义3若一种数值方法对于固定的，当时有，则称该方法是收敛的。显然，数值方法收敛是指，对单步法(1)有下述收敛性定理：

定理1假设单步法(1)具有阶精度，且增量函数关于满足利普希茨条件

(2)又设初值准确，即，则其整体截断误差为

(3)证明：设以表示取用公式(1)求得的结果,即(4)则为局部截断误差，由于所给的方法具有阶精度，按定义2，存在定数，使又由(4)和(1)，有利用条件(2)，有：从而有即对整体截断误差成立下列递推关系：反复递推，可得：证明：设以表示取用公式(1)求得的结果,即(4)则为局部截断误差，由于所给的方法具有阶精度，按定义2，存在定数，使又由(4)和(1)，有利用条件(2)，有：从而有即对整体截断误差成立下列递推关系：反复递推，可得：再注意到当时，，最终得到下列估计式：由此可以断定，如果初值是准确的，则(3)成立。

推论：对于一个阶的显式单步法，若微分方程的右端函数关于满足利普希茨条件，且初值是精确的，则显式欧拉法，改进欧拉法和龙格－库塔方法是收敛的。定理1表明，时单步法收敛，并且当是初值问题的解，且具有阶精度时，有展开式：再注意到当时，，最终得到下列估计式：由此可以断定，如果初值是准确的，则(3)成立。

推论：对于一个阶的显式单步法，若微分方程的右端函数关于满足利普希茨条件，且初值是精确的，则显式欧拉法，改进欧拉法和龙格－库塔方法是收敛的。定理1表明，时单步法收敛，并且当是初值问题的解，且具有阶精度时，有展开式：所以的充要条件是，而，于是可给出如下的定义：定义4单步法的增量函数称为该单步法与初值问题相容。

以上讨论表明阶方法当与初值问题相容，反之，相容的方法至少是一阶的。于是，由定理1可知，线性单步方法收敛的充分必要条件是该方法是相容的。绝对稳定性与绝对稳定域

所以的充要条件是，而，于是可给出如下的定义：定义4单步法的增量函数称为该单步法与初值问题相容。

以上讨论表明阶方法当与初值问题相容，反之，相容的方法至少是一阶的。于是，由定理1可知，线性单步方法收敛的充分必要条件是该方法是相容的。绝对稳定性与绝对稳定域定义5若一种数值方法在节点值上有大小为的扰动，而于以后各节点上产生的偏差均不超过，则称该方法是稳定的。例4考察初值问题其准确解是一个按指数曲线衰减的很快的函数，如图：绝对稳定性与绝对稳定域

定义5若一种数值方法在节点值上有大小为的扰动，而于以后各节点上产生的偏差均不超过，则称该方法是稳定的。例4考察初值问题其准确解是一个按指数曲线衰减的很快的函数，如图：

用欧拉方法解该方程，得到：若取，则欧拉法的具体形式为：若取，则欧拉法的具体形式为：显然，前一形式计算不稳定，后一形式计算稳定。再考察后退欧拉方法，可以得到，当时，计算是稳定的。具体数据见书上356页表9－4。例题表明：稳定性不但与方法有关，也与步长有关。当然也与有关，为了只考察数值方法本身，通常只检验将数值方法用于解模型方程的稳定性，模型方程为：，其中为复数。

对于一般的方程，将方程中的在解域内某点作泰勒展开并局部线性化：忽略高阶项，令和对上式显然，前一形式计算不稳定，后一形式计算稳定。再考察后退欧拉方法，可以得到，当时，计算是稳定的。具体数据见书上356页表9－4。例题表明：稳定性不但与方法有关，也与步长有关。当然也与有关，为了只考察数值方法本身，通常只检验将数值方法用于解模型方程的稳定性，模型方程为：，其中为复数。

对于一般的方程，将方程中的在解域内某点作泰勒展开并局部线性化：忽略高阶项，令和对上式作变量代换，得到：

这就是模型方程。首先研究欧拉方程的稳定性。模型方程的欧拉公式为：假设在节点上有一扰动，它的传播使节点值产生大小为的扰动值，假设用按欧拉公式得出的计算过程不再有新的误差，则扰动值满足可见扰动值满足原来的差分方程，这样如果差分方程的解是不增长的，即有：

这就是模型方程。首先研究欧拉方程的稳定性。模型方程的欧拉公式为：假设在节点上有一扰动，它的传播使节点值产生大小为的扰动值，假设用按欧拉公式得出的计算过程不再有新的误差，则扰动值满足可见扰动值满足原来的差分方程，这样如果差分方程的解是不增长的，即有：则它就是稳定的。这一结论对于以下研究的其它方法也适用。显然，为了保证，只要，这在的复平面上是一个以(-1,0)为圆心，1为半径的单位圆我们称其为欧拉法的绝对稳定域，一般情形可有如下的定义。

定义6单步法用于模型，若得到的解，