第七章常微分方程数值解法.

1、1计算方法计算方法 吴筑筑编吴筑筑编2本章主要内容：本章主要内容：3引言：引言：可求出方程可求出方程y=1+ex的通解为的通解为 y=x+ex+c, 将初值条件将初值条件 x=0, y=2 代入得代入得 2=1+c, 故故 c=1,所以初值问题的解为所以初值问题的解为 y=x+ex+1 dxeyx)1 (求解初值问题求解初值问题 1(0)2xyey=+= 4引言：引言： 00)(,),(yxybaxyxfy本章解决的问题：本章解决的问题：一阶常微分方程的初值问题一阶常微分方程的初值问题5引言：引言：l若方程若方程 y=f(x,y)的右端函数的右端函数f(x,y)在闭矩形域在闭矩形域R：x0ax

2、 x0a， y0by y0b上满足：上满足：（1） f(x,y)在在R上连续上连续,（2）在）在R上关于上关于y满足满足Lipschitz(李普希兹李普希兹)条件条件即存在常数即存在常数L，对，对R上任意点均有以下不等式成立：上任意点均有以下不等式成立： |f(x,y1 )f(x,y2 )|L|y1y2|, xa,b, y1,y2R 则上述初值问题存在唯一的连续可微的则上述初值问题存在唯一的连续可微的解函数解函数 y = y(x)。6引言：引言：又如初值问题又如初值问题 12(0)0yxyy =-= 可求出它的解为可求出它的解为 220 xxtyee dt-=但要进一步计算指定点的函数值，还需

3、要用数值但要进一步计算指定点的函数值，还需要用数值积分方法。有些微分方程的解是隐函数，例如积分方法。有些微分方程的解是隐函数，例如要求函数值还需要解超越方程。应用中所处理的要求函数值还需要解超越方程。应用中所处理的微分方程往往很复杂且大都得不出一般解，所以微分方程往往很复杂且大都得不出一般解，所以一般用数值解法。一般用数值解法。1x yyex+=7引言：引言：数值解法：数值解法： 给定节点给定节点a=x0 x1xn=b, 将初值问题离散将初值问题离散化为差分方程化为差分方程,求出解函数求出解函数 y(x) 在这些点的近似在这些点的近似值值y1 ,y2 ,yn 。 所求得的近似值称为所求得的近似

4、值称为数值解数值解。 本章中总假定步长本章中总假定步长h为定值，节点为定值，节点xi=x0+ihi=1，2，38l 7.1.1 欧拉法及其截断误差欧拉法及其截断误差l 7.1.2 改进的欧拉法及预测改进的欧拉法及预测-校正公式校正公式 7.1 欧拉法和改进的欧拉法欧拉法和改进的欧拉法97.1.1 欧拉法及其截断误差欧拉法及其截断误差 00)(),(yxyyxfy初值问题：初值问题：1、欧拉公式的构造思想：用差商代替导数、欧拉公式的构造思想：用差商代替导数 012,nx x xx鬃设设 等距，步长为等距，步长为1,0,1 , ,1iixxhin+- =鬃 ()( )( )( )()()( )(

5、),y xhy xy xy xhy xh yxy xh f x yh+- +=+令令x=xi , x+h=xi+1 , y(xi )yi ，y(xi+1 ) yi+1 ,初值问初值问题离散化为题离散化为 )(, 2 , 1 , 0,),(001xyyiyxfhyyiiii(欧拉公式欧拉公式) 107.1.1 欧拉法及其截断误差欧拉法及其截断误差例例 取步长取步长 h=0.1，用欧拉法求解初值问题用欧拉法求解初值问题 解解:00( , ),0,1f x yxyxy=+=(0)1yxyy=+= y1=y0+h f(x0 ,y0 )=1+0.1(0 + 1 )=1.1y2=y1+h f(x1 ,y1

6、 )=1.1 + 0.1(0.1 + 1.1 )=1.22y3=y2+hf(x2 ,y2 )=1.22+0.1 (0.2+1.22)=1.362y10=y9+h f(x9 ,y9 )=y9+0.1(x9 + y9 )=3.18748 )(, 2 , 1 , 0,),(001xyyiyxfhyyiiii117.1.1 欧拉法及其截断误差欧拉法及其截断误差 )(, 2 , 1 , 0,),(001xyyiyxfhyyiiii2、欧拉公式几何意义、欧拉公式几何意义:用折线代替曲线计算解函数的近似值。用折线代替曲线计算解函数的近似值。 1100( ,)(),0,1,2,()iiiiiiyyf x yx

7、xiyy x127.1.1 欧拉法及其截断误差欧拉法及其截断误差3、数值公式的误差来源。、数值公式的误差来源。 (1)局部截断误差局部截断误差（简称（简称截断误差截断误差）：假设）：假设 yi=y(xi )是准确的是准确的 ，计算，计算yi+1所产生的误差所产生的误差 y(xi+1 ) - yi+1 若局部截断误差可以表示为若局部截断误差可以表示为O(hk+1), k为正整数，为正整数，则称公式是则称公式是k阶公式阶公式。（2）由于实际上）由于实际上yi不是准确值，因此它的误差会不是准确值，因此它的误差会传播下去。传播下去。 （3）实际计算时，每一步都可能产生舍入误差。）实际计算时，每一步都可

8、能产生舍入误差。 137.1.1 欧拉法及其截断误差欧拉法及其截断误差4、欧拉公式的截断误差是、欧拉公式的截断误差是O(h2)，公式是公式是1 阶的。阶的。 因为因为1( ,)( )( )iiiiiiyyh f x yy xh yx+=+211()( ) ( ) ( )2iiiy xy xyx hyh（泰勒公式）（泰勒公式） 两式相减，由设两式相减，由设 yi=y(xi ) ，有，有 221112iiy xyyhO h2( )11( )( ) ( ) ( )() ( )()2!nniiiiy xy xyxyxxyxxn147.1.2 改进的欧拉法及预测改进的欧拉法及预测-校正公式校正公式对微分

9、方程对微分方程y=f(x,y) 两边求两边求xi 到到xi+1 的定积分，有的定积分，有11()( )( , ( )iixiixy xy xf t y t dt+-=利用梯形公式计算积分，有利用梯形公式计算积分，有 1111()( ) ( , ( )(, ()2iiiiiiiixxy xy xf x y xf xy x+-+1、改进的欧拉公式的构造、改进的欧拉公式的构造157.1.2 改进的欧拉法及预测改进的欧拉法及预测-校正公式校正公式将将y(xi ) 、y(xi+1 )分别用分别用yi 、yi+1 代替，构造相应的代替，构造相应的数值公式：数值公式：（改进的欧拉公式）（改进的欧拉公式） )

10、(, 2 , 1 , 0,),(),(200111xyyiyxfyxfhyyiiiiii167.1.2 改进的欧拉法及预测改进的欧拉法及预测-校正公式校正公式2、改进的欧拉公式的截断误差为、改进的欧拉公式的截断误差为O(h3)因而改进的欧拉法是二阶的。因而改进的欧拉法是二阶的。 177.1.2 改进的欧拉法及预测改进的欧拉法及预测-校正公式校正公式3、 改进的欧拉法的具体使用格式。改进的欧拉法的具体使用格式。 改进的欧拉法是隐式公式改进的欧拉法是隐式公式 ,计算时常用迭代法。一计算时常用迭代法。一般每一步先由欧拉公式计算出般每一步先由欧拉公式计算出yi+1 的初始值的初始值yi+1(0),再迭

11、代计算再迭代计算yi+1。 , 2 , 1 , 0),(),(2),()(11)1(1)0(1kyxfyxfhyyyxfhyykiiiiikiiiii当满足当满足 (1)( )11|kkiiyye+-时，取时，取 (1)11.kiiyy+=可证明当可证明当f(x,y)满足一定条件时，迭代是收敛的。满足一定条件时，迭代是收敛的。 )(, 2 , 1 , 0,),(),(200111xyyiyxfyxfhyyiiiiii187.1.2 改进的欧拉法及预测改进的欧拉法及预测-校正公式校正公式改进的欧拉法的预测校正公式改进的欧拉法的预测校正公式 ( )1( )111( ,) ( ,)(,)20,1,2

12、,piiiipiiiiiiyyh f x yhyyf x yf xyi可证明预测校正公式的截断误差也为可证明预测校正公式的截断误差也为 O(h3)。197.1.2 改进的欧拉法及预测改进的欧拉法及预测-校正公式校正公式例例 取步长取步长h=0.2,用改进的欧拉法的预测校正公用改进的欧拉法的预测校正公式求解初值问题的数值解式求解初值问题的数值解y1 , y2 .(0)1yxyy =+= 解解 00( , ),0,1f x yxy xy=+=( )1( )1110.2()0.21.20.1()()piiiiiipiiiiiiyyxyxyyyxyxy+=+=+=+预测预测- -校正公式具体是校正公式

13、具体是 207.1.2 改进的欧拉法及预测改进的欧拉法及预测-校正公式校正公式()1000.21.21.2pyxy=+=( )1000110.1()()10.1(010.21.2)1.24pyyxyxy=+=( )2110.21.20.20.21.2 1.241.528pyxy=+=+=11.24y =21.5768y =( )2111220.1()()1.240.1(0.21.240.41.528)1.5768pyyxyxy=+=+=21设：改用后差商设：改用后差商 替代方程中的导数项替代方程中的导数项 7.1.2 向后向后 (隐式隐式)欧拉公式欧拉公式 111()iiiy xy xy xh

14、111,iiiiyyhf xy可以得到向后欧拉公式可以得到向后欧拉公式l这是隐式欧拉格式，也是一阶方法，精度与欧拉这是隐式欧拉格式，也是一阶方法，精度与欧拉公式相当。计算公式相当。计算yi+1通常用迭代法通常用迭代法:(0)1(1)( )111,0,1,2,iiiikkiiiiyyhf x yyyhf xyk227.1.2 两步欧拉公式两步欧拉公式设改用中心差商设改用中心差商 替代方程替代方程 中的导数项中的导数项 ，再离散化，即可，再离散化，即可导出下列格式导出下列格式无论是显式欧拉公式还是隐式欧拉公式，它们都无论是显式欧拉公式还是隐式欧拉公式，它们都是单步法，其特点是计算时只用到前一步的信

15、息是单步法，其特点是计算时只用到前一步的信息yi，而该格式却调用了前面两步的信息而该格式却调用了前面两步的信息yi-1,yi ，两步欧，两步欧拉格式因此而得名。拉格式因此而得名。两步欧拉格式具有更高的精度，它是二阶方法。两步欧拉格式具有更高的精度，它是二阶方法。1112iiy xy xh,nnnyxf xy x112,iiiiyyhfx y23引言：引言：（回顾）（回顾） 00)(,),(yxybaxyxfy本章解决的问题：本章解决的问题：一阶常微分方程的初值问题一阶常微分方程的初值问题247.1 欧拉法和改进的欧拉法欧拉法和改进的欧拉法预测校正公式预测校正公式 ()1()111( ,) (

16、,)(,)20,1,2,piiiipiiiiiiyyh f x yhyyf x yf xyi改进的欧拉公式改进的欧拉公式 )(, 2 , 1 , 0,),(),(200111xyyiyxfyxfhyyiiiiii )(, 2 , 1 , 0,),(001xyyiyxfhyyiiii欧拉公式欧拉公式 257.1 欧拉法和改进的欧拉法欧拉法和改进的欧拉法两步欧拉公式两步欧拉公式112,nnnnyyhfxy向后欧拉公式向后欧拉公式111,nnnnyyhfxy267.2 龙格龙格-库塔法库塔法(R-K法法)l7.2.1 二阶龙格二阶龙格-库塔公式库塔公式l7.2.2 四阶龙格四阶龙格-库塔公式库塔公式

17、27引言：引言：公式构造思想公式构造思想：从泰勒公式出发，寻找更高阶的：从泰勒公式出发，寻找更高阶的数值公式。数值公式。 例如，泰勒公式计算到二阶可得例如，泰勒公式计算到二阶可得 231()( )( )( )()2!y xhy xy x hyx hO h+=+( )( , ( )( , ( )( )( , ( )( , ( )( , ( )xyy xf x y xdf x y xyxfx y xfx y xf x y xdx = =+令令ixx=则则1ixhx+=略去余项，得出一个二阶的数值公式为略去余项，得出一个二阶的数值公式为21( ,)( ,)( ,)( ,)2iiiixiiyiiiih

18、yyf x y hfx yfx yf x y+=+因因28引言：引言：理论上按此方式可以得到更高阶的公式。但需理论上按此方式可以得到更高阶的公式。但需要计算复合函数的高阶导数，使算法复杂而不实要计算复合函数的高阶导数，使算法复杂而不实用。用。 龙格龙格库塔的思想（间接地运用泰勒公式）：库塔的思想（间接地运用泰勒公式）：利用利用y(x)在若干个点上的函数值和导数值，作出一在若干个点上的函数值和导数值，作出一个适当的线性组合，使这个线性组合按个适当的线性组合，使这个线性组合按h展开后的展开后的泰勒公式与泰勒公式与y(x+h)的泰勒公式有较多的项达到一致，的泰勒公式有较多的项达到一致，从而得出较高阶

19、的数值公式。从而得出较高阶的数值公式。 29引言：引言：R阶龙格阶龙格库塔（库塔（ Runge-Kutta)法的一般形式：法的一般形式：11122111()( ,)(,),(2,3, ),(1,2, ),(2,3, ),(2,3, ,2,3,1)iirriissisisjjjsssjyyh cKc Kc KKf x yKf xa h yhb Ksrcsr asrb sr js其中都是待定系数307.2.1 二阶龙格二阶龙格-库塔公式库塔公式设想一个有二阶精度的数值公式形状为设想一个有二阶精度的数值公式形状为 a, b为待定系数。为待定系数。1(, ()iiiiyybhf xah y xah+=

20、+仍令仍令x=xi，则，则x+h=xi+1。如果能找出。如果能找出a,b,使得使得 3()( ) ()()y xhy xbh yxahO h略去余项就可得到上面所希望的近似计算公式了。略去余项就可得到上面所希望的近似计算公式了。 因此考虑因此考虑 )( )()(haxyhbxyhxyhT 在在h=0处求泰勒公式得处求泰勒公式得 ( )23(0)(0)(0)()2TT hTThhO h=+由于由于(0)0T=317.2.1 二阶龙格二阶龙格-库塔公式库塔公式( )2()()()()Thyxhabyxaha bhyxahabyxah =+-+-+-+( )()()()Thy xhby xahabh

21、yxah=+-+-+ )( )()(haxyhbxyhxyhT (0)( )( )Ty xby x =-(0)(12)( )Tab yx=-由由T(h)的泰勒公式的泰勒公式( )23(0)(0)(0)()2TT hTThhO h=+327.2.1 二阶龙格二阶龙格-库塔公式库塔公式( )231(1)( )()( )()2T hb yx hab yx hO h=-+-+为使为使T(h)=O(h3),令令 101 20bab-=-=，解出解出 1 2,1ab=，得，得 32,2hOhxyhxfhxyhxyhT 整理得整理得337.2.1 二阶龙格二阶龙格-库塔公式库塔公式利用利用( )( )( )

22、2,2222hhhhf xy xyxyxyxO h骣骣骣琪珑+=+=+鼢珑鼢珑桫桫桫3()( )(,( )( , ( )()22hhy xhy xh f xy xf x y xO h+=+可以推出可以推出 取取x=xi并略去并略去O(h3)便得到二阶龙格库塔公式便得到二阶龙格库塔公式 1(,( ,)22iiiiiihhyyh f xyf x y+=+或或12112( ,),(,)22iiiiiihhkf x ykf xykyyhk+=+=+347.2.1 二阶龙格二阶龙格-库塔公式库塔公式3()( )(,( )( , ( )()22hhy xhy xh f xy xf x y xO h+=+可

23、以推出可以推出 取取x=xi并略去并略去O(h3)便得到二阶龙格库塔公式便得到二阶龙格库塔公式 1(,( ,)22iiiiiihhyyh f xyf x y+=+或或12112( ,),(,)22iiiiiihhkf x ykf xykyyhk+=+=+ 32,2hOhxyhxfhxyhxyhT 357.2.2 四阶龙格四阶龙格-库塔公式库塔公式仿照上述的讨论，可导出仿照上述的讨论，可导出四阶龙格库塔公式四阶龙格库塔公式: : 121324311234( ,),(,)22(,),(,)22(22)6iiiiiiiiiihhkf x ykf xykhhkf xykkf xh yhkhyykkkk

24、+=+=+=+=+ 36例例 取步长取步长h=0.2，用四阶龙格库塔公式求下面，用四阶龙格库塔公式求下面初值问题的数值解。初值问题的数值解。 7.2.2 四阶龙格四阶龙格-库塔公式库塔公式 2,(01)(0)1yyx yxy =-= 21(,)22iihhkf xyk=+解解 00( , )2,0,1,0.2f x yyx yxyh=-=由公式得由公式得 2001(2,2)(0.1,1.1)0.91818kf xhyhkf=+=100(,)1kf xy=1(,)iikf xy=377.2.2 四阶龙格四阶龙格-库塔公式库塔公式11234(22)6iihyykkkk+=+3002(2,2)(0.

25、1,1.09182)0.90864kf xhyhkf=+=4003(,)(0.2,1.18173)0.84324kf xh yhkf=+=18323. 16)22(2 . 0432101 kkkkyy32(,)22iihhkf xyk=+43(,)iikf xh yhk=+21yx=+数值解数值解yi与准确解与准确解y(xi)的对照见表的对照见表 准确解是准确解是00.20.40.60.81.011.183231.341671.483281.612511.7321411.183221.341641.483241.612451.73205xiyiy(xi)387.3 线性多步法线性多步法7.3.

26、1 四阶阿达姆斯四阶阿达姆斯(Adams)外插公式外插公式7.3.2 四阶阿达姆斯四阶阿达姆斯(Adams)内插公式内插公式7.3.0 多步法的概念多步法的概念7.3.3 初始出发值的计算初始出发值的计算7.3.4 阿达姆斯预测阿达姆斯预测-校正公式校正公式397.3.0 多步法的概念多步法的概念单步法单步法-计算计算yi+1时只使用时只使用yi的值。的值。多步法多步法-计算计算yi+1时使用前面的时使用前面的k个个y值，即由值，即由yi-k+1 , yi-k+2 , ,yi-1, yi计算计算yi+1。(k=1,2,)线性多步法线性多步法-计算计算yi+1的公式由的公式由yi-k+1 , y

27、i-k+2 , ,yi-1, yi 的线的线性组合表达。性组合表达。407.3.1 四阶阿达姆斯外插公式四阶阿达姆斯外插公式设想用设想用yi-3 ,yi-2 ,yi-1 ,yi 的值计算的值计算yi+1。为方便讨论由。为方便讨论由 (3 ), (2 ), (), ( )y xhy xhy xhy x-出发计算出发计算y(x+h), 由初值问题的方程由初值问题的方程y=f(x,y(x)两两边从边从x到到x+h积分积分, 可得到等价的积分方程可得到等价的积分方程 hxxdssysfxyhxy)(,()()(417.3.1 四阶阿达姆斯外插公式四阶阿达姆斯外插公式设想运用数值积分方法，取设想运用数值

28、积分方法，取 x-3h, x-2h, x- h, x 为插值基点，做为插值基点，做 f(s,y(s) 的三次拉格朗日插值，的三次拉格朗日插值，用它近似计算上式的积分用它近似计算上式的积分 。这样得到的数值积分公式是这样得到的数值积分公式是f(s,y(s)在在4个插值基个插值基点处的函数值的线性组合。点处的函数值的线性组合。 3)2()()()()(3210hxybhxybhxybxybhxyhxy 由于由于 f(x-ih,y(x-ih)= y(x-ih)，所得到的计算所得到的计算 y(x+h) 的近似公式形为：的近似公式形为：427.3.1 四阶阿达姆斯外插公式四阶阿达姆斯外插公式( )(4)

29、2345(0)(0)(0)(0)(0)()2!3!4!TTTT hTThhhhO h=+( )(0)ThT=LL为达到四阶精度，希望确定参数为达到四阶精度，希望确定参数b0 ,b1 ,b2 ,b3使满足使满足 )(3)2()()()()(53210hOhxybhxybhxybxybhxyhxy 运用在运用在 h=0 处的泰勒公式得处的泰勒公式得 0123( )()( )( ) () (2 ) (3 )T hy xhy xh b yxb yxhb yxhb yxh(0)0T=( )(0)T hT =LL( )(0)ThT =LL(4)(4)( )(0)ThT=LL437.3.1 四阶阿达姆斯外插

30、公式四阶阿达姆斯外插公式为达到四阶精度，希望确定参数为达到四阶精度，希望确定参数b0 ,b1 ,b2 ,b3使满足使满足 )(3)2()()()()(53210hOhxybhxybhxybxybhxyhxy 运用在运用在 h=0 处的泰勒公式得处的泰勒公式得 3(4)4533112299411(2) ( )()( )()62224632bbbbbbyx hyx hO h+-+ )3( )2( )( )()()(3210hxybhxybhxybxybhxyhxy 201231231(1) ( )(23 ) ( )2bbbbyx hbbbyx h=-+447.3.1 四阶阿达姆斯外插公式四阶阿达姆

31、斯外插公式0123123331122110,23029941120 ,062224632bbbbbbbbbbbbb-=+=-=+= 01235559373,2424248bbbb= -= -)( )3( 9)2( 37)( 59)( 5524)()(5hOhxyhxyhxyxyhxyhxy 为使误差等于为使误差等于O(h5)，令令h, h2,h3,h4 的系数为的系数为0，得方程组：得方程组：求得求得代入前面的公式得代入前面的公式得457.3.1 四阶阿达姆斯外插公式四阶阿达姆斯外插公式令令x=xi 并记并记 2222 (2 )(2 , (2 )(2 ,)(,)iiiiiiiiyxhf xh

32、y xhf xh yf xyf-=-= ( )( , ( )( ,)iiiiiiyxf x y xf x yf=1111 ()(, ()(,)(,)iiiiiiiiyxhf xh y xhf xh yf xyf-=-=3333 (3 )(3 , (3 )(3 ,)(,)iiiiiiiiyxhf xh y xhf xh yf xyf-=-=1( ),(),iiiy xyy xhy+467.3.1 四阶阿达姆斯外插公式四阶阿达姆斯外插公式略去余项，得到略去余项，得到四阶阿达姆斯外插公式四阶阿达姆斯外插公式 ：1123(5559379)24iiiiiihyyffff+-=+-+-这是显式公式。公式的

33、截断误差为这是显式公式。公式的截断误差为O(h5)。 477.3.2 四阶阿达姆斯内插公式四阶阿达姆斯内插公式把四阶阿达姆斯外插公式中使用的把四阶阿达姆斯外插公式中使用的yi-3 ,yi-2 ,yi-1 ,yi改为改为 ：yi-2 ,yi-1 ,yi ,yi+1 ,经类似的推导可得近似公式：经类似的推导可得近似公式： )( )2( )( )( )()(3210hxybhxybhxybxybhxyhxy 0123123331122110,20241120,062224636bbbbbbbbbbbbb-=+-=-=+-= 确定待定系数的方程组为确定待定系数的方程组为487.3.2 四阶阿达姆斯内插公式四阶阿达姆斯内插公式012319 24,5 24,1 24,3 8bbbb= -=解得解得得到四阶阿达姆斯内插公式得到四阶阿达姆斯内插公式 1112(9195)24iiiiiihyyffff+-=+-+这是一个隐式公式，截断误差也是这是一个隐式公式，截断误差也是O(h5)。 497.3.3 初始出发值的计算初始出发值的计算 阿达姆斯公式的特点是计算公式简单，只需简单阿达姆斯公式的特点是计算公式简单，只需简单的算术运算，计算量少，结