(*---------- Package that constructs embedded Runge-Kutta pairs of orders 8(7) ----------*) BeginPackage["T87`"]; Clear["T87`*"]; T87::usage="T87[c2,c5,c6,c7,c8,c10,c11,a87,w13,ww12,ww13] returns the coefficients matrices a,w,ww,c of RK embedded pair of orders 8(7)" Begin["`Private`"]; Clear["T87`Private`*"]; T87[c2_,c5_,c6_,c7_,c8_,c10_,c11_,a87_,w13_,ww12_,ww13_]:= Module[ {c, c3, c4, c9, w, ww, e, a, w1, w6, w7, w8, w9, w10, w11, w12, ww1, ww6, ww7, ww8, ww9, ww10, ww11, ae, ac, ac2, ac3, cc, ii, baci, bbaci, a32, a43, a53, a54, a64, a65, a74, a75, a76, a84, a85, a86, a94, a95, a96, a97, a98, a104, a105, a106, a107, a108, a109, a114, a115, a116, a117, a118, a119, a1110, a124, a125, a126, a127, a128, a129, a1210, a1211, a134, a135, a136, a137, a138, a139, a1310, a1311, a21, a31, a41, a51, a61, a71, a81, a91, a101, a111, a121, a131}, c = {0, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11, 1, 1}; w = {w1, 0, 0, 0, 0, w6, w7, w8, w9, w10, w11, w12, w13}; ww = {ww1, 0, 0, 0, 0, ww6, ww7, ww8, ww9, ww10, ww11, ww12, ww13}; e = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}; a = {{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}, {a31, a32, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {a41, 0, a43, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {a51, 0, a53, a54, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {a61, 0, 0, a64, a65, 0, 0, 0, 0, 0, 0, 0, 0}, {a71, 0, 0, a74, a75, a76, 0, 0, 0, 0, 0, 0, 0}, {a81, 0, 0, a84, a85, a86, a87, 0, 0, 0, 0, 0, 0}, {a91, 0, 0, a94, a95, a96, a97, a98, 0, 0, 0, 0, 0}, {a101, 0, 0, a104, a105, a106, a107, a108, a109, 0, 0, 0, 0}, {a111, 0, 0, a114, a115, a116, a117, a118, a119, a1110, 0, 0, 0}, {a121, 0, 0, a124, a125, a126, a127, a128, a129, a1210, a1211, 0, 0}, {a131, 0, 0, a134, a135, a136, a137, a138, a139, a1310, a1311, 0, 0}}; c9 = (14*c6^2*(7*c7^2*c8 + c7*(7*c8^2 - 12*c8 + 1) + c8) + c6*(14*c7^2*(7*c8^2 - 12*c8 + 1) - 7*c7*(24*c8^2 - 33*c8 + 4) + 14*c8^2 - 28*c8 + 3) + 14*c7^2*c8 + c7*(14*c8^2 - 28*c8 + 3) + 3*c8)/ (2*(7*c6^2*(7*c7^2*(15*c8^2 - 10*c8 + 2) - 2*c7*(35*c8^2 - 26*c8 + 6) + 14*c8^2 - 12*c8 + 3) - 7*c6*(2*c7^2*(35*c8^2 - 26*c8 + 6) - c7*(52*c8^2 - 42*c8 + 11) + 12*c8^2 - 11*c8 + 3) + 7*c7^2*(14*c8^2 - 12*c8 + 3) - 7*c7*(12*c8^2 - 11*c8 + 3) + 21*c8^2 - 21*c8 + 6)); c3 = (3*c2)/2; c4 = (9*c2)/4; a32 = (9*c2)/8; a43 = (27*c2)/16; {w1, w6, w7, w8, w9, w10, w11, w12} = Solve[{w.e == 1, w.c == 1/2, w.c^2 == 1/3, w.c^3 == 1/4, w.c^4 == 1/5, w.c^5 == 1/6, w.c^6 == 1/7, w.c^7 == 1/8}, {w1, w6, w7, w8, w9, w10, w11, w12}][[1, All, 2]]; {ww1, ww6, ww7, ww8, ww9, ww10, ww11} = Solve[{ww.e == 1, ww.c == 1/2, ww.c^2 == 1/3, ww.c^3 == 1/4, ww.c^4 == 1/5, ww.c^5 == 1/6, ww.c^6 == 1/7}, {ww1, ww6, ww7, ww8, ww9, ww10, ww11}][[1, All, 2]]; ae = a.e - c; ac = a.c - c^2/2; ac2 = a.c^2 - c^3/3; ac3 = a.c^3 - c^4/4; cc = DiagonalMatrix[c]; ii = IdentityMatrix[13]; baci = w.(a + cc - ii); bbaci = ww.(a + cc - ii); {a53, a54, a64, a65, a74, a75, a76, a84, a85, a86, a94, a95, a96, a97, a98, a104, a105, a106, a107, a108, a109, a114, a115, a116, a117, a118, a119, a1110, a124, a125, a126, a127, a128, a129, a1210, a1211, a134, a135, a136, a137, a138, a139, a1310, a1311} = Solve[ Join[(w.(cc - ii).a)[[4 ;; 5]], (w.(cc - ii).(cc - ii).a)[[4 ;; 5]], ac[[5 ;; 12]], ac2[[5 ;; 12]], ac3[[7 ;; 13]], baci[[4 ;; 10]], bbaci[[4 ;; 8]], {(ww.(cc - ii).a)[[4]], -(1/35) + w.(c a.c^4), -(1/40) + w.(c^2 a.c^4), -(1/48) + w.(c a.c^5), -(1/35) + ww.(c a.c^4)}] == Array[0 &, 44], {a53, a54, a64, a65, a74, a75, a76, a84, a85, a86, a94, a95, a96, a97, a98, a104, a105, a106, a107, a108, a109, a114, a115, a116, a117, a118, a119, a1110, a124, a125, a126, a127, a128, a129, a1210, a1211, a134, a135, a136, a137, a138, a139, a1310, a1311}][[1, All, 2]]; {a21, a31, a41, a51, a61, a71, a81, a91, a101, a111, a121, a131} = Simplify[Solve[ae[[2 ;; 13]] == Array[0 &, 12], {a21, a31, a41, a51, a61, a71, a81, a91, a101, a111, a121, a131}][[1, All, 2]]]; Return[{a,w,ww,c}]]; (*----------------------------------------------------------------------------------------*) (*----------- coefficients of T8(7) pair for quadtruple precision computations -----------*) (*--------------------- given for use with build-in function NDSolve ---------------------*) T87bvec= {959469921003535/20735873900418433,0,0,0,0,83661087663817387/226096222469839182, 228743606234324881/883020026679163794,3544120671195926375/8063503515187523, 164403934540876/64548125027903185,1872154679941434671/50440600905843744, -3908844507545666995/8324248434152054,-402658040159189839/58491143516062232,16491/120125}; T87evec= {10839870895445/185203278486104297,0,0,0,0,70876466420204/75470597442438275, -16496614726651/75468957694148719,54859577937538923405/14355606386435513, 17895137500075704819/2362101251104030,-969063659258770673/19161965088242370, -39193463899416576680/3455993185375433,-973424986199410385/155592419305288032,16491/120125}/10; T87cvec= {3102/110773,41448895555141/353624691619188,41448895555141/235749794412792, 49442/119883,51187/105369,61011/376738,77114/79499,147909751614626799/152923788158104127, 74279/78046,72043/74409,1,1}; T87amat={ {3102/110773}, {-17033458900934993/132978864382888258,17659313382611255/71989792689293837}, {41448895555141/942999177651168,0,41448895555141/314333059217056}, {33544131897542527/99303639017753176,0,-123806032279621065/100880451772826828, 80881552191452041/62126727673226683}, {3901178494518027/70202052982346435,0,0,12244602153330846/48744104078022083, 11363782051482252/63479278340035273}, {7281184019796491/108906123149933189,0,0,8912953764743186/75237479424494327, -1193193435755019/24043824215671157,3001381510813201/114340525306552991}, {-297808918551351805/103302384399153762,0,0,-2387409947307450796/38235137422988677, -320655295147743895/172685972706995386,266830735262229145/73369592821183637,8174527/126711}, {-312230898179118543/111335375555652709,0,0,-5921685522031592717/97516557935639304, -122516042059134140/66440638491697461,143089054978597281/39930960285352934, 1966780853930863533/31340008936176199,27204097600957/30119714219091834}, {-497327926559154029/208366132906665209,0,0,-2070519061247416919/40105304012179956, -139926368413626755/79789745208684688,436822604663916242/133157501626893287, 4951999978536596383/92678477827402881,-1662171172972759/32043786293542537, 320510318790859/5467452906511140}, {-267997292446794835/94625648159795289,0,0,-1326916430444389167/21635054137957163, -50510473210813287/27322222661367848,680595213260915461/188925642391189177, 1090597603926315985/17207867085312708,-818226826952911/56758278493554744, 794276136679319/44163223221855014,-495594365453263/165024671142376612}, {-286074472550848766/70568381571246193,0,0,-2666282586603439301/29766446888618900, -394981932622811234/181671027945865139,354437914440687571/72293255173230666, 1737172167669457231/18855481952627537,-1908527156826626453/17978177470082379, 14359180611877865064/20075894067162869,-1863006586402493967/31715262582627044, -5146117877451253921/9346764321565133}, {-2286460617615599450/148215689608432541,0,0,-21511651826330234931/52669819756106150, -949790098629780736/69310896259636617,2488552272190713800/64326656295428697, 14577683994864478388/35463253730030943,-34626716477448076238/6579786536866391, 267076469802229885930/7436961774107587,-15666088518007151408/5323429123670105, -39614246945332388915/1429199330541022,0}}; T87Coefficients[8,p_]:=N[{T87amat,T87bvec,T87cvec,T87evec},p]; (*----------------------------------------------------------------------------------------*)