CH01常微方程數(shù)值解1.1-1.7解析

1、主要內(nèi)容第第3 3章章 拋物型方程的有限差分方法拋物型方程的有限差分方法第第1 1章章 常微分方程初值問(wèn)題數(shù)值解法常微分方程初值問(wèn)題數(shù)值解法第第2 2章章 橢圓型方程的有限差分方法橢圓型方程的有限差分方法第第4 4章章 雙曲型方程的有限差分方法雙曲型方程的有限差分方法第第5 5章章 非線(xiàn)性雙曲型守恒律方程的差分方法非線(xiàn)性雙曲型守恒律方程的差分方法第第6 6章章 有限元方法簡(jiǎn)介有限元方法簡(jiǎn)介第第1 1章章 常微分方程初值問(wèn)題的常微分方程初值問(wèn)題的數(shù)值解法數(shù)值解法 1.11.1 引言引言 1.2 1.2 歐拉法歐拉法( (Euler方法方法) 1.3 1.3 梯形法、隱式格式的迭代計(jì)算梯形法、隱式

2、格式的迭代計(jì)算 1.4 1.4 一般單步法、一般單步法、Runge-Kutta格式格式 1.5 1.5 線(xiàn)性多步法線(xiàn)性多步法 1.6 1.6 誤差的事后估計(jì)法、步長(zhǎng)的自動(dòng)選擇誤差的事后估計(jì)法、步長(zhǎng)的自動(dòng)選擇 1.7 1.7 高階常微分方程高階常微分方程( (組組) )的數(shù)值方法的數(shù)值方法 1.1 1.1 引言引言 目標(biāo)在于給出解在一些離散點(diǎn)上的近似值。目標(biāo)在于給出解在一些離散點(diǎn)上的近似值。 本章研究常微分方程初值問(wèn)題的主要數(shù)值本章研究常微分方程初值問(wèn)題的主要數(shù)值解法，包括基本方法和基本理論問(wèn)題。解法，包括基本方法和基本理論問(wèn)題。1.2 1.2 歐拉法歐拉法( (Euler方法方法) ) 1.2

3、.1 1.2.1 歐拉方法歐拉方法 )(|),(000 xyyXxxyxfdxdyxx)2 . 1 () 1 . 1 (考慮常微分方程初值問(wèn)題注：在后面的討論中，我們總認(rèn)為這個(gè)初值問(wèn)題的解存在、唯一且連續(xù)依賴(lài)于初值條件，即初值問(wèn)題(1.1)，(1.2)是適定的。),(yxfdxdy)(,()()(1nnnnxyxhfxyxy 將解的存在區(qū)間 N 等分，得到N 個(gè)小區(qū)間。 任取一個(gè)小區(qū)間 ,1nnxx，由原方程 1),()()(1nnxxnndxyxfxyxy得),(yxf在區(qū)間 ,1nnxx上，用 在點(diǎn) nx上的值來(lái) 代替),(yxf，得到NxXxxhnn01為步長(zhǎng)。其中圖圖1.1 0)(,(

4、xyxfnx1nxx.)( ! 21)( )()(21nnnnxyhxhyxyxy)(,()( nnnxyxfxy)(,()()(1nnnnxyxhfxyxy在上式中分別用ny和來(lái)代替1ny)(nxy)(1nxy和并由n 的任意性，得到)(,.2 , 1 , 0),(001xyynyxhfyynnnn(1.3)這就是歐拉公式。歐拉公式亦可由Taylor 展式得到在上式中分別用ny和來(lái)代替1ny)(nxy)(1nxy和),(1nnnnyxhfyy則得),()(11nnnnnnyxfxxyy一般而言，并不要求步長(zhǎng)相等，則有(1.4)0圖圖1.2nxy0y3210 xxxx xyy x幾何意義 例例

5、 1.11.1 以的數(shù)值解，并與精確解為步長(zhǎng)，用歐拉法求初值問(wèn)題1 . 0h1)0(yyxedxdyx比較。xexxy)2(21)(2 (1)、計(jì)算格式本身不能準(zhǔn)確描述原來(lái)的方程 誤差的產(chǎn)生：誤差的產(chǎn)生：(2)、計(jì)算機(jī)本身引入的誤差 計(jì)算機(jī)輸出的是歐拉方程的近似解 ,而不是精確解 。因此ny)(nxy)()()(nnnnnnyyxyyxyy(1.5) 可見(jiàn)，為了使計(jì)算得到的解 是 的好的近似，我們要求：ny)(nxy (1)歐拉方法的精確解 是微分方程精確解 的很好近似。ny)(nxy(2) 是 好的近似。nyny問(wèn)題(1)稱(chēng)為格式的收斂性問(wèn)題。問(wèn)題(2)稱(chēng)為格式的穩(wěn)定性問(wèn)題。1.2.2 收斂

6、性研究收斂性研究0hxnhxxn0)(xyyn 所謂收斂性問(wèn)題,就是研究時(shí),要求,。111)(nnnyxy整體截?cái)嗾`差整體截?cái)嗾`差局部截?cái)嗾`差局部截?cái)嗾`差nnnnRyxye1*11)()(,()(1*nnnnxyxhfxyy這里1),()(,()()(11nnxxnnnnnnyxhfydxxyxfxyyxy(1.6)即ThTh 定理定理 1.11.1221|MhRn(1.12)| )(|max0 xyMXxx 其中h。為步長(zhǎng)，XxCxyy,)(02nR的局部截?cái)嗾`差，則歐拉方法滿(mǎn)足 假定Back其中 R 為局部截?cái)嗾`差的上界。 定理定理 1.21.2 設(shè) f(x,y) 關(guān)于 y 滿(mǎn)足Lipsc

7、hitz 條件，L為相應(yīng)的Lipschitz常數(shù)，則歐拉方法的整體截?cái)酀M(mǎn)足n) 1(|)(0)(00LxXLxXnehLRe(1.13)誤差Th1.4由定理1.1，1.2，可得 定理定理 1.31.3 設(shè) f(x,y) 關(guān)于 y 滿(mǎn)足Lipschitz 條件，L為相應(yīng)的Lipschitz常數(shù)，(1.14)XxCxyy,)(02)(00 xyy )(nxy且當(dāng)h0，并有估計(jì)式，則歐拉方法的解一致收斂到初值問(wèn)題(1.1)，(1.2)的解) 1(2|)(0)(00 xXLxXLneLMheny)(00 xyy 00) 1(2|)(0 xXLneMLh)(|hOn如果，即，由此有(1.15)即 歐拉方

8、法的整體截?cái)嗾`差與歐拉方法的整體截?cái)嗾`差與h h 同階，由同階，由 的的表達(dá)式可知，表達(dá)式可知， ，這說(shuō)明局部截?cái)嗾`差比整，這說(shuō)明局部截?cái)嗾`差比整體截?cái)嗾`差高一階。體截?cái)嗾`差高一階。nR)(2hORn我們稱(chēng)歐拉方法為一階格式。我們稱(chēng)歐拉方法為一階格式。1.2.3 1.2.3 穩(wěn)定性研究穩(wěn)定性研究 前已指出歐拉方法的穩(wěn)定性問(wèn)題是決定歐拉法在利用計(jì)算機(jī)能否得到精確解的關(guān)鍵問(wèn)題，只有穩(wěn)定的算法才可能是有用的算法。定義1.100,zy如果存在正常數(shù) c 及 ，使對(duì)任意初始值0h01),(zzxhfzznnnn與01),(yyxhfyynnnn，由nnzy ,計(jì)算所得之解滿(mǎn)足估計(jì)式0000;0|xXnh

9、hhzyczynn則稱(chēng)歐拉方法穩(wěn)定歐拉方法穩(wěn)定。注意注意: 這里 分別是以 為初值得到的精確值，毫無(wú)舍入誤差，因此這里穩(wěn)定性定義式對(duì)初值的穩(wěn)定性，即研究初值誤差在計(jì)算過(guò)程中的傳遞問(wèn)題。nnzy ,00, zy 定理定理1.4 1.4 在定理在定理1.21.2的條件下，歐拉方法是穩(wěn)的條件下，歐拉方法是穩(wěn)定的。定的。Th1.2 由定理1.2，我們看到如初始誤差 ，則整體截?cái)嗾`差的階完全由局部截?cái)嗾`差的階決定，事實(shí)上，若局部截?cái)嗾`差階為 ，則整體截?cái)嗾`差階為 。因此為了提高數(shù)值算法的精度，往往從提高局部截?cái)嗾`差的階入手，這也時(shí)構(gòu)造高精度差分方程數(shù)值方法的主要依據(jù)。00)(1phO)(phO1.3 梯

10、形法、隱式格式的迭代計(jì)算梯形法、隱式格式的迭代計(jì)算dxxyxfdxyxfnnnnxxxxnn11)(,(),()(,(nnxyxhf 在歐拉方法的推導(dǎo)過(guò)程，用矩形公式近似計(jì)算積分1),()()(1nnxxnndxyxfxyxy)(,(2)(1nnnnxyxfxx)(.(11nnxyxf若用梯形公式近似計(jì)算積分，則0圖圖1.31nxnxxyxf,yxf,)(,(2)(),(11nnnnxxxyxfxxdxyxfnn)(.(11nnxyxf因此有)(,()(,()(21)()(1111nnnnxnnnxyxfxyxfxxxyxy),(),(21111nnnnnnyxfyxfhyy(1.16)這是一

11、個(gè)隱式格式。梯形公式局部截?cái)嗾`差：首先給出后退歐拉公式的局部截?cái)嗾`差：)(,()()(111nnnnxyxhfxyxy后退歐拉公式假定 和解 充分光滑，則yxf, xy*111)(nnnyxye)(2132hOh 而歐拉公式的局部截?cái)嗾`差為*111)(nnnyxye)(2132hOh 梯形公式的局部截?cái)嗾`差為)(,(21)()()(11*11nnnnnnnxyxfhxyxyyxye),(*11nnyxf)(3hO 類(lèi)似于歐拉法，對(duì)梯形法和后退歐拉法也可平行地建立它們的整體截?cái)嗾`差的階分別為 和 ，以及格式的收斂性和穩(wěn)定性等定理。)(2hO)(hO),(),(2111nnnnnnyxfyxfhy

12、y(1.16)如何求解 ，采用迭代法，其格式如下： 1ny初始猜測(cè))0(1)(11)1(1),(),(2nnnpnnnpnyyxfyxfhyy(1.18)迭代法的收斂性前已指出，梯形法是一個(gè)隱式格式12Lh(1.19)迭代法收斂的充分條件。 當(dāng) ，有下面的預(yù)報(bào)-校正格式：0p也稱(chēng)為改進(jìn)歐拉公式改進(jìn)歐拉公式。校正格式預(yù)報(bào)格式),(),(2),()0(111)0(1nnnnnnnnnnyxfyxfhyyyxhfyy(1.20)當(dāng)然也可迭代多次： 校正格式預(yù)報(bào)格式),(),(2),()(11)1(1)0(1pnnnnnpnnnnnyxfyxfhyyyxhfyy(1.21) 當(dāng)步長(zhǎng) 取得適當(dāng)小，用預(yù)報(bào)

13、格式(歐拉法)已能算出比較好的近似值，故迭代收斂很快，通常只需迭代二三次就可滿(mǎn)足精度要求，如果迭代多次仍不收斂，說(shuō)明步長(zhǎng)過(guò)大，必須減少步長(zhǎng) ，再進(jìn)行計(jì)算。 hh 梯形法較之歐拉法提高了精度，但增加了迭代次數(shù)，因此增加了計(jì)算工作量。 例例 1.21.2 試用預(yù)報(bào)校正格式(1.20)解初值問(wèn)題1|1 , 010 xyxxyy(1.21)1 . 0h取 。從計(jì)算結(jié)果可以看出，歐拉法精度較低，預(yù)報(bào)校正格式精度有所改善，大約精確到3為有效數(shù)字。 1.4 1.4 一般單步法、一般單步法、 格式格式KuttaRunge 前面，我們研究了歐拉法和梯形法，它們有一個(gè)共同的特點(diǎn)，即在格式中只包括 的值，或者說(shuō)由

14、，僅使用 的值計(jì)算出 的值，這種格式稱(chēng)為單步格式，下面研究一般單步法。11,nnnnyxyx1nnxxny1ny1.4.1 1.4.1 一種構(gòu)造單步法的方法一種構(gòu)造單步法的方法泰勒級(jí)數(shù)法泰勒級(jí)數(shù)法 設(shè)初值問(wèn)題00|),(yyyxfdxdyxx的解 階可微，將 在 點(diǎn)展開(kāi)為泰勒級(jí)數(shù)，有 ) 1()(qxy xy0 x)()(!.)(! 2)()()(10)(02000 qqqhOxyqhxyhxyhxyhxy(1.22)由方程可得),()(000yxfxyyyxfyxfyxfxyyx ),(),(),()(因此),(),(),()(0000000yxfyxfyxfxyyx )()(2)()(2y

15、xyyyxyxxfffffffffxydxdxy ,.,xxyxfff 式中， 都是相對(duì)于變量的偏導(dǎo)數(shù)。于是式 (1.22)可寫(xiě)成 )(),(,()()(10000qhOhxyxhxyhxy其中)(,(11110000| )(,(!1(),(,(xyxqjjjjxyxfdxdhjhxyx舍去 ，可得)(1qhO),(0001hyxhyy),(1112hyxhyy),(1hyxhyynnnn稱(chēng) 為一般單步法，顯然局部截?cái)嗾`差 ),(1hyxhyynnnn*)(nnnyxye),(,()()(111hxyxhxyxynnnnqjnjjjnnxydxdhjxyxy111)(!1)()()(1qhO所

16、以局部截?cái)嗾`差為 ，在式(1.23)中令 ，即得歐拉法。)(1qhO1q 為任意關(guān)于 的函數(shù)，其對(duì)于微分方程 ),(,(hxyx),(,(hxyx),(yxfdxdy的解 滿(mǎn)足 )(xy)(),(,()()(1qhOhxyxhxyhxy(1.24)且 為使上式成立的最大整數(shù)，則稱(chēng)q,.2 , 1 , 0),(1nhyxhyynnnn(1.25)為 階單步法，歐拉法為一階單步法，泰勒級(jí)數(shù)法式(1.23)為 階單步法。qq 定義定義 1.2 1.2 給出單步法),(1hyxhyynnnn1.4.2 1.4.2 一般單步法基本理論一般單步法基本理論 定義定義 1.31.3 如果， 則稱(chēng)單步法 為與初

17、值問(wèn)題(1.3)相容的。 ),()0 ,(yxfyx),(1hyxhyynnnn)(),()()(1qhOhyxhxyhxy從而)0 ,()()(lim0yxhxyhxyh由單步法的定義得:因此有如下定義: 定理定理 1.51.5 如果 對(duì)于 ， 以及所有實(shí)數(shù) 滿(mǎn)足 條件，則單步法(1.23)穩(wěn)定。滿(mǎn)足：的解nnzy ,),(hyxXxx000hh yLipschitz|00zyczynn01),(zhzxhzznnnn01),(yhyxhyynnnn欲使定理成立要證明存在常數(shù) , 對(duì)0, 00hc 定理定理 1.61.6 如果 對(duì)于 ， 以及所有實(shí)數(shù) 關(guān)于 滿(mǎn)足 條件，則 收斂的充要條件是格

18、式相容，即滿(mǎn)足 。),(hyxXxx000hh yhyx,Lipschitz),(1hyxhyynnnn),()0 ,(yxfyx由相容性可以得到格式的收斂性: 定理定理 1.71.7 在定理1.5的條件下，如果局部截?cái)嗾`差 為 ，則單步法 的整體截?cái)嗾`差 滿(mǎn)足nnR)(1qhO),(1hyxhyynnnnnnnyxy)() 1(|)(0)(001xXLqxXLneLcheee(1.26)特別若 ,則 ，整體截?cái)嗾`差比局部截?cái)嗾`差低一階。00)(qnhO 關(guān)于單步法的整體截?cái)嗾`差 ，有:1.4.3 1.4.3 格式格式KuttaRunge 從前面討論可見(jiàn)，構(gòu)造高階單步法的關(guān)鍵在于構(gòu)造 ，使),

19、(hyx)(),(,()()(11qnnnnhOhxyxhxyxy中的局部截?cái)嗾`差階盡可能高。前面我們利用泰勒級(jí)數(shù)法構(gòu)造了一個(gè)歐拉法，這時(shí) ，局部截?cái)嗾`差與 同階，這是一個(gè)一階格式。為了要求 ，利用泰勒級(jí)數(shù)法得到一個(gè)二階格 1),(),(,()()(21qhOhxyxhfxyxynnnn)(2hO2q),(),(),(2),(21nnnnynnxnnnnyxfyxfyxfhyxhfyy(1.27)這時(shí)我們有 )(,()()(1nnnnxyxhfxyxy)()(,()(,()(,(232hOxyxfxyxfxyxfhnnnnynnx 格式(1.27)計(jì)算過(guò)程中要求函數(shù) 的二個(gè)偏導(dǎo)數(shù) 在 處的值，

20、比較麻煩，可以預(yù)計(jì)，利用泰勒級(jí)數(shù)法推導(dǎo)出的高階格式需要求更多的偏導(dǎo)數(shù)值，計(jì)算繁復(fù)。那么是否可以避免計(jì)算偏導(dǎo)數(shù)，而得到高階單步格式的 呢？分析梯形法的預(yù)報(bào)校正格式(1.20) ),(yxfyxff ,)(,(nnxyx),(hyx校正格式預(yù)報(bào)公式),(),(2),()0(111)0(1nnnnnnnnnnyxfyxfhyyyxhfyy 二級(jí)二階 方法。一般而言，二級(jí)二階 格式可以寫(xiě)成KuttaRunge KuttaRunge ),(),()(12122122111KhbyhaxfKyxfKKcKchyynnnnnn(1.28) 適當(dāng)選擇參數(shù) ，使局部截?cái)嗾`差 21221,bacc)(,()()(

21、111nnnnnxyxfchxyxyR)(,()(,(2122nxnnxyxfhbxyhaxfc)(3hO由)(,(2)(,()(21nnxnnnnxyxfhxyxhfxyR)() )(,()(,(3hOxyxfxyxfnnnny)(,()(,()(21nnnnnxyxfhcxyxfhcxy)() )(,()(,()(,(3212hOxyxfhbxyxfhaxyxfnnnnynnx 因此要求滿(mǎn)足212112212221cbcacc(1.29) (1) 取 ，則 ,即得二級(jí)二階 法211c1,212122bacKuttaRunge ),(),()(2121211hKyhxfKyxfKKKhyyn

22、nnnnn(1.30)01c21, 12122bac)21,21(),(12121hKyhxfKyxfKhKyynnnnnn(1.31)這是一個(gè)含有四個(gè)參數(shù)、三個(gè)方程的方程組,因此由一個(gè)自由參數(shù)，解答不唯一。改進(jìn)歐拉公式變形歐拉公式 (2)令 ， ，由此得算式為 (3) 取 ， ，則有411c32,432122bac)32,32(),()3(4121211hKyhxfKyxfKKKhyynnnnnn(1.32) 根據(jù)同樣的思想,可以構(gòu)造更高階精度的Runge-Kutta 方法。),(),(),()(232131331212213322111hKbhKbyhaxfKhKbyhaxfKyxfKKc

23、KcKchyynnnnnnnn(1.33) 三級(jí)三階 Runge-Kutta法一般可以寫(xiě)成：616131)(31)(3121)(2112132332232323132212323133212223322232313212332232bbcbacbbcbcbbacbacacacbbcbcacacccc(1.34)(41hORn故要求，必須有(1) 令 ，則 ，6131 cc1,21,64322aac1, 2,21313221bbb故有 三級(jí)三階算法Kutta)2,()2,2(),()4(62131213211hKhKyhxfKhKyhxfKyxfKKKKhyynnnnnnnn(1.36)(2)令

24、 ，解得 ，32,3132aa43, 0,313221ccb0,32,4131321bbc，故有 三級(jí)三階 算法HeumKR)32,32()3,3(),()3(4213121311hKhKyhxfKhKyhxfKyxfKKKhyynnnnnnnn(1.37)式(1.36)，(1.37)是三級(jí)三階 格式。KuttaRunge 格式的局部截?cái)嗾`差為 。)(4hO可以設(shè)計(jì)四級(jí)四階 格式KuttaRunge ),(),(),(),()(3432421414423213133121221443322111hKbhKbhKbyhaxfKhKbhKbyhaxfKhKbyhaxfKyxfKKcKcKcKchy

25、ynnnnnnnnnn(1.38) 為了達(dá)到四級(jí)四階格式，可得13個(gè)參數(shù)滿(mǎn)足11個(gè)方程241121)(81)(61)(4131211433224432342224322234433422432323433422432233443333222442332224433224321434241432313212bbacbabacbacababacbaacbabacbacacacacacacacacacacccccbbbabbaba(1) 經(jīng)典四級(jí)四階 格式取定 ，則得：KuttaRunge 2132 aa),()2,2()2,2(),()122(6342312143211hKyhxfKhKyhxfK

26、hKyhxfKyxfKKKKKhyynnnnnnnnnn(1.39)這是最為著名的經(jīng)典四級(jí)四階 格式。KuttaRunge 格式的局部截?cái)嗾`差為 。)(5hO1.5 1.5 線(xiàn)性多步法線(xiàn)性多步法 前面利用歐拉法或梯形法求未知函數(shù)在 的近似值 ，基本思想是在積分表達(dá)式1nxx1nydxxyxfxyxynnxxnn1)(,()()(1中被積函數(shù) ，用水平直線(xiàn)或連接 兩點(diǎn)的直線(xiàn) 代替。然而為了近似 中的曲線(xiàn) 也可用多點(diǎn)插值曲線(xiàn)，如我們利用 插值可得經(jīng)過(guò) 的曲線(xiàn) 為1,),(nnxxxyxf)(,(),(nnxyxfyxf)(,(,(nnnxyxfx)(,(,(111nnnxyxfx)(,(),(1n

27、nnxyxfhxxyxf)(,(11nnnxyxfhxx1,nnxx)(,(xyxfLagrange)(,(,(222nnnxyxfx)(,(,( ,111nnnxyxfx)(,( ,nnxyx)(,(2xyxL)(,(2)()(,(22212nnnnxyxfhxxxxxyxL)(,()(1122nnnnxyxfhxxxx)(,(2)(212nnnnxyxfhxxxx(1.50)用它近似 中的 ，則得到積分近似值1,nnxx)(,(xyxfdxxyxLxyxynnxxnn1)(,()()(21)(,(125)(22nnnxyxfhxy)(,(1223)(,(121611nnnnxyxfxyxf

28、像歐拉格式一樣，我們可得近似求解格式如下：),(125221nnnnyxfhyy),(1223),(121611nnnnyxfyxf(1.51)與歐拉格式 不同之處在于增加了包括 , 的兩項(xiàng),即由計(jì)算 ，與歐拉格式僅由前面一點(diǎn) 計(jì)算 的這種單步法不同，格式(1.51)稱(chēng)為多步法。 現(xiàn)在，設(shè)已給出常微分方程初值問(wèn)題),(1nnnnyxhfyy),(22nnyxf),(11nnyxf),(),(),(1122nnnnnnyxyxyx1ny),(nnyx1ny00|),(yyyxfdxdyxx的解 在 處的近似值 ，或者說(shuō)給出表頭)(xynxxx,.,10nyyy,.,10 xynxxxx210ny

29、yyy210研究如何由表頭給出 處的近似值 。 根據(jù)Xhnxxn) 1(011ny1)()()(1nnxxnndxxyxyxy(1.52)和 ，利用表頭的值和 插值法求出 的近似表達(dá)式，再利用(1.52)就得到 的近似值。 具體操作如下： 用 表示用 處 的值構(gòu)造出的 的 插值多項(xiàng)式，用 表示相應(yīng)的插值余項(xiàng)，即從而),()(yxfxyLagrange)(xy)(1nxy)(,xLknknnnxxx,.,1)(),.,(),(1knnnxyxyxy)(xyLagrangeknr,)()()(,xrxLxyknkn11)()()()(,1nnnnxxknxxknnndxxrdxxLxyxy舍去余項(xiàng) ，并用 代替 ，則可得 的近似值的表達(dá)式1)(,nnxxknkndxxrRjy)(jxy)(1nxydxLyynnxxknnn1*,1(1.53) 為 中 的值用 代替， 為局部截?cái)嗾`差。上述格式中，被插值點(diǎn) 不包括在插值基點(diǎn)所決定的最大區(qū)間 內(nèi)，故稱(chēng)為外插公式。由常微分方程， 根據(jù) ， ， ， 得插值公式為 ，由于插值基點(diǎn)為等距，令 ，利用牛頓后插公式，有)(*,xLkn)(,xLkn)(jxyjyknR,)(1nnxxxxn