(*---------- Package that constructs embedded Runge-Kutta pairs of orders 9(8) ----------*) (*------------ Use rational numbers in the input, otherwise modify properly -------------*) BeginPackage["T98`"]; Clear["T98`*"]; T98::usage="T98[d2,d5,d9,d10,d11,d13,d14,a116,w16,ww15,ww16] returns the coefficients matrices a,w,ww,d of RK embedded pair of orders 9(8)" Begin["`Private`"]; Clear["T98`Private`*"]; T98[d2_,d5_,d9_,d10_,d11_,d13_,d14_,a116_,ww15_,ww16_,w16_]:= Module[{d,d3,d4,d6,d7,d8,d12,w,ww,e,a,w1,w6,w7,w8,w9,w10,w11,w12,w13,w14,w15,ww1,ww6,ww7,ww8, ww9,ww10,ww11,ww12,ww13,ww14,ae,ad,ad2,ad3,ad4,dd,ii,wadi,wwadi,a32,a43,a53,a54,a64, a65,a74,a75,a76,a84,a85,a86,a87,a94,a95,a96,a97,a98,a104,a105,a106,a107,a108,a109, a114,a115,a117,a118,a119,a1110,a124,a125,a126,a127,a128,a129,a1210,a1211,a134,a135, a136,a137,a138,a139,a1310,a1311,a1312,a146,a147,a148,a149,a1410,a1411,a1412,a1413, a156,a157,a158,a159,a1510,a1511,a1512,a1513,a1514,a166,a167,a168,a169,a1610,a1611, a1612,a1613,a1614,a21,a31,a41,a51,a61,a71,a81,a91,a101,a111,a121,a131,a141,a151, a161,p,pp,x}, a = {{0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}, {a21,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}, {a31,a32,0,0,0,0,0,0,0,0,0,0,0,0,0,0}, {a41,0,a43,0,0,0,0,0,0,0,0,0,0,0,0,0}, {a51,0,a53,a54,0,0,0,0,0,0,0,0,0,0,0,0}, {a61,0,0,a64,a65,0,0,0,0,0,0,0,0,0,0,0}, {a71,0,0,a74,a75,a76,0,0,0,0,0,0,0,0,0,0}, {a81,0,0,0,0,a86,a87,0,0,0,0,0,0,0,0,0}, {a91,0,0,0,0,a96,a97,a98,0,0,0,0,0,0,0,0}, {a101,0,0,0,0,a106,a107,a108,a109,0,0,0,0,0,0,0}, {a111,0,0,0,0,a116,a117,a118,a119,a1110,0,0,0,0,0,0}, {a121,0,0,0,0,a126,a127,a128,a129,a1210,a1211,0,0,0,0,0}, {a131,0,0,0,0,a136,a137,a138,a139,a1310,a1311,a1312,0,0,0,0}, {a141,0,0,0,0,a146,a147,a148,a149,a1410,a1411,a1412,a1413,0,0,0}, {a151,0,0,0,0,a156,a157,a158,a159,a1510,a1511,a1512,a1513,a1514,0,0}, {a161,0,0,0,0,a166,a167,a168,a169,a1610,a1611,a1612,a1613,a1614,0,0}}; d = {0,d2,d3,d4,d5,d6,d7,d8,d9,d10,d11,d12,d13,d14,1,1}; w = {w1,0,0,0,0,0,0,w8,w9,w10,w11,w12,w13,w14,w15,w16}; ww = {ww1,0,0,0,0,0,0,ww8,ww9,ww10,ww11,ww12,ww13,ww14,ww15,ww16}; e = {1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1}; d7 = d9*(2*Sqrt[6]/15+4/5); d8 = 4*d9/3; d6 = d8*(4*d7-3*d8)/(2*(3*d7-2*d8)); d4 = d6*(4*d5-3*d6)/(2*(3*d5-2*d6)); d3 = 2*d4/3; a32 = d3^2/(2*d2); a43 = d4^2/(2*d3); a53 = 3*d5^2*(3*d4-2*d5)/(4*d4^2); a54 = -d5^2*(d4-d5)/d4^2; ae = a.e-d; ad = a.d-d^2/2; ad2 = a.d^2-d^3/3; ad3 = a.d^3-d^4/4; ad4 = a.d^4-d^5/5; {a64,a65,a74,a75,a76,a86,a87,a96,a97,a98} = Simplify[Solve[Join[ad[[6;;9]],ad2[[6;;9]],ad3[[7;;9;;2]]] == Array[0&,10], {a64,a65,a74,a75,a76,a86,a87,a96,a97,a98}]][[1,All,2]]; p = x*(x-d8)*(x-d9)*(x-d10)*(x-d11); pp = (x-d12)*p; d12 = Simplify[ Solve[ Integrate[p*(1-x)^3/6, {x,0,1}]*Integrate[pp*(1-x), {x,0,1}] -Integrate[p*(1-x)^2/2,{x,0,1}]*Integrate[pp*(1-x)^2/2,{x,0,1}] == 0,d12][[1,1,2]]]; {w8,w9,w10,w11,w12,w13,w14,w15} = Solve[{w.d==1/2, w.d^2==1/3, w.d^3==1/4, w.d^4==1/5, w.d^5==1/6,w.d^6==1/7,w.d^7==1/8,w.d^8==1/9},{w8,w9,w10,w11,w12,w13,w14,w15}][[1,All,2]]; {ww8,ww9,ww10,ww11,ww12,ww13,ww14} = Solve[{ww.d==1/2,ww.d^2==1/3,ww.d^3==1/4,ww.d^4==1/5, ww.d^5==1/6,ww.d^6==1/7,ww.d^7==1/8},{ww8,ww9,ww10,ww11,ww12,ww13,ww14}][[1,All,2]]; dd = DiagonalMatrix[d];ii=IdentityMatrix[16]; wadi = w.(a+dd-ii); wwadi = ww.(a+dd-ii); {a1514,a1614} = Solve[{wadi[[14]] == 0, wwadi[[14]] == 0}, {a1514,a1614}][[1,All,2]]; {a106,a107,a108,a109,a117,a118,a119,a1110} = Simplify[ Solve[ Join[ad[[10;;11]], ad2[[10;;11]],ad3[[10;;11]],ad4[[10;;11]]]==Array[0&,8], {a106,a107,a108,a109,a117,a118,a119,a1110}][[1,All,2]]]; {a126,a136,a146} = Simplify[ Solve[{(w.(dd-ii).a)[[6]] == 0,(ww.(dd-ii).a)[[6]] == 0, (w.(dd-ii).(dd-d14*ii).a)[[6]]==0},{a126,a136,a146}][[1,All,2]]]; {a156,a166} = Solve[{wadi[[6]] == 0, wwadi[[6]] == 0},{a156,a166}][[1,All,2]]; {a127,a137,a147,a157,a167} = Solve[ {(w.(dd-ii).a)[[7]] == 0, (ww.(dd-ii).a)[[7]] == 0, (w.(dd-ii).(dd-d14*ii).a)[[7]] == 0, wadi[[7]] == 0, wwadi[[7]] == 0}, {a127,a137,a147,a157,a167}][[1,All,2]]; {a128,a129,a1210,a1211} = Simplify[ Solve[{ad[[12]] == 0, ad2[[12]] == 0, ad3[[12]] == 0, ad4[[12]]==0},{a128,a129,a1210,a1211}][[1,All,2]]]; a1413 = Simplify[ Solve[ w.(dd-ii).a.(dd-d10*ii).(dd-d9*ii).(dd-d8*ii).(dd-d12*ii).(dd-d11*ii).d== Integrate[(x-1)*Integrate[(x-d10)*(x-d9)*(x-d8)*(x-d12)*(x-d11)*x, {x, 0, x}], {x, 0, 1}], a1413]][[1,1,2]]; {a1513,a1613} = Solve[ {wadi[[13]] == 0, wwadi[[13]] == 0}, {a1513,a1613}][[1,All,2]]; {a1312,a1412} = Solve[ {ww.(dd-ii).a.(dd-d9*ii).(dd-d8*ii).(dd-d11*ii).(dd-d10*ii).d == Integrate[(x-1)*Integrate[(x-d9)*(x-d8)*(x-d11)*(x-d10)*x, {x, 0, x}], {x, 0, 1}], w.(dd-ii).(dd-d12*ii).a.(dd-d10*ii).(dd-d9*ii).(dd-d8*ii).(dd-d11*ii).d == Integrate[(x-1)*(x-d12)*Integrate[(x-d10)*(x-d9)*(x-d8)*(x-d11)*x, {x, 0, x}],{x, 0, 1}]}, {a1312,a1412}][[1,All,2]]; {a138,a139,a1310,a1311,a148,a149,a1410,a1411,a158,a159,a1510,a1511,a1512, a168,a169,a1610,a1611,a1612} = Simplify[ Solve[ Join[ ad[[13;;16]], ad2[[13;;16]], ad3[[13;;16]], ad4[[13;;16]], {w.(d*a.d^5)-1/48,ww.(d*a.d^5)-1/48}] == Array[0&,18], {a138,a139,a1310,a1311,a148,a149,a1410,a1411,a158,a159,a1510,a1511,a1512, a168,a169,a1610,a1611,a1612}][[1,All,2]]]; {a21,a31,a41,a51,a61,a71,a81,a91,a101,a111,a121,a131,a141,a151,a161} = Simplify[ Solve[ ae[[2;;16]] == Array[0&,15], {a21,a31,a41,a51,a61,a71,a81,a91,a101,a111,a121,a131,a141,a151,a161}][[1,All,2]]]; {w1,ww1} = Solve[ {w.e==1, ww.e==1}, {w1, ww1}][[1,All,2]]; Return[{a,w,ww,d}]; End[]; EndPackage[]]; (*----------------------------------------------------------------------------------------*) (*----------- coefficients of T9(8) pair for quadtruple precision computations -----------*) (*--------------------- given for use with build-in function NDSolve ---------------------*) T98cvec= {3/184,1005991230711768/26022242251007185,1003144597305563/17299071752613603, 3230/3269,3634874590315107/41460143603477948,840156096102871/4026814668841799,703/2847, 703/3796,2135/9119,1531/2485,15434796193306477/18528760750836549,1999/2000,6594/6595,1,1}; T98amat={{3/184}, {-433678958387722/60461976373487901,800463031060200/17465287845577487}, {377246711516913/26022242251007185,0,1003144597305563/23065429003484804}, {5232938947227754278/42414207887323079,0,-18576057673326955826/47337109241323507, 9722802843766070268/36006158955671251}, {1274983010403228/59501240010408527,0,0,1653987885575953/24968953451991825, 70447607619/36673234034080855}, {904721379210022/5952931408395025,0,0,-23381308838657835/41704208630353019, 2066826469901/11925118158788070,6895820618639281/11173940517165586}, {703/25623,0,0,0,0,2453249100410123/19386166204217813,3412083729453488/36711207832699739}, {13357/485888,0,0,0,0,2029374314414594/16090898152406377,1287982913423758/28873883427687143, -6327/485888}, {965192994961021/35149230657258373,0,0,0,0,4127541420931126/32665647335015649, 3583837594502223/41664019148508749,-246496779813877/27175788804018514, 28747455070668/8549877477842473}, {63799716039495390/22202823855545599,0,0,0,0,-43523/2744,46511948996933887/19326005113913521, 4190224756003065953/25772283742746511,4202436480827105093/56735195078669575, -6231363253366161136/27638382159529459}, {-2708709512927940344/43559765744627999,0,0,0,0,9732832383152222769/28307424515467810, -1467124908252536684/25216366962527245,-123168960431653229617/37767835273866068, -45716261318046911311/29725039922507753,254588505817700334157/55659194024546537, 154479208885061992/61881010299032267}, {7815889611024124839/10273038606706643,0,0,0,0,-132865466831736468380/31621404109692561, 17703281196398756191/24701300463848133,566474708394637174745/14201552570678562, 307285346872690815991/16345778788818343,-2053676792180155828764/36714291312129259, -78247540885035338/2948693499382457,18260998947401049/15112168351518079}, {21763183378499925071/28449009966588093,0,0,0,0,-131859313072947008828/31210661767011655, 26844746596644780532/37252237041563455,1253469865929367878849/31252862654010364, 463799691295259099617/24536685538656591,-4389603511184674820033/78045747665396128, -1128181715104199933/42276540809151831,25502481921637862/20992853044100563, -9559594005944/54330273344414971}, {14646506801530214193/19100557410937333,0,0,0,0,-100013811078882007482/23616683892616855, 51418842743714417428/71184230362278145,449729579268755668521/11186469186029341, 830367055814329487149/43824976316629463,-1510129653054461421096/26785760772223265, -1401877877271794513/52404745175389990,16736569978314491/13745172576768074, -962836086693/28355588058859954,-9941175715160/45586581308747583}, {5917471745174541109/8344408594020065,0,0,0,0,-107320434306581771105/27401651799311177, 5093729026420265343/7626561482131418,2060040620258320258974/55399901444516827, 707106302890477200174/40350224587566143,-1164323447229089993686/22328563534132263, -831410876325262488/33593954156280041,11423107778999224/9908693206689583, 165211034625068/39884924026051713,-38387832271169/18010889018302554,0}}; T98bvec={922536264951379/23804051189793426,0,0,0,0,0,0,225220497265063051/38410941914692596, 36148972660529108/25187799949529035,-198426619086731179/28886469242593278, 6156310435334661/21948640431420694,4487700950738411/30126758819889407, 29127285852344652967/26541231024453067,-98614101020366183246/27482176375033227, 126863537044611133977/50939017338016615,1421/3078}; T98evec= {32020052252709/30901738203108092,0,0,0,0,0,0,-5097456320246635/16138439327639183, -14836286301942601/220086346297809410,8712584401648310/23055283999904231, 201598787147503/16251552598785487,-775224267694528/25963200824081993, 37505781412258941025/34357113542750541,-70601687975788697909/19705466275422783, 451638954762572460757303/181291962706001132785,0}/10; T98Coefficients[9, p_] := N[{T98amat, T98bvec, T98cvec, T98evec}, p]; (*----------------------------------------------------------------------------------------*) (*------------------ run inhomogeneous for tol=10^-16,10^-17,...10^-24 -------------------*) Needs["DifferentialEquations`NDSolveProblems`"]; Needs["DifferentialEquations`NDSolveUtilities`"]; system = NDSolveProblem[{{y1'[t] == y2[t], y2'[t] == -100*y1[t] + 99*Sin[t]}, {y1[0] == 1, y2[0] == 11}, {y1[t], y2[t]}, {t, 0, 20*Pi}, {}, {}, {}}]; refsol = {1, 11}; TabulateResults[labels_List, names_List, data_List] := DisplayForm[ FrameBox[ GridBox[Apply[{labels, ##} &, MapThread[Prepend, {data, names}]], RowLines -> True, ColumnLines -> True]]] /; SameQ[Length[names], Length[data]]; T98 = {"ExplicitRungeKutta", "Coefficients" -> T98Coefficients, "DifferenceOrder" -> 9, "StiffnessTest" -> False}; V98 = {"ExplicitRungeKutta", "DifferenceOrder" -> 9, "StiffnessTest" -> False}; names = {"T9(8)", "V98"}; methods = {T98, V98}; labels = {"Method", "Steps", "Cost", "Error"}; resul1 = Table[{0, 0}, {uu, 16, 24}]; resul2 = resul1; Table[data = CompareMethods[system, refsol, methods, WorkingPrecision -> 33, AccuracyGoal -> uu, PrecisionGoal -> uu]; resul1[[uu - 15, ;;]] = data[[1, 2 ;; 3]]; resul2[[uu - 15, ;;]] = data[[2, 2 ;; 3]]; TabulateResults[labels, names, data], {uu, 16, 24}] (*----------------------------------------------------------------------------------------*)