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

1、第一節(jié) 問題的提出含有一個或多個導(dǎo)數(shù)及其函數(shù)的方程式稱為微分方程，在工程中常遇到求解微分方程的問題。很多微分方程的解不能用初等函數(shù)來表示，有時即使能夠用解析式表示其解，但計算量太大而不實用(表達式過于復(fù)雜)。需要用數(shù)值方法來求解，一般只要求得到若干個點上的近似值或者解的簡單的近似表達式(精度要求滿足即可)。 常微分方程（一個自變量）微分方程 偏微分方程（一個以上自變量）一一. 定解問題定解問題實際中求解常微分方程的所謂定解問題有兩類：初值問題和邊值問題定解指已知因變量和/或其導(dǎo)數(shù)在某些點上是已知的(約束條件)。1. 邊解問題約束條件為已知，在自變量的任一非初值上，已知函數(shù)值和/或其導(dǎo)數(shù)值，如

2、y ” = f (x,y,y) 求解y y(a)=，y(b)=常?？梢詫⑦吔鈫栴}轉(zhuǎn)化為初值問題求解。2. 初值問題約束條件為在自變量的初值上已知函數(shù)值，如：y= f (x,y) dy/dx = f (x,y) xx0 y(x0) = y0 y(x0) = y0 求解y(x)，以滿足上述兩式，即在 ax0 x1 xnb上的y(xi)的近似值 yi (i=0,1,2,n)。通常取等距節(jié)點，即h= xi+1-xi ，有 xi = x0+ih(i=0,1,2,n)初值問題的數(shù)值解法特點：按節(jié)點順序依次推進，由已知的y0 , y1 , , yi ，求出yi+1，這可以通過遞推公式得到。 單步法：利用前一

3、個單步的信息 (一個點)，在y=f(x) 上找下 一點yi， 有歐拉法， 龍格庫格法。 初值問題的常見解法 預(yù)測校正法：多步法，利用一個以 上的前點信息求f(x)上 的下一個yi，常用迭 代法，如改進歐拉法， 阿當姆斯法。第二節(jié) 歐拉法y= f (x,y) 求解計算y(xi) (i=1,2,n)的近似值yiy(x0) = y0 一、歐拉折線法一、歐拉折線法1. 解一階方程初值問題的幾何意義 y= f (x,y) (x,y)R，R=axb，-y+也即解y = y(x)所表示的積分曲線 y = f(x,y) dx+c上，每一點(x，y)的切線斜率等于函數(shù)f(x,y)在該點的值，即： y ( xk)

4、 = f (xk，yk)y(x0)=y0，表示積分曲線從P0（x0，y0）點出發(fā)且在P0（x0，y0）的切線斜率為f（x0，y0），故設(shè)想積分曲線y = y (x)在x = x0 附近可用切線近似代替。2. 歐拉折線法在 （ x0，y0 ）點上的切線方程為y=y0+ f (x0 , y0)(x-x0)設(shè)方程與x = x1交點P1（x1，y1）縱坐標y1取為y(x1)的近似值，則有y(x1) y1 = y0+ f (x0,y0)(x1-x0) = y0+ h f (x0,y0)同理 ：在(x1，y1 )上的切線方程 y = y1+ f (x1,y1)(x-x1)與x=x2交點P2(x2,y2)的

5、縱坐標y2取為y(x2)的近似值有y(x2) y2 = y1+ f (x1,y1)(x2-x1) = y1+ h f (x1,y1)上述的h = x2x1= x1x0過Pi(xi,yi)作斜率為f(xi,yi)的切線方程與x = xi+1的交點，有歐拉公式 y(xi+1) yi+1 = yi+ h f (xi，yi) (i=0,1,n-1) xi = x0 +i h (i=1,2,n)例1：求解初值問題y=-2xy 0 x1.8 y(0)=1取h=0.1解： f(x,y) =-2xy， x0=0， y(x0)=1， h=0.1 所以y(xi+1) yi+1 = yi+ h f (xi,yi)

6、= yi +0.1( -2xiyi) = (1-0.2xi)yi 3. 誤差若y(xi) = yi，則y(xi+1) yi+1 0由泰勒展開式在xi上展開有y(xi+1)= y(xi + h)= y(xi) + y (xi) h+(h2/2!) y”()而yi+1 = yi+ h f (xi,yi) = y(xi) + h f(xi, y(xi) = y(xi)+h y(xi)所以y(xi+1) yi+1 = (h2/2) y”() = O(h2)這是局部截斷誤差(因為認為yi = y(xi)實際上可能yi y(xi)，從而有誤差的積累，所以歐拉方法雖然簡單，但精度不夠，很少采用，只是說明這個

7、方法，有助理解其他方法。二、改進歐拉法二、改進歐拉法1. 梯形公式 y= f(x,y) 在xi , xi+1上的積分有即現(xiàn)被積函數(shù)中y (x)未知，用數(shù)值積分(矩形公式)11),(iiiixxxxdxyxfdxy1)(,()()(1iixxiidxxyxfxyxy11)(,()(,()()(1iiiixxxxiiidxxyxfydxxyxfxyxyhxyxfdxxyxfiixxii)(,()(,(1用y(xi) yi ，則y(xi+1) yi+1 = yi + h f (xi , yi)若用梯形公式求積分，則所以實際上是?。?xi，yi）及（xi+1，yi+1）兩點的平均斜率作為在（xi，yi

8、）點的斜率代入歐拉公式中得到上式的?，F(xiàn)有yi+1在公式右邊，如用歐拉公式求出y( xi+1)的一個粗糙近似值yi+1（預(yù)測值），代入梯形公式中得yi+1（校正值）。)(,()(,(2)(,(111iiiixxxyxfxyxfhdxxyxfii)(,()(,(2)(1111iiiiiiixyxfxyxfhyyxyyi+1（校正值） 較yi+1 （預(yù)測值）準確，即有改進的歐拉公式（預(yù)測-校正公式） yi+1 = yi + hf (xi，yi) yi+1= yi + h/2f (xi , y(xi), f (xi+1 , y( xi+1)為便于編程，改為 yp= yi + h f（xi，yi） yi

9、+1= yi + h/2（k1 + k2） yc = yi + h f（xi+1，yp） 或 k1 = f（xi，yi） yi+1= 1/2（ yp + yc ） k2 = f（xi+1，yi + h k1 ）此外為提高精度，可迭代求yi+1，即有 yi+1（0） = yi + hf（xi，yi） yi+1（k+1） = yi + h/2f（xi，yi ）， f（xi+1，yi+1 （k） ）直到 yi+1（k+1） - yi+1（k） 時，取y（ xi+1） yi+1（k+1）例2 用改進歐拉法求例1中的初值問題解： yp= yi + hf（xi，yi）=（10.2 xi ） yi yc =

10、 yi + hf（xi+1，yp）= yi 0.2（ xi +0.1） yp yi+1= 1/2（ yp + yc ） 計算見P144 表7-2-22. 誤差 y(xi+1)= y(xi+h)= y(xi) +h y( xi)+(h2/2!) y”(xi)+ k1= f（xi，yi）= f（xi，y （ xi ）= y（ xi）k2 = f（xi+1，yi + hk1 ）= f（xi +h ， y （ xi ）+ h k1 ）= f（xi，y （ xi ）+h f（xi，y （ xi ）/x + h k1 f（xi，y （ xi ）/y += y（ xi）+h f（xi，y （ xi ）/x+

11、 y（ xi） f（xi，y （ xi ）/y+O(h2)= y（ xi）+h y”（ xi）+ O(h2)所以：yi+1 = yi + h/2 (k1 + k2) = y(xi)+h/2y(xi)+ y(xi)+h y”(xi)+O(h3) = y(xi)+hy(xi)+h2/2 y”(xi)+O(h3)y(xi +1)yi+1 =O(h3 )局部截斷誤差O(hp+1)時，方法具有p階精度。因此歐拉公式有一階精度，而改進歐拉公式有二階精度。3.改進歐拉法N-S圖及程序讀入 x,y讀入數(shù)據(jù) h,xn輸出x,yxxnyl=y+h f(x,y)x = x+h輸出x,y結(jié)束定義函數(shù) f(x,y)=.

12、y=y+0.5*h*f(x,y)+f(x+h,yl)第三節(jié) 龍格庫塔方法一、基本思想一、基本思想根據(jù)微分中值定理有：y(xi +h)= y(xi+1)y (xi)/h (0 1，h= xi+1 - xi)即f (xi +h, y(xi +h) = y(xi+1)yi/h 是平均斜率y(xi+1)= yi +hf (xi +h , y(xi +h) = yi +hk*記 k* = f (xi +h , y(xi +h)而計算k*的方法是歐拉公式中 k*= f (xi , yi )精度低改進歐拉公式中 k*=0.5*f (xi , yi )+ f (xi+1, yi+1) =0.5*k1+ f (

13、xi +1 , yi +hk1) =0.5*k1 + k2 多取了一個點使精度得以提高。如果設(shè)法在xi ， xi+1上多預(yù)報幾個點的斜率值，然后將它們的加權(quán)平均值作為k* = f (xi +h , y(xi +h)的近似值，則有可能構(gòu)造出更高精度的計算公式，這就是龍格庫塔方法的基本思想。二二 、二階龍格庫塔公式、二階龍格庫塔公式在xi ， xi+1內(nèi)有 xi+l = xi +lh（0l 1）取k*=1 k1 + 2 k2 (是xi, xi+1兩個點的斜率值k1, k2加權(quán)平均值）則 yi+1 = yi +hk* = yi +h（1 k1 + 2 k2 ）用歐拉公式，取k1= f（ xi ， y

14、i ）而k2 = f（ xi+l， yi+l）= f（ xi +l， yi + lhk1 ） yi+1= yi + h（ 1 k1 + 2 k2 ） k1 = f（ xi ， yi ） k2 = f（ xi +l， yi + lhk1 ）現(xiàn)計算1 ， 2 及l(fā)。若要使上述公式具有二階精度，即y(xi +1)yi+1 = O(h3)設(shè)y(xi) = yi ，展開y(xi +1)= y (xi)+h y(xi)+ h2/2 y”(xi)+O(h3) k1= f (xi , yi ) = y(xi)k2 = f (xi +l, yi+lhk1) = f(xi , yi)+ lhfx(xi , yi)

15、+ fy(xi , yi)+ f (xi , yi)+O(h2) = y(xi)+ l h y”(xi)+O(h2)yi+1= y(xi)+ h(1 +2) y(xi)+ h22 l y”(xi)+O(h3)具有二階精度，只需 1 + 2 = 1 三個未知量，兩個方程 2 l = 1/2（1）當l=1時，即xi +l =xi +1，解得1 = 2 = 1 /2 yi+1= yi + h/2（ k1 + k2 ） k1 = f（ xi ， yi ） k2 = f（ xi +l， yi + hk1 ）即為改進歐拉公式（2）取l=1/2，則1 =0，2 =1，有變形的歐拉公式 yi+1= yi +

16、h k2 k1 = f（ xi ， yi ） k2 = f（ xi +h/2， yi + h/2 k1 ）三、高階龍格庫塔公式三、高階龍格庫塔公式在xi ， xi+1上除xi ， xi+1外再增加一點xi +m = xi +mh （l m 1） 設(shè)xi +m 處的斜率為k3 ，則：k*= 1k1 +2 k2 +3 k3 yi+1= yi+ hk*= yi+ h (1k1 +2 k2 +3 k3)k1= f（ xi ， yi ） k2 = f（xi +lh， yi+lhk1） k3 = f（xi +m， yi+m）= f（xi +mh， yi+mh（1k1 + 2k2 ）得 yi+1= yi +

17、 h（ 1 k1 + 2 k2 + 3 k3 ） k1 = f（ xi ， yi ） k2 = f（ xi +lh， yi +lhk1 ） k3 = f（xi +mh， yi+mh（1k1 + 2k2 ）利用泰勒展開方法，選取1 , 2 , 3 , l , m及 1 , 2使上述公式具有三階精度O(h4)。得 1 + 2 =1 1 + 2 + 3 =1有7個未知數(shù)，5個方程， 2 l+ 3 m=1/2可得一族解，所得公式 2 l2+ 3 m2=1/3為三階龍格庫塔公式 3 lm 2 =1/61. 庫塔公式取l=1/2, m=1, 則1 =1/6, 2 =2/3, 3 =1/6, 1 =-1,2

18、 =2得庫塔公式 yi+1= yi + h/6（ k1 + 4k2 + k3 ） k1 = f（ xi ， yi ） k2 = f（ xi +h/2， yi +h/2k1 ） k3 = f（xi +h， yi - hk1+ 2hk2 ）2. 四階龍格庫塔公式在xi ， xi+1上用四個點的ki (i=1,2,3,4)做k*的加權(quán)平均值，可得一族經(jīng)典四階龍格庫塔公式。 yi+1= yi + h/6（ k1 + 2k2 +2 k3 + k4 ） k1 = f（ xi ， yi ） k2 = f（ xi +h/2， yi +h/2k1 ） k3 = f（ xi +h/2 ， yi+ h/2k2 ）

19、k4 = f（ xi +h， yi+ hk3）例3：用四階龍格庫塔方法求解y=2xy 0 x1.8 y(0)=1取h=0.2局部截斷誤差O(h5)，精度提高。解： yi+1= yi + 0.2/6（ k1 + 2k2 +2 k3 + k4 ） k1 = -2 xi yi k2 = f（xi +0.1， yi +0.1k1 ）=-2（xi +0.1）（yi +0.1k1 ） k3 = f（xi +0.1， yi+ 0.1k2 ）=-2（xi +0.1）（yi +0.1k2 ） k4 = f（ xi +0.2，yi+ 0.2k3）=-2（xi +0.2）（yi +0.2k3 ）計算結(jié)果見P150表

20、7-3-1，在表7-3-2中對比了幾種方法的誤差。表7-3-3說明再提高階數(shù)沒有多大意義。四階是兼顧了精度和計算量的理想公式。四、四、 四階龍格四階龍格-庫塔公式的幾何意義。庫塔公式的幾何意義。xixi+hABECDR SPi（ xi， yi ）hk1斜率為f（ xi， yi ）的PA，AM=hki，PA中點R(xi+h/2,yi+h/2ki)的斜率作PB, BM=hk2 ，有PB中點S (xi+h/2,yi+h/2k2)的斜率作PC，CM=hk3，C點(xi+h,yi+hk3)的斜率作PD，DM=hk4，取EM=hk*=h/6*(k1+2k2+2k3+k4 )y(xi+h)= yi+h =

21、yi + 1/6*(k1+2k2+2k3 +k4 )Mhk2hk3hk4五、五、(步長的選擇步長的選擇)自動定步長四階龍格庫塔自動定步長四階龍格庫塔設(shè)從xi出發(fā) ，先以h為步長，用四階龍格庫塔公式求得y(xi+1)的近似值yi+1(h) ，則y(xi+1)yi+1(h)ch5再以h/2為步長，兩次計算后(先yi+h/2，再算yi+1)可得yi+1 y(xi+1)yi+1(h/2)c(h/2)5 = c/16* h5有 y(xi+1)yi+1(h/2) yi+1(h/2) yi+1(h) /15即利用事后誤差估計法 | yi+1(h/2) yi+1(h) |，判斷是否停止計算。六、四階龍格庫塔法

22、N-S圖讀入x0, y0讀入數(shù)據(jù) h讀入 xnx0 xnk1, k2, k3, k4y0= y0+ h/6(k1+2k2+2k3+k4)x0 = x0+h結(jié)束定義函數(shù) f(x, y)輸出k1, k2, k3, k4 , x0 , y04 線形多步法線形多步法l單步法中，求解yi+1時，實際已求出了y0，y1yi 。利用多個y值，期望得到較高精度的yi+1ly（xi+1）= y （xi） + xi+1 f（ x ， y （x） ） dxl xil用更精確的求積方法，可用更高次的插值多次來代f（x，y），選取的插值節(jié)點不同，則解法不同。l一 阿當姆斯內(nèi)插公式及誤差l除取xi ， xi+1兩個節(jié)點外

23、，其他節(jié)點取在xi ， xi+1時，f 的值未知，故選取xi ， xi+1以外的節(jié)點為xi-1 ， xi-2等。l現(xiàn)取xi-1 ， xi-2 ， xi ， xi+1四個插值節(jié)點，由于計算yi+1時， yi-2 ， yi-1 ， yi ，已求解出來，插值多項式lL3（x）=（x- xi ）（x- xi-1）（x- xi-2）/（xi+1- xi ）（ xi+1 - xi-1）（ xi+1 - xi-2） f（ xi+1 ， y （ xi+1 ） ）+ （x- xi+1 ）（x- xi-1）（x- xi-2）/（xi- xi+1 i ）（ xi - xi-1）（ xi - xi-2） f（ xi

24、， y （ xi ） ）+ （x- xi+1 ）（x- xi）（x- xi-2）/（xi-1- xi+1 ）（ xi-1 - xi）（ xi-1 - xi-2） f （ xi-1 ， y （ xi-1 ） ）+ （x- xi+1 ）（x- xi）（x- xi-1）/（xi-2- xi ）（ xi-2 - xi-1）（ xi-2 - xi+1） f （ xi-2 ， y （ xi-2 ） ）lf（ x ， y （x） ） L3（x）l所以yi+1 y （ xi） + xi+1 L3（x）dxl xil等距節(jié)點時， xi+1 - xi = xi - xi-1 = xi+1 - xi-2 = h，設(shè)

25、yi= f（ xi ， y （xi） ）l令x= xi+th 0t 1lyi+1 = yi + h/6 f（ xi+1 ， y （xi+1） ） 1t（t-1）（t+2）dt l 0l-h/2 f（ xi ， y （xi） ） 1（t-1）（t+1）（t+2）dtl 0l+h/2 f（ xi-1 ， y （xi-1） ） 1（t-1）t（t+2）dtl 0l-h/6 f（ xi-2， y （xi-2） ） 1（t-1）t（t+1）dtl 0l=y（ xi ）+h/249 yi+1+19 yi-5 yi-1+ yi-2 l阿當姆斯內(nèi)插公式是四階公式ly（ xi+1）- yi+1 =（-19/72

26、0）h5yi（5）+ =0（h5）l內(nèi)插公式中由于選取了xi+1 點為節(jié)點，從而使公式成為隱式，若取xi-3 ， xi-2 ， xi-1 ， xi 作為插值節(jié)點，即：l i ilL3（x）= (x- xj )/(xk-xj)f(xj -y(xj )l k=i-3 j=I-3jkl所以y(xj +1)= yi+1 = xi+1 L3（x）dxl xil= yi +h/2455 yi+59 yi-1+37 yi-2-9 yi-3 l外推公式ly(xj +1)- yi+1 = （25/720）h5yi（5）+ =0（h5）l計算f（x，y）的次數(shù)，算yi+1 時 yi-1 ，yi-2， yi-3 已

27、在計算yi時得到，故只要算yi就可，故無需算四次，故計算變化四階龍格-庫塔公式，而0（h5）相同，但需用其它方法先行計算y1 ， y2 ， y3 后才能用該法。依次計算y4 ，y5 ，所以不具有自啟動性，此外h要改變也很麻煩。l若將外推公式作預(yù)測式，內(nèi)插公式作校正，可得：l 阿當姆斯預(yù)測校正式lyj +1（0）= yi +h/2455 yi+59 yi-1+37 yi-2-9 yi-3 lyj +1= yi +h/249f(xj +1 ， yj +1（0)+19 yi+5 yi-1- yi-2 l例4：ly=-2xy 0 x1.2l ly(0)=1l取h=0.2l解：四階龍格-庫塔求出， y1， y2 ， y3 ，由 y0 ， y1 ， y2 ， y3阿當姆斯預(yù)測校正公式。計算見p159 表7-6，7-7l一 一階方程組lu=（x，u，v）, u(x 0)= u0l lv=(x，u，v) , v(x 0)= v0 初值解l l記y=（u，v）T，y0=（ u0 ， v0 ）T F=（ ， ）Tl l y= F（ x ， y ）l l y(x 0)= y0l用龍格庫塔解l yi+1= yi + h/6（ k1 + 2k2 +2 k3 + k4 ）l k1 = f（ xi ， yi ）l k2 = f（xi +h/2， yi + h/2k1 ）l k3 = f