第6章 常微分方程数值解法

绪论在工程和科学计算中，所建立的各种常微分方程的初值或边值问题，除很少几类的特殊方程能给出解析解，绝大多数的方程是很难甚至不可能给出解析解的，其主要原因在于积分工具的局限性。因此，人们转向用数值方法去解常微分方程，并获得相当大的成功，讨论和研究常微分方程的数值解法是有重要意义的。6.1初值问题的Euler方法,...)2,1,0()(,)(0,...),2,1,0()(),(00nyxyyxxyhnnhaxyxyyxfdxdynnnnn即上的近似值在点去获得解值解是指通过某种方法假定为常数。该式的数为步长，一般总记问题设一阶常微分方程初值初值问题的Euler方法。算出出发，逐次公式，它可以从这就是显式的）（的近似值，则有表示以于是该式可离散为：代替散化，并用方法首先将微分算子离为实现这一目标，...,,`Euler)1(,...2,1,0),()())(,()()(,Euler321010yyyynyxhfyyxyyxyxfhxyhxyxxnnnnnnnnnnn初值问题的Euler方法。才能得到步要解函数方程的不同在于，它每算一方法，它与显式公式或向后这就是隐式的）（的近似值，则有表示以于是该式可离散为：代替如果用11111101)2(Euler`Euler)2(,...2,1,0),()())(,()()(,nnnnnnnnnnnnynyxhfyyxyyxyxfhxyhxyxx初值问题的Euler方法)1(11)()1(1)(11)1(1)0(1111,||3,...2,1,0)],(),([2),(,...2,1,0)],(),([2knnknknknnnnnknnnnnnnnnnnnyyyykyxfyxfhyyyxhfyyynyxfyxfhyy取时当）（）（时常用以下迭代式：计算称为梯形公式。）（平均值的结果，则得如果取以上两式的算术初值问题的Euler方法收敛。）产生的序列时，则由（即满足常数。如果步长为条件，满足对变量设函数定理...)2,1,0}({32,120LipschitzLipschitzy),(1.1.6)(1kyLhhLhLyxfkn初值问题的Euler方法。故有由假设知）有）和（证明：由式（1)1(11)0(111)(11)(11111)(k1n1nlim,0)2(lim:||)2(......||2|),(),(|2|yy|32nknkkknnkknnknnnnyyhLyyhLyyhLyxfyxfh初值问题的Euler方法。并取即校正算法构成一类预估迭代一次一般只由于迭代工作量较大计算对于)c(11)(11)(`)(11)],(),([2),(,,,,)2(nnpnnnnncnnnnpnnyyyxfyxfhyyyxhfyyy初值问题的Euler方法))],(,()([2,Euler,...)2,1,0(),(),()(2111121211nnnnnnnnnnnnnnyxhfyxfyxfhyynkyhxhfkyxhfkkkfyy亦可写成方法该式称为改进上式还常写成初值问题的Euler方法进行比较。精确解并与法求解法和改进试分别用设初值问题例xyyyxydxdy21,EulerEuler1)0(21.1.6初值问题的Euler方法计算结果如下表所示：法：改进的法：上结果，此时计算解：取,...)2,1,0())1.0(2(1.0)2(1.0)(21Euler,...)2,1,0()2(1.0]1,0[,1.011212111nkyxkykyxykkkyynyxyyyEulerxhnnnnnnnnnnnnnxEuler法y改进的Euler法y精确解01.0000001.0000001.0000000.11.0000001.0959091.0954450.21.1918181.1840971.1832160.31.2774381.2662011.2649110.41.3582131.3433601.3416410.51.4351331.4164021.4142140.61.5089661.4859561.4832400.71.5803381.5525141.5491930.81.6497831.6164751.6124520.91.7177791.6781661.6733201.01.7847701.7378671.7320516.1.2误差概述有关，称为增量函数。与函数为而隐式单步法一般形式显式单步法一般形式为),(),,,,(),,(1111yxfhyxyxhyyhyxhyynnnnnnnnnn误差概述上的局部截断误差。称为点则由单步法计算出值计算没有误差，即的步在点果第上的整体截断误差。如称为在点则的值或隐式逐步计算，得出发，由单步法显式从初值定义11~111~11111100,)(,),()(,)(6.1.1nnnnnnnnnnnnnnxyxyTyxyyxnxyxyeyxyxy误差概述称为为局部截断误差主))(,(则非零项)())((阶精度。如果则称该算法具有误差为截断如果给果给定的算法的21..6定义121111pnnppnnnpnhxyxghOhx,yxgTp)O(hT误差概述)()(2)O()(y2Euler式对于)O(h)(2Euler对于显式则精度相对的越高。越大，误差阶，一个算法，局部截断一般431321321hOxyhThxhTxyhTｐnnnnnn：对于梯形公式：法隐：法说来6.1.3数值稳定性分析稳值不大，则，则称该算法的舍入误的舍入误差。如果算则称该算法是数值稳定||||，得值方法又算，如果结果为计算，则时有一舍入误差假设1111数，也成绝对稳定且通过某种数实际计算nnnnn~nnn~nnρρρyyρyyρy数值稳定性分析定义6.1.3若某数值算法的绝对稳定性区域包含hλ平面上的左半平面Re(hλ)＜0，则称该方法是A隐式Euler法是A6.2Runge-Kutta方法。这里仍假定，使局部截断误差，，，适当选择参数可设为方法启发，更一般算式受改进的)(),()(,...)2,1,0(),(),(Euler31112112122111nnnnnnnnnnnxyyhOyxyTnkyhxhfkyxhfkkkyyRunge-Kutta方法)()),(),(),(()()(),(),(),(:Taylor323122hOyxfyxfyxfhxyhhOyxfhkyxfhyxhfknnynnnnxnnnynnxnn展开式由二元函数)()],(),(),([)()()(3222211hOhyxfyxfyxfxyhxyynnynnnnxnnn于是Runge-Kutta方法一。自由参数，即解答不唯程，因此有一个由于四个参数，三个方展式相比较得与21211:Taylor2221方法。这是改进的此时算式为可得取Euler),(),()(21,1,21,21)1(12121121kyhxhfkyxhfkkkyynnnnnnRunge-Kutta方法.K-R)21,21(),(,21,1,0)2(1212121方法这是二阶此时算式为可得取kyhxhfkyxhfkkyynnnnnn方法。这也是二阶（又有算式可得取K-R)32,32(),()341,32,43,41)3(12121121kyhxhfkyxhfkkkyynnnnnn6.2.2四阶Runge-Kutta方法4342414323132123432421414323213133121221443322111,...)2,1,0(),(),(),(),(KR，，约定分方程初值问题可设为方法类似，解一阶常微与二阶nkkkyhxhfkkkyhxhfkkyhxhfkyxhfkkkkkyynnnnnnnnnn四阶Runge-Kutta方法稳定。但不是该公式结构整齐算式求得经典四阶类似以上推导过程，可A,),()21,21()21,21(),()22(61KR332312143211kyhxhfkkyhxhfkkyhxhfkyxhfkkkkkyynnnnnnnnnn6.2.3R-K法的稳定性nnnnyhhhhkyhhhkyhhkyhkyxy))(41)(21)(())(41)(21())(21(dd,K-R4324323221有应用于实验方程法为例以经典R-K法的稳定性nnnnnhhhhyhhhhkkkkyy]!4)(!3)(!2)()(1[:]!4)(!3)(!2)()(1[)22(61432143243211相应的误差为于是R-K法的稳定性.,,,,078.2:,1|!4)(!3)(!2)()(1|:|,|||4321稳定性必须很小才能保证算法即限制很大步长的绝对值较大时当因此可得稳定性区间为负实数时当于是得稳定性区域为由稳定性要求hhyfhhhhhnn6.2.5隐式R-K法.,,,.,A,,Euler..,||,,,K-R11又能保证数值稳定精度提高既可以构造高阶隐式方法因此限制较小即对步长定的稳它们均是但从稳定性分析知迭代求解是它们的缺点需要多次法和梯形方法计算隐式格式的向后否则数值不稳定有较大的限制对步长较大时当但从稳定性分析知使用方便可以直接计算出法的优点是从显式hyhyfyynnn隐式R-K法2221212111222121222121111122111),(),(KR约定方法形式为：设二阶隐式kkyhxhfkkkyhxhfkkkyynnnnnn隐式R-K法61)()(31211Taylor),(Taylor),(22221121221111222211221121121的同次幂系数得：比较展式的在再反复带入上式，并与，展开求得的做二元函数在点，将hyxykyxkknnnjnn隐式R-K法)32,32(),()3(41K-R320)1(2221121121kyhxhfkkkyxhfkkkyynnnnnn方法，得二阶三级隐式，令隐式R-K法)41)6341(,)6321(())6341(4,)6321(()(21K-R63216321)2(21221121121kkyhxhfkkkyhxhfkkkyynnnnnn方法得二阶四级隐式，令6.3线形多步法单步法主要依据yn的信息去计算yn+1。线性多步法是想依据yn，yn-1…，yn-r(r≥1)的信息去计算yn+1。考虑到线性组合较为方便，因此，线性多步法一般形式可设为.,,0).,...,1,0,1)(,(,)...(