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

1计算方法华中科技大学数学与统计学院第六章常微分方程初值问题的数值解法计算方法课程组2§6.1基本离散方法§6.2Runge-Kutta方法§6.3线性多步法§6.4收敛性与稳定性§6常微分方程数值解法3考虑一阶常微分方程的初值问题:例如：,[,]()0101y=x+yxy其解析解为：12xyxe§6.1基本离散方法0)(],[),(yaybaxyxfdxdy4[,]xdtx01012-tye但是,只有一些特殊类型的微分方程问题能够得到用解析表达式表示的函数解，而大量的微分方程问题很难得到其解析解。因此，只能依赖于数值方法去获得微分方程的数值解。例如：,[,]()01012-xy=exy其解析解为：很难得到其解析解5例如：,[,]()0101y=x+yxy其解析解为12xyxe只有一些特殊类型的微分方程问题能够得到用解析表达式表示的函数解，而大量的微分方程问题很难得到其解析解。因此，只能依赖于数值方法去获得微分方程的数值解。要计算出解函数y(x)在一系列节点a=x0x1…xn=b处的近似值),...,1()(nixyyii通常取节点间距为步长，通常采用等距节点，即取hi=h(常数)。它适合计算机求解,应用广泛,具有应用价值。)1,...,0(1nixxhiii6Taylor级数复习函数在点作Taylor级数展开：()yyx0xx这里x,x0都可以是任意一点。00000(n)n001()()()()()()2!()()2yx=yx+yxx-x+yxx-x1++yxx-x+n!7Taylor级数复习()yyx0xx()()()()()()()!2i+1iiiii2n!nnhhyxyx+hyxhyxyxyx()()()()()()()!2i-1ii()2nniiih-hyxyx-hyxhyxyxyxn!则：()()()()()!2i+1iiii()2n(n+1)ihhyxyxhyxhyxyxyxn!()()()()()!2i-1iiii()()2n(n+1)ih-hyxyx-hyxhyxyxyxn!(n)n00000001()()()()()()()()2!21yx=yx+yxx-x+yxx-x++yxx-x+n!86.1.2Euler公式利用等距分割，数值微分来代替导数项，建立差分格式。ibaxim-=1、向前差商公式)(''2)(')()(1nnnnyhxyhxyxy)(''2))(,()()(1nnnnnyhxyxfhxyxy)(''2))(,()()(21nnnnnyhxyxhfxyxy所以，可以构造差分方程),(1nnnnyxhfyy称为局部截断误差。显然，这个误差在逐步计算过程中会传播，积累。因此还要估计这种积累9定义:在假设yi=y(xi)，即第i步计算是精确的前提下，考虑的截断误差Ri=y(xi+1)yi+1称为局部截断误差。定义:若某算法的局部截断误差为O(hp+1)，则称该算法有p阶精度。收敛性：)(''2))(,()()(21nnnnnyhxyxhfxyxy),(1nnnnyxhfyy考察局部误差的传播和积累102、向后差商公式)(''2)(')()(11nnnnyhxyhxyxy)(''2))(,()()(111nnnnnyhxyxfhxyxy)(''2))(,()()(2111nnnnnyhxyxhfxyxy),(111nnnnyxhfyy是隐格式，要迭代求解)0(1)(11)1(1),(nknnnknyyxhfyy可以由向前差商公式求出113、中心差商公式)(''2)(')()(111nnnnyhxyhxyxy),(1111nnnnyxhfyy是多步，2阶格式，该格式不稳定12对微分方程],[,),('baxyxfdxdyy积分有:1),(nnxxyxfdxdy)(21hOen类似，可以算出其误差估计式：2阶的方法所以，有)],(),([2111nnnnnnyxfyxfhyy是个隐式的方法，要用迭代法求解11()()(,())nnxnnxyxyxfxyxdx)('''12))](,())(,([2)()(2111fhxyxfxyxfhxyxynnnnnn局部截断误差4、梯形公式135、欧拉公式的改进：•隐式欧拉法向后差商近似导数hxyxyxy)()()(011x0x1))(,()(1101xyxfhyxy)1,...,0(),(111niyxfhyyiiii由于未知数yi+1同时出现在等式的两边，不能直接得到，故称为隐式欧拉公式，而前者称为显式欧拉公式。14•中点欧拉公式（欧拉二步法）中心差商近似导数hxyxyxy2)()()(021x0x2x1))(,(2)()(1102xyxfhxyxy1,...,1),(211niyxfhyyiiii假设，则可以导出即中点公式具有2阶精度。)(),(11iiiixyyxyy)()(311hOyxyRiii需要2个初值y0和y1来启动递推过程，这样的算法称为双步法/*double-stepmethod*/，而前面的三种算法都是单步法/*single-stepmethod*/。方法显式欧拉隐式欧拉梯形公式中点公式简单精度低稳定性最好精度低,计算量大精度提高计算量大精度提高,显式多一个初值,可能影响精度15•改进欧拉法Step1:先用显式欧拉公式作预测，算出Step2:再将代入隐式梯形公式的右边作校正，得到)1,...,0(),(,),(211niyxfhyxfyxfhyyiiiiiiii1111(,)[(,)(,)]2nnnnnnnnnnyyhfxyhyyfxyfxy1ny1ny16Euler方法、隐式Euler方法、梯形方法与单步法计算公式的显式单步法对应关系),,(1hyxhyykkkk隐式单步法),,,(11hyyxhyykkkkk),(),,(kkkkyxfhyx),(1kkkkkyxfhyy显式Euler方法隐式Euler方法),,(111kkkkkyxfhyy)],(),([2111kkkkkkkyxfyxfhyy梯形方法(隐式)),(),,,(111kkkkkyxfhyyx),(21),(21),,,(111kkkkkkkyxfyxfhyyx6.1.3总结17算例：分别用Euler公式和改进的Euler公式求解：,()()201100yxxyxy1.0h取步长，计算y(0.5)的近似值解：欧拉公式：.(),,,,2i1iiii001i0190yyxxyy改进的Euler公式：.()./(),,,,yyxxyyyxxyxxy2iiiii+122i+1iiiii+1i+1i+1001012i0190y12-xy=xx+-e精确解：18算例分别用显式Euler方法，梯形方法和预估－校正Euler方法初值问题解：取h=0.1，),(1kkkkkyxfhyy)1(kkkkxyhy1),(xyyxfkkkkkhxhyh)1(101100109kyk0101.,..kkkhxkhk===,1xydxdy10x1)0(y(1)Euler方法为:续19算例分别用显式Euler方法，梯形方法和预估－校正Euler方法解初值问题解：取h=0.1，1),(xyyxf,1xydxdy10x1)0(y梯形方法为：)],(),([2111kkkkkkkyxfyxfhyy)]1()1[(211kkkkkkxyxyhy).1(1.0)1(,1.0,1.01khkxkkhxhkkkkk10110521191kyykk续20算例分别用显式Euler方法，梯形方法和预估－校正Euler方法解初值问题解：取h=0.1，1),(xyyxf,1xydxdy10x1)0(y梯形方法为：)],(),([2111kkkkkkkyxfyxfhyy)]1()1[(211kkkkkkxyxyhy).1(1.0)1(,1.0,1.01khkxkkhxhkkkkk10110521191kyykk1.00095.0905.01kyykk预估－校正Euler方法：续21Euler方法梯形方法预估－校正方法0.01.0000000.01.0000000.01.0000000.00.11.0000004.8×10－31.0047627.5×10－51.0050001.6×10－40.21.0100008.7×10－31.0185941.4×10－41.0190252.9×10－40.31.0290001.2×10－21.0406331.9×10－41.0412184.0×10－40.41.0561001.4×10－21.0700962.2×10－41.0708004.8×10－40.51.0904901.6×10－21.1062782.5×10－41.1070765.5×10－40.61.1314411.7×10－21.1485372.7×10－41.1494045.9×10－40.71.1782971.8×10－21.1962952.9×10－41.1972106.2×10－40.81.2304671.9×10－21.2490193.0×10－41.2499756.5×10－40.91.2874201.9×10－21.3062643.1×10－41.3072286.6×10－41.01.3486781.9×10－21.3675733.1×10－41.3685146.6×10－4kxky)(kkxyy)(kkxyy)(kkxyykyky数值例子表明，梯形方法和预估－校正Euler方法比显式Euler方法有更好的精度。续22基本思想:根据微分中值定理有：()()()[]11i+iii+yxyxyξξx,xh()()()1i+iyxyxhyξ因此()(),01iiyxhyxθhθ()(,())iiiyxhfxθhyxθh*()iyxhK*(,())iiKfxθhyxθh这里表示在区间上函数[,]1ii+xx()yyx的平均斜率。§6.2龙格-库塔方法建立高精度的单步递推格式。236.2Runge－Kutta法由Taylor展开)()!1()(!)(')()()1(1)(1nkknkknnnykhxykhxhyxyxy记为1nhT),()('yxfxy'),(),()(''yyxfyxfxyyx)('''xy所以，可以构造格式),(),(),(!2),(21nnnnynnxnnnnyxfyxfyxfhyxhfyy这种格式使用到了各阶偏导数，使用不便。从另一个角度看，11)),(,,()()(nnnnnhTfxyxhhFxyxy取(x,y)及其附近的点做线性组合，表示F，问题就好办了。当然，要求此时的展开精度相同。这种方法称为Runge－Kutta法24在(x,y)处展开有)(),(),(),(),(),(),,,(221221hOyxfyxhfb