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

/*numericalmethodsforordinarydifferentialequations*/

考慮一階常微分方程的初值問題/*initial-valueproblem*/:只要f(x,y)在[a,b]r1上連續(xù)，且關(guān)于y

滿足lipschitz

條件，即存在與x,y無關(guān)的常數(shù)l

使對任意定義在[a,b]上的y1(x)和y2(x)都成立，則上述ivp存在唯一解。要計(jì)算出解函數(shù)y(x)在一系列節(jié)點(diǎn)a=x0<x1<…<xn=b

處的近似值節(jié)點(diǎn)間距為步長，通常采用等距節(jié)點(diǎn)，即取hi=h

(常數(shù))。第9章常微分方程數(shù)值解§1

引言§2歐拉方法/*euler’smethod*/

歐拉公式：x0x1向前差商近似導(dǎo)數(shù)記為定義在假設(shè)yi=y(xi)，即第

i

步計(jì)算是精確的前提下，考慮的截?cái)嗾`差ri=y(xi+1)

yi+1稱為局部截?cái)嗾`差/*localtruncationerror*/。定義若某算法的局部截?cái)嗾`差為o(hp+1)，則稱該算法有p

階精度。

歐拉法的局部截?cái)嗾`差：歐拉法具有1階精度。ri

的主項(xiàng)/*leadingterm*/亦稱為歐拉折線法

/*euler’spolygonalarcmethod*/

例1：解：h=0.2,xi=1+ih，精確解為：y=x2-2xy1=y0+hf(x0,y0)=-1y2=y1+hf(x1,y1)=-0.9333

y3=y2+hf(x2,y2)=-0.8y4=y3+hf(x3,y3)=-0.6y5=y4+hf(x4,y4)=-0.3333y6=y5+hf(x5,y5)=0可以看出誤差隨著計(jì)算在積累。xky(xk)ykek1.2-0.96-10.041.4-0.84-0.93330.09331.6-0.64-0.80.161.8-0.36-0.60.242.00-0.33330.33332.20.4400.44步長h=0.2§1euler’smethod隱式歐拉法/*impliciteulermethod*/向后差商近似導(dǎo)數(shù)x0x1))(,()(1101xyxfhyxy+)1,...,0(),(111-=+=+++niyxfhyyiiii由于未知數(shù)yi+1

同時出現(xiàn)在等式的兩邊，不能直接得到，故稱為隱式/*implicit*/

歐拉公式，而前者稱為顯式/*explicit*/歐拉公式。一般先用顯式計(jì)算一個初值，再迭代求解。隱式歐拉法的局部截?cái)嗾`差：即隱式歐拉公式具有1階精度。

歐拉公式的改進(jìn)：§1euler’smethod/*trapezoidformula*/—顯/隱式兩種算法的平均注：的確有局部截?cái)嗾`差，即梯形公式具有2

階精度，比歐拉方法有了進(jìn)步。但注意到該公式是隱式公式，計(jì)算時不得不用到迭代法。中點(diǎn)歐拉公式/*midpointformula*/中心差商近似導(dǎo)數(shù)x0x2x1假設(shè)，則可以導(dǎo)出即中點(diǎn)公式具有2階精度。需要2個初值y0和y1來啟動遞推過程，這樣的算法稱為雙步法/*double-stepmethod*/，而前面的三種算法都是單步法/*single-stepmethod*/。梯形公式方法§1euler’smethod顯式歐拉隱式歐拉梯形公式中點(diǎn)公式簡單精度低穩(wěn)定性最好精度低,計(jì)算量大精度提高計(jì)算量大精度提高,顯式多一個初值,可能影響精度can’tyougivemeaformulawithalltheadvantagesyetwithoutanyofthedisadvantages?doyouthinkitpossible?well,callmegreedy…ok,let’smakeitpossible./*modifiedeuler’smethod*/step1:

先用顯式歐拉公式作預(yù)測，算出),(1iiiiyxfhyy+=+step2:再將代入隱式梯形公式的右邊作校正，得到1+iy)],(),([2111+++++=iiiiiiyxfyxfhyy注：此法亦稱為預(yù)測-校正法/*predictor-correctormethod*/?？梢宰C明該算法具有2階精度，同時可以看到它是個單步遞推格式，比隱式公式的迭代求解過程簡單。后面將看到，它的穩(wěn)定性高于顯式歐拉法?！?euler’smethod改進(jìn)歐拉法例2：解：梯形公式為：可以預(yù)測校正方法來求：解yk+1出來比較困難，遇到的是一個二次方程，即：xkyky(xk)ek0.11.09596911.09544510.0004640.21.18409661.18321600.0008810.31.26620141.26491110.0012900.41.34336021.34164080.0017190.51.41640191.41421360.0021880.61.47595561.48323970.0027160.71.55251411.54919330.0033210.81.61647481.61245150.0040230.91.67816641.67332010.0047461.01.73786741.73205080.005817/*runge-kuttamethod*/建立高精度的單步遞推格式。單步遞推法的基本思想是從(xi,yi)點(diǎn)出發(fā)，以某一斜率沿直線達(dá)到(xi+1

,yi+1

)點(diǎn)。歐拉法及其各種變形所能達(dá)到的最高精度為2階。

考察改進(jìn)的歐拉法，可以將其改寫為：斜率一定取k1k2的平均值嗎？步長一定是一個h

嗎？§3龍格-庫塔法§2runge-kuttamethod首先希望能確定系數(shù)1、2、p，使得到的算法格式有2階精度，即在的前提假設(shè)下，使得

step1:將k2在(xi,yi)

點(diǎn)作taylor展開將改進(jìn)歐拉法推廣為：),(),(][12122111phkyphxfkyxfkkkhyyiiiiii++==++=+llstep2:將k2代入第1式，得到§2runge-kuttamethodstep3:將yi+1與y(xi+1)在xi點(diǎn)的泰勒展開作比較要求，則必須有：這里有個未知數(shù)，個方程。32存在無窮多個解。所有滿足上式的格式統(tǒng)稱為2階龍格-庫塔格式。注意到，就是改進(jìn)的歐拉法。q:

為獲得更高的精度，應(yīng)該如何進(jìn)一步推廣？其中i

(i=1,…,m)，i

(i=2,…,m)

和ij

(i=2,…,m;j=1,…,i1

)

均為待定系數(shù)，確定這些系數(shù)的步驟與前面相似?！?runge-kuttamethod)...,(......),(),(),(]...[1122112321313312122122111--++++++=+++=++==++++=mmmmmmimiiiiiimmiihkhkhkyihxfkhkhkyhxfkhkyhxfkyxfkkkkhyybbbabbaballl

最常用為四級4階經(jīng)典龍格-庫塔法

/*classicalrunge-kuttamethod*/

：例3：其精確解為：用經(jīng)典龍格－庫塔方法，h=0.1(1)求y1，此時解：以下計(jì)算用表格列出：

1.11.20027141.20027111.21.40181581.40181531.31.60523931.60523871.41.81079361.71079301.52.01856282.01856211.62.22854702.22854631.72.44070382.44070301.82.65496932.65496851.92.87126972.87126892.03.08952793.0895271附：三階龍格－庫塔方法:共有八個參數(shù)：1,

2,

3,

1,

2,

21,

31,

32，1、三階龍格－庫塔公式格式如下：