2Runge-Kutta方法

一、高阶显式单步法的构造方法§2龙格-库塔（Runge-Kutta）法显式单步法的一般形式：1(,,)nnnnyyhxyh构造高阶方法，即如何确定增量函数(,,)xyh使得该方法的局部截断误差的阶数尽量高。111()()pnnyxyOh为尽可能大的整数p1、Taylor级数方法设满足初值问题：000(,)();dyfxydxyxyxx()()yx2112()()()()!nnnnnhyxyyxhyxyx则111()()pnnyxyOh1()()()(,,)!pppnnnnhyxOhyhxyhp12()(,,)()()()!!pphhxyhyxyxyxp若取的各阶导数的计算()yx()(,())dyxfxyxdx22()(,())......dyxfxyxdx()(,())yxfxyx(,())(,())()xyxyxdfxyxfffdx(,())Fxyx2222(,())(,())(())xxxyyyyxyxyxdfxyxfffffffffdx(,())(,())(,())yGxyxfxyxFxyx以此类推，可以求出的各阶导数，则有()yx223(,())(,,)(,)(,)()!!yxyxhhxyhfxyFxyGfF(,,)(,())xyhfxyx即为Euler方法2、Runge-Kutta方法基本思想：利用在某些特殊点上的函数值的线性组合来构造高阶单步法的增量函数(,)fxyN级（阶）Runge-Kutta方法的一般形式：1(,,)nnnnyyhxyh1(,,)Niiixyhck其中1(,)kfxy1123(,),,,...,iiiijjjkfxahyhbkiN11iiijjabN=1：Euler方法当N1时，适当选取式中参数，使该方法的阶数尽量高前式中参数的求解方法利用二元函数的Taylor展开公式：01(,)()(,)!iifxhykhkfxyixy将每个展开，然后代入ik1(,,)nnnnyyhxyh从而得到关于的方程组，解之即得。,,iijiabc223(,())(,,)(,)(,)()!!yxyxhhxyhfxyFxyGfF与Taylor级数法中的增量函数比较幂次相同项的系数h二级方法：N=21(,)kfxy22211(,)kfxahyhbk2211(,)[]xyfxyahfhbkf2222322211211122[()()]()xxxyyyahfabhkfbhkfOh2(,)[]xyfxyahfff223222()[(]()xxxyyyahfffffOh2232212222(,,)()()cahxyhccfcahFGOh223(,,)()!yhhxyhfFGfF比较关于的幂次相同项的系数，得到方程组：h1222112ccca方程组有无穷多解：二级方法有无穷多种常见的3种二级方法：中点法（修正的Euler法）1221012,,cca取122(,(,))nnnnnnhhyyhfxyfxyRunge-Kutta二级方法122112,cca取12[(,)(,(,))]nnnnnnnnhyyfxyfxhyhfxy3()OhHeun（休恩）二级方法2232212222(,,)()()cahxyhccfcahFGOh223(,,)()!yhhxyhfFGfF2h要求项的系数尽量相同得到方程组：121;cc2212;ca22213ca122132443,,cca1223433[(,)(,(,))]nnnnnnnnhhyyfxyfxhyfxy三级方法：N=3类似于N=2的推导方法，可得到1231;ccc223312;caca22223313;caca323216cab4()Oh常见的2种三阶方法：Kutta三阶方法112346()nnhyykkk1(,);nnkfxy2122(,)nnhhkfxyk3122(,)nnkfxhyhkhkHeun三阶方法（见教材）四级方法：N=45()Oh局部截断误差常见的2种四阶方法：经典Runge-Kutta方法Kutta四阶方法（见教材）11234226()nnhyykkkk1(,)nnkfxy2122(,)nnhhkfxyk3222(,)nnhhkfxyk43(,)nnkfxhyhk解：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二、相容性/*Consistency*/()()(,(),)yxhyxhxyxh0[()()...]()[(,(),)...]yxhyxyxhxyx0[()(,(),)]...hyxxyx1()pOh100()(,(),)pyxxyx对于阶方法：1(,,)nnnnyyhxyhp()若方法（**）的增量函数满足：2Def0(,,)(,)xyfxy则称该方法与初值问题（*）相容。设方法（**）与初值问题（*）相容，且满足L-条件，则该方法（**）是收敛的，即当固定，时nxx0h()nnyyx1(,,)nnnnyyxyhh0(,(),)()nnnxyxyx()(,())nnnyxfxyx再由相容性得：上式说明：当时，方法（**）趋于原微分方程0h本章讨论的数值方法都是与原初值问题相容的三、绝对稳定性/*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法的绝对稳定区间