八 常微分方程初值问题的数值解2

12'()11'''''()()()()()()2!!()(,),()(,)(,)(),,nppnnnnnnxypTayloryxhhyyxyxhyxyxyxPyxfxyyxfxyfxyfxy若用阶多项式近似函数有：其中。但由于公式中各阶偏导数计算复杂，不实用。§6.5Runge-Kutta(龙格-库塔)方法在高精度的单步法中,应用最广泛的是Runge-Kutta(龙格-库塔)方法一、Runge-Kutta法的基本思想（1）'(0)(0)(0)''(1)(1)(1)'''(2)(2)(2)()(1);;;;2,3,jjjjyffffyffxyffyffxyffyffjxy一般地有111112121(,)11()22(,)(,)nnnnnnnnnnEulerEuleryyhKEulerKfxyyyhKKEulerKfxyKfxhyhK如果将公式与改进公式写成下列形式：公式改进公式Runge-Kutta法的基本思想（2）11(,)()(,)(,)nnfxyyxyfxyfxy以上两组公式都使用函数在某些点上的值的线性组合来计算的近似值。Euler公式：每步计算一次的值，为一阶方法。改进Euler公式：需计算两次的值，二阶方法。Runge-Kutta法的基本思想（3）(,)(,)()nnnfxyxyTayloryxxTaylor于是可考虑用函数在若干点上的函数值的线性组合来构造近似公式，构造是要求近似公式在处的展开式与解在处的展开式的前面几项重合，从而使近似公式达到所需要的阶数。即避免求偏导，又提高了方法的精度，此为R-K方法的基本思想。11111(,)(,)(2,3,)pnniiinniininijjjyyhcKKfxyKfxahyhbKip二、二阶龙格－库塔方法11111RK(,)(,)(2,3,),(,)()pnniiinniininijjjiijinnnyyhcKKfxyKfxahyhbKipabcxyTayloryxxTaylor一般地，方法设近似公式为其中，都是参数，确定它们的原则是使近似公式在处的展开式与在处的展开式的前面项尽可能多地重合。11122122211()2(,)(,)nnnnnnyyhcKcKKfxyKfxahyhbK当p时，近似公式为112221'122'321(,)[(,)(,(,))]{(,)[(,)(,)(,)(,)]}()nnnnnnnnnnnnnnnxnnynnnnxyTayloryyhcfxycfxahyhbfxyyhcfxycfxyahfxyhbfxyfxyOh上式在处的展开式为12''232221()(,)[(,)(,)(,)]()nnnxnnynnnnyccfxyhcafxybfxyfxyhOh12'312''3()()()()()()2(,)[(,)(,)(,)]()2nnnnnnnnnxnnynnnnyxxTaylorhyxyxhyxyxOhyfxyhhfxyfxyfxyOh在处的展开式为'1122'321{(,)[(,)(,)(,)(,)]}()nnnnnnxnnynnnnyyhcfxycfxyahfxyhbfxyfxyOh122222111/21/2),cccacbO3有无穷多组解，每一组解得一近似公式，局部截断误差均为(h这些方法统称二阶方法。122211121211,1,2()/2(,)(,)nnnnnnccabEuleryyhKKKfxyKfxhyhK取此为改进公式。近似公式为122211212110,1,,2(,)(2,2)nnnnnnccabyyhKKfxyKfxhyhK取此为常用的二阶公式，称为中点公式。／／三、三阶龙格－库塔方法11231213123(4)6(,)(,)22(2,)nnnnnnnnpRKRKhyyKKKKfxyhhKfxyKKfxhyhKhK类似地，对，即三个点，通过更复杂的计算，可导出三阶公式。常用的三阶公式为：四、四阶龙格－库塔方法1123412132434(22)6(,)(,)22(,)22(,)nnnnnnnnnnpRKRKhyyKKKKKfxyhhKfxyKhhKfxyKKfxhyhK对，即四个点，可导出四阶公式。常用的四阶公式为：解：201()dyxydxyy01(,)x例2：用经典的Runge-Kutta方法求解下列初值问题。01.h经典的四阶Runge-Kutta公式：11234226()nnhyykkkk12;nnnxkyy2112222()nnnhxhkykhyk4332()nnnxhkyhkyhk3222222();nnnhxhkykhyk0.10.20.30.40.51.09541.18321.26491.34161.4142nxny0.60.70.80.91.01.48321.54921.61251.67331.7321nxny同保留5位的精确值完全一致：0.10.20.30.40.51.09541.18321.26491.34161.4142nxny0.60.70.80.91.01.48321.54921.61251.67331.7321nxny21yx二、高阶和隐式Runge-Kutta方法注：对于显式N级R-K方法，最多只能得到N级方法；N1,2,3,45,6,78,910,11,…NN-1N-2()pN已经证明N级R-K方法的阶具有下列关系：()pN2N若要得到N阶以上方法，则使用N级隐式R-K方法N级隐式R-K方法的一般形式：11Nnniiiyyhck11(,);,...,NininijjjkfxahyhbkiNN级隐式R-K法可以达到2N阶(1)一级二阶的隐式中点方法：11nnyyhk1122(,)nnhkhkfxy(2)二级四阶的隐式R-K方法：1212()nnhkkyy112131326446[(),()]nnhkkfxhyhk221131326464[(),()]nnhkkfxhyhk三、变步长方法基本思想：根据精度自动地选择步长对于经典Runge-Kutta方法：012,,,nStep1:设从出发，以为步长，经过一步计算得到nxh511()()hnnyxyChStep2:取为步长，再从出发，经过两步计算得到nx2h521122()()()hnnhyxyC21111116(/)()()()hnnhnnyxyyxy2111116(/)()(())()hhnnnnyxyyxy221111115(/)(/)()()[]hhhnnnnyxyyy211(/)()||hhnnyy记如果，则将步长折半进行计算，直到为止此时取为最终结果；21(/)hny如果，则将步长加倍进行计算，直到为止此时将步长折半一次计算，得到的为最终结果。一、收敛性/*Convergence*/§3单步法的收敛性、相容性和绝对稳定性1(,,)nnnnyyhxyh对于初值问题的一种000(,)();dyfxydxyxyxx()1Def单步法产生的近似解，如果对于任一固定的，均有，0nxxnh则称该单步法是收敛的。0lim()nnhyyx类似地可以定义隐式单步法、多步法（§4）的收敛性31.Th设初值问题（*）对应的下列单步法是阶的，1(,,)nnnnyyhxyhp且函数满足对的Lipschitz条件，即存在常数y0L121212|(,,)(,,)|||,xyhxyhLyyyy则该单步法是收敛的，且()()pnnyxyOh证明：()nnneyxy记由截断误差的定义11()()(,(),)nnnnnyxyxhxyxhT11[(,(),)(,,)]nnnnnnneehxyxhxyhT因为单步法是阶的：p000,hhh满足11||pnTCh11||||||pnnneehLeCh||ne其中11,phLCh212||||||nnneee3231||()ne2101||||(...)nnnee10001||exp[()]||{exp[()]}pnnneLxxeChLLxx101{exp[()]}pnChLLxx()pneOh00()h二、绝对稳定性/*AbsoluteStibility*/计算过程中产生的舍入误差对计算结果的影响首先以Euler公式为例，来讨论一下舍入误差的传播:1(,)nnnnyyhfxy设实际计算得到的点的近似函数值为，nnnyynxny其中为精确值，为误差n1(,)nnnnyyhfxy111nnnyy11[(,)(,)][(,)]nnnnnnynnhfxyfxyhfx如果，则误差是不增的，故可认为是稳定的11||yhf例如：对于初值问题0()yyyxa精确解为0xxyae而实际求解的初值问题为0()yyyxaa精确解为0()xxyaae在处的误差为nx0nxxae可见误差随着的增加呈指数函数增长nx如果初值问题为0()yyyxa精确解为0xxyae实际求解的初值问题为0()yyyxaa精确解为0()xxyaae在处的误差为nx0nxxae可见误差随着的增加呈指数函数递减nx当时，微分方程是不稳定的；0yf而时，微分方程是稳定的。0yf上面讨论的稳定性，与数值方法和方程中有关f实验方程：0,Re()yyC1(,,)nnnnyyhxyh对单步法应用实验方程，1()nnyEhy3Defh如果，当时，则称该1()Eh单步法是绝对稳定的，在复平面上复变量满足1()Eh的区域，称为该单步法的绝对稳定域，它与实轴的交集称为绝对稳定区间。1111()()pnnnTyxyOh若单步法是阶的，则p由实验方程可得：1()exp()nnyxyh1exp()()()pnnyhEhyOh()exp()Ehh11()()nnnnyyhyhy例3：分别求Euler法和经典的R-K法的绝对稳定区间。1()Ehh解：Euler公式：1(,)nnnnyyhfxy将其应用于实验方程绝对稳定域：11h1120hh当时，R绝对稳定区间：20(,)经典的R-K公式：11234226()nnhyykkkk1(,)nnkfxy211222(,)()nnnhhhkfxykyk322222(,)()nnnhhhkfxykyk433(,)()nnnkfxhyhkyhk22()nhy23224()nhhy3nyhk2341234()()()()!!!hhhEhh22334411234()!!!nnhhhyhyny当时，R绝对稳定区间：27850(.,)23411234()()()()!!!hhhEhh