数值计算方法第六章.

第六章常微分方程初值问题的数值解法6.1欧拉方法6.2计算公式的误差分析6.3龙格－库塔方法6.4向一阶方程组与高阶方程的推广问题的提出的数值解求初值问题)16(],[)(),(000bxxyxyyxfy数值求解方法的近似值数值逐个求解出节点上的函生成节点作等距分割对区间)(,),(),()2(),,,1,0(,],[)1(21000NixyxyxyNxbhniihxxbx6.1欧拉方法6.1.1欧拉公式与改进欧拉公式算法：Nxbhniihxxbxi000),,,1,0(,],[生成节点作等距分割对区间处的近似值在处的精确值在iiiixxyyxxyxy)()()(111))(,()()())(,()(),(],[11iiiiiixxiixxxxiidxxyxfxyxydxxyxfdxxyyxfyxx即：积分上对在)46())(,()()(11iixxiidxxyxfxyxy化为积分方程问题从而把微分方程问题转选择不同的数值积分公式来求近似值就得到初值问题的各种数值解法1))(,(iixxdxxyxf1.欧拉公式))(,()()(1iiiixyxfhxyxy得，时由当)46())(,())(,(1iixxxyxfhdxxyxfiioxy),(yxfY0x1x1ixbxnixh)26()1,,1,0()(),(001Nixyyyxfhyyiiii由此可建立计算公式这称为欧拉公式2.后退欧拉公式oxy),(yxfY0x1x1ixbxnixh得，时由当)46())(,())(,(111iixxxyxfhdxxyxfii)36()1,,1,0()(),(00111Nixyyyxfhyyiiii由此可建立另一公式))(,()()(111iiiixyxfhxyxy这称为后退欧拉公式后退欧拉公式是一个隐式公式，通常采用迭代法求解。比较并与精确解的数值解xxexxyxyyxey)2(21)(,]1,0[1)0(2例6.1以h=0.1为步长，用欧拉法求常微分方程初值问题)9,2,1,0(1)(01iyyexhyyixiiii解：由欧拉公式得xiyiy(xi)y(xi)-yi01100.10.9000000.9093629.362x10-3………..….….立表(见表6-1（书121页）)6.1.2梯形公式与改进欧拉公式3.梯形公式)56()()],(),([2))](,())(,([2)()())](,())(,([2))(,(00111111111xyyyxfyxfhyyxyxfxyxfhxyxyxyxfxyxfhdxxyxfiiiiiiiiiiiiiiiixxii从而导出梯形公式得由---梯形公式也是隐式单步法公式oxy),(yxfY0x1x1ixbxnixh图1梯形公式用梯形公式计算时，通常取欧拉公式的解作为迭代初值进行迭代计算，即采用下式),1,0()],(),([2),()(11)1(1)0(1kyxfyxfhyyyxfhyykiiiiikiiiii)1(1)(1)2()3(1)2()2(1)2()1(1)2()0(1)1(0kkyyyyyyy(1)(2))1(2)(2)2()3(2)2()2(2)2()1(2)2()0(2)1()0(1kkyyyyyyy11)(1)1(1)(1)(,||,iikikikiyxyyyyk的近似值值作为函数取时使得迭代计算到某个)1,,1,0()2()]~,(),([2)1(),(~,11111)1(1Niyxfyxfhyyyxfhyyyyiiiiiiiiiiii则得到显式计算公式代替取直接简单地4.改进欧拉公式这称为改进欧拉公式NNyyyyyyyyy)2(3)2(3)1(2)2(2)1(1)2(1)1(0~~~~例6.2仍取步长h=0.1，采用改进欧拉法重新计算例6.1的常微分方程初值问题。)9,1,0(1)]~()[(2)(~011111iyyexyexhyyyexhyyixiixiiiixiiiiii(见表6-2（书125页）)这时改进欧拉公式为解xiyiy(xi)y(xi)-yi01100.10.9095240.9093621.363x10-4………..….….立表6.2计算公式的误差分析定义6.1若yi+1是yi=y(xi)从计算得到的近似解，则称y(xi+1)－yi+1为所用公式的局部截断误差，简称为截断误差。定理6.1若单步法yi+1=yi+h(xi,yi,h)的局部截断误差为O(hp+1)，且增量函数(x,y,h)关于y满足李普希兹条件，即存在常数L0，使对成立不等式yy~,|~||),~,(),,(|yyLhyxhyx则其整体截断误差y(xi)－yi=O(hp)截断误差的估计(基本假设：yi=y(xi))设y(x)C3[x0,b],则)76()()(2)()()()(321hOxyhxyhxyhxyxyiiiii（1）对欧拉公式，有)86()()()(2)()()())(,()(232111hOhOxyhyxyxyhxyxyxfhxyyiiiiiiiii故因此，欧拉公式的局部截断误差为O(h2)（2）对后退欧拉公式，有)96()()()(2)())(,()(23211111hOhOxyhyxyxyxfhxyyiiiiiii故因此，后退欧拉公式的局部截断误差为O(h2)（3）对梯形公式，注意到其公式可改写为)],([)],([21111iiiiiiiyxfhyyxhfyy故由式（6-8）和（6-9）得)106()()]()(2[)]()(2[21)],([)()],([)(21)(33232111111hOhOxyhhOxyhyxfhyxyyxfhyxyyxyiiiiiiiiiiii因此，梯形公式的局部截断误差为O(h3)（4）对改进欧拉公式，有)]~,(),([2)()()(~1111iiiiiiiiiyxfyxfhxyyxyhxyy而由，故有),(),(),(yxfyyxfyyxfyyx得)()(2)()()()()()(][))(,())()(,()~,(32122))(,(11hOxyhxyhxyyhOxyhxyhOfyhfhxyxfxyhxyhxfyxfiiiiiixyxyxiiiiiiiii所以与式（6-7）比较得y(xi+1)－yi+1=O(h3)因此，改进欧拉公式的局部截断误差为O(h3)定义6.2若一种求解常微分方程初值问题的数值计算方法的局部截断误差为O(hp+1)，则称该方法为p阶精度，或称该方法为p阶方法。由此定义知，欧拉方法与后退欧拉方法为一阶精度，梯形法与改进欧拉方法为二阶精度。6.3龙格-库塔方法由中值定理，有),(,))(,()()()()(111iiiiiixxyfhyxxxyxy因此，以上介绍的各种单步法本质上都是对平均斜率f(,y())进行近似，龙格-库塔据之提出了适当选取若干点上的斜率值作近似以构造高精度计算公式的方法，其基本思想是基于泰勒展式的待定系数法。6.3.1二阶Rung-Kutta公式问题：建立二阶精度的计算格式形为待定系数为这里babhKyahxfKyxfKKKhyyiiiiii,,,)126(),(),()(2112122111在y(xi)=yi的假设下，有)(][)()(][))(,())(,()())(,(2))(,(2))(,(1121hOfybfahxyhOfbhKfahxyxfbhKxyahxfKxyxyxfKiiiixyxyxixyxyxiiiiiii)(][)()()(3))(,(222211hOfybfahxyhxyyiixyxyxiii故解)(][)()()(3))(,(222211hOfybfahxyhxyyiixyxyxiii变形欧拉公式根据格式为二阶精度，即y(xi+1)－yi+1=O(h3)比较两式系数得)136(2/12/112221ba而)(]2121[)()()()()2/()()()(3))(,(2321hOfyfhxyhxyhOxyhxyhxyxyiixyxyxiiiiii系数满足(6-13)的形为(6-12)计算格式统称为二阶R-K公式。当令1=1/2时，解得2=1/2，a=b=1，即为改进欧拉公式。若令1=0，解得2=1，a=b=1/2，则得另一计算公式)146()2/,2/(),(12121hKyhxfKyxfKhKyyiiiiii6.3.2四阶R-K公式每一步需计算的f值的个数1234567n8精度阶1234456n-21965年，Butcher研究发现显式R-K公式的精度与需要组合的斜率值的个数具有如下关系可见，超过四阶精度的R-K公式效率并不高，实际计算通常选用如下四阶格式)156(),()2/,2/()2/,2/(),()22(6342312143211hKyhxfKhKyhxfKhKyhxfKyxfKKKKKhyyiiiiiiiiii经典R-K公式这时经典R-K公式为例6.3取步长h=0.2，采用经典R-K法计算例6.1的常微分方程初值问题。)()()2()2()2()2()22(63)(42)2/(31)2/(2143211hKyehxKKhyehxKKhyehxKyexKKKKKhyyihxiihxiihxiixiiiiiii取h=0.2计算得到表6-4(书133页）。与例6.1和例6.2比较可见，用经典R-K法计算得到的解比用欧拉法和改进欧拉法所得到的解精确得多。解6.3.3步长的自动选择对于p阶精度的计算格式，当取步长为h时，记为从y(xi)计算得到的y(xi+1)(xi+1=xi+h)的近似解，则有][1hiy)(121)(21)()(~,)~(2~2)()()(][1]2/[1]2/[11][11]2/[111]2/[111][11hihiphiiphiihiiphiiphiiyyyxyyxyyxyCChChCyxyCChyxy所以故有较小时为常数为常数为便于进行事后误差估计，实际计算时通常采用步长减半算法。记，则对给定的精度要求，可根据按如下方式调整步长：||][1]2/[1hihiyy（1）若，则把步长逐次减半计算，直至为止，这时最终得到的解即为满足精度要求的近似解。（2）若，则把步长逐次加倍计算，直至为止，这时取前一次步长计算所得到的解作为满足精度要求的近似解。