剛性微分方程組隱式龍格庫塔方法

1、 畢業(yè)設(shè)計(jì)（論文）畢 業(yè) 設(shè) 計(jì)題 目：剛性系統(tǒng)的隱式RK方法學(xué) 院： 數(shù)理學(xué)院 專業(yè)名稱： 信息與計(jì)算科學(xué) 學(xué) 號(hào)： 201241210127 學(xué)生姓名： 丁楠 指導(dǎo)教師： 汪玉霞 2016年05月15日摘要本文主要介紹單步隱式Runge Kutta方法,簡(jiǎn)要的介紹了Gauss型隱式Runge Kutta方法、Radau型隱式Runge Kutta方法和Lobatto型隱式Runge Kutta方法。并利用這些基本的隱式Runge Kutta方法來對(duì)剛性方程組進(jìn)行數(shù)值求解，并將隱式Runge Kutta方法與顯式經(jīng)典Runge Kutta方法求解的結(jié)果進(jìn)行對(duì)比，說明兩種數(shù)值解法的優(yōu)缺點(diǎn)。關(guān)鍵

2、詞：剛性系統(tǒng) 隱式Runge Kutta方法 單步方法 Newton迭代法AbstractThis paper mainly introduces the Implicit Runge-Kutta Methods and a simple description of Gauss implicit Runge-Kutta method ,Radau implicit Runge-Kutta method and Lobatto implicit Runge-Kutta method . These basic Implicit Runge-Kutta methods are used to s

3、olve the stiff equations. These implicit Runge-Kutta methods iare compared with the classical explicit Runge-Kutta method. This paper explain the advantages and disadvantages of the two kind of numerical methods.Keywords: Stiff system Implicit Runge-Kutta method One step method Newton iterative meth

4、od目錄1.緒論11.1剛性方程11.2隱式RK方法的研究意義21.3 RK方法的研究現(xiàn)狀32.單步RK方法的收斂性和穩(wěn)定性52.1單步RK方法的收斂性52.2單步RK方法的穩(wěn)定性63.三類隱式RK方法83.1引言83.2 Gauss型隱式RK方法93.3 Radau型隱式RK方法103.4 Lobatto型隱式RK方法114隱式RK方法的實(shí)現(xiàn)134.1非線性系統(tǒng)的改進(jìn)134.2簡(jiǎn)化的Newton迭代法135數(shù)值實(shí)驗(yàn)與結(jié)果分析15參考文獻(xiàn)18附錄211. 緒論1.1剛性方程對(duì)于一般的線性常系數(shù)系統(tǒng)y'=Ay+(t)A為m×m的矩陣，特征值為i(i=1,2,m)。定義123 若

5、一個(gè)系統(tǒng)滿足（1）Rei<0, i=1,2,m（2）maxiReiminiRei=R1其中R為剛性比，則這個(gè)系統(tǒng)稱為剛性系統(tǒng)。定義227 若線性系統(tǒng)y'=Ay x0,T或非線性系統(tǒng)y'=f(x,y) x0,T的矩陣A或Jacobi矩陣fy的特征值i滿足max1imRei1則其是剛性的。定理1（解的存在性與唯一性）（1）對(duì)于所有x,yD，函數(shù)fx,y是連續(xù)的；（2）對(duì)于任何x,y,(x,y*)D，存在常數(shù)L，是函數(shù)滿足fx,y-f(x,y*)Ly-y*則初值問題y=fx,y axbya= 有唯一解。其中y=(y1,y2,ym)T，D=x,y|axb,-<a<b&

6、lt;-<yi<,i=1,2,m 。其中L被稱為Lipschitz常數(shù)定義3 如果一個(gè)常微分系統(tǒng)的Lipschitz常數(shù)L很大(大于20)，則它是剛性的。1.2隱式RK方法的研究意義在常微分方程及常微分方程組的數(shù)值解法中，Runge Kutta方法是目前應(yīng)用最為廣泛的數(shù)值解法之一，同時(shí)又具有誤差小，精度高的特點(diǎn)。盡管顯式Runge Kutta方法能夠非常準(zhǔn)確、快速的給出大部分常微分方程組的數(shù)值解。但是在化學(xué)、自動(dòng)控制電力系統(tǒng)等領(lǐng)域中，會(huì)出現(xiàn)一些病態(tài)的常微分方程組，也就是剛性方程組。剛性方程組對(duì)于數(shù)值解法的穩(wěn)定性要求苛刻，比如方程組y1=-0.01y1-99.99y2, y10=2y

7、2=-100y2, y20=1將其表示為矩陣形式：y1y2=-0.01-99.990-100y1y2令A(yù)=-0.01-99.990-100發(fā)現(xiàn)A特征值為：1=-100，2=-0.01，剛性比s=12=100001。方程組的解為：y1x=e-100x+e-0.01xy2x=e-100x 快瞬態(tài) 慢瞬態(tài)解由快瞬態(tài)和慢瞬態(tài)兩部分構(gòu)成。由于慢瞬態(tài)的部分，y1x衰減變得十分緩慢。當(dāng)自變量變到x=391時(shí)，函數(shù)值還未下降到初值的1%，求解區(qū)間至少取為(0,391)。另一方面，由于快瞬態(tài)的部分，y2x衰減的非?？欤虼瞬介L要取得非常小。從絕對(duì)穩(wěn)定性的方面來看，如果用四階顯式經(jīng)典RK方法求解，絕對(duì)穩(wěn)定區(qū)間要求

8、h(-2.78,0)，則要求h<0.0278。這樣，在(0,391)上就要計(jì)算14 065步，計(jì)算量巨大，因此計(jì)算區(qū)間(0,391)內(nèi)的解時(shí)，舍入誤差積累會(huì)特別嚴(yán)重。例如取求解區(qū)間為0,1，用不同步長h來計(jì)算y11和y21的值。利用四階顯式經(jīng)典RK方法求解如下：hy11y210.042.9802322e+172.9802322e+170.029.9004983e-011.3929556e-240.019.9004983e-012.5300364e-430.0019.9004983e-013.7204130e-440.00019.9004983e-013.7200760e-44真值9.90

9、04983e-013.7200760e-44很顯然，保留八位有效數(shù)字的情況下，要保持良好的精度，步長要取得非常小，這就增加了計(jì)算量。而隨著步長的加大，誤差也會(huì)越來越大。當(dāng)步長加大到絕對(duì)穩(wěn)定與之外時(shí)（即h>0.0278），計(jì)算結(jié)果就完全不可信了。對(duì)于剛性方程組，顯式方法已經(jīng)遠(yuǎn)遠(yuǎn)不能勝任了，一般采用絕對(duì)穩(wěn)定性更好的方法（如隱式Runge Kutta方法）進(jìn)行求解，本文采用單步隱式Runge Kutta方法，而對(duì)于隱式方法中的級(jí)值得求解，本文采用Newton迭代法。1.3 RK方法的研究現(xiàn)狀研究基于標(biāo)準(zhǔn)模型方程的Runge Kutta方法的常見形式為：yi+1=yi+hj=1rjkj k1=f

10、xi,yi kj=f(xi+ih,yi+hl=1j-1jlkl) (j=2,3,r) （顯式）yi+1=yi+hj=1rjkj kj=f(xi+ih,yi+hl=1rjlkl) (j=1,2,3,r) （隱式）因?yàn)镽unge Kutta方法是比較成熟的常微分方程數(shù)值解法。所以如今主要是對(duì)于經(jīng)典的Runge Kutta方法進(jìn)行完善和擴(kuò)充，在一定的條件下，提高級(jí)數(shù)以提高精度?；蛘呤菍unge Kutta方法與某些領(lǐng)域結(jié)合使用。在1994年，費(fèi)景高1給出了一種顯式Runge Kutta并行方法，從而實(shí)現(xiàn)Runge Kutta方法在多處理機(jī)上的應(yīng)用。1997年，Enenkel和Jackson2實(shí)現(xiàn)了

11、Runge Kutta方法的對(duì)角隱式并行改進(jìn)。1999年，廖文遠(yuǎn)和李慶揚(yáng)5給出了一類求解剛性常微分方程的半隱式多步Runge Kutta方法。2000年，張誠堅(jiān)和余洪兵3針對(duì)非線性延遲系統(tǒng)構(gòu)造了一類并行預(yù)校算法，給出其算法的局部誤差估計(jì)，數(shù)值實(shí)驗(yàn)表明該算法是有效的，且具一定的可比性。2003年，李愛雄4通過對(duì)傳統(tǒng)單支方法的計(jì)算格式進(jìn)行改造，得到了解DDEs的兩類單支并行算法單支并行預(yù)校算法和單支并行算法，并對(duì)方法的收斂性和穩(wěn)定性做出了分析。2008年，龐麗君和朱永忠6給出了一類隨機(jī)微分方程Runge Kutta方法的指數(shù)穩(wěn)定性。2.單步RK方法的收斂性和穩(wěn)定性2.1單步RK方法的收斂性對(duì)于常微

12、分初值問題 y=fx,y axbya= (1)的單步顯式公式 yi+1=yi+h(xi,yi,h) i=0,1,n-1y0= (2)局部截?cái)嗾`差可以表示為 Ri+1=yxi+1-yxi+hxi,yi,h i=0,1,n-1 3定理216：設(shè)y(x)為(1)的解，yii=0n為(2)的解。如果：（1）存在常數(shù)c，使得Ri+1chp+1 (i=0,1,2,n-1)（2）存在a>0，使得maxx,yD0ah(x,y,h)yL其中D=x,y|axb,yx-yyx+記c0=cLeLb-a-1，則當(dāng)hmina,pc0 時(shí)，有E(h)c0hp證:由(3)得yxi+1=yxi+hxi,y(xi),h+R

13、i+1 (i=0,1,n-1) (4)將(4)與(2)相減yxi+1-yxi=yxi-yi+hxi,yxi,h-xi,yi,h+ Ri+1 i=0,1,n-1 由yx0-y0=0知道，在i=0時(shí)，yxi-yic0hp成立?，F(xiàn)在假設(shè)0ik-1時(shí)也是成立的。在hpc0時(shí)有：yxi-yic0pc0p= (i=0,1,n-1)進(jìn)而 xi,yxi,h-xi,yi,h =xi,i,hyyxi-yiLyxi-yi 0ik-1其中i是yxi和yi之間的數(shù)。于是定理結(jié)合條件與(4)式，可以得到 yxi+1-yi+1yxi-yi+h xi,yxi,h-xi,yi,h+ Ri+1yxi-yi+Lhyxi-yi+ch

14、p+1=1+Lhyxi-yi+chp+1 0ik-1 從而y(xk-yk)1+Lhkyx0-y0+1+Lhk-11+Lh-1chp+1因?yàn)閥x0-y0=0及1+LhkeLkheL(b-a)，得到y(tǒng)xk-ykcLeLb-a-1hp=c0hp所以當(dāng)i=k時(shí)yxi-yic0hp也是成立的。（證畢）對(duì)于單步顯式格式而言，當(dāng)f(x,y)和f(x,y)y在D內(nèi)連續(xù)時(shí)，那么定理1的條件（2）總是滿足的。所以滿足上述條件的單步顯式公式，只要也滿足相容性條件x,yx,0=f(x,y(x)那么它一定是收斂的。2.2單步RK方法的穩(wěn)定性定義4 對(duì)于初值問題(1),若yii=0n是由(2)得到，zii=0n是下面擾動(dòng)

15、問題的解：zi+1=zi+hxi,zi,h+i+1 (i=0,1,2,n-1)z0=+0 如果存在正常數(shù)C，0及h0，使所有的(0,0，h(0,h0，當(dāng)max0ini時(shí)，有max0inyi-ziC就稱(2)是穩(wěn)定的或者零穩(wěn)定的。定理316 在定理1的條件下，單步顯式公式(2)是穩(wěn)定的。3.三類隱式RK方法3.1引言對(duì)于初值問題(1)，隱式Runge Kutta方法的一般形式為yi+1=yi+hj=1rbjkj (5) kj=fxi+cih,yi+hl=1rajlkl j=1,2,3,r (6)其中，xi=x0+ih，n=0,1,2,，為數(shù)軸上的離散點(diǎn)列；h為步長，yi為解y(xi)的近似值；c

16、1，c2，cr稱為隱式Runge Kutta方法的步長；b1，b2，br為權(quán)系數(shù)。令A(yù)=ajl，j=1,2,3,r。稱A為系數(shù)矩陣，因此上式可以簡(jiǎn)化表示為c1c2a11a1rar1arrb1br使用Butcher陣的記號(hào)，上表可以表示為cAbT可以看到，如果A是一個(gè)主對(duì)角元素均為零的下三角矩陣，那么上表就可以表示一個(gè)顯式Runge Kutta方法。如果A是一個(gè)主對(duì)角線非零的下三角矩陣，相應(yīng)的Runge Kutta方法就是半隱式方法，(5)式等號(hào)左右就含有級(jí)值k1，k2，kr，計(jì)算ki時(shí)要解一個(gè)包含r個(gè)未知量ki的非線性方程組。如果A是滿足Aij(ij)不全為零，則相應(yīng)的方法就是隱式方法。若 系

17、統(tǒng)是m維的，對(duì)于隱式Runge-Kutta方法實(shí)現(xiàn)的每一個(gè)遞推步，均需要求解一個(gè)m×r維的非線性方程組，一般用牛頓迭代法求解。例如01200001200-120162316這是Kutta得到的3級(jí)3階顯式Runge Kutta公式。而012100014140010162316與012100052413-124162316162316分別是3級(jí)4階的半隱式Runge Kutta公式和隱式Runge Kutta公式。要求出具體的Runge Kutta公式，一般有兩種辦法。一種是將(6)在點(diǎn)xi,yi處進(jìn)行泰勒展開，并且代入到(5)中，再與y(xi+h)在xi處的泰勒展開式進(jìn)行對(duì)比，從而確

18、定參數(shù)c，b，A。另一種方法就是將微分方程轉(zhuǎn)化為等價(jià)的積分方程，用數(shù)值積分得到表達(dá)式，再進(jìn)行對(duì)比得到參數(shù)?；诤笠环N方法，可以構(gòu)造得到以下三種不同的隱式Runge-Kutta方法。3.2 Gauss型隱式RK方法17設(shè)c1，c2，cs為Ps2c-1=0的根，這里Ps(t)是0,1上的s次Legendre多項(xiàng)式，0<ci<1，i=1,s。求s級(jí)的2s階的Gauss型Runge Kutta公式的參數(shù)c，b，A需要以下步驟:（1）求出s次的Legendre多項(xiàng)式Ps2c-1=0的s個(gè)零點(diǎn)c1，c2，cs；（2）由線性方程組j=1sbjcjk-1=1k k=1,s確定系數(shù)bj，j=1,s；

19、（3）計(jì)算系數(shù)矩陣A=CW-1在這個(gè)基礎(chǔ)上，就給出了s=1,2,5，p=2s的一系列Gauss型隱式Runge Kutta方法。定理4910 如果一個(gè)數(shù)值方法應(yīng)用于方程y'=y,其中為復(fù)常數(shù)，則對(duì)于所有，Re<0和對(duì)于固定步長h>0，當(dāng)n時(shí)，得到的數(shù)值解趨于零，則稱這種方法是A穩(wěn)定的。定理5910 s級(jí)Gauss型隱式Runge Kutta方法是2s階的，其穩(wěn)定函數(shù)是s,sPade近似且是A穩(wěn)定的。以下列出s=3，p=6的Gauss型隱式Runge Kutta方法之一：12-15101212+151053629-1515536-1530536+152429536-15245

20、36+153029+1515536518495183.3 Radau型隱式RK方法17這種方法的參數(shù)c，b，A需要下面幾個(gè)步驟求得：（1）求多項(xiàng)式Pst-Ps-1(t)的零點(diǎn)c1，c2，cs，并指定c1>0，cs=1；（2）計(jì)算系數(shù)bk，bk=01Pst-Ps-1(t)t-ckPs'ck-Ps-1'(ck)dt k=1,s（3）計(jì)算系數(shù)矩陣A，A=CW-1例如s=3，p=5的RadauA型Runge Kutta方法之一為：4-6104+610188-76360296-16961800-2+36225296+1696180088+76360-2-3622516-63616+

21、6361916-63616+63619定理631 s級(jí)RadauA型Runge Kutta方法和s級(jí)RadauA型Runge Kutta方法是2s-1階的，穩(wěn)定函數(shù)是(s-1,s)次對(duì)角Pade近似。這兩種都是A穩(wěn)定的。3.4 Lobatto型隱式RK方法17這種方法的參數(shù)c，b，A需要下面幾個(gè)步驟求得：（1）求多項(xiàng)式Pst-Ps-2(t)的零點(diǎn)c1，c2，cs，并指定c1=0，cs=1；（2）計(jì)算系數(shù)bk，bk=01Pst-Ps-2(t)t-ckPs'ck-Ps-2'(ck)dt k=1,s（3）計(jì)算系數(shù)矩陣A，A=D-1(WT)-1N-CTD列出4級(jí)6階RadauA型隱式R

22、unge Kutta方法05-5105+5101011+5120025-512011-512011225+135120512025-1351200-1+512025+5120512-1-5120112112512512112定理731 s級(jí)LobattoA，B，C型隱式Runge Kutta方法是2s-2階的，LobattoA和B型Runge Kutta方法的方法的穩(wěn)定函數(shù)是(s-1,s-1)對(duì)角Pade近似，LobattoC型隱式Runge Kutta方法是(s-2,s)次對(duì)角Pade近似。所以這些方法都是A穩(wěn)定的。4隱式RK方法的實(shí)現(xiàn)4.1非線性系統(tǒng)的改進(jìn)對(duì)于s級(jí)隱式隱式Runge Kut

23、ta方法Yi=yn+hj=1saijfxn+cjh,Yj i=1,s (7)yn+1=yn+hj=1sbjfxn+cjh,Yj (8)令 zi=Yi-y (9)于是(7)變?yōu)?zi=j=1saijfxn+cjh,yn+zj (10)當(dāng)?shù)玫降?10)解z1，z1時(shí)，由顯式(8)可以很快得到解yn+1。若隱式Runge Kutta方法的系數(shù)矩陣式非奇異的，那么(10)可以寫為z1zs=Ahf(xn+c1h,yn+z1)hf(xn+csh,yn+zs) (11)所以(8)可以看成與yn+1=yn+i=1sdizi等價(jià)的，其中d1,ds=b1,bsA-1 (12)比如s=3，p=5的RadauA型Ru

24、nge Kutta方法中，d=0,0,1。4.2簡(jiǎn)化的Newton迭代法對(duì)于一般的非線性微分方程，系統(tǒng)(10)必須用迭代的方法求解。本文采用Newton迭代法。Newton迭代法應(yīng)用于系統(tǒng)(10)，每次迭代都需要一個(gè)系數(shù)矩陣I-ha11fyxn+c1h,yn+z1ha1sfyxn+csh,yn+zs-has1fyxn+c1h,yn+z1hassfyxn+csh,yn+zs的線性方程組。定義5 若s階矩陣B=bij，m階矩陣A=aij，且有BA=b11Ab1sAbs1AbssA則稱運(yùn)算是B與A的Kronecker積。為了簡(jiǎn)化(10)的計(jì)算，我們用近似值Jfyxn,yn代替Jacobi矩陣fyxn

25、+cih,yn+zi的值，于是方程(10)的簡(jiǎn)化Newton迭代法變?yōu)镮-hAJZk=-Zk+h(AI)F(Z(k) Z(k+1)=Z(k)+Z(k) (13)這里Z(k+1)=(z1k,zs(k)T是解的第k次近似，Zk=Z1k,Zs(k)T是增量，F(xiàn)(Z(k)是FZk=fxn+c1h,yn+z1(k),fxn+csh,yn+zs(k)T的縮寫。每一次的迭代需要s次由函數(shù)f的求值和一個(gè)NS維線性方程組的求解。因?yàn)榫仃嘔-hAJ對(duì)于所有的迭代是相同的，所以其LU分解在一個(gè)積分步內(nèi)只要進(jìn)行一次，進(jìn)而減少了計(jì)算時(shí)間。5數(shù)值實(shí)驗(yàn)與結(jié)果分析問題：y1'=10-7×y1+sinx y1

26、0=1y2'=10-3×y2+cosx y20=1 寫成矩陣形式就是：y1'y2'=10-70010-3y1y2+sinxcosx很明顯該方程組的剛性比R=10-310-7=1041,因此是個(gè)剛性方程組。取步長為0.1，在區(qū)間0,1內(nèi)分別用四級(jí)四階經(jīng)典顯式Runge Kutta方法yi+1=yi+16k1+2k2+2k3+k4k1=fxi,yi k2=fxi+12h,yi+12hk1 k3=fxi+23h,yi+23hk2 k4=fxi+h,yi+hk3 和二級(jí)四階隱式Runge Kutta方法yi+1=yi+12k1+k2 k1=f(xi+3-36h,yi+

27、14k1h+3-2312k2h)k2=f(xi+3+36h,yi+14k2h+3+2312k1h)進(jìn)行求解。得到結(jié)果如下表一：真值結(jié)果xi y1 y2 0.0000000e+00 1.0000000e+00 1.0000000e+00 1.0000000e-01 1.0049958e+00 1.0999384e+00 2.0000000e-01 1.0199334e+00 1.1988893e+00 3.0000000e-01 1.0446635e+00 1.2958649e+00 4.0000000e-01 1.0789390e+00 1.3898974e+00 5.0000000e-01

28、1.1224175e+00 1.4800481e+00 6.0000000e-01 1.1746644e+00 1.5654174e+00 7.0000000e-01 1.2351579e+00 1.6451531e+00 8.0000000e-01 1.3032934e+00 1.7184598e+00 9.0000000e-01 1.3783901e+00 1.7846058e+00 1.0000000e+00 1.4596978e+00 1.8429313e+00表二：隱式結(jié)果xi y1 y2 0.0000000e+00 1.0000000e+00 1.0000000e+00 1.000

29、0000e-01 1.0049958e+00 1.0999384e+00 2.0000000e-01 1.0199334e+00 1.1988893e+00 3.0000000e-01 1.0446635e+00 1.2958649e+00 4.0000000e-01 1.0789390e+00 1.3898974e+00 5.0000000e-01 1.1224175e+00 1.4800481e+00 6.0000000e-01 1.1746644e+00 1.5654173e+00 7.0000000e-01 1.2351579e+00 1.6451531e+00 8.0000000e-

30、01 1.3032934e+00 1.7184598e+00 9.0000000e-01 1.3783901e+00 1.7846058e+00 1.0000000e+00 1.4596978e+00 1.8429313e+00表三：顯式結(jié)果xi y1 y2 0.0000000e+00 1.0000000e+00 1.0000000e+00 1.0000000e-01 1.0049958e+00 1.0999336e+00 2.0000000e-01 1.0199334e+00 1.1988802e+00 3.0000000e-01 1.0446635e+00 1.2958521e+00 4.

31、0000000e-01 1.0789390e+00 1.3898814e+00 5.0000000e-01 1.1224175e+00 1.4800297e+00 6.0000000e-01 1.1746645e+00 1.5653972e+00 7.0000000e-01 1.2351579e+00 1.6451319e+00 8.0000000e-01 1.3032934e+00 1.7184382e+00 9.0000000e-01 1.3783901e+00 1.7845846e+00 1.0000000e+00 1.4596978e+00 1.8429111e+00很明顯的看到，在保

32、留八位小數(shù)的情況下，二級(jí)四階隱式Runge Kutta方法與真值無異，能夠精確到小數(shù)點(diǎn)后七位以上。而經(jīng)典Runge Kutta只能精確到小數(shù)點(diǎn)后四位。因此對(duì)于剛性方程組，相同步長下，隱式Runge Kutta方法的精度比顯式Runge Kutta方法高得多。并且由于步長小，在實(shí)際的求解區(qū)間里面，涉及的遞推次數(shù)少，從而函數(shù)的計(jì)算量就小。參考文獻(xiàn)1 費(fèi)景高. 并行顯式 RungeKutta 公式的實(shí)現(xiàn)J. 計(jì)算機(jī)工程與設(shè)計(jì), 1994 (5): 34-40.2 Enenkel R F, Jackson K R. DIMSEMs-diagonally implicit single-eigenval

33、ue methods for the numerical solution of stiff ODEs on parallel computersJ. Advances in Computational Mathematics, 1997, 7(1-2): 97-133.3 張誠堅(jiān), 余紅兵. 非線性延遲系統(tǒng)的一類并行預(yù)校算法J. 華中理工大學(xué)學(xué)報(bào), 2000, 28(11): 104-106.4 李愛雄, 劉偉豐, 張誠堅(jiān), 等. 求解 DDEs 的多導(dǎo)龍格庫塔方法的漸近穩(wěn)定性 Asymptotic stability of multi-derivative Runge-Kutta meth

34、od for delay differential equationsJ. 華中科技大學(xué)學(xué)報(bào) (自然科學(xué)版), 2002, 6: 108-110.5 廖文遠(yuǎn), 李慶揚(yáng). 一類求解剛性常微分方程的半隱式多步 RK 方法J. 清華大學(xué)學(xué)報(bào): 自然科學(xué)版, 1999, 39(6): 1-4.6 龐立君, 朱永忠. 類隨機(jī)微分方程 Runge. Kutta 方法的指數(shù)穩(wěn)定性J. 河海大學(xué)學(xué)報(bào) (自然 科學(xué)版), 2008, 36(3).7 Curtiss C F, Hirschfelder J O. Integration of stiff equationsJ. Proc. Nat. Acad. S

35、ci, 1952, 38(235): 1.8 Gear C W. 常微分方程初值問題的數(shù)值解法J. 1978.9 Butcher J C. Implicit runge-kutta processesJ. Mathematics of Computation, 1964, 18(85): 50-64.10 Ehle B L. High order A-stable methods for the numerical solution of systems of DE'sJ. BIT Numerical Mathematics, 1968, 8(4): 276-278.11 Burrag

36、e K, Butcher J C, Chipman F H. An implementation of singly-implicit Runge-Kutta methodsJ. BIT Numerical Mathematics, 1980, 20(3): 326-340.12 Shampine L F. Implementation of implicit formulas for the solution of ODEsJ. SIAM Journal on Scientific and Statistical Computing, 1980, 1(1): 103-118.13 Lindb

37、erg B. On smoothing and extrapolation for the trapezoidal ruleJ. BIT Numerical Mathematics, 1971, 11(1): 29-52.14 Lindberg B. IMPEX 2-a procedure for solution of systems of stiff differential equationsR. CM-P00069411, 1973.15 Bader G, Deuflhard P. A semi-implicit mid-point rule for stiff systems of

38、ordinary differential equationsJ. Numerische Mathematik, 1983, 41(3): 373-398.16 孫志忠, 袁慰平, 聞?wù)鸪? 數(shù)值分析M. 東南大學(xué)出版社, 2002.17 劉德貴, , 費(fèi)景高. 剛性大系統(tǒng)數(shù)字仿真方法M. 河南科學(xué)技術(shù)出版社, 1996.18 Cockburn B, Shu C W. The RungeKutta discontinuous Galerkin method for conservation laws V: multidimensional systemsJ. Journal of Comput

39、ational Physics, 1998, 141(2): 199-224.19 Monovasilis T, Kalogiratou Z, Simos T E. Exponentially fitted symplectic rungeKuttaNyström methodsJ. Appl. Math. Inf. Sci, 2013, 7(1): 81-85.20 Hochbruck M, Pazur T. Implicit Runge-Kutta Methods and Discontinuous Galerkin Discretizations for Linear Maxw

40、ell's EquationsJ. SIAM Journal on Numerical Analysis, 2015, 53(1): 485-507.21 Papadopoulos D F, Simos T E. A modified Runge-Kutta-Nyström method by using phase lag properties for the numerical solution of orbital problemsJ. Applied Mathematics & Information Sciences, 2013, 7(2): 433-437

41、.22 Zhong X, Shu C W. A simple weighted essentially nonoscillatory limiter for RungeKutta discontinuous Galerkin methodsJ. Journal of Computational Physics, 2013, 232(1): 397-415.23 Lambert J D, Lambert J D. Computational methods in ordinary differential equationsM. London: Wiley, 1973.24 Zhu J, Zho

42、ng X, Shu C W, et al. RungeKutta discontinuous Galerkin method using a new type of WENO limiters on unstructured meshesJ. Journal of Computational Physics, 2013, 248: 200-220.25 Jayakumar T, Maheshkumar D, Kanagarajan K. Numerical solution of fuzzy differential equations by Runge-Kutta method of ord

43、er fiveJ. International Journal of Applied Mathematical Science, 2012, 6: 2989-3002.26 Dimarco G, Pareschi L. Asymptotic preserving implicit-explicit Runge-Kutta methods for nonlinear kinetic equationsJ. SIAM Journal on Numerical Analysis, 2013, 51(2): 1064-1087.27 Miranker W L, Miranker A. Numerica

44、l Methods for Stiff Equations: And Singular Perturbation ProblemsM. Springer Science & Business Media, 2001.28 Hadi M, Anderson M, Husein A. Using 4th order Runge-Kutta method for solving a twisted Skyrme string equationC/THE 4TH INTERNATIONAL CONFERENCE ON THEORETICAL AND APPLIED PHYSICS (ICTAP

45、) 2014. AIP Publishing, 2016, 1719: 030005.29 Jimenez J C, Sotolongo A, Sanchez-Bornot J M. Locally linearized runge kutta method of dormand and princeJ. Applied Mathematics and Computation, 2014, 247: 589-606.30 Jimenez J C, Sotolongo A, Sanchez-Bornot J M. Locally linearized runge kutta method of

46、dormand and princeJ. Applied Mathematics and Computation, 2014, 247: 589-606.31 Wanner G, Hairer E. Solving ordinary differential equations IIM. Springer-Verlag, Berlin, 1991.附錄緒論程序：a=0;b=1; %求解區(qū)間f1=(t,x,y)(-0.01*x-99.99*y);f2=(t,x,y)(-100*y);h=0.0001; %步長T=zeros(1,(b-a)/h)+1);X=zeros(1,(b-a)/h)+1);

47、Y=zeros(1,(b-a)/h)+1);T=a:h:b;X(1)=2; %初值Y(1)=1;for j=1:(b-a)/h k11=feval(f1,T(j),X(j),Y(j); k12=feval(f2,T(j),X(j),Y(j); k21=feval(f1,T(j)+h/2,X(j)+h/2*k11,Y(j)+h/2*k12); k22=feval(f2,T(j)+h/2,X(j)+h/2*k11,Y(j)+h/2*k12); k31=feval(f1,T(j)+h/2,X(j)+h/2*k21,Y(j)+h/2*k22); k32=feval(f2,T(j)+h/2,X(j)+h

48、/2*k21,Y(j)+h/2*k22); k41=feval(f1,T(j)+h,X(j)+h*k31,Y(j)+h*k32); k42=feval(f2,T(j)+h,X(j)+h*k31,Y(j)+h*k32); X(j+1)=X(j)+h*(k11+2*k21+2*k31+k41)/6; Y(j+1)=Y(j)+h*(k12+2*k22+2*k32+k42)/6; R=T' X' Y' endsave 0.0001結(jié)果.txt R ascii真值程序：y1=exp(-100)+exp(-0.01);y2=exp(-100);R=1 y1 y2;save 真值結(jié)果.txt R ascii第五章程序：syms x;a=0;b=1; h=0.1; f1=dsolve('Dy=0.0000001*y+sin(x)','y(0)=1','x');f2=dsolve('Dy=0.001*y+cos(x)','y(0)=1','x');X=a:h:b;Y1=double(subs(f1,x,X);Y2=double(subs(f2,x,X);R1=X' Y1' Y2'xlswrite('結(jié)果',R1,'真值')f