常微分方程初值問題數(shù)值解的實現(xiàn)和分析—四階Rungekutta方法與預估校正算法畢業(yè)論文

1、. . . . 1 / 22數(shù)值分析課程設計常微分方程初值問題數(shù)值解的實現(xiàn)和分析四階 Runge-kutta 方法與預估-校正算法. . . . I / 22常微分方程初值問題數(shù)值解的實現(xiàn)和分析四階 Runge-kutta 方法與預估-校正算法摘 要求解常微分方程的初值問題，Euler 方法，改進的 Euler 方法與梯形方法精度比較低，所以本文構造高精度單步的四級 Runge-kutta 方法與高精度的多步預估校正算法與其 Matlab 編程來實現(xiàn)對常微分方程初值問題的求解，使在求解常微分方程時，對以前積分方法的收斂速度與精度都有了很高的提高。關鍵詞：Runge-kutta 方法，Adams

2、 方法，預估校正算法，Matlab目 錄1.前言 12. 幾個簡單的數(shù)值積分法 12.1 Runge-kutta 方法 12.1.1 Runge-kutta 方法的應用 52.2 預估校正算法 72.2.1 Adams 數(shù)值積分方法簡介與預估校正算法 72.2.2 預估校正算法的應用 123. 結果分析 15總結 16參考文獻 17英文原文和中文翻譯 181 英文原文 182 中文翻譯 19. . . . 1 / 221.前言常微分方程的初值問題是微分方程定解問題的一個典型代表，以下面的例子介紹常微分方程初值問題數(shù)值解的基本思想和原理。例 1.1 一重量垂直作用于彈簧所引起的震蕩，當運動阻力與

3、速度的平方成正 比時，可借助如下二階常微分方程描述若令和，則上述二階常微分方程可化成等價的一階常微分方程組類似于例 1.1，對于 m階常微分方程其中。 若定義可得如下等價的一階常微分方程組我們知道多數(shù)常微分方程主要靠數(shù)值解法。所謂數(shù)值解法，就是尋求解在一系列離散節(jié)點上的近似值。相鄰兩個節(jié)點之間的間距稱為步長1。2. 幾個簡單的數(shù)值積分法2.1 Runge-kutta 方法 Runge 在 1985 年提出了一種基于 Euler 折線法的新的數(shù)值方法，此后這種新的數(shù)值. . . . 2 / 22方法又經(jīng)過其同胞 K.Heun 和 Kutta 的努力2，發(fā)展完善成為后世所稱的 Runge-kutt

4、a方法。 Runge 重要發(fā)現(xiàn)的靈感主要來源于把 Euler 方法應用于初值問題 (2.1.1)和 f 不顯含 y 時的初值問題 (2.1.2)之間的類比。Runge 觀察到，當把 Euler 方法應用到(2.1.2)型初值問題時，得到差分方程事實上，這是在計算積分問題 (2.1.3)時采用了左矩形法則即用高為，寬為 h 的矩形代替了(2.1.3)式中的積分值。類似的，當把 Euler 方法應用到初值問題(2.1.1)時，得到差分方程這是在計算積分問題 (2.1.4)時采用了左矩形法則. . . . 3 / 22顯然(2.1.3)和(2.1.4)式右側數(shù)值積分的精度決定著數(shù)值解的精度。此時，R

5、unge 發(fā)現(xiàn)，這個矩形公式的逼近程度并不高。如果在計算積分問題(2.1.4)時采用中點法則或梯形法則，則數(shù)值解的精度會更高。采用中點法則和梯形法則分別代替(2.1.4)式中的積分會得到 (2.1.5) (2.1.6)(2.1.5)和(2.1.6)是兩個隱式的方程，它們在早期發(fā)現(xiàn) Runge-kutta 方法過程中扮演了一個重要角色。Runge 的出發(fā)點是：分別用一個 Euler 步代替(2.1.5)中未知的和(2.1.6)式右側未知的值，這樣 Runge 得到如下的中點格式：和梯形格式：這兩個方程具有很好的幾何解釋。它們是有折線構成的，這些折線假定微分方程所確定的斜率在前面的點上已經(jīng)計算出來

6、了。與他的后繼者一樣，Runge 用 Taylor 展開和系數(shù)對比說明上述兩個方法的階為 2. Runge 意識到，用中點法和梯形法得出的兩個方法的階數(shù)還不夠高，所以，他設想，如果使用比中點法和梯形法精度更高的Simpson 法則得到的方法階數(shù)會提高，如果用 M 和 T 分別表示用中點法和梯形法算得的數(shù)值積分。Simpson 法則可以寫成眾所周知的形式 S=M+(T-M)/3。Taylor 展開表明，如果 f 依賴于. . . . 4 / 22y，則這個表達形式只是 2 階的，接著 Runge 發(fā)現(xiàn)，通過重復使用 Euler 步驟對梯形法則做適當?shù)恼{整，會使格式 =M+( -M)/3 成為 3

7、 階方法，他還把他的方法與其 Taylor展開式拓展到微分方程組3。五年后，Heun 在其 1900 年的文章中評論 Runge 的方法時說，Runge 獲得的上述方法是歸納性的而且是令人費解的，他主使用更具一般性的 Gauss 求積公式其中于是可以把一般的 Gauss 格式擴充為 (2.1.7)把(2.1.7)式的右端進行二元 Taylor 展開并與的 Taylor 展開式的對應的系數(shù)比較，適當選取參數(shù)使方法具有盡可能高的精度。上述算法公式中的系數(shù)希望如此確定，使得(2.1.7)的 Taylor 展開式中所有 h 冪次不超過 p 的那些項與中的相應項的系數(shù)相等2。假定 p=4 并取 s=3，

8、用類似的推導，可建立下述常用的顯型四級 Runge-kutta 方法：. . . . 5 / 22截斷誤差為，當右端函數(shù) f 不依賴于 y 時，上述公式可簡化為2.1.1 Runge-kutta 方法的應用 用四級 Runge-kutta 方法求初值問題 =y2 e-x ，y(1)=1 在區(qū)間1,2上的數(shù)值解(取 h=0.1) (1)先求出常微分方程的精確解，再用四階 Runge-kutta 方法進行求解。Y1=dsolve(Dy=y2*exp(-x),y(1)=1,x) 運行后的精確解為：Y1=1/(1/exp(x) - 1/exp(1) + 1)現(xiàn)給出常用的四階 Runge-kutta 方

9、法求解常微分方程的 Matlab 主程序如下：function X,Y=Rungek(funfcn,x0,b,y0,h)x=x0;y=y0;n=fix(b-x0)/h); %求步數(shù)i=1;X=zeros(n,1);Y=zeros(n,1);X(i)=x0;Y(i)=y0; %賦初值for i=2:n k1=feval(funfcn,x,y); k2=feval(funfcn,x+h/2,y+h*k1/2); k3=feval(funfcn,x+h/2,y+h*k2/2); k4=feval(funfcn,x+h,y+h*k3); y=y+h/6*(k1+2*k2+2*k3+k4); Y(i)=

10、y; x=x+h;. . . . 6 / 22 X(i)=x;end %按照龍格-庫塔方法進行求解X,Y1=1./(1./exp(X) 1./exp(1) + 1), %精確解plot(X,Y,mh,X,Y1,bo) %繪圖gridlegend(用四階龍哥庫塔方法計算的數(shù)值解,精確解 y(x), %圖形說明wcha=abs(Y-Y1), %截斷誤差自定義函數(shù)：function z=funfcn(x,y)z=y.2.*exp(-x);在 matlab 工作窗口輸入如下程序：x0=1,b=2;y0=1;h=0.1;X,Y=Rungek(funfcn,x0,b,y0,h)運行后得：X = 1.000

11、0 1.1000 1.2000 1.3000 1.4000 1.5000 1.6000 1.7000 1.8000 1.9000Y1 = 1.0000 1.0363 1.0714 1.1054 1.1380 1.1692 1.1990. . . . 7 / 22 1.2273 1.2540 1.2793wcha = 1.0e-007 * 0 0.0111 0.0290 0.0518 0.0777 0.1054 0.1338 0.1620 0.18960.215911.11.21.31.41.51.61.71.81.9211.051.11.151.21.251.31.35 值 值 值 值 值 值

12、 值 值 值 值 值 值 值 值 值值 值 值 y=(x)2.2 預估校正算法2.2.1 Adams 數(shù)值積分方法簡介與預估校正算法 構造高階精度格式的一個重要途徑是，將方程改寫成如下的積分形式 (2.2.1.1). . . . 8 / 22然后以適當?shù)牟钪刀囗検浇铺娲鲜街械谋环e函數(shù)，這時后面將要介紹的Adams 差值方法的出發(fā)點和思想。選取 k+1 個插值結點：。假定 y(x)在上述結點上的近似值已被求出，記為，并令，它被視為在 x=上的近似值。利用上述結點和近似值作 k 次插值多項式用它近似(2.2.1.1)中的被積函數(shù)，得到近似公式 (2.2.1.2)將由此公式確定的當做的近似值。在

13、上述插值中，被插值點 x位于包含所有插值結點的最小區(qū)間的外部，故稱(2)為外插公式。容易看出，當 k=0 時，(2.2.1.2)就是向前 Euler 公式。為了得到(2.2.1.2)的具體計算公式，令x=，并將在區(qū)間寫成 Newton 向后插值公式的形式其中代表 j 階向前差分算子。引進記號則有 (2.2.1.3)將表達式(2.2.1.3)代入(2.2.1.2)得到. . . . 9 / 22 (2.2.1.4)其中系數(shù) (2.2.1.5)它們與 k 和 n 無關。例如，由于所以外插方法(2.2.1.4)和(2.2.1.5)的截斷誤差為 (2.2.1.6)已知差分可用函數(shù)值的線性組合表示，故又

14、可將外插公式(2.2.1.4)改寫成 (2.2.1.7)其中系數(shù)已造表，應用時可查。Adams 外插公式表K公式0,1. . . . 10 / 2223仍然從積分關系式(2.2.1.1)出發(fā)，將其中的被積函數(shù)改用以為插值結點的 k 次 Lagrange 插值多項式近似。這里，由于被插值點 x位于包含所有插值結點的最小區(qū)間的部，故稱相應的積分公式為插公式1。用類似的方法可得到 Adams 插公式表如下Adams 插公式表K公式0123在求解常微分方程初值問題的數(shù)值積分法中，還有一類高階精度方法是多步方法。一般來說，隱式多步方法（如 Adams 插法）比相應的顯式多步方法精度高，并有較好的穩(wěn)定性質

15、。然而，隱式方法計算的每一步需要求解非線性方程，所以它的計算比顯式多步方法要復雜得多。為了克服隱式多步方法計算方面的困難，人們經(jīng)常隱式多步法和顯式多步法結合起來使用，即預估校正的算法。下面介紹預估和校正的具體做法和計算公式：隱式 k 步方法可寫成如下形式： (2.2.1.8)其中，可認為是已知的，只是關于的非. . . . 11 / 22線性方程。通常采用迭代法求解，計算公式如下：(2.2.1.9)容易證明，當 h 適當小并且恰當?shù)剡x取初值時，上述迭代過程是收斂的。顯而易見，隱式方法(2.2.1.8)的每一步的計算量由(2.2.1.9)的迭代次數(shù)決定，因此選取好的初始值是十分重要的。一種非常自

16、然的做法就是取為用相應的k 步外插法計算的值，即令 (2.2.1.10)由(2.2.1.10)和(2.2.1.9)構成的算法被稱為預估校正算法，其中(2.2.1.10)為預估公式(記為 P)，(2.2.1.9)為校正公式(記為 C)。由于這里的預估值比較精確，因此校正中所需的迭代次數(shù)一般不會很多。如果出現(xiàn)迭代收斂慢的情況，應縮小步長在進行計算。預估校正算法的計算步驟：首先用預估公式(2.2.1.10)得到，然后計算，緊接著使用校正公式(2.2.1.9)求出，這樣完成了一步校正。對重復上述過程，可得。如此循環(huán)，可得進行 N 次校正的預估校正算法的計算公式為P :E : . . . . 12 /

17、22C : (2.2.1.11)將上述過程簡記為 P(EC)N,這里 E 代表計算 f 的函數(shù)值。最后介紹一個提高近似解精度的有效方法：Richardson 外推法。用表示步長 h 按某一數(shù)值方法求得的近似解。假定該方法具有 p 階精度，此時將關于 h 展開有(即有) (2.2.1.12)這也就是說，。據(jù)此，當步長 h 縮小一半時，相應的近似解滿足 (2.2.1.13)將(2.2.1.13)乘以后減去(2.2.1.12)，可得由此可知，當按如下公式 (2.2.1.14)計算，并將所得作為的近似值，那么其精度要比至少要提高一階，即有。通常稱(2.2.1.14)為 Richardson 外推公式。

18、2.2.2 預估校正算法的應用 用下述預估校正方法 . . . . 13 / 22yn+4yn+3 =(三步四階 Adams 外插方法)yn+4yn+3 =(三步四階 Adams 插方法)計算出題(1)初值問題在區(qū)間1,2上的數(shù)值解(取 h=0.1)。現(xiàn)給出預估校正算法的主程序如下：function k,X,Y,wucha,P=dAdamsyx(fcn,x0,b,y0,h)x=x0;y=y0;p=128;n=fix(b-x0)/h); %求步數(shù)if n5 returnendX=zeros(p,1);Y=zeros(p,length(y);f=zeros(p,1);k=1;X(k)=x;Y(k,

19、:)=y;for k=2:4 c1=1/6;c2=2/6;c3=2/6;c4=1/6;a2=1/2; x1=x+h/2;x2=x+h/2;x3=x+h;k1=feval(fcn,x,y); y1=y+h*k1/2;x=x+h;k2=feval(fcn,x1,y1); y2=y+h*k2/2;k3=feval(fcn,x2,y2); y3=y+h*k3;x=x+h;k4=feval(fcn,x3,y3); y=y+h*(c1*k1+c2*k2+c3*k3+c4*k4);X(k)=x;Y(k,:)=y;end %四階 Runge-kutta 算法X,Y;f=feval(fcn,X(1:4),Y(1

20、:4);f=ffor k=4:n X(k+1)=X(1)+h*k;f(k)=feval(fcn,X(k),Y(k); p=Y(k)+(h/24)*(f(k-3:k)*-9 37 -59 55); f=f(2) f(3) f(4) feval(fcn,X(k+1),p), Y(k+1)=Y(k)+(h/24)*(f*1 -5 19 9); f(4)=feval(fcn,X(k+1),Y(k+1);k=k+1; end %預估校正算法for k=1:n. . . . 14 / 22wucha(k+1)=norm(Y(k+1)-Y(k); %每次迭代的截斷誤差endX=X(1:n+1);Y=Y(1:

21、n+1,:);n=1:n+1,wucha=wucha(1:n,:);P=n,X,Y,wucha;自定義函數(shù)：function z=fcn(x,y) z=y.2.*exp(-x);在 matlab 工作窗口輸入：x0=1;b=2;y0=1;h=0.1;subplot(3,1,1)k,X,Y1,wucha1,P1=dAdamsyx(fcn,x0,b,y0,h)y=1./(1./exp(X) - 1./exp(1) + 1);plot(X,Y1,mh,X,y,bo),gridlegend(用 adams 預測-校正公式計算的數(shù)值解,精確解 y=1/(1/exp(x) - 1/exp(1) + 1)w

22、uch1=abs(y-Y1),P1,y,wuch1title(dy/dx=y2*exp(-x),y(1)=1 在1,2的數(shù)值解和精確解的圖形)運行后得到：k =10X = 1.0000 1.2000 1.4000 1.6000 1.4000 1.5000 1.6000 1.7000 1.8000 1.9000 2.0000. . . . 15 / 22Y1 = 1.0000 1.0363 1.0680 1.0955 1.1217 1.1533 1.1833 1.2113 1.2379 1.2629 1.2865wucha1 = 0 0.0363 0.0317 0.0275 0.0262 0.0

23、316 0.0300 0.0280 0.0266 0.0251 0.0235P1 = 1.0000 1.0000 1.0000 0 2.0000 1.2000 1.0363 0.0363 3.0000 1.4000 1.0680 0.0317 4.0000 1.6000 1.0955 0.0275 5.0000 1.4000 1.1217 0.0262 6.0000 1.5000 1.1533 0.0316 7.0000 1.6000 1.1833 0.0300 8.0000 1.7000 1.2113 0.0280 9.0000 1.8000 1.2379 0.0266 10.0000 1.

24、9000 1.2629 0.0251 11.0000 2.0000 1.2865 0.0235wuch1 = 0 0.0352 0.0700 0.1036 0.0163. . . . 16 / 22 0.0159 0.0157 0.0160 0.0162 0.01640.016511.11.21.31.41.51.61.71.81.9211.21.4dy/dx=y2*exp(-x),y(1)=1值 1,2值 值 值 值 值 值 值 值 值 值 值 值 adams值 值 -值 值 值 值 值 值 值 值 值 值值 值 值 y=1/(1/exp(x) - 1/exp(1) + 1)3. 結果分析由

25、四階 Runge-kutta 方法計算的結果和 Adams 預估校正算法計算的結果進行對比可知，前一種方法計算的數(shù)值解與精確解在數(shù)據(jù)點處吻合的相當好，后一種方法則沒有前一種方法精確。總 結這次的數(shù)值分析課程設計我感覺很算是有收獲的。此次課程設計不僅使我認識到了 Matlab 的強大用途與涉與的領域之廣，重要的是通過對 Runge-kutta 數(shù)值積分方法和 Adams 積分方法的相關研究，使我對數(shù)值積分方法又有了新的理解和更深的認識。綜合課程設計讓我對常微分方程初值問題的數(shù)值積分方法得到鞏固和進一步的加深，對已有知識有了更進一步的理解和認識。再次，雖然我在課程設計中碰到了很多. . . . 1

26、7 / 22的問題，但通過查閱相關書籍，資料，并通過自己鉆研，特別是經(jīng)過小組成員的討論和反復的糾正，我得到了很大的幫助，不僅給了我思路上的開闊，還讓我認識到了自己對以前所學知識的不足方面。而且，通過這次課程設計，我也發(fā)現(xiàn)了自身的很多不足之處，在以后的學習中，我會不斷的完善自我，不斷進取，使自己在此方面有一個好的發(fā)展。參考文獻1 黃明游,果忱.數(shù)值分析M.:高等教育,2008. 2 林立軍，郭松云.常微分方程數(shù)值解法-Runge-kutta 方法的歷史淺析，師大學學報(自然科學版)，第 26 卷第 2 期，2003. 3 黃鐸,蘭平,王鳳.數(shù)值分析M.:科學,2005.英文原文和中文翻譯1 英文

27、原文PLOT Linear plot. PLOT(X,Y) plots vector Y versus vector X. If X or Y is a matrix,then the vector is plotted versus the rows or columns of the matrix,whichever line up. If X is a scalar and Y is a vector, length(Y) disconnected points are plotted.PLOT(Y) plots the columns of Y versus their index.I

28、f Y is complex, PLOT(Y) is equivalent to PLOT(real(Y),imag(Y).In all other uses of PLOT, the imaginary part is ignored.Various line types, plot symbols and colors may be obtained with PLOT(X,Y,S) where S is a character string made from one elementfrom any or all the following 3 columns: b blue . poi

29、nt - solid g green o circle : dotted r red x x-mark -. dashdot c cyan + plus - dashed m magenta * star (none) no line y yellow s square k black d diamond v triangle (down) triangle (up) triangle (right) p pentagram h hexagram For example, PLOT(X,Y,c+:) plots a cyan dotted line with a plus at each da

30、ta point; PLOT(X,Y,bd) plots blue diamond at each data point but does not draw any line. PLOT(X1,Y1,S1,X2,Y2,S2,X3,Y3,S3,.) combines the plots defined bythe (X,Y,S) triples, where the Xs and Ys are vectors or matrices and the Ss are strings. For example, PLOT(X,Y,y-,X,Y,go) plots the data twice, wit

31、h asolid yellow line interpolating green circles at the data points. The PLOT command, if no color is specified, makes automatic use ofthe colors specified by the axes ColorOrder property. The defaultColorOrder is listed in the table above for color systems where thedefault is blue for one line, and

32、 for multiple lines, to cycle through the first six colors in the table. For monochrome systems, PLOT cycles over the axes LineStyleOrder property. If you do not specify a marker type, PLOT uses no marker. If you do not specify a line style, PLOT uses a solid line.PLOT(AX,.) plots into the axes with

33、 handle AX. PLOT returns a column vector of handles to lineseries objects, one handle per plotted line. The X,Y pairs, or X,Y,S triples, can be followed by parameter/value pairs to specify additional properties of the lines. For example, PLOT(X,Y,LineWidth,2,Color,.6 0 0) will create a plot with a dark red line width of 2 points.Backwards compatibility PLOT(v6,.) creates line objects instead of lineseries objects for compatibility with MATLAB 6.5 and earlier.See also plottools, semilogx, semilogy, loglog, plotyy, plot3, grid, title, xlabel, ylabel, axis, axes, hold, legend, subplot, scat