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

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

2、Runge-Kutta方法單步方法Newto"迭代法AbstractTliispapermainlyintroducestheImplicitRunge-KuttaMethodsandasimpledesaiptionofGaussimplicitRunge-Kuttamethod,RadauimplicitRunge-KuttamethodandLobattoimplicitRunge-Kuttamethod.TliesebasicIinplicitRunge-Kuttamethodsareusedtosolvethestiffequations.TlieseimplicitRun

3、ge-KuttamethodsiarecomparedwithtlieclassicalexplicitRunge-KuttamethodTliispaperexplaintlieadvantagesanddisadvantagesofthetwokindofnumericalmethods.Keywords:StiffsystemliiiplicitRuiige-KuttamethodOnestepmethodNewtoniterativemethod目錄1 .緒論11.1 剛性方程11.2 隱式RK方法的研究意義21.3 RK方法的研究現(xiàn)狀32 .單步RK方法的收斂性和穩(wěn)定性52.1 單步

4、RK方法的收斂性52.2 單步RK方法的穩(wěn)定性63 .三類隱式RK方法83.1 弓I言83.2 Gauss型隱式RK方法93.3 Radau型隱式RK方法103.4 Lobatto型隱式RK方法114隱式RK方法的實現(xiàn)1336 非線性系統(tǒng)的改進(jìn)1336 簡化的Newton迭代法135數(shù)值實驗與結(jié)果分析15參考文獻(xiàn)18附錄21.緒論剛性方程對于一般的線性常系數(shù)系統(tǒng)y'=4y+3(t)4為mxm的矩陣，特征值為Z(i=1,2,，m)。定義13若一個系統(tǒng)滿足RetV。,i=1,2,mmaxiReXiZmiTii|Z?eAJ=R»1其中R為剛性比，則這個系統(tǒng)稱為剛性系統(tǒng)。定義ZB1若

5、線性系統(tǒng)y'=AyxG0,T或非線性系統(tǒng)y'=f(x,y)xwo,T的矩陣4或Jacobi矩陣。/dy的特征值人滿足maxReAt»1則其是剛性的。定理1(解的存在性與唯一性)(1)對于所有(y)wD,函數(shù)f(x,y)是連續(xù)的：(2)對于任何Q,y),(x,y*)WD,存在常數(shù)L,是函數(shù)滿足-faWly-y*ll則初值問題(y=/Uy)(a<x<b)ly(a)=rj有唯一解。其中y=O"2,),，D=(x,y)|a<x<b,-co<a<b<8；8<yt<a,Z=1,2,，力。其中N被稱為/必*4/2常數(shù)定

6、義3如果一個常微分系統(tǒng)的Lipschitz常數(shù)L很大(大于20),則它是剛性的。隱式RK方法的研究意義在常微分方程及常微分方程組的數(shù)值解法中，Runge-Kats方法是目前應(yīng)用最為廣泛的數(shù)值解法之一，同時乂具有誤差小，精度高的特點。盡管顯式Hunge-Kutta方法能夠非常準(zhǔn)確、快速的給出大部分常微分方程組的數(shù)值解。但是在化學(xué)、自動控制電力系統(tǒng)等領(lǐng)域中，會出現(xiàn)一些病態(tài)的常微分方程組，也就是剛性方程組。剛性方程組對于數(shù)值解法的穩(wěn)定性要求苛刻，比如方程組(yi=-O.Olyi-99.99為，71(0)=2ly2=-100y2,y2(0)=1將其表示為矩陣形式:令yiyf2-o.oio-99.991

7、ryq-100電1-0.01-99.991剛性比s=M=10000I人21-0-100J發(fā)現(xiàn)4特征值為：而=-100,A2=-0.01,方程組的解為：)=e100x+e_001xy2(xl=je-100快臟態(tài)慢瞬態(tài)解由快瞬態(tài)和慢瞬態(tài)兩部分構(gòu)成。由于慢瞬態(tài)的部分，yi(x)衰減變得十分緩慢。當(dāng)口變量變到x=391時，函數(shù)值還未下降到初值的1%，求解區(qū)間至少取為(0,391)。另一方面，由于快瞬態(tài)的部分，以。)衰減的非?？?，因此步長要取得非常小。從絕對穩(wěn)定性的方面來看，如果用四階顯式經(jīng)典RK方法求解，絕對穩(wěn)定區(qū)間要求放”-2.78,0),則要求h<0.0278o這樣，在(0,391)上就要計

8、算14065步，計算量巨大，因此計算區(qū)間(0,391)內(nèi)的解時，舍入誤差積累會特別嚴(yán)重。例如取求解區(qū)間為0,1,用不同步長九來計算力(1)和及(1)的值。利用四階顯式經(jīng)典RK方法求解如下：h0.04yi(i)y2(i)0.022.9802322e+179.9004983e-012.9802322e+17l.3929556e-240.019.9004983e-012.53003646430.0019.9004983e-013.7204130e-440.00019.9004983&013.7200760e-44真值9.9004983e-013.7200760e-44很顯然，保留八位有效數(shù)字

9、的情況下，要保持良好的精度，步長要取得非常小，這就增加了計算量。而隨著步長的加大，誤差也會越來越大。當(dāng)步長加大到絕對穩(wěn)定與之外時（即力0.0278）,計算結(jié)果就完全不可信了。對于剛性方程組，顯式方法已經(jīng)遠(yuǎn)遠(yuǎn)不能勝任了，一般采用絕對穩(wěn)定性更好的方法（如隱式Hunge-Kutta方法）進(jìn)行求解，本文采用單步隱式Rimge-Kutm方法，而對于陶式方法中的級值得求解，本文采用Newto/i迭代法。RK方法的研究現(xiàn)狀研究基于標(biāo)準(zhǔn)模型方程的Riwge-Kutta方法的常見形式為:yt+i=yt+h勺；=i,的=/(“)JTkj=f(xt+/1也+h'ujiki)(J=2,3,r)(顯式)rM+i

10、=M+九2%勺，r角%=+入彼,+h'內(nèi)的)(j=1,2,3,r)i=i（隱式）因為Rimge-Kutta方法是比較成熟的常微分方程數(shù)值解法。所以如今主要是對于經(jīng)典的Runge-Kutta方法進(jìn)行完善和擴充，在一定的條件下，提高級數(shù)以提高精度?；蛘呤菍imge-Kuaa方法與某些領(lǐng)城結(jié)合使用。在1994年，費景高給出了一種顯式Kunge-Kutta并行方法，從而實現(xiàn)Hunge-Kutta方法在多處理機上的應(yīng)用。1997年，Eneiikel和Jackson實現(xiàn)了Rimge-方法的對角隱式并行改進(jìn)。1999年，廖文遠(yuǎn)和李慶揚給出了一類求解剛性常微分方程的半隱式多步Mmge-Kutm方法。

11、2000年，張誠堅和余洪兵針對非線性延遲系統(tǒng)構(gòu)造了一類并行預(yù)校算法，給出其算法的局部誤差估計，數(shù)值實驗表明該算法是有效的，且具一定的可比性。2003年，李愛雄通過對傳統(tǒng)單支方法的計算格式進(jìn)行改造，得到了解DDEs的兩類單支并行算法單支并行預(yù)校算法和單支并行算法，并對方法的收斂性和穩(wěn)定性做出了分析。2008年，龐麗君和朱永忠給出了一類隨機微分方程Runge-Kutta方法的指數(shù)穩(wěn)定性。.單步RK方法的收斂性和穩(wěn)定性單步RK方法的收斂性對于常微分初值問題(y=/(y)(a<%<b)ly(a)=r的單步顯式公式仍+i=%+h(p&,y"九)(i=o,l,-,n-1)H部

12、截斷誤差可以表示為&+1=yto+i)-tyM+九3(%,見，九)(i=o,L-1)(3)定理2161：設(shè)y(x)為Q)的解，殖之。為(2)的解。如果：(1)存在常數(shù)c,使得因+11Wc/iP+i(i=0,l,2,.,n-l)(2)存在q>0,使得0<p(x,y,/i)max<L(x»)皿dyQ<a<h其中。0=(x,y)|as入,OyQ)。，y,y(x)十。記c。=1，則當(dāng)hWm加卜時，有Eg)<cQhP證:由(3)得y(x(+i)=y(%i)+廂(孫yQ)九)+4+iG=0,L,汽-1)(4)將(4)與(2)相減y®+i)-y

13、(%()=y®)-+九3(卬丫(%)九)一+Ri+i(i=。，1，n1)由yGo)%=。知道，在i=o時，|y(xj-yjwJ臚成立?，F(xiàn)在假設(shè)0<i<k-1時也是成立的。在hWP區(qū)時有：Coly()-yilJJ=a(i=0,1,-,n-1)進(jìn)而即(孫y®)，九)一W(MM，h)ld(p%,h)/、=一一(yGr)wL|y(/一%)l0<i<k-i其中q是y(%)和yi之間的數(shù)。于是定理結(jié)合條件與(4)式，可以得到1丫(修+1)-%+11w|y(Xi)-yd+h|<p(%i，y(%i),A)-<p(Xj,yi,h)|+|/?l+1|Mly(

14、x()-y(|+Lhy(xt)-yt+c/ip+1=(1+L/i)|y(xt)y114-chP+i0<i<k1從而.(1+Lh)k1lyQk-%)1V(1+Lh)kD-%l+7-ch1(1+Ln)1因為y(xo)-y0=。及(1+Lh»<eLkhWeb-a),得到lv(4)-ykyeL(b'a)-lhp=cohpL所以當(dāng)i=%時|曠(即)一yjWc()/iP也是成立的。(證畢)對于單步顯式格式而言，當(dāng)fQ,y)和竽kED。內(nèi)連續(xù)時，那么定理1的條件(2)總是滿足的。所以滿足上述條件的單步顯式公式，只要也滿足相容性條件w(3(x),0)=/Q,yQ)那么它一定

15、是收斂的。單步RK方法的穩(wěn)定性定義4對于初值問題(1),若%乜0是由(2)得到，區(qū)?=0是下面擾動問題的解:仔+1=4+九/(，/，九)+(SH1(i=0,1,2,，九一1)Izo=+6。如果存在正常數(shù)C,%及九0,使所有的£G(0,2,hE(0,hQ,當(dāng)maxoq夕141<£時,有maxIVizA<Cs0<i<n就稱(2)是穩(wěn)定的或者零穩(wěn)定的。定理3n司在定理1的條件下，單步顯式公式(2)是穩(wěn)定的。.三類隱式RK方法引言對于初值問題Q)，隱式Runge-Kutta方法的一般形式為ryt+1=yi+hbjkj(5)j=1rkj=f(勺+Cihtyt+

16、引的)(j=1,2,3,r)(6)其中，距=而+訪，n=0,12,為數(shù)軸上的離散點列;h為步長,%為解y(%)的近似值:cvc2».»稱為隱式/hinge-KSa方法的步長;bvb2,.»%為權(quán)系數(shù)。令4=(即)，j=l,2,3,.r。稱4為系數(shù)矩陣，因此上式可以簡化表示為C1。1?bi.br使用Btche廠陣的記號，上表可以表示為可以看到，如果4是一個主對角元素均為零的下三角矩陣，那么上表就可以表示一個顯式Rimge-Kaa方法。如果4是一個主對角線非零的下三角矩陣，相應(yīng)的Rzmge-Katta方法就是半隱式方法，(5)式等號左右就含有級值的，依，%,計算加時要

17、解一個包含r個未知量匕的非線性方程組。如果4是滿足4)(iWj)不全為零，則相應(yīng)的方法就是隱式方法。若系統(tǒng)是m維的，對于隱式Runge-Kutta方法實現(xiàn)的每一個遞推步，均需要求解一個mxr維的非線性方程組，一般用牛頓迭代法求解。例如012024160001-002-120121636這是Kutta得到的3級3階顯式Rimge-Kutta公式。而001144011263000151224311263121636分別是3級4階的半隱式Rioige-Kutta公式和隱式Runge-Kutta公式。要求出具體的Rimge-Kutta公式，一般有兩種辦法。一種是將（6）在點（環(huán)）處進(jìn)行泰勒展開，并且代

18、入到（5）中，再與yaIh）在/處的泰勒展開式進(jìn)行對比，從而確定參數(shù)c,44。另一種方法就是將微分方程轉(zhuǎn)化為等價的積分方程，用數(shù)值積分得到表達(dá)式，再進(jìn)行對比得到參數(shù)?；诤笠环N方法，可以構(gòu)造得到以下三種不同的隱式Runge-Kutta方法。3.2Gauss型隱式RK方法設(shè)C1，。2，Cs為己（2c-1）=0的根，這里鳥是0,1上的s次Legendre多項式，OVqV1,i=l,，S。求s級的2s階的Gauss型Range-Kiztta公式的參數(shù)c,b,4需要以下步驟：（1）求出s次的Legend廠e多項式a（2c1）=0的s個零點q,c2».»cs；（2）由線性方程組工4$

19、-1=1/kk=1,s自i確定系數(shù)與，j=L,s：(3)計算系數(shù)矩陣A=CW'1在這個基礎(chǔ)上，就給出了s=1,2,5,p=2s的一系列Gauss型隱式-Kutta方法。定理4【刈如果一個數(shù)值方法應(yīng)用于方程y'=4y,其中A為更常數(shù)，則對于所有A,Re/lV0和對于固定步長九>0,當(dāng)ms時，得到的數(shù)值解趨于零，則稱這種方法是力穩(wěn)定的。定理5網(wǎng)""s級Gaass型隱式Runge-Kutta方法是2s階的，其穩(wěn)定函數(shù)是(s,s)Pad近似且是4穩(wěn)定的。1V152-0"121V152+0以卜,列出s=3,p=6的Gauss型隱式Range-Kutta

20、方法之一:52VK5VK369-TT36305vB25715"4,3624936245V152V155+-+36309153654518918Radau型隱式RK方法這種方法的參數(shù)c,從4需要下面幾個步驟求得:(1)求多項式R(t)Ps_i(t)的零點Ci，Q，Cs，并指定Cl>0,Cs=l：(2)計算系數(shù)比,I-/y£Jo("CQE'Q)-*©)(3)計算系數(shù)矩陣4,A=CIV-180016-V636-例如s=3,p=5的Radau/ASlRunge-Kiztta方法之一為:4-V6104+y/6101一2+3仍-225-2-3V6-22

21、5-1988-776-360296+169點296-169后180088+7點-360-16+V616-V616+V613636-9定理6311s級Radau/4型Runge-Kutta方法和s級RadauHA型Rmge-Kiztta方法是2s-1階的，穩(wěn)定函數(shù)是（s-l,s）次對角Pad近似。這兩種都是4穩(wěn)定的。Lobatto型隱式RK方法這種方法的參數(shù)c,b,4需要下面幾個步驟求得：（1）求多項式R（t）-尸$一2«）的零點C1，。2，Cs，并指定Ci=0,Cs=1：（2）計算系數(shù)比，廣G（t）_Ps_2（t）以=Jo«M）因3）？一2©）"&quo

22、t;=（3）計算系數(shù)矩陣/,A=DT(m)T(N-C)T。列出4級6階Rada”/4型隱式Runge-Kutta方法05yS105+V5101000011+乃25-V525-13x/5-14-V512012012012011-V525+13425+遍-1-V51201201201201551訪H訪運1551定理7叫級Lobatto川A,川B,川C型隱式Runge-K"q方法是2s-2階的，LobattolHA和川B型.Runge-Kutta方法的方法的穩(wěn)定函數(shù)是(s-l,s1)對角Pad近似，LobattoU/C型隱式Runge-Kutta方法是(s2,s)次對角PadG近似。所以這

23、些方法都是4穩(wěn)定的。4隱式RK方法的實現(xiàn)4.1非線性系統(tǒng)的改進(jìn)對于s級隱式隱式Ruzige-Kutta方法s匕=yn+hWf&+q九普)J=1syn+1=ynth'bjfGn十gh2)j=l(7)Zl=Yt-y(9)于是(7)變?yōu)?=W5jf(xn+Cjhfyn+Zj)六1Q0)當(dāng)?shù)玫降腝0)解Zi，Zi時，由顯式(8)可以很快得到解為+1。若隱式Range-Kutta方法的系數(shù)矩陣式非奇異的，那么(10)可以寫為降1=A/if(/+csh9ynzs)(ID所以(8)可以看成與%+i=yn+244i=l等價的，其中31,4)=Si，,瓦)4T(12)比如s=3,p=5的Rada

24、”HARunge-Kutta方法中，d=(0,0,1)4.2簡化的Newton迭代法對于一般的非線性微分方程，系統(tǒng)Q0)必須用迭代的方法求解。本文采用Newto九迭代法。Newtori迭代法應(yīng)用于系統(tǒng)(10),每次迭代都需要一個系數(shù)矩陣dfdfI-han(xn+Cihtyn4-zjhals(xn+cshfyn+zs)*dfdfhasi(xn+cYhfyn+zi)-/iass(%n+cshfyn+zs)dydy的線性方程組。定義5若s階矩陣8=(加m階矩陣4=(%),且有。1遇比54B®A=Ji4遇bssA.則稱運算8是B與4的Kronecker積。為了簡化(10)的計算，我們用近似值

25、9?(如加iJacobidf/dy(xn+cth,yn+%)的值，于是方程Q0)的簡化Newt。幾迭代法變?yōu)?Z-hA0/)AZ=-Z")+h(A0/)F(Z)Z(k+D=Z(")+AZ(13)這里Z(k+1)=(zf),29)7是解的第k次近似，AZg=(Z，)(Z,)y是增量，/?(2的)是(2(切)=。(5”1勺兒叫Izj"),"(如I/,為Iz$)的縮寫。每一次的迭代需要s次由函數(shù)r的求值和一個NS維線性方程組的求解。因為矩陣(/M/)對于所有的迭代是相同的，所以其LU分解在一個積分步內(nèi)只要進(jìn)行一次，進(jìn)而減少了計算時間。5數(shù)值實驗與結(jié)果分析問題：

26、(y=10-7Xy】+sinWy/O)=1y2=103xy2+cos(x)y2(O)=1寫成矩陣形式就是：Mrio-70iRi+QMx)yj-1oi。一3八)JcosM)很明顯該方程組的剛性比R=J=104»1,因此是個剛性方程組。取步長為0.1,在區(qū)間0,1內(nèi)分別用四級四階經(jīng)典顯式Rinige-Kutta方法1%+1=%+工*1+2A2+2k3+”4)o的=fd%)<6=f("),+加i)B=f(即+|九%+|h&)土4=f(%+htyt4-hk3)和二級四階隱式Rimge-Kutta方法1%+1=%+萬(的+攵2)3-x/313-2>/3，的=f(M

27、+72+1的九十行一心九)o41Z+V313+2V32=f(左+72%+加入+kih)O1Z進(jìn)行求解。得到結(jié)果如下表一：真值結(jié)果%加20.0000000e+001.0000000e+001.0000000e+001.00000000-011.0049958P+001.0999384P+002.0000000e-011.0199334e+001.1988893e+001.0446635e+001.2958649e+003.0000000e-014.0000000e-01S.OOOOOOOe-Ol6.0000000e-017.0000000e-018.0000000e-019.0000000e-0

28、11.0000000e+001.0789390e+001.1224175e+001.1746644e+001.2351579e+001.3032934e+001.3783901e+001.4596978e+001.3898974e+001.4800481e+001.5654174e+001.6451531e+001.7184598e+001.7846058e+001.8429313e+00表二:隱式結(jié)果0.0000000e+00l.OOOOOOOe-Ol2.0000000e-013.0000000e-014.0000000e-015.0000000e-016.0000000e-017.0000

29、000e-018.0000000e-019.0000000e-01l.OOOOOOOe+OO1.0000000e+001.0049958e+001.0199334e+001.0446635e+001.0789390e+001.1224175e+001.1746644e+001.2351579e+001.3032934e+001.3783901e+001.4596978e+001.0000000e+001.0999384e+001.1988893e+001.2958649e+001.3898974e+001.4800481e+001.5654173e+001.6451531C+001.71845

30、98e+001.7846058e+001.8429313e+00表三：顯式結(jié)果0.0000000e+001,0000000e+001.0000000e+001.0000000e-011.0049958e+001.0999336e+00Z0000000e-011.0199334e+001.1988802e+003.0000000e-011.0446635e+001.2958521e+004.0000000e-011.0789390e+001.3898814e+005.0000000e-011.1224175e+001.4800297e+006.0000000e-011.1746645e+001.

31、5653972e+007.0000000e-011.2351579e+001.6451319e+008.0000000e-011.3032934e+001.7184382e+009.0000000e-011.3783901e+001.7845846e+001.0000000e+001.4596978e+001.8429111e+00很明顯的看到，在保留八位小數(shù)的情況下，二級四階隱式Rimge-Kutta方法與真值無異，能夠精確到小數(shù)點后七位以上。而經(jīng)典Rimge-Kutta只能精確到小數(shù)點后四位。因此對于剛性方程組，相同步長下，隱式Range-Kutta方法的精度比顯式Rimge-Kia”方法

32、高得多。并且由于步長小，在實際的求解區(qū)間里面，涉及的遞推次數(shù)少，從而函數(shù)的計算量就小。參考文獻(xiàn)囚費景高,并行顯式Runge-Kutta公式的實現(xiàn)J.計算機工程與設(shè)計，1994(5):34-40.2EneiikelRF,JacksonKR.DIMSEMs-diagonallyimplicitsingle-eigenvaluemethodsforthemunericalsolutionofstiffODEsonparallelcomputersJ.AdvancesinComputationalMathematics,1997,7(1-2):97-133.3張誠堅，余紅兵.非線性延遲系統(tǒng)的一類并行預(yù)

33、校算法J.華中理工大學(xué)學(xué)報,2000,28(11):104-106.4李愛雄，劉偉豐，張誠堅，等.求解DDEs的多導(dǎo)龍格庫塔方法的漸近穩(wěn)定性Asymptoticstabilityofmulti-derivativeRmige-KuttamethodfordelaydifferentialequationsJ.華中科技大學(xué)學(xué)報(自然科學(xué)版),2002,6:108T10.5廖文遠(yuǎn)，李慶揚.一類求解剛性常微分方程的半隱式多步RK方法J.清華大學(xué)學(xué)報：自然科學(xué)版,1999,39(6):1-4.龐立君，朱永忠.類隨機微分方程Runge.Kutta方法的指數(shù)穩(wěn)定性J.河海大學(xué)學(xué)報(自然科學(xué)版),2008,

34、36(3).7CiirtissCF,HirsclifelderJO.IiitegiationofstiffeqiiationsJ.Proc.Nat.Acad.Sci,1952,38(235):1.8GearCW.常微分方程初值問題的數(shù)值解法J.1978.ButcherJC.Implicitninge-kuttaprocessesJ.MathematicsofComputation,1964,18(85):50-64.10ElileBL.HigliorderA-stablemethodsforthenumencalsolutionofsystemsofDEfsJBITNumericalMathe

35、matics,1968,8(4):276-278.llBuiTageK,ButcherJC,CliipmanFH.Animplementationofsingly-implicitRunge-KuttamethodsJ.BITNumericalMathematics,1980.20(3):326-340.12ShampineLF.IiiiplementatioiiofimplicitfonnulasfortliesolutionofODEsJ.SIAMJounialonScientificandStatisticalComputing,1980,1(1):103-118.LindbergB.O

36、nsmootliingandextiapolationforthetrapezoidalmleJ.BITNumericalMathematics,1971,11(1):29-52.LindbergB.IMPEX2-aprocedureforsolutionofsystemsofstiffdifferential18equationsR.CM-P000694U,1973.BaderG.DeuflhardP.Asemi-implicitmid-pointmleforstiffsystemsofordinarydifferentialeqiiationsJ.NumerischeMathematik,

37、1983,41(3):373-398.16孫志忠，袁慰平，聞?wù)鸪?數(shù)值分析M.東南大學(xué)出版社,2002.17劉德貴,，費景高.剛性大系統(tǒng)數(shù)字仿真方法M.河南科學(xué)技術(shù)出版社,1996.CockburnB,ShuCW.TlieRunge-KuttadiscontinuousGalerkinmethodforconseivationlawsV:multidimensionalsystemsJ.JournalofComputationalPhysics,1998,141(2):199-224.MonovasilisT,KalogiratouZ,SimosTE.Exponentiallyfittedsy

38、inplecticiunge-Kutta-NystidmmethodsJ.Appl.Math.hif.Sci,2013,7(1):81-85.HochbnickM,PazurT.ImplicitRunge-KuttaMethodsandDiscontinuousGalerkinDiscretizationsforLinearMaxwell'sEqiiationsJ.SIAMJournalonNmnericalAnalysis,2015,53(1):485-507.PapadopoulosDF,SimosTE.AmodifiedRunge-Kutta-Nystiommethodbyusi

39、ngphaselagpropeitiesfortlienumericalsolutionoforbitalproblemsJ.AppliedMathematics&InfonnationSciences,2013,7(2):433-437.ZhongX,ShuCW.AsimpleweiglitedessentiallynonoscillatoiylimiterforRungeKuttadiscontinuousGalerkinmethodsJ.JournalofComputationalPhysics,2013,232(1):397-415.LambertJD,LambertJD.Co

40、mputationalmethodsinordinarydifferentialequationsM.London:Wiley,1973.ZhuJ,ZhongX,ShuCW,etal.Runge-KuttadiscontinuousGalerkinmethodusinganewtypeofWENOlimitersonunstnichiredmeshesJ.JounialofComputationalPhysics,2013,248:200-220.Jayakumar工MaheshkiimarD,KanagarajanK.Numericalsolutionoffiizzydifferential

41、equationsbyRunge-KuttamethodoforderfiveJ.IiitemationalJounialofAppliedMathematicalScience,2012,6:2989-3002.DimarcoG,PareschiL.Asymptoticpreservingimplicit-explicitRmige-KuttamethodsfornonlinearkineticequationsJ.SIAMJournalonNumericalAnalysis,2013,51(2):1064-1087.MirankerWL,MiraiikerA.NumericalMethod

42、sforStiffEquations:AndSingularPerturbationProblemsM.SpringerScience&BusinessMedia,2001.HadiM,AndersonM.HuseinA.Using4thorderRunge-KuttamethodforsolvingatwistedSkynnestringeqiiationC/THE4THINTERNATIONALCONFERENCEONTHEORETICALANDAPPLIEDPHYSICS(ICTAP)2014.AIPPublisliing,2016,1719:030005.JimenezJC,S

43、otolongoA,Sanchez-BomotJM.LocallylinearizedrungekiittamethodofdonnandandprinceJ.AppliedMathematicsandComputation,2014,247:589-606.JimenezJC,SotolongoA,Sanchez-BomotJM.LocallylinearizedrungekiittamethodofdonnandandprinceJ.AppliedMathematicsandComputation,2014,247:589-606.31WaimerG,HaiierE.Solvingoidi

44、narydifferentialequationsnM.Springer-Verlag,Berlin,1991.附錄緒論程序：a=O;b=l;%求解區(qū)間fl=(tc,y)(-0.01*x-99.99*y);f2=(t,y)(-100*y);h=0.0001;%)步長T=zeros(l,(b-a)/li)+l);X=zeros(l,(b-a)/h)+l);Y=zeros(l,(b-a)/h)+l);T=a:h:b;X(l)=2;%初值Y(l)=l;forj=l:(b-a)/hkll=feval(fl,TG)(j),Y(j);kl2=feval(f2,T(j),X(j),Y(j);k21=feva

45、l(fl,T(j)+h/2,X(j)+h/2*kll,Y(j)-Hi/2*kl2);k22=feval(f2,T(j)+h/2,X(j)+h/2*kll,Y也2*1<12);k31=feval(fl,T(j)+h/2,X(j)-Hi/2*k21,Yg)+h/2*k22);k32=feval(f2,TO)-Hi/2,xg)-Hi/2*k2i,Y0)+h/2*k22);k4l=feval(fl,T(j)+h,X(j)+h*k31,Y(j)+h*k32);k42=feval(f2,TG)+h,X(j)+h*k31,Y(j)+h*k32);X(j+l)=X(j)-Hi*(kll+2*k21+2*

46、k31+k41)/6;Y(j+l)=Y(j)-Hi*(kl2+2*k22+2*k32+k42)/6;R=TXfY1;endsave0.0001結(jié)果.txtR-ascii真值程序：yl=exp(-l00)+exp(-0.01);y2=exp(-100);R=iyiy2;save真值結(jié)果.txtR-ascii第五章程序：syinsx;a=0;b=l;h=0.1;fl=dsolve(Dy=0.00000014：y+sin(x)','y(0)=l,x,)f2=dsolve(Dy=0.001*y+cos(xy,y(0)=l；|x);X=a:h:b;Yl=doiible(subs(fl,x

47、,X);Y2=double(subs(f2,x,X);Rl=pCYY2;xlswriteC結(jié)果JR"真值)fonnatlongfl=(x,yl)(0.0000001*yl-l-sin(x);f2=(x,y2)(0.001*y2+cos(x);a=0;b=l;11=0.1;X=zeros(l,(ba)/h+l)；Yl=zeros(l,(b-a)/h+l);Y2=zeros(l,(b-a)/h+l);X=a:h:b;Yl(l)=l；Y2(l)=l;syinsklk2;forj=l:(b-a)/lixlO=l1;x20=l1;tol=le-10;xl=xl0(l);x2=x20(l);yl

48、=xl0(2);y2=x20(2);Kl=fal(klJc2)(kl-feval(fl(j)+h*(3-sqrt(3)/6),YlG)-Hi*(l/4)*kl+h*(3-2*sqit(3)/12)*k2),xl,yl);K2=feval(klJc2)(k2-feval(fl(j)-Hi*(3+sqrt(3)/6),Yl(j)-Hi*(3+2*sqrt(3)/12)*kl+h*(l/4)*k2),xl,yl);%迭代函數(shù)Gl=feval(klJc2)(klfeval(f2X(j)+h*(3sqrt(3)/6),Y2(j)+h*(l/4)*kl+h*(3-2*sqrt(3)/12)*k2),x2,

49、y2);G2=feval(kl,k2)(k2-feval(I2<(j)-Hi*(3-bs<pl(3)/6),Y2g)+h*(3-l-2+scpt(3)/12)*kH-h*(l/4)*k2)K2,y2);%迭代函數(shù)F1=K1K2;F2=G1G2;Jl=double(subs(jacobian(kl-feval(fl,X0)+h*(3-sqrt(3)/6),Yl(j)+h*(l/4)*kl+h*(3-2*sqrt(3)/12)*k2);k2-feval(fl(j)+h*(3+sqrt(3)/6),Yl(j)44i*(3+2*sqit(3)/12)*kl-Hi*(l/4)*k2)J,klk2),kl>k2,xl,yl);J2=double(siibs(jacobian(kl-feval(f2,X(j)4-h*(3-sqit(3)/6),Y2(j)-Hi,*!(l/4)*kl-Hi*(3-2*sqrt(3)/12)*k2);k2-feval(f2(j)+h*(3+scp-t(3)/6),Y2g)-Hi*(34-2*scp