计算方法第7章常微分方程

1、第七章第七章 常微分方程数值解常微分方程数值解本章研究一阶常微分方程初值问题的数值解。常微分方程分为线性,非线性两类,线性微分方程是非线性微分的特例.而高阶可化为一阶方程组.若将方程组中的未知量写成一个向量,那么本章所涉及的方程表示为( ,) , (7.1)( )dyf x yxa bdxy ah=近似方法有两类:一类是近似解析方法.如级数解法,逐次逼近法等;另一类就是数值解法,它可以给出解在一些离散点上的近似值.nxxxxx210nyyyyy210在使用数值法求解之前,需要考虑解的存在性和唯一性问题. 有如下定理:数值解数值解定理定理7.1 设f(x,y)在D=(x,y)|ax b,yR上有

2、定义且连续,同时满足如下李普希茨(Lipschitz)条件则对 a, b, R,初值问题(7.1)在a,b上存在唯一的连续可微解 y (x).(7.2)中L称为李普希茨常数.( , )( , ),( , ),( , ),0 (7.2)f x yf x yL yyx yD x yDL-挝 +%0 x0y求解之前,应检验是否满足Lipschitz条件.如果f(x,y)对y可导,且( ,)f x yy有界,则令(, )( ,)maxx yDf x yLy=则问题(7.2)一定满足Lipschitz条件.因为：( , )( ,)( ,)()f xf x yf x yyyL yyyx-=-要计算出数值解

3、要计算出数值解, ,实际上就是求对应于一系列已知节点实际上就是求对应于一系列已知节点 a a = = x x0 0 x x1 1 x xn n= = b b处的函数值处的函数值y y0 0, , y y1,1, , y yn n通常采用通常采用等距节点等距节点，即取即取 hi = h (常数常数)。记理论解为y(xn),数值解为yn, fn=f (xn, yn),()(, ()nnny xf xy x=01nxxxLLL01nyyyLLLixiy数值解数值解1(0,.,1)iiihxxin+=-=-节点间距节点间距称为步长，称为步长，7.2 7.2 欧拉欧拉(Euler)(Euler)方法方法

4、1.泰勒函数展开对y(x)在xn点泰勒展开,并令x=xn+1得:21n1()()()(), (,) (7.3)2!nnnnnnhy xy xhy xyxxxx+=+略去误差项略去误差项,有有假设y(x)充分连续可微.下面是导出欧拉(Euler)方法的三种途径.2. 数值微分由导数的定义知:当h充分小时有:1()()(, () (7.5)nnnny xy xhf xy x+1()()()(, (). (7.7)nnnnny xy xy xf xy xh+-=即同样可得公式(7.5):1()()(, ()nnnny xy xhf xy x+3 数值积分在xn, xn+1上对 积分得( , ( )y

5、f x y x=11()()( ,( )nnxnnxy xy xf x y xdx+=+利用数值积分的左矩形公式得1( ,( )(,()nnxnnxf x y xdxhf xy x+即同样可得公式(7.5):1()()(, ()nnnny xy xhf xy x+nx1nx+在公式(7.5)中分别用yn, yn+1近似代替y(xn),y(xn+1)得:10(,) n=1,2, (7.6)nnnnyyhf xyyh+=+= L(7.6)称为欧拉折线法,简称欧拉(Euler)法.欧拉法欧拉法例1（p198）：用欧拉法求解初值问题2(0)1xyyyy =-= 在0,1上的数值解,取h=0.1,并与精

6、确解比较.精确解为:12yx=+因为2( ,)xf x yyy=-由欧拉法求解公式得:102()1, n=1,2,10nnnnnxyyh yyy+=+-= L解解:计算结果见计算结果见p1987.2.2 隐式公式及改进的欧拉方法隐式公式及改进的欧拉方法利用数值积分的右矩形公式得111( , ( )(, ()nnxnnxf x y x dxhf xy x+进一步得1110(,) n=1,2, (7.10)nnnnyyhf xyyh+=+= L这是一个隐式方法.nx1nx+称为称为隐式欧拉法隐式欧拉法利用数值积分的梯形公式得得:111( , ( ) (, ()(, ()2nnxnnnnxhf x

7、y x dxf xy xf xy x+1110 (,)(,)2 n=1,2, (7.11)nnnnnnhyyf xyf xyyh+=+=L这也是一个隐式方法,精度比欧拉法好，但不便计算。一般用下式计算：通常采用的算法:11011(,), (7.13) (,)(,),2,0,1,2,nnnnnnnnnnyhf xyhyyf xyf xynyyh+=+=+=L称作称作预估预估-校正法校正法,也称也称改进的欧拉法改进的欧拉法改进的欧拉法也可以写成:例2(p200) : 用预估预估-校正法校正法(7.14)求解例1中的初值问题.1121210(),2(,), (7.14)(,), n=0,1,2,.n

8、nnnnnhyyKKKf xyKf xh yhKyh+ =+=+=L=-= 2(0)1xyyyy与（与（7.13）等价）等价哟哟设则依公式(7.14)得:2( , )xf x yyy=-1211112022()()21, n=1,2,.nnnnnnnnxKyyxhKyhKyhKhyyKKy+=-+=+-+=+=L计算结果见计算结果见p201解解:定义：假设7.2.3 误差分析误差分析()nnyy x=则称误差为局部截断误差。11()nny xy+-1. 欧拉公式的局部截断误差局部截断误差分析 设yn=y(xn),则有此时欧拉公式可写成1()()nnnyy xhy x+=+对精确解y(x),由泰

9、勒公式得231()()()()()2nnnnhy xy xhy xyxO h+=+所以欧拉公式的局部截断误差为:23211()()()()2nnnhy xyyxO hO h+-=+=就是1ny+(*)()()(,)nnnnny xf x y xf xy=2.改进的欧拉公式的改进的欧拉公式的局部截断误差局部截断误差分析分析由于()(, ()nnny xf xy x=所以1(,)(, ()()nnnnnKf xyf xy xy x=()(, ()(, ()(, ()nxnnnnynnyxfxy xf xy xfxy x=+211212(,)(, () = (, ()(, ()(, ()() = (

10、)()()nnnnnnxnnynnnnKf xh yhKf xh y xhKf xy xhfxy xhK fxy xO hy xhyxO h=+=+将K1, K2代入公式(7.14)中得231()()()()2nnnnhyy xhy xyxO h+=+，两边对x求导： (*)(*)- (*)得到得到改进的欧拉公式的局部截断误差为类似地,梯形公式的局部截断误差为311()()nny xyO h+-=34311()()()()12nnnhy xyyxO hO h+-= -+=如果单步差分公式的局部截断误差为O(hp+1),则称该公式为p阶方法,这里p为正整数.由此知欧拉法是一阶法,改进的欧拉法和梯

11、形法为二阶方法.7.3 龙格-库塔(Runge-Kutta)方法7.3.1龙格-库塔方法的构造 为了得到高阶的单步显式差分公式,我们可以对y(x)在xn处进行泰勒展开:2( )112(1)( )1()()()()()()2!= ()(,)(, ()(, ()()2!pppnnnnnpppnnnnnnnhhy xy xhy xyxyxO hphhy xhf xyfxy xfxy xO hp+=+LL于是可得差分公式2(1)( )10()=(,)(,)(,)2! n=1,2,.ppnnnnnnnnhhy xyhf xyfxyfxypyh+= LL其中:(1)( , )( , ( )( , )( ,

12、 ) ( , )xydfx yf x y xfx yfx y f x ydx=+2(2)22( , )( , ( )()().dfx yf x y xdxfffffxyyxy=抖抖=+抖抖由于涉及到高阶复合函数的导数,计算较困难.下面介绍如何化解此困难.以2阶方法为例设11122121()(,),(,)nnnnnnyyhKKKf xyKf xh yhKllab+=+=+ 对K2的右端在(xn,yn)处做一阶泰勒展开22( ,)( ,)( ,)( ,)()nnxnnnnynnKf x yh f x yh f x y f x yO hab=+从而21122232() (,)(,)(,) (,)()

13、nnnnxnnynnnnyyhf xyhfxyhfxyf xyO hlll al b+=+而y(xn+1)在xn处的泰勒展开式为23123()()()()()2()(,)(, ()2 (, () (,)() (7.18)nnnnnnnxnnynnnnhy xy xhy xyxO hhy xhf xyfxy xfxy xf xyO h+=+=+21122232()(,)(,)(,)(,)() (7.17)nnnnxnnynnnnyyhf xyhfxyhfxyf xyO hlll al b+=+令yn=y(xn),比较(7.16)和(7.18),得1222111, , 22llalbl+=这是一个

14、有三个方程四个未知数的方程组.有无穷多解,所以有无穷多个二阶龙格-库塔方法.令121, = =12llab=得改进的欧拉法.令1210,1, = =2llab=得中点公式,即12121(,)1(,)2nnnnnnyyhKKf xyKf xh yhK+=+=+令 得休恩公式，即 高阶龙格-库塔公式可作类似地推导。下面是常用的三阶，四阶龙格-库塔公式。12322,443llab=112121(3)4(,)22(,)33nnnnnnhyyKKKf xyKf xh yhK+=+=+ 三阶龙格-库塔公式1123121312(4)6(,)11(,)22(,2)nnnnnnnnhyyKKKKf xyKf xh yhKKf xh yhKhK+=+=+=+-+标准四阶龙格-库塔公式112341213243(22)6(,)11(,)2211(,)22(,)nnnnnnnnnnhyyKKKKKf xyKf xh yhKKf xh