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

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

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

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

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

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

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

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

8、法。注：可以證明該算法具有 2 階精度，同時(shí)可以看到它是個(gè)單步遞推格式（只迭代一次） ，比隱式梯形公式的迭代求解過(guò)程簡(jiǎn)單。它的穩(wěn)定性高于顯式歐拉法。腳標(biāo)用 i四、改進(jìn)歐拉法（預(yù)報(bào)-校正法）Step 1: 先用顯式歐拉公舉例：舉例：第六章常微分方程初值問(wèn)題的數(shù)值解法-課件xEuler法y改進(jìn)的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改進(jìn)的Euler法y精確解01.0000006.2 龍格 - 庫(kù)塔法一、 泰勒級(jí)數(shù)法 龍格庫(kù)塔(Runge-Kut ta)法(簡(jiǎn)稱(chēng)為R-K方法)是一類(lèi)高精度的一步法，這類(lèi)方法與泰勒級(jí)數(shù)法有著密 切的關(guān)系。 設(shè)有初值問(wèn)題由 泰勒展開(kāi)式 從理論上講，只要解y(x)有任意階導(dǎo)數(shù)，泰勒展開(kāi)

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

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

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

13、如何進(jìn)一步推廣？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)典龍格-庫(kù)塔法其中i ( i = 1, , m )，i ( i 四階經(jīng)典龍格-庫(kù)塔

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

15、須算Ki 的個(gè)數(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 */定義 若某算法對(duì)于任意固定的 x = xi = x0 + i h

20、，當(dāng) h0 ( 同時(shí) i ) 時(shí)有 yi y( xi )，則稱(chēng)該算法是收斂的。 例：就初值問(wèn)題 考察歐拉顯式格式的收斂性。解：該問(wèn)題的精確解為 歐拉公式為對(duì)任意固定的 x = xi = i h ，有3 收斂性與穩(wěn)定性 /* Convergency and 穩(wěn)定性 /* Stability */例：考察初值問(wèn)題 在區(qū)間0, 0.5上的解。分別用歐拉顯、隱式格式和改進(jìn)的歐拉格式計(jì)算數(shù)值解。0.00.10.20.30.40.5精確解改進(jìn)歐拉法 歐拉隱式歐拉顯式 節(jié)點(diǎn) 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