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.数值解法的误差、阶与绝对稳定性单步法的一般形式可以表示成11(,,,)kkkkkyyhxyyh显式：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方法是一阶方法.例1试求梯形方法的阶和局部截断误差主项.解该单步公式的局部截断误差是131111()()((,())(,()))2kkkkkkkhTyxhyxfxyxfxyx1()()(()())2kkkkhyxhyxyxyx23()()()(()()23!2kkkkkhhhyxhyxyxyxyx24()())()2kkhyxhyxOh34()().12khyxOh故局部截断误差主项是3()12khyx，方法是二阶的.定义3设某种数值方法在ky上大小为的扰动，于以后各()nynk上产生的偏差均不超过，则称该132数值方法是稳定的。通常用试验方程'yy（为复数）来讨论求解数值方法绝对稳定性.Euler方法稳定性将Euler公式用于试验方程'yy，得到1(1)kkkkyyhyhy设计算ky时有误差,k则有11(1)()kkkkyhy得1(1)kkh要想1kk，只须11h，因此Euler方法在11h时是绝对稳定的，其绝对稳定域为复平面h上以(-1,0)为中心的单位圆盘.绝对稳定区间为20.h6.2Runge-Kutta方法13311111,,(2,3,,)mkkiiikkrrkrkrjjjyyhcKKfxyKfxahyhbKrm称为m级R-K方法.增量函数是()1,,,,miiixyhcKxyh构造过程134以2m来说明Runge-Kutta方法的构造方法和过程，对一般的Runge-Kutta方法可类似处理.2m的Runge-Kutta公式为11122kkyyhcKcK式中1,kkKfxy，22211,kkKfxahyhbK.由,yfxy，可得(),(),(),()xyyxfxyxfxyxfxyx在kx处做Taylor展开，有2322!(,())(,())2!kkkkkkkxkkyxyxhyxyxhhOhhyxhfxyxfxyx3(,())(,())ykkkkfxyxfxyxOh对,,kkxyxh在(,())kkxyx做二元Taylor展开，有12,,(,())(,())kkkkkkxyxhcfxyxcfxyx21(,())(,())ykkkkhbfxyxfxyx22(,())()xkkahfxyxOh12()(,())kkccfxyx221(,())(,())ykkkkcbfxyxfxyx22(,())()xkkafxyxhOh由1,,kkkkkTyxhyxhxyxh，有1351122211(,())()(,())2kkkxkkTccfxyxhcafxyx232211()(,())(,())()2ykkkkcbfxyxfxyxhOh选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故1411122200111limlim12hnhhyyehhhh211ln112201limhhhheeh231116201limhOhheeh