1、第八章常微分方程数值解法8.5.1 一阶方程组和高阶方程一阶方程组和高阶方程 考虑一阶常微分方程组的初值问题：考虑一阶常微分方程组的初值问题：，TNyyyyxy0020100 那么，（那么，（8.5.1）式可以写成）式可以写成 。，00yxyyxfyi （8.5.2）若将其中的为知函数、方程的右端项都表示成向量形式：若将其中的为知函数、方程的右端项都表示成向量形式：，TNTNffffyyyy2121 初始条件表示成初始条件表示成 。，NiyxyyyyxfyiiNii2100021 （8.5.1）8.5 一阶方程组的数值解法一阶方程组的数值解法第八章常微分方程数值解法可见，式（可见，式（8.5.
2、2）在形式上与一个方程的初值问题一样。关于一个方程的初值）在形式上与一个方程的初值问题一样。关于一个方程的初值为体的数值方法均适用于方程组。相应的理论问题也可类似地讨论。下面仅写为体的数值方法均适用于方程组。相应的理论问题也可类似地讨论。下面仅写出两种数值方法作说明。出两种数值方法作说明。梯形方法：梯形方法：，1112 nnnnnnyxfyxfhyy或表达为或表达为 ，Niyxfyxfhyynninniniin212111 其中其中 是第是第 个因变量个因变量 在节点在节点 处的近似值，相应地，处的近似值，相应地，。niyi xyinx nniyxf，nNnnniyyyxf，21经典经典R-K
3、方法：方法：，43211226KKKKhyynn 第八章常微分方程数值解法其中其中 ，34231212222hKyhxfKKhyhxfKKhyhxfKyxfKnnnnnnnn 或表达为或表达为，NiKKKKhyyiiiiniin2122643211其中其中 。，NnNnnniiNnNnnniiNnNnnniinNnnniihKyhKyhKyhxfKKhyKhyKhyhxfKKhyKhyKhyhxfKyyyxfK33223114222221131122111221122222222 第八章常微分方程数值解法对于高阶方程，可把它转化为一阶方程组。例如，考察下列对于高阶方程，可把它转化为一阶方程组。
4、例如，考察下列 m 阶微分方程。阶微分方程。，110,001mkyxyyyyxfykkmm（8.5.3）只要引进新的变量只要引进新的变量，1321 mmyyyyyyyy则可将则可将 m 阶方程（阶方程（8.5.3）化为如下的一阶方程组）化为如下的一阶方程组 。，mkyyyyyxfymiyykxkmmii21121102110 （8.5.4）因此，可用求解方程组形式的方法来求解（因此，可用求解方程组形式的方法来求解（8.5.4）。）。第八章常微分方程数值解法 先考虑两个简单的初值问题。先考虑两个简单的初值问题。问题问题1：问题问题2：。，3200sinxcosx2sinx22112vuvuvu
5、。，3200sinxcosx999sinx299999812vuvuvu这两个问题有同样的解这两个问题有同样的解 。xxexvxuxcossin1128.5.2 刚性方程组刚性方程组第八章常微分方程数值解法 采用四阶经典采用四阶经典R-K方法来计算上面的两个问题，以相同的误差要求来自动选方法来计算上面的两个问题，以相同的误差要求来自动选取步长，计算从取步长，计算从 x=0 到到 x=10。第一个问题可用相当大的步长，而第二个问题。第一个问题可用相当大的步长，而第二个问题能能使用的步长小到难以接受。如果改用某种低阶隐式公式，那么这两个问题均可用使用的步长小到难以接受。如果改用某种低阶隐式公式，那
6、么这两个问题均可用较大的步长，计算出大致符合要求的解来。上述显示出来的现象称刚性。问题较大的步长，计算出大致符合要求的解来。上述显示出来的现象称刚性。问题2是刚性的，问题是刚性的，问题1是非刚性的。由于这两个问题的解是相同的，因此这种现象不是非刚性的。由于这两个问题的解是相同的，因此这种现象不是问题的解的作用，而是方程组的一种特性所引起的。基于这个事实，较为正确是问题的解的作用，而是方程组的一种特性所引起的。基于这个事实，较为正确的应称之为刚性方程组而不是刚性问题。的应称之为刚性方程组而不是刚性问题。考虑方程组的通解。对于问题考虑方程组的通解。对于问题1，方程组的系数矩阵特征值为，方程组的系数
7、矩阵特征值为 和和 ，其通解为，其通解为1132 ，xxeexvxuxxcossin1111321 第八章常微分方程数值解法其中其中 为任意常数为任意常数21，xxeexvxuxxcossin998111100021 其中其中 为任意常数。对于问题为任意常数。对于问题2，方程组的系数矩阵的特征值为，方程组的系数矩阵的特征值为 和和 ，其通解为，其通解为111000221，数值计算中出现的现象可以用稳定性来解释。两个问题的特征值都是实的，因数值计算中出现的现象可以用稳定性来解释。两个问题的特征值都是实的，因此可只考虑稳定区间。经典此可只考虑稳定区间。经典 R-K 方法的绝对稳定区间近似为（方法的
8、绝对稳定区间近似为（-2.785，0）。）。对于问题对于问题1，如果，如果 或或 时，可以是稳定的。对于问时，可以是稳定的。对于问题题2，要求，要求 或或 ，才能保证稳定。由上，才能保证稳定。由上可以看出，一定精度范围内，可以看出，一定精度范围内，h 完全由绝对稳定性决定。完全由绝对稳定性决定。0785.23，h928.0h0785.21000，h002785.0h第八章常微分方程数值解法 由通解（由通解（8.5.6）可见，当）可见，当 时，（时，（8.5.6）右边的第一项和第二项都趋）右边的第一项和第二项都趋于零，这两项瞬态解。趋于零的快慢取决于特征值的大小。显然，第二项很快趋于零，这两项瞬
9、态解。趋于零的快慢取决于特征值的大小。显然，第二项很快趋于零，此项称为快瞬态解，而第一项称为慢瞬态解。（于零，此项称为快瞬态解，而第一项称为慢瞬态解。（8.5.6）右边的第三项称为）右边的第三项称为稳态解。实际计算表明稳态解。实际计算表明 ，当第二项很快趋于零以后，当第二项很快趋于零以后，要使整个方程组的解趋于，要使整个方程组的解趋于稳稳态解，必须由第一项来决定计算终止与否。因此计算中，要在一个很长的区间上态解，必须由第一项来决定计算终止与否。因此计算中，要在一个很长的区间上处处用小步长来计算，这就是处处用小步长来计算，这就是刚性现象刚性现象。由（由（8.5.5）和（）和（8.5.6）可见，计
10、算的）可见，计算的步步数与量数与量 有关。有关。一般地，考虑非齐次常系数方程组一般地，考虑非齐次常系数方程组 x12/xAyy ，（8.5.7）其中，其中，为常系数矩阵。设为常系数矩阵。设 A 有不同的特征值，有不同的特征值，它们对应的特征向量它们对应的特征向量 ，则方程组（，则方程组（8.5.7）的通解为）的通解为 mmmRARy，mkCk，21mkCumk，21 ，mkxkxexyk1 第八章常微分方程数值解法其中其中 为任意常数，为任意常数，为（为（8.5.7）的特解。）的特解。假定特征值假定特征值 的实部为负，即的实部为负，即mkk，21k。，mkk210Re现设现设 A 的特征值按其
11、实部的绝对值大小排列：的特征值按其实部的绝对值大小排列：mReReRe21当当 时，时，kmkkxkuek1趋向于零，此项称为趋向于零，此项称为瞬态解瞬态解，而，而 则称为则称为稳态解稳态解。如果。如果 大，那么对应大，那么对应项项 当当 x 增加快速衰减，此项称为增加快速衰减，此项称为快瞬态解快瞬态解。如果。如果 小，那么对小，那么对应的项应的项 当当 x 增加时衰减慢，称其为增加时衰减慢，称其为慢瞬态解慢瞬态解。kRekxkuekkRekxkuek第八章常微分方程数值解法当我们计算稳态解时，必须求到当我们计算稳态解时，必须求到 可以忽略为止，所以可以忽略为止，所以 越越小，计算的区间越长。
12、另一方面，为使小，计算的区间越长。另一方面，为使 均在绝对稳定性均在绝对稳定性区域内，显然，区域内，显然，的值大时，必须采用很小的步长的值大时，必须采用很小的步长 h。因此，引入微分方。因此，引入微分方程程组（组（8.5.7）的）的刚性比刚性比111uex 1Re mkhk，21 m Re ，1ReRe mS （8.5.8）那么，我们似乎可以用刚性比来描述刚性方程组，那么，我们似乎可以用刚性比来描述刚性方程组，即方程组（即方程组（8.5.7）中）中 A 的的全部特征值有负的实部并刚性比全部特征值有负的实部并刚性比 S 是大的，那么（是大的，那么（8.5.7）是刚性的。）是刚性的。上述描述性定义
13、有时也会发生一些不妥，比如，该定义为能包含实际问上述描述性定义有时也会发生一些不妥，比如，该定义为能包含实际问题中常常出现的特征值实部为小的正数或等于零的情况。因此，我们引入下题中常常出现的特征值实部为小的正数或等于零的情况。因此，我们引入下面的定义。面的定义。第八章常微分方程数值解法 定义定义8.6 当具有有限的绝对稳定区域的数值方法应用到一个任意初始条件当具有有限的绝对稳定区域的数值方法应用到一个任意初始条件的的方程组时，如果在求解区间上必须用非常小的步长，则称此方程组在该区间上是方程组时，如果在求解区间上必须用非常小的步长，则称此方程组在该区间上是刚性的刚性的。刚性方程组有其自身的特点，一般显式方法难于应用。梯形方法、隐刚性方程组有其自身的特点，一般显式方法难于应用。梯形方法、隐 式式Euler法对法对 h不不 限制，可适用于一定类型的刚性方程组的求解。这里，我们不详限制，可适用于一定类型的刚性方程组的求解。这里，我们不详细讨论刚性方程的求解。细讨论刚性方程的求解。
