常微分方程的數值解法

常微分方程的數值解法1第1頁，共33頁，2023年，2月20日，星期四上一頁下一頁返回

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

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

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

滿足Lipschitz

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

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

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

上一頁下一頁返回

3第3頁，共33頁，2023年，2月20日，星期四第一節(jié)

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

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

4第4頁，共33頁，2023年，2月20日，星期四例1

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

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

5第5頁，共33頁，2023年，2月20日，星期四其部分結果見下表

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

上一頁下一頁返回

6第6頁，共33頁，2023年，2月20日，星期四歐拉公式的幾何意義：x0P0x1P1x2P2xnPn幾何意義：用折線近似代替方程的解曲線,因而也稱Euler方法為折線法.上一頁下一頁返回

7第7頁，共33頁，2023年，2月20日，星期四二、后退的歐拉公式也用一階差商逼近導數令yn+1為y(xn+1)的近似值，則可得稱為后退Euler公式已知yn時,必須通過解方程才能求出yn＋1

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

Euler公式和后退Euler公式都是由yn去計算yn+1，因此，稱它們?yōu)閱尾椒āＩ弦豁撓乱豁摲祷?/p>

8第8頁，共33頁，2023年，2月20日，星期四定義在假設yi=y(xi)，即第

i

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

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

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

Euler方法的精度

其中：

上一頁下一頁返回

9第9頁，共33頁，2023年，2月20日，星期四所以，Euler方法具有1階精度。將在點處一階Taylor展開上一頁下一頁返回

10第10頁，共33頁，2023年，2月20日，星期四所以，后退的Euler方法也具有1階精度。將在點處一階Taylor展開隱式Euler方法的精度

上一頁下一頁返回

11第11頁，共33頁，2023年，2月20日，星期四—顯、隱式兩種算法的平均

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

12第12頁，共33頁，2023年，2月20日，星期四四、改進的歐拉公式Step1:

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

13第13頁，共33頁，2023年，2月20日，星期四為了便于編程,常將改進的歐拉公式寫為：上一頁下一頁返回

14第14頁，共33頁，2023年，2月20日，星期四例2用改進的歐拉法解例1中的初值問題.解：取步長h=0.1,

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

15第15頁，共33頁，2023年，2月20日，星期四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第16頁，共33頁，2023年，2月20日，星期四例3對下面的初值問題解

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

17第17頁，共33頁，2023年，2月20日，星期四計算結果見下表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第18頁，共33頁，2023年，2月20日，星期四第二節(jié)

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

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

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

上一頁下一頁返回

20第20頁，共33頁，2023年，2月20日，星期四以m=2的情形為例說明建立R-K公式的方法.其局部截斷誤差為：上一頁下一頁返回

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

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

22第22頁，共33頁，2023年，2月20日，星期四常用的標準四階R－K公式（經典R-K方法）最常用的四階標準R－K公式（經典R-K方法）為：上一頁下一頁返回

23第23頁，共33頁，2023年，2月20日，星期四例用四階標準R-K公式解初值問題

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

24第24頁，共33頁，2023年，2月20日，星期四已知

上一頁下一頁返回

25第25頁，共33頁，2023年，2月20日，星期四同理可計算得具體結果見下表至少具有四位有效數字.比較:上節(jié)用改進的Euler公式計算,取h=0.1，最多具有四位有效數字

。

上一頁下一頁返回

26第26頁，共33頁，2023年，2月20日，星期四改進的Euler公式每前進一步只要計算兩次f值,而4階R-K公式每前進一步要計算四次f值,但改進的Euler法的步長比4階R-K法的小一半,兩者計算總量差不多.

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

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

取小。上一頁下一頁返回

27第27頁，共33頁，2023年，2月20日，星期四第三節(jié)

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

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

(同時i)時有yi

y(xi

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

，有上一頁下一頁返回

28第28頁，共33頁，2023年，2月20日，星期四穩(wěn)定性

/*Stability*/例：考察初值問題在區(qū)間[0,0.5]上的解。分別用歐拉法、隱式歐拉法和改進的歐拉格式計算數值解。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第29頁，共33頁，2023年，2月20日，星期四定義若某算法在計算過程中任一步產生的誤差在以后的計算中都逐步衰減，則稱該算法是絕對穩(wěn)定的.一般分析時為簡單起見，只考慮試驗方程常數l<0，可以是復數當步長取為h

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

絕對穩(wěn)定，z

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

比算法B

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

的