第六章常微分方程初值問題的數(shù)值解法-課件

1、常微分方程初值問題的數(shù)值解法第6章常微分方程初值問題的數(shù)值解法第6章引言在實際問題中，常需要求解微分方程(如發(fā)電機運動方程)。只有簡單的和典型的微分方程可以求出解析解，而在實際問題中的微分方程往往無法求出解析解。在高等數(shù)學(xué)中我們見過以下常微分方程：-(1)-(2)一階常微分方程引言在實際問題中，常需要求解微分方程(如發(fā)電機運動方程)。只-(3)(1),(2)式稱為初值問題,(3)式稱為邊值問題-(4)另外,在實際應(yīng)用中還經(jīng)常需要求解常微分方程組: 本課程主要研究問題一階常微分方程(1)的數(shù)值解法,我們首先介紹初值問題(1)的解存在的條件-(3)(1),(2)式稱為初值問題,(定理 只要 f (

2、x, y) 連續(xù)，且關(guān)于 y 滿足 Lipschitz 條件，即存在與 x, y 無關(guān)的常數(shù) L 使對任意定義在 a, b 上的 y1(x) 和 y2(x) 都成立，則初值問題（1）存在唯一解。（通常采用等距節(jié)點）對于問題(1)要求它的數(shù)值解定理 只要 f (x, y) 連續(xù)，且關(guān)于 y 滿足 常微分方程數(shù)值解公式的推導(dǎo) 求初值問題數(shù)值解的方法是步進法，即從已知的初值y0出發(fā)，通過一定的計算求y1 ，然后由y1或y0和y1求出y2 ，依次計算到y(tǒng)n ，即在計算出yk后計算yk+1 ，這時有單步法：計算yk+1時，只利用yk多步法：計算yk+1時，用到y(tǒng)k, yk-1, yk-2,常微分方程數(shù)值

3、解公式的主要推導(dǎo)方法泰勒展開利用差商利用數(shù)值積分法常微分方程數(shù)值解公式的推導(dǎo) 求初值問題數(shù)值解的方法是1、泰勒展開的求解思路：將 按泰勒級數(shù)展開用 的近似值 代入上式右端，記所得結(jié)果為，則得到數(shù)值解序列的計算公式：1、泰勒展開的求解思路：將 按泰勒2、化導(dǎo)數(shù)為差商的求解方法思路：若在點 處的導(dǎo)數(shù)用差商來近似代替，如向前差商則微分方程初值問題化為將近似號改為等號，精確解 改為近似解 ，得2、化導(dǎo)數(shù)為差商的求解方法思路：若在點 處的導(dǎo)數(shù)用差商來3、數(shù)值積分的求解思路：如果將微分方程 化為， 然后在各小區(qū)間 上對其兩邊進行積，即如用矩形數(shù)值積分公式可得：3、數(shù)值積分的求解思路：如果將微分方程 以上三

4、種方法推導(dǎo)出同一個數(shù)值求解公式:這個數(shù)值公式稱為歐拉(Euler)公式。以上三種方法推導(dǎo)出同一個數(shù)值求解公式:這個數(shù)值公式稱為歐拉(6.1 歐拉方法一、 歐拉格式：x0 x1向前差商近似導(dǎo)數(shù)記為歐拉公式幾何意義 用一條通過初始點的折線近似表示解曲線 ,亦稱為歐拉折線法 ，或稱為矩形法。一般形式1、顯式歐拉公式6.1 歐拉方法一、 歐拉格式：x0 x1向前差商近似導(dǎo)數(shù)第六章常微分方程初值問題的數(shù)值解法-課件在假設(shè) yk = y(xk)，即第 k 步計算是精確的前提下，考慮的截斷誤差 Rk = y(xk+1) yk+1 稱為局部截斷誤差 。定義 若某算法的局部截斷誤差為O(hp+1)，則稱該算法有

5、p 階精度。定義 歐拉法的局部截斷誤差：歐拉法具有 1 階精度。局部截斷誤差和階數(shù)在假設(shè) yk = y(xk)，即第 k 步計算是精確的2、隱式歐拉格式向后差商近似導(dǎo)數(shù)x0 x1)(,()(1101xyxfhyxy+)1,.,0(),(111-=+=+nkyxfhyykiik由于未知數(shù) yk+1 同時出現(xiàn)在等式的兩邊，不能直接得到，故稱為隱式 歐拉公式，而前者稱為顯式 歐拉公式。一般先用顯式計算一個初值，再迭代求解。 隱式歐拉法的局部截斷誤差：即隱式歐拉公式具有 1 階精度。2、隱式歐拉格式向后差商近似導(dǎo)數(shù)x0 x1)(,()(110二、兩步歐拉格式（中點公式）中心差商近似導(dǎo)數(shù)x0 x2x1假

6、設(shè) ,則可以導(dǎo)出即兩步歐拉格式具有 2 階精度。該方法需要2個初值 y0和 y1來啟動遞推過程，這樣的算法稱為雙步法。二、兩步歐拉格式（中點公式）中心差商近似導(dǎo)數(shù)x0 x2x1假設(shè)三、 梯形公式 顯、隱式兩種算法的平均注：有局部截斷誤差 ， 即梯形公式具有2 階精度，比歐拉方法有了進步。但注意到該公式是隱式公式，計算時不得不用到迭代法，不易求解。 對歐拉法進行改進，用梯形公式計算右側(cè)積分，即計算公式三、 梯形公式 顯、隱式兩種算法的平均注：有局部截斷誤差 梯形格式算法計算步驟：先用(1)式計算出 處 。再用(2)式反復(fù)進行迭代，得到計算公式-(1)-(2)類似地得到 用 控制迭代次數(shù)， 為允許

7、誤差。把滿足誤差要求的 作為 的近似值。 梯形格式算法計算步驟：先用(1)式計算出 處 方 法顯式歐拉法 隱式歐拉法 梯形公式 中點公式 簡單 精度低 穩(wěn)定性最好 精度低, 計算量大 精度提高 計算量大 精度提高, 顯式 多一個初值, 可能影響精度 不同方法比較方 法顯式歐拉法 隱式歐拉法 梯形公式 中點公式四、改進歐拉法（預(yù)報-校正法）Step 1: 先用顯式歐拉公式作預(yù)報，算出),(1iiiiyxfhyy+=+Step 2: 再將 代入隱式梯形公式的右邊作校正，得到1+iy),(),(2111+=iiiiiiyxfyxfhyy它可表示為嵌套形式表示為平均化形式此法稱為預(yù)報-校正法，是顯式算

8、法。注：可以證明該算法具有 2 階精度，同時可以看到它是個單步遞推格式（只迭代一次） ，比隱式梯形公式的迭代求解過程簡單。它的穩(wěn)定性高于顯式歐拉法。腳標(biāo)用 i四、改進歐拉法（預(yù)報-校正法）Step 1: 先用顯式歐拉公舉例：舉例：第六章常微分方程初值問題的數(shù)值解法-課件xEuler法y改進的Euler法y精確解01.0000001.0000001.0000000.11.0000001.0959091.0954450.21.1918181.1840971.1832160.31.2774381.2662011.2649110.41.3582131.3433601.3416410.51.435133

9、1.4164021.4142140.61.5089661.4859561.4832400.71.5803381.5525141.5491930.81.6497831.6164751.6124520.91.7177791.6781661.6733201.01.7847701.7378671.732051xEuler法y改進的Euler法y精確解01.0000006.2 龍格 - 庫塔法一、 泰勒級數(shù)法 龍格庫塔(Runge-Kut ta)法(簡稱為R-K方法)是一類高精度的一步法，這類方法與泰勒級數(shù)法有著密 切的關(guān)系。 設(shè)有初值問題由 泰勒展開式 從理論上講，只要解y(x)有任意階導(dǎo)數(shù)，泰勒展開

10、方法就可以構(gòu)造任意階求yk+1公式。但由于計算這些導(dǎo)數(shù)是非常復(fù)雜的，所以這種方法實際上不能用來解初值問題。 6.2 龍格 - 庫塔法一、 泰勒級數(shù)法 設(shè)有初值問題二、龍格庫塔法的基本思路等價于：（積分中值定理） R-K方法基本思想：用 在幾個不同點的加權(quán)平均值（線性組合）來代替準(zhǔn)確的 的值，構(gòu)造近似公式。再把近似公式與解的泰勒展開 式進行比較，使前面的若干項相同，從而使近似公式達到一定的階數(shù)。 這樣龍格庫塔法保留了泰勒級數(shù)展開法的高階局部截斷誤差，又避免了高階導(dǎo)數(shù)的計算。 我們先分析歐拉法 與預(yù)估校正法。設(shè)有初值問題二、龍格庫塔法的基本思路等價于：（積分中值定理第六章常微分方程初值問題的數(shù)值解

11、法-課件推廣這種單步法稱為Runge-Kutta方法,簡記為R-K公式.Ki為某些點上的斜率，或f(x,y)在某些點上的值。推廣這種單步法稱為Runge-Kutta方法,簡記為R-K公三、 二階龍格 - 庫塔法目標(biāo)：建立高精度的單步遞推格式。單步遞推法的基本思想是從 ( xi , yi ) 點出發(fā)，以某一斜率沿直線達到 ( xi+1 , yi+1 ) 點。歐拉法及其各種變形所能達到的最高精度為2階。 考察改進的歐拉法，可以將其改寫為：斜率一定取K1 K2 的平均值嗎？ 步長一定是一個h 嗎？ 腳標(biāo)用 i三、 二階龍格 - 庫塔法目標(biāo)：建立高精度的單步遞推格式。首先希望能確定系數(shù) 1、2、p，使

12、得到的算法格式有2階精度，即在 的前提假設(shè)下，使得 Step 1: 將 K2 在 ( xi , yi ) 點作 Taylor 展開將改進歐拉法推廣為：),(),(12122111phKyphxfKyxfKKKhyyiiiiii+=+=+llStep 2: 將 K2 代入yi+1表達式，得到首先希望能確定系數(shù) 1、2、p，使得到的算法格式有2階精Step 3: 將 yi+1 與 y( xi+1 ) 在 xi 點的泰勒展開作比較要求 ，則必須有：這里有 個未知數(shù)， 個方程。32存在無窮多個解。所有滿足上式的格式統(tǒng)稱為2階龍格 - 庫塔格式。注意到， 就是改進的歐拉法。 Q: 為獲得更高的精度，應(yīng)該

13、如何進一步推廣？Step 3: 將 yi+1 與 y( xi+1 ) 在 x其中i ( i = 1, , m )，i ( i = 2, , m ) 和 ij ( i = 2, , m; j = 1, , i1 ) 均為待定系數(shù)，確定這些系數(shù)的步驟與前面相似。其解不唯一。).,(.),(),(),(.1122112321313312122122111-+=+=+=+=mm mmmmimiiiiiimmiihKhKhKyhxfKhKhKyhxfKhKyhxfKyxfKKKKhyybbbabbaballl最常用為四階經(jīng)典龍格-庫塔法其中i ( i = 1, , m )，i ( i 四階經(jīng)典龍格-庫塔

14、法公式四、 四階龍格 - 庫塔法用四個f函數(shù)值的線性組合得到四階龍格 - 庫塔法。經(jīng)典龍格-庫塔法公式具有四階精度，因此可取大步長。四階經(jīng)典龍格-庫塔法公式四、 四階龍格 - 庫塔法用四個f注： 龍格-庫塔法的主要運算在于計算 Ki 的值，即計算 f 的值。Butcher 于1965年給出了計算量與可達到的最高精度階數(shù)的關(guān)系：753可達到的最高精度642每步須算Ki 的個數(shù)高于四階時每步計算量增加較多，但精度提高不快，因此使用的比較少。由于龍格-庫塔法的導(dǎo)出基于泰勒展開，故精度主要受解函數(shù)的光滑性影響。對于光滑性不太好的解，最好采用低階算法而將步長h 取小。注：753可達到的最高精度642每步

15、須算Ki 的個數(shù)高于四Lab 15. Runge-Kutta (Order Four) Use Runge-Kutta method of order four to approximate the solution of the initial-value problem, , and (1)You are supposed to write a functionvoid RK4 ( double (*f)( ), double a, double b, double y0, int n, FILE *outfile)to approximate the solution of Proble

16、m (1) with y= f, x in a, b, and the initial value of y being y0. Output the approximating values of y on the n+1 equally spaced grid points from a to b to outfile.InputThere is no input file. Instead, you must hand in your function in a *.h file. The rule of naming the *.h file is the same as that o

17、f naming the *.c or *.cpp files.Output ( represents a space) For each test case, you are supposed to print n+1 lines, and each line is in the following format: fprintf( outfile, “%8.4f%12.8en”, x, y );Lab 15. Runge-Kutta (Order FoSample Judge Program#include #include #include 98115001_15.hdouble f0(

18、double x, double y) return (yx*x+1.0); void main( ) FILE *outfile = fopen(out.txt, w);int n;double a, b, y0;a = 0.0; b = 2.0; y0 = 0.50; n = 10;RK4(f0, a, b, y0, n, outfile);fprintf(outfile, n);fclose(outfile);Sample Output ( represents a space) 0.00005.00000000e0010.20008.29293333e0010.40001.214076

19、21e+0000.60001.64892202e+0000.80002.12720268e+0001.00002.64082269e+0001.20003.17989417e+0001.40003.73234007e+0001.60004.28340950e+0001.80004.81508569e+0002.00005.30536300e+000Sample Judge ProgramSample Out3 收斂性與穩(wěn)定性 /* Convergency and Stability */ 收斂性 /* Convergency */定義 若某算法對于任意固定的 x = xi = x0 + i h

20、，當(dāng) h0 ( 同時 i ) 時有 yi y( xi )，則稱該算法是收斂的。 例：就初值問題 考察歐拉顯式格式的收斂性。解：該問題的精確解為 歐拉公式為對任意固定的 x = xi = i h ，有3 收斂性與穩(wěn)定性 /* Convergency and 穩(wěn)定性 /* Stability */例：考察初值問題 在區(qū)間0, 0.5上的解。分別用歐拉顯、隱式格式和改進的歐拉格式計算數(shù)值解。0.00.10.20.30.40.5精確解改進歐拉法 歐拉隱式歐拉顯式 節(jié)點 xi 1.00002.0000 4.00008.0000 1.6000101 3.2000101 1.00002.5000101 6.25001021.56251023.90631039.76561041.00002.50006.25001.56261013.90631019.76561011.00004.97871022.47881031.234110