數(shù)值分析知識要點

1、第五章 常微分方程數(shù)值解/* Numerical Methods for Ordinary Differential Equations */ 待求解的問題：一階常微分方程的初值問題 /* Initial-Value Problem */: 解的存在唯一性（“常微分方程”理論）：只要 f (x, y) 在a, b R1 上連續(xù)，且關(guān)于 yLipschitz 條件，即存在與 x, y 無關(guān)的常數(shù) L 使對任意定義在 a, b 上的 y1(x) 和 y2(x) 都成立，則上述IVP存在唯一解。1解析解法：（常微分方程理論）只能求解極少一類常微分方程；實際中給定的問題不一定是解析表達(dá)式，而是函數(shù)表，

2、無法用解析解法。如何求解計算解函數(shù) y(x) 在一系列節(jié)點 a = x0 x1 xn= b 處的近似值節(jié)點間距 為步長，通常采用等距節(jié)點，即取 hi = h (常數(shù))。數(shù)值解法: 求解所有的常微分方程步進(jìn)式：根據(jù)已知的或已求出的節(jié)點上的函數(shù)值計算當(dāng)前節(jié)點上的函數(shù)值，一步一步向前推進(jìn)。因此只需建立由已知的或已求出的節(jié)點上的函數(shù)值求當(dāng)前節(jié)點函數(shù)值的遞推公式即可。234-Eulers Method1 歐拉方法 /* Eulers Method */1 Eulers MethodTaylor展開法5幾何意義亦稱為歐拉折線法 /* Eulers polygonal arc method*/ 幾何直觀是幫

3、助我們尋找解決一個問題的思路的好辦法哦6定義在假設(shè) yn = y(xn)，即第 n 步計算是精確的前提下，考慮公式或方法本身帶來的誤差: Rn = y(xn+1) yn+1 , 稱為局部截斷誤差 /* local truncation error */。說明 顯然，這種近似有一定誤差，而且步長越大，誤差越大，如何估計這種誤差y(xn+1) yn+1 ？1 Eulers Method7截斷誤差: 實際上，y(xn) yn， yn 也有誤差，它對yn+1的誤差也有影響，見下圖。但這里不考慮此誤差的影響，僅考慮方法或公式本身帶來的誤差，因此稱為方法誤差或截斷誤差。局部截斷誤差的分析：由于假設(shè)yn =

4、 y(xn) ，即yn準(zhǔn)確，因此分析局部截斷誤差時將y(xn+1) 和 yn+1都用點xn上的信息來表示，工具：Taylor展開。 歐拉法的局部截斷誤差：Rn+1 的主項/* leading term */1 Eulers Method8定義若某算法的局部截斷誤差為O(hp+1)，則稱該算法有p 階精度。 歐拉法具有 1 階精度。在xn點用一階向前差商近似一階導(dǎo)數(shù) 在第二章討論牛頓插值公式時介紹了差商的概念和性質(zhì)，各階差商可以近似各階導(dǎo)數(shù)，具有不同的精度，且可以用函數(shù)值來表示。上一章中數(shù)值微分的方法之一就是用差商近似導(dǎo)數(shù)Eulers method1 Eulers Method91 Eulers

5、 Method 歐拉公式的改進(jìn)：隱式歐拉法或后退歐拉法 /* implicit Euler method or backward Euler method*/xn+1點向后差商近似導(dǎo)數(shù)隱式或后退歐拉公式10由于未知數(shù) yn+1 同時出現(xiàn)在等式的兩邊，故稱為隱式 /* implicit */ 歐拉公式，而前者稱為顯式 /* explicit */ 歐拉公式。隱式公式不能直接求解，一般需要用Euler顯式公式得到初值，然后用Euler隱式公式迭代求解。因此隱式公式較顯式公式計算復(fù)雜，但穩(wěn)定性好（后面分析）。收斂性1 Eulers Method11 見上圖， 顯然，這種近似也有一定誤差，如何估計這種

6、誤差y(xn+1) yn+1 ？方法同上，基于Taylor展開估計局部截斷誤差。但是注意，隱式公式中右邊含有f(xn+1 , yn +1 ) ,由于yn +1不準(zhǔn)確，所以不能直接用y (xn+1)代替f(xn+1 , yn +1 ) 設(shè)已知曲線上一點 Pn (xn , yn ),過該點作弦線，斜率為(xn+1 , yn +1 ) 點的方向場f(x,y)方向,若步長h充分小，可用弦線和垂線x=xn+1的交點近似曲線與垂線的交點。幾何意義xnxn+1PnPn+1xyy(x)1 Eulers Method12 隱式歐拉法的局部截斷誤差：1 Eulers Method131 Eulers Method

7、14 隱式歐拉法的局部截斷誤差：即隱式歐拉公式具有 1 階精度。1 Eulers Method15比較尤拉顯式公式和隱式公式及其局部截斷誤差顯式公式隱式公式1 Eulers Method16 若將這兩種方法進(jìn)行算術(shù)平均，即可消除誤差的主要部分/*leading term*/而獲得更高的精度,稱為梯形法1 Eulers Method 梯形公式 /* trapezoid formula */ 顯、隱式兩種算法的平均注：的確有局部截斷誤差 ， 即梯形公式具有2 階精度，比歐拉方法有了進(jìn)步。但注意到該公式是隱式公式，計算時不得不用到迭代法，其迭代收斂性與歐拉公式相似。17梯形法的迭代計算和收斂性收斂性

8、1 Eulers Method18梯形法的簡化計算 迭代計算量大，且難以預(yù)測迭代次數(shù)。為了控制計算量，通常只迭代一次就轉(zhuǎn)入下一點的計算。用顯式公式作預(yù)測，梯形公式作校正，得到如下預(yù)測校正系統(tǒng)，也稱為改進(jìn)尤拉法： 改進(jìn)歐拉法 /* modified Eulers method */Step 1: 先用顯式歐拉公式作預(yù)測，算出),(1nnnnyxfhyy+=+Step 2: 再將 代入隱式梯形公式的右邊作校正，得到1+ny),(),(2111+=nnnnnnyxfyxfhyy1 Eulers Method19注：此法亦稱為預(yù)測-校正法 /* predictor-corrector method *

9、/?？梢宰C明該算法具有 2 階精度，同時可以看到它是個單步遞推格式，比隱式公式的迭代求解過程簡單。后面將看到，它的穩(wěn)定性高于顯式歐拉法。1 Eulers Method其它形式20幾何解釋xnxn+1ABPn+1=(A+B)/2尤拉法后退尤拉法梯形法1 Eulers Method21令x=x1,得Another point of view對右端積分采用左矩形、右矩形、梯形積分公式，即可得尤拉顯式、隱式、梯形公式1 Eulers Method22 中點歐拉公式 /* midpoint formula */中心差商近似導(dǎo)數(shù)x0 x2x1假設(shè) ，則可以導(dǎo)出即中點公式也具有 2 階精度，且是顯式的。需要

10、2個初值 y0和 y1來啟動遞推過程，這樣的算法稱為雙步法 /* double-step method */，而前面的三種算法都是單步法 /* single-step method */。1 Eulers Method23幾何解釋xnxn+1尤拉法后退尤拉法中點法xn-1令x=x2,得Another point of view對右端積分采用中矩形公式即得中點公式1 Eulers Method24 預(yù)測-校正-改進(jìn)系統(tǒng)中點法具有二階精度，且是顯式的，與梯形公式精度相匹配，用中點公式作預(yù)測，梯形公式作校正，得到如下預(yù)測校正系統(tǒng)：校正誤差約為預(yù)測誤差的1/41 Eulers Method25預(yù)測誤差

11、和校正誤差的事后誤差估計式利用上兩式可以估計預(yù)測值和校正值與準(zhǔn)確值的誤差，可以期望，利用這兩個誤差分別作預(yù)測值和校正值的補償，有可能提高精度。 設(shè)pn,cn分別為第n步的預(yù)測值和校正值，即此時cn+1未知，故用pn -cn代替1 Eulers Method26 預(yù)測-校正-改進(jìn)公式注：利用該算法計算yn+1時，需要1 Eulers Method27summary28兩個預(yù)測-校正系統(tǒng)尤拉兩步法和梯形公式構(gòu)成的預(yù)測-校正-改進(jìn)系統(tǒng)尤拉公式和梯形公式構(gòu)成的預(yù)測-校正系統(tǒng)293031例 32HW: p.201 #1-5證明中點法和梯形公式的精度為2階332 龍格 - 庫塔法 /* Runge-Kut

12、ta Method */建立高精度的單步遞推格式:在改進(jìn)尤拉法和尤拉兩步法預(yù)測-校正系統(tǒng)中,預(yù)測公式都是單步法,如果預(yù)測誤差很小,則通過校正后得到的近似值誤差會更小,因此需要研究高精度的單步法. 1. Taylor級數(shù)法IVP:設(shè)其解為y=y(x)由Taylor展開，有(1)2 Runge-Kutta Method34(2)2 Runge-Kutta Method35要使公式具有p階精度，則在(1)式中截取前p+1項，用(2)式計算各階導(dǎo)數(shù)，即得下面Taylor公式：局部截斷誤差(3)2 Runge-Kutta Method362.Taylor公式(3)表面上看形式簡單，但具體構(gòu)造時往往很困難

13、，因為按(2)式求導(dǎo)，這一過程可能很復(fù)雜。因此通常不直接用Taylor公式，而借鑒其思想提出其它公式。1. 由此看出，一種方法具有p階精度公式對不超過p次的多項式準(zhǔn)確成立（局部截斷誤差為0）。這一等價條件也可以用來判斷一種方法的精度。2 Runge-Kutta Method37單步遞推法的基本思想是從 ( xn , yn ) 點出發(fā)，以某一斜率沿直線達(dá)到 ( xn+1 , yn+1 ) 點。歐拉法及其各種變形所能達(dá)到的最高精度為2階。 2. RungeKutta Method由微分中值定理，有k*稱為區(qū)間xn, xn+1上的平均斜率，只要知道平均斜率，就可計算y(xn+1).因此只要對平均斜率

14、提供一種近似算法，則由(4)式可導(dǎo)出一種相應(yīng)的求解公式。(4) 2 Runge-Kutta Method38例2 Runge-Kutta Method39 由此看出，改進(jìn)的尤拉公式用xn與xn+1兩個節(jié)點的斜率的算術(shù)平均作為平均斜率， xn+1點的斜率通過已知信息yn來預(yù)測。 考察改進(jìn)的歐拉法，可以將其改寫為：斜率一定取K1， K2 的平均值嗎？步長一定是h 嗎？即第二個節(jié)點一定是xn+1嗎？2 Runge-Kutta Method402 Runge-Kutta Method首先希望能確定系數(shù) 1、2、p，使得到的算法格式有2階精度，即在 的前提假設(shè)下，使得 Step 1: 將 K2 在 (

15、xi , yi ) 點作 Taylor 展開將改進(jìn)歐拉法推廣為：),(),(12122111phKyphxfKyxfKKKhyyiiiiii+=+=+llStep 2: 將 K2 代入第1式，得到 2階RungeKutta Method412 Runge-Kutta MethodStep 3: 將 yi+1 與 y( xi+1 ) 在 xi 點的泰勒展開作比較要求 ，則必須有：這里有 個未知數(shù)， 個方程。32存在無窮多個解。所有滿足上式的格式統(tǒng)稱為2階龍格 - 庫塔格式。注意到， 就是改進(jìn)的歐拉法；p=1/2, 1=0, 2=1, 變形尤拉公式。42Q: 為獲得更高的精度，應(yīng)該如何進(jìn)一步推廣？

16、改進(jìn)的Euler 公式推廣為二階Runge-Kutta公式帶來這樣的啟示：若在xn, xn+1上多預(yù)測幾個點的斜率值，然后將它們的算術(shù)平均作為平均斜率，則有可能構(gòu)造出具有更高精度的計算公式。-Runge-Kutta方法的基本思想。注：二階Runge-Kutta公式用多算一次函數(shù)值f 的辦法避開了二階Taylor級數(shù)法所要計算的f 的導(dǎo)數(shù)。在這種意義上，可以說Runge-Kutta方法實質(zhì)上是Taylor級數(shù)法的變形。2 Runge-Kutta Method43其中i ( i = 1, , m )，i ( i = 2, , m ) 和 ij ( i = 2, , m; j = 1, , i1 )

17、 均為待定系數(shù)，確定這些系數(shù)的步驟與前面相似。 2 Runge-Kutta Method).,(.),(),(),(.1122112321313312122122111-+=+=+=+=mm mmmmimiiiiiimmiihKhKhKyhxfKhKhKyhxfKhKyhxfKyxfKKKKhyybbbabbaballl 高階RungeKutta Method4445 Gill公式：4階經(jīng)典龍格-庫塔公式的一種改進(jìn)2 Runge-Kutta Method 最常用為四級4階經(jīng)典龍格-庫塔法 /* Classical Runge-Kutta Method */ ：462 Runge-Kutta M

18、ethod注： 龍格-庫塔法的主要運算在于計算 Ki 的值，即計算 f 的值。Butcher 于1965年給出了計算量與可達(dá)到的最高精度階數(shù)的關(guān)系：753可達(dá)到的最高精度642每步須算Ki 的個數(shù) 由于龍格-庫塔法的導(dǎo)出基于泰勒展開，故精度主要受解函數(shù)的光滑性影響。對于光滑性不太好的解，最好采用低階算法而將步長h 取小。4748492 Runge-Kutta Method 變步長的RungeKutta MethodQ: 由局部截斷誤差可以看出，步長 h 越小，局部截斷誤差越?。坏介L減小，在一定求解范圍（區(qū)間）內(nèi)要完成的步數(shù)就增加了，步數(shù)增加會引起計算量增大，導(dǎo)致舍入誤差積累。因此要選取適當(dāng)?shù)?/p>

19、步長。選擇步長時要考慮兩個問題： 1.如何衡量和檢驗計算結(jié)果的精度？ 2.如何根據(jù)所獲得的精度處理步長？5051HW: p.201 #6-8523 單步法的收斂性與穩(wěn)定性 /* Convergency and Stability */ 前面介紹了兩大類微分方程數(shù)值解法：一類是用差商近似導(dǎo)數(shù)得到的尤拉系列公式，另一類是基于平均斜率概念的RungeKutta公式。基本思想都是通過某種離散化手續(xù)，將微分方程轉(zhuǎn)化為差分方程（代數(shù)方程）來求解。 Q1. 這種轉(zhuǎn)化是否合理？要看差分問題的解yn當(dāng)h0時是否收斂到微分方程的解y(xn),即是否成立 yn y(xn)， h0. -收斂性問題 Q2. 實際計算時

20、，由于舍入誤差的影響，差分方程的解本身也有誤差，這種誤差在計算過程中會不會擴(kuò)大 -穩(wěn)定性問題3 Convergency and Stability53 收斂性 /* Convergency */定義 若某算法對于任意固定的 x = xi = x0 + i h，當(dāng) h0 ( 同時 i ) 時有 yi y( xi )，則稱該算法是收斂的。 例：就初值問題 考察歐拉顯式格式的收斂性。解：該問題的精確解為 歐拉公式為對任意固定的 x = xi = i h ，有3 Convergency and Stability54顯式單步法的收斂性(1)3 Convergency and Stability55而整

21、體截斷誤差為證明: 設(shè) 表示當(dāng)yn =y(xn)時， 由公式(1)求得的結(jié)果，即則局部截斷誤差為(2)(2)-(1),得3 Convergency and Stability56從而3 Convergency and Stability573 Convergency and Stability58 Euler法的收斂性 ： (x,y)=f(x,y),故當(dāng)f(x,y)滿足Lipschitz條件時，尤拉法收斂；3 Convergency and Stability593 Convergency and Stability 穩(wěn)定性 /* Stability */例：考察初值問題 在區(qū)間0, 0.5上

22、的解。分別用歐拉顯、隱式格式和改進(jìn)的歐拉格式計算數(shù)值解。0.00.10.20.30.40.5精確解改進(jìn)歐拉法 歐拉隱式歐拉顯式 節(jié)點 xi 1.00002.0000 4.00008.0000 1.6000101 3.2000101 1.00002.5000101 6.25001021.56251023.90631039.76561041.00002.50006.25001.56261013.90631019.76561011.00004.97871022.47881031.23411046.14421063.0590107What is wrong ?! An Engineer complai

23、ns: Math theorems are so unstable that a small perturbation on the conditions will cause a crash on the conclusions!603 Convergency and Stability613 Convergency and Stability定義若某算法在計算過程中任一步產(chǎn)生的誤差在以后的計算中都逐步衰減，則稱該算法是絕對穩(wěn)定的 /*absolutely stable */。一般分析某算法的穩(wěn)定性時，為簡單起見，只考慮模型方程或試驗方程 /* test equation */常數(shù)，可以是復(fù)

24、數(shù)當(dāng)步長取為 h 時，將某算法應(yīng)用于上式，并假設(shè)只在初值產(chǎn)生誤差 ，則若此誤差以后逐步衰減，就稱該算法相對于 絕對穩(wěn)定， 的全體構(gòu)成絕對穩(wěn)定區(qū)域。我們稱算法A 比算法B 穩(wěn)定，就是指 A 的絕對穩(wěn)定區(qū)域比 B 的大。hl h=h623 Convergency and Stability633 Convergency and Stability643 Convergency and Stability653 Convergency and Stability例：隱式龍格-庫塔法而顯式 1 4 階方法的絕對穩(wěn)定區(qū)域為其中2階方法 的絕對穩(wěn)定區(qū)域為0ReImgk=1k=2k=3k=4-1-2-3-1

25、23ReImg無條件穩(wěn)定注：一般來說，隱式歐拉法的絕對穩(wěn)定性比同階的顯式法的好。66小結(jié)67 尤拉兩步公式與尤拉單步公式相比，使用兩個節(jié)點上的已知信息將精度提高一階?？梢栽O(shè)想，計算y(xn+1)時，充分利用前面已經(jīng)求出的節(jié)點上的 y 及 y 值的線性組合來近似y(xn+1)，精度會大大提高。-線性多步法的基本思想。 構(gòu)造線性多步法有多種途徑，這里介紹兩種： 基于數(shù)值積分的構(gòu)造方法； 基于Taylor展開的構(gòu)造方法。4 線性多步法 /* Multistep Method */68).(.110111101rnrnnnrnrnnnffffhyyyy-+-+=bbbbaaa線性多步法的通式可寫為：當(dāng)

26、 10 時，為隱式公式; 1=0 則為顯式公式。 基于數(shù)值積分的構(gòu)造法將 在 上積分，得到只要近似地算出右邊的積分 ，則可通過 近似y(xn+1) 。而選用不同近似式 I，可得到不同的計算公式。例如利用左矩形積分公式得到尤拉公式；梯形積分公式得到梯形公式。69一般地，利用插值原理所建立的一系列數(shù)值積分方法也可以導(dǎo)出解微分方程的一系列計算公式。運用插值方法的關(guān)鍵在于選取合適的插值節(jié)點。假設(shè)已構(gòu)造出f(x,y(x)的插值多項式Pr(x),則704 Multistep Method 亞當(dāng)姆斯顯式公式 /* Adams explicit formulae */利用r+1 個節(jié)點上的被積函數(shù)值構(gòu)造 r

27、階牛頓后插多項式, 有Newton插值余項/*Adams 顯式公式 */外推技術(shù) /* extrapolation */71table 5-6, p.181jj0111/225/1233/8實際計算時，將差分展開 ijrj=i72局部截斷誤差為：注：Br 與yn+1 計算公式中 fn , , fnk 各項的系數(shù) 均可查表得到 。 見下表。7310123rfnfn1fn2fn3Br74例：常用的是 r = 3 的4階亞當(dāng)姆斯顯式公式75Adams顯式公式用 作為插值節(jié)點，在求積區(qū)間xn, xn+1上用插值函數(shù)Nr(x)近似f(x,y(x)，而xn+1不在插值節(jié)點內(nèi)，因此是一個外推的過程。雖然公式

28、是顯式的，便于計算，但效果并不理想，比如穩(wěn)定性較差等。因此改用通過r+1個節(jié)點的插值多項式Nr(x)近似f(x,y(x)，由于xn+1是其中一個插值點，因此是內(nèi)插多項式，但導(dǎo)出的公式是隱式的。764 Multistep Method 亞當(dāng)姆斯隱式公式 /* Adams implicit formulae */利用r+1 個節(jié)點上的被積函數(shù)值 fn+1 , fn , , fnr+1 構(gòu)造r 階牛頓前插多項式。與顯式多項式完全類似地可得到一系列隱式公式。局部截斷誤差為：7710123kfi+1fifi1fi2Br小于Br注： 與yn+1 計算公式中 fn+1 , , fn+1k 各項的系數(shù) 均可查

29、表得到 。 見下表。78常用的是 k = 3 的4階亞當(dāng)姆斯隱式公式較同階顯式穩(wěn)定例：794 Multistep Method 亞當(dāng)姆斯預(yù)測-校正系統(tǒng) /* Adams predictor-corrector system */Step 1: 用Runge-Kutta 法計算前 k 個初值；Step 2: 用Adams 顯式計算預(yù)測值；Step 3: 用同階Adams 隱式計算校正值。注意：三步所用公式的精度必須相同。通常用經(jīng)典Runge-Kutta 法配合4階Adams 公式。4階Adams隱式公式的截斷誤差為4階Adams顯式公式的截斷誤差為Predicted value pn+1Corr

30、ected value cn+180Modified value mn+1預(yù)測值與校正值的事后誤差估計81Adams顯隱預(yù)測-校正-改進(jìn)系統(tǒng)824 Multistep Method Adams 4th-Order predictor-corrector AlgorithmTo approximate the the solution of the initial-value problemAt (N+1) equally spaced numbers in the interval a, b.Input: endpoints a, b; integer N; initial value y0

31、.Output: approximation y at the (N+1) values of x.Step 1 Set h = (b a) / N ; x0 = a; y0 = y0; Output ( x0, y0 );Step 2 For i = 1, 2, 3 Compute yi using classical Runge-Kutta method; Output ( xi , yi );Step 3 For i = 4, , N do steps 4-10Step 5 ; /* predict */ Step 6 ; /* modify */Step 7 ; /* correct

32、*/Step 8 ; /* modify the final value */Step 9 Output ( xi+1 , yi+1 ); Step 10 For j = 0, 1, 2, 3 Set xi = xi+1 ; yi = yi+1 ; /* Prepare for next iteration */Step 11 STOP. 應(yīng)為( ci+1 pi+1 ), 但因ci+1 尚未算出，只好用( ci pi )取代之。834 Multistep Method 基于泰勒展開的構(gòu)造法).(.110111101rnrnnnrnrnnnffffhyyyy-+-+=bbbbaaa 將通式中的右

33、端各項 yn1, , ynk ; fn+1, fn1, , fnk (yn+1, , ynk)分別在 xn 點作泰勒展開，與精確解 y(xn+1) 在 xn 點的泰勒展開作比較。通過令同類項系數(shù)相等，得到足以確定待定系數(shù)0, , r ; 1, 0, , r 的等式，則可構(gòu)造出線性多步法的公式。當(dāng) 10 時，為隱式公式; 1=0 則為顯式公式。84局部截斷誤差為85例：設(shè))(3322110221101-+=nnnnnnnnyyyyhyyyybbbbaaa確定式中待定系數(shù)0, 1, 2, 0, 1, 2, 3, 使得公式具有4階精度。4 Multistep Method解：/* y(xn) = y

34、n */個未知數(shù)個方程75864 Multistep Method 令 1 = 2 = 0Adams 顯式公式 以 yn+1 取代 yn3，并取 1 = 2 = 0Adams 隱式公式 以 yn3 取代 yn3 ，則可導(dǎo)出另一組4 階顯式算法，其中包含了著名的米爾尼 /* Milne */ 公式(4步4階顯式公式)其局部截斷誤差為注：上式也可通過數(shù)值積分導(dǎo)出，即將 在區(qū)間 上積分，得到 再過 做 f 的2次插值多項式即可。取 1 = 1, 2 = 0得到辛甫生 /* Simpson */ 公式與Milne 公式匹配使用,構(gòu)成預(yù)測-校正系統(tǒng)辛甫生 /* Simpson */ 公式（2步4階）在區(qū)

35、間xn1, xn+1上積分，并用Simpson數(shù)值積分公式來近似積分項，亦可得此Simpson公式。87其局部截斷誤差為Milne-Simpson預(yù)測校正系統(tǒng)的局部截斷誤差(14/450.31, -1/90 -0.01)都比同價的Adams公式 (251/720 0.34, -19/720 -0.02)小，計算量又比4階Runge-Kutta公式少。缺點是, 校正公式是弱穩(wěn)定的，即對某些問題來說，其穩(wěn)定性不好。Hamming對Simpson公式的穩(wěn)定性作了改進(jìn)，得到Hamming公式。884 Multistep Method Milne-Simpson 系統(tǒng)的缺點是穩(wěn)定性差，為改善穩(wěn)定性，考慮

36、另一種隱式校正公式：要求公式具有4 階精度。通過泰勒展開，可得到 個等式，從中解出 個未知數(shù)，則有 個自由度。561取 1 = 1 得Simpson 公式哈明 /* Hamming */ 用1 的不同數(shù)值進(jìn)行試驗，發(fā)現(xiàn)當(dāng)1 = 0 時，公式的穩(wěn)定性較好，即：其局部截斷誤差為注：哈明公式不能用數(shù)值積分方法推導(dǎo)出來。89Milne-Hamming預(yù)測-校正-改進(jìn)系統(tǒng)905 微分方程組與高階方程 /* Systems of Differential Equations and Higher-Order Equations */ 一階微分方程組IVP的一般形式為：=)(,.),(,()(.)(,.),

37、(,()(1111xyxyxfxyxyxyxfxymmmm初值0002020101)(,.,)(,)(mmyxyyxyyxy=將問題記作向量形式，令：前述所有公式皆適用于向量形式。91 高階微分方程5 Systems of DEs and Higher-Order Equations=-10)1(1000)1()()(,.,)(,)(),.,(nnnnaxyaxyaxyyyyxfy化作一階微分方程組求解。引入新變量初值條件為：926 邊值問題的數(shù)值解 /* Boundary-Value Problems */2 階常微分方程邊值問題 打靶法 /* shooting method */先猜測一個初始斜率 y (a) = s，通過解初值問題y(b) = (s)找出s*使得(s*) = ，即把問題轉(zhuǎn)化為求方程 (s) = 0 的根。yx0abyx()b斜率=s0j()s0斜率=s1j()s1每計算一個(s) 都必須解一個ODE.93 有限差分法 /* finite di