第7章常微分方程及方程組的數(shù)值解法

常微分方程的定解問題一般分為初值問題和邊值問題，而邊值問題常?？梢曰癁槌踔祮栴}來求解，所以本章只討論初值問題。設(shè)在X—Y平面有：當(dāng)

求。在工程問題中，只有極少有解析解，大多數(shù)只能求數(shù)值解。即在區(qū)域：求y(xi)的近似解：。

(7-1)§7.1常微分方程的數(shù)值法本節(jié)研究常用的數(shù)值方法。先說明幾個符號的含義：一、Taylor級數(shù)展開法把yi+1在yi處展開成Taylor級數(shù)，則有因為所以

O(h2)為h2的同階無窮小，略去高階項得如下的近似關(guān)系：

§7.1常微分方程的數(shù)值法這就是常微分方程數(shù)值解法中的向前歐拉(Euler)公式。二、數(shù)值微分法這種方法就是利用差商近似代替導(dǎo)數(shù)，即(7-2)§7.1常微分方程的數(shù)值法略去二階導(dǎo)數(shù)項，用代替，得利用它們求問題(7-1)式數(shù)值的方法，分別稱為歐拉(Euler)法，后退歐拉法，中點法。三、數(shù)值積分法將(7-1)式在區(qū)間內(nèi)積分得(7-3)§7.1常微分方程的數(shù)值法對上式右端的積分使用左端矩形積分公式，則有：于是(7-1)式的數(shù)值解法但如果上式右端的積分使用的是梯形積分公式，則有從而(7-1)式的數(shù)值解法也改為這就是改進(jìn)的歐拉公式。(7-4)(7-5)§7.1常微分方程的數(shù)值法四、幾個基本概念（1）單步法：總是由前一點yi推算后一點yi+1，只需要一個點的值。除此之外，還有多步法，如Adams法。（2）顯示法：yi+1只出現(xiàn)等式一邊。若同時出現(xiàn)在右邊，則為隱式法。（3）截項誤差：當(dāng)h趨于零時，它是h2的同階無窮小，（7-2-1）誤差階數(shù)為2。如果y(x)是一次式，則O(h2)0,此時歐拉法是精確的（直線逼近直線），同時也說明歐拉法是一階代數(shù)精度?！?.2一階常微分方程的數(shù)值法讓我們研究(7-1)式形式的一階常微分方程的數(shù)值求解方法。7.2.1歐拉折線法一、算法分析二、幾何意義參見圖7-1中的折線f01f11f21f31，歐拉折線法命名由此而來。容易看出，歐拉折線法的數(shù)值解離理論曲線越來越遠(yuǎn)。折線f01f11f21f31為后面要討論的歐拉二次法。(7-6)§7.2一階常微分方程的數(shù)值法圖7-1歐拉折線法及改進(jìn)的歐拉法的幾何意義§7.2一階常微分方程的數(shù)值法三、通用程序框圖與程序設(shè)計subroutineeulercompf(x0,y0)x0=ax=x0+hdoi=1,n+1y=y0+h*f(x0,y0)enddoendsubroutineeulerx=x0;y0=yH=(b-a)/(n+1)§7.2一階常微分方程的數(shù)值法10subroutineeuler(a,b,n,y0,y)20x0=a30h=(b-a)/(n+1)40doi=1,n+150x=x0+h60y=y0+h*f(x0,y0)70x0=x;y0=y80enddo90endsubroutineeuler§7.2一一階常微分方方程的數(shù)值法法7.2.2改改進(jìn)的歐拉拉折線法（梯梯形法）一、算法分析析在歐拉法中，，泰勒級數(shù)是是取前兩項。。下面看看取取三項時是什什么樣子。為了求yi+1，在yi處進(jìn)行泰勒級級數(shù)展開。依依據(jù)前面所講講的，可以寫寫成：因為所以§7.2一一階常微分方方程的數(shù)值法法這就是歐拉二二次逼近法的的遞推式，其其局部誤差為為O(h3)?？梢?，比歐歐拉法精度高高（隱式）。。在上式中，，等號的右邊邊也含有yi+1,要出yi+1一般要用迭代代法。為了簡簡便，用歐拉拉法先求出y*i+1，作為yi+1的預(yù)報值。于于是，歐拉二二次逼近法的的算法可表示示為：二、通用程序序框圖與程序序設(shè)計(7-7)§7.2一一階常微分方方程的數(shù)值法法subroutinemeulercompf(x0,y0)x0=a;y=y0x=x0+hdoi=1,n+1y1=y0+h*f(x0,y0)enddoendsubroutinemeulerx=x0;y0=yH=(b-a)/(n+1)y=y0+h*(f(x0,y0)+f(x,y1))/2compf(x,y1)§7.2一一階常微分方方程的數(shù)值法法10subroutinemeuler(a,b,n,y0,y)20x0=a;y=y030h=(b-a)/(n+1)40doi=1,n+150x=x0+h60y1=y0+h*f(x0,y0)70y=y0+h*(f(x0,y0)+f(x,y1))/270x0=x;y0=y80enddo90endsubroutinemeuler§7.2一一階常微分方方程的數(shù)值法法7.2.2Runge-Kutta方法一、算法分析析一階常微分初初值問題的求求解，在小區(qū)區(qū)間上歸結(jié)為為如何選擇去去逼近該區(qū)間間曲線的直線線。如果直線線已經(jīng)確定，，已知yi，就很容易求求出yi+1。依據(jù)Lagrange微分中值定定理，有：其中，h=xi+1－xi;0<<1。同時有y’(xi+h)=f(x,y)=f(xi+h,y(xi+h))，這就是曲線線y=y(x)在xi+h處切線的斜率率。如果已確定，則斜斜率也就確定定。有了斜率率和已知點（xi,yi）坐標(biāo)，就能確確定直線：§7.2一一階常微分方方程的數(shù)值法法上式中的f為為小區(qū)間曲線線的平均斜率率。剩下的問問題就是如何何確定該平均均斜率？斜率率的不同選擇擇方法或算法法，決定不同同的計算格式式。如：用起起始點的斜率率，則Ki=y’(xi)是曲線在xi點的斜率，用用此斜率代替替平均斜率。。上式就是歐歐拉算法。如如：用區(qū)間上上兩端點的斜斜率的平均值值，則其中，Ki=f(xi,yi)，Ki+1=f(xi+1,yi+1)，用兩點斜率率的算術(shù)平均均值代替平均均斜率。此處處的yi+1用歐拉法作預(yù)預(yù)報值，則(7-8)(7-9)§7.2一一階常微分方方程的數(shù)值法法這就是歐拉二二次逼近法的的算法式。（（7-8）式比比（7-9））式的計算精度度要高，這說說明增加f(x,y)或或斜率的計算算次數(shù)可以提提高截斷的階階數(shù)，從而提提高了精度。。這就是Runge-Kutta法的基本思想想。將不同點點的斜率進(jìn)行行線性組合，，其一般的形形式為：只要適當(dāng)選取取1、2、a、b的值，就能保保證局部截斷斷誤差為O(h3)。其中一個特特例是1=1/2=2，a=b=1,此時就變成了了（7-9）式，亦稱二階階Runge-Kutta法。四階Runge-Kutta法就是在小區(qū)區(qū)間上計算四四個斜率，爾爾后進(jìn)行線性性組合，一般般形式為：§7.2一一階常微分方方程的數(shù)值法法適當(dāng)選擇待定定系數(shù)，可得得到截斷誤差差O(h5)。四階龍格格—庫塔法的的常用算法為為：§7.2一一階常微分方方程的數(shù)值法法二、通用程序序框圖與程序序設(shè)計(7-10)§7.2一一階常微分方方程的數(shù)值法法subroutinerkuttax0=a;y=y0x1=x0+(i-1)*hdoi=1,nx=x1+hh(j)enddoendsubroutineykuttayy=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/2doj=1,4ff=f(x,y)compf(x,y)y=y1+ff*hh(j)y=yyenddo§7.2一一階常微分方方程的數(shù)值法法10subroutinerkutta(a,b,n,y0,yy)15dimensionhh(5)20x0=a;y=y030h=(b-a)/n40hh(1)=h/2;hh(2)=hh(1)50hh(3)=h;hh(4)=h;hh(5)=h/260doi=1,n70x1=x0+(i-1)*h80x=x1;y1=y;yy=y90doj=1,4100ff=f(x,y)110x=x1+hh(j)120yy=yy+ff*hh(j+1)/3130y=y1+ff*hh(j)140enddo150y=yy160enddo170endsubroutinerkutta§7.3一階常微分方方程組

的數(shù)數(shù)值解法我們所要解決決的實際問題題，往往是由由多個常微分分方程組成的的聯(lián)立方程組組，所以本節(jié)節(jié)我們討論一一階常微分方方程組的數(shù)值值解法。一階常微分方方程組的一般般形式為：一、算法概述述我們可以把求求單個常微分分方程的四階階Runge-Kutta法的遞推計算公公式（7.10）式推廣到求解常微分分方程組的情情形，即：(7-11)§7.3一階常微分方方程組

的數(shù)數(shù)值解法(7-12)§7.3一階常微分方方程組

的數(shù)數(shù)值解法在求時時，必須把把所有的（（））都都先求出來，，因為括號內(nèi)內(nèi)的通式為；；同同理，求時時，，必須把所有有的都都先求出來來，因為括號號內(nèi)的通式為為；；求時時，必須須把所有的都都先先求出來，因因為括號內(nèi)的的通式為。。二、通用程序序框圖和程序序設(shè)計subroutinehrkuttax0=ax1=x0+(i-1)*hdoi=1,nx=x1+hh(k)enddoendsubroutinehrkuttayy(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/2dok=1,4callfun(x,y,f)y(l)=y1(l)+f(l)*hh(k)y(j)=yy(j)enddodoj=1,my1(j)=y(j);yy(j)=y(j)enddodol=1,menddodoj=1,menddoready(i),i=1,m§7.3一階常微分方方程組

的數(shù)數(shù)值解法10subroutineorkutta(a,b,n,m,yy)20dimensionhh(5),y(m),y1(m),f(m)30&yy(mm)40x0=a50doj=1,m60read(5,’f12.5’)y(j)70enddo80h=(b-a)/n90hh(1)=h/2;hh(2)=hh(1)100hh(3)=h;hh(4)=h;hh(5)=h/2110doi=1,n120x1=x0+(i-1)*h130x=x1140doj=1,m150y1(j)=y(j);yy(j)=y(j)160enddo170dok=1,4180callfun(x,y,f)190x=x1+hh(k)200dol=1,m210yy(l)=yy(l)+f(l)*hh(k+1)/3220y(l)=y1(l)+f(l)*hh(k)230enddo240enddo250doj=1,m260y(j)=yy(j)270enddo280enddo290endsubroutineorkutta§7.4高階常微分方方程（組）的的數(shù)值解法法一、算法分析析關(guān)于高階常微微分方程（組組）初值問題題的數(shù)值解法法，一般可以以引入新的變變元，把它們們化為一階常常微分方程組組的初值問題題用Runge-Kutta法來求解。一般般地，對于初初值問題：則（7.11）式可被化化為等價的初初值問題：(7-11)§7.4高階常微分方方程（組）的的數(shù)值解法法則可用前面介介紹的任何一一種方法求解解（7.12）式。(7-12)§7.4高階常微分方方程（組）的的數(shù)值解法法例7-1某振動過程的的動態(tài)模型如如下，求其在在區(qū)間[0，，0.2]內(nèi)內(nèi)的數(shù)值解，，取h=0.2。解：這是一個個二階常微分分方程組的初初值問題，令令,可將原方程組組化為如下的的等價形式：：§7.4高階常微分方方程（組）的的數(shù)值解法法解：這是二、通用程序序框圖與程序序設(shè)計subroutinehrkuttax0=a;y=y0x1=x0+(i-1)*hdoi=1,nx=x1+hh(k)enddoendsubroutinehrkuttayy(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/2dok=1,4callfun(x,y,f)y(l)=y1(l)+f(l)*hh(k)y(j)=yy(j)enddodoj=1,my1(j)=y(j);yy(j)=y(j)enddodol=1,menddodoj=1,menddo§7.4高階常微分方方程（組）的的數(shù)值解