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

1、第六章常微分方程的數(shù)值解法6.0 引言6.1 算法構造的主要途徑6.2 Runge-Kutta Method算法6.3 線性多步法6.4 線性多步法的一般形式6.5 算法的穩(wěn)定性、收斂性6.0 引言1主要考慮如下的一階常微分方程初值問題的求解:微分方程的解就是求一個函數(shù)y=y(x)，使得該函數(shù)滿足微分方程并且符合初值條件。2. 例如微分方程:xy-2y=4x ；初始條件: y(1)=-3。于是可得一階常微分方程的初始問題。顯然函數(shù)y(x)=x2-4x滿足以上條件,因而是該初始問題的微分方程的解。3. 但是,只有一些特殊類型的微分方程問題能夠得到用解析表達式表示的函數(shù)解，而大量的微分方程問題很難

2、得到其解析解，有的甚至無法用解析表達式來表示。因此，只能依賴于數(shù)值方法去獲得微分方程的數(shù)值解。4微分方程的數(shù)值解：設微分方程問題的解y(x)的存在區(qū)間是a,b，初始點x0=a，將a,b進行劃分得一系列節(jié)點x0 , x1,.,xn，其中a= x0 x1 xn =b。y(x)的解析表達式不容易得到或根本無法得到，我們用數(shù)值方法求得y(x)在每個節(jié)點xk的近似值y(xk)，即yy(xk)，這樣y0, y1,.,yn稱為微分方程的數(shù)值解。如圖所示:ab x0 x1 x2 . xn-1 xn6.1 算法構造的主要途徑1 歐拉公式1.1 構造的思想：微分方程初值問題: 利用差商代替一階導數(shù)，即，則。于是，

3、可求出y(x1)的近似值y1，同樣地，可利用x1處的微分方程可得：一般地，利用在xn處的微分方程可得：此式稱為歐拉公式。1.2 幾何意義：對于微分方程y=2(x+1),其通解是y=(x+1)2+c,是一個曲線族，當給定初值條件y(0)=2，其特解為y=(x+1)2+1。如圖所示：由y(0)=2，過該曲線上一點(0,2)作曲線的切線，其斜率，切線為：,因此可計算出y1,如此，可根據(jù)：，故歐拉法又稱歐拉折線法。1.3 算例：例：解：h=0.2 , xi=1+ih y1= y0+hf(x0，y0)=-1+0.2y2= y1+hf(x1，y1)=-1+0.2y3= y2+hf(x2，y2)= -0.9

4、333+0.2y4= y3+hf(x3，y3)=0.8+0.2y5= y4+hf(x4，y4)=0.6+0.2y6= y5+hf(x5，y5)=0.3333+0.2精確解為：y=x2-2xxky(xk)ykek1.2-0.96-10.041.4-0.84-0.93330.09331.6-0.64-0.80.161.8-0.36-0.60.242.00-0.33330.33332.20.4400.44可以看出誤差隨著計算在積累。1.4 Euler法的特點和誤差迭代格式：特點：(1)單步方法；(2)顯式格式；(3)局部截斷誤差為。局部截斷誤差：當時，由按照歐拉方法計算來的的誤差稱為局部截斷誤差。即

5、，是局部截斷誤差。如：歐拉法得：因此，局部截斷誤差是。2改進Euler法2.1方法構造在微分方程初值問題,對其從到進行定積分得:將右端的定積分用梯形公式來進行近似計算。用和來分別代替和得計算格式:這就是改進歐拉方法。2.2 顯式格式和隱式格式在歐拉式中每一步計算已知,直接用格式可以計算出,此類格式稱為顯式格式。而在改進歐拉方法中在每一步計算中是未知,待求的,未知量在中這是一個方程,如是非線形或超越函數(shù),此方程是無法直接解出來(要依靠迭代法才能解出)。這類格式稱為隱式格式。2.3 算例例: 用改進歐拉方法求解。解：解得:注意：由于,是線形函數(shù)可以從隱式格式中解出。問題的精確解是歐拉方法誤差3預測

6、校正方法由于改進Euler法是隱式格式，無法從格式中直接求出必須要解方程。下面用預測校正方法來求隱式格式中的。預測值:校正值:此式相當于對隱式格式求時采用迭代的方法,用歐拉格式得到的作為初始值迭代公式迭代一次而已,此公式代入后得:如改寫成平均的形式為:6.2 龍格-庫塔法Runge-Kutta Method1 龍格-庫塔法的思想1.1考慮微分方程的初值問題:根據(jù)微分中值定理有：，其中01。于是即我們稱為y(x)在區(qū)間xk, xk+1上的平均斜率，記作K，其中，是存在但是未知的。因此，如何對平均斜率K進行近似計算，相應地就得到一種近似公式，或稱為微分方程的一種計算格式。1.2 例如：用f(xk,

7、yk)作為平均斜率K的近似值就得到歐拉格式：。用作為平均斜率K的近似值就得到比歐拉格式高一階精度的格式，即，改進歐拉格式的預測-校正方法:1.3 啟發(fā)(Motivative)能否在二維平面中xxk, xk+1, yyk, yk+1上多找一些f(x,y)點，在這些點上作函數(shù)值的平均數(shù)并以此作為平均斜率K的近似值。由于自由度的增加，使得的p能夠提高，從而達到提高精度的目標，這就是龍格-庫塔法的基本思想。1.4 R-K的一般形式在二維平面區(qū)域上取N個點，得到公式:，其中Kj是二元函數(shù)f(x,y)在這N個點上的值。其中是待定常數(shù)。2 二階龍格-庫塔法2.1計算格式: 應用Taylor公式求參數(shù)：。將f

8、(x,y)在(xk , yk)上展開又因以上兩式代入中，又由得比較兩邊得四個未知量,三個方程,有無窮多組解，但局部截斷誤差均為,是二階精度的。例如:(1) 取,則,就是改進歐拉公式的預測-校正方法(或二階龍格-庫塔法)。(2) 取，則，則,得二階中點法：3 三階龍格庫塔方法3.1 計算格式: 參數(shù)有,共有八個參數(shù)。將f和均在點Tayor展開,有利用,按展開式代入(高于的不計)得:以上展開式代入中比較得:(1)(2)(3)(4)(5)(6)(7)(8)上述八個方程有六個獨立方程。八個未知量,六個方程有無窮多組解,但其截斷誤差均為,具有三階精度。例如:，得:此方法稱為三階Kutta方法。例如:,得

9、:此方法稱為三階Heun法。4 四階龍格庫塔方法4.1 顯式四階龍格-庫塔方法的一般形式是:利用與三階龍格-庫塔類似的處理方法，經過推導，可得古典四階龍格-庫塔算法:其特點是：截斷誤差為，階數(shù)為四階精度，為顯式四階龍格庫塔方法。計算均要計算四個f的值。一般一階常微分方程初值問題的解均用四階龍格-庫塔方法進行計算，其精度滿足實際問題的精度要求。數(shù)值例子：例：初值問題用四階古典RungeKutta方法，。0.21.18322931.18321600.41.34168031.34164080.61.48328381.48323970.81.61251721.61245151.01.73214631.

10、7320508例：6.3線性多步法1 單步方法和多步方法單步方法：歐拉方法、改進歐拉方法、龍格-庫塔方法均是單步方法，即在每一步計算時，只要前面一個值已知的條件下就可以計算出。單步方法特點：(1)可以自成系統(tǒng)進行直接計算，因為初始條件只有一個已知，由可以計算，不必借助于其它方法，這種單步方法是自開始的。(2)如果格式簡單如歐拉方法，則只有一階精度，如果提高精度，則計算很復雜，如RungeKutta方法。(3)公式的構造推導也很復雜。多步方法：利用前面已知計算出來的，由已經計算好的個值來計算，這樣可以提高算法的精度，該方法稱為多步方法，利用k個值計算，稱為k步方法。多步方法的特點：(1)因初始條

11、件只有一個，運用多步方法要借助高階的單步方法來開始。例如，已知用單步的四階Runge-Kutta方法計算，再計算，再由計算，用單步方法有后運用四階的四步方法，由計算；由計算；由計算；一直下去，可以用多步方法，并且始終達到四階精度。(2)多步方法比較簡單，只要在這四個點的函數(shù)值的線性組合，而且每步中后三個函數(shù)值下一步還可使用。2顯式Adams方法：2.1 構造的思想考慮微分方程初值問題將微分方程在上積分，2.2 顯式Adams方法若已知來計算，簡記，用的拉格朗日插值多項式代替f，截斷部分，用等距步長，上面積分很簡單，得到的方法就是顯式四階Adams方法。以上可以看到該方法的局部截斷誤差是因而是四

12、階精度的。例如：解：取，首先用四階RungeKutta方法來起步，計算出，,下面不必用RungeKutta方法，而開始用四階Adams方法。(1)、求(2)、求只要補算(3)、求只要補算現(xiàn)列表看用Adams方法求出的誤差，精解為0.81.61142311.61245151.0210-31.01.72984031.73205082.2110-31.21.84066161.84390893.2410-32.3 隱式Adams方法用作為插值結點，由于也是插值結點，必帶來從而導致是隱式格式。用插值多項式來代替積分中的得：截掉得近似公式：得：，從而得四階隱式Adams方法。因，而是未知的，故這是隱式格式

13、。隱式格式的解法用預測-校正法：用顯式格式作為預測值，再用隱式格式來校正。預測值：校正值：6.4 線性多步法的一般形式微分方程初值問題: 1 k步線性法的一般結構是常數(shù)線性多步法 , k步法不同時為零 k步法 ,隱式結構;否則為顯式結構。2 k步p階線性法的推導由: 定義算子：由Taylor展開右端有代入上式：其中，因此，若有次可微，令,則因而算子對應的線性多步法為p階k步方法。3算例例：考慮2步法記于是：故一般的二步法為：同時可以算出：當時，上述公式為三階的。當時，上述公式為四階的。6.5 收斂性與穩(wěn)定性1 定義收斂性：若某算法對于任意固定的x = xn= x0 + n h，當h0 ( 同時

14、n ) 時,有yny( xn)，則稱該算法是收斂的。1.2 算例例：對于初值問題考察歐拉顯式格式的收斂性。解：該問題的精確解為歐拉公式為即：對任意固定的x = xn = n h，有由2 穩(wěn)定性例：考察初值問題在區(qū)間0, 0.5上的解。分別用歐拉顯式、隱式格式和改進的歐拉格式計算數(shù)值解。節(jié)點xi歐拉顯式歐拉隱式改進歐拉法精確解0.00.10.20.30.40.51.0000-2.00004.0000-8.0000 1.6000101-3.20001011.00002.500010-16.250010-21.562510-23.906310-39.765610-41.00002.50006.250

15、01.56261013.90631019.76561011.00004.978710-22.478810-31.234110-46.144210-63.059010-7誤差包括截斷誤差(算法理論誤差)和舍入誤差兩個部分。后者由計算機字長等決定，屬于穩(wěn)定性問題。若某算法在計算過程中任一步產生的誤差在以后的計算中都逐步衰減，則稱該算法是絕對穩(wěn)定的討論數(shù)值方法的穩(wěn)定性，通常只考慮試驗方程當步長取為h時，將某算法應用于上式，并假設只在初值產生誤差，則若此誤差以后逐步衰減，就稱該算法相對于絕對穩(wěn)定，的全體構成絕對穩(wěn)定區(qū)域。有時，我們稱算法A 比算法B 穩(wěn)定，就是指 A 的絕對穩(wěn)定區(qū)域比 B 的絕對穩(wěn)定區(qū)域大。2.1 算例：考察顯式歐拉法算法的穩(wěn)定性解：由得：于是：因此，要保證初始誤