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

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

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

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

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

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

=

y(x)在x1,x2,x3

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

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

yn

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

=

f

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

=0，y0

=1，

f(x,y)=-y

+x

+1顯式歐拉法：12第十二頁(yè)，共六十五頁(yè)，2022年，8月28日例：（續(xù)）隱式歐拉法：化簡(jiǎn)得：梯形公式歐拉法：13第十三頁(yè)，共六十五頁(yè)，2022年，8月28日計(jì)算結(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第十四頁(yè)，共六十五頁(yè)，2022年，8月28日局部截?cái)嗾`差為了簡(jiǎn)化分析某常微分方程數(shù)值算法的誤差，現(xiàn)假設(shè)yn

=

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

為該常微分方程數(shù)值算法的局部截?cái)嗾`差如果某個(gè)常微分方程數(shù)值算法的局部截?cái)嗾`差可表示為

O(h

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

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

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

y(xn)誤差改進(jìn)歐拉法歐拉法梯形法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可見(jiàn)，改進(jìn)歐拉法的誤差數(shù)量級(jí)與梯形法大致相同，而比歐拉法小得多。

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

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

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

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

在xn處的值yn=

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

=

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

(xn)-

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

y(xn+1)

的近似值變步長(zhǎng)的龍格-庫(kù)塔方法（續(xù)）以上每步的截?cái)嗾`差約為cn(h/2)5，于是兩步的局部截?cái)嗾`差為：于是：整理得：43第四十三頁(yè)，共六十五頁(yè)，2022年，8月28日變步長(zhǎng)的龍格-庫(kù)塔方法（續(xù)）記：，給定的精度要求為eD>e，反復(fù)將步長(zhǎng)折半計(jì)算，直至D<e，取最終得到的作為y(xn+1)的近似值。D<e，反復(fù)將步長(zhǎng)加倍計(jì)算，直至D>e，再將步長(zhǎng)折半一次計(jì)算，最終得到符合精度要求的y(xn+1)的近似值。44第四十四頁(yè)，共六十五頁(yè)，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時(shí)，近似解是否收斂到精確解 ，它應(yīng)當(dāng)是一個(gè)固定節(jié)點(diǎn)，因此h

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

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

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

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

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

0且n

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

yn

為單步法的近似解yn

的整體截?cái)?/p>

誤差。單步法收斂h

0且n

時(shí)，yn

的整體截?cái)嗾`差

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

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

滿足：Lipschitz條件：初值y0

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

yn準(zhǔn)確的前提下求得的近似值為：算法精度為

p階，局部截?cái)嗾`差：49第四十九頁(yè)，共六十五頁(yè)，2022年，8月28日50第五十頁(yè)，共六十五頁(yè)，2022年，8月28日即：若初值是準(zhǔn)確的，則e0=0

，從而整體截?cái)嗾`差為：y

=e

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

0且n

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

y0

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

處的數(shù)值解

yn

的擾動(dòng)，而在以后各節(jié)點(diǎn)ym

(m>n)上產(chǎn)生的擾動(dòng)為，如果：?jiǎn)尾椒ǖ姆€(wěn)定性（續(xù)）定義：設(shè)在節(jié)點(diǎn)xn

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

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