北京交通大学(数字分析研究生课程)6微分方程数值解

120第6章常微分方程初值问题数值解法6.1问题的描述和基本概念1、常微分方程初值问题一般形式0(,)()yfxyyay式中(,)fxy已知，0()yay称为初值条件.初值问题的数值方法和数值解求函数()yyx在若干离散点kx上的近似值(0,1,)kyk的方法称为初值问题的数值方法，而称(0,1,)kyk为初值问题的数值解.1212.建立数值解法的思想与方法微分方程初值问题的数值解法是用离散化方法将初值问题化为差分方程后再求解的方式.设节点为011nnaxxxx距离1kkkhxx称为步长.求数值解一般是从0y开使逐次顺序求出12,,yy.初值问题的解法有单步法和多步法两种：单步法：计算1ky时只用到ky一个值；多步法：计算1ky时要用1,,,kkklyyy多个值。数值解法还有显格式和隐格式之分。122微分方程离散化方法主要有数值微分法，数值积分法和Taylor展开法1)数值微分法由'()(,())kkkyxfxyx，用数值微分的2点前差公式代替'()kyx，得近似离散化方程11()()'()(,())kkkkkkkyxyxyxfxyxxx记1kkhxx，做()kkyyx，“”，得差分方程1(,)kkkkyyfxyh即1(,)kkkkyyhfxy（Euler公式）由初值条件0()yya及Euler公式可求出数值解12,,,,nyyy.Euler公式是显式单步法.1232）数值积分法在1[,]kkxx上对'(,)yfxy两边取定积分，得111()()'(,())kkkkxxkkxxyxyxydxfxyxdx右端积分用左矩形公式（数值积分公式）得1()()(,())kkkkyxyxhfxyx于是得到求初值问题的Euler方法1(,)kkkkyyhfxy124右端积分用右矩形公式（数值积分公式）得111()()(,())kkkkyxyxhfxyx于是得到求初值问题的后退Euler方法1+1+1(,)kkkkyyhfxy后退Euler方法是隐式的.125右端积分用梯形公式（数值积分公式）得近似离散化方程：111()()[(,())(,())]2kkkkkkhyxyxfxyxfxyx于是得到求初值问题的梯形方法111[(,)(,)]2kkkkkkhyyfxyfxy该公式是隐式单步法.1263）Taylor展开法因为初值问题中函数(,)fxy是已知函数，由(,)yfxy，可以计算''y，'''y，…，于是有函数()yyx在kx处的Taylor展式212()()'()''()2!()(,())[(,())]2!kkkkkkkkxxhyxyxhyxyxhdyxhfxyxfxyxdx取上式右端前若干项，得近似离散化方程.例如取前两项有1()()(,())kkkkyxyxhfxyx于是又得到Euler公式：1(,)kkkkyyhfxy.1273.数值解法的误差、阶与绝对稳定性单步法数学描述为111(,,,,)kkkkkkyyhxxyyh显式：1(,,)kkkkyyhxyh其中(,,)xyh称为增量函数.128显式单步法的一些概念定义1称111()kkkeyxy为单步法在节点1kx的整体截断误差，而称11()()(,(),)kkkkkTyxyxhxyxh为在1kx点的局部截断误差。()kyx表示解()yx在kx的值，是准确值，没有误差；ky表示由数值解公式得出()kyx的近似值，是数值解，有截断误差.129局部截断误差1kT的理解假设在计算()kyx时没有误差（()kkyyx）下,计算出的1ky（1()(,(),)kkkkyyxhxyxh）与1()kyx的误差111kkkTyxy(计算一步的误差）.定义2如果数值解法的局部截断误差为11()PkTOh则称该方法具有p阶精度或该方法是p阶方法.方法的阶越高，方法越好.130局部截断误差的主项如果某方法是p阶方法，11()PkTOh按h可展为1121()(,())()PPPkkkTOhgxyxhOh则称1(,())Pkkgxyxh为局部截断误差的主项.在同阶方法中，局部截断误差的主项越小，方法越好.对Euler方法1(,)kkkkyyhfxy，有1()()(,())kkkkkTyxhyxhfxyx将()kyxh在kx点展开，有2()()'()''()2!kkkkhyxhyxhyxyx2()(,())''()2!kkkkhyxhfxyxyx故有231()().2kkyxThOhEuler方法是一阶方法.131例1试求梯形方法的阶和局部截断误差主项.解该单步公式的局部截断误差是111()()((,())(,()))2kkkkkkkhTyxhyxfxyxfxyx1()()(()())2kkkkhyxhyxyxyx22()()()(()()23!2kkkkkhhhyxhyxyxyxyx24()())()2kkhyxhyxOh34()().12khyxOh故局部截断误差主项是3()12khyx，方法是二阶的.132定义3设某种数值方法在ky上大小为的扰动，于以后各()nynk上产生的偏差均不超过，则称该数值方法是稳定的。通常用试验方程'yy（为复数）来讨论求解数值方法绝对稳定性.Euler方法稳定性将Euler公式用于试验方程'yy，得到1(1)kkkkyyhyhy设计算ky时有误差,k则有11(1)()kkkkyhy得1(1)kkh要想1kk，只须11h，因此Euler方法在11h时是绝对稳定的，其绝对稳定域为复平面h上以-1为中心的单位圆盘.绝对稳定区间为20.h1336.2Runge-Kutta方法11111,,(2,3,,)mkkiiikkrrkrkrjjjyyhcKKfxyKfxahyhbKrm称为m级R-K方法.增量函数是1,,miiixyhcK（x,y,h）134构造过程以2m来说明Runge-Kutta方法的构造方法和过程，对一般的Runge-Kutta方法可类似处理.2m的Runge-Kutta公式为11122kkyyhcKcK式中1,kkKfxy，22211,kkKfxahyhbK.引进符号,fxyxf，,xxfxyxf，,()yyfxyxf，…由,yfxy，可得xyyxfff，…在kx处做Taylor展开，有23232!2!kkkkkxyyxyxhyxyxhhOhhyxhffffOh对,,kkxyxh在(,())kkxyx做二元Taylor展开，有212221,,kkxyxyxhcfcfahfhbffOh21222221xyccfcafcbffhOh由1,,kkkkkTyxhyxhxyxh，有231122222111122kxyTccfhcafcbffhOh135选1222221111,0,022cccacb有局部截断误差3TOh，这样可得到二阶Runge-Kutta公式.取20ct，则式（6.13）的解为11ct，22112abt取不同的t可得出不同的二阶Runge-Kutta公式.如取12t时，得到改进的Euler公式1,,,2kkkkkkkkhyyfxyfxhyhfxy1t时，得到中点公式1(,(,)).22kkkkkkhhyyhfxyfxy136经典Runge-Kutta公式112341213243226,,22,22,kkkkkkkkkkhyyKKKKKfxyhhKfxyKhhKfxyKKfxhyhK四阶方法.137例1设初值问题为100yyy00.4x分别用Euler方法（0.025h），改进Euler方法（0.05h）和经典Runge-Kutta方法（0.1h）计算。解Euler方法计算格式（0.025h）为10.0251kkkyyy改进的Euler方法计算格式（0.05h）为1121210.02510.051kkkkyyKKKyKyK138经典Runge-Kutta方法计算格式（0.1h）为1123412132431(22)6010.0510.0510.11kkkkkkyyKKKKKyKyKKyKKyK它们的初值00y，计算结果及准确解列于下表kxEuler方法改进Euler方法经典R-K法)(kxy000000.10.0963120.0951230.095162500.095162580.20.1833480.1811930.181269100.181269250.30.2620010.2590850.259181580.259181780.40.3330790.3295630.329679710.32967995139例2给定初值问题0(,)()yfxyaxbyay1）分析求解公式122[(,)3(,(,))]433mmmmmmmmhyyfxyfxhyhfxy的局部截断误差，指出它是几阶公式；2）证明用上面求解公式计算初值问题'01(0)1yyxy的数值解1nmmy成立极限2011lim6nhyyhe本题中的节点是等距节点，h为步长，n为由节点分割区间[a,b]的份数.140解由题意有,,0,1,,mbahxamhmnn1）局部截断误差1122()()[(,())3(,()(,()))]433mmmmmmmmmhTyxyxfxyxfxhyxhfxyx1322()()()(,()())4433mmmmmmhyxyxyxhfxhyxhyx将1mmyxyxh在mx点做Taylor展开到3h项，将22(,()())33mmmfxhyxhyx在,mmxyx点做二元Taylor展开到2h项，则有3431,()()6mmymmhTyxfxyxOhOh得所用公式是2阶的.2）显然所给初值问题的准确解为xyxe。由给出的数值解计算公式有2111111222021[3()]1432111111222nnnnnnnnhhyyyyhyhhyhhyhhhh故1122200111limlim12hnhhyyehhhh211ln112201limhhhheeh141231116201limhOhheeh2311620lim1hOhheeh1234201lim6hehOhOhh16e式中221111ln1ln1