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

1、HUBEIFOIJYTECMNICUNIVERSITY畢業(yè)設計題目：剛性系統(tǒng)的隱式RK方法學院：數(shù)理學院專業(yè)名稱：信息與計算科學學號：201241210127學生姓名：丁楠指導教師：汪玉霞摘要本 文 主 要 介 紹 單 步 隱 式 ？ ? 方 彷 簡 要 的 介 紹 了 ？ ? 舞 蹴 ?式?、？?？?？式？?？?珈?？和 ？ ？ ？ ？ ？ ？ ？ ? ？ ? 并 利 用 這 些 基 本 的 隱式？?？對剛性方程組進行數(shù)值求解，并將隱式？?？顯式經(jīng)典？?？解的結(jié)果進行對比，說明兩種數(shù)值解法的優(yōu)缺點。關(guān)鍵詞：剛性系統(tǒng)隱式？?單步方法？?？AbstractThispapermainlyintro

2、ducestheImplicitRunge-KuttaMethodsandasimpledescriptionofGaussimplicitRunge-Kuttamethod,RadauimplicitRunge-KuttamethodandLobattoimplicitRunge-Kuttamethod.ThesebasicImplicitRunge-Kuttamethodsareusedtosolvethestiffequations.TheseimplicitRungeKuttamethodsiarecomparedwiththeclassicalexplicitRunge-Kuttam

3、ethod.Thispaperexplaintheadvantagesanddisadvantagesofthetwokindofnumericalmethods.Keywords:StiffsystemImplicitRunge-KuttamethodOnestepmethodNewtoniterativemethod目錄1 .緒論 1.1.1 剛性方程 1.1.2 隱式 RK 方法的研究意義 21.3 RK 方法的研究現(xiàn)狀 3.2 .單步 RK 方法的收斂性和穩(wěn)定性 52.1 單步 RK 方法的收斂性.5.2.2 單步 RK 方法的穩(wěn)定性.6.3 .三類隱式 RK 方法 8.3.1 引言

4、8.3.2 Gauss 型隱式 RK 方法 9.3.3 Radau 型隱式 RK 方法 1.03.4 Lobatto 型隱式 RK 方法 1.14 隱式 RK 方法的實現(xiàn) 134.1 非線性系統(tǒng)的改進 134.2 簡化的 Newton 迭代法 1.35 數(shù)值實驗與結(jié)果分析 15參考文獻.18附錄 211.緒論1.1 剛性方程對于一般的線性常系數(shù)系統(tǒng)?=?（?）?必？X?勺矩陣，特征值為？決?？1,2,?,?）0定義 123若一個系統(tǒng)滿足(1)?0,?1,2,?,?(2)?儂?靄?？?？1其中？叫剛性比，則這個系統(tǒng)稱為剛性系統(tǒng)。定義 227若線性系統(tǒng)_?=?0,?或非線性系統(tǒng)?=?(?)?0,?

5、的矩陣？城？?？?？?？?特征值？到足?11V?則其是剛性的。定理 1（解的存在性與唯一性）（1）對于所有（?？e?函數(shù)？?是連續(xù)的；（2）對于任何（?,（?）e?存在常數(shù)？是函數(shù)滿足I?-?（?）IV?那則初值問題?)(?=?有唯一解。其 中 ？ ?= （ ?,?,?,?）?,?= （ ?|?0?-oo?kOO；-OO?0.0278）,計算結(jié)果就完全不可信了。對于剛性方程組，顯式方法已經(jīng)遠遠不能勝任了，一般采用絕對穩(wěn)定性更好的方法（如隱式？?方法） 進行求解， 本文采用單步隱式？?？而對于隱式方法中的級值得求解，本文采用？?？。1.3 RK 方法的研究現(xiàn)狀研究基于標準模型方程的？?？常見形式

6、為：?+1=?+?E?=1?=?-1?=?（?+?,?+?E?（?=2,3,？?）?=1（顯式）?+1=?+?E?=1?=?（?+?,?+?E?（?=1,2,3,？?）?=1（隱式）因為？?如海？比較成熟的常微分方程數(shù)值解法。所以如今主要是對于經(jīng)典的？?雷？行完善和擴充，在一定的條件下，提高級數(shù)以提高精度。或者是將？?燃四？某些領域結(jié)合使用。在 1994 年，費景高1給出了一種顯式？?？法，從而實現(xiàn)？?？多處理機上的應用。1997年，Enenkel 和 JacksorF 實現(xiàn)了？?？?？對角隱式0.010.0010.00019.9004983e-019.9004983e-012.5300364

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

8、RK 方法的收斂性對于常微分初值問題?(?7c?公？)?=?的單步顯式公式?+1=?+?(?)(?=0,1,?,?1)1?=?局部截斷誤差可以表示為?+1=?+1)-?+?)定理 216:設?(?)(1)的解，?嫉?=(2)的解。如果:(1)存在常數(shù)？使得(?=0,1,?,?1)|?+1|(2)存在？?0,使得?歹+1(=?0,1,2,?,?1)?(西)|?(?,法切?10v?v?.其中？?=(?|?-?0?&?+?CfCCCC一,.?記?0?=力?-?1-1,則當？&?,7時，有0?(?)?(?證:由(3)得?+1)=?+?(?)+?+1伴？0,1,?,?-1)將(4)與(2)相減?+1)-

9、?=?-?汾？?務？)-?)+?+1(?=0,1,?,?1)由？?)-?=0 知道，在？0 時，|?-?/0?(?城立。?.現(xiàn)在假設 00?實？1 時也是成立的。在？&、七時有：?|?7?在？)-?7?)|?-?趙其中？戲?和？彩間的數(shù)。于是定理結(jié)合條件與(4)式，可以得到|?+i)-?%j|?-?+?|?,?)-?)|+|?+i|?-?+?-?+?與+1=(1+?|?-?+?+10?c?1從而(1+?-1|?(?)|(1+?)-?|+(1+?-1?因為？?)-?=0 及(1+?0?？?(?-?),得到?|?-?w?-?-1?=?所以當？時|?-?務&?也是成立的。(證畢)對于單步顯式格式而言

10、，當？?(?和?:翳翟？?內(nèi)連續(xù)時，那么定理 1 的條件(2)總是滿足的。所以滿足上述條件的單步顯式公式，只要也滿足相容性條件?),0)=?(?(?)那么它一定是收斂的。2.2 單步 RK 方法的穩(wěn)定性定義 4 對于初值問題(1)，若?4?=是由得到，?金?=思下面擾動問題的解:?+1=?必？?)+?+1(?0,1,2,?,?1)進而_?|?-?&?(?(9=?%?)1?(?)|CK?c?-1?=?+?如果存在正常數(shù)？很？o,使所有的？?e(0,?的，?(0,?0,當？?富嚼?時,有?|?-?”？0V?1就稱(2)是穩(wěn)定的或者零穩(wěn)定的。定理 316在定理 1 的條件下，單步顯式公式(2)是穩(wěn)定

11、的3 .三類隱式RK方法3.1 引言對于初值問題(1),隱式？?勃四?一股形式為?+1=?+?E?(5)?=1?=?+?,?+?E?(?1,2,3,?(6)?=1其中，？?=?+?=0,1,2,?,為數(shù)軸上的離散點列；？為步長，？的解？?(?見勺近似值；？?，？?，？?稱為隱式？?？步長；？?，？，？劾權(quán)系數(shù)。令？?=(?1,2,3，？稱？為系數(shù)矩陣，因此上式可以簡化表示為?1?1?？?？使用？?舞?必，上表可以表示為可以看到，如果？是一個主對角元素均為零的下三角矩陣，那么上表就可以表示一個顯式 ？ ? ？ 如 果 ？ 用 一 個 主 對 角 線 非 零 的 下 三 角 矩 陣 ， 相 應的?

12、是半隱式方法，(5)式等號左右就含有級值？，?,?,?計算？弼寸要解一個包含 r 個未知量？?非線性方程組。如果？圖滿足？?箕？*全為零,則相應的方法就是隱式方法。若系統(tǒng)是 m 維的，對于隱式 RungeKutta 方法實現(xiàn)的每一個遞推步，均需要求解一個 mxr 維的非線性方程組，一般用牛頓迭代法求解。例如0120012-1這是 Kutta 得到的 3 級 3 階顯式？?W?而00011044010211216360121000511243-24121636121636分別是 3 級 4 階的半隱式？?隱式？?要求出具體的？?,？一般有兩種辦法。一種是將(6)在點(？?處進行泰勒展開，并且代入

13、到(5)中，再與？?(?*?)在？班的泰勒展開式進行對比，從而確定參數(shù)?另一種方法就是將微分方程轉(zhuǎn)化為等價的積分方程，用數(shù)值積分得到表達式，再進行對比得到參數(shù)?；诤笠环N方法，可以構(gòu)造得到以下三種不同的隱式Runge-Kutta 方法。3.2Gauss 型隱式 RK 方法17設?，？，？欲/2?1)=0 的根，這里留?0,1上的？：？?項式，0?0,當？?一oo時， 得到的數(shù)值解趨于零， 則稱這種方法是？?急定的。定理 5910?皴？?？?方/？2?O 勺， 其穩(wěn)定函數(shù)是(?浮?以且是？?急定的。以下列出？?=3,?=6的？?？?？?？?雄四？一：524555369-1536-305V1525

14、V1536+24936-245V152等536+-30-9+行36545189183.3 Radau 型隱式 RK 方法17這種方法的參數(shù)？弟要下面幾個步驟求得:(1)求多項式解?？-?-1(?)勺零點？1?,?,，?并指定？1?0,?=1;(2)計算系數(shù)？為囹?？-?-1(?)(?制?(?)-?-1(?劉?=?,?,?(3)計算系數(shù)矩陣?1V152-記12_1巫2+70-1?=/0?=?例如？笑 3,?=5 的？?？?第期為:16-v616+v6136369631?和？徼？?2?1 階的，穩(wěn)定函數(shù)是（?-1,?冰對角？?這兩種都是？?急定的。3.4 Lobatto 型隱式 RK 方法17這種

15、方法的參數(shù)？弟要下面幾個步驟求得:求多項式？?（?-?-2（?的零點？?，？?，?并指定？1?=0,?=1；計算系數(shù)？為?=?1(?3-1(?-?IKf列出 4 級 6 階？?隱式？?121212124-v6104+610188-7V6360296+169V6180016-淺36296-169V6-2+3V6180088+7V636016+v636225-2-3V622519(1)(2)?=倒?-?-2(?)?(?)?(?)-?-2(?)=?,?,?(3)計算系數(shù)矩陣?-70_5-v510_5+v5101011+v5120_11-v5120112025-v5120_25+13V51205120

16、25-13V5120_25+v51205120-1+v5120_-1-v5120112定理 731?田？型隱式？?？?？?？2?2 階的，？?他？?時?理?少游？方法的穩(wěn)定函數(shù)是（?？1,?1）對角?W?M?（?2,?次對角？?所以這些方法都是？?急定的。4隱式RK方法的實現(xiàn)4.1 非線性系統(tǒng)的改進對于？皴隱式隱式？?？?=?+?(?+?,?1,?,?(7)?=1?+1=?+?E?(?+?,冽(8)?=1令?=?-?于是變?yōu)?=?E?(?+?,?+初(10)?=1當?shù)玫降?10)解?,？， ？寸， 由顯式(8)可以很快得到解？?+1。 若隱式？?方法的系數(shù)矩陣式非奇異的，那么(10)可以寫為?

17、+?,?+?)?=?(11)?(?+?,?+?所以(8)可以看成與?+1=?+E?=1等價的，其中(?,?,?=(?,?(12)比如？?=3,?=5 的？?？?,?=(0,0,1)。4.2 簡化的 Newton 迭代法對于一般的非線性微分方程，系統(tǒng)（10）必須用迭代的方法求解。本文采用?。？?？應用于系統(tǒng)（10）,每次迭代都需要一個系數(shù)矩陣?1而尹+?,?+?)?為務?+?,?+?-?1還必?+?1?,?2+?2?r巧劣?%+?,?+?分I,I的線性方程組。定義 5 若?階矩陣？?=(?人？確夕!陣?=(?且有?B?A=?則稱運算？是？為?為了簡化（10）的計算，我們用近似值?2K 我?）代替

18、？?硒?+?,?+?少的值，于是方程（10）的簡化？?曬?？變?yōu)??=-?(?+?(?)?(?)/zI?夕?+1)=?)+?(?)?=(?,?,?)?建解白向第 k 次近似,?=(?),?,&?)?是增量，？?(涔)是？?a?)=(?(?+?,?+?，?),？，？?(?+?,?+?)的縮寫。每一次的迭代需要？?由函數(shù)？?勺求值和一個？?維線性方程組的求解。 因為矩陣 （？?寸于所有的迭代是相同的，所以其？?？解在一個積分步內(nèi)只要進行一次，進而減少了計算時間。(13)這里？夕?+1）5數(shù)值實驗與結(jié)果分析10-7x?+?(0)=110-3X?+?(0)=1?)03)(1)+(?)10-?2?)10

19、-3很明顯該方程組的剛性比？?=1071012?1,因此是個剛性方程組。1+313+23?=?(?+6?,邙4?+12?)進行求解。得到結(jié)果如下表一：真值結(jié)果問題：寫成矩陣形式就是:?=?=,2?(QQ) 2(10-7(0取步長為 0.1,在區(qū)間0,1內(nèi)分別用四級四階經(jīng)典顯式？?1?+1=?海 6(?+2?+2?+?)?=?11?=?7?計 2?,?+2?)?=?7?&+2?,?+2?)33?=?+?,?+?)和二級四階隱式？?？?？?1?+1=?+2(?+?)3-v313-23?=?(?+6?,?計 4?+12?)0.0000000e+001.0000000e-012.0000000e-01

20、3.0000000e-011.0000000e+001.0049958e+001.0199334e+001.0446635e+001.0000000e+001.0999384e+001.1988893e+001.2958649e+00?4.0000000e-015.0000000e-016.0000000e-017.0000000e-018.0000000e-019.0000000e-011.0000000e+001.0789390e+001.1224175e+001.1746644e+001.2351579e+001.3032934e+001.3783901e+001.4596978e+00

21、1.3898974e+001.4800481e+001.5654174e+001.6451531e+001.7184598e+001.7846058e+001.8429313e+00表二：隱式結(jié)果?%?0.0000000e+001.0000000e-012.0000000e-013.0000000e-014.0000000e-015.0000000e-016.0000000e-017.0000000e-018.0000000e-019.0000000e-011.0000000e+001.0000000e+001.0049958e+001.0199334e+001.0446635e+001.07

22、89390e+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.6451531e+001.7184598e+001.7846058e+001.8429313e+00表二：顯式結(jié)果?0.0000000e+001.0000000e-012.0000000e-011.0000000e+001.0049958e

23、+001.0199334e+001.0000000e+001.0999336e+001.1988802e+003.0000000e-014.0000000e-015.0000000e-016.0000000e-017.0000000e-018.0000000e-019.0000000e-011.0000000e+001.0446635e+001.0789390e+001.1224175e+001.1746645e+001.2351579e+001.3032934e+001.3783901e+001.4596978e+001.2958521e+001.3898814e+001.4800297e+

24、001.5653972e+001.6451319e+001.7184382e+001.7845846e+001.8429111e+00很明顯的看到，在保留八位小數(shù)的情況下，二級四階隱式？?胡渴真值無異，能夠精確到小數(shù)點后七位以上。而經(jīng)典？?？確到小數(shù)點后四位。因此對于剛性方程組，相同步長下，隱式？?？精度比顯式？?得多。并且由于步長小，在實際的求解區(qū)間里面，涉及的遞推次數(shù)少，從而函數(shù)的計算量就小。參考文獻1費景高.并行顯式 Runge-Kutta 公式的實現(xiàn)J.計算機工程與設計,1994(5):34-40.2EnenkelRF,JacksonKR.DIMSEMs-diagonallyimpli

25、citsingle-eigenvaluemethodsforthenumericalsolutionofstiffODEsonparallelcomputersJ.AdvancesinComputationalMathematics,1997,7(1-2):97-133.3張誠堅，余紅兵.非線性延遲系統(tǒng)的一類并行預校算法J.華中理工大學學報，2000,28(11):104-106.4李愛雄，劉偉豐，張誠堅，等.求解 DDEs 的多導龍格庫塔方法的漸近穩(wěn)定性Asymptoticstabilityofmulti-derivativeRunge-Kuttamethodfordelaydifferen

26、tialequationsJ.華中科技大學學報(自然科學版)，2002,6:108-110.5廖文遠，李慶揚.一類求解剛性常微分方程的半隱式多步 RK 方法J.清華大學學報:自然科學版,1999,39(6):1-4.6龐立君，朱永忠.類隨機微分方程 Runge.Kutta 方法的指數(shù)穩(wěn)定性J.河海大學學報(自然科學版),2008,36(3).7 CurtissCF,HirschfelderJO.IntegrationofstiffequationsJ.Proc.Nat.Acad.Sci,1952,38(235):1.8 GearCW.常微分方程初值問題的數(shù)值解法J.1978.9 Butcher

27、JC.Implicitrunge-kuttaprocessesJ.MathematicsofComputation,1964,18(85):50-64.10 EhleBL.HighorderA-stablemethodsforthenumericalsolutionofsystemsofDEsJ.BITNumericalMathematics,1968,8(4):276-278.11 BurrageK,ButcherJC,ChipmanFH.Animplementationofsingly-implicitRunge-KuttamethodsJ.BITNumericalMathematics,

28、1980,20(3):326340.12 ShampineLF.ImplementationofimplicitformulasforthesolutionofODEsJ.SIAMJournalonScientificandStatisticalComputing,1980,1(1):103-118.13 LindbergB.OnsmoothingandextrapolationforthetrapezoidalruleJ.BITNumericalMathematics,1971,11(1):2952.14 LindbergB.IMPEX2-aprocedureforsolutionofsys

29、temsofstiffdifferentialequationsR.CM-P00069411,1973.15 BaderG,DeuflhardP.Asemi-implicitmid-pointruleforstiffsystemsofordinarydifferentialequationsJ.NumerischeMathematik,1983,41(3):373398.16孫志忠，袁慰平,聞震初.數(shù)值分析M.東南大學出版社，2002.17劉德貴，費景高.剛性大系統(tǒng)數(shù)字仿真方法M.河南科學技術(shù)出版社,1996.18 CockburnB,ShuCW.TheRunge-Kuttadiscontin

30、uousGalerkinmethodforconservationlawsV:multidimensionalsystemsJ.JournalofComputationalPhysics,1998,141(2):199224.19 MonovasilisT,KalogiratouZ,SimosTE.ExponentiallyfittedsymplecticrungeKutta-Nystr?mmethodsJ.Appl.Math.Inf.Sci,2013,7(1):885.20 HochbruckM,PazurT.ImplicitRunge-KuttaMethodsandDiscontinuou

31、sGalerkinDiscretizationsforLinearMaxwellsEquationsJ.SIAMJournalonNumericalAnalysis,2015,53(1):485-507.21 PapadopoulosDF,SimosTE.AmodifiedRunge-Kutta-Nystr?mmethodbyusingphaselagpropertiesforthenumericalsolutionoforbitalproblemsJ.AppliedMathematics&InformationSciences,2013,7(2):433437.22 ZhongX,ShuCW

32、.AsimpleweightedessentiallynonoscillatorylimiterforRungeKuttadiscontinuousGalerkinmethodsJ.JournalofComputationalPhysics,2013,232(1):397-415.23 LambertJD,LambertJD.ComputationalmethodsinordinarydifferentialequationsM.London:Wiley,1973.24 ZhuJ,ZhongX,ShuCW,etal.Runge-KuttadiscontinuousGalerkinmethodu

33、singanewtypeofWENOlimitersonunstructuredmeshesJ.JournalofComputationalPhysics,2013,248:20Q220.25 JayakumarT,MaheshkumarD,KanagarajanK.NumericalsolutionoffuzzydifferentialequationsbyRunge-KuttamethodoforderfiveJ.InternationalJournalofAppliedMathematicalScience,2012,6:29893002.26 DimarcoG,PareschiL.As

34、ymptoticpreservingimplicit-explicitRunge-KuttamethodsfornonlinearkineticequationsJ.SIAMJournalonNumericalAnalysis,2013,51(2):1064-1087.27 MirankerWL,MirankerA.NumericalMethodsforStiffEquations:AndSingularPerturbationProblemsM.SpringerScience&BusinessMedia,2001.28 HadiM,AndersonM,HuseinA.Using4thorde

35、rRunge-KuttamethodforsolvingatwistedSkyrmestringequationC/THE4THINTERNATIONALCONFERENCEONTHEORETICALANDAPPLIEDPHYSICS(ICTAP)2014.AIPPublishing,2016,1719:030005.29 JimenezJC,SotolongoA,SanchezBornotJM.LocallylinearizedrungekuttamethodofdormandandprinceJ.AppliedMathematicsandComputation,2014,247:589-6

36、06.30 JimenezJC,SotolongoA,SanchezBornotJM.LocallylinearizedrungekuttamethodofdormandandprinceJ.AppliedMathematicsandComputation,2014,247:589-606.31 WannerG,HairerE.SolvingordinarydifferentialequationsIIM.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);

37、h=0.0001;%步長T=zeros(1,(b-a)/h)+1);X=zeros(1,(b-a)/h)+1);Y=zeros(1,(b-a)/h)+1);T=a:h:b;X(1)=2;%初值Y(1)=1;forj=1:(b-a)/hk11=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/

38、2*k21,Y(j)+h/2*k22);k32=feval(f2,T(j)+h/2,X(j)+h/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=TXY;endsave0.0001 結(jié)果.txtRascii真值程序：y1=exp(-100)+exp(-0.01);y2=exp(-

39、100);R=1y1y2;save 真值結(jié)果.txtR-ascii第五章程序：symsx;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=XY1Y2;xlswhte(結(jié)果,R1,真值)formatlongf1=(x,y1)(0.0000001*y1+sin(x);f2=(x,y2)(0.001*y2+cos(x);a=0;b=1;h=0.1;X

40、=zeros(1,(b-a)/h+1);Y1=zeros(1,(b-a)/h+1);Y2=zeros(1,(b-a)/h+1);X=a:h:b;Y1(1)=1;Y2(1)=1;symsk1k2;forj=1:(b-a)/hx10=11;x20=11;tol=1e-10;x1=x10(1);x2=x20(1);y1=x10(2);y2=x20(2);K1=feval(k1,k2)(k1-feval(f1,X(j)+h*(3-sqrt(3)/6),Y1(j)+h*(1/4)*k1+h*(3-2*sqrt(3)/12)*k2),x1,y1);K2=feval(k1,k2)(k2-feval(f1,X

41、(j)+h*(3+sqrt(3)/6),Y1(j)+h*(3+2*sqrt(3)/12)*k1+h*(1/4)*k2),x1,y1);%迭代函數(shù)G1=feval(k1,k2)(k1-feval(f2,X(j)+h*(3-sqrt(3)/6),Y2(j)+h*(1/4)*k1+h*(3-2*sqrt(3)/12)*k2),x2,y2);G2=feval(k1,k2)(k2-feval(f2,X(j)+h*(3+sqrt(3)/6),Y2(j)+h*(3+2*sqrt(3)/12)*k1+h*(1/4)*k2),x2,y2);%迭代函數(shù)F1=K1K2;F2=G1G2;J1=double(subs(

42、jacobian(k1-feval(f1,X(j)+h*(3-sqrt(3)/6),Y1(j)+h*(1/4)*k1+h*(3-2*sqrt(3)/12)*k2);k2-feval(f1,X(j)+h*(3+sqrt(3)/6),Y1(j)+h*(3+2*sqrt(3)/12)*k1+h*(1/4)*k2),k1k2),k1,k2,x1,y1);J2=double(subs(jacobian(k1-feval(f2,X(j)+h*(3-sqrt(3)/6),Y2(j)+h*(1/4)*k1+h*(3-2*sqrt(3)/12)*k2);k2-feval(f2,X(j)+h*(3+sqrt(3)/6),Y2(j)+h*(3+2*sqrt(3)/12)*k1+h*(1/4)*k2),k1k2),k1,k2,x2,y2);x11=x10-F1/J1;x21=x20-F2/J2;while(norm(x11-x10)tol)x10=x11;x1=x10;y1=x10(2);K1=feval(k1,k2)(k1-feval(f1,X(j)+h*(3-sqrt(3)/6),Y1(j)+h*(1/4)*k1+h*(3-2*sqrt(3)/12)