微分方程數(shù)值解

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

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

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

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

5、方法具有p階精度或該方法是p階方法。方法的階越高，方法越好。l 主局部截?cái)嗾`差或局部截?cái)嗾`差的主項(xiàng)如果某方法是p階方法，按可展為則稱為主局部截?cái)嗾`差。在同階方法中，主局部截?cái)嗾`差越小，方法越好。例6.1常微分方程初值問題的單步法為試求其局部截?cái)嗾`差主項(xiàng)并回答它是幾階精度的？解 該單步公式的局部截?cái)嗾`差是故局部截?cái)嗾`差主項(xiàng)是，方法是一階的。l 求階p的另一方法因?yàn)槿サ粝聵?biāo)，有若將在x點(diǎn)展開有則知該方法的階是p。例如，對Euler方法，有那么將在x點(diǎn)展開，有故有因此，Euler方法是一階方法。定義6.3 設(shè)用某種數(shù)值方法求初值問題（6.1）在任意節(jié)點(diǎn)的數(shù)值解時(shí)，滿足，則稱該數(shù)值方法是絕對穩(wěn)定的。這

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

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

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

9、取如下兩次計(jì)算格式由來計(jì)算（改進(jìn)的Euler方法）：預(yù)測 （6.9）校正 （6.10）預(yù)測校正公式。易證它也是二階方法。6.5 Runge-Kutta方法理論上，由Taylor展開法可以構(gòu)造出解初值問題（6.1）的高階數(shù)值方法，但這涉及到要計(jì)算的高階導(dǎo)數(shù)，因此很不方便。本節(jié)的Runge-Kutta方法從函數(shù)本身著手，避開計(jì)算的高階導(dǎo)數(shù)，來構(gòu)造高階數(shù)值方法。基本思想將微分方程初值問題轉(zhuǎn)化為積分方程問題，再對積分方程中的定積分使用待定的m點(diǎn)插值型求積公式構(gòu)造高階的函數(shù)展開模式以獲得高階數(shù)值方法。1. 構(gòu)造原理將初值問題的方程在上積分對定積分用m個(gè)點(diǎn)的插值型求積公式取，并對中的做展開，得解初值問題的

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

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

12、unge-Kutta方法計(jì)算格式（）為它們的初值，計(jì)算結(jié)果及準(zhǔn)確解列于下表Euler方法改進(jìn)Euler方法經(jīng)典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）分析求解公式的局部截?cái)嗾`差，指出它是幾階公式；2）證明用上面求解公式計(jì)算初值問題的數(shù)值解

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

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

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

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