精選第8章 常微分方程初值問題的數(shù)值解法

1、1 第8章 常微分方程初值問 題的數(shù)值解法 ?基礎(chǔ)知識(shí) ?歐拉方法 ?龍格庫塔方法 2 8.1 基礎(chǔ)知識(shí) 一、問題的提出 很多實(shí)際問題都需要求解常微分方程。例如單擺問題。 常微分方程分為線性常微分方程和非線性常微分方程，又可以分為一階 常微分方程和高階常微分方程。通過變量的替換，可以把高階常微分方 程轉(zhuǎn)化為一階常微分方程再求解。對于一階常微分方程組，可以寫成向 量形式的單個(gè)方程，求解方法與一階常微分方程相似。因此本章只討論 一階常微分方程的初值問題： ? ? ? ? ? ? ? 00 )( , ),( yxy bxayxf dx dy 目前在常微分方程理論中，只能求出某些特殊類型常微分方程的解

2、析解， 對大部分常微分方程，用解析方法求出常微分方程的精確解非常困難， 甚至不存在解的解析表達(dá)式。為滿足工程實(shí)踐的需要，常常用數(shù)值解法 求常微分方程的近似解。 在本章中假設(shè)討論的一階常微分方程的初值問題的解y(x)存在、唯一且存在、唯一且 足夠光滑，方程本身是穩(wěn)定的，即精確解y(x)連續(xù)且依賴于初始值及右 端函數(shù)。 3 8.1 基礎(chǔ)知識(shí) 二、數(shù)值解法 一階常微分方程初值問題的數(shù)值解法的主要思想，是對區(qū)間a,b上的節(jié)點(diǎn) a=x 0 x1xnxn+1b 建立y(xn)的近似值yn的某一遞推格式，利用初值y0和已計(jì)算出的和已計(jì)算出的y1,y2, yk-1遞推出yk，并且用這一方法反復(fù)遞推，依次得到y(tǒng)

3、k+1,yk+2,yn。這一 求解方法稱為步進(jìn)式求解，相鄰2個(gè)節(jié)點(diǎn)的距離稱為步長，記為hi=xi+1xi。 為便于計(jì)算，常取成等距節(jié)點(diǎn)，稱為定步長，這時(shí)把步長記為 h。 一階常微分方程初值問題的數(shù)值解法有多種分類方法。 多步法不能自行啟動(dòng)，必須先用單步法計(jì)算出y1,y2,yk-1，才能啟動(dòng)一 個(gè)r步的多步法。 一種分類方法為： 單步法：每一輪遞推只用到前面一輪的遞推結(jié)果，遞推格式為： yk=yk-1hT(x k-1,yk-1) 多步法：每一輪遞推要用到前面多輪遞推的結(jié)果，遞推格式為： yk=yk-1hT(x k-r,yk-r,xk-r+1 ,yk-r+1, xk-1,yk-1)，其中r1。 4

4、 8.1 基礎(chǔ)知識(shí) 二、數(shù)值解法（續(xù)） 另一種分類方法為： 顯式方法：遞推公式的右端都是已知量，可以直接計(jì)算出遞推的結(jié) 果，遞推格式為：yk=yk-1hT(x k-r,yk-r,xk-r+1 ,yk-r+1, xk-1,yk-1) 隱式方法：遞推公式左端的未知量也出現(xiàn)在公式的右端，遞推格式 為： yk=yk-1hT(x k-r,yk-r,xk-r+1 ,yk-r+1, xk,yk) 隱式方法的遞推公式其實(shí)是一個(gè)方程。解方程的運(yùn)算量可能較大，為避免 解方程，常采用預(yù)測校正系統(tǒng)。 預(yù)測校正系統(tǒng)：每一輪遞推包括預(yù)測和校正這2個(gè)步驟。先用顯式方 法計(jì)算出yk，作為迭代的初值，這一過程稱為預(yù)測；再把隱式

5、方法的遞推 公式作為迭代公式，把預(yù)測值yk代入迭代公式右端進(jìn)行迭代，這一過程稱 為校正。在校正時(shí)往往迭代1次或幾次，校正值的精度就會(huì)有大幅提高。 一階常微分方程初值問題的數(shù)值解法一般是對連續(xù)的初值問題進(jìn)行離散化 處理，把微分方程轉(zhuǎn)化為代數(shù)方程來求解。常用的離散化方法有： 基 于數(shù)值微分的離散化方法， 基于數(shù)值積分的離散化方法， 基于泰勒 展開的離散化方法。 5 8.2 歐拉方法 一、顯式歐拉法 在一階常微分方程初值問題的數(shù)值解法中，顯式歐拉（Euler ）法是最簡 單的一種。顯式歐拉法有明顯的幾何含義，缺點(diǎn)是精度不高。對于一階 常微分方程的初值問題： ? ? ? ? ? ? ? 00 )( ,

6、 ),( yxy bxayxf dx dy 顯式歐拉法的遞推公式為：yk=yk-1hf(xk-1,yk-1)，k=1,2,3,。 顯式歐拉法每一輪遞推只用到前面一輪遞推的結(jié)果，因此它是單步法。 由基于數(shù)值微分的離散化方法、基于數(shù)值積分的離散化方法、基于泰 勒展開的離散化方法都可以推導(dǎo)出顯式歐拉法的遞推公式。 用顯式歐拉法求解一階常微分方程初值問題的過程，就是以已知的 (x0,y0) 作為起點(diǎn)，代入顯式歐拉法的遞推公式的右端，計(jì)算出y(x)在x1處的近似 值y1；再以(x1,y1)作為起點(diǎn)，用顯式歐拉法的遞推公式計(jì)算出y(x)在x2處 的近似值y2；。 6 8.2 歐拉方法 一、顯式歐拉法（續(xù)）

7、 顯式歐拉法的遞推公式是斜率為f(xk-1,yk-1)，經(jīng)過點(diǎn)(xk-1,yk-1)的直線方程。 上述顯式歐拉法遞推過程的幾何含義，就是用曲線y(x)在點(diǎn)(x0,y0)處的切 線段代替y(x)在區(qū)間x 0,x1內(nèi)的曲線段；再把曲線段的終點(diǎn)(x1,y(x1)近似 為切線段的終點(diǎn)(x1,y1)，把曲線y(x)在點(diǎn)(x1,y1)處切線的斜率處切線的斜率f(x1,y(x1)近 似為f(x1,y1) ，做曲線y(x)在點(diǎn)(x1,y1)處的切線，并在區(qū)間x 1,x2內(nèi)用近似 的切線段代替y(x)的曲線段；。歐拉法用一系列折線近似地代替曲線 y(x)，因此它又稱為歐拉折線法。 定義：假設(shè)某一步遞推的起點(diǎn)是精

8、確的，即yi-1=y(xi-1)，那么這一步遞推 的截?cái)嗾`差R i=y(xi)yi稱為局部截?cái)嗾`差。 定義：若某算法的局部截?cái)嗾`差為O(h p+1)，則稱此算法有p階精度。 顯式歐拉法具有一階精度。 顯式歐拉法在步長過大時(shí)誤差較大；在步長較小時(shí)需要多步遞推，可能 出現(xiàn)誤差積累的現(xiàn)象。由于顯式歐拉法的精度不高，因此在實(shí)際應(yīng)用中 用的較少。 7 8.2 歐拉方法 顯式歐拉法的算法 輸入求解區(qū)間的邊界a,b，步長h，起點(diǎn)x0=a處的縱坐標(biāo)y0。 遞推的次數(shù)n=(b-a)/h。 for(i=0;in;i+) /循環(huán)1次完成1輪遞推 xi+1=xi+h; yi+1=yi+h*f(xi,yi); for(

9、i=0;i=n;i+) /輸出計(jì)算結(jié)果 輸出xi,yi。 8 8.2 歐拉方法 顯式歐拉法對應(yīng)的程序 #include #define MAXSIZE 50 double f(double x,double y); void main(void) double a,b,h,xMAXSIZE,yMAXSIZE; long i,n; printf( 請輸入求解區(qū)間a,b：); scanf(%lf,%lf, printf( 請輸入步長h：); scanf(%lf, n=(long)(b-a)/h); x0=a; printf( 請輸入起點(diǎn)x0=%lf處的縱坐標(biāo)y0：,x0); scanf(%lf,

10、for(i=0;in;i+) xi+1=xi+h; yi+1=yi+h*f(xi,yi); printf( 計(jì)算結(jié)果為：計(jì)算結(jié)果為：); for(i=0;i=n;i+) printf( x%ld=%lf,y%ld=%lf,i,xi,i,yi); double f(double x,double y) return(); /* 計(jì)算并返回函數(shù)值f(x,y)*/ 9 8.2 歐拉方法 二、歐拉方法的變形 方法一：隱式歐拉法（后退的歐拉法） 遞推公式為： yk=yk-1hf(x k,yk)，k=1,2,3,。 用隱式歐拉法求解一階常微分方程初值問題的過程，就是以已知的 (x0,y0)作為起點(diǎn)，計(jì)算出

11、y(x)在x1處的近似值y1；再以(x1,y1)作為起點(diǎn)， 計(jì)算出y(x)在x2處的近似值y2；，象這樣反復(fù)地遞推。 顯然，隱式歐拉法是隱式方法，遞推公式的右端有未知量yk。隱式歐拉 法需要求解1個(gè)方程。為避免解方程，常用顯式歐拉法的計(jì)算結(jié)果作為迭 代的初值yk(0)，把隱式歐拉法的遞推公式作為迭代公式反復(fù)迭代，得到 迭代序列yk(0),yk(1),yk(2),。如果步長足夠小，那么迭代序列收斂于yk。 隱式歐拉法具有一階精度。 隱式歐拉法精度不高，計(jì)算復(fù)雜，用的比較少。 隱式歐拉法由點(diǎn)(x0,y0)遞推到點(diǎn)(x1,y1)的幾何含義，是把曲線y(x)在點(diǎn) (x1,y1)處的切線平行移動(dòng)，移動(dòng)到

12、經(jīng)過點(diǎn)(x0,y0) ，在區(qū)間x 0,x1內(nèi)用此 直線段代替y(x)的曲線段。 10 8.2 歐拉方法 二、歐拉方法的變形 方法二：梯形公式法 遞推公式為： yk=yk-1(h/2)*(f(xk-1,yk-1)+f(xk,yk)，k=1,2,3,。 用梯形公式法求解一階常微分方程初值問題的過程，就是以已知的 (x0,y0)作為起點(diǎn)，計(jì)算出y(x)在x1處的近似值y1；再以(x1,y1)作為起點(diǎn)， 計(jì)算出y(x)在x2處的近似值y2；，象這樣反復(fù)地遞推。 梯形公式法具有2階精度。 顯然，梯形公式法是隱式方法，需要求解方程。為避免解方程，常用 顯式歐拉法的計(jì)算結(jié)果作為迭代的初值yk(0)，把梯形公

13、式法的遞推公 式作為迭代公式反復(fù)迭代，得到迭代序列yk(0),yk(1),yk(2),。如果步 長h足夠小，那么迭代序列收斂于yk。 梯形公式法由點(diǎn)(x0,y0)遞推到點(diǎn)(x1,y1)的幾何含義，是經(jīng)過點(diǎn)(x0,y0)做一 直線段，在區(qū)間x 0,x1內(nèi)用此直線段代替y(x)的曲線段。此直線段的斜 率等于曲線y(x)在點(diǎn)(x0,y0)處的切線斜率和曲線y(x)在點(diǎn)(x1,y1)處的切線 斜率的平均值。類似地，由點(diǎn)(xk-1,yk-1)遞推出點(diǎn)(xk,yk)，k=1,2,3,。 11 8.2 歐拉方法 二、歐拉方法的變形 方法三：中點(diǎn)歐拉方法（兩步歐拉方法 ） 遞推公式為： yk+1=yk-12h

14、f(xk,yk)，k=1,2,3,。 中點(diǎn)歐拉方法是雙步法，需要2個(gè)初值y0和y1才能啟動(dòng)遞推過程。一般 先用單步法由點(diǎn)(x0,y0)計(jì)算出(x1,y1)，再用中點(diǎn)歐拉方法反復(fù)地遞推。 中點(diǎn)歐拉方法由點(diǎn)(x0,y0)和(x1,y1)遞推出點(diǎn)(x2,y2)的幾何含義，是經(jīng)過 點(diǎn)(x0,y0)做一直線段，在區(qū)間x 0,x2內(nèi)用此直線段代替y(x)的曲線段。 此直線段的斜率等于曲線y(x)在點(diǎn)(x1,y1)處的切線斜率。類似地，由 點(diǎn)(xk-1,yk-1)和點(diǎn)(xk,yk)遞推出點(diǎn)(xk+1,yk+1) ，其中k=1,2,3,。 中點(diǎn)歐拉方法具有2階精度。 12 8.2 歐拉方法 二、歐拉方法的變形

15、 上述各種方法的比較： 由拉格郎日中值定理，在區(qū)間(xk-1,xk)內(nèi)必定存在k，使得： y(xk) y(xk-1)=y( k)(xkxk-1) y(xk)=y(xk-1)+hy( k)=y(xk-1)+hf(k,y(k)， 如果知道k，代入這個(gè)遞推公式，那么遞推過程得到的序列y0,y1, y2,沒有誤差。求k往往很困難，因此常用一個(gè)易求的值近似地代替 k 。顯式歐拉法、隱式歐拉法、梯形公式法、中點(diǎn)歐拉方法的區(qū)別是 對k 、y(k)的近似方法不同。 ?顯式歐拉法把k近似為區(qū)間x k-1,xk的起點(diǎn)xk-1； ?隱式歐拉法把k近似為區(qū)間x k-1,xk的終點(diǎn)xk； ?梯形公式法把處k的導(dǎo)數(shù)y(k

16、)近似為區(qū)間x k-1,xk的起點(diǎn)和終點(diǎn)導(dǎo)數(shù) 的平均值； ?中點(diǎn)歐拉方法考察的區(qū)間為x k-1,xk+1，k(xk-1,xk+1)，k被近似為區(qū) 間x k-1,xk+1的中點(diǎn)xk。 13 8.2 歐拉方法 三、改進(jìn)的歐拉法 改進(jìn)的歐拉法是一種預(yù)測校正方法，它的每一輪遞推包括預(yù)測和校正 這2個(gè)步驟：先用顯式歐拉公式計(jì)算出y*k， y*k=yk-1hf(xk-1,yk-1) 這一步稱為“預(yù)測”；再用梯形公式迭代一次， yk=yk-1(h/2)*(f(x k-1,yk-1)+f(xk,y*k) 這一步稱為“校正”。改進(jìn)的歐拉法精度比顯式歐拉法高，不需要解方 程，是一種更實(shí)用的方法。 14 8.2 歐

17、拉方法 改進(jìn)的歐拉法的算法 輸入求解區(qū)間的邊界a,b，步長h，起點(diǎn)x0=a處的縱坐標(biāo)y0。 遞推的次數(shù)n=(b-a)/h。 for(i=0;in;i+) /循環(huán)1次完成1輪遞推 xi+1=xi+h; yi+1=yi+h*f(xi,yi); / 預(yù)測 yi+1=yi+h*(f(xi,yi)+f(xi+1,yi+1)/2; / 校正 for(i=0;i=n;i+) /輸出計(jì)算結(jié)果 輸出xi,yi。 15 8.2 歐拉方法 改進(jìn)的歐拉法對應(yīng)的程序 #include #define MAXSIZE 50 double f(double x,double y); void main(void) doub

18、le a,b,h,xMAXSIZE,yMAXSIZE; long i,n; printf( 請輸入求解區(qū)間a,b：); scanf(%lf,%lf, printf( 請輸入步長h：); scanf(%lf, n=(long)(b-a)/h); x0=a; printf( 請輸入起點(diǎn)x0=%lf處的縱坐標(biāo)y0：,x0); scanf(%lf, for(i=0;in;i+) xi+1=xi+h; yi+1=yi+h*f(xi,yi); yi+1=yi+h*(f(xi,yi)+f(xi+1,yi+1)/2; printf( 計(jì)算結(jié)果為：); for(i=0;i=n;i+)printf( x%ld=%

19、lf,y%ld=%lf,i,xi,i,yi); double f(double x,double y) return(); /* 計(jì)算并返回函數(shù)值f(x,y)*/ 16 8.3 龍格庫塔方法 一、泰勒展開方法 顯式歐拉法、隱式歐拉法、梯形公式法、中點(diǎn)歐拉方法的缺陷是精度較 低。本節(jié)介紹高精度的一步法。 直接用泰勒公式展開是一種高階顯式一步法。在由點(diǎn)(xk-1,yk-1)遞推出點(diǎn) (xk,yk)時(shí)，把y(xk)在x=xk-1處做泰勒展開： y(xk)=y(xk-1)+hy(x k-1)+(h 2/2!)y(x k-1)+(h n/n!)y(n)(xk-1)+R n(x) 其中余項(xiàng)Rn(x)=(h

20、 n+1 /(n+1)!)y(n+1)()，(xk-1,xk) 泰勒公式中的各階導(dǎo)數(shù)依次為：y(x)=f，y(x)=f x+fy f， 以已知的(x0,y0)作為起點(diǎn)，依次地用泰勒公式由點(diǎn)(xk-1,yk-1)遞推出點(diǎn) (xk,yk)，k=1,2,3,，完成求解的過程。 如果y(x)在考察的區(qū)間內(nèi)具有直到(n+1)階導(dǎo)數(shù)，那么泰勒公式可以 計(jì)算到前(n+1)項(xiàng)，此時(shí)泰勒展開方法具有階精度。 從理論上講，如果y(x)在考察的區(qū)間內(nèi)足夠光滑，那么泰勒展開方 法可以具有任意階精度。但是高階泰勒展開方法計(jì)算量大，求 f(x,y(x)的高階導(dǎo)數(shù)困難，因此泰勒展開方法并不實(shí)用。 17 8.3 龍格庫塔方法

21、 二、龍格庫塔法基本思想 龍格庫塔（RungeKutta ）法，簡稱RK方法，是一種高階顯式一步 法，而且不需要計(jì)算導(dǎo)數(shù)。Runge首先提出了間接使用泰勒展開式的方法： 用f(x,y)在一些點(diǎn)上函數(shù)值的線性組合代替y(x)的各階導(dǎo)數(shù)，構(gòu)造yn+1的表 達(dá)式，比較這個(gè)表達(dá)式與y(xn+1)在x=xn處泰勒展開式，確定yn+1的表達(dá)式中 的系數(shù)，使yn+1的表達(dá)式與y(xn+1)泰勒展開式前面的若干項(xiàng)相等，從而具 有對應(yīng)泰勒展開式的精度階次，這就是龍格庫塔法的主要思想。 p級(jí)龍格庫塔公式的一般形式為： yn+1=yn+h ，n=1,2,3,。 ? ? ? p i iik c 1 其中k1=f(xn

22、,yn)，ki= ，i=2,3,p。 ),( 1 1 ? ? ? ? i j jijnin kbhyhaxf 上式中ci、ai、b ij是與具體常微分方程初值問題以及n、h無關(guān)的常數(shù)。 無關(guān)的常數(shù)。k i 是f(x,y)在某些點(diǎn)處函數(shù)值，ci是在線性組合求yn+1時(shí)k i的“權(quán)”，ai和bij用 來確定k i在f(x,y)上的位置。 以p級(jí)龍格庫塔公式作為遞推公式，以已知的(x0,y0)作為起點(diǎn)，依次地 由點(diǎn)(xn,yn)遞推出點(diǎn)(xn+1,yn+1)，n=0,1,2,，完成求解的過程。 18 8.3 龍格庫塔方法 二、龍格庫塔法基本思想（續(xù)） 確定常數(shù)ci、ai、bij的原則是使龍格庫塔公式

23、與泰勒公式前面的盡可能 多的項(xiàng)相等。如果p級(jí)龍格庫塔公式等于泰勒公式的前(q+1)項(xiàng)，那么這 個(gè)龍格庫塔公式具有q階精度，稱此公式為p級(jí)q階龍格庫塔公式。 一級(jí)一階龍格庫塔公式就是顯式歐拉公式。 二級(jí)龍格庫塔公式與泰勒公式前面3項(xiàng)相等，具有二階精度。這個(gè)方 程組有3個(gè)方程，4個(gè)變元，因此它有1個(gè)自由參數(shù)，二級(jí)二階龍格庫 塔公式有無窮多個(gè)。改進(jìn)的歐拉法（預(yù)測校正法）和中點(diǎn)歐拉方法 （兩步歐拉方法）是二級(jí)龍格庫塔公式的特例。 標(biāo)準(zhǔn)（經(jīng)典）龍格庫塔公式為： yn+1=yn+(h/6)(k 1+2k2+2k3+k4)， 其中k1=f(xn,yn)， k2=f(xn+h/2,yn+k1h/2)， k3=f(xn+h/2,yn+k2h/2)， k4=f(xn