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

2、初值上，已知函數(shù)值和/或其導(dǎo)數(shù)值，如 y ” = f (x,y,y) 求解y y(a)=，y(b)= 常?？梢詫⑦吔鈫?wèn)題轉(zhuǎn)化為初值問(wèn)題求解。,2. 初值問(wèn)題 約束條件為在自變量的初值上已知函數(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é)點(diǎn)，即h= xi+1-xi ，有 xi = x0+ih(i=0,1,2,n) 初值問(wèn)題的數(shù)值解法特點(diǎn)：按節(jié)點(diǎn)順序依次推進(jìn)，由已知的y0 , y1 , , yi ，求出

3、yi+1，這可以通過(guò)遞推公式得到。,單步法：利用前一個(gè)單步的信息 (一個(gè)點(diǎn))，在y=f(x) 上找下 一點(diǎn)yi， 有歐拉法， 龍格庫(kù)格法。 初值問(wèn)題的常見解法 預(yù)測(cè)校正法：多步法，利用一個(gè)以 上的前點(diǎn)信息求f(x)上 的下一個(gè)yi，常用迭 代法，如改進(jìn)歐拉法， 阿當(dāng)姆斯法。,第二節(jié) 歐拉法 y= f (x,y) 求解計(jì)算y(xi) (i=1,2,n)的近似值yi y(x0) = y0 一、歐拉折線法 1. 解一階方程初值問(wèn)題的幾何意義 y= f (x,y) (x,y)R，R=axb，-y+ 也即解y = y(x)所表示的積分曲線 y = f(x,y) dx+c 上，每一點(diǎn)(x，y)的切線斜率等

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

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

6、， h=0.1 所以y(xi+1) yi+1 = yi+ h f (xi,yi) = yi +0.1( -2xiyi) = (1-0.2xi)yi 計(jì)算見P142 表7-2-1。,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) 這是局部截?cái)嗾`差(因?yàn)?/p>

7、認(rèn)為yi = y(xi) 實(shí)際上可能yi y(xi)，從而有誤差的積累，所以歐拉方法雖然簡(jiǎn)單，但精度不夠，很少采用，只是說(shuō)明這個(gè)方法，有助理解其他方法。,二、改進(jìn)歐拉法 1. 梯形公式 y= f(x,y) 在xi , xi+1上的積分有 即 現(xiàn)被積函數(shù)中y (x)未知，用數(shù)值積分(矩形公式),用y(xi) yi ，則 y(xi+1) yi+1 = yi + h f (xi , yi) 若用梯形公式求積分，則 所以 實(shí)際上是?。?xi，yi）及（xi+1，yi+1）兩點(diǎn)的平均斜率作為在（xi，yi）點(diǎn)的斜率代入歐拉公式中得到上式的。 現(xiàn)有yi+1在公式右邊，如用歐拉公式求出y( xi+1)的一個(gè)

8、粗糙近似值yi+1（預(yù)測(cè)值），代入梯形公式中得yi+1（校正值）。,yi+1（校正值） 較yi+1 （預(yù)測(cè)值）準(zhǔn)確，即有改進(jìn)的歐拉公式（預(yù)測(cè)-校正公式） 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+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+

9、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） 時(shí)，取y（ xi+1） yi+1（k+1）,例2 用改進(jìn)歐拉法求例1中的初值問(wèn)題 解： yp= yi + hf（xi，yi）=（10.2 xi ） yi yc = yi + hf（xi+1，yp）= yi 0.2（ xi +0.1） yp yi+1= 1/2（ yp + yc ） 計(jì)算見P144 表7-2-2 2. 誤差 y(xi+1)= y(xi+h)= y(xi) +h y( xi)+(h

10、2/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+ 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

11、”(xi)+O(h3) = y(xi)+hy(xi)+h/2 y”(xi)+O(h3) y(xi +1)yi+1 =O(h3 ) 局部截?cái)嗾`差O(hp+1)時(shí)，方法具有p階精度。 因此歐拉公式有一階精度，而改進(jìn)歐拉公式有二階精度。,3.改進(jìn)歐拉法N-S圖及程序,讀入 x,y,讀入數(shù)據(jù) h,xn,輸出x,y,xxn,yl=y+h f(x,y),x = x+h,輸出x,y,結(jié)束,定義函數(shù) f(x,y)=.,y=y+0.5*h*f(x,y)+f(x+h,yl),第三節(jié) 龍格庫(kù)塔方法 一、基本思想 根據(jù)微分中值定理有： y(xi +h)= y(xi+1)y(xi)/h (0 1，h= xi+1 - x

12、i) 即 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) 而計(jì)算k*的方法是 歐拉公式中 k*= f (xi , yi )精度低 改進(jìn)歐拉公式中 k*=0.5*f (xi , yi )+ f (xi+1, yi+1) =0.5*k1+ f (xi +1 , yi +hk1) =0.5*k1 + k2 多取了一個(gè)點(diǎn)使精度得以提高。,如果設(shè)法在xi ， xi+1上多預(yù)報(bào)幾個(gè)點(diǎn)的斜率值，然后將它們的加權(quán)平均值作為k* = f (x

13、i +h , y(xi +h)的近似值，則有可能構(gòu)造出更高精度的計(jì)算公式，這就是龍格庫(kù)塔方法的基本思想。 二 、二階龍格庫(kù)塔公式 在xi ， xi+1內(nèi)有 xi+l = xi +lh（0l 1） 取k*=1 k1 + 2 k2 (是xi, xi+1兩個(gè)點(diǎn)的斜率值k1, k2加權(quán)平均值） 則 yi+1 = yi +hk* = yi +h（1 k1 + 2 k2 ） 用歐拉公式，取k1= f（ xi ， yi ）而 k2 = f（ xi+l， yi+l）= f（ xi +l， yi + lhk1 ） yi+1= yi + hf（ 1 k1 + 2 k2 ） k1 = f（ xi ， yi ） k2

14、 = f（ xi +l， yi + lhk1 ）,現(xiàn)計(jì)算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)+ 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

15、 l y”(xi)+O(h3),具有二階精度，只需 1 + 2 = 1 三個(gè)未知量，兩個(gè)方程 2 l = 1/2 （1）當(dāng)l=1時(shí)，即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 ）即為改進(jìn)歐拉公式 （2）取l=1/2，則1 =0，2 =1，有變形的歐拉公式 yi+1= yi + h k2 k1 = f（ xi ， yi ） k2 = f（ xi +h/2， yi + h/2 k1 ）,三、高階龍格庫(kù)塔公式 在xi ， xi+1上除xi ， xi+1

16、外再增加一點(diǎn) 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 + h（ 1 k1 + 2 k2 + 3 k3 ） k1 = f（ xi ， yi ） k2 = f（ xi +lh， yi +lhk1 ） k3 = f（xi +mh， yi+mh

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

18、 ） k3 = f（xi +h， yi - hk1+ 2hk2 ）,2. 四階龍格庫(kù)塔公式 在xi ， xi+1上用四個(gè)點(diǎn)的ki (i=1,2,3,4)做k*的加權(quán)平均值，可得一族經(jīng)典四階龍格庫(kù)塔公式。 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 ） k4 = f（ xi +h， yi+ hk3） 例3：用四階龍格庫(kù)塔方法求解 y=2xy 0 x1.8 y(0)=1取h=0.2,局部截?cái)嗾`差O(h5)，精度提高。,

19、解： 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 ） 計(jì)算結(jié)果見P150表7-3-1，在表7-3-2中對(duì)比了幾種方法的誤差。表7-3-3說(shuō)明再提高階數(shù)沒有多大意義。四階是兼顧了精度和計(jì)算量的理想公式。,四、 四階龍格-庫(kù)塔公式的幾

20、何意義。,xi,xi+h,A,B,E,C,D,R,S,Pi（ xi， yi ）,hk1,斜率為f（ xi， yi ）的PA，AM=hki，PA中點(diǎn)R(xi+h/2,yi+h/2ki)的斜率作PB, BM=hk2 ，有PB中點(diǎn)S (xi+h/2,yi+h/2k2)的斜率作PC，CM=hk3，C點(diǎn)(xi+h,yi+hk3)的斜率作PD，DM=hk4，取EM=hk* =h/6*(k1+2k2+2k3+k4 ) y(xi+h)= yi+h = yi + 1/6*(k1+2k2+2k3 +k4 ),M,hk2,hk3,hk4,五、自動(dòng)定步長(zhǎng)四階龍格庫(kù)塔 設(shè)從xi出發(fā) ，先以h為步長(zhǎng)，用四階龍格庫(kù)塔公式求

21、得y(xi+1)的近似值yi+1(h) ，則 y(xi+1)yi+1(h)ch5 再以h/2為步長(zhǎng)，兩次計(jì)算后(先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 即利用事后誤差估計(jì)法 | yi+1(h/2) yi+1(h) |，判斷是否停止計(jì)算。,六、四階龍格庫(kù)塔法N-S圖,讀入x0, y0,讀入數(shù)據(jù) h,讀入 xn,x0xn,k1, k2, k3, k4,y0= y0+ h/6(k1+2k2+2k3+k4),x0 = x0+h,結(jié)束,定義函數(shù) f(

22、x, y),輸出k1, k2, k3, k4 , x0 , y0,4 線形多步法,單步法中，求解yi+1時(shí)，實(shí)際已求出了y0，y1yi 。利用多個(gè)y值，期望得到較高精度的yi+1 y（xi+1）= y （xi） + xi+1 f（ x ， y （x） ） dx xi 用更精確的求積方法，可用更高次的插值多次來(lái)代f（x，y），選取的插值節(jié)點(diǎn)不同，則解法不同。 一 阿當(dāng)姆斯內(nèi)插公式及誤差 除取xi ， xi+1兩個(gè)節(jié)點(diǎn)外，其他節(jié)點(diǎn)取在xi ， xi+1時(shí)，f 的值未知，故選取xi ， xi+1以外的節(jié)點(diǎn)為xi-1 ， xi-2等。 現(xiàn)取xi-1 ， xi-2 ， xi ， xi+1四個(gè)插值節(jié)點(diǎn)，由

23、于計(jì)算yi+1時(shí)， yi-2 ， yi-1 ， yi ，已求解出來(lái)，插值多項(xiàng)式,L3（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 ， y （ xi ） ）+ （x- xi+1 ）（x- xi）（x- xi-2）/（xi-1- xi+1 ）（ xi-1 - xi）（ xi-1 - xi-2） f （ xi-1

24、， 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 ） ） f（ x ， y （x） ） L3（x） 所以yi+1 y （ xi） + xi+1 L3（x）dx xi 等距節(jié)點(diǎn)時(shí)， xi+1 - xi = xi - xi-1 = xi+1 - xi-2 = h，設(shè)yi= f（ xi ， y （xi） ） 令x= xi+th 0t 1,yi+1 = yi + h/6 f（ xi+1 ， y （xi+1） ） 1t（t-1）（t+2）dt 0 -

25、h/2 f（ xi ， y （xi） ） 1（t-1）（t+1）（t+2）dt 0 +h/2 f（ xi-1 ， y （xi-1） ） 1（t-1）t（t+2）dt 0 -h/6 f（ xi-2， y （xi-2） ） 1（t-1）t（t+1）dt 0 =y（ xi ）+h/249 yi+1+19 yi-5 yi-1+ yi-2 阿當(dāng)姆斯內(nèi)插公式是四階公式 y（ xi-1）- yi+1 =（-19/720）h5yi（5）+ =0（h5）,二 阿當(dāng)姆斯外推公式及誤差,內(nèi)插公式中由于選取了xi+1 點(diǎn)為節(jié)點(diǎn)，從而使公式成為隱式，若取xi-3 ， xi-2 ， xi-1 ， xi 作為插值節(jié)點(diǎn)，即：

26、 i i L3（x）= (x- xj )/(xk-xj)f(xj -y(xj ) k=i-3 j=I-3jk 所以y(xj +1)= yi+1 = xi+1 L3（x）dx xi = yi +h/2455 yi+59 yi-1+37 yi-2-9 yi-3 外推公式 y(xj +1)- yi+1 = （25/720）h5yi（5）+ =0（h5）,計(jì)算f（x，y）的次數(shù)，算yi+1 時(shí) yi-1 ，yi-2， yi-3 已在計(jì)算yi時(shí)得到，故只要算yi就可，故無(wú)需算四次，故計(jì)算變化四階龍格-庫(kù)塔公式，而0（h5）相同，但需用其它方法先行計(jì)算y1 ， y2 ， y3 后才能用該法。依次計(jì)算y4

27、，y5 ，所以不具有自啟動(dòng)性，此外h要改變也很麻煩。 若將外推公式作預(yù)測(cè)式，內(nèi)插公式作校正，可得： 阿當(dāng)姆斯預(yù)測(cè)校正式 yj +1（0）= yi +h/2455 yi+59 yi-1+37 yi-2-9 yi-3 yj +1= yi +h/249f(xj +1 ， yj +1（0)+19 yi+5 yi-1- yi-2 ,例4： y=-2xy 0 x1.2 y(0)=1 取h=0.2 解：四階龍格-庫(kù)塔求出， y1， y2 ， y3 ，由 y0 ， y1 ， y2 ， y3阿當(dāng)姆斯預(yù)測(cè)校正公式。計(jì)算見p159 表7-6，7-7,5 一階方程組與高階方程,一 一階方程組 u=（x，u，v）, u(x 0)= u0 v=(x，u，v) , v(x 0)= v0 初值解 記y=（u，v）T，y0=（ u0 ， v0 ）T F=（ ， ）T y= F（ x ， y ） y(x 0)= y0 用龍格庫(kù)塔解, yi+1= yi + h/6（ k1 + 2k2 +2 k3 + k4 ） k1 = f（ xi ， yi ） k2 = f（xi +h/2， yi + h/