第九章 常微分方程的數(shù)值解法 ppt課件

1、第九章第九章 常微分方程數(shù)值解法常微分方程數(shù)值解法第一節(jié)第一節(jié) Euler方法方法第三節(jié)第三節(jié) 單步法的收斂性和穩(wěn)定性單步法的收斂性和穩(wěn)定性第二節(jié)第二節(jié) Runge-Kutta方法方法上一頁上一頁 下一頁下一頁 返回返回 上一頁上一頁 下一頁下一頁 返回返回 本章介紹求解微分方程數(shù)值解的基本思想和方法本章介紹求解微分方程數(shù)值解的基本思想和方法. 含有自變量、未知函數(shù)和它的一階導(dǎo)數(shù)和高階導(dǎo)數(shù)的方程. 常微分方程常微分方程 它是描述運動、變化規(guī)律的重要數(shù)學(xué)方法之一，它是描述運動、變化規(guī)律的重要數(shù)學(xué)方法之一，分為兩類：分為兩類：1. 初值問題初值問題 即給出未知函數(shù)及導(dǎo)數(shù)在初始點的值；即給出未知函數(shù)

2、及導(dǎo)數(shù)在初始點的值；2. 邊值問題邊值問題 即給出未知函數(shù)及或它的某些導(dǎo)數(shù)即給出未知函數(shù)及或它的某些導(dǎo)數(shù)在區(qū)間兩個端點的值在區(qū)間兩個端點的值 。 考慮一階常微分方程的初值問題考慮一階常微分方程的初值問題 :0)(,),(yaybaxyxfy只要只要 f (x, y) 在在a, b R1 上連續(xù)，且關(guān)于上連續(xù)，且關(guān)于 y 滿足滿足 Lipschitz 條件，即存在與條件，即存在與 x, y 無關(guān)的常數(shù)無關(guān)的常數(shù) L 使使對任意定義在對任意定義在 a, b 上的上的 y1(x) 和和 y2(x) 都成立，則上述問題都成立，則上述問題解存在唯一解。解存在唯一解。| ),(),(|2121yyLyxf

3、yxf 所謂數(shù)值解法就是要計算出初值問題的解函數(shù)所謂數(shù)值解法就是要計算出初值問題的解函數(shù) y(x) 在一系列在一系列離散點離散點 a = x0 x1 xN= b上的近似值：上的近似值： y0 , y1 ,yN .節(jié)點間距節(jié)點間距 為步長，為步長，) 1,., 0(1 nixxhiii通常采用等距節(jié)點，即取通常采用等距節(jié)點，即取 hi = h (常數(shù)常數(shù))。 yn 稱為問題的數(shù)值解稱為問題的數(shù)值解.數(shù)值解所滿足的離散方程統(tǒng)稱為差分格式數(shù)值解所滿足的離散方程統(tǒng)稱為差分格式. 上一頁上一頁 下一頁下一頁 返回返回 第一節(jié)第一節(jié) 歐拉方法歐拉方法一、一、 歐拉公式歐拉公式點點列列出出方方程程式式在在令

4、令等等份份分分成成將將求求解解區(qū)區(qū)間間nnxNnnhaxNabhNba. ), 1 , 0(,(*)(,()(nnnxyxfxy ,)()()()(1hxyxyxyxynnnn即即用用一一階階差差商商近近似似代代替替，將將令令yn為為y(xn)的近似值，將上式代入的近似值，將上式代入(*)式可得式可得1, 2 , 1 , 0 ),(1Nnyxhfyynnnn此式稱為歐拉此式稱為歐拉(Euler)公式公式.,210Nyyyy值值問問題題的的數(shù)數(shù)值值解解即即可可由由此此式式依依次次求求出出初初已已知知nnnyxyT)( 為Euler方法的局部截斷誤差.上一頁上一頁 下一頁下一頁 返回返回 例例1

5、用歐拉公式解初值問題用歐拉公式解初值問題 1)0(10,2yxyxyy解：解： 取步長取步長 h=0.1, 歐拉公式的具體形式為：歐拉公式的具體形式為：)2(1 . 0),(1nnnnnnnnyxyyyxfhyy由由此此式式可可得得已已知知其其中中, 1),10, 1 , 0(1 . 00ynnnhaxn 1 . 11 . 01)2(00001 yxyhyy191818. 1)1 . 12 . 01 . 1 ( 1 . 01 . 1)2(11112 yxyhyy依次計算可得依次計算可得 y3，y4，y5，y6，y7，y8，y9，y10上一頁上一頁 下一頁下一頁 返回返回 其部分結(jié)果見下表其部分

6、結(jié)果見下表 可見可見Euler方法的計算結(jié)果精度不太高。方法的計算結(jié)果精度不太高。 xy21 準(zhǔn)準(zhǔn)確確解解上一頁上一頁 下一頁下一頁 返回返回 出出發(fā)發(fā)的的一一條條曲曲線線表表示示從從的的解解方方程程),()()(),(0000yxPxyyaybxayxfy 歐拉公式的幾何意義：歐拉公式的幾何意義：x0P0 x1P1),)(,()(),(0000000000 xxyxfyxxxyyyyxP 的的切切線線方方程程過過點點,),(000101hyxfyyhxxx的的交交點點的的縱縱坐坐標(biāo)標(biāo)為為與與直直線線 x2P2xnPn2111212111111111,),(),)(,(),(),(yhyxfy

7、yhxxxxxyxfyyyxfyxP是是歐歐拉拉公公式式求求出出的的正正好好的的交交點點的的縱縱坐坐標(biāo)標(biāo)為為與與的的直直線線為為且且斜斜率率為為過過點點 .1y的的正正好好是是由由歐歐拉拉公公式式求求出出幾何意義：幾何意義：用折線近似代替方程的解曲線用折線近似代替方程的解曲線,因而也稱因而也稱Euler方法為折線法方法為折線法.上一頁上一頁 下一頁下一頁 返回返回 二、二、 后退的歐拉公式后退的歐拉公式)(,()(1111點點列列出出方方程程在在nnnnxyxfxyx 也用一階差商逼近導(dǎo)數(shù)也用一階差商逼近導(dǎo)數(shù)hxyxyxynnn)()()(11令令yn+1為為y(xn+1)的近似值，則可得的近

8、似值，則可得1, 2 , 1 , 0 ),(111Nnyxhfyynnnn稱為后退稱為后退Euler公式公式 知知 yn時時,必須通過解方程才能求出必須通過解方程才能求出yn1 ,這樣的公式稱為這樣的公式稱為隱式公式隱式公式, 而而Euler公式為顯式公式公式為顯式公式. Euler公式和后退公式和后退Euler公式都是由公式都是由yn去計算去計算yn+1 ，因而，稱，因而，稱它們?yōu)閱尾椒?。它們?yōu)閱尾椒āＩ弦豁撋弦豁?下一頁下一頁 返回返回 定義定義 在假設(shè)在假設(shè) yi = y(xi)，即第，即第 i 步計算是精確的前提下，考步計算是精確的前提下，考慮的截斷誤差慮的截斷誤差 Ti+1 = y(

9、xi+1) yi+1 稱為局部截斷誤差。稱為局部截斷誤差。定義定義 若某算法的局部截斷誤差為若某算法的局部截斷誤差為O(hp+1)，則稱該算法有，則稱該算法有p 階精度。階精度。顯然顯然, p越大越大, 精度越高精度越高.三、三、 局部截斷誤差與方法的階局部截斷誤差與方法的階（將準(zhǔn)確解（將準(zhǔn)確解)(xy代入公式的左、右兩端，其左端與右端之差）代入公式的左、右兩端，其左端與右端之差） Euler方法的精度方法的精度 )(,()()()(111nnnnnnnxyxhfxyhxyyxyT其中：其中：)()(,(nnnxyxyxf 上一頁上一頁 下一頁下一頁 返回返回 所以，所以， Euler方法具有

10、方法具有 1 階精度。階精度。),( , ! 2)()()()(12 nnnnnxxhyxyhxyhxy )(hxynnx將將在點在點處一階處一階TaylorTaylor展開展開)()()(2hOxyhxynn)()()()()()(221hOxyhxyhOxyhxyTnnnnn上一頁上一頁 下一頁下一頁 返回返回 所以，所以， 后退的后退的Euler方法也具有方法也具有 1 階精度。階精度。 , ! 2)()()()()(2111hyxyhxyhxyxynnnn )(nxy1nx將將在點在點處一階處一階TaylorTaylor展開展開 )()()(211hOxyhxynn 隱式隱式Euler

11、方法的精度方法的精度 )( )()()()()(2121111hOxyhhOxyhxyxyTnnnnn上一頁上一頁 下一頁下一頁 返回返回 顯、隱式兩種算法的平均顯、隱式兩種算法的平均),(),(2111 nnnnnnyxfyxfhyy 歐拉公式的改進歐拉公式的改進)()(1233hOxyhTn 其局部誤差為：其局部誤差為：此公式具有此公式具有2 2階精度階精度. .稱平均公式或梯形公式稱平均公式或梯形公式梯形公式可由下迭代式計算：其中迭代初值是梯形公式可由下迭代式計算：其中迭代初值是Euler公式提供公式提供. ),1 ,0(),(),(2),(111101kyxfyxfhyyyxhfyyk

12、nnnnnknnnnn上一頁上一頁 下一頁下一頁 返回返回 四、改進的歐拉公式四、改進的歐拉公式Step 1: 先用顯式歐拉公式作預(yù)測，算出先用顯式歐拉公式作預(yù)測，算出),(1iiiiyxfhyy Step 2: 再將再將 代入隱式梯形公式的右邊作校正，得到代入隱式梯形公式的右邊作校正，得到1 iy),(),(2111 iiiiiiyxfyxfhyy注：此法亦稱為預(yù)測注：此法亦稱為預(yù)測-校正法校正法 。 可以證明該算法具有可以證明該算法具有 2 階精度，同時可以看到它是個單階精度，同時可以看到它是個單步遞推格式，比隱式公式的迭代求解過程簡單。它的精步遞推格式，比隱式公式的迭代求解過程簡單。它的

13、精度高于顯式歐拉法。度高于顯式歐拉法。 )1,., 0(),(,),(211 niyxfhyxfyxfhyyiiiiiiii上一頁上一頁 下一頁下一頁 返回返回 為了便于編程為了便于編程, 常將改進的歐拉公式寫為：常將改進的歐拉公式寫為：cpnpnncnnnpyyyyxhfyyyxhfyy21),(),(11上一頁上一頁 下一頁下一頁 返回返回 例例2 用改進的歐拉法解例用改進的歐拉法解例1中的初值問題中的初值問題.1)0(10,2yxyxyy解：取步長解：取步長 h=0.1, 改進歐拉法的具體改進歐拉法的具體 形式為形式為2/)()2()2(11cpnpnpncnnnnpyyyyxyhyyy

14、xyhyy 095909. 1)091818. 11 . 1(21)(21091818. 1)1 . 12 . 01 . 1( 1 . 01)2(1 . 0 1 . 1)01( 1 . 01)2(1 . 01100000 yyyyxyyyyxyyycpppcp具體計算過程如下具體計算過程如下上一頁上一頁 下一頁下一頁 返回返回 xn改進的歐拉法改進的歐拉法誤差誤差xn改進的歐拉法改進的歐拉法誤差誤差0100.61.4859560.0027160.21.1840960.0000880.81.6164760.0040240.41.3433600.0017191.01.7378690.005818

15、184096.1)(21180942.1)2(1 .0187250.1)2(1 .02211111cpppcpyyyyxyyyyxyyy依次計算可得依次計算可得 y3，y4，y5，y6，y7，y8，y9，y10其部分結(jié)果見下表其部分結(jié)果見下表 上一頁上一頁 下一頁下一頁 返回返回 例例3 對下面的初值問題對下面的初值問題1)0(10,yxyy解解 （1取步長取步長 h=0.1, 歐拉方法的具體公式為歐拉方法的具體公式為nnnnyyhyy9.0)(1 905. 02/ )(91. 0)(1 . 0 9 . 0)(1 . 01ncpnnpncnnnpyyyyyyyyyyyy（2取步長取步長 h=0

16、.1, 改進的歐拉方法的具體公式為改進的歐拉方法的具體公式為取步長取步長h=0.1，分別用，分別用Euler方法、改進的方法、改進的Euler方法求數(shù)值解。方法求數(shù)值解。上一頁上一頁 下一頁下一頁 返回返回 計算結(jié)果見下表計算結(jié)果見下表Euler方法方法改進的改進的Euler方法方法xnynyn0.10.20.30.40.50.60.70.80.91.00.900 0000.810 0000.729 0000.656 1000.590 4900.531 4410.478 2970.430 4670.387 4210.348 6790.905 0000.819 0250.741 2180.670

17、 8020.607 0760.549 4040.497 2100.449 9750.407 2280.368 541上一頁上一頁 下一頁下一頁 返回返回 .,),(),(11拉拉公公式式高高精精度度比比歐歐值值代代替替平平均均斜斜率率兩兩點點的的斜斜率率的的平平均均相相當(dāng)當(dāng)于于取取點點hkyxyxnnnn第二節(jié)第二節(jié) 龍格龍格 - 庫塔法庫塔法基本思想基本思想考察改進的歐拉法，可以將其改寫為：考察改進的歐拉法，可以將其改寫為：2/)( ),( ),(2111121kkhyyhkyxfkyxfknnnnnn斜率斜率一定取一定取k1 k2 的平均值嗎？的平均值嗎？步長一定是一個步長一定是一個h 嗎

18、？嗎？根根據(jù)據(jù)微微分分中中值值定定理理的的準(zhǔn)準(zhǔn)確確解解為為初初值值問問題題設(shè)設(shè),)(),()(0 yaybxayxfyxy),(),(,()()()()(11 nnnnnnnnnxxyhfxyhyxyxy .,)()(,(1上上的的平平均均斜斜率率在在區(qū)區(qū)間間為為解解曲曲線線 nnnnxxxyyf 只要能對平均斜率提供一種近似算法只要能對平均斜率提供一種近似算法,就能得到一種對應(yīng)的差分格式就能得到一種對應(yīng)的差分格式.,),(),(精精度度較較低低平平均均斜斜率率近近似似代代替替處處的的斜斜率率點點歐歐拉拉公公式式相相當(dāng)當(dāng)于于取取一一個個nnnnyxfyx上一頁上一頁 下一頁下一頁 返回返回 .

19、,1公公式式可可能能構(gòu)構(gòu)造造出出精精度度更更高高的的代代替替平平均均斜斜率率用用他他們們的的加加權(quán)權(quán)平平均均值值上上多多預(yù)預(yù)測測幾幾個個點點的的斜斜率率在在區(qū)區(qū)間間 nnxx例如取例如取 m 個點的斜率構(gòu)造如下形式的公式個點的斜率構(gòu)造如下形式的公式 )(),(),(),(),(2211111,1123213133121221mmnnmmmmnmnmnnnnnnkckckchyyhkbhkbyhaxfkhkbhkbyhaxfkhkbyhaxfkyxfk,的的常常數(shù)數(shù)無無關(guān)關(guān)是是與與其其中中nfcbaiiji該公式稱為該公式稱為m級龍格庫塔級龍格庫塔(Runge-Kutta)公式公式,簡稱簡稱R-

20、K公式公式.求解：只需將公式的局部截斷誤差在求解：只需將公式的局部截斷誤差在xn點進行點進行Taylor展開展開,令其令其前面盡可能多的項為前面盡可能多的項為0, 便可導(dǎo)出便可導(dǎo)出ai，bij，ci所滿足的方程組所滿足的方程組,即即可從中求出這些系數(shù)可從中求出這些系數(shù). 上一頁上一頁 下一頁下一頁 返回返回 以以 m=2 的情形為例說明建立的情形為例說明建立R-K公式的方法公式的方法.公公式式級級建建立立令令KRaba 2,212 )(),(),(22111121kckchyyahkyahxfkyxfknnnnnn.,21待待定定其其中中cca其局部截斷誤差為：其局部截斷誤差為：)()()()

21、(2211111kckchxyxyyxyTnnnn點點進進行行泰泰勒勒展展開開得得在在nx)()(,(1nnnxyxyxfk )()()(,()()(,(22hOxyyfxfahxyxfxyahxyahxfknnnnnn )(6)(2)()()(321nnnnnxyhxyhxyhxyxy),(),(),(),(),(),()(yxfyxfyxfdxdyyxfyxfyxfdxdxyyxyx )()()(2hOxyahxynn 上一頁上一頁 下一頁下一頁 返回返回 因此有：因此有：)()(21)()(132221hOxyhacxyhccTnn .,2項項系系數(shù)數(shù)為為零零與與必必須須的的階階數(shù)數(shù)盡盡

22、可可能能高高要要使使hhT 211221accc即滿足：即滿足： 而對于而對于h3h3，若將，若將k2k2的的TaylorTaylor展開式多取一項展開式多取一項, ,會發(fā)現(xiàn)會發(fā)現(xiàn)h3h3項的系項的系數(shù)不可能為數(shù)不可能為0. 0. 而對于上式有無窮多個解而對于上式有無窮多個解, ,它的每一組解都給出了一個局部截它的每一組解都給出了一個局部截斷誤差為斷誤差為 的二級的二級R-KR-K公式公式, ,即二階即二階R-KR-K公式公式. .)(3hOT 當(dāng)取當(dāng)取1,2121acc時時, ,二階二階R-KR-K公式就是改進的公式就是改進的EulerEuler公式公式 這里有這里有 個未知個未知數(shù)，數(shù)，

23、個方程。個方程。32上一頁上一頁 下一頁下一頁 返回返回 )22( ),( ),( ),( ),(43216134222312221kkkkyyhkyhxfkkyxfkkyxfkyxfkhnnnnhnhnhnhnnn 常用的標(biāo)準(zhǔn)四階常用的標(biāo)準(zhǔn)四階RK公式經(jīng)典公式經(jīng)典R-K方法）方法）最常用的四階標(biāo)準(zhǔn)最常用的四階標(biāo)準(zhǔn)RK公式經(jīng)典公式經(jīng)典R-K方法為：方法為：上一頁上一頁 下一頁下一頁 返回返回 例例 用四階標(biāo)準(zhǔn)用四階標(biāo)準(zhǔn)R-K公式解初值問題公式解初值問題 1)0(10 ,2yxyxyy解：取解：取 h=0.2 h=0.2 ，四階標(biāo)準(zhǔn)，四階標(biāo)準(zhǔn)R-KR-K法的具體格式如下法的具體格式如下: : )

24、22(6 )(2)(),(2)2(2)2()2,2(2)2(2)2()2,2( /2),(432113334222311121kkkkhyyhkyhxhkyhkyhxfkkhyhxkhykhyhxfkkhyhxkhykhyhxfkyxyyxfknnnnnnnnnnnnnnnnnnnnnn上一頁上一頁 下一頁下一頁 返回返回 10y5 , 1 , 0 , 2 . 0nnxn知知 183229. 1)22(6 84324. 033849. 0181728. 1)(2)(90864. 018318. 0091818. 12)2(2)2( 91818. 01 . 12 . 01 . 12)2(2- )

25、2( 1012-4321013003042002031001020001kkkkhyyhkyhxhkykkhyhxkhykkhyhxkhykyxyk上一頁上一頁 下一頁下一頁 返回返回 同理可計算得同理可計算得， , ,5432yyyy具體結(jié)果見下表具體結(jié)果見下表 至少具有四位有效數(shù)字至少具有四位有效數(shù)字. 比較比較: :上節(jié)用改進的上節(jié)用改進的EulerEuler公式計算公式計算, ,取取h=0.1h=0.1，最多具有四位有，最多具有四位有效數(shù)字效數(shù)字 。 上一頁上一頁 下一頁下一頁 返回返回 改進的改進的EulerEuler公式每前進一步只要計算兩次公式每前進一步只要計算兩次f f 值值,

26、 ,而而4 4階階R-KR-K公式每前進一步要計算四次公式每前進一步要計算四次f f 值值, ,但改進的但改進的EulerEuler法的步長比法的步長比4 4階階R-KR-K法的小一半法的小一半, ,兩者計算總量差不多兩者計算總量差不多. . 而4階R-K法的效果要比改進的Euler法好. 由于龍格由于龍格-庫塔法的導(dǎo)出基于泰勒展開，庫塔法的導(dǎo)出基于泰勒展開，故精度主要受解函數(shù)的光滑性影響。對故精度主要受解函數(shù)的光滑性影響。對于光滑性不太好的解，最好采用低階算于光滑性不太好的解，最好采用低階算法而將步長法而將步長h 取小。取小。上一頁上一頁 下一頁下一頁 返回返回 第三節(jié)第三節(jié) 單步法的收斂性

27、與穩(wěn)定性單步法的收斂性與穩(wěn)定性 收斂性收斂性 /* Convergency */定義定義 若某算法對于任意固定的若某算法對于任意固定的 x = xi = x0 + i h，當(dāng)，當(dāng) h0 ( 同時同時 i ) 時有時有 yi y( xi )，則稱該算法是收斂的。，則稱該算法是收斂的。 例：就初值問題例：就初值問題 考察歐拉顯式格式的收斂性?？疾鞖W拉顯式格式的收斂性。0)0(yyyy 解：該問題的精確解為解：該問題的精確解為 xeyxy 0)( 歐拉公式為歐拉公式為iiiiyhyhyy)1 (1 0)1 (yhyii 對任意固定的對任意固定的 x = xi = i h ，有，有iixhhxihyh

28、yy )1()1(/10/0 ehhh /10)1(lim)(0ixxyeyi 上一頁上一頁 下一頁下一頁 返回返回 穩(wěn)定性穩(wěn)定性 /* Stability */例：考察初值問題例：考察初值問題 在區(qū)間在區(qū)間0, 0.5上的解。上的解。分別用歐拉法、隱式歐拉法和改進的歐拉格式計算數(shù)值解。分別用歐拉法、隱式歐拉法和改進的歐拉格式計算數(shù)值解。 1)0()(30)(yxyxy0.00.10.20.30.40.5精確解精確解改進歐拉法改進歐拉法隱式歐拉法隱式歐拉法歐拉法歐拉法 節(jié)點 xixey30 1.00002.0000 4.00008.0000 1.6000101 3.2000101 1.0000

29、2.5000101 6.25001021.56251023.90631039.76561041.00002.50006.25001.56261013.90631019.76561011.00004.97871022.47881031.23411046.14421063.0590107What is wrong ?!上一頁上一頁 下一頁下一頁 返回返回 定義定義 若某算法在計算過程中任一步產(chǎn)生的誤差在以后的計若某算法在計算過程中任一步產(chǎn)生的誤差在以后的計算中都逐步衰減，則稱該算法是絕對穩(wěn)定的算中都逐步衰減，則稱該算法是絕對穩(wěn)定的. 一般分析時為簡單起見，只考慮試驗方程一般分析時為簡單起見，只考慮試驗方程yy 常數(shù)常數(shù)l0，可以是復(fù)數(shù)可以是復(fù)數(shù) 當(dāng)步長取為當(dāng)步長取為 h 時，將某算法應(yīng)用于上式