第四章常微分方程數(shù)值解法

第四章常微分方程數(shù)值解法1第一頁，共六十五頁，2022年，8月28日常微分方程數(shù)值解法引言（常微分方程數(shù)值解法概述）顯式歐拉法、隱式歐拉法、二步歐拉法局部截斷誤差與精度改進的歐拉方法龍格-庫塔方法收斂性與穩(wěn)定性簡述一階常微分方程組與高階常微分方程2第二頁，共六十五頁，2022年，8月28日引言一階常微分方程初值問題：微分方程初始條件定理：若f(x,y)在某閉區(qū)域R：

上連續(xù)，且在R域內(nèi)滿足李普希茲(Lipschitz)條件，即存在正數(shù)L，使得對于R域內(nèi)的任意兩值y1,y2，下列不等式成立：

則上述初值問題的連續(xù)可微的解y(x)存在并且唯一。3第三頁，共六十五頁，2022年，8月28日引言（續(xù)）實際生產(chǎn)與科研中，除少數(shù)簡單情況能獲得初值問題的初等解（用初等函數(shù)表示的解）外，絕大多數(shù)情況下是求不出初等解的。有些初值問題即便有初等解，也往往由于形式過于復(fù)雜而不便處理。實用的方法是在計算機上進行數(shù)值求解：即不直接求y(x)的顯式解，而是在解所存在的區(qū)間上，求得一系列點xn(n

=0,1,2,…)上解的近似值。4第四頁，共六十五頁，2022年，8月28日歐拉（Euler）方法方法一化導(dǎo)數(shù)為差商的方法由于在逐步求解的過程中，y(xn)的準確值無法求解出來，因此用其近似值代替。為避免混淆，以下學(xué)習(xí)簡記：y(xn)：待求函數(shù)y(x)在xn處的精確函數(shù)值yn：待求函數(shù)y(x)在xn處的近似函數(shù)值5第五頁，共六十五頁，2022年，8月28日代入初值問題表達式可得：根據(jù)y0

可以一步步計算出函數(shù)y

=

y(x)在x1,x2,x3

x4,…上的近似值y1,y2,y3,y4,

…常微分方程數(shù)值解是一組離散的函數(shù)值數(shù)據(jù)，它的精確表達式很難求解得到，但可以進行插值計算后用插值函數(shù)逼近y(x)6第六頁，共六十五頁，2022年，8月28日歐拉方法（續(xù)）方法二數(shù)值積分法同樣以近似值

yn

代替精確值y(xn)可得：將微分方程y

=

f

(x,y)在區(qū)間[xn,xn+1]上積分：7第七頁，共六十五頁，2022年，8月28日歐拉方法的幾何意義xy08第八頁，共六十五頁，2022年，8月28日隱式歐拉法在數(shù)值積分法推導(dǎo)中，積分的近似值取為積分區(qū)間寬度與右端點處的函數(shù)值乘積，即：這樣便得到了隱式歐拉法：含有未知的函數(shù)值隱式歐拉法沒有顯式歐拉法方便9第九頁，共六十五頁，2022年，8月28日二步歐拉法在數(shù)值積分法推導(dǎo)中，積分區(qū)間寬度選為兩步步長，即積分區(qū)間為：[xn-1,xn+1]，則：以y(x)在xn-1,xn上的近似值代替精確值可得：需要前兩步的計算結(jié)果中矩形公式10第十頁，共六十五頁，2022年，8月28日梯形公式歐拉法在數(shù)值積分法中，如果用梯形公式近似計算f(x,y)在區(qū)間[xn,xn+1]上的積分，即：用近似值代替精確值可得梯形公式歐拉法：上式右端出現(xiàn)了未知項，可見梯形法是隱式歐拉法的一種；實際上，梯形公式歐拉法是顯式歐拉法與隱式歐拉法的算術(shù)平均。11第十一頁，共六十五頁，2022年，8月28日例用顯式歐拉法、隱式歐拉法、梯形法求解初值問題：取h=0.1，計算到x=0.5，并與精確解進行比較解：由已知條件可得：h=0.1，x0

=0，y0

=1，

f(x,y)=-y

+x

+1顯式歐拉法：12第十二頁，共六十五頁，2022年，8月28日例：（續(xù)）隱式歐拉法：化簡得：梯形公式歐拉法：13第十三頁，共六十五頁，2022年，8月28日計算結(jié)果：xn顯式法yn隱式法yn梯形法yn精確解y(xn)0.011110.11.0000001.0090911.0047621.0048370.21.0100001.0264461.0185941.0197310.31.0290001.0513151.0406331.0408180.41.0561001.0830141.0700971.0703200.51.0904901.1209221.1062781.106531本題的精確解為：顯式法隱式法梯形法14第十四頁，共六十五頁，2022年，8月28日局部截斷誤差為了簡化分析某常微分方程數(shù)值算法的誤差，現(xiàn)假設(shè)yn

=

y(xn)，即在前一步y(tǒng)n準確的前提下，估計：稱上述誤差Tn+1

為該常微分方程數(shù)值算法的局部截斷誤差如果某個常微分方程數(shù)值算法的局部截斷誤差可表示為

O(h

p+1)，則稱該數(shù)值算法的精度是p

階歐拉法的精度為一階；二步歐拉法的精度為二階；梯形法的精度為二階。15第十五頁，共六十五頁，2022年，8月28日泰勒展開法如果初值問題中的f

(x,y)充分可微，則可將y(xn+1)在點xn處展開：如果只保留線性項，忽略h2及以上各項，則：顯式歐拉公式16第十六頁，共六十五頁，2022年，8月28日局部截斷誤差的分析利用泰勒公式展開，比較各算法與展開式的前幾項將y(xn+1)在xn點處用泰勒公式展開：顯式歐拉法的局部截斷誤差：歐拉法1階精度17第十七頁，共六十五頁，2022年，8月28日補充：二元函數(shù)微分中值定理18第十八頁，共六十五頁，2022年，8月28日y(xn+1)在xn點處展開：隱式歐拉法：1階精度19第十九頁，共六十五頁，2022年，8月28日分別將y(xn+1),y(xn-1)在xn點處用泰勒公式展開：二步歐拉法的局部截斷誤差：二步歐拉法：2階精度20第二十頁，共六十五頁，2022年，8月28日梯形公式歐拉法：y(xn+1)在xn點處展開：2階精度21第二十一頁，共六十五頁，2022年，8月28日各種歐拉法的比較方法精度評述顯式歐拉法1最簡單，精度低隱式歐拉法1不便計算，穩(wěn)定性好二步歐拉法2需要兩步初值，且第2個初值只能由其它方法給出，可能對后面的遞推精度有影響梯形公式歐拉法2精度有所提高，但為隱式，需要迭代求解，計算量大22第二十二頁，共六十五頁，2022年，8月28日改進的歐拉法從上述例子可以看到，梯形法由于具有二階精度，其局部截斷誤差比顯式歐拉法和隱式歐拉法小，但梯形法實質(zhì)上是一種隱式算法顯式歐拉法是一個顯式算法，雖然計算量較小，但是精度不高綜合兩種方法的長處，可以先用顯式歐拉法求出y(xn+1)的一個粗略近似值，然后用它代入梯形法公式的右端，用梯形法計算y(xn+1)的較為精確的近似值。23第二十三頁，共六十五頁，2022年，8月28日改進的歐拉法（續(xù)）按照上述思想，可以建立如下預(yù)報-校正系統(tǒng)：按以上兩式求解常微分方程的算法稱為改進的歐拉法，它還可以表示為：嵌套形式平均化形式2階精度24第二十四頁，共六十五頁，2022年，8月28日用改進歐拉法求上例所述的初值問題并與歐拉法和梯形法比較誤差的大小。解：采用改進歐拉法的嵌套形式：25第二十五頁，共六十五頁，2022年，8月28日計算結(jié)果xn改進歐拉法yn精確解

y(xn)誤差改進歐拉法歐拉法梯形法0.11.0050001.0048371.610-44.810-37.510-50.21.0192051.0197312.910-48.710-31.410-40.31.0412181.0408184.010-41.210-21.910-40.41.0708021.0703204.810-41.410-22.210-40.51.1070761.1065315.510-41.610-22.510-4可見，改進歐拉法的誤差數(shù)量級與梯形法大致相同，而比歐拉法小得多。

26第二十六頁，共六十五頁，2022年，8月28日改進的歐拉法的意義改進的歐拉法的平均化形式y(tǒng)(xn+1)在點xn處的一階展開式為：27第二十七頁，共六十五頁，2022年，8月28日改進的歐拉法的幾何意義P0xyRhQ斜率:k2S斜率:k128第二十八頁，共六十五頁，2022年，8月28日龍格-庫塔(Runge-Kutta)方法改進的歐拉法（2階精度）y(xn+1)在點xn處的一階泰勒展開式為：顯式歐拉法（1階精度）29第二十九頁，共六十五頁，2022年，8月28日龍格-庫塔方法（續(xù)）顯式歐拉法用一個點的值k1作為k*的近似值改進的歐拉公式用二個點的值k1和k2的平均值作為k*近似值；改進的歐拉法比顯式歐拉法精度高；在[xn,xn+1]內(nèi)多預(yù)報幾個點的ki

值，并用其加權(quán)平均值作為k*的近似值從而構(gòu)造出具有更高精度的計算公式，這就是龍格-庫塔方法的基本思想。30第三十頁，共六十五頁，2022年，8月28日二階龍格-庫塔方法以k1和k2的加權(quán)平均來近似取代k*為分析局部截斷誤差，令yn=y(xn)，由泰勒公式得：31第三十一頁，共六十五頁，2022年，8月28日補充：二元泰勒展開式32第三十二頁，共六十五頁，2022年，8月28日用二元泰勒公式展開將k1,k2代入

中可得：33第三十三頁，共六十五頁，2022年，8月28日二階龍格-庫塔方法（續(xù)）2階精度34第三十四頁，共六十五頁，2022年，8月28日四個未知變量，只有三個方程，有無窮多組解每組解的構(gòu)成的龍格-庫塔方法均為二階二階龍格-庫塔方法即為改進的歐拉方法變形的歐拉法中點方法35第三十五頁，共六十五頁，2022年，8月28日三階龍格-庫塔方法三階龍格-庫塔方法是用三個值k1,k2,k3的加權(quán)平均來近似取代k*要使三階龍格-庫塔方法具有三階精度，必須使其局部截斷誤差為O(h4)將k1,k2,k3代入yn+1的表達式中，在(xn,

yn)處用二元泰勒公式展開，與y(xn+1)在xn處的泰勒展開式比較36第三十六頁，共六十五頁，2022年，8月28日三階龍格-庫塔方法（續(xù)）類似二階龍格-庫塔方法的推導(dǎo)過程，8個待定系數(shù)c1,c2,c3,a2,a3,b21,b31,b32應(yīng)滿足：8個未知參數(shù)，6個方程，有無窮多組解庫塔公式37第三十七頁，共六十五頁，2022年，8月28日四階龍格-庫塔方法類似可以推出四階龍格-庫塔公式，常用的有：標(biāo)準四階龍格-庫塔公式38第三十八頁，共六十五頁，2022年，8月28日四階龍格-庫塔方法（續(xù)）吉爾（Gill）公式4階以上龍格-庫塔公式的計算量太大，并且精度不一定提高，有時反而會降低，因此實際應(yīng)用中一般選用四階龍格-庫塔已足可滿足精度要求。39第三十九頁，共六十五頁，2022年，8月28日用經(jīng)典四階龍格-庫塔方法求解前例的初值問題，并與改進歐拉法、梯形法在x5=0.5處比較其誤差大小解：采用經(jīng)典四階龍格-庫塔公式：40第四十頁，共六十五頁，2022年，8月28日四階R-K方法的精度比二階方法高得多精確解為：R-K方法的誤差：改進歐拉法的誤差：梯形法的誤差：41第四十一頁，共六十五頁，2022年，8月28日變步長的龍格-庫塔方法設(shè)y(xn)

在xn處的值yn=

y(xn)，當(dāng)xn+1

=

xn+h時 y(xn+1)的近似值為，由于四階R-K方法的精度為4階，故局部截斷誤差為：用四階R-K方法求解初值問題精度較高，但要從理論上給出誤差|y

(xn)-

yn|的估計式則比較困難；那么應(yīng)如何判斷計算結(jié)果的精度以及如何選擇合適的步長h？通常是通過不同步長在計算機上的計算結(jié)果進行近似估計。42第四十二頁，共六十五頁，2022年，8月28日若以h/2為步長，從xn出發(fā)，經(jīng)過兩步計算，得到

y(xn+1)

的近似值變步長的龍格-庫塔方法（續(xù)）以上每步的截斷誤差約為cn(h/2)5，于是兩步的局部截斷誤差為：于是：整理得：43第四十三頁，共六十五頁，2022年，8月28日變步長的龍格-庫塔方法（續(xù)）記：，給定的精度要求為eD>e，反復(fù)將步長折半計算，直至D<e，取最終得到的作為y(xn+1)的近似值。D<e，反復(fù)將步長加倍計算，直至D>e，再將步長折半一次計算，最終得到符合精度要求的y(xn+1)的近似值。44第四十四頁，共六十五頁，2022年，8月28日單步法的收斂性顯式單步法可統(tǒng)一寫成：增量函數(shù)，僅依賴于函數(shù)f，且僅僅是xn,yn,h的函數(shù)求y

=

y(x)求y(xn)，xn=x0+

nh離散化求yny(xn)某種數(shù)值方法h

0時，近似解是否收斂到精確解 ，它應(yīng)當(dāng)是一個固定節(jié)點，因此h

0時應(yīng)同時附帶n

45第四十五頁，共六十五頁，2022年，8月28日單步法的收斂性（續(xù)）對于p

階的常微分方程數(shù)值算法，當(dāng)h

0，n時，是否yn+1y(xn+1)？p

階算法的局部截斷誤差為：顯然：局部截斷誤差的前提假設(shè)是：局部截斷誤差0并不能保證算法收斂46第四十六頁，共六十五頁，2022年，8月28日單步法的收斂性（續(xù)）定義：若求解某初值問題的單步數(shù)值法，對于固定的當(dāng)h

0且n

時，它的近似 解趨向于精確解y(xn)，即：則稱該單步法是收斂的。定義：稱y(xn)-

yn

為單步法的近似解yn

的整體截斷

誤差。單步法收斂h

0且n

時，yn

的整體截斷誤差

047第四十七頁，共六十五頁，2022年，8月28日單步法的收斂性（續(xù)）收斂性定理若某單步法滿足以上條件，則該方法是收斂的則該單步法的整體截斷誤差為：若單步法具有p

階精度，且增量函數(shù)關(guān)于y

滿足：Lipschitz條件：初值y0

是準確的48第四十八頁，共六十五頁，2022年，8月28日假設(shè)在前一步

yn準確的前提下求得的近似值為：算法精度為

p階，局部截斷誤差：49第四十九頁，共六十五頁，2022年，8月28日50第五十頁，共六十五頁，2022年，8月28日即：若初值是準確的，則e0=0

，從而整體截斷誤差為：y

=e

x為單調(diào)增函數(shù)，當(dāng)時當(dāng)h

0且n

時，則51第五十一頁，共六十五頁，2022年，8月28日單步法的穩(wěn)定性在討論單步法收斂性時一般認為數(shù)值方法本身的計算過程是準確的，實際上并非如此：初始值y0有誤差d=

y0

-y(x0)后續(xù)的每一步計算均有舍入誤差這些初始和舍入誤差在計算過程的傳播中是逐步衰減的還是惡性增長就是數(shù)值方法的穩(wěn)定性問題52第五十二頁，共六十五頁，2022年，8月28日定義：若一種數(shù)值方法在節(jié)點xn

處的數(shù)值解

yn

的擾動，而在以后各節(jié)點ym

(m>n)上產(chǎn)生的擾動為，如果：單步法的穩(wěn)定性（續(xù)）定義：設(shè)在節(jié)點xn

處用數(shù)值算法得到的理想數(shù)值解為yn，而實際計算得到的近似解為，稱差值：為第n步的數(shù)值解的擾動。則稱該數(shù)值方法是穩(wěn)定的。53第五十三頁，共六十五頁，2022年，8月28日單步法的穩(wěn)定性（續(xù)）歐拉法：由于函數(shù)f

(x,y)的多樣性，數(shù)值穩(wěn)定性的分析相當(dāng)復(fù)雜，通常只研究模型方程考