常微分方程的差分方法

1、安工大2004 4.1數(shù)值分析簡明教程第三章 常微分方程的差分方法1 歐拉方法2 改進的歐拉方法3 龍格-庫塔方法4 亞當姆斯方法5 收斂性與穩(wěn)定性6 方程組與高階方程的情形7 邊值問題安工大2004 4.2數(shù)值分析簡明教程引言 科學技術當中常常需要求解常微分方程的定解問題。這類問題的最簡單的形式，是本章著重要考察的一階方程的初值問題： 本章中我們假定右函數(shù)適當光滑以保證初值問題解的存在唯一。雖然求解常微分方程有各種各樣的解析方法，但求解從實際問題中歸結出來的微分方程要靠數(shù)值解法。 差分法是一類重要的數(shù)值方法，這類方法是要尋求 離散節(jié)點上的近似解 ，相鄰節(jié)點間距 稱為步長步長。初值問題的各種差

2、分方法都采用“步進式”，即求解過程順著節(jié)點排列的次序一步一步地向前推進。描述這類算法，只要給出從已知信息 計算 的遞推公式，這類計算格式統(tǒng)稱為差分格差分格 式式。00,yf x yy xy12nxxx12,ny yy 1nnhxx12,nnnyyyny| ),(),(:|*yyLyxfyxfLipshitz條件且滿足安工大2004 4.3數(shù)值分析簡明教程歐拉格式 微分方程的本質(zhì)特征是方程中含有導數(shù)項，這也是它難于求解的癥結所在。數(shù)值解法的第一步就是設法消除其導數(shù)值，這項手續(xù)稱為離散化。離散化。實現(xiàn)離散化的基本途徑就用差商代替導數(shù)。譬如，若在點 列出方程，并用差商 代替 ，結果有設用 的近似值

3、代入上式右端，記所得結果為 ，這樣導出的計算公式就是眾所周知的歐拉（歐拉（Euler）格式）格式，若初值 是已知的，則依據(jù)上式即可逐步算出數(shù)值解 。nx1nny xy xhnyx1,nnnny xy xhf xy xny xny1ny1,0,1,2,nnnnyyhf xyn12,y y 0y安工大2004 4.4數(shù)值分析簡明教程安工大2004 4.5數(shù)值分析簡明教程安工大2004 4.6數(shù)值分析簡明教程萊昂哈德歐拉（Leonhard Euler，1707-1783）18世紀最優(yōu)秀的數(shù)學家，也是歷史上最偉大的數(shù)學家之一，被稱為“分析的化身”。萊昂哈德歐拉（Leonhard Euler，1707-

4、1783），1707年出生在瑞士的巴塞爾城，小時候他就特別喜歡數(shù)學，不滿10歲就開始自學代數(shù)學。這本書連他的幾位老師都沒讀過，可小歐拉卻讀得津津有味，遇到不懂的地方，就用筆作個記號，事后再向別人請教。13歲就進巴塞爾大學讀書，這在當時是個奇跡，曾轟動了數(shù)學界。小歐拉是這所大學，也是整個瑞士大學校園里年齡最小的學生。在大學里得到當時最有名的數(shù)學家微積分權威約翰伯努利（Johann Bernoulli，1667-1748年）的精心指導，并逐漸與其建立了深厚的友誼。約翰伯努利后來曾這樣稱贊青出于藍而勝于藍的學生：“我介紹高等分析時，他還是個孩子，而你將他帶大成人?！眱赡旰蟮南奶?，歐拉獲得巴塞爾大學的

5、學士學位，次年，歐拉又獲得巴塞爾大學的哲學碩士學位。1725年，歐拉開始了他的數(shù)學生涯。1783年9月18日，在不久前才剛計算完氣球上升定律的歐拉，在興奮中突然停止了呼吸，煙斗從手中落下，口里喃喃地說：“我要死了”，歐拉終于“停止了生命和計算”，享年76歲。歐拉生活、工作過的三個國家：瑞士、俄國、德國，都把歐拉作為自己的數(shù)學家，為有他而感到驕傲 . 安工大2004 4.7數(shù)值分析簡明教程他從19歲開始發(fā)表論文，直到76歲，半個多世紀寫下了浩如煙海的書籍和論文可以說歐拉是科學史上最多產(chǎn)的一位杰出的數(shù)學家，據(jù)統(tǒng)計他那不倦的一生，共寫下了886本書籍和論文（七十余卷，牛頓全集八卷，高斯全集十二卷），

6、其中分析、代數(shù)、數(shù)論占40%，幾何占18%，物理和力學占28%，天文學占11%，彈道學、航海學、建筑學等占3%，彼得堡科學院為了整理他的著作，足足忙碌了四十七年。到今幾乎每一個數(shù)學領域都可以看到歐拉的名字，從初等幾何的歐拉線，多面體的歐拉定理，立體解析幾何的歐拉變換公式，四次方程的歐拉解法到數(shù)論中的歐拉函數(shù)，微分方程的歐拉方程，級數(shù)論的歐拉常數(shù)，變分學的歐拉方程，復變函數(shù)的歐拉公式等等，數(shù)也數(shù)不清他對數(shù)學分析的貢獻更獨具匠心，無窮小分析引論一書便是他劃時代的代表作，當時數(shù)學家們稱他為分析學的化身安工大2004 4.8數(shù)值分析簡明教程19世紀偉大數(shù)學家高斯（Gauss，1777-1855年）曾說

7、：研究歐拉的著作永遠是了解數(shù)學的最好方法歐拉的一生，是為數(shù)學發(fā)展而奮斗的一生，他那杰出的智慧，頑強的毅力，孜孜不倦的奮斗精神和高尚的科學道德，永遠是值得我們學習的歐拉在數(shù)學、物理、天文、建筑以至音樂、哲學方面都取得了輝煌的成就。在數(shù)學的各個領域，常常見到以歐來命名的公式、定理、和重要常數(shù)。課本上常見的如（1736年），i（1777年），e（1748年），sin和cos（1748年），tg（1753年），x（1755年），（1755年），f(x)（1734年）等，都是他創(chuàng)立并推廣的。歌德巴赫猜想也是在他與歌德巴赫的通信中提出來的。歐拉還首先完成了月球繞地球運動的精確理論，創(chuàng)立了分析力學、剛體力學

8、等力學學科，深化了望遠鏡、顯微鏡的設計計算理論。安工大2004 4.9數(shù)值分析簡明教程歐拉不但重視教育，而且重視人才。當時法國的拉格朗日只有19歲，而歐拉已48歲。拉格朗日與歐拉通信討論“等周問題”，歐拉也在研究這個問題。后來拉格朗日獲得成果，歐拉就壓下自己的論文，讓拉格朗日首先發(fā)表，使他一舉成名 。安工大2004 4.10數(shù)值分析簡明教程七橋問題七橋問題 七橋問題七橋問題Seven Bridges Problem18世紀著名古典數(shù)學問題之一。在哥尼斯堡的一個公園里，有七座橋?qū)⑵绽赘駹柡又袃蓚€島及島與河岸連接起來(如圖)。問是否可能從這四塊陸地中任一塊出發(fā)，恰好通過每座橋一次，再回到起點？歐拉

9、于1736年研究并解決了此問題，他把問題歸結為如下右圖的“一筆畫”問題，證明上述走法是不可能的。給出了連通網(wǎng)絡可一筆畫的充要條件是它們是連通的，且奇頂點它們是連通的，且奇頂點(通過此點弧通過此點弧的條數(shù)是奇數(shù)的條數(shù)是奇數(shù))的個數(shù)為的個數(shù)為0或或2。安工大2004 4.11數(shù)值分析簡明教程1733年，年僅26歲的歐拉擔任了彼得堡科學院數(shù)學教授1735年，歐拉解決了一個天文學的難題（計算彗星軌道），這個問題經(jīng)幾個著名數(shù)學家?guī)讉€月的努力才得到解決，而歐拉卻用自己發(fā)明的方法，三天便完成了然而過度的工作使他得了眼病，并且不幸右眼失明了，這時他才28歲1741年歐拉應普魯士彼德烈大帝的邀請，到柏林擔任科學

10、院物理數(shù)學所所長，直到1766年，后來在沙皇喀德林二世的誠懇敦聘下重回彼得堡，不料沒有多久，左眼視力衰退，最后完全失明不幸的事情接踵而來，1771年彼得堡的大火災殃及歐拉住宅，帶病而失明的64歲的歐拉被圍困在大火中，雖然他被別人從火海中救了出來，但他的書房和大量研究成果全部化為灰燼了沉重的打擊，仍然沒有使歐拉倒下，他發(fā)誓要把損失奪回來歐拉完全失明以后，雖然生活在黑暗中，但仍然以驚人的毅力與黑暗搏斗，憑著記憶和心算進行研究，直到逝世，竟達17年之久 安工大2004 4.12數(shù)值分析簡明教程安工大2004 4.13數(shù)值分析簡明教程歐拉格式的精度 為簡化分析，人們常在 為準確即 的前提下估計誤差 ，

11、這種誤差稱為局部截斷誤差局部截斷誤差。 如果一種數(shù)值方法的局部截斷誤差為 ， 則稱它的的精度是 階的，或稱之為 階方法。 歐拉格式僅為一階方法。nynnyy x11nny xyp1pO hp)( 21)()( )(),( 21)( )()(,2)(2111211yhyxyEulerxhyxyyyhxhyxyxyxxhTaylorxxynnnnnnnnnnn公式為展開式階點的在安工大2004 4.14數(shù)值分析簡明教程隱式歐拉格式 設改用后差商 替代方程 中的導數(shù)項 ，再離散化，即可導出下列格式該格式右端含有未知的 ，它實際上是個關于 的函數(shù)方程。故稱該格式隱式歐拉格式隱式歐拉格式。隱式歐拉格式也

12、是一階方法,與歐拉格式相當。11nny xy xh111,nnnyxf xy x1nyx111,nnnnyyhf xy1ny1ny)(21)(),( )()(21)( )()(),(21)( )()(2)(211112112111yhyxyxhyxyyyhxhyxyxyyhxhyxyxyTaylorxxynnnnnnnnnnnn展開式階點的在安工大2004 4.15數(shù)值分析簡明教程兩步歐拉格式1 設改用中心差商 替代方程 中的導數(shù)項 ，再離散化，即可導出下列格式 無論是顯式歐拉格式還是隱式歐拉格式，它們都是單步法，其特點是計算時只用到前一步的信息 ，而該格式卻調(diào)用了前面兩步的信息 ，兩步歐拉格

13、式兩步歐拉格式因此而得名。 兩步歐拉格式具有更高的精度，它是二階方法。1112nny xy xh,nnnyxf xy x112,nnnnyyhf xyny1,nnyy安工大2004 4.16數(shù)值分析簡明教程兩步歐拉格式2.)()()( 2)()()( 2)()(),()(21)( )()(),()(21)( )()(,2)(3111131132132111公式為二階方法所以兩步公式為兩步展開式階點的在EulerhOyxyEulerxhyxyyhOxhyxyxyhOxyhxhyxyxyhOxyhxhyxyxyxxxxhTaylorxxynnnnnnnnnnnnnnnnnnnnn安工大2004 4

14、.17數(shù)值分析簡明教程微分方程的數(shù)值解haxbxxaxbxaxyxybxxaxnahbxxxxahnhbyxxahhbbaxhyyahxhbbaxhyyhbbaxhyynnhxxyynbaxhyyEulerbxaxxyybaxynnnnnnnnnnnnnnnnnn21)21(21)(21) 121 ()()(3)(2)()(0,0, 0, 0)0(0),(,21)(:0)0(,:12121212122311120010012公式由易知其準確解設微分方程安工大2004 4.18數(shù)值分析簡明教程梯形公式1 設將方程 的兩端從 到 求積分，即得，顯然，只要能近似的算出其中的積分項，我們就可以得到計算

15、 的差分格式。 若我們用梯形法計算積分項：再離散化，即可得如下計算公式與梯形求積公式相呼應的這一差分格式稱為梯形格式梯形格式。,yf x ynx1nx 11,nnxnnxy xy xfx y xdx 111,2nnxnnnnxhf x y xf xy xf xy x111,2nnnnnnhyyf xyf xy1ny x安工大2004 4.19數(shù)值分析簡明教程梯形公式2),(),(2),(),(,)(,()(,(2)()(),(,()(,(2)(,(,),(,()()()(,()(,(),(,()()()(,()(,()()()( )(,(),(,()( 11111111111111111111

16、11nnnnnnnnnnnnnnnnnnnnxxnnnnnnxxnnnnnnxxnnxxxxyxfyxfhyyxyxyyyxyxfxyxfhxyxyxyxfxyxfhdxxyxfEulerxyxhfxyxyxyxhfdxxyxfEulerxyxhfxyxyxyxhfdxxyxfxyxydxxydxxyxfxyxfxynnnnnnnnnn即得梯形公式替換分別用式若使用較準確的梯形公公式即得隱式公式即得則設安工大2004 4.20數(shù)值分析簡明教程改進的歐拉格式 先用歐拉法求得一個初步的近似值，記為 ，稱之為預報值，然后用它替代梯形法右端的 再直接計算 ，得到校正值 ，這樣建立的預報校正系統(tǒng)稱為改進

17、的歐拉格式改進的歐拉格式： 預報 校正 它有下列平均化形式： 實踐表明，改進的歐拉格式明顯改善了精度。1ny1ny1ny111,2nnnnnnhyyfxyfxy11,/2pnnncnnpnpcyyhf xyyyhf xyyxy1nf),(1nnnnyxhfyy安工大2004 4.21數(shù)值分析簡明教程安工大2004 4.22數(shù)值分析簡明教程龍格-庫塔法的設計思想 考察差商 ，根據(jù)微分中值定理，存在點 ,利用所給方程 得我們稱 為區(qū)間 上的平均斜率，這樣只要對平均斜率 提供一種算法，相應地我們便導出一種計算格式。 龍格庫塔（RungeKutta）方法設計思想就是設法在 內(nèi)多預報幾個點的斜率值，然后

18、把它們加權平均作為平均斜率，以期望構造出更高精度的計算格式。1nny xy xh yf 11,nnnny xy xhfyxx， ,Kfy1,nnxxK1,nnxx安工大2004 4.23數(shù)值分析簡明教程二階龍格-庫塔方法1 隨意考察區(qū)間 內(nèi)一點 ， 用兩個點 的斜率 的加權平均代替平均斜率 ,于是我們就得到如下計算格式: (類似于改進的歐拉格式)其中有兩個待定參數(shù) ， 適當選取它們的值，就可使上述格式有較高的精度。若 ，該格式是二階的 ，故統(tǒng)稱滿足這一條件的一族格式為二階龍格庫塔格式二階龍格庫塔格式。特別地，當 時，上述格式即為改進的歐拉格式改進的歐拉格式，如果取 ，則上述格式稱為變形的歐拉格

19、式，亦稱為中點格式中點格式。1,nnxx,01npnxxphp,nnpxx12,K KK1121211,nnnnn pnyyhKKKf xyKf xyphK, p12p11,2p1,12p安工大2004 4.24數(shù)值分析簡明教程二階龍格-庫塔方法221),()( )( )()1()()( )( )(),()( )(),(),(),(),(),()( ),(),(),(2),(21),(),(),(),(22),(322112221112122phOxyphxhyxyKKhyyhOxphyxyhOdxyxdfphxyhOyxfphKyxphfyxfphKyphxfphKyxfKxyyxfKbaf

20、sbarsfbafrbasfbarfbafsbrafTayloryxfnnnnnnnnnnnnynnxnnnnnpnnnnyyxyxxyx又展開式階元的安工大2004 4.25數(shù)值分析簡明教程二階龍格-庫塔方法3)2,(),()1(,1,21),(),()(2)1(,21, 121)()( 21)( )()(2)(12121211121211211321KhyxfKyxfKhKyyKuttaRungeKKhyyphKyxfKyxfKKKhyyEulerKKhyypphOxyhxhyxyxyTaylorxynpnnnnnnnnpnnnnnnnnnnn公式為二階時當公式為改進的時當展開式階的安工大

21、2004 4.26數(shù)值分析簡明教程作業(yè)P.124第9題安工大2004 4.27數(shù)值分析簡明教程三階龍格-庫塔方法 為了進一步提高精度，我們可以考慮用三個點 的斜率值 加權平均得出平均斜率 的近似值，其中，于是就可以構造所謂的三階龍格庫塔格式三階龍格庫塔格式 ，下列庫塔格式庫塔格式是其中的一種： ,nnpn qxxx123,K KKK,01,1npnn qnxxphqxxqh pq1123121/2131124/6,/2,2nnnnnnnnyyh KKKKf xyKf xyhKKf xyhKK安工大2004 4.28數(shù)值分析簡明教程四階龍格-庫塔方法 繼續(xù)上述過程，我們可以導出四階龍格庫塔格式四

22、階龍格庫塔格式，下列經(jīng)典格式是其中的一種： 值得注意的是，龍格庫塔法的推導基于泰勒展開法，因而它要求解具有較好的光滑性。如果解的光滑性差，則該方法得到的解反而不好。11234121/2131/2241322/6,/2,/2,nnnnnnnnnnyyh KKKKKf xyKf xyhKKf xyhKKf xyhK安工大2004 4.29數(shù)值分析簡明教程變步長的龍格-庫塔方法 同積分的數(shù)值計算一樣，微分方程的數(shù)值解法也需要選擇步長。同樣，我們可以采取步長加倍或折半的辦法選擇步長，即通過檢查步長折半前后的兩種計算結果的偏差：來判斷選取的步長是否合適，具體可以分為兩種情況來處理：對于給定精度 ，若 ，

23、則反復將步長折半進行計算直到 為止，取步長折半后的“新值”作為結果；相反的，反復將步長加倍直到 ，取步長加倍前的“老值”作為結果。 這種通過步長加倍或折半的手續(xù)處理步長的方法稱為變步長變步長方方法法。 211hhnnyy安工大2004 4.30數(shù)值分析簡明教程亞當姆斯格式1 亞當姆斯（亞當姆斯（Adams)方法方法的設計思想是充分利用計算 之前已得到一系列節(jié)點 上的斜率值來減少計算量。譬如，我們可以用 兩點的斜率的加權平均作為區(qū)間 上的平均斜率，于是可設計出如下二階亞當姆斯格式: Adams方法亦成為線性多步方法亦成為線性多步法法類似的，可導出如下三階和四階亞當姆斯格式：1ny1,nnxx1,

24、nnx x1,nnxx113/2nnnnyyhyy11223165/12nnnnnyyhyyy11235559379/24nnnnnnyyhyyyy安工大2004 4.31數(shù)值分析簡明教程亞當姆斯格式2)3(2,21)()( )( )()()( )( )( 1)( )()( 21)( )()(2)()( ),()( ),()1(1132121132111121211nnnnnnnnnnnnnnnnnnnnnnnnnyyhyyAdamshOxyhxhyxyyhOxhyxyxyTaylorxxyhOxyhxhyxyxyTaylorxyxyyxfKxyyxfKKKhyy格式二階展開式階的在展開式階的

25、假設安工大2004 4.32數(shù)值分析簡明教程亞當姆斯格式3格式即法均可以看作為線性多步梯形公式公式兩步公式隱式公式AdamsyyhyyxfyxfhyyhyyyxhfyyEulerhyyyxhfyyEulerhyyyxhfyyEulernnnnnnnnnnnnnnnnnnnnnnnnnnn,)(2),(),(2,2),(2,),(,),(,111111111111安工大2004 4.33數(shù)值分析簡明教程隱式亞當姆斯格式1 同樣，我們也可導出如下隱式的二階、三階和四階亞當姆斯格式：11/2nnnnyyh yy11158/12nnnnnyyhyyy11129195/24nnnnnnyyhyyyy安工

26、大2004 4.34數(shù)值分析簡明教程隱式亞當姆斯格式2)(2,21)()()( )()()(21)( )()(2)()( ),()( ),()1(,113213211112121111111nnnnnnnnnnnnnnnnnnnnnnnnnnnnnyyhyyAdamshOxyhxhyxyyhOxyhxhyxyxyTaylorxyxyyxfKxyyxfKKKhyyxxxxxxxxxAdams格式是梯形公式二階隱式展開式階的最簡單的形式為即采用內(nèi)插過程的數(shù)值的數(shù)值來預報也可以由是外推過程的數(shù)值的數(shù)值來預報格式是由二階安工大2004 4.35數(shù)值分析簡明教程作業(yè))51623(12:3211nnnnn

27、yyyhyyAdams格式階試證明下式為安工大2004 4.36數(shù)值分析簡明教程證明3階Adams格式)()(6)(2)( )()()2*12521*1216)()2*1251216)()12512161223)( )()51623(12)()(2)(2)( ),()(2)()( ,4)3(324)3(322113)3(223)3(21hOxyhxyhxhyxyhOxyhxyhxhyxyyyyhyyhOxyhxhyxyyhOxyhxhyxyyTaylornnnnnnnnnnnnnnnnnnnnn展開式由安工大2004 4.37數(shù)值分析簡明教程亞當姆斯預報-校正系統(tǒng) 仿照改進的歐拉格式的構造方法

28、，將顯式和隱式兩種亞當姆斯格式相匹配，可構成下列亞當姆斯預報校正系統(tǒng)亞當姆斯預報校正系統(tǒng)： 預報 校正11231115559379/24,nnnnnnnnnyyhyyyyyfxy11231119195/24,nnnnnnnnnyyhyyyyyf xy安工大2004 4.38數(shù)值分析簡明教程改進的亞當姆斯預報-校正系統(tǒng) 我們可以方便地估計出亞當姆斯預報校正系統(tǒng)的截斷誤差，從而依據(jù)這種估計將該系統(tǒng) 就可改進為如下精度更高的計算方案： 預報 改進 校正 改進1123111115559379/24251270,nnnnnnnnnnnnnpyhyyyympcpmfxm111211111119195/24

29、19270,nnnnnnnnnnnnncyhmyyyyccpyf xy安工大2004 4.39數(shù)值分析簡明教程記住的公式),()2,2()2,2(),()22(63142312143211hKyxfKKhyhxfKKhyhxfKyxfKKKKKhyynnnnnnnnnn),(),(2),(1111nnnnnnnnnnyxfyxfhyyyxhfyy校正值預報值),(1nnnnyxhfyy安工大2004 4.40數(shù)值分析簡明教程收斂性問題1 在用差分格式求解微分方程時我們要考慮差分格式的收斂性。我們稱差分格式是收斂的，如果對任意固定的 ，數(shù)值解 當 （同時 ）時趨于準確解 。 以下我們研究歐拉方法

30、的收斂性。我們記 ，記 為 關于 的李普希茲常數(shù)，經(jīng)反復遞推，可得其中 為常數(shù)。 若初始 是準確的，即 ，則當 時，有 。 這說明歐拉方法是收斂收斂的。0nxxnhny0h n ny xnnney xyLfy01TLTLnCee eehL,C T0y00e 0h 0ne 安工大2004 4.41數(shù)值分析簡明教程收斂性問題2nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnehLChyxyeyxyeehLyxyhLeyyyyLyxfyxfLipshitzyxhfyyxyxhfxyyyyyxyyxyChyxyCxhyxyyyhyxyyhxhyxyxyxyyy)1 (|

31、)(|,)(|,)1 (|)(| ),(),(:|),(),(,()(|)(|)(|)(| ,),( )(),( 21|)(|),( 21)( )()(,.)(,211111*11111111211121121其中條件又存在常數(shù)其中即整體截斷誤差下面討論整體收斂性時的局部截斷誤差準確即假設前面討論局部收斂性安工大2004 4.42數(shù)值分析簡明教程收斂性問題3.),0( , 0, 0,.) 1()1 (,1,! 2)1 ( 1)1(1)1 ()1 ()1 ()1 (1 ,)1 ()1 (1 )1 (,)1 (0000020002012022202121格式收斂即即準確又為常數(shù)設Eulerheey

32、eeeLCheeehLTnhxxehLexxeeeehLhLhLChehLhLhLCheehLhLCheehLCheehLChenTLTLnTLnhLnnhLxnnnnnnn安工大2004 4.43數(shù)值分析簡明教程收斂性問題4., 1)0(,0,0,0,22,)221 ()221 ()221 (,0)22(, 1)0(,22)(2)(2, 1)0(,0,0,122)22(221111的精確解此即為設任給梯形格式考慮yxyyeyhxhxhhehhhhhhyhxnxnhxhhyyyhhyyyhyyyhyyyxyyxnhhhxhhhhhxnnnnnnnnnnnn安工大2004 4.44數(shù)值分析簡明教程穩(wěn)定