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

Chapter6

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

§6.0Introduction§6.1歐拉方法

(Euler’sMethod)§6.2龍格-庫(kù)塔法

(Runge-KuttaMethod)§6.3收斂性與穩(wěn)定性

(ConvergencyandStability)§6.4線性多步法

(MultistepMethods)§6.5微分方程組與高階方程

(SystemsofDifferentialEquationsandHigher-OrderEquations)

§6.6邊值問(wèn)題的數(shù)值解

(Boundary-ValueProblems)

考慮一階常微分方程的初值問(wèn)題/*Initial-ValueProblem*/:只要f(x,y)在[a,b]

R1上連續(xù)，且關(guān)于y滿足Lipschitz

條件，即存在與x,y無(wú)關(guān)的常數(shù)L使對(duì)任意定義在[a,b]上的y1(x)和y2(x)都成立，則上述IVP存在唯一解。要計(jì)算出解函數(shù)y(x)在一系列節(jié)點(diǎn)a=x0<x1<…<xn=b

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

(常數(shù))?！?.0Introduction

歐拉公式：x0x1向前差商近似導(dǎo)數(shù)記為定義在假設(shè)yi=y(xi)，即第i步計(jì)算是精確的前提下，考慮的截?cái)嗾`差Ri=y(xi+1)

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

歐拉法的局部截?cái)嗾`差：歐拉法具有1階精度。Ri的主項(xiàng)/*leadingterm*/亦稱(chēng)為歐拉折線法

/*Euler’spolygonalarcmethod*/

§6.1歐拉方法

(Euler’sMethod)例6.1

用歐拉法求初值問(wèn)題當(dāng)h=0.02時(shí)在區(qū)間[0,0.10]上的數(shù)值解。方程真解：〔參考程序〕nxnyny(xn)

n=y(xn)-

yn001.00001.0000010.020.98200.98250.000520.040.96500.96600.001030.060.94890.95030.001440.080.93370.93540.001750.100.91920.92120.0020歐拉公式的變形與改進(jìn)：

隱式歐拉法/*implicitEulermethod*/向后差商近似導(dǎo)數(shù)x0x1))(,()(1101xyxfhyxy+

)1,...,0(),(111-=+=+++niyxfhyyiiii由于未知數(shù)yi+1

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

歐拉公式，而前者稱(chēng)為顯式/*explicit*/歐拉公式。一般先用顯式計(jì)算一個(gè)初值，再迭代求解。

隱式歐拉法的局部截?cái)嗾`差：即隱式歐拉公式具有1階精度。Hey!Isn’ttheleadingtermofthelocaltruncationerrorofEuler’smethod?Seemsthatwecanmakeagooduseofit…§6.1歐拉方法

(Euler’sMethod)

梯形公式/*trapezoidformula*/—顯、隱式兩種算法的平均注：的確有局部截?cái)嗾`差，即梯形公式具有2階精度，比歐拉方法有了進(jìn)步。但注意到該公式是隱式公式，計(jì)算時(shí)不得不用到迭代法，其迭代收斂性與歐拉公式相似。

中點(diǎn)歐拉公式/*midpointformula*/中心差商近似導(dǎo)數(shù)x0x2x1假設(shè)，則可以導(dǎo)出即中點(diǎn)公式具有2階精度。需要2個(gè)初值y0和y1來(lái)啟動(dòng)遞推過(guò)程，這樣的算法稱(chēng)為雙步法/*double-stepmethod*/，而前面的三種算法都是單步法/*single-stepmethod*/?！?.1歐拉方法

(Euler’sMethod)方法

顯式歐拉隱式歐拉梯形公式中點(diǎn)公式簡(jiǎn)單精度低穩(wěn)定性最好精度低,計(jì)算量大精度提高計(jì)算量大精度提高,顯式多一個(gè)初值,可能影響精度Can’tyougivemeaformulawithalltheadvantagesyetwithoutanyofthedisadvantages?Doyouthinkitpossible?Well,callmegreedy…OK,let’smakeitpossible.§6.1歐拉方法

(Euler’sMethod)改進(jìn)歐拉法/*modifiedEuler’smethod*/Step1:先用顯式歐拉公式作預(yù)測(cè)，算出),(1iiiiyxfhyy+=+Step2:再將代入隱式梯形公式的右邊作校正，得到1+iy)],(),([2111+++++=iiiiiiyxfyxfhyy注：此法亦稱(chēng)為預(yù)測(cè)-校正法/*predictor-correctormethod*/?？梢宰C明該算法具有2階精度，同時(shí)可以看到它是個(gè)單步遞推格式，比隱式公式的迭代求解過(guò)程簡(jiǎn)單。后面將看到，它的穩(wěn)定性高于顯式歐拉法。§6.1歐拉方法

(Euler’sMethod)歐拉法的截?cái)嗾`差：

改進(jìn)歐拉法的截?cái)嗾`差：【例6.2】

在區(qū)間[0,1.5]上，取h=0.1，求解〔參考程序〕。本題的精確解為，可用來(lái)檢驗(yàn)近似解的精確程度。計(jì)算結(jié)果如下表：

xn

歐拉法yn迭代一次改進(jìn)歐拉法yn準(zhǔn)確解01110.11.11.0959091.0954450.21.1918181.1840961.1832160.31.2774381.2602011.2649110.41.3582131.3433601.3416410.51.4351331.4161021.4142140.61.5089661.4829561.4832400.71.5803381.5525151.5491930.81.6497831.6164761.6124520.91.7177791.6781681.6733201.01.7847701.7378691.7320511.11.851181.7958221.7888541.21.9174641.8522421.8439091.31.9840461.9073231.8973671.42.0514041.9612531.9493591.52.1200522.0142072.000000建立高精度的單步遞推格式。單步遞推法的基本思想是從(xi,yi)點(diǎn)出發(fā)，以某一斜率沿直線達(dá)到(xi+1

,yi+1

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

考察改進(jìn)的歐拉法，可以將其改寫(xiě)為：斜率一定取K1K2的平均值嗎？步長(zhǎng)一定是一個(gè)h嗎？§6.2龍格-庫(kù)塔法

(Runge-KuttaMethod)首先希望能確定系數(shù)

1、

2、p，使得到的算法格式有2階精度，即在的前提假設(shè)下，使得

Step1:將K2在(xi,yi)點(diǎn)作Taylor展開(kāi)將改進(jìn)歐拉法推廣為：),(),(][12122111phKyphxfKyxfKKKhyyiiiiii++==++=+llStep2:將K2代入第1式，得到§6.2Runge-KuttaMethodsStep3:將yi+1與y(xi+1)在xi點(diǎn)的泰勒展開(kāi)作比較要求，則必須有：這里有個(gè)未知數(shù)，個(gè)方程。32存在無(wú)窮多個(gè)解。所有滿足上式的格式統(tǒng)稱(chēng)為2階龍格-庫(kù)塔格式。注意到，就是改進(jìn)的歐拉法。Q:為獲得更高的精度，應(yīng)該如何進(jìn)一步推廣？§6.2Runge-KuttaMethods其中

i

(i=1,…,m)，

i

(i=2,…,m)

和

ij

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

)

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

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

/*ClassicalRunge-KuttaMethod*/：§6.2Runge-KuttaMethods注：

龍格-庫(kù)塔法的主要運(yùn)算在于計(jì)算Ki

的值，即計(jì)算f

的值。Butcher于1965年給出了計(jì)算量與可達(dá)到的最高精度階數(shù)的關(guān)系：753可達(dá)到的最高精度642每步須算Ki的個(gè)數(shù)

由于龍格-庫(kù)塔法的導(dǎo)出基于泰勒展開(kāi)，故精度主要受解函數(shù)的光滑性影響。對(duì)于光滑性不太好的解，最好采用低階算法而將步長(zhǎng)h

取小。HW:p.127#1,2,3,5§6.2Runge-KuttaMethods【例6.3】

用Runge-Kutta(龍格―庫(kù)塔)方法解初值問(wèn)題y’=x2–y(0≤x≤1)y(0)=1.〔參考程序〕 解：

取h=0.1，

【作業(yè)】用Runge-Kutta方法解例6.2，與那里的結(jié)果比較。

收斂性/*Convergency*/定義若某算法對(duì)于任意固定的x=xi=x0+ih，當(dāng)h0

(同時(shí)i

)時(shí)有yi

y(xi

)，則稱(chēng)該算法是收斂的。例：就初值問(wèn)題考察歐拉顯式格式的收斂性。解：該問(wèn)題的精確解為歐拉公式為對(duì)任意固定的x=xi=ih，有

§6.3收斂性與穩(wěn)定性

(ConvergenceandStability)

穩(wěn)定性/*Stability*/例：考察初值問(wèn)題在區(qū)間[0,0.5]上的解。分別用歐拉顯、隱式格式和改進(jìn)的歐拉格式計(jì)算數(shù)值解。0.00.10.20.30.40.5精確解改進(jìn)歐拉法

歐拉隱式歐拉顯式

節(jié)點(diǎn)xi

1.0000

2.00004.0000

8.00001.6000101

3.2000101

1.00002.5000101

6.25001021.56251023.90631039.76561041.00002.50006.25001.56261013.90631019.76561011.00004.97871022.47881031.23411046.14421063.0590107Whatiswrong??!AnEngineercomplains:"Maththeoremsaresounstablethatasmallperturbationontheconditionswillcauseacrashontheconclusions!"§6.3ConvergenceandStability定義若某算法在計(jì)算過(guò)程中任一步產(chǎn)生的誤差在以后的計(jì)算中都逐步衰減，則稱(chēng)該算法是絕對(duì)穩(wěn)定的/*absolutelystable*/。一般分析時(shí)為簡(jiǎn)單起見(jiàn)，只考慮試驗(yàn)方程/*testequation*/常數(shù)，可以是復(fù)數(shù)當(dāng)步長(zhǎng)取為h時(shí)，將某算法應(yīng)用于上式，并假設(shè)只在初值產(chǎn)生誤差，則若此誤差以后逐步衰減，就稱(chēng)該算法相對(duì)于絕對(duì)穩(wěn)定，的全體構(gòu)成絕對(duì)穩(wěn)定區(qū)域。我們稱(chēng)算法A比算法B穩(wěn)定，就是指A的絕對(duì)穩(wěn)定區(qū)域比B的大。hlh=h§6.3ConvergenceandStability例：考察顯式歐拉法由此可見(jiàn)，要保證初始誤差

0以后逐步衰減，必須滿足：0-1-2ReImg例：考察隱式歐拉法可見(jiàn)絕對(duì)穩(wěn)定區(qū)域?yàn)椋?10ReImg注：一般來(lái)說(shuō)，隱式歐拉法的絕對(duì)穩(wěn)定性比同階的顯式法的好。§6.3ConvergenceandStability例：隱式龍格-庫(kù)塔法而顯式1~4階方法的絕對(duì)穩(wěn)定區(qū)域?yàn)槠渲?階方法的絕對(duì)穩(wěn)定區(qū)域?yàn)?ReImgk=1k=2k=3k=4-1-2-3---123ReImg無(wú)條件穩(wěn)定§6.3ConvergenceandStability用若干節(jié)點(diǎn)處的y及y’值的線性組合來(lái)近似y(xi+1)。)...(...110111101kikiiikikiiiffffhyyyy--+---+++++++++=bbbbaaa其通式可寫(xiě)為：

基于數(shù)值積分的構(gòu)造法將在上積分，得到只要近似地算出右邊的積分，則可通過(guò)近似y(xi+1)。而選用不同近似式Ik，可得到不同的計(jì)算公式?！?.4線性多步法

(MultistepMethods)當(dāng)10時(shí)，為隱式公式;1=0那么為顯式公式。

亞當(dāng)姆斯顯式公式/*Adamsexplicitformulae*/利用k+1個(gè)節(jié)點(diǎn)上的被積函數(shù)值構(gòu)造k階牛頓后插多項(xiàng)式,有Newton插值余項(xiàng)/*顯式計(jì)算公式*/局部截?cái)嗾`差為：例：k=1時(shí)有§6.4MultistepMethods注：一般有，其中Bk與yi+1計(jì)算公式中fi,…,fi

k

各項(xiàng)的系數(shù)均可查表得到。10123kfifi1fi2fi3…Bk…………………Misprintonp.204常用的是k=3的4階亞當(dāng)姆斯顯式公式§6.4MultistepMethods

亞當(dāng)姆斯隱式公式/*Adamsimplicitformulae*/利用k+1個(gè)節(jié)點(diǎn)上的被積函數(shù)值fi+1

,fi,…,fi

k+1

構(gòu)造k階牛頓前插多項(xiàng)式。與顯式多項(xiàng)式完全類(lèi)似地可得到一系列隱式公式，并有，其中與fi+1

,fi,…,fi

k+1的系數(shù)亦可查表得到。~~10123kfi+1fifi1fi2…Bk…………………~常用的是k=3的4階亞當(dāng)姆斯隱式公式小于Bk較同階顯式穩(wěn)定§6.4MultistepMethods

亞當(dāng)姆斯預(yù)測(cè)-校正系統(tǒng)/*Adamspredictor-correctorsystem*/Step1:用Runge-Kutta法計(jì)算前k

個(gè)初值；Step2:用Adams顯式計(jì)算預(yù)測(cè)值；Step3:用同階Adams隱式計(jì)算校正值。注意：三步所用公式的精度必須相同。通常用經(jīng)典Runge-Kutta法配合4階Adams公式。Hey!LookatthelocaltruncationerroroftheexplicitandimplicitAdamsmethods:andDon’tyouthinkthere’ssomethingyoucando?4階Adams隱式公式的截?cái)嗾`差為4階Adams顯式公式的截?cái)嗾`差為當(dāng)h充分小時(shí)，可近似認(rèn)為

i

i

，則：Predictedvaluepi+1Modifiedvaluemi+1Correctedvalueci+1Modifiedfinalvalueyi+1外推技術(shù)/*extrapolation*/§6.4MultistepMethods

Adams4th-Orderpredictor-correctorAlgorithmToapproximatethethesolutionoftheinitial-valueproblemAt(N+1)equallyspacednumbersintheinterval[a,b].Input:endpointsa,b;integerN;initialvaluey0

.Output:approximationyatthe(N+1)valuesofx.Step1Seth=(b

a)/N;x0=a;y0=y0;Output(x0,y0);Step2Fori=1,2,3ComputeyiusingclassicalRunge-Kuttamethod;Output(xi,yi);Step3Fori=4,…,Ndosteps4-10

Step5;/*predict*/

Step6;/*modify*/

Step7

;/*correct*/

Step8;/*modifythefinalvalue*/

Step9

Output(xi+1

,yi+1

);

Step10Forj=0,1,2,3Setxi

=xi+1

;yi

=yi+1

;/*Preparefornextiteration*/Step11STOP.應(yīng)為(ci+1

pi+1

),但因ci+1

尚未算出，只好用(ci

pi)取代之?！?.4MultistepMethods

基于泰勒展開(kāi)的構(gòu)造法)...(...110111101kikiiikikiiiffffhyyyy--+---+++++++++=bbbbaaa

將通式中的右端各項(xiàng)yi1,…,yi

k;fi+1,fi1,…,fi

k

分別在

xi點(diǎn)作泰勒展開(kāi)，與精確解

y(xi+1)在xi點(diǎn)的泰勒展開(kāi)作比較。通過(guò)令同類(lèi)項(xiàng)系數(shù)相等，得到足以確定待定系數(shù)

0,…,

k;

1,

0,…,

k

的等式，則可構(gòu)造出線性多步法的公式?！?.4MultistepMethods例：設(shè))(3322110221101-----+

+

+

+

+++=iiiiiiiiyyyyhyyyybbbbaaa確定式中待定系數(shù)

0,

1,

2,

0,

1,

2,

3,

使得公式具有4階精度。解：/*y(xi)=yi*/個(gè)未知數(shù)個(gè)方程75§6.4MultistepMethods

令

1=

2=0Adams

顯式公式

以y

i+1取代y

i1，并取

1=

2=0Adams

隱式公式

以yi3取代y

i3，則可導(dǎo)出另一組4階顯式算法，其中包含了著名的米爾尼

/*Milne*/公式其局部截?cái)嗾`差為注：上式也可通過(guò)數(shù)值積分導(dǎo)出，即將在區(qū)間上積分，得到再過(guò)做f的插值多項(xiàng)式即可。取

1=1,

2=0得到辛甫生

/*Simpson*/公式與Milne公式匹配使用

辛甫生

/*Simpson*/公式在區(qū)間[xi1,xi+1]上積分，并用Simpson數(shù)值積分公式來(lái)近似積分項(xiàng)，亦可得此Simpson公式。§6.4MultistepMethods

Milne-Simpson系統(tǒng)的缺點(diǎn)是穩(wěn)定性差，為改善穩(wěn)定性，考慮另一種隱式校正公式：要求公式具有4階精度。通過(guò)泰勒展開(kāi)，可得到個(gè)等式，從中解出個(gè)未知數(shù)，那么有個(gè)自由度。561取

1=1得Simpson公式哈明

/*Hamming*/用

1的不同數(shù)值進(jìn)行試驗(yàn)，發(fā)現(xiàn)當(dāng)

1=0

時(shí)，公式的穩(wěn)定性較好，即：其局部截?cái)嗾`差為注：Hamming公式不能用數(shù)值積分方法推導(dǎo)出來(lái)。§6.4MultistepM