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

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

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

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

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

滿足Lipschitz

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

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

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

上一頁下一頁返回

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

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

為Euler方法的局部截斷誤差.上一頁下一頁返回

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

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

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

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

可見Euler方法的計算結果精度不太高。

上一頁下一頁返回

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

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

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

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

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

i

步計算是精確的前提下，考慮的截斷誤差Ti+1=y(xi+1)

yi+1稱為局部截斷誤差。定義若某算法的局部截斷誤差為O(hp+1)，則稱該算法有p

階精度。顯然,p越大,精度越高.三、局部截斷誤差與方法的階（將準確解代入公式的左、右兩端，其左端與右端之差）

Euler方法的精度

其中：

上一頁下一頁返回

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

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

上一頁下一頁返回

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

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

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

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

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

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

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

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

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

上一頁下一頁返回

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

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

17第十七頁，共三十三頁，2022年，8月28日計算結果見下表Euler方法改進的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é)

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

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

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

上一頁下一頁返回

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

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

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

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

23第二十三頁，共三十三頁，2022年，8月28日例用四階標準R-K公式解初值問題

解：取h=0.2，四階標準R-K法的具體格式如下:上一頁下一頁返回

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

上一頁下一頁返回

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

。

上一頁下一頁返回

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

而4階R-K法的效果要比改進的Euler法好.

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

取小。上一頁下一頁返回

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

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

/*Convergency*/定義若某算法對于任意固定的x=xi=x0+ih，當h0

(同時i)時有yi

y(xi

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

，有上一頁下一頁返回

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

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

節(jié)點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日定義若某算法在計算過程中任一步產(chǎn)生的誤差在以后的計算中都逐步衰減，則稱該算法是絕對穩(wěn)定的.一般分析時為簡單起見，只考慮試驗方程常數(shù)l<0，可以是復數(shù)當步長取為h

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

絕對穩(wěn)定，z

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

比算法B

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