微分方程數(shù)值解

1、第6章 常微分方程初值問題數(shù)值解法本章探討常微分方程特解的常用數(shù)值方法的構造和原理，主要介紹求常微分方程初值問題的常用方法和有關知識。重點論述Euler方法、Runge-Kutta方法和線性多步法的原理、構造、局部截斷誤差和穩(wěn)定性等內容。6.1 實際案例工程技術里某些振動問題可以表示為單擺的運動，其運動規(guī)律的微分方程為：怎樣求出其特解？該微分方程不能用通常的解析方法來求解！怎樣解不能用解析方法求解的微分方程特解問題是本章要解決的問題。6.2 問題的描述和基本概念1、常微分方程初值問題l 一般形式 （6.1）式中已知，稱為初值條件。l 初值問題的數(shù)值方法和數(shù)值解求函數(shù)在若干離散點上的近似值的方法

2、稱為初值問題的數(shù)值方法，而稱為初值問題的數(shù)值解。2. 建立數(shù)值解法的思想與方法微分方程初值問題的數(shù)值解法是用離散化方法將初值問題化為差分方程后再求解的方式。設區(qū)間a,b上的一組節(jié)點為距離稱為步長。求數(shù)值解一般是從開使逐次順序求出。初值問題的解法有單步法和多步法兩種：l 單步法：計算時只用到一個值；l 多步法：計算時要用多個值。數(shù)值解法還有顯格式和隱格式之分。l 微分方程離散化方法主要有數(shù)值微分法，數(shù)值積分法和Taylor展開法。l 初值問題化為差分方程的方法1. 用離散方法去掉方程中的導數(shù)得到近似離散化方程；2. 在近似離散化方程中用代替；3. 在近似離散化方程中將近似號“”用等號“=”代替。

3、1) 數(shù)值微分法由初值問題（6.1）有，用數(shù)值微分的2點前差公式代替，得近似離散化方程記，做，“”，得差分方程寫容易計算的形式 （6.2）（Euler公式）由初值條件及式（6.2）可求出（6.1）的數(shù)值解。公式（6.2）是顯式單步法。2）數(shù)值積分法在上對兩邊取定積分，得對右端積分采用梯形公式（數(shù)值積分公式）得近似離散化方程：于是得到求初值問題（6.1）的梯形方法該公式是隱式單步法。3）Taylor展開法因為初值問題中函數(shù)是已知函數(shù)，由，可以計算，于是有函數(shù)在處的Taylor展式取上式右端前若干項，得近似離散化方程。例如取前兩項有于是又得到Euler公式：6.3數(shù)值解法的誤差、階與絕對穩(wěn)定性單步

4、法數(shù)學描述為顯式： （6.4）隱式：（6.5）其中稱為增量函數(shù)，與有關。l 顯式單步法的一些概念定義6.1 設是問題（6.1）的解，是經過式（6.4）求出的的計算解，則稱為單步法（6.4）在節(jié)點的整體截斷誤差，而稱 （6.6）為在點的局部截斷誤差。表示解在的值，是準確值，沒有誤差；表示由數(shù)值解公式得出的近似值，是數(shù)值解，有截斷誤差；表示用計算機計算給出的計算解，有舍入誤差。l 局部截斷誤差的理解假設在計算時沒有誤差（）下，由式（6.4）計算出的（）與的誤差。整體截斷誤差還要加上與的誤差?？疾斐踔祮栴}解法的優(yōu)劣，引入階的概念。定義6.2 如果初值問題（6.1）對某種數(shù)值解法的局部截斷誤差為則稱該

5、方法具有p階精度或該方法是p階方法。方法的階越高，方法越好。l 主局部截斷誤差或局部截斷誤差的主項如果某方法是p階方法，按可展為則稱為主局部截斷誤差。在同階方法中，主局部截斷誤差越小，方法越好。例6.1常微分方程初值問題的單步法為試求其局部截斷誤差主項并回答它是幾階精度的？解 該單步公式的局部截斷誤差是故局部截斷誤差主項是，方法是一階的。l 求階p的另一方法因為去掉下標，有若將在x點展開有則知該方法的階是p。例如，對Euler方法，有那么將在x點展開，有故有因此，Euler方法是一階方法。定義6.3 設用某種數(shù)值方法求初值問題（6.1）在任意節(jié)點的數(shù)值解時，滿足，則稱該數(shù)值方法是絕對穩(wěn)定的。這

6、里是計算機計算時得出的計算解的舍入誤差，。通常用試驗方程 （為復數(shù)）來討論求解初值問題的數(shù)值方法絕對穩(wěn)定性；對具體初值問題，可取。穩(wěn)定性常與步長h有關。定義6.4 某方法在試驗方程中絕對穩(wěn)定的復平面范圍稱為該方法的絕對穩(wěn)定域，它與復平面實軸的交稱為該方法的絕對穩(wěn)定區(qū)間。絕對穩(wěn)定域包含復平面左半平面的方法稱為是A-穩(wěn)定的。絕對穩(wěn)定域越大，對應的方法絕對穩(wěn)定性越好。6.4 Euler方法的有關問題用Euler公式，求解初值問題（6.1）數(shù)值解的方法稱為Euler方法。1 )Euler方法的幾何意義Euler方法常稱為折線法。2) Euler方法的誤差設為的計算解，滿足其中。則有Euler方法的局部

7、截斷誤差；Euler方法的總體截斷誤差由，,有因為對任意m都有，可得由k的任意性，可知Euler方法的總體截斷誤差為說明當h足夠小時，由Euler方法計算所得數(shù)值解可以很好地逼近準確解，從而Euler方法是收斂的。3) Euler方法穩(wěn)定性將Euler公式用于試驗方程，得到 （6.8）設計算時有舍入誤差，則有得要想，只須，因此Euler方法在時是絕對穩(wěn)定的，其絕對穩(wěn)定域為復平面上以-1為中心的單位圓盤。由圖形可知Euler方法的絕對穩(wěn)定區(qū)間為。-1-2圖6.3圖6.4若指定是負實數(shù)，則有當步長h滿足時，可保證Euler方法的計算絕對穩(wěn)定。例6.2 ，取用Euler方法求其數(shù)值解并與其準確解作比

8、較。 解 本題準確解為，Euler法計算公式為直接計算結果與對應準確值的誤差有如下表xk0.0250.050.0750.1yk-1.52.25-3.3755.0625y(xk)- yk1.58-2.243.37-5.06可見計算結果波動大，不穩(wěn)定，計算結果不可用。原因：，h的穩(wěn)定性范圍，而本題的h=0.0250.02，故計算結果不穩(wěn)定。4) 改進的Euler方法易證梯形方法（6.3）的絕對穩(wěn)定域為復平面上的部分，它是復平面的整個左半平面（圖6.4）。因此是A-穩(wěn)定的。其絕對穩(wěn)定域比Euler方法大，是二階方法，故梯形方法比Euler方法好，但美中不足的是梯形方法是隱式方法。為避免迭代，實用中采

9、取如下兩次計算格式由來計算（改進的Euler方法）：預測 （6.9）校正 （6.10）預測校正公式。易證它也是二階方法。6.5 Runge-Kutta方法理論上，由Taylor展開法可以構造出解初值問題（6.1）的高階數(shù)值方法，但這涉及到要計算的高階導數(shù)，因此很不方便。本節(jié)的Runge-Kutta方法從函數(shù)本身著手，避開計算的高階導數(shù)，來構造高階數(shù)值方法?；舅枷雽⑽⒎址匠坛踔祮栴}轉化為積分方程問題，再對積分方程中的定積分使用待定的m點插值型求積公式構造高階的函數(shù)展開模式以獲得高階數(shù)值方法。1. 構造原理將初值問題的方程在上積分對定積分用m個點的插值型求積公式取，并對中的做展開，得解初值問題的

10、Runge-Kutta一般公式形式（6.11）利用Taylor展開式對這些參數(shù)適當選擇就可以構造高階方法了。Runge-Kutta方法（6.11）的增量函數(shù)是要構造p階數(shù)值方法，可選擇參數(shù)使局部截斷誤差即可。2 構造過程以來說明Runge-Kutta方法的構造方法和過程，對一般的Runge-Kutta方法可類似處理。的Runge-Kutta公式為 （6.12）式中 ，。其增量函數(shù)為引進符號， ，由，可得，在處做Taylor展開，有對在做二元Taylor展開，有由，有選 （6.13）有局部截斷誤差，這樣從式（6.13）中解出，代入到式（6.12）中，即得到一組二階Runge-Kutta公式。由于

11、式（6.13）中有3個方程4個參數(shù)，若令一個參數(shù)為自由變量，如取，則式（6.13）的解為，取不同的可得出不同的二階Runge-Kutta公式。如取時，得到改進的Euler公式3. Runge-Kutta方法的階與級的關系計算函數(shù)值的次數(shù)稱為級。l Runge-Kutta方法的級與階對應關系：計算的次數(shù) 1 2 3 4 5 6 7對應方法的最高階 1 2 3 4 4 5 6l 經典Runge-Kutta公式 （6.14）例 設初值問題為分別用Euler方法（），改進Euler方法（）和經典Runge-Kutta方法（）計算。解 Euler方法計算格式（）為改進的Euler方法計算格式（）為經典R

12、unge-Kutta方法計算格式（）為它們的初值，計算結果及準確解列于下表Euler方法改進Euler方法經典R-K法 000000.10.096 3120.095 1230.095 162500.095 162 580.20.183 3480.181 1930.181 269 100.181 269 250.30.262 0010.259 0850.259 181 580.259 181 780.40.333 0790.329 5630.329 679 710.329 679 95例6.5給定初值問題1）分析求解公式的局部截斷誤差，指出它是幾階公式；2）證明用上面求解公式計算初值問題的數(shù)值解

13、成立極限 本題中的節(jié)點是等距節(jié)點，h為步長，n為由節(jié)點分割區(qū)間a,b的份數(shù)。解 由題意有1）局部截斷誤差將在點做Taylor展開到項，將在點做二元Taylor展開到項，則有得所用公式是2階的。2）顯然所給初值問題的準確解為。由給出的數(shù)值解計算公式有故式中6.6線性多步法線性多步法的一般計算格式為 （6.15）式中均為常數(shù),，為等距節(jié)點，步長為h。若不同時為零，計算一個需要用到的前個值,式（6.15 ）稱為線性步法。當n>1時就稱為線性多步方法。構造解初值問題（6.1）的線性多步法有基于數(shù)值積分方法和Taylor展開方法兩種手段。1. 基于數(shù)值積分的構造方法1) 構造原理與分析將初值問題在

14、積分，m為正整數(shù)，得然后選擇的不同插值函數(shù)代替并積分之，即可得到計算公式。l Adams方法給出構造過程。已初值問題的解在個點上的數(shù)值解為，記，考慮數(shù)表構造的次Lagrange插值多項式，在式（6.16）中取=1，有因為式中 。記，有于是可得n步Adams顯示（外插）公式 （6.17）其局部截斷誤差為式中的定積分計算利用了變換。常用的Adams公式1、4階Adams顯式公式局部截斷誤差 2、4階Adams隱式公式局部截斷誤差 例6.7 構造3階Adams顯示公式。解 取三個相鄰的等距節(jié)點對應的數(shù)據做二次Lagrange插值多項式注意到，作變換，有故有于是得到3階Adams外推公式其局部截斷誤差

15、為2. 基于Taylor展開的構造方法l 先給出線性多步法的計算公式模式l 然后對局部截斷誤差表達式在處作Taylor展開以確定公式的系數(shù)。例 6.8 設求初值問題的線性3步公式具有如下形式h為步長，試求系數(shù)a,b,c使該公式的階數(shù)盡可能高，并寫出其局部截斷誤差。解 局部截斷誤差為因為有3個待定系數(shù)，選擇截斷誤差的前3項的系數(shù)為零，得關于待定系數(shù)的線性方程組該方程組有唯一解 ，用此值代入第四項的系數(shù)，有故有當時，所給公式的階數(shù)達到最高，其值為3，對應的局部截斷誤差為6.7步長的自動選取問題數(shù)值解法步長h不變，會產生求出的各節(jié)點上的數(shù)值的解具有不同的精度的情況！要使所有數(shù)值解都滿足給定的精度，必

16、須選用較小的步長h，但這又會產生計算量的增加和舍入誤差的累積！因此最好的方式是步長可以變化。通常計算機上多采用根據誤差估計來自動選步長值的方式處理這個問題。l 自動選步長方法的原理和過程設某種方法的階為p>0, 設為步長。記表示由計算一步算出的在的數(shù)值解，有如將步長折半為，用表示由計算兩個半步長算出的在的數(shù)值解，有假設，則有整理后，有記，則在給定精度后，若有，就有說明用步長折半前后在的兩次數(shù)值解之差的大小可判別在處數(shù)值解的精度！步長自動選取的算法1. 先選定某個步長，然后計算；2. 若，（說明精度足夠，但為防止精度太高導致求數(shù)值解的計算量太大，做步長擴大處理）將加倍為，計算，若還有，再將加倍為，類似做下去直到某個，則選定步長；3. 若，（說明精度不夠，做步長縮小處理，以提高精度）將折半為，計算，若還有，再將折半為，類似做下去知道某個，則選定步長。6.8 一階微