常微分方程数值解法课件

常微分方程及方程組的數值解法

常微分方程的定解問題一般分為初值問題和邊值問題，而邊值問題常常可以化為初值問題來求解，所以本章只討論初值問題。設在X—Y平面有：當

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

(7-1)§7.1常微分方程的數值法

本節研究常用的數值方法。先說明幾個符號的含義：一、Taylor級數展開法把yi+1在yi處展開成Taylor級數，則有因為所以

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

§7.1常微分方程的數值法

這就是常微分方程數值解法中的向前歐拉(Euler)公式。二、數值微分法這種方法就是利用差商近似代替導數，即(7-2)§7.1常微分方程的數值法略去二階導數項，用代替，得利用它們求問題(7-1)式數值的方法，分別稱為歐拉(Euler)法，後退歐拉法，中點法。三、數值積分法將(7-1)式在區間內積分得(7-3)§7.1常微分方程的數值法對上式右端的積分使用左端矩形積分公式，則有：於是(7-1)式的數值解法但如果上式右端的積分使用的是梯形積分公式，則有從而(7-1)式的數值解法也改為這就是改進的歐拉公式。(7-4)(7-5)§7.1常微分方程的數值法四、幾個基本概念（1）單步法：總是由前一點yi推算後一點yi+1，只需要一個點的值。除此之外，還有多步法，如Adams法。（2）顯示法：yi+1只出現等式一邊。若同時出現在右邊，則為隱式法。（3）截項誤差：當h趨於零時，它是h2的同階無窮小，（7-2-1）誤差階數為2。如果y(x)是一次式，則O(h2)

0,此時歐拉法是精確的（直線逼近直線），同時也說明歐拉法是一階代數精度。§7.2一階常微分方程的數值法讓我們研究(7-1)式形式的一階常微分方程的數值求解方法。7.2.1歐拉折線法一、演算法分析二、幾何意義參見圖7-1中的折線f01f11f21f31，歐拉折線法命名由此而來。容易看出，歐拉折線法的數值解離理論曲線越來越遠。折線f01f11f21f31為後面要討論的歐拉二次法。(7-6)§7.2一階常微分方程的數值法圖7-1歐拉折線法及改進的歐拉法的幾何意義§7.2一階常微分方程的數值法三、通用程式框圖與程式設計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一階常微分方程的數值法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一階常微分方程的數值法7.2.2改進的歐拉折線法（梯形法）一、演算法分析在歐拉法中，泰勒級數是取前兩項。下麵看看取三項時是什麼樣子。為了求yi+1，在yi處進行泰勒級數展開。依據前面所講的，可以寫成：因為所以

§7.2一階常微分方程的數值法這就是歐拉二次逼近法的遞推式，其局部誤差為O(h3)。可見，比歐拉法精度高（隱式）。在上式中，等號的右邊也含有yi+1,要出yi+1一般要用迭代法。為了簡便，用歐拉法先求出y*i+1，作為yi+1的預報值。於是，歐拉二次逼近法的演算法可表示為：二、通用程式框圖與程式設計(7-7)§7.2一階常微分方程的數值法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一階常微分方程的數值法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一階常微分方程的數值法7.2.2Runge-Kutta方法一、演算法分析一階常微分初值問題的求解，在社區間上歸結為如何選擇去逼近該區間曲線的直線。如果直線已經確定，已知yi，就很容易求出yi+1。依據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）座標，就能確定直線：§7.2一階常微分方程的數值法上式中的f為社區間曲線的平均斜率。剩下的問題就是如何確定該平均斜率？斜率的不同選擇方法或演算法，決定不同的計算格式。如：用起始點的斜率，則Ki=y’(xi)是曲線在xi點的斜率，用此斜率代替平均斜率。上式就是歐拉演算法。如：用區間上兩端點的斜率的平均值，則其中，Ki=f(xi,yi)，Ki+1=f(xi+1,yi+1)，用兩點斜率的算術平均值代替平均斜率。此處的yi+1用歐拉法作預報值，則

(7-8)(7-9)§7.2一階常微分方程的數值法

這就是歐拉二次逼近法的演算法式。（7-8）式比（7-9）式的計算精度要高，這說明增加f(x,y)或斜率的計算次數可以提高截斷的階數，從而提高了精度。這就是Runge-Kutta法的基本思想。將不同點的斜率進行線性組合，其一般的形式為：只要適當選取

1、

2、a、b的值，就能保證局部截斷誤差為O(h3)。其中一個特例是

1=1/2=

2，a=b=1,此時就變成了（7-9）式，亦稱二階Runge-Kutta法。四階Runge-Kutta法就是在社區間上計算四個斜率，爾後進行線性組合，一般形式為：§7.2一階常微分方程的數值法

適當選擇待定係數，可得到截斷誤差O(h5)。四階龍格—庫塔法的常用演算法為：§7.2一階常微分方程的數值法二、通用程式框圖與程式設計(7-10)§7.2一階常微分方程的數值法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一階常微分方程的數值法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一階常微分方程組

的數值解法

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

(7-11)§7.3一階常微分方程組

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

的數值解法

在求時，必須把所有的（）都先求出來，因為括弧內的通式為；同理，求時，必須把所有的都先求出來，因為括弧內的通式為；求時，必須把所有的都先求出來，因為括弧內的通式為。二、通用程式框圖和程式設計

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一階常微分方程組

的數值解法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高階常微分方程（組）

的數值解法

一、演算法分析關於高階常微分方程（組）初值問題的數值解法，一般可以引入新的變元，把它們化為一階常微分方程組的初值問題用Runge-Kutta法來求解。一般地，對於初值問題：

則（7.11）式可被化為等價的初值問題：(7-11)§7.4高階常微分方程（組）

的數值解法則可用前面介紹的任何一種方法求解（7.12）式。(7-12)§7.4高階常微分方程（組）

的數值解法例7-1某振動過程的動態模型如下，求其在區間[0，0.2]內的數值解，取h=0.2。解：這是一個二階常微分方程組的初值問題，令,可將原方程組化為如下的等價形式：

§7.4高階常微分方程（組）

的數值解法解：這是二、通用程式框圖與程式設計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高階常微分方程（組）

的數值解法10subroutinerkutta(a,b,n,m,y,yy)20