常微分方程的數(shù)值解法

常微分方程的數(shù)值解法1第一頁，共三十三頁，2022年，8月28日上一頁下一頁返回

本章介紹求解微分方程數(shù)值解的基本思想和方法.

含有自變量、未知函數(shù)和它的一階導(dǎo)數(shù)和高階導(dǎo)數(shù)的方程.常微分方程它是描述運(yùn)動(dòng)、變化規(guī)律的重要數(shù)學(xué)方法之一，分為兩類：1.初值問題即給出未知函數(shù)及導(dǎo)數(shù)在初始點(diǎn)的值；2.邊值問題即給出未知函數(shù)及（或）它的某些導(dǎo)數(shù)在區(qū)間兩個(gè)端點(diǎn)的值。2第二頁，共三十三頁，2022年，8月28日

考慮一階常微分方程的初值問題:只要f(x,y)在[a,b]R1上連續(xù)，且關(guān)于y

滿足Lipschitz

條件，即存在與x,y無關(guān)的常數(shù)L

使對(duì)任意定義在[a,b]上的y1(x)和y2(x)都成立，則上述問題解存在唯一解。所謂數(shù)值解法就是要計(jì)算出初值問題的解函數(shù)y(x)在一系列離散點(diǎn)a=x0<x1<…<xN=b上的近似值：y0,y1,……yN.節(jié)點(diǎn)間距為步長(zhǎng)，通常采用等距節(jié)點(diǎn)，即取hi=h

(常數(shù))。{yn}稱為問題的數(shù)值解.數(shù)值解所滿足的離散方程統(tǒng)稱為差分格式.

上一頁下一頁返回

3第三頁，共三十三頁，2022年，8月28日第一節(jié)

歐拉方法一、歐拉公式令yn為y(xn)的近似值，將上式代入(*)式可得此式稱為歐拉(Euler)公式.

為Euler方法的局部截?cái)嗾`差.上一頁下一頁返回

4第四頁，共三十三頁，2022年，8月28日例1

用歐拉公式解初值問題解：取步長(zhǎng)h=0.1,歐拉公式的具體形式為：依次計(jì)算可得………

y3，y4，y5，y6，y7，y8，y9，y10上一頁下一頁返回

5第五頁，共三十三頁，2022年，8月28日其部分結(jié)果見下表

可見Euler方法的計(jì)算結(jié)果精度不太高。

上一頁下一頁返回

6第六頁，共三十三頁，2022年，8月28日歐拉公式的幾何意義：x0P0x1P1x2P2xnPn幾何意義：用折線近似代替方程的解曲線,因而也稱Euler方法為折線法.上一頁下一頁返回

7第七頁，共三十三頁，2022年，8月28日二、后退的歐拉公式也用一階差商逼近導(dǎo)數(shù)令yn+1為y(xn+1)的近似值，則可得稱為后退Euler公式已知yn時(shí),必須通過解方程才能求出yn＋1

,這樣的公式稱為隱式公式,而Euler公式為顯式公式.

Euler公式和后退Euler公式都是由yn去計(jì)算yn+1，因此，稱它們?yōu)閱尾椒?。上一頁下一頁返?/p>

8第八頁，共三十三頁，2022年，8月28日定義在假設(shè)yi=y(xi)，即第

i

步計(jì)算是精確的前提下，考慮的截?cái)嗾`差Ti+1=y(xi+1)

yi+1稱為局部截?cái)嗾`差。定義若某算法的局部截?cái)嗾`差為O(hp+1)，則稱該算法有p

階精度。顯然,p越大,精度越高.三、局部截?cái)嗾`差與方法的階（將準(zhǔn)確解代入公式的左、右兩端，其左端與右端之差）

Euler方法的精度

其中：

上一頁下一頁返回

9第九頁，共三十三頁，2022年，8月28日所以，Euler方法具有1階精度。將在點(diǎn)處一階Taylor展開上一頁下一頁返回

10第十頁，共三十三頁，2022年，8月28日所以，后退的Euler方法也具有1階精度。將在點(diǎn)處一階Taylor展開隱式Euler方法的精度

上一頁下一頁返回

11第十一頁，共三十三頁，2022年，8月28日—顯、隱式兩種算法的平均

歐拉公式的改進(jìn)其局部誤差為：此公式具有2階精度.稱平均公式或梯形公式梯形公式可由下迭代式計(jì)算：其中迭代初值是Euler公式提供.上一頁下一頁返回

12第十二頁，共三十三頁，2022年，8月28日四、改進(jìn)的歐拉公式Step1:

先用顯式歐拉公式作預(yù)測(cè)，算出),(1iiiiyxfhyy+=+Step2:再將代入隱式梯形公式的右邊作校正，得到1+iy)],(),([2111+++++=iiiiiiyxfyxfhyy注：此法亦稱為預(yù)測(cè)-校正法?？梢宰C明該算法具有2階精度，同時(shí)可以看到它是個(gè)單步遞推格式，比隱式公式的迭代求解過程簡(jiǎn)單。它的精度高于顯式歐拉法。上一頁下一頁返回

13第十三頁，共三十三頁，2022年，8月28日為了便于編程,常將改進(jìn)的歐拉公式寫為：上一頁下一頁返回

14第十四頁，共三十三頁，2022年，8月28日例2用改進(jìn)的歐拉法解例1中的初值問題.解：取步長(zhǎng)h=0.1,

改進(jìn)歐拉法的具體形式為具體計(jì)算過程如下上一頁下一頁返回

15第十五頁，共三十三頁，2022年，8月28日xn改進(jìn)的歐拉法誤差xn改進(jìn)的歐拉法誤差0100.61.4859560.0027160.21.1840960.0000880.81.6164760.0040240.41.3433600.0017191.01.7378690.005818依次計(jì)算可得………

y3，y4，y5，y6，y7，y8，y9，y10其部分結(jié)果見下表

上一頁下一頁返回

16第十六頁，共三十三頁，2022年，8月28日例3對(duì)下面的初值問題解

（1）取步長(zhǎng)h=0.1,歐拉方法的具體公式為（2）取步長(zhǎng)h=0.1,改進(jìn)的歐拉方法的具體公式為取步長(zhǎng)h=0.1，分別用Euler方法、改進(jìn)的Euler方法求數(shù)值解。上一頁下一頁返回

17第十七頁，共三十三頁，2022年，8月28日計(jì)算結(jié)果見下表Euler方法改進(jìn)的Euler方法xnynyn0.10.20.30.40.50.60.70.80.91.00.9000000.8100000.7290000.6561000.5904900.5314410.4782970.4304670.3874210.3486790.9050000.8190250.7412180.6708020.6070760.5494040.4972100.4499750.4072280.368541上一頁下一頁返回

18第十八頁，共三十三頁，2022年，8月28日第二節(jié)

龍格-庫塔法基本思想考察改進(jìn)的歐拉法，可以將其改寫為：斜率一定取k1k2的平均值嗎？步長(zhǎng)一定是一個(gè)h

嗎？只要能對(duì)平均斜率提供一種近似算法,就能得到一種對(duì)應(yīng)的差分格式.上一頁下一頁返回

19第十九頁，共三十三頁，2022年，8月28日例如取m個(gè)點(diǎn)的斜率構(gòu)造如下形式的公式該公式稱為m級(jí)龍格－庫塔(Runge-Kutta)公式,簡(jiǎn)稱R-K公式.求解：只需將公式的局部截?cái)嗾`差在xn點(diǎn)進(jìn)行Taylor展開,令其前面盡可能多的項(xiàng)為0,便可導(dǎo)出ai，bij，ci所滿足的方程組,即可從中求出這些系數(shù).

上一頁下一頁返回

20第二十頁，共三十三頁，2022年，8月28日以m=2的情形為例說明建立R-K公式的方法.其局部截?cái)嗾`差為：上一頁下一頁返回

21第二十一頁，共三十三頁，2022年，8月28日因此有：而對(duì)于h3，若將k2的Taylor展開式多取一項(xiàng),會(huì)發(fā)現(xiàn)h3項(xiàng)的系數(shù)不可能為0.而對(duì)于上式有無窮多個(gè)解,它的每一組解都給出了一個(gè)局部截?cái)嗾`差為的二級(jí)R-K公式,即二階R-K公式.當(dāng)取時(shí),二階R-K公式就是改進(jìn)的Euler公式

這里有個(gè)未知數(shù)，個(gè)方程。32上一頁下一頁返回

22第二十二頁，共三十三頁，2022年，8月28日常用的標(biāo)準(zhǔn)四階R－K公式（經(jīng)典R-K方法）最常用的四階標(biāo)準(zhǔn)R－K公式（經(jīng)典R-K方法）為：上一頁下一頁返回

23第二十三頁，共三十三頁，2022年，8月28日例用四階標(biāo)準(zhǔn)R-K公式解初值問題

解：取h=0.2，四階標(biāo)準(zhǔn)R-K法的具體格式如下:上一頁下一頁返回

24第二十四頁，共三十三頁，2022年，8月28日已知

上一頁下一頁返回

25第二十五頁，共三十三頁，2022年，8月28日同理可計(jì)算得具體結(jié)果見下表至少具有四位有效數(shù)字.比較:上節(jié)用改進(jìn)的Euler公式計(jì)算,取h=0.1，最多具有四位有效數(shù)字

。

上一頁下一頁返回

26第二十六頁，共三十三頁，2022年，8月28日改進(jìn)的Euler公式每前進(jìn)一步只要計(jì)算兩次f值,而4階R-K公式每前進(jìn)一步要計(jì)算四次f值,但改進(jìn)的Euler法的步長(zhǎng)比4階R-K法的小一半,兩者計(jì)算總量差不多.

而4階R-K法的效果要比改進(jìn)的Euler法好.

由于龍格-庫塔法的導(dǎo)出基于泰勒展開，故精度主要受解函數(shù)的光滑性影響。對(duì)于光滑性不太好的解，最好采用低階算法而將步長(zhǎng)h

取小。上一頁下一頁返回

27第二十七頁，共三十三頁，2022年，8月28日第三節(jié)

單步法的收斂性與穩(wěn)定性收斂性

/*Convergency*/定義若某算法對(duì)于任意固定的x=xi=x0+ih，當(dāng)h0

(同時(shí)i)時(shí)有yi

y(xi

)，則稱該算法是收斂的。例：就初值問題考察歐拉顯式格式的收斂性。解：該問題的精確解為歐拉公式為對(duì)任意固定的x=xi=ih

，有上一頁下一頁返回

28第二十八頁，共三十三頁，2022年，8月28日穩(wěn)定性

/*Stability*/例：考察初值問題在區(qū)間[0,0.5]上的解。分別用歐拉法、隱式歐拉法和改進(jìn)的歐拉格式計(jì)算數(shù)值解。0.00.10.20.30.40.5精確解改進(jìn)歐拉法隱式歐拉法歐拉法

節(jié)點(diǎn)xi

1.00002.00004.00008.00001.6000101

3.2000101

1.00002.5000101

6.25001021.56251023.90631039.76561041.00002.50006.25001.56261013.90631019.76561011.00004.97871022.47881031.23411046.14421063.0590107Whatiswrong??!上一頁下一頁返回

29第二十九頁，共三十三頁，2022年，8月28日定義若某算法在計(jì)算過程中任一步產(chǎn)生的誤差在以后的計(jì)算中都逐步衰減，則稱該算法是絕對(duì)穩(wěn)定的.一般分析時(shí)為簡(jiǎn)單起見，只考慮試驗(yàn)方程常數(shù)l<0，可以是復(fù)數(shù)當(dāng)步長(zhǎng)取為h

時(shí)，將某算法應(yīng)用于上式，并假設(shè)在初值產(chǎn)生誤差，則若此誤差以后逐步衰減，就稱該算法相對(duì)于z=lh

絕對(duì)穩(wěn)定，z

的全體構(gòu)成絕對(duì)穩(wěn)定區(qū)域。我們稱算法A

比算法B

穩(wěn)定，就是指A的絕對(duì)穩(wěn)定區(qū)域比B