數(shù)值分析71常微分方程初值問題數(shù)值解引言

1、第第7章章 常微分方程初值問題數(shù)值解法常微分方程初值問題數(shù)值解法 7.1 引言引言 7.2 歐拉方法與改進(jìn)的歐拉方法歐拉方法與改進(jìn)的歐拉方法 7.3 龍格龍格- -庫塔方法庫塔方法7.1 引引 言言本章將考察本章將考察一階方程的初值問題一階方程的初值問題 )2 . 1(.)()1 . 1(),(00yxyyxfy 我們知道，只要我們知道，只要 f(x, y) 適當(dāng)光滑適當(dāng)光滑譬如關(guān)于譬如關(guān)于 y 滿滿足足利普希茨利普希茨(Lipschitz)條件條件理論上就可以保證初值問題的解理論上就可以保證初值問題的解 yf(x) 存在并且唯一存在并且唯一.),(),(yyLyxfyxf 雖然求解常微分方程

2、有各種各樣的雖然求解常微分方程有各種各樣的解析方法解析方法，但，但解析方法只能用來求解一些特殊類型的方程，實(shí)際問解析方法只能用來求解一些特殊類型的方程，實(shí)際問題中歸結(jié)出來的微分方程主要靠題中歸結(jié)出來的微分方程主要靠數(shù)值解法數(shù)值解法. 所謂所謂數(shù)值解法數(shù)值解法, 就是尋求解就是尋求解y(x)在一系列離散節(jié)點(diǎn)在一系列離散節(jié)點(diǎn) 121nnxxxx上的近似值上的近似值 y1, y2,yn, yn+1,. 相鄰兩個(gè)節(jié)點(diǎn)的間距相鄰兩個(gè)節(jié)點(diǎn)的間距hn=xn+1- -xn 稱為稱為步長步長. 假定假定 hi=h (i=1,2,)為為定數(shù)定數(shù), 這時(shí)這時(shí)節(jié)點(diǎn)為節(jié)點(diǎn)為 xi=x0+ih (i=0,1,2,) (等

3、距節(jié)點(diǎn)等距節(jié)點(diǎn)). 初值問題的初值問題的數(shù)值解法數(shù)值解法有個(gè)有個(gè)基本特點(diǎn)基本特點(diǎn): 都采取都采取“步步進(jìn)式進(jìn)式”，即求解過程順著節(jié)點(diǎn)排列的次序一步一步地，即求解過程順著節(jié)點(diǎn)排列的次序一步一步地向前推進(jìn)向前推進(jìn). 首先，要對微分方程離散化，建立求解數(shù)值解的首先，要對微分方程離散化，建立求解數(shù)值解的遞推公式遞推公式. 一類是計(jì)算一類是計(jì)算 yn+1 時(shí)只用到前一點(diǎn)的值時(shí)只用到前一點(diǎn)的值 yn，稱，稱為為單步法單步法. 另一類是用到另一類是用到 yn+1前面前面 k 點(diǎn)的值點(diǎn)的值 yn, yn-1, yn-k+1，稱為，稱為 k步法步法. 其次，要研究公式的其次，要研究公式的局部截?cái)嗾`差局部截?cái)嗾`差

4、和和階階，數(shù)值解，數(shù)值解 yn 與精確解與精確解 y(xn) 的的誤差估計(jì)誤差估計(jì)及及收斂性收斂性，還有遞推公式的還有遞推公式的計(jì)算穩(wěn)定性計(jì)算穩(wěn)定性等問題等問題.7.2 歐拉方法與改進(jìn)的歐拉方法歐拉方法與改進(jìn)的歐拉方法7.2.1 歐拉法與后退歐拉法歐拉法與后退歐拉法 在在xoy平面上，微分方程平面上，微分方程(1.1)式式的解的解y=f(x)稱作稱作它的它的積分曲線積分曲線，積分曲線積分曲線上一點(diǎn)上一點(diǎn)(x, y)的切線斜率等的切線斜率等于函數(shù)于函數(shù)f(x, y)的值的值. 如果按如果按 f(x, y) 在在xoy平面上建立一平面上建立一個(gè)個(gè)方向場方向場，那么，那么，積分曲線積分曲線上每一點(diǎn)的

5、切線方向均上每一點(diǎn)的切線方向均與方向場在該點(diǎn)的方向相一致與方向場在該點(diǎn)的方向相一致.基于上述幾何解釋，我們從初始點(diǎn)基于上述幾何解釋，我們從初始點(diǎn)P0(x0, y0)出發(fā)出發(fā),先依方向場在該點(diǎn)的方向推進(jìn)到先依方向場在該點(diǎn)的方向推進(jìn)到x=x1上一點(diǎn)上一點(diǎn)P1，然后，然后再從再從P1點(diǎn)依方向場在該點(diǎn)的方向推進(jìn)到點(diǎn)依方向場在該點(diǎn)的方向推進(jìn)到 x=x2 上一點(diǎn)上一點(diǎn)P2 , 循環(huán)前進(jìn)做出一條循環(huán)前進(jìn)做出一條折線折線P0 P1 P2. 一般地，設(shè)已做出該折線的頂點(diǎn)一般地，設(shè)已做出該折線的頂點(diǎn)Pn, ,過過Pn(xn, yn)依方向場的方向再推進(jìn)到依方向場的方向再推進(jìn)到Pn+1(xn+1, yn+1)，顯然

6、兩個(gè)，顯然兩個(gè)頂點(diǎn)頂點(diǎn)Pn, ,Pn+1的坐標(biāo)有關(guān)系的坐標(biāo)有關(guān)系斜斜率率),(111nnnnnnnnyxfhyyxxyy ) 1 .2 (),(1nnnnyxhfyy 這就是著名的這就是著名的( (顯式顯式) )歐拉歐拉( (Euler) )公式公式. . 若初值若初值y0已已知，則依公式知，則依公式(2.1)可逐次逐步算出各點(diǎn)數(shù)值解可逐次逐步算出各點(diǎn)數(shù)值解. .即即),(0001yxhfyy ),(1112yxhfyy 例例1 用歐拉公式求解初值問題用歐拉公式求解初值問題)2 . 2(. 1)0(),10(2 yxyxyy 解解 取步長取步長h=0.1，歐拉公式的具體形式為，歐拉公式的具體形

7、式為)2(1nnnnnyxyhyy 其中其中xn=nh=0.1n (n=0,1,10), 已知已知y0 =1, 由此式可得由此式可得191818. 1)1 . 12 . 01 . 1 ( 1 . 01 . 1)2(1 . 11 . 01)2(1111200001 yxyhyyyxyhyy依次計(jì)算下去，依次計(jì)算下去，部分計(jì)算結(jié)果部分計(jì)算結(jié)果見下表見下表. xy21 與準(zhǔn)確解與準(zhǔn)確解 相比，可看出歐拉公式的計(jì)算結(jié)相比，可看出歐拉公式的計(jì)算結(jié)果精度很差果精度很差. xn 歐拉公式數(shù)值解歐拉公式數(shù)值解yn準(zhǔn)確解準(zhǔn)確解y(xn) 誤差誤差 0.2 0.4 0.6 0.8 1.0 1.191818 1.3

8、58213 1.508966 1.649783 1.784770 1.183216 1.341641 1.483240 1.612452 1.732051 0.008602 0.016572 0.025726 0.037331 0.052719 歐拉公式具有明顯的幾何意義歐拉公式具有明顯的幾何意義, , 就是就是用折線近似用折線近似代替方程的解曲線代替方程的解曲線，因而常稱公式，因而常稱公式(2.1)為為歐拉折線法歐拉折線法. .( )yy xxynx1nxnp1np1np x 還可以通過幾何直觀來考察歐拉方法的精度還可以通過幾何直觀來考察歐拉方法的精度. .假假設(shè)設(shè)yn=y(xn), ,即頂

9、點(diǎn)即頂點(diǎn)Pn落在積分曲線落在積分曲線y=y(x)上，那么，上，那么，按歐拉方法做出的折線按歐拉方法做出的折線PnPn+1便是便是y=y(x)過點(diǎn)過點(diǎn)Pn的切線的切線. .從圖形上看從圖形上看, ,這這樣定出的頂點(diǎn)樣定出的頂點(diǎn)Pn+1偏離偏離了原來的積分曲線，可了原來的積分曲線，可見歐拉方法是見歐拉方法是相當(dāng)粗糙相當(dāng)粗糙的的. . 為了分析計(jì)算公式的精度，通?？捎脼榱朔治鲇?jì)算公式的精度，通?？捎锰├照归_泰勒展開將將y(xn+1)在在xn處展開，則有處展開，則有).,()(2),()()(2)()()()(1221 nnnnnnnnnnnnxxyhyxhfxyyhxyhxyhxyxy 在在yn=y

10、(xn)的前提下，的前提下，f(xn,yn )=f(xn,y(xn)=y ( (xn n) ). .于是于是可得歐拉法可得歐拉法(2.1)的的公式誤差公式誤差為為)3 . 2(),(2)(2)(2211nnnnxyhyhyxy 稱為此方法的稱為此方法的局部截?cái)嗾`差局部截?cái)嗾`差. . 如果對方程如果對方程(1.1)從從xn到到xn+1積分，得積分，得) 4 . 2(.)(,()()(11 nnxxnndttytfxyxy右端積分用右端積分用左矩形公式左矩形公式 hf(xn,y(xn) 近似，再以近似，再以 yn代代替替 y(xn)，yn+1代替代替 y(xn+1)也得到歐拉公式也得到歐拉公式(2

11、.1)，局部，局部截?cái)嗾`差也是截?cái)嗾`差也是(2.3).稱為稱為( (隱式隱式) )后退的歐拉公式后退的歐拉公式. . 如果右端積分用如果右端積分用右矩形公式右矩形公式 hf(xn+1,y(xn+1)近似，近似，則得到另一個(gè)公式則得到另一個(gè)公式) 5 . 2 (),(111 nnnnyxhfyy 后退的歐拉公式與歐拉公式有著后退的歐拉公式與歐拉公式有著本質(zhì)的區(qū)別本質(zhì)的區(qū)別, 后后者是關(guān)于者是關(guān)于yn+1的一個(gè)直接計(jì)算公式，這類公式稱作是的一個(gè)直接計(jì)算公式，這類公式稱作是顯式的顯式的；前者公式的；前者公式的右端含有未知的右端含有未知的yn+1，它實(shí)際上，它實(shí)際上是關(guān)于是關(guān)于yn+1的一個(gè)函數(shù)方程的

12、一個(gè)函數(shù)方程, ,這類方程稱作是這類方程稱作是隱式的隱式的. . 顯式顯式與與隱式隱式兩類方法各有特點(diǎn)，考慮到數(shù)值穩(wěn)兩類方法各有特點(diǎn)，考慮到數(shù)值穩(wěn)定性等其他因素，人們有時(shí)需要選用定性等其他因素，人們有時(shí)需要選用隱式隱式方法，但方法，但使用使用顯式顯式算法遠(yuǎn)比算法遠(yuǎn)比隱式隱式方便方便. .注：注：隱式方程通常用隱式方程通常用迭代法迭代法求解，而迭代過程的實(shí)求解，而迭代過程的實(shí)質(zhì)是質(zhì)是逐步逐步顯式化顯式化. . 設(shè)用歐拉公式設(shè)用歐拉公式),() 0(1nnnnyxhfyy 給出迭代初值給出迭代初值 ，用它代入，用它代入(2.5)式的式的右端，使之轉(zhuǎn)右端，使之轉(zhuǎn)化為顯式，直接計(jì)算得化為顯式，直接計(jì)算

13、得) 0(1 ny),() 0(11) 1 (1 nnnnyxhfyy然后再用然后再用 代入代入(2.5)式，又有式，又有) 1 (1 ny).,() 1 (11) 2(1 nnnnyxhfyy如此反復(fù)進(jìn)行，得如此反復(fù)進(jìn)行，得) 6 . 2 ()., 1 , 0(),()(11) 1(1 kyxhfyyknnnkn如果如果 f(x, y)對對y滿足滿足Lipschitz條件條件(1.3). 由由(2.6)減減(2.5)得得),(),(11)(111) 1(1 nnknnnknyxfyxfhyy.1)(1 nknyyhL由此可知，只要由此可知，只要hL1，迭代法，迭代法(2.6)就收斂到解就收斂

14、到解. .7.2.2 梯形方法梯形方法 為得到比為得到比歐拉法精度高的計(jì)算公式，在等歐拉法精度高的計(jì)算公式，在等式式(2.4) 右端積分用梯形求積公式近似右端積分用梯形求積公式近似, 并用并用yn代替代替y(xn), yn+1代替代替y(xn+1)，則得，則得 ) 7 . 2(,),(),(2111 nnnnnnyxfyxfhyy稱為稱為矩形方法矩形方法. 矩形方法矩形方法是是隱式單步法隱式單步法，用迭代法求解，同后，用迭代法求解，同后退的歐拉方法一樣，仍用退的歐拉方法一樣，仍用歐拉法歐拉法提供迭代初值，則提供迭代初值，則矩形迭代公式矩形迭代公式為為 ) 8 . 2()., 1 , 0(),(

15、),(2);,()(11) 1()0(11 kyxfyxfhyyyxhfyyknnnnnknnnnn為了分析迭代過程的收斂性為了分析迭代過程的收斂性, 將將(2.7)與與(2.8)相減相減, 得得),(),(2)(1111) 1(11knnnnknnyxfyxfhyy 于是有于是有,2)(11) 1(11knnknnyyhLyy 使得使得, 12 hL則當(dāng)則當(dāng)k時(shí)有時(shí)有 , 這說明迭代過程這說明迭代過程(2.8)是收斂是收斂的的.1)(1 nknyy7.2.3 改進(jìn)的歐拉公式改進(jìn)的歐拉公式我們看到，梯形方法雖然提高了精度，但其算法我們看到，梯形方法雖然提高了精度，但其算法復(fù)雜，在應(yīng)用迭代公式復(fù)

16、雜，在應(yīng)用迭代公式(2.8)進(jìn)行實(shí)際計(jì)算時(shí)，每迭代進(jìn)行實(shí)際計(jì)算時(shí)，每迭代一次，都要重新計(jì)算函數(shù)一次，都要重新計(jì)算函數(shù)f(x, y )的值，而迭代又要反的值，而迭代又要反復(fù)進(jìn)行若干次，計(jì)算量很大，而且往往難以預(yù)測復(fù)進(jìn)行若干次，計(jì)算量很大，而且往往難以預(yù)測. 為為了控制計(jì)算量，通常只迭代一兩次就轉(zhuǎn)入下一步的計(jì)了控制計(jì)算量，通常只迭代一兩次就轉(zhuǎn)入下一步的計(jì)算，這算，這就簡化了算法就簡化了算法. .具體地說，我們先用歐拉公式求得一個(gè)初步的近具體地說，我們先用歐拉公式求得一個(gè)初步的近似值似值 ，稱之為，稱之為預(yù)測值預(yù)測值，此預(yù)測值，此預(yù)測值 的精度可能很的精度可能很差，再用梯形公式差，再用梯形公式(2.

17、7)將它校正一次，即按將它校正一次，即按(2.8)式迭式迭代一次，這個(gè)結(jié)果代一次，這個(gè)結(jié)果稱之為稱之為校正值校正值. .1 ny1 ny這樣建立的這樣建立的預(yù)測預(yù)測校正系統(tǒng)校正系統(tǒng)通常稱為通常稱為改進(jìn)的歐拉公式改進(jìn)的歐拉公式:或表為下列或表為下列平均化形式平均化形式 ).(21),(),(11cpnpnncnnnpyyyyxhfyyyxhfyy( (2.9) ), ),(1nnnnyxhfyy 預(yù)測預(yù)測.),(),(2111 nnnnnnyxfyxfhyy校正校正例例2 用改進(jìn)的歐拉法解用改進(jìn)的歐拉法解例例1中的初值問題中的初值問題(2.2). . 解解 仍取步長仍取步長h=0.1, ,改進(jìn)的

18、歐拉公式為改進(jìn)的歐拉公式為 ).(21),2(),2(11cpnpnpncnnnnpyyyyxyhyyyxyhyy部分計(jì)算結(jié)果部分計(jì)算結(jié)果見下表見下表 xnyn 誤差誤差 xnyn 誤差誤差00.20.4 1. 1.184096 1.34336000.0000880.0017190.60.81.01.4859561.6164761.7378690.0027160.0040240.05818同例同例1中的歐拉法的計(jì)算結(jié)果比較，明顯改善了精度中的歐拉法的計(jì)算結(jié)果比較，明顯改善了精度. 例例3 3 (兩種方法的精度比較兩種方法的精度比較)用歐拉公式和改進(jìn)的歐拉公式解下述初值問題，并與用歐拉公式和改進(jìn)

19、的歐拉公式解下述初值問題，并與其準(zhǔn)確解其準(zhǔn)確解y=e- -x+x進(jìn)行比較進(jìn)行比較. . 1) 0(6 . 00, 1yxxyy解解 取步長取步長h=0.1，xk=kh (k=0,1,6). .用兩種方法用兩種方法進(jìn)行計(jì)算對應(yīng)結(jié)果及絕對誤差見下表進(jìn)行計(jì)算對應(yīng)結(jié)果及絕對誤差見下表 xn 歐拉公式歐拉公式 改進(jìn)歐拉公式改進(jìn)歐拉公式 yn 誤差誤差 yn 誤差誤差 0.0 0.1 0.2 0.3 0.4 0.5 0.6 1.000000 1.000000 1.010000 1.029000 1.056100 1.090490 1.13144104.810-38.710-3 1.210-21.410-2

20、1.610-21.810-2 1 1.005000 1.019025 1.041218 1.070802 1.107076 1.14940401.610-42.910-44.010-44.810-45.510-4 5.910-47.2.3 單步法的局部截?cái)嗾`差與階單步法的局部截?cái)嗾`差與階 初值問題初值問題(1.1),(1.2)的的單步法單步法可用可用一般形式一般形式表示為表示為11( , ),(2.10)nnnnnyyhx y yh其中多元函數(shù)其中多元函數(shù) 與與f(x, y )有關(guān)，當(dāng)有關(guān)，當(dāng) 含有含有yn+1時(shí)，方法時(shí)，方法是是隱式隱式的，若不含的，若不含yn+1則為則為顯式方法顯式方法，所

21、以，所以顯式單步顯式單步法法可表示為可表示為 (x, y, h)稱為稱為增量函數(shù)增量函數(shù)，例如對歐拉法，例如對歐拉法(2.1)有有1( , , ),(2.11)nnnnyyh x y h. ),(),(yxfhyxnn 它的局部截?cái)嗾`差已由它的局部截?cái)嗾`差已由(2.3) 給出給出, 對一般顯式單步對一般顯式單步法則可如下定義法則可如下定義. . 定義定義1 設(shè)設(shè)y(x)是初值問題是初值問題(1.1),(1.2)的的準(zhǔn)確解準(zhǔn)確解, 稱稱11()( )( , ( ), ).(2.12)nnnnnTy xy xh x y xh為顯式單步法為顯式單步法( (2.11) )的的局部截?cái)嗾`差局部截?cái)嗾`差.

22、 . Tn+1之所以稱為之所以稱為局部的局部的，是假設(shè)在，是假設(shè)在 xn前各步前各步?jīng)]沒有誤差有誤差. .當(dāng)當(dāng) yn= =y(xn) 時(shí)，計(jì)算一步，則有時(shí)，計(jì)算一步，則有.),(,()()(),()()(11111 nnnnnnnnnnnThxyxhxyxyhyxhyxyyxy 所以，局部截?cái)嗾`差可理解為用方法所以，局部截?cái)嗾`差可理解為用方法(2.11)計(jì)算計(jì)算一步一步的誤差，也即公式的誤差，也即公式(2.11)中用準(zhǔn)確解中用準(zhǔn)確解y(x)代替數(shù)代替數(shù)值解產(chǎn)生的公式誤差值解產(chǎn)生的公式誤差. 根據(jù)定義根據(jù)定義, 顯然顯然歐拉法的局部歐拉法的局部截?cái)嗾`差截?cái)嗾`差為為).()(2)()()()(,()()(3211hOxyhxyhxyhxyxyxhfxyx