chapt常微分方程數(shù)值解法

1、第第7章章 常微分方程及方程組的常微分方程及方程組的數(shù)值解法數(shù)值解法 常微分方程的定解問(wèn)題一般分為初值問(wèn)題和邊值問(wèn)題，而邊值問(wèn)題常?？梢曰癁槌踔祮?wèn)題來(lái)求解，所以本章只討論初值問(wèn)題。設(shè)在XY平面有：當(dāng) 求 。在工程問(wèn)題中，只有極少有解析解，大多數(shù)只能求數(shù)值解。即在區(qū)域：求y(xi)的近似解： 。 ( , )dyf x ydx(7-1)00()y xy( )yy x012,nax x xxb012,nyy yy7.1 常微分方程的數(shù)值法 本節(jié)研究常用的數(shù)值方法。先說(shuō)明幾個(gè)符號(hào)的含義：一、Taylor級(jí)數(shù)展開(kāi)法把yi+1在yi處展開(kāi)成Taylor級(jí)數(shù)，則有因?yàn)樗?O(h2)為h2的同階無(wú)窮小，略去

2、高階項(xiàng)得如下的近似關(guān)系： 111111(,)()iiiiiiiff xyyy xxxh21111()()2iiiiiiiiyyxxyxxy1( ,),iiiiiiyf x yfxxh 21()iiiyyhfO h 7.1 常微分方程的數(shù)值法 這就是常微分方程數(shù)值解法中的向前歐拉(Euler)公式。二、數(shù)值微分法這種方法就是利用差商近似代替導(dǎo)數(shù)，即100()iiiyyhfy xy 111211()( )( )( )2()( )()( )2()()( )( )26iiiiiiiiiiiiy xy xhy xyhy xy xhy xyhy xy xhy xyh(7-2)7.1 常微分方程的數(shù)值法略去

3、二階導(dǎo)數(shù)項(xiàng)，用 代替 ，得利用它們求問(wèn)題(7-1)式數(shù)值的方法，分別稱(chēng)為歐拉(Euler)法，后退歐拉法，中點(diǎn)法。三、數(shù)值積分法將(7-1)式在區(qū)間 內(nèi)積分得iy( )iy x11111100( ,)(,)2( ,)()iiiiiiiiiiiiyyhf x yyyhf xyyyhf x yy xy(7-3)1 ,iix x111d( ,)diiiixxiiiixxyyyf x yx7.1 常微分方程的數(shù)值法對(duì)上式右端的積分使用左端矩形積分公式，則有：于是(7-1)式的數(shù)值解法但如果上式右端的積分使用的是梯形積分公式，則有從而(7-1)式的數(shù)值解法也改為這就是改進(jìn)的歐拉公式。1( ,)diixi

4、iixf x yxhf100()iiiyyhfy xy (7-4)111( ,)d2iixiiiixf x yxh ff110012()iiiyyh ffy xy(7-5)7.1 常微分方程的數(shù)值法四、幾個(gè)基本概念（1）單步法：總是由前一點(diǎn)yi推算后一點(diǎn)yi+1，只需要一個(gè)點(diǎn)的值。除此之外，還有多步法，如Adams法。（2）顯示法：yi+1只出現(xiàn)等式一邊。若同時(shí)出現(xiàn)在右邊，則為隱式法。（3）截項(xiàng)誤差：當(dāng)h趨于零時(shí)，它是h2的同階無(wú)窮小，（7-2-1）誤差階數(shù)為2。如果y(x)是一次式，則O(h2)0,此時(shí)歐拉法是精確的（直線逼近直線），同時(shí)也說(shuō)明歐拉法是一階代數(shù)精度。 223()/2!/3!O

5、 hyhyh7.2 一階常微分方程的數(shù)值法讓我們研究(7-1)式形式的一階常微分方程的數(shù)值求解方法。7.2.1 歐拉折線法一、算法分析二、幾何意義參見(jiàn)圖7-1中的折線f01f11f21f31，歐拉折線法命名由此而來(lái)。容易看出，歐拉折線法的數(shù)值解離理論曲線越來(lái)越遠(yuǎn)。折線f01f11f21f31為后面要討論的歐拉二次法。10100(1)( ,)()iiiiiixxhxihyyhf x yyy x (7-6)7.2 一階常微分方程的數(shù)值法( )yf x01f11f21f31f02f12f22f32f1x2x3x4x5xxy圖7-1 歐拉折線法及改進(jìn)的歐拉法的幾何意義7.2 一階常微分方程的數(shù)值法三、

6、通用程序框圖與程序設(shè)計(jì)subroutine eulercomp f(x0,y0)x0=ax=x0+hdo i=1,n+1y=y0+h*f(x0,y0)end doend subroutine eulerx=x0;y0=yH=(b-a)/(n+1)7.2 一階常微分方程的數(shù)值法10 subroutine euler(a,b,n,y0,y)20 x0=a30 h=(b-a)/(n+1)40 do i=1,n+150 x=x0+h60 y=y0+h*f(x0,y0)70 x0=x;y0=y80 end do90 end subroutine euler7.2 一階常微分方程的數(shù)值法7.2.2 改進(jìn)的

7、歐拉折線法（梯形法）一、算法分析 在歐拉法中，泰勒級(jí)數(shù)是取前兩項(xiàng)。下面看看取三項(xiàng)時(shí)是什么樣子。 為了求yi+1，在yi處進(jìn)行泰勒級(jí)數(shù)展開(kāi)。依據(jù)前面所講的，可以寫(xiě)成：因?yàn)樗?23123/2()( ,)( ,)/2()iiiiiiiiiyyhyhyO hdf x yyhf x yhO hdx 11( ,)/( (,)()/( )iiiiiidf x ydxf xyf x yhO h111/2 ( ,)(,)iiiiiiyyhf x yf xy7.2 一階常微分方程的數(shù)值法這就是歐拉二次逼近法的遞推式，其局部誤差為O(h3)?？梢?jiàn)，比歐拉法精度高（隱式）。在上式中，等號(hào)的右邊也含有yi+1, 要出

8、yi+1一般要用迭代法。為了簡(jiǎn)便，用歐拉法先求出y*i+1，作為yi+1的預(yù)報(bào)值。于是，歐拉二次逼近法的算法可表示為：二、通用程序框圖與程序設(shè)計(jì)1*1*11100( ,) ( ,)(,)/2()iiiiiiiiiiiixxhyyhf x yyyhf x yf xyyy x(7-7)7.2 一階常微分方程的數(shù)值法subroutine meulercomp f(x0,y0)x0=a;y=y0 x=x0+hdo i=1,n+1y1=y0+h*f(x0,y0)end doend subroutine meulerx=x0;y0=yH=(b-a)/(n+1)y=y0+h*(f(x0,y0)+f(x,y1

9、)/2comp f(x,y1)7.2 一階常微分方程的數(shù)值法10 subroutine meuler(a,b,n,y0,y)20 x0=a;y=y030 h=(b-a)/(n+1)40 do i=1,n+150 x=x0+h60 y1=y0+h*f(x0,y0)70 y=y0+h*(f(x0,y0)+f(x,y1)/270 x0=x;y0=y80 end do90 end subroutine meuler7.2 一階常微分方程的數(shù)值法7.2.2 Runge-Kutta方法一、算法分析一階常微分初值問(wèn)題的求解，在小區(qū)間上歸結(jié)為如何選擇去逼近該區(qū)間曲線的直線。如果直線已經(jīng)確定，已知yi，就很容易

10、求出yi+1。依據(jù)Lagrange微分中值定理，有：其中，h=xi+1xi; 01。同時(shí)有 y(xi+h)=f(x,y)=f(xi+h,y(xi+h)，這就是曲線y=y(x)在xi+h處切線的斜率。如果已確定，則斜率也就確定。有了斜率和已知點(diǎn)（xi,yi）坐標(biāo)，就能確定直線： 1()( )()iiiy xy xhy xh1()(, (), ,iiiiiiyyxxf xh y xhx x7.2 一階常微分方程的數(shù)值法上式中的f為小區(qū)間曲線的平均斜率。剩下的問(wèn)題就是如何確定該平均斜率？斜率的不同選擇方法或算法，決定不同的計(jì)算格式。如：用起始點(diǎn)的斜率，則Ki=y(xi)是曲線在xi點(diǎn)的斜率，用此斜率

11、代替平均斜率。上式就是歐拉算法。如：用區(qū)間上兩端點(diǎn)的斜率的平均值，則其中，Ki=f(xi,yi)，Ki+1=f(xi+1,yi+1)，用兩點(diǎn)斜率的算術(shù)平均值代替平均斜率。此處的yi+1用歐拉法作預(yù)報(bào)值，則 1(, ()iiiiyyhf xh y xh 1( , ( )iiiiiiyyhf x y xyhK(7-8)1111/2 ( , ( )(, () /2 ()iiiiiiiiiyyhf x y xf xy xyhKK(7-9)7.2 一階常微分方程的數(shù)值法 這就是歐拉二次逼近法的算法式。（7-8）式比（7-9）式的計(jì)算精度要高，這說(shuō)明增加f(x,y)或斜率的計(jì)算次數(shù)可以提高截?cái)嗟碾A數(shù)，從而

12、提高了精度。這就是Runge-Kutta法的基本思想。將不同點(diǎn)的斜率進(jìn)行線性組合，其一般的形式為： 只要適當(dāng)選取1、2、a、b的值，就能保證局部截?cái)嗾`差為O(h3)。其中一個(gè)特例是1=1/2=2，a=b=1,此時(shí)就變成了（7-9）式，亦稱(chēng)二階Runge-Kutta法。四階Runge-Kutta法就是在小區(qū)間上計(jì)算四個(gè)斜率，爾后進(jìn)行線性組合，一般形式為： 11( , ( ),(,)iiiiiiiKf x y xKf xyhK11122121()( ,),(,)iiiiiiyyhKKKf x yKf xah ybhK7.2 一階常微分方程的數(shù)值法 適當(dāng)選擇待定系數(shù)，可得到截?cái)嗾`差O(h5)。四階龍

13、格庫(kù)塔法的常用算法為： 1112233441211132213243415263()( ,)(,)(,)(,)iiiiiiiiiiyyhKKKKKf x yKf xa h ybhKKf xa h yb hKb hKKf xa h yb hKb hKb hK 7.2 一階常微分方程的數(shù)值法二、通用程序框圖與程序設(shè)計(jì)112341213243(22)/6( ,)(/2,/2)(/2,/2)(,/2)iiiiiiiiiiyyhKKKKKf x yKf xhyhKKf xhyhKKf xh yhK (7-10)7.2 一階常微分方程的數(shù)值法subroutine rkuttax0=a;y=y0 x1=x0

14、+(i-1)*hdo i=1,nx=x1+hh(j)end doend subroutine ykuttayy=yy+ff*hh(j+1)/3h=(b-a)/nx=x1;y1=y;yy=yhh(1)=h/2;hh(2)=hh(1);hh(3)=h;hh(4)=h;hh(5)=h/2do j=1,4ff=f(x,y)comp f(x,y)y=y1+ff*hh(j)y=yyend do7.2 一階常微分方程的數(shù)值法10 subroutine rkutta(a,b,n,y0,yy)15 dimension hh(5)20 x0=a;y=y030 h=(b-a)/n40 hh(1)=h/2;hh(2)

15、=hh(1)50 hh(3)=h;hh(4)=h;hh(5)=h/260 do i=1,n70 x1=x0+(i-1)*h80 x=x1;y1=y;yy=y90 do j=1,4100 ff=f(x,y)110 x=x1+hh(j)120 yy=yy+ff*hh(j+1)/3130 y=y1+ff*hh(j)140 end do150 y=yy160 end do170 end subroutine rkutta7.3 一階常微分方程組的數(shù)值解法 我們所要解決的實(shí)際問(wèn)題，往往是由多個(gè)常微分方程組成的聯(lián)立方程組，所以本節(jié)我們討論一階常微分方程組的數(shù)值解法。 一階常微分方程組的一般形式為：一、算法

16、概述 我們可以把求單個(gè)常微分方程的四階Runge-Kutta法的遞推計(jì)算公式（7.10）式推廣到求解常微分方程組的情形，即： 123,00( ,)()jjjjnjjdyyfx y yyyydxyyx (7-11)7.3 一階常微分方程組的數(shù)值解法1,11,2,3,21,1,12,2,1,1,1,31,1,22,2,2,2,2,41,1( ,)(,)22222(,)22222(,iijjiiiij in ijjiiij ijn injjiiij ijn injiixxhkfx yyyyyhhhhhkfxykykykykhhhhhkfxykykykykkf xh yhk,32,2,3,3,3,1,

17、1,2,3,4,00,)1(23)6()ij ijn inj ij ijjjjjjyhkyhkyhkyyh kkkkyyx(7-12)7.3 一階常微分方程組的數(shù)值解法 在求 時(shí)，必須把所有的 （ ）都先求出來(lái)，因?yàn)槔ㄌ?hào)內(nèi)的通式為 ；同理，求 時(shí)，必須把所有的 都先求出來(lái)，因?yàn)槔ㄌ?hào)內(nèi)的通式為 ；求 時(shí)，必須把所有的 都先求出來(lái)，因?yàn)槔ㄌ?hào)內(nèi)的通式為 。二、通用程序框圖和程序設(shè)計(jì) ,1jk1,2,3,jn,12j ijhyk,2jk,3jk,2jk,22j ijhyk,4jk,3jk,32j ijhyksubroutine hrkuttax0=ax1=x0+(i-1)*hdo i=1,nx=x1+

18、hh(k)end doend subroutine hrkuttayy(l)=yy(l)+f(l)*hh(k+1)/3h=(b-a)/nx=x1hh(1)=h/2;hh(2)=hh(1);hh(3)=h;hh(4)=h;hh(5)=h/2do k=1,4call fun(x,y,f)y(l)=y1(l)+f(l)*hh(k)y(j)=yy(j)end dodo j=1,my1(j)=y(j);yy(j)=y(j)end dodo l=1,mend dodo j=1,mend doread y(i), i=1,m7.3 一階常微分方程組的數(shù)值解法10 subroutine orkutta(a,b

19、,n,m,yy)20 dimension hh(5),y(m),y1(m),f(m) 30 &yy(mm)40 x0=a50 do j=1,m60 read(5,f12.5) y(j)70 end do80 h=(b-a)/n90 hh(1)=h/2;hh(2)=hh(1)100 hh(3)=h;hh(4)=h;hh(5)=h/2110 do i=1,n120 x1=x0+(i-1)*h130 x=x1140 do j=1,m150 y1(j)=y(j);yy(j)=y(j)160 end do170 do k=1,4180 call fun(x,y,f)190 x=x1+hh(k)2

20、00 do l=1,m210 yy(l)=yy(l)+f(l)*hh(k+1)/3220 y(l)=y1(l)+f(l)*hh(k)230 end do240 end do250 do j=1,m260 y(j)=yy(j)270 end do280 end do290 end subroutine orkutta7.4 高階常微分方程（組）的數(shù)值解法 一、算法分析 關(guān)于高階常微分方程（組）初值問(wèn)題的數(shù)值解法，一般可以引入新的變?cè)阉鼈兓癁橐浑A常微分方程組的初值問(wèn)題用Runge-Kutta法來(lái)求解。一般地，對(duì)于初值問(wèn)題： 則（7.11）式可被化為等價(jià)的初值問(wèn)題： ( )( )(1)(1)(1

21、)(1)( , ,)(1,2,3, )nnjnniiiid yyf x y y yyydxdyysindx(7-11)(1)iiyy(1,2,3, )in7.4 高階常微分方程（組）的數(shù)值解法則可用前面介紹的任何一種方法求解（7.12）式。 1223111230( ,)()(1,2,3, )jjnnnjniiyyyyyyyyyf x y yyyyy xsin (7-12)7.4 高階常微分方程（組）的數(shù)值解法例 7-1 某振動(dòng)過(guò)程的動(dòng)態(tài)模型如下，求其在區(qū)間0，0.2內(nèi)的數(shù)值解，取h=0.2。解：這是一個(gè)二階常微分方程組的初值問(wèn)題，令 , 可將原方程組化為如下的等價(jià)形式： 11221212120

22、.15278.78110.5700.15110.57110.570(0)(0)0(0)(0)1yyyyyyyyyy13yy 24yy 7.4 高階常微分方程（組）的數(shù)值解法解：這是二、通用程序框圖與程序設(shè)計(jì)132431241212341858.533737.133737.133737.133(0)(0)0(0)(0)1yyyyyyyyyyyyyy subroutine hrkuttax0=a;y=y0 x1=x0+(i-1)*hdo i=1,nx=x1+hh(k)end doend subroutine hrkuttayy(l)=yy(l)+f(l)*hh(k+1)/3h=(b-a)/nx=x1hh(1)=h/2;hh(2)=hh(1);hh(3)=h;hh(4)=h;hh(5)=h/2do k=1,4call fun(x,y,f)y(l)=y1(l)+f(l)*hh(k)y(j)=yy(j)end dodo j=1,my1(j)=y(j);yy(j)=y(j)end dodo l=1,mend dodo j=1,mend do7.4 高階常微分方程（組）的數(shù)值解法10 subroutine rkutta(a,b,