第4章常微分方程數值解-ppt課件

1、第第4章章 常微分方程數值解常微分方程數值解n4.1 微分方程在化工中的運用微分方程在化工中的運用 n4.2 歐拉歐拉Euler公式公式 n4.3 龍格龍格-庫塔方法庫塔方法 n4.4 常微分方程組的數值解法常微分方程組的數值解法 n4.5 程序例如及運用程序例如及運用 4.1 微分方程在化工中的運用微分方程在化工中的運用 微分方程在化工中運用的簡單而又典型的例子是套管式換熱器的穩(wěn)態(tài)溫度分布。首先作以下假設：1、套管內側為液體，其溫度只隨套管的長度改動而改動，忽略溫度的徑向變化；套管環(huán)隙為蒸汽，其溫度在任何位置均為恒定值，可以為是飽和蒸汽的溫度。2、忽略套管內側流體的縱向熱傳導。3、在整個套管

2、長度方向上，總傳熱系數K不變。4.1 微分方程在化工中的運用微分方程在化工中的運用0)()(222dldtdltuCrtTrdlKtuCrPWp蒸汽入口流體入口，u，t0冷凝液出口流體出口，u，tL圖4-1 套管式換熱器溫度分布表示圖流入的熱量+傳入的熱量-流出的熱量=0)( 2tTrCuKdldtWP4-14-24-34.1 微分方程在化工中的運用微分方程在化工中的運用n另一個在化工中常見的微分方程是物料冷卻過程的數學模型，其模型可用下式表示： )(0TTkdtdT4-3在微分方程中我們稱自變量函數只需一個的微分方程為常微分方程，自變量函數個數為兩個或兩個以上的微分方程為偏微分方程。給定微分

3、方程及其初始條件，稱為初值問題；給定微分方程及其邊境條件，稱為邊值問題。4.1 微分方程在化工中的運用微分方程在化工中的運用n在化工模擬中主要碰到的是常微分方程的初值問題：)( ,)(),()( 0bxayayyxfxy)( ,)(),(0bxayayyxfdxdy4-5或對于大多數常微分方程的初值問題, 只能計算它的數值解。常微分方程初值問題的數值解就是求y(x)在求解區(qū)間a,b上各個分點序列xn，n =1,2,m的數值解yn。在計算中商定y(xn)表示常微分方程準確解的值，yn表示y(xn)的近似值。4. 2 歐拉歐拉Euler公式公式n4.2.1 向前歐拉公式向前歐拉公式n4.2.2 向

4、后歐拉公式向后歐拉公式 n4.2.3 中心歐拉公式中心歐拉公式n4.2.4 梯形公式梯形公式 4.2.1 向前歐拉公式向前歐拉公式),(1nnnnyxhfyy4-6下式為計算近似值的向前歐拉公式：圖4.1 歐拉折線法圖4-2 歐拉折線法幾何表示圖4.2.1 向前歐拉公式向前歐拉公式 實例實例n例例4.1：假定某物體的溫度：假定某物體的溫度w因自熱而產生的熱量可因自熱而產生的熱量可以使物體在每秒鐘內以以使物體在每秒鐘內以4%的速度增長，同時該物體的速度增長，同時該物體由于散熱可使其溫度在每秒種內下降由于散熱可使其溫度在每秒種內下降100k，那么物，那么物體溫度隨時間變化的微分方程：體溫度隨時間變

5、化的微分方程： n t以秒為單位以秒為單位) n 分別以初始溫分別以初始溫x(0)=1500k，y(0)=2500k，z(0)=3500k用歐拉公式預測用歐拉公式預測24秒后的物體溫度趨勢。秒后的物體溫度趨勢。10004. 0wdtdw4.2.1 向前歐拉公式向前歐拉公式 實例實例n解：n w0分別以x0=1500，y0=2500，z0=3500代入。計算結果見表4-1。1,10004. 1)10004. 0(1hwwhwwnnnnn xn yn zn n xn yn zn 1 2 3 4 5 6 7 8 9 10 11 12 1460 1418.4 1375.14 1330.14 1283.

6、35 1234.68 1184.07 1131.43 1076.69 1019.76 960.546 898.968 2500 2500 2500 2500 2500 2500 2500 2500 2500 2500 2500 2500 3540 3581.6 3624.86 3669.86 3716.65 3765.32 3815.93 3868.57 3923.31 3980.24 4039.45 4101.03 13 14 15 16 17 18 19 20 21 22 23 24 834.926 768.324 699.056 627.019 552.1 474.183 393.151

7、 308.877 221.232 130.081 35.2845 -63.3042 2500 2500 2500 2500 2500 2500 2500 2500 2500 2500 2500 2500 4165.07 4231.68 4300.94 4372.98 4447.9 4525.82 4606.85 4691.12 4778.77 4869.92 4964.72 5063.3 05101520250500100015002000250030003500400045005000wt圖4-3 三種初始值的溫度變化曲線 表4-14.2.1 向前歐拉公式向前歐拉公式 實例實例n從表4-1可

8、以看到當自熱引起物體溫度升高的速度小于散熱引起溫度下降的速度，物體的溫度隨時間而逐漸減少：當自熱引起物體溫度升高的速度與散熱引起溫度下降的速度平衡時，物體的溫度堅持不變；當自熱引起物體溫度升高的速度大于散熱引起溫度下降的速度，物體的溫度隨時間而增長。在圖4-3中L1，L2，L3分別表示初始值3500，2500和1500的三條溫度變化趨勢曲線。 h充分小時，以上迭代收斂。記 ，那么 h充分小時，可保證 ，其中L為李普希茲條件。 4.2.2 向后歐拉公式向后歐拉公式),(111nnnnyxhfyy)0(1ny給定精度直到)(1)1(1)(11)1(1)0(1, 2 , 1 , 0,),(),(kn

9、knknnnknnnnnyykyxhfyyyxhfyy向后歐拉公式：4-6式(4-7)是yn+1的非線性方程，即隱式歐拉公式，用迭代法求得yn+1。初始值 由向前歐拉公式提供。),()(1yxhfyynn最簡單的迭代公式為：),()( 1yxhfyny1),(1hLyxhfny4.2.3 中心歐拉公式中心歐拉公式ny(x)的在x=x1處的中心差商式：n 又 ，可得到y(tǒng)(x2)的近似值y2計算公式：n 類似地，可得到計算y(xn+1)近似值yn+1的計算公式：n n 公式(4-8)稱為中心格式。按公式(4-8)，需求知道yn-1, yn的值才干求得yn+1的值。因此，要先用其它公式計算出y1，再

10、用中心格式算出y2, y3,。y1可用向前歐拉公式計算，為提高精度，也可用向后歐拉公式計算。hxyxyxy2)()()( 021)(,()( 111xyxfxy),(21102yxhfyy),(211nnnnyxhfyy(4-8)4.2.4 梯形公式梯形公式n梯形公式：),(),(2111nnnnnnyxfyxfhyy),(),(2),(1111nnnnnnnnnnyxfyxfhyyyxhfyy),(,(),(211nnnnnnnnyxhfyxfyxfhyy(4-9)梯形公式也是隱式格式。n用顯式的歐拉公式和隱式的梯形公式給出的一次預估-校正公式：(4-10)上式也稱為改良的歐拉公式，它可合并

11、成：4.2.4 梯形公式梯形公式 實實例例n例4.2：請用預估-校正公式改良的歐拉公式解右面初值問題：0.40.0 , 1)0(2xyydxdyn解：1 . 0, 10hy)(22)(12)1(12)0(1knnnknnnnyyhyyhyyy 用下面的迭代公式，對每個點迭代4次，k=1,2,3,4。 4.2.4 梯形公式梯形公式 實實例例n nx ny )(nxy )(nnxyy 1 2 3 4 0.1 0.2 0.3 0.4 1.1118 1.2520 1.4311 1.6763 1.1111 1.2500 1.4326 1.6667 0.0007 0.0020 0.0095 0.0004

12、xy11該方程的準確解是 計算結果如表4-2所示。表4-2 計算結果4.3 龍格龍格-庫塔方法庫塔方法 本課程主要研討其實踐運用,，故直接給出各類龍格-庫塔公式。),(),()(2121211hkyhxfkyxfkkkhyynnnnnn)2,2(),(12121khyhxfkyxfkhkyynnnnnnn1、二階龍格-庫塔其中c1=0,c2=1,a=1/2,b=1/2。其中c1=1/2, c2=1/2, a=1, b=1或(4-11)(4-12)4.3 龍格龍格-庫塔方法庫塔方法n2、三階龍格、三階龍格-庫塔公式庫塔公式 )2,()2,2(),()4(6) 1 (2131213211hkhky

13、hxfkkhyhxfkyxfkkkkhyynnnnnnnn2312131132,3231,31,34 )2(hkyhxfkhkyhxfkyxfkkkhyynnnnnnnn23121321143,4321,21,4329 ) 3(hkyhxfkhkyhxfkyxfkkkkhyynnnnnnnn(4-15)(4-14)(4-13)4.3 龍格龍格-庫塔方法庫塔方法n3、四階龍格、四階龍格庫塔公式庫塔公式 342312143211,21,2121,21,226 ) 1 (hkyhxfkhkyhxfkhkyhxfkyxfkkkkkhyynnnnnnnnnn321421312143211,31,3231

14、,31,338 )2(hkhkhkyhxfkhkhkyhxfkhkyhxfkyxfkkkkkhyynnnnnnnnnn(4-17)(4-16)4.3 龍格龍格-庫塔方法庫塔方法 實例實例n例例4.3：用四階龍格：用四階龍格庫塔公式庫塔公式(4-16)求解下面初求解下面初值問題值問題8 . 01 . 0 , 1) 0 (cos2xyxydxdy2 . 0cos2 . 01 . 0cos1 . 01 . 0cos1 . 0cos2262 . 02342232122143211nnnnnnnnnnxkykxkykxkykxykkkkkyyn nx ny nxy nnxyy 1 2 3 4 0.2 1

15、.24789 1.24792 0.00003 0.4 1.63762 1.63778 0.00016 0.6 2.29618 2.29696 0.00078 0.8 3.53389 3.53802 0.00413 n解：取步長h=0.2，計算公式為：表4-3 計算結果4.3 龍格龍格-庫塔方法庫塔方法 步長的步長的選擇選擇n下面以四階龍格-庫塔方法為例，闡明如何自動選擇步長，使計算結果滿足給定精度的要求。n 設從節(jié)點xn出發(fā)，先以h為步長，利用四階龍格-庫塔公式方法經過一步計算得y(xn+1)的近似值，記為 ，由于公式的部分截斷誤差是y(h5)，故有 n 當h不大時，c可近似地看作常數。然后將

16、步長h對折，即取h/2為步長，從出發(fā)經過兩步計算求y(xn+1)的近似值，記為 ，每一步計算的部分截斷誤差為c(h/2)5，于是就有 )(1hny5)(11)(chyxyhnn5)2/(11)2/()(hcyxyhnn)2/(1hny(4-18)4.3 龍格龍格-庫塔方法庫塔方法 步長的步長的選擇選擇n把它與4-18式相比，可得：161)()()(11)2/(11hnnhnnyxyyxy)()()(11151)2/(11hnnhnnyxyyxy經整理可得：4.3 龍格龍格-庫塔方法庫塔方法 步長的步長的選擇選擇n這闡明以 作為y(xn+1) 的近似值，其誤差可用先后兩次計算結果之差來表示，因此

17、，只需調查n 能否成立。假設成立，那么可將 作為y(xn+1)的近似值；假設不成立，那么將步長再次對折進展計算，直到不等式成立為止，并取最后的 作為計算結果。以上方法就是計算過程中自動選擇步長的方法，也稱為變步長方法。)2/(1hny)2/(11)(hnnyxy)2/(1hny)2/(1hny4.3 龍格龍格-庫塔方法庫塔方法 步長的步長的選擇選擇用一個例子闡明：由VB選用規(guī)范四階龍格-庫塔方法計算得 因此合理選擇步長既能保證精度又能減少計算量。 1)0( yyxdxdy4.4 常微分方程組的數值解法常微分方程組的數值解法n4.4.1 一階常微分方程組的數值解法一階常微分方程組的數值解法n4.

18、4.2 高階常微分方程數值方法高階常微分方程數值方法 4.4.1 一階常微分方程組的數值解法一階常微分方程組的數值解法n將由m個一階方程組成的常微分方程初值問題：btaayayayyyytfdtdyyyytfdtdyyyytfdtdymmmmmmm)()()(),(),(),(22112121222111)(),(aYytFdtdYmmmmmmyytfyytfyytfytFtytytytY211121121,),(),(),(),(,)()()()(向量方式：其中：(4-19)(4-20)4.4.1 一階常微分方程組的數值解法一階常微分方程組的數值解法n下面以兩個方程組為例，給出相應的計算公式

19、。n常微分方程組： n歐拉公式： n預估校正公式： btazazyayzytfdtdzzytfdtdy )()(),(),(00),(),(11nnnnnnnnnnzythgzzzythfyy),(),(),(),(2),(),(1111111111nnnnnnnnnnnnnnnnnnnnnnnnnnzytgzytfzytgzytfhzyzyzytgzytfhzyzy(4-21)(4-22)4.4.1 一階常微分方程組的數值解法一階常微分方程組的數值解法n四階龍格庫塔公式：)2(3)1(3)2(3)1(3)2(4)1(44)2(2)1(2)2(2)1(2)2(3)1(33)2(1)1(1)2(

20、1)1(1)2(2)1(22)2(1)1(11)2(4)1(4)2(3)1(3)2(2)1(2)2(1)1(11143211,2,2,22,2,22,2,22,2,2),(),(226226hkzhkyhtghkzhkyhtfkkKkhzkhyhtgkhzkhyhtfkkKkhzkhyhtgkhzkhyhtfkkKzytgzytfkkKkkkkkkkkhzyzyKKKKhYYnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn(4-23)4.4.1 一階常微分方程組的數值解法一階常微分方程組的數值解法 實例實例n例4.4：兩種微生物，其數量分別是u=u(t)，v=v(t)，t的單位為分

21、，其中一種微生物以吃另一種微生為生，兩種微生物的增長函數如以下常微分方程組所示，預測3分鐘后這一對微生物的數量。2.1)0(6.1)0(001.015106.045.020109.0vuuvvvdtdvuvuudtdu4.4.1 一階常微分方程組的數值解法一階常微分方程組的數值解法 實例實例n解：記uvvvvuguvuuvuf001. 015106. 0),(45. 020109. 0),(11111111,2,nnnnnnnnnnnnnnnnnnnnvugvufvugvufhvuvuvugvufhvuvut(m) )(tu )(tv 1 1.6 1.2 2 0.28766 1.33103 3 0.0291308 1.47776 4 0.00021727 1.63947 用歐拉預估校正公式(4-22)表4-4 計算結果4.4.2 高階常微分方程數值方法高階常微分方程數值方法n以三階常微分方程為例闡明高階常微分方程的數值計算步驟。 )()()()(,)()2() 1 ()0(3btaayayayyyytfdxtyd )()()()(3221tydttdytydttdy)2(3)1(2)0(132133221)()()()(),()