第八章常微分方程数值解

1、第八章 常微分方程数值解由常微分方程理论可知，我们只能求一些特殊类型的常微分方程。而实际上许多常微分方程求解非常困难。本章主要讨论一阶常微分方程的初值问题： （81）从理论上讲只要方程中的连续且关于满足李普希兹（Lipschitz）条件，即存在常数L，使则常微分方程存在唯一解。微分方程数值解：就是求微分方程的解在一系列离散节点处的近似值（i=1，2，n） 称为由到的步长，通常取为常数h。求数值解，首先将微分方程离散化，常用方法有：（1） 用差商代替微商若用向前差商代替微商，即 （i=1，2，n）则得即（2） 数值积分法利用数值积分法左矩形公式 =可得同样算法（3） 用泰勒（Taylor）公式得

2、离散化计算公式§1 欧拉（Euler）方法1.1欧拉方法对一阶微分方程（81），等分区间为份，则 由以上讨论可知，无论用一阶向前差商，还是用数值积分法左矩形公式，或者用泰勒公式取前两项都可得到同样的离散化计算公式代入初值则得到数值算法： （i=1，2，n1） （82）称其为欧拉方法。几何上欧拉方法就是用一条折线近似表示曲线。（如图）： P0 Pi Pi+1 O 1.2欧拉方法的误差估计定义1 局部截断误差：假设为准确值，用某数值算法计算产生的误差，称为该数值算法的局部截断误差。定义2 整体截断误差：准确解与数值解的误差，。设有二阶导数，由泰勒公式有：= 所以 = ，（83）当h充分小

3、时，欧拉方法的局部截断误差与h2是同阶无穷小，称其为一阶方法。定义3 如果一数值解法的局部截断误差为，则称该算法为阶算法。1.3 改进的欧拉方法由微分方程数值解的三种基本构造方法知，若取不同的差商（如向后），不同的数值积分公式（如梯形公式），以及泰勒公式取前三项、四项等可得不同的算法。如果用梯形公式计算积分： （84）且 = （85）由于此方程为的隐式方程，不易求解。一般将其与欧拉方法联合使用。可得算法 （86）（k=0，1，2，；i=1，2，n1）实际计算中，当比较小时，常取一次迭代后的近似值为，于是有改进的欧拉方法 （i=0,1,2,n-1）例1 用欧拉方法和改进的欧拉方法求微分方程 的数

4、值解（取h=0.1）。 解：由欧拉方法（82），得数值计算公式=0.1× 计算结果如表8-1 由改进的欧拉方法（86），得数值计算公式 计算结果如表8-2表8-1xi 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 yi 1.0000 1.0069 1.0208 1.0391 1.0628 1.0923 1.1269 1.1643 误差 0.0000 0.0037 0.0077 0.0993 0.0120 0.0151 0.0189 0.0222 表8-2xi 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 yi 1.0000 1.0033 1.013

5、2 1.0292 1.0506 1.0773 1.1079 1.1422 误差 0.0000 0.0000 0.0000 0.0002 0.0020 0.0101 0.0122 0.0053 例2 用欧拉法、改进欧拉法求微分方程数值解（h=0.1）。 解 由欧拉方法（82），得数值计算公式=0.1由改进的欧拉方法（86），得数值计算公式 计算结果如下表8-3（Euler） （改进Euler） 0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.90001.00001.0000 1.1000 1.2009 1.3043 1.41

6、18 1.5246 1.6443 1.7722 1.9100 2.0591 2.2213 1.0000 1.1005 1.2027 1.3083 1.4190 1.5362 1.6614 1.7965 1.9430 2.1030 2.2782 00.00050.00180.00400.00720.01150.01720.02430.03310.04390.0569§2 龙格库塔（RungeKutta）方法2.1 泰勒展开法由于欧拉方法为一阶方法，为了提高算法的阶，有必要讨论更高阶的方法。在泰勒展开式中取更多的项，如取p1项可得p阶算法。其中）可用复合函数求导法则计算。如p=2时得二阶

7、泰勒方法2.2 龙格库塔法 为了避免计算高阶导数龙格库塔方法利用某些点处的值的线性组合构造计算公式，使其按泰勒公式展开后与初值问题解的泰勒展开式比较，有尽可能多的项相同。 龙格库塔法的一般形式为： （87）下面以二阶龙格库塔法为例说明龙格库塔法的构造过程。 （88）将K2在处按泰勒公式展开，则 = =另一方面=所以为使局部截断误差的阶尽可能高，应使 方程组有无穷多组解，取定参数则得到许多具体的二阶龙格库塔公式。如取 =， =1 则得龙格库塔法即改进的欧拉公式。同理可得更高阶的龙格库塔法。常用的四阶（经典）龙格库塔法： （89）例3 用二阶龙格库塔法求例1（ ，h取0.1）的数值解。解 由（89

8、）得计算得结果见表8-4 表8-4_xi 0.0 0.1 0.2 0.3 0.4 0.5 yi 1.000000 1.003322 1.013159 1.029124 1.050718 1.077217误差 0.000000 0.000002 0.000008 0.000011 0.000024 0.000031 _§3阿达姆斯（Adams）方法将初值问题写成等价形式记，当k取不同的值以及对用不同的插值多项式近似时得到不同的数值算法。3.1阿达姆斯外插公式上式中取k=0，并取为节点，作函数的三次插值多项式 =其插值余项为 则=积分后得阿达姆斯外插公式 （n=3，4，5） （810）其

9、局部截断误差为 Rn= （811）由于积分区间为，插值区间为，因此称其为阿达姆斯外插公式。3.2 阿达姆斯内插公式同样取k=0，并取为节点，作函数的三次插值多项式可得 （812） 余项Rn （813）由于阿达姆斯内插公式是隐式公式，故用它计算时也需用迭代法。例4 用四阶阿达姆斯外插公式求初值问题 的数值解，取h=0.1。 解 由（810）得 = （n=3，4，5） 计算结果见表8-5 表8-5_xi 0.0 0.1 0.2 0.3 0.4 0.5 yi 0 0.004837 0.018731 0.040818 0.070323 0.106535误差 2.87×10-6 4.82

10、15;10-6_§4 收敛性与稳定性4.1收敛性 定义4 如果一数值方法对任意固定的点，当h=时有，则称该方法是收敛的。定理1 如果关于y满足李普希兹（Lipschitz）条件，即存在常数L，使，且有界，则欧拉方法的整体截断误差满足其中M=。（证略）由整体截断误差估计式知，当初值误差为0时有当时，欧拉方法收敛。4.2 稳定性一个即使收敛的数值方法，由于初值一般带有误差，且计算过程中不断产生舍入误差，随着误差的传播，对计算结果可能产生很大的影响。定义5 设用某数值方法计算时，所得实际结果为，且由误差引起以后各节点处（j>i）的误差为，如果总有，则称该算法是绝对稳定的。由于稳定性的

11、讨论比较复杂，常用试验方程=(其中为常数)，并把能使某一数值方法绝对稳定的的允许取值范围称为该方法的绝对稳定域。对试验方程应用欧拉算法得=（1h）设上有误差，它的传播使产生误差i1，设y*i=按欧拉公式计算得y*i1=时不产生新的误差，则 y*i1=（1h）y*i 两式相减，可得i1=（1h）为保证误差在以后的计算中不至增大，应选h使 即 （814）当h满足上式时，则称欧拉方法是条件稳定的。上式称为稳定条件。对一般方程，可近似地取以便判断稳定性，并确定计算时的步长。例5 对初值问题 求欧拉方法的稳定条件。解： 因为=-20 所以取步长为 =0.1时欧拉方法是条件稳定的。§5 方程组与

12、高阶方程的数值解法5.1 一阶方程组的数值解法设初值问题由单个方程的经典龙格库塔法可得其中 （）（）（）（）5.2 高阶方程的数值解法对于高阶常微分方程初值问题，原则上总可转化为一阶方程组来来解。例如，二阶常微分方程初值问题作变换z= 则方程化为一阶方程组则可用一阶方程组的数值解法来计算。习 题 八1 用欧拉方法和改进欧拉方法求初值问题的数值解（取h=0.5），并将计算结果与准确解比较。2 对于初值问题 （x>0）分别导出欧拉方法和改进欧拉方法的表达式。3 用四阶龙格库塔方法求初值问题：的数值解。4 用阿达姆斯外插法求的数值解。取h=0.1。5 试求向后欧拉方法的绝对稳定区域。6 对初值问题用欧拉方法求数值解，分别取h=0.1，h=0.01，将所得结果与精确 比较，并用稳定条件说明比较结果。7 用改进欧