# JuliaDiffEq/DiffEqDevTools.jl

Switch branches/tags
Nothing to show
Fetching contributors…
Cannot retrieve contributors at this time
7869 lines (7541 sloc) 427 KB
 #= Classic Verner 8(9) excluded since it's just hard to write from the paper. If anyone wants to add them I'd accept them. =# """ Gauss-Legendre Order 2. """ function constructGL2(T::Type = Float64) c = [1//2] A = [1//2] α = [1] A = map(T,A) α = map(T,α) c = map(T,c) return(ImplicitRKTableau(A,c,α,2)) end """ Gauss-Legendre Order 4. """ function constructGL4(T::Type = Float64) c = [(3-sqrt(3))/6;(3+sqrt(3))/6] A = [ 1/4 (3-2*sqrt(3))/12; (3+2*sqrt(3))/12 1/4 ] α = [1/2;1/2] A = map(T,A) α = map(T,α) c = map(T,c) return(ImplicitRKTableau(A,c,α,4)) end """ Gauss-Legendre Order 6. """ function constructGL6(T::Type = Float64) c = [(5-sqrt(15))/10;1/2;(5+sqrt(15))/10]; A = [ 5/36 (10-3*sqrt(15))/45 (25-6*sqrt(15))/180; (10+3*sqrt(15))/72 2/9 (10-3*sqrt(15))/72; (25+6*sqrt(15))/180 (10+3*sqrt(15))/45 5/36 ]; α = [5/18;4/9;5/18]; A = map(T,A) α = map(T,α) c = map(T,c) return(ImplicitRKTableau(A,c,α,6)) end """ Heun's Order 2 method. """ function constructHeun(T::Type = Float64) A = [0 0 1 0] c = [0;1] α = [1//2;1//2] αEEst = [1;0] A = map(T,A) α = map(T,α) c = map(T,c) αEEst = map(T,αEEst) return(ExplicitRKTableau(A,c,α,2,αEEst=αEEst,adaptiveorder=1)) end """ Ralston's Order 2 method. """ function constructRalston(T::Type = Float64) A = [0 0 2//3 0] c = [0;2//3] α = [1//4;3//4] A = map(T,A) α = map(T,α) c = map(T,c) return(ExplicitRKTableau(A,c,α,2)) end """ Euler's method. """ function constructEuler(T::Type = Float64) A = Matrix{T}(undef,1,1) A[1] = 0 c = [0] α = [1] A = map(T,A) α = map(T,α) c = map(T,c) return(ExplicitRKTableau(A,c,α,1)) end """ Kutta's Order 3 method. """ function constructKutta3(T::Type = Float64) A = [0 0 0 1//2 0 0 -1 2 0] c = [0;1//2;1] α = [1//6;2//3;1//6] A = map(T,A) α = map(T,α) c = map(T,c) return(ExplicitRKTableau(A,c,α,3)) end """ Classic RK4 method. """ function constructRK4(T::Type = Float64) A = [0 0 0 0 1//2 0 0 0 0 1//2 0 0 0 0 1 0] c = [0;1//2;1//2;1] α = [1//6;1//3;1//3;1//6] A = map(T,A) α = map(T,α) c = map(T,c) return(ExplicitRKTableau(A,c,α,4)) end """ Classic RK4 3/8's rule method. """ function constructRK438Rule(T::Type = Float64) A = [0 0 0 0 1//3 0 0 0 -1//3 1 0 0 1 -1 1 0] c = [0;1//3;2//3;1] α = [1//8;3//8;3//8;1//8] A = map(T,A) α = map(T,α) c = map(T,c) return(ExplicitRKTableau(A,c,α,4)) end """ Explicit SSP method of order 2 using 2 stages. """ function constructSSPRK22(T::Type = Float64) A = [0 0 1 0] c = [0; 1] α = [1//2; 1//2] A = map(T,A) α = map(T,α) c = map(T,c) return(ExplicitRKTableau(A,c,α,2)) end """ Explicit SSP method of order 3 using 3 stages. """ function constructSSPRK33(T::Type = Float64) A = [0 0 0 1 0 0 1//4 1//4 0] c = [0; 1; 1//2] α = [1//6; 1//6; 2//3] A = map(T,A) α = map(T,α) c = map(T,c) return(ExplicitRKTableau(A,c,α,3)) end """ Explicit SSP method of order 3 using 4 stages. """ function constructSSPRK43(T::Type = Float64) A = [0 0 0 0 1//2 0 0 0 1//2 1//2 0 0 1//6 1//6 1//6 0] c = [0; 1//2; 1; 1//2] α = [1//6; 1//6; 1//6; 1//2] A = map(T,A) α = map(T,α) c = map(T,c) return(ExplicitRKTableau(A,c,α,3)) end """ Explicit SSP method of order 4 using 10 stages. """ function constructSSPRK104(T::Type = Float64) A = [0 0 0 0 0 0 0 0 0 0 1//6 0 0 0 0 0 0 0 0 0 1//6 1//6 0 0 0 0 0 0 0 0 1//6 1//6 1//6 0 0 0 0 0 0 0 1//6 1//6 1//6 1//6 0 0 0 0 0 0 1//15 1//15 1//15 1//15 1//15 0 0 0 0 0 1//15 1//15 1//15 1//15 1//15 1//6 0 0 0 0 1//15 1//15 1//15 1//15 1//15 1//6 1//6 0 0 0 1//15 1//15 1//15 1//15 1//15 1//6 1//6 1//6 0 0 1//15 1//15 1//15 1//15 1//15 1//6 1//6 1//6 1//6 0] c = [0; 1//6; 1//3; 1//2; 2//3; 1//3; 1//2; 2//3; 5//6; 1] α = [1//10; 1//10; 1//10; 1//10; 1//10; 1//10; 1//10; 1//10; 1//10; 1//10] A = map(T,A) α = map(T,α) c = map(T,c) return(ExplicitRKTableau(A,c,α,4)) end """ Implicit Euler Method """ function constructImplicitEuler(T::Type = Float64) A = Matrix{T}(undef,1,1) A[1] = 1 c = [1] α = [1] A = map(T,A) α = map(T,α) c = map(T,c) return(ImplicitRKTableau(A,c,α,1)) end """ Order 2 Midpoint Method """ function constructMidpointRule(T::Type = Float64) A = Matrix{T}(undef,1,1) A[1] = 1//2 c = [1//2] α = [1] A = map(T,A) α = map(T,α) c = map(T,c) return(ImplicitRKTableau(A,c,α,2)) end """ Order 2 Trapezoidal Rule (LobattoIIIA2) """ function constructTrapezoidalRule(T::Type = Float64) A = [0 0 1//2 1//2] c = [0;1] α = [1//2;1//2] αEEst = [1;0] A = map(T,A) α = map(T,α) c = map(T,c) αEEst = map(T,αEEst) return(ImplicitRKTableau(A,c,α,2,αEEst=αEEst,adaptiveorder=1)) end """ LobattoIIIA Order 4 method """ function constructLobattoIIIA4(T::Type = Float64) A = [0 0 0 5//24 1//3 -1//24 1//6 2//3 1//6] c = [0;1//2;1] α = [1//6;2//3;1//6] αEEst = [-1//2;2;-1//2] A = map(T,A) α = map(T,α) c = map(T,c) αEEst = map(T,αEEst) return(ImplicitRKTableau(A,c,α,4,αEEst=αEEst,adaptiveorder=3)) end """ LobattoIIIB Order 2 method """ function constructLobattoIIIB2(T::Type = Float64) A = [1//2 0 1//2 0] c = [0;1] α = [1//2;1//2] αEEst = [1;0] A = map(T,A) α = map(T,α) c = map(T,c) αEEst = map(T,αEEst) return(ImplicitRKTableau(A,c,α,2,αEEst=αEEst,adaptiveorder=1)) end """ LobattoIIIB Order 4 method """ function constructLobattoIIIB4(T::Type = Float64) A = [1//6 -1//6 0 1//6 1//3 0 1//6 5//6 0] c = [0;1//2;1] α = [1//6;2//3;1//6] αEEst = [-1//2;2;-1//2] A = map(T,A) α = map(T,α) c = map(T,c) αEEst = map(T,αEEst) return(ImplicitRKTableau(A,c,α,4,αEEst=αEEst,adaptiveorder=3)) end """ LobattoIIIC Order 2 method """ function constructLobattoIIIC2(T::Type = Float64) A = [1//2 -1//2 1//2 1//2] c = [0;1] α = [1//2;1//2] αEEst = [1;0] A = map(T,A) α = map(T,α) c = map(T,c) αEEst = map(T,αEEst) return(ImplicitRKTableau(A,c,α,2,αEEst=αEEst,adaptiveorder=1)) end """ LobattoIIIC Order 4 method """ function constructLobattoIIIC4(T::Type = Float64) A = [1//6 -1//3 1//6 1//6 5//12 -1//12 1//6 2//3 1//6] c = [0;1//2;1] α = [1//6;2//3;1//6] αEEst = [-1//2;2;-1//2] A = map(T,A) α = map(T,α) c = map(T,c) αEEst = map(T,αEEst) return(ImplicitRKTableau(A,c,α,4,αEEst=αEEst,adaptiveorder=3)) end """ LobattoIIIC* Order 2 method """ function constructLobattoIIICStar2(T::Type = Float64) A = [0 0 1 0] c = [0;1] α = [1//2;1//2] A = map(T,A) α = map(T,α) c = map(T,c) return(ImplicitRKTableau(A,c,α,2)) end """ LobattoIIIC* Order 4 method """ function constructLobattoIIICStar4(T::Type = Float64) A = [0 0 0 1//4 1//4 0 0 1 0] c = [0;1//2;1] α = [1//6;2//3;1//6] αEEst = [-1//2;2;-1//2] A = map(T,A) α = map(T,α) c = map(T,c) αEEst = map(T,αEEst) return(ImplicitRKTableau(A,c,α,4,αEEst=αEEst,adaptiveorder=3)) end """ LobattoIIID Order 2 method """ function constructLobattoIIID2(T::Type = Float64) A = [1//2 1//2 -1//2 1//2] c = [0;1] α = [1//2;1//2] A = map(T,A) α = map(T,α) c = map(T,c) return(ImplicitRKTableau(A,c,α,2)) end """ LobattoIIID Order 4 method """ function constructLobattoIIID4(T::Type = Float64) A = [1//6 0 -1//6 1//12 5//12 0 1//2 1//3 1//6] c = [0;1//2;1] α = [1//6;2//3;1//6] αEEst = [-1//2;2;-1//2] A = map(T,A) α = map(T,α) c = map(T,c) αEEst = map(T,αEEst) return(ImplicitRKTableau(A,c,α,4,αEEst=αEEst,adaptiveorder=3)) end """ RadauIA Order 3 method """ function constructRadauIA3(T::Type = Float64) A = [1//4 -1//4 1//4 5//12] c = [0;2//3] α = [1//4;3//4] A = map(T,A) α = map(T,α) c = map(T,c) return(ImplicitRKTableau(A,c,α,3)) end """ RadauIA Order 5 method """ function constructRadauIA5(T::Type = Float64) A = [1//9 (-1-sqrt(6))/18 (-1+sqrt(6))/18 1//9 11//45+7*sqrt(6)/360 11//45-43*sqrt(6)/360 1//9 11//45+43*sqrt(6)/360 11//45-7*sqrt(6)/360] c = [0;3//5-sqrt(6)/10;3//5+sqrt(6)/10] α = [1//9;4//9+sqrt(6)/36;4//9-sqrt(6)/36] A = map(T,A) α = map(T,α) c = map(T,c) return(ImplicitRKTableau(A,c,α,5)) end """ RadauIIA Order 3 method """ function constructRadauIIA3(T::Type = Float64) A = [5//12 -1//12 3//4 1//4] c = [1//3;1] α = [3//4;1//4] A = map(T,A) α = map(T,α) c = map(T,c) return(ImplicitRKTableau(A,c,α,3)) end """ RadauIIA Order 5 method """ function constructRadauIIA5(T::Type = Float64) sq6 = sqrt(6) A = [11//45-7sq6/360 37//225-169sq6/1800 -2//225+sq6/75 37//225+169sq6/1800 11//45+7sq6/360 -2//225-sq6/75 4//9-sq6/36 4//9+sq6/36 1//9] c = [2//5-sq6/10;2/5+sq6/10;1] α = [4//9-sq6/36;4//9+sq6/36;1//9] A = map(T,A) α = map(T,α) c = map(T,c) return(ImplicitRKTableau(A,c,α,5)) end """ Runge-Kutta-Fehlberg Order 4/5 method. """ function constructRKF5(T::Type = Float64) A = [0 0 0 0 0 0 1//4 0 0 0 0 0 3//32 9//32 0 0 0 0 1932//2197 -7200//2197 7296//2197 0 0 0 439//216 -8 3680//513 -845//4104 0 0 -8//27 2 -3544//2565 1859//4104 -11//40 0] c = [0;1//4;3//8;12//13;1;1//2] α = [16//135;0;6656//12825;28561//56430;-9//50;2//55] αEEst = [25//216;0;1408//2565;2197//4104;-1//5;0] A = map(T,A) α = map(T,α) c = map(T,c) αEEst = map(T,αEEst) return(ExplicitRKTableau(A,c,α,5,αEEst=αEEst,adaptiveorder=4)) end """ Runge's First Order 5 method """ function constructRungeFirst5(T::Type=Float64) A = zeros(T,6,6) c = zeros(T,6) α = zeros(T,6) c[2]=1//5 c[3]=2//5 c[4]=1 c[5]=3//5 c[6]=4//5 A[2,1]=1//5 A[3,1]=0 A[3,2]=2//5 A[4,1]=9//4 A[4,2]=-5 A[4,3]=15//4 A[5,1]=-63//100 A[5,2]=9//5 A[5,3]=-13//20 A[5,4]=2//25 A[6,1]=-6//25 A[6,2]=4//5 A[6,3]=2//15 A[6,4]=8//75 A[6,5]=0 α[1]=17//144 α[2]=0 α[3]=25//36 α[4]=1//72 α[5]=-25//72 α[6]=25//48 A = map(T,A) α = map(T,α) c = map(T,c) return(ExplicitRKTableau(A,c,α,5)) end """ Cassity's Order 5 method """ function constructCassity5(T::Type=Float64) A = zeros(T,6,6) c = zeros(T,6) α = zeros(T,6) c[2]=1//7 c[3]=5//14 c[4]=9//14 c[5]=6//7 c[6]=1 A[2,1]=1//7 A[3,1]=-367//4088 A[3,2]=261//584 A[4,1]=41991//2044 A[4,2]=-2493//73 A[4,3]=57//4 A[5,1]=-108413//196224 A[5,2]=58865//65408 A[5,3]=5//16 A[5,4]=265//1344 A[6,1]=-204419//58984 A[6,2]=143829//58984 A[6,3]=171//202 A[6,4]=2205//404 A[6,5]=-432//101 α[1]=1//9 α[2]=7//2700 α[3]=413//810 α[4]=7//450 α[5]=28//75 α[6]=-101//8100 A = map(T,A) α = map(T,α) c = map(T,c) return(ExplicitRKTableau(A,c,α,5)) end """ Lawson's 5th order scheme An Order Five Runge Kutta Process with Extended Region of Stability, J. Douglas Lawson, Siam Journal on Numerical Analysis, Vol. 3, No. 4, (Dec., 1966) pages 593-597 """ function constructLawson5(T::Type=Float64) A = zeros(T,6,6) c = zeros(T,6) α = zeros(T,6) c[2]=1//12 c[3]=1//4 c[4]=1//2 c[5]=3//4 c[6]=1 A[2,1]=1//12 A[3,1]=-1//8 A[3,2]=3//8 A[4,1]=3//5 A[4,2]=-9//10 A[4,3]=4//5 A[5,1]=39//80 A[5,2]=-9//20 A[5,3]=3//20 A[5,4]=9//16 A[6,1]=-59//35 A[6,2]=66//35 A[6,3]=48//35 A[6,4]=-12//7 A[6,5]=8//7 α[1]=7//90 α[2]=0 α[3]=16//45 α[4]=2//15 α[5]=16//45 α[6]=7//90 A = map(T,A) α = map(T,α) c = map(T,c) return(ExplicitRKTableau(A,c,α,5)) end """ Luther and Konen's First Order 5 Some Fifth-Order Classical Runge Kutta Formulas, H.A.Luther and H.P.Konen, Siam Review, Vol. 3, No. 7, (Oct., 1965) pages 551-558. """ function constructLutherKonen5(T::Type = Float64) A = zeros(T,6,6) c = zeros(T,6) α = zeros(T,6) c[2]=1/2 c[3]=1/2-1/10*5^(1/2) c[4]=1/2 c[5]=1/2+1/10*5^(1/2) c[6]=1 A[2,1]=1/2 A[3,1]=1/5 A[3,2]=3/10-1/10*5^(1/2) A[4,1]=1/4 A[4,2]=1/4 A[4,3]=0 A[5,1]=1/20-1/20*5^(1/2) A[5,2]=-1/5 A[5,3]=1/4+3/20*5^(1/2) A[5,4]=2/5 A[6,1]=1/4*5^(1/2)-1/4 A[6,2]=1/2*5^(1/2)-1/2 A[6,3]=5/4-1/4*5^(1/2) A[6,4]=-2 A[6,5]=5/2-1/2*5^(1/2) α[1]=1/12 α[2]=0 α[3]=5/12 α[4]=0 α[5]=5/12 α[6]=1/12 A = map(T,A) α = map(T,α) c = map(T,c) return(ExplicitRKTableau(A,c,α,5)) end """ Luther and Konen's Second Order 5 Some Fifth-Order Classical Runge Kutta Formulas, H.A.Luther and H.P.Konen, Siam Review, Vol. 3, No. 7, (Oct., 1965) pages 551-558. """ function constructLutherKonen52(T::Type = Float64) A = zeros(T,6,6) c = zeros(T,6) α = zeros(T,6) c[2]=2/5 c[3]=1/2 c[4]=1 c[5]=1/2-1/10*15^(1/2) c[6]=1/2+1/10*15^(1/2) A[2,1]=2/5 A[3,1]=3/16 A[3,2]=5/16 A[4,1]=1/4 A[4,2]=-5/4 A[4,3]=2 A[5,1]=3/20-1/100*15^(1/2) A[5,2]=-1/4 A[5,3]=3/5-2/25*15^(1/2) A[5,4]=-1/100*15^(1/2) A[6,1]=-3/20-1/20*15^(1/2) A[6,2]=-1/4 A[6,3]=3/5 A[6,4]=3/10-1/20*15^(1/2) A[6,5]=1/5*15^(1/2) α[1]=0 α[2]=0 α[3]=4/9 α[4]=0 α[5]=5/18 α[6]=5/18 A = map(T,A) α = map(T,α) c = map(T,c) return(ExplicitRKTableau(A,c,α,5)) end """ Luther and Konen's Third Order 5 Some Fifth-Order Classical Runge Kutta Formulas, H.A.Luther and H.P.Konen, Siam Review, Vol. 3, No. 7, (Oct., 1965) pages 551-558. """ function constructLutherKonen53(T::Type = Float64) A = zeros(T,6,6) c = zeros(T,6) α = zeros(T,6) c[2]=3/25 c[3]=5/18 c[4]=45/89 c[5]=3/5-1/10*6^(1/2) c[6]=3/5+1/10*6^(1/2) A[2,1]=3/25 A[3,1]=-85/1944 A[3,2]=625/1944 A[4,1]=610425/1409938 A[4,2]=-961875/1409938 A[4,3]=532170/704969 A[5,1]=7411/37500-673/18750*6^(1/2) A[5,2]=0 A[5,3]=27621/228125*6^(1/2)-6561/228125 A[5,4]=1180229/2737500-126736/684375*6^(1/2) A[6,1]=-5351/62500-7087/281250*6^(1/2) A[6,2]=0 A[6,3]=2736423/1140625+786753/1140625*6^(1/2) A[6,4]=73736589/86687500+101816534/195046875*6^(1/2) A[6,5]=-30448/11875-12903/11875*6^(1/2) α[1]=1/9 α[2]=0 α[3]=0 α[4]=0 α[5]=4/9+1/36*6^(1/2) α[6]=4/9-1/36*6^(1/2) A = map(T,A) α = map(T,α) c = map(T,c) return(ExplicitRKTableau(A,c,α,5)) end """ S.N. Papakostas and G. PapaGeorgiou higher error more stable A Family of Fifth-order Runge-Kutta Pairs, by S.N. Papakostas and G. PapaGeorgiou, Mathematics of Computation,Volume 65, Number 215, July 1996, Pages 1165-1181. """ function constructPapakostasPapaGeorgiou5(T::Type = Float64) A = zeros(T,7,7) c = zeros(T,7) α = zeros(T,7) αEEst = zeros(T,7) c[2]=64//315 c[3]=115//381 c[4]=762//935 c[5]=25//28 c[6]=1 c[7]=1 A[2,1]=64//315 A[3,1]=480815//6193536 A[3,2]=462875//2064512 A[4,1]=344904825219069//345923838700000 A[4,2]=-2360077908867//601606676000 A[4,3]=40439332863108//10810119959375 A[5,1]=12078745127989699//5009699157786624 A[5,2]=-791781731775//81669668864 A[5,3]=39297175833216951//4694461413969824 A[5,4]=-10508413393960625//54097233987826176 A[6,1]=2251421737440863//828701767536000 A[6,2]=-39895842357//3782730880 A[6,3]=34564628685305112534//3916944972468643375 A[6,4]=12051135495733565//36492943960723917 A[6,5]=-26808346215168//82592376030125 A[7,1]=2405713//26289000 A[7,2]=0 A[7,3]=63896466577779//141024193000600 A[7,4]=454128848141375//589615117674696 A[7,5]=-1359311744//2892576375 A[7,6]=256979//1656648 α[1]=2405713//26289000 α[2]=0 α[3]=63896466577779//141024193000600 α[4]=454128848141375//589615117674696 α[5]=-1359311744//2892576375 α[6]=256979//1656648 αEEst[1]=1818563883019//20194131951000 αEEst[2]=0 αEEst[3]=5513862498202899713//12036555896794210600 αEEst[4]=324806515311046773125//452918159177876804664 αEEst[5]=-126112324722496//317422653663375 αEEst[6]=137695258717//1272569071032 αEEst[7]=1//42 A = map(T,A) α = map(T,α) c = map(T,c) αEEst = map(T,αEEst) return(ExplicitRKTableau(A,c,α,5,αEEst=αEEst,adaptiveorder=4,fsal=true)) end """ S.N. Papakostas and G. PapaGeorgiou less stable lower error Strictly better than DP5 A Family of Fifth-order Runge-Kutta Pairs, by S.N. Papakostas and G. PapaGeorgiou, Mathematics of Computation,Volume 65, Number 215, July 1996, Pages 1165-1181. """ function constructPapakostasPapaGeorgiou52(T::Type = Float64) A = zeros(T,7,7) c = zeros(T,7) α = zeros(T,7) αEEst = zeros(T,7) c[2]=35//159 c[3]=42//131 c[4]=131//143 c[5]=21//22 c[6]=1 c[7]=1 A[2,1]=35//159 A[3,1]=7476//85805 A[3,2]=20034//85805 A[4,1]=2438549411//1983961980 A[4,2]=-3707256508//716430715 A[4,3]=25077455105//5158301148 A[5,1]=105337889067//64388030080 A[5,2]=-1698584121//245755840 A[5,3]=6869523776931//1096562558080 A[5,4]=-927215289//26981535520 A[6,1]=67512025387//32454592380 A[6,2]=-20051384//2293935 A[6,3]=10587214001321//1373901639516 A[6,4]=731293420//8209319229 A[6,5]=-144610048//1077690663 A[7,1]=669707//6932520 A[7,2]=0 A[7,3]=2215522905683//4570867891800 A[7,4]=349043981//116904400 A[7,5]=-2234144//575505 A[7,6]=9363//7120 α[1]=669707//6932520 α[2]=0 α[3]=2215522905683//4570867891800 α[4]=349043981//116904400 α[5]=-2234144//575505 α[6]=9363//7120 αEEst[1]=2243660497//23535905400 αEEst[2]=0 αEEst[3]=7589131232781673//15518096492661000 αEEst[4]=1104461697911//396890438000 αEEst[5]=-6925033984//1953839475 αEEst[6]=3529851//3021550 αEEst[7]=1//112 A = map(T,A) α = map(T,α) c = map(T,c) αEEst = map(T,αEEst) return(ExplicitRKTableau(A,c,α,5,αEEst=αEEst,adaptiveorder=4,fsal=true)) end """ Runge–Kutta pairs of orders 5(4) using the minimal set of simplifying assumptions, by Ch. Tsitouras, TEI of Chalkis, Dept. of Applied Sciences, GR34400, Psahna, Greece. """ function constructTsitouras5(T::Type = Float64) A = zeros(T,7,7) c = zeros(T,7) α = zeros(T,7) αEEst = zeros(T,7) c[2] = convert(T,161//1000) c[3] = convert(T,327//1000) c[4] = convert(T,9//10) c[5] = convert(T,big".9800255409045096857298102862870245954942137979563024768854764293221195950761080302604") c[6] = convert(T,1) c[7] = convert(T,1) A[2,1] = convert(T,161//1000) A[3,1] = convert(T,big"-.8480655492356988544426874250230774675121177393430391537369234245294192976164141156943e-2") A[3,2] = convert(T,big".3354806554923569885444268742502307746751211773934303915373692342452941929761641411569") A[4,1] = convert(T,big"2.897153057105493432130432594192938764924887287701866490314866693455023795137503079289") A[4,2] = convert(T,big"-6.359448489975074843148159912383825625952700647415626703305928850207288721235210244366") A[4,3] = convert(T,big"4.362295432869581411017727318190886861027813359713760212991062156752264926097707165077") A[5,1] = convert(T,big"5.325864828439256604428877920840511317836476253097040101202360397727981648835607691791") A[5,2] = convert(T,big"-11.74888356406282787774717033978577296188744178259862899288666928009020615663593781589") A[5,3] = convert(T,big"7.495539342889836208304604784564358155658679161518186721010132816213648793440552049753") A[5,4] = convert(T,big"-.9249506636175524925650207933207191611349983406029535244034750452930469056411389539635e-1") A[6,1] = convert(T,big"5.861455442946420028659251486982647890394337666164814434818157239052507339770711679748") A[6,2] = convert(T,big"-12.92096931784710929170611868178335939541780751955743459166312250439928519268343184452") A[6,3] = convert(T,big"8.159367898576158643180400794539253485181918321135053305748355423955009222648673734986") A[6,4] = convert(T,big"-.7158497328140099722453054252582973869127213147363544882721139659546372402303777878835e-1") A[6,5] = convert(T,big"-.2826905039406838290900305721271224146717633626879770007617876201276764571291579142206e-1") A[7,1] = convert(T,big".9646076681806522951816731316512876333711995238157997181903319145764851595234062815396e-1") A[7,2] = convert(T,1//100) A[7,3] = convert(T,big".4798896504144995747752495322905965199130404621990332488332634944254542060153074523509") A[7,4] = convert(T,big"1.379008574103741893192274821856872770756462643091360525934940067397245698027561293331") A[7,5] = convert(T,big"-3.290069515436080679901047585711363850115683290894936158531296799594813811049925401677") A[7,6] = convert(T,big"2.324710524099773982415355918398765796109060233222962411944060046314465391054716027841") α[1] = convert(T,big".9646076681806522951816731316512876333711995238157997181903319145764851595234062815396e-1") α[2] = convert(T,1//100) α[3] = convert(T,big".4798896504144995747752495322905965199130404621990332488332634944254542060153074523509") α[4] = convert(T,big"1.379008574103741893192274821856872770756462643091360525934940067397245698027561293331") α[5] = convert(T,big"-3.290069515436080679901047585711363850115683290894936158531296799594813811049925401677") α[6] = convert(T,big"2.324710524099773982415355918398765796109060233222962411944060046314465391054716027841") αEEst[1]= convert(T,big".9468075576583945807478876255758922856117527357724631226139574065785592789071067303271e-1") αEEst[2]= convert(T,big".9183565540343253096776363936645313759813746240984095238905939532922955247253608687270e-2") αEEst[3]= convert(T,big".4877705284247615707855642599631228241516691959761363774365216240304071651579571959813") αEEst[4]= convert(T,big"1.234297566930478985655109673884237654035539930748192848315425833500484878378061439761") αEEst[5]= convert(T,big"-2.707712349983525454881109975059321670689605166938197378763992255714444407154902012702") αEEst[6]= convert(T,big"1.866628418170587035753719399566211498666255505244122593996591602841258328965767580089") αEEst[7]= convert(T,1//66) A = map(T,A) α = map(T,α) c = map(T,c) αEEst = map(T,αEEst) return(ExplicitRKTableau(A,c,α,5,αEEst=αEEst,adaptiveorder=4,fsal=true)) end """ An Efficient Runge-Kutta (4,5) Pair by P.Bogacki and L.F.Shampine Computers and Mathematics with Applications, Vol. 32, No. 6, 1996, pages 15 to 28 """ function constructBogakiShampine5(T::Type = Float64) A = zeros(T,8,8) c = zeros(T,8) α = zeros(T,8) αEEst = zeros(T,8) αEEst2 = zeros(T,8) c[2]=1//6 c[3]=2//9 c[4]=3//7 c[5]=2//3 c[6]=3//4 c[7]=1 c[8]=1 A[2,1]=1//6 A[3,1]=2//27 A[3,2]=4//27 A[4,1]=183//1372 A[4,2]=-162//343 A[4,3]=1053//1372 A[5,1]=68//297 A[5,2]=-4//11 A[5,3]=42//143 A[5,4]=1960//3861 A[6,1]=597//22528 A[6,2]=81//352 A[6,3]=63099//585728 A[6,4]=58653//366080 A[6,5]=4617//20480 A[7,1]=174197//959244 A[7,2]=-30942//79937 A[7,3]=8152137//19744439 A[7,4]=666106//1039181 A[7,5]=-29421//29068 A[7,6]=482048//414219 A[8,1]=587//8064 A[8,2]=0 A[8,3]=4440339//15491840 A[8,4]=24353//124800 A[8,5]=387//44800 A[8,6]=2152//5985 A[8,7]=7267//94080 α[1]=587//8064 α[2]=0 α[3]=4440339//15491840 α[4]=24353//124800 α[5]=387//44800 α[6]=2152//5985 α[7]=7267//94080 α[8]=0 αEEst[1]=6059//80640 αEEst[2]=0 αEEst[3]=8559189//30983680 αEEst[4]=26411//124800 αEEst[5]=-927//89600 αEEst[6]=443//1197 αEEst[7]=7267//94080 αEEst2[1]=2479//34992 αEEst2[2]=0 αEEst2[3]=123//416 αEEst2[4]=612941//3411720 αEEst2[5]=43//1440 αEEst2[6]=2272//6561 αEEst2[7]=79937//1113912 αEEst2[8]=3293//556956 return(ExplicitRKTableau(A,c,α,5,αEEst=αEEst,adaptiveorder=4)) end """ Explicit Runge-Kutta Pairs with One More Derivative Evaluation than the Minimum, by P.W.Sharp and E.Smart, Siam Journal of Scientific Computing, Vol. 14, No. 2, pages. 338-348, March 1993. """ function constructSharpSmart5(T::Type = Float64) A = zeros(T,7,7) c = zeros(T,7) α = zeros(T,7) αEEst = zeros(T,7) αEEst2 = zeros(T,7) c[2]=16//105 c[3]=8//35 c[4]=9//20 c[5]=2//3 c[6]=7//9 c[7]=1 A[2,1]=16//105 A[3,1]=2//35 A[3,2]=6//35 A[4,1]=8793//40960 A[4,2]=-5103//8192 A[4,3]=17577//20480 A[5,1]=347//1458 A[5,2]=-7//20 A[5,3]=3395//10044 A[5,4]=49792//112995 A[6,1]=-1223224109959//9199771214400 A[6,2]=1234787701//2523942720 A[6,3]=568994101921//3168810084960 A[6,4]=-105209683888//891227836395 A[6,5]=9//25 A[7,1]=2462504862877//8306031988800 A[7,2]=-123991//287040 A[7,3]=106522578491//408709510560 A[7,4]=590616498832//804646848915 A[7,5]=-319138726//534081275 A[7,6]=52758//71449 α[1]=1093//15120 α[2]=0 α[3]=60025//190992 α[4]=3200//20709 α[5]=1611//11960 α[6]=712233//2857960 α[7]=3//40 αEEst[1]=84018211//991368000 αEEst[2]=0 αEEst[3]=92098979//357791680 αEEst[4]=17606944//67891005 αEEst[5]=3142101//235253200 αEEst[6]=22004596809//70270091500 αEEst[7]=9//125 A = map(T,A) α = map(T,α) c = map(T,c) αEEst = map(T,αEEst) return(ExplicitRKTableau(A,c,α,5,αEEst=αEEst,adaptiveorder=4)) end """ constructBogakiShampine3() Constructs the tableau object for the Bogakai-Shampine Order 2/3 method. """ function constructBogakiShampine3(T::Type = Float64) A = [0 0 0 0 1//2 0 0 0 0 3//4 0 0 2//9 1//3 4//9 0] c = [0;1//2;3//4;1] α = [2//9;1//3;4//9;0] αEEst = [7//24;1//4;1//3;1//8] A = map(T,A) α = map(T,α) c = map(T,c) αEEst = map(T,αEEst) return(ExplicitRKTableau(A,c,α,3,αEEst=αEEst,adaptiveorder=2)) end """ constructCashKarp() Constructs the tableau object for the Cash-Karp Order 4/5 method. """ function constructCashKarp(T::Type = Float64) A = [0 0 0 0 0 0 1//5 0 0 0 0 0 3//40 9//40 0 0 0 0 3//10 -9//10 6//5 0 0 0 -11//54 5//2 -70//27 35//27 0 0 1631//55296 175//512 575//13824 44275//110592 253//4096 0] c = [0;1//5;3//10;3//5;1;7//8] α = [37//378;0;250//621;125//594;0;512//1771] αEEst = [2825//27648;0;18574//48384;13525//554296;277//14336;1//4] A = map(T,A) α = map(T,α) c = map(T,c) αEEst = map(T,αEEst) return(ExplicitRKTableau(A,c,α,5,αEEst=αEEst,adaptiveorder=4)) end """ Runge-Kutta-Fehberg Order 4/3 """ function constructRKF4(T::Type = Float64) c = [0;1//4;4//9;6//7;1] A = [0 0 0 0 0 1//4 0 0 0 0 4//81 32//81 0 0 0 57//98 -432//343 1053//686 0 0 1//6 0 27//52 49//156 0 ] α = [43//288;0;243//416;343//1872;1//12] αEEst = [1//6;0;27//52;49//156;0] A = map(T,A) α = map(T,α) c = map(T,c) αEEst = map(T,αEEst) return(ExplicitRKTableau(A,c,α,4,αEEst=αEEst,adaptiveorder=3)) end """ Butcher's Third Order 6 On Runge-Kutta Processes of High Order, by J. C. Butcher, Journal of the Australian Mathematical Society, Vol. 4, (1964), pages 179 to 194 """ function constructButcher63(T::Type = Float64) c = [0;1/2;2/3;1/3;5/6;1/6;1] A = [0 0 0 0 0 0 0 1//2 0 0 0 0 0 0 2//9 4//9 0 0 0 0 0 7//36 2//9 -1//12 0 0 0 0 -35//144 -55//36 35//48 15//8 0 0 0 -1//360 -11//36 -1//8 1//2 1//10 0 0 -41//260 22//13 43//156 -118//39 32//195 80//39 0 ] α = [13//200;0;11//40;11//40;4//25;4//25;13//200] A = map(T,A) α = map(T,α) c = map(T,c) return(ExplicitRKTableau(A,c,α,6)) end """ Butcher's First Order 6 method On Runge-Kutta Processes of High Order, by J. C. Butcher, Journal of the Australian Mathematical Society, Vol. 4, (1964), pages 179 to 194 """ function constructButcher6(T::Type = Float64) A = zeros(T,7,7) c = zeros(T,7) α = zeros(T,7) c[2]=1//2-1//10*5^(1//2) c[3]=1//2+1//10*5^(1//2) c[4]=1//2-1//10*5^(1//2) c[5]=1//2+1//10*5^(1//2) c[6]=1//2-1//10*5^(1//2) c[7]=1 A[2,1]=1//2-1//10*5^(1//2) A[3,1]=-1//10*5^(1//2) A[3,2]=1//2+1//5*5^(1//2) A[4,1]=7//20*5^(1//2)-3//4 A[4,2]=1//4*5^(1//2)-1//4 A[4,3]=3//2-7//10*5^(1//2) A[5,1]=1//12-1//60*5^(1//2) A[5,2]=0 A[5,3]=1//6 A[5,4]=7//60*5^(1//2)+1//4 A[6,1]=1//60*5^(1//2)+1//12 A[6,2]=0 A[6,3]=3//4-5//12*5^(1//2) A[6,4]=1//6 A[6,5]=-1//2+3//10*5^(1//2) A[7,1]=1//6 A[7,2]=0 A[7,3]=-55//12+25//12*5^(1//2) A[7,4]=-7//12*5^(1//2)-25//12 A[7,5]=5-2*5^(1//2) A[7,6]=5//2+1//2*5^(1//2) α[1]=1//12 α[2]=0 α[3]=0 α[4]=0 α[5]=5//12 α[6]=5//12 α[7]=1//12 A = map(T,A) α = map(T,α) c = map(T,c) return(ExplicitRKTableau(A,c,α,6)) end """ Butcher's Second Order 6 method On Runge-Kutta Processes of High Order, by J. C. Butcher, Journal of the Australian Mathematical Society, Vol. 4, (1964), pages 179 to 194 """ function constructButcher62(T::Type = Float64) A = zeros(T,7,7) c = zeros(T,7) α = zeros(T,7) c[2]=1//3 c[3]=2//3 c[4]=1//3 c[5]=1//2 c[6]=1//2 c[7]=1 A[2,1]=1//3 A[3,1]=0 A[3,2]=2//3 A[4,1]=1//12 A[4,2]=1//3 A[4,3]=-1//12 A[5,1]=-1//16 A[5,2]=9//8 A[5,3]=-3//16 A[5,4]=-3//8 A[6,1]=0 A[6,2]=9//8 A[6,3]=-3//8 A[6,4]=-3//4 A[6,5]=1//2 A[7,1]=9//44 A[7,2]=-9//11 A[7,3]=63//44 A[7,4]=18//11 A[7,5]=0 A[7,6]=-16//11 α[1]=11//120 α[2]=0 α[3]=27//40 α[4]=27//40 α[5]=-4//15 α[6]=-4//15 α[7]=11//120 A = map(T,A) α = map(T,α) c = map(T,c) return(ExplicitRKTableau(A,c,α,6)) end """ Verner Order 5/6 method A Contrast of a New RK56 pair with DP56, by Jim Verner, Department of Mathematics. Simon Fraser University, Burnaby, Canada, 2006. """ function constructVerner6(T::Type = Float64) A = zeros(T,8,8) c = zeros(T,8) α = zeros(T,8) αEEst = zeros(T,8) c[2]=1//7 c[3]=2//9 c[4]=3//7 c[5]=2//3 c[6]=3//4 c[7]=1 c[8]=1 A[2,1]=1//7 A[3,1]=4//81 A[3,2]=14//81 A[4,1]=291//1372 A[4,2]=-27//49 A[4,3]=1053//1372 A[5,1]=86//297 A[5,2]=-14//33 A[5,3]=42//143 A[5,4]=1960//3861 A[6,1]=-267//22528 A[6,2]=189//704 A[6,3]=63099//585728 A[6,4]=58653//366080 A[6,5]=4617//20480 A[7,1]=10949//6912 A[7,2]=-69//32 A[7,3]=-90891//68096 A[7,4]=112931//25920 A[7,5]=-69861//17920 A[7,6]=26378//10773 A[8,1]=1501//19008 A[8,2]=-21//88 A[8,3]=219519//347776 A[8,4]=163807//926640 A[8,5]=-417//640 A[8,6]=1544//1539 A[8,7]=0 α[1]=79//1080 α[2]=0 α[3]=19683//69160 α[4]=16807//84240 α[5]=0 α[6]=2816//7695 α[7]=1//100 α[8]=187//2800 αEEst[1]=763//10800 αEEst[2]=0 αEEst[3]=59049//197600 αEEst[4]=88837//526500 αEEst[5]=243//4000 αEEst[6]=12352//38475 αEEst[7]=0 αEEst[8]=2//25 A = map(T,A) α = map(T,α) c = map(T,c) αEEst = map(T,αEEst) return(ExplicitRKTableau(A,c,α,6,adaptiveorder=5,αEEst=αEEst)) end """ Dormand-Prince Order 5//6 method P.J. Prince and J. R. Dormand, High order embedded Runge-Kutta formulae, Journal of Computational and Applied Mathematics . 7 (1981), pp. 67-75. """ function constructDormandPrince6(T::Type = Float64) A = zeros(T,8,8) c = zeros(T,8) α = zeros(T,8) αEEst = zeros(T,8) c[2]=1//10 c[3]=2//9 c[4]=3//7 c[5]=3//5 c[6]=4//5 c[7]=1 c[8]=1 A[2,1]=1//10 A[3,1]=-2//81 A[3,2]=20//81 A[4,1]=615//1372 A[4,2]=-270//343 A[4,3]=1053//1372 A[5,1]=3243//5500 A[5,2]=-54//55 A[5,3]=50949//71500 A[5,4]=4998//17875 A[6,1]=-26492//37125 A[6,2]=72//55 A[6,3]=2808//23375 A[6,4]=-24206//37125 A[6,5]=338//459 A[7,1]=5561//2376 A[7,2]=-35//11 A[7,3]=-24117//31603 A[7,4]=899983//200772 A[7,5]=-5225//1836 A[7,6]=3925//4056 A[8,1]=465467//266112 A[8,2]=-2945//1232 A[8,3]=-5610201//14158144 A[8,4]=10513573//3212352 A[8,5]=-424325//205632 A[8,6]=376225//454272 A[8,7]=0 α[1]=61//864 α[2]=0 α[3]=98415//321776 α[4]=16807//146016 α[5]=1375//7344 α[6]=1375//5408 α[7]=-37//1120 α[8]=1//10 αEEst[1]=821//10800 αEEst[2]=0 αEEst[3]=19683//71825 αEEst[4]=175273//912600 αEEst[5]=395//3672 αEEst[6]=785//2704 αEEst[7]=3//50 αEEst[8]=0 A = map(T,A) α = map(T,α) c = map(T,c) αEEst = map(T,αEEst) return(ExplicitRKTableau(A,c,α,6,adaptiveorder=5,αEEst=αEEst)) end """ Sharp-Verner Order 5/6 method Completely Imbedded Runge-Kutta Pairs, by P. W. Sharp and J. H. Verner, SIAM Journal on Numerical Analysis, Vol. 31, No. 4. (Aug., 1994), pages. 1169 to 1190. """ function constructSharpVerner6(T::Type = Float64) A = zeros(T,9,9) c = zeros(T,9) α = zeros(T,9) αEEst = zeros(T,9) c[2]=1//12 c[3]=2//15 c[4]=1//5 c[5]=8//15 c[6]=2//3 c[7]=19//20 c[8]=1 c[9]=1 A[2,1]=1//12 A[3,1]=2//75 A[3,2]=8//75 A[4,1]=1//20 A[4,2]=0 A[4,3]=3//20 A[5,1]=88//135 A[5,2]=0 A[5,3]=-112//45 A[5,4]=64//27 A[6,1]=-10891//11556 A[6,2]=0 A[6,3]=3880//963 A[6,4]=-8456//2889 A[6,5]=217//428 A[7,1]=1718911//4382720 A[7,2]=0 A[7,3]=-1000749//547840 A[7,4]=819261//383488 A[7,5]=-671175//876544 A[7,6]=14535//14336 A[8,1]=85153//203300 A[8,2]=0 A[8,3]=-6783//2140 A[8,4]=10956//2675 A[8,5]=-38493//13375 A[8,6]=1152//425 A[8,7]=-7168//40375 A[9,1]=53//912 A[9,2]=0 A[9,3]=0 A[9,4]=5//16 A[9,5]=27//112 A[9,6]=27//136 A[9,7]=256//969 A[9,8]=-25//336 α[1]=53//912 α[2]=0 α[3]=0 α[4]=5//16 α[5]=27//112 α[6]=27//136 α[7]=256//969 α[8]=-25//336 αEEst[1]=617//10944 αEEst[2]=0 αEEst[3]=0 αEEst[4]=241//756 αEEst[5]=69//320 αEEst[6]=435//1904 αEEst[7]=10304//43605 αEEst[8]=0 αEEst[9]=-1//18 A = map(T,A) α = map(T,α) c = map(T,c) αEEst = map(T,αEEst) return(ExplicitRKTableau(A,c,α,6,adaptiveorder=5,αEEst=αEEst,fsal=true)) end """ Verner 1991 Second Order 5/6 method Some Ruge-Kutta Formula Pairs, by J.H.Verner, SIAM Journal on Numerical Analysis, Vol. 28, No. 2 (April 1991), pages 496 to 511. """ function constructVerner9162(T::Type = Float64) A = zeros(T,9,9) c = zeros(T,9) α = zeros(T,9) αEEst = zeros(T,9) c[2]=1//8 c[3]=1//6 c[4]=1//4 c[5]=1//2 c[6]=3//5 c[7]=4//5 c[8]=1 c[9]=1 A[2,1]=1//8 A[3,1]=1//18 A[3,2]=1//9 A[4,1]=1//16 A[4,2]=0 A[4,3]=3//16 A[5,1]=1//4 A[5,2]=0 A[5,3]=-3//4 A[5,4]=1 A[6,1]=134//625 A[6,2]=0 A[6,3]=-333//625 A[6,4]=476//625 A[6,5]=98//625 A[7,1]=-98//1875 A[7,2]=0 A[7,3]=12//625 A[7,4]=10736//13125 A[7,5]=-1936//1875 A[7,6]=22//21 A[8,1]=9//50 A[8,2]=0 A[8,3]=21//25 A[8,4]=-2924//1925 A[8,5]=74//25 A[8,6]=-15//7 A[8,7]=15//22 A[9,1]=11//144 A[9,2]=0 A[9,3]=0 A[9,4]=256//693 A[9,5]=0 A[9,6]=125//504 A[9,7]=125//528 A[9,8]=5//72 α[1]=11//144 α[2]=0 α[3]=0 α[4]=256//693 α[5]=0 α[6]=125//504 α[7]=125//528 α[8]=5//72 αEEst[1]=1//18 αEEst[2]=0 αEEst[3]=0 αEEst[4]=32//63 αEEst[5]=-2//3 αEEst[6]=125//126 αEEst[7]=0 αEEst[8]=-5//63 αEEst[9]=4//21 A = map(T,A) α = map(T,α) c = map(T,c) αEEst = map(T,αEEst) return(ExplicitRKTableau(A,c,α,6,adaptiveorder=5,αEEst=αEEst)) end """ Verner 1991 First Order 5/6 method Some Ruge-Kutta Formula Pairs, by J.H.Verner, SIAM Journal on Numerical Analysis, Vol. 28, No. 2 (April 1991), pages 496 to 511. """ function constructVerner916(T::Type = Float64) A = zeros(T,9,9) c = zeros(T,9) α = zeros(T,9) αEEst = zeros(T,9) c[2]=1//8 c[3]=4//9-4//45*10^(1//2) c[4]=2//3-2//15*10^(1//2) c[5]=9//16 c[6]=1//2 c[7]=9//10 c[8]=1 c[9]=1 A[2,1]=1//8 A[3,1]=-268//405+92//405*10^(1//2) A[3,2]=448//405-128//405*10^(1//2) A[4,1]=1//6-1//30*10^(1//2) A[4,2]=0 A[4,3]=1//2-1//10*10^(1//2) A[5,1]=11547//32768+405//16384*10^(1//2) A[5,2]=0 A[5,3]=-18225//32768-5103//16384*10^(1//2) A[5,4]=12555//16384+2349//8192*10^(1//2) A[6,1]=19662371//51149376+441281//12787344*10^(1//2) A[6,2]=0 A[6,3]=-3786045//5683264-252663//710408*10^(1//2) A[6,4]=1570556745//1821486112+290041461//910743056*10^(1//2) A[6,5]=-41227072//512292969+1374464//512292969*10^(1//2) A[7,1]=-154207593//369412160-1829424339//11544130000*10^(1//2) A[7,2]=0 A[7,3]=2659895739//1847060800+653855409//1154413000*10^(1//2) A[7,4]=-349492176711//591982986400-359784638379//1479957466000*10^(1//2) A[7,5]=153920585664//92497341625+311066673408//462486708125*10^(1//2) A[7,6]=-1944//1625-6804//8125*10^(1//2) A[8,1]=70594945601//21406013856+21473424323//21406013856*10^(1//2) A[8,2]=0 A[8,3]=-794525145//88090592-249156075//88090592*10^(1//2) A[8,4]=866290968775//254097312624+256998959765//254097312624*10^(1//2) A[8,5]=-15964196472448//1286367645159-5039429245312//1286367645159*10^(1//2) A[8,6]=17017//1116+5075//1116*10^(1//2) A[8,7]=42875//90396+16625//90396*10^(1//2) A[9,1]=31//324-37//4860*10^(1//2) A[9,2]=0 A[9,3]=0 A[9,4]=37435//69228-3235//69228*10^(1//2) A[9,5]=-1245184//1090341+9699328//16355115*10^(1//2) A[9,6]=71//54-74//135*10^(1//2) A[9,7]=625//486-250//729*10^(1//2) A[9,8]=-23//21+37//105*10^(1//2) α[1]=31//324-37//4860*10^(1//2) α[2]=0 α[3]=0 α[4]=37435//69228-3235//69228*10^(1//2) α[5]=-1245184//1090341+9699328//16355115*10^(1//2) α[6]=71//54-74//135*10^(1//2) α[7]=625//486-250//729*10^(1//2) α[8]=-23//21+37//105*10^(1//2) αEEst[1]=5//54-2//135*10^(1//2) αEEst[2]=0 αEEst[3]=0 αEEst[4]=2390//17307+2290//17307*10^(1//2) αEEst[5]=40960//121149+262144//605745*10^(1//2) αEEst[6]=2//27-64//135*10^(1//2) αEEst[7]=0 αEEst[8]=150029//443709-236267//2218545*10^(1//2) αEEst[9]=2411//126774+1921//63387*10^(1//2) A = map(T,A) α = map(T,α) c = map(T,c) αEEst = map(T,αEEst) return(ExplicitRKTableau(A,c,α,6,adaptiveorder=5,αEEst=αEEst)) end """ From Verner's Website """ function constructVernerRobust6(T::Type = Float64) A = zeros(T,9,9) c = zeros(T,9) α = zeros(T,9) αEEst = zeros(T,9) c[2]=9//50 c[3]=1//6 c[4]=1//4 c[5]=53//100 c[6]=3//5 c[7]=4//5 c[8]=1 c[9]=1 A[2,1]=9//50 A[3,1]= 29//324 A[3,2]= 25//324 A[4,1]= 1//16 A[4,2]=0 A[4,3]=3//16 A[5,1]= 79129//250000 A[5,2]=0 A[5,3]=-261237//250000 A[5,4]=19663//15625 A[6,1]= 1336883//4909125 A[6,2]=0 A[6,3]=-25476//30875 A[6,4]=194159//185250 A[6,5]= 8225//78546 A[7,1]=-2459386//14727375 A[7,2]=0 A[7,3]=19504//30875 A[7,4]=2377474//13615875 A[7,5]=-6157250//5773131 A[7,6]=902//735 A[8,1]=2699//7410 A[8,2]=0 A[8,3]=-252//1235 A[8,4]=-1393253//3993990 A[8,5]=236875//72618 A[8,6]=-135//49 A[8,7]=15//22 A[9,1]=11//144 A[9,2]=0 A[9,3]=0 A[9,4]=256//693 A[9,5]=0 A[9,6]=125//504 A[9,7]=125//528 A[9,8]=5//72 α[1]=11//144 α[2]=0 α[3]=0 α[4]=256//693 α[5]=0 α[6]=125//504 α[7]=125//528 α[8]=5//72 αEEst[1]=28//477 αEEst[2]=0 αEEst[3]=0 αEEst[4]=212//441 αEEst[5]=-312500//366177 αEEst[6]=2125//1764 αEEst[7]=0 αEEst[8]=-2105//35532 αEEst[9]=2995//17766 A = map(T,A) α = map(T,α) c = map(T,c) αEEst = map(T,αEEst) return(ExplicitRKTableau(A,c,α,6,adaptiveorder=5,αEEst=αEEst)) end """ From Verner's Website """ function constructVernerEfficient6(T::Type = Float64) A = zeros(T,9,9) c = zeros(T,9) α = zeros(T,9) αEEst = zeros(T,9) c[2]=3//50 c[3]=1439//15000 c[4]=1439//10000 c[5]=4973//10000 c[6]=389//400 c[7]=1999//2000 c[8]=1 c[9]=1 A[2,1]=3//50 A[3,1]=519479//27000000 A[3,2]=2070721//27000000 A[4,1]=1439//40000 A[4,2]=0 A[4,3]=4317//40000 A[5,1]=109225017611//82828840000 A[5,2]=0 A[5,3]=-417627820623//82828840000 A[5,4]=43699198143//10353605000 A[6,1]=-8036815292643907349452552172369//191934985946683241245914401600 A[6,2]=0 A[6,3]=246134619571490020064824665//1543816496655405117602368 A[6,4]=-13880495956885686234074067279//113663489566254201783474344 A[6,5]=755005057777788994734129//136485922925633667082436 A[7,1]=-BigInt(1663299841566102097180506666498880934230261)//BigInt(30558424506156170307020957791311384232000) A[7,2]=0 A[7,3]=130838124195285491799043628811093033//631862949514135618861563657970240 A[7,4]=-BigInt(3287100453856023634160618787153901962873)//BigInt(20724314915376755629135711026851409200) A[7,5]=2771826790140332140865242520369241//396438716042723436917079980147600 A[7,6]=-1799166916139193//96743806114007800 A[8,1]=-BigInt(832144750039369683895428386437986853923637763)//BigInt(15222974550069600748763651844667619945204887) A[8,2]=0 A[8,3]=818622075710363565982285196611368750//3936576237903728151856072395343129 A[8,4]=-BigInt(9818985165491658464841194581385463434793741875)//BigInt(61642597962658994069869370923196463581866011) A[8,5]=BigInt(31796692141848558720425711042548134769375)//BigInt(4530254033500045975557858016006308628092) A[8,6]=-14064542118843830075//766928748264306853644 A[8,7]=-1424670304836288125//2782839104764768088217 A[9,1]=382735282417//11129397249634 A[9,2]=0 A[9,3]=0 A[9,4]=5535620703125000//21434089949505429 A[9,5]=13867056347656250//32943296570459319 A[9,6]=626271188750//142160006043 A[9,7]=-51160788125000//289890548217 A[9,8]=163193540017//946795234 α[1]=382735282417//11129397249634 α[2]=0 α[3]=0 α[4]=5535620703125000//21434089949505429 α[5]=13867056347656250//32943296570459319 α[6]=626271188750//142160006043 α[7]=-51160788125000//289890548217 α[8]=163193540017//946795234 αEEst[1]=124310637869885675646798613//2890072468789466426596827670 αEEst[2]=0 αEEst[3]=0 αEEst[4]=265863151737164990361330921875//1113197271463372303940319369579 αEEst[5]=3075493557174030806536302953125//6843749922042323876546949699876 αEEst[6]=67798000008733879813263055//29532792147666737550036372 αEEst[7]=-1099436585155390846238326375//15055706496446408859196167 αEEst[8]=26171252653086373181571802//368794478890732346033505 αEEst[9]=1//30 A = map(T,A) α = map(T,α) c = map(T,c) αEEst = map(T,αEEst) return(ExplicitRKTableau(A,c,α,6,adaptiveorder=5,αEEst=αEEst,fsal=true)) end function constructVerner7(T::Type = Float64) A = zeros(T,10,10) c = zeros(T,10) α = zeros(T,10) αEEst = zeros(T,10) c[2] = 1//200 c[3] = 49//450 c[4] = 49//300 c[5] = 911//2000 c[6] = 3480084980//5709648941 c[7] = 221//250 c[8] = 37//40 c[9] = 1 c[10] = 1 A[2,1] = 1//200 A[3,1] = -4361//4050 A[3,2] = 2401//2025 A[4,1] = 49//1200 A[4,3] = 49//400 A[5,1] = 2454451729//3841600000 A[5,3] = -9433712007//3841600000 A[5,4] = 4364554539//1920800000 A[6,1] = -BigInt(6187101755456742839167388910402379177523537620)//BigInt(2324599620333464857202963610201679332423082271) A[6,3] = BigInt(27569888999279458303270493567994248533230000)//BigInt(2551701010245296220859455115479340650299761) A[6,4] = -BigInt(37368161901278864592027018689858091583238040000)//BigInt(4473131870960004275166624817435284159975481033) A[6,5] = BigInt(1392547243220807196190880383038194667840000000)//BigInt(1697219131380493083996999253929006193143549863) A[7,1] = 11272026205260557297236918526339//1857697188743815510261537500000 A[7,3] = -48265918242888069//1953194276993750 A[7,4] = 26726983360888651136155661781228//1308381343805114800955157615625 A[7,5] = -2090453318815827627666994432//1096684189897834170412307919 A[7,6] = BigInt(1148577938985388929671582486744843844943428041509)//BigInt(1141532118233823914568777901158338927629837500000) A[8,1] = BigInt(1304457204588839386329181466225966641)//BigInt(108211771565488329642169667802016000) A[8,3] = -1990261989751005//40001418792832 A[8,4] = BigInt(2392691599894847687194643439066780106875)//BigInt(58155654089143548047476915856270826016) A[8,5] = -BigInt(1870932273351008733802814881998561250)//BigInt(419326053051486744762255151208232123) A[8,6] = BigInt(1043329047173803328972823866240311074041739158858792987034783181)//BigInt(510851127745017966999893975119259285040213723744255237522144000) A[8,7] = -311918858557595100410788125//3171569057622789618800376448 A[9,1] = BigInt(17579784273699839132265404100877911157)//BigInt(1734023495717116205617154737841023480) A[9,3] = -18539365951217471064750//434776548575709731377 A[9,4] = BigInt(447448655912568142291911830292656995992000)//BigInt(12511202807447096607487664209063950964109) A[9,5] = -BigInt(65907597316483030274308429593905808000000)//BigInt(15158061430635748897861852383197382130691) A[9,6] = BigInt(273847823027445129865693702689010278588244606493753883568739168819449761)//BigInt(136252034448398939768371761610231099586032870552034688235302796640584360) A[9,7] = BigInt(694664732797172504668206847646718750)//BigInt(1991875650119463976442052358853258111) A[9,8] = -19705319055289176355560129234220800//72595753317320295604316217197876507 A[10,1] = -511858190895337044664743508805671//11367030248263048398341724647960 A[10,3] = 2822037469238841750//15064746656776439 A[10,4] = -BigInt(23523744880286194122061074624512868000)//BigInt(152723005449262599342117017051789699) A[10,5] = BigInt(10685036369693854448650967542704000000)//BigInt(575558095977344459903303055137999707) A[10,6] = -BigInt(6259648732772142303029374363607629515525848829303541906422993)//BigInt(876479353814142962817551241844706205620792843316435566420120) A[10,7] = 17380896627486168667542032602031250//13279937889697320236613879977356033 α[1] = 96762636172307789//2051985304794103980 α[4] = 312188947591288252500000//1212357694274963646019729 α[5] = 13550580884964304000000000000//51686919683339547115937980629 α[6] = BigInt(72367769693133178898676076432831566019684378142853445230956642801)//BigInt(475600216991873963561768100160364792981629064220601844848928537580) α[7] = 1619421054120605468750//3278200730370057108183 α[8] = -66898316144057728000//227310933007074849597 α[9] = 181081444637946577//2226845467039736466 αEEst[1] = 117807213929927//2640907728177740 αEEst[4] = 4758744518816629500000//17812069906509312711137 αEEst[5] = 1730775233574080000000000//7863520414322158392809673 αEEst[6] = BigInt(2682653613028767167314032381891560552585218935572349997)//BigInt(12258338284789875762081637252125169126464880985167722660) αEEst[7] = 40977117022675781250//178949401077111131341 αEEst[10] = 2152106665253777//106040260335225546 A = map(T,A) α = map(T,α) c = map(T,c) αEEst = map(T,αEEst) return(ExplicitRKTableau(A,c,α,7,adaptiveorder=6,αEEst=αEEst)) end """ Verner Efficient 8 """ function constructVerner8(T::Type = Float64) A = zeros(T,13,13) c = zeros(T,13) α = zeros(T,13) αEEst = zeros(T,13) c[2] = 1//20 c[3] = 341//3200 c[4] = 1023//6400 c[5] = 39//100 c[6] = 93//200 c[7] = 31//200 c[8] = 943//1000 c[9] = 7067558016280//7837150160667 c[10] = 909//1000 c[11] = 47//50 c[12] = 1 c[13] = 1 A[2,1] = 1//20 A[3,1] = -7161//1024000 A[3,2] = 116281//1024000 A[4,1] = 1023//25600 A[4,3] = 3069//25600 A[5,1] = 4202367//11628100 A[5,3] = -3899844//2907025 A[5,4] = 3982992//2907025 A[6,1] = 5611//114400 A[6,4] = 31744//135025 A[6,5] = 923521//5106400 A[7,1] = 21173//343200 A[7,4] = 8602624//76559175 A[7,5] = -26782109//689364000 A[7,6] = 5611//283500 A[8,1] = -1221101821869329//690812928000000 A[8,4] = -125//2 A[8,5] = -1024030607959889//168929280000000 A[8,6] = 1501408353528689//265697280000000 A[8,7] = 6070139212132283//92502016000000 A[9,1] = -BigInt(1472514264486215803881384708877264246346044433307094207829051978044531801133057155)//BigInt(1246894801620032001157059621643986024803301558393487900440453636168046069686436608) A[9,4] = -BigInt(5172294311085668458375175655246981230039025336933699114138315270772319372469280000)//BigInt(124619381004809145897278630571215298365257079410236252921850936749076487132995191) A[9,5] = -BigInt(12070679258469254807978936441733187949484571516120469966534514296406891652614970375)//BigInt(2722031154761657221710478184531100699497284085048389015085076961673446140398628096) A[9,6] = BigInt(780125155843893641323090552530431036567795592568497182701460674803126770111481625)//BigInt(183110425412731972197889874507158786859226102980861859505241443073629143100805376) A[9,7] = BigInt(664113122959911642134782135839106469928140328160577035357155340392950009492511875)//BigInt(15178465598586248136333023107295349175279765150089078301139943253016877823170816) A[9,8] = BigInt(10332848184452015604056836767286656859124007796970668046446015775000000)//BigInt(1312703550036033648073834248740727914537972028638950165249582733679393783) A[10,1] = -BigInt(29055573360337415088538618442231036441314060511)//BigInt(22674759891089577691327962602370597632000000000) A[10,4] = -20462749524591049105403365239069//454251913499893469596231268750 A[10,5] = -180269259803172281163724663224981097//38100922558256871086579832832000000 A[10,6] = BigInt(21127670214172802870128286992003940810655221489)//BigInt(4679473877997892906145822697976708633673728000) A[10,7] = BigInt(318607235173649312405151265849660869927653414425413)//BigInt(6714716715558965303132938072935465423910912000000) A[10,8] = 212083202434519082281842245535894//20022426044775672563822865371173879 A[10,9] = -BigInt(2698404929400842518721166485087129798562269848229517793703413951226714583)//BigInt(469545674913934315077000442080871141884676035902717550325616728175875000000) A[11,1] = -BigInt(2342659845814086836951207140065609179073838476242943917)//BigInt(1358480961351056777022231400139158760857532162795520000) A[11,4] = -996286030132538159613930889652//16353068885996164905464325675 A[11,5] = -26053085959256534152588089363841//4377552804565683061011299942400 A[11,6] = BigInt(20980822345096760292224086794978105312644533925634933539)//BigInt(3775889992007550803878727839115494641972212962174156800) A[11,7] = BigInt(890722993756379186418929622095833835264322635782294899)//BigInt(13921242001395112657501941955594013822830119803764736) A[11,8] = BigInt(161021426143124178389075121929246710833125)//BigInt(10997207722131034650667041364346422894371443) A[11,9] = BigInt(300760669768102517834232497565452434946672266195876496371874262392684852243925359864884962513)//BigInt(4655443337501346455585065336604505603760824779615521285751892810315680492364106674524398280000) A[11,10] = -31155237437111730665923206875//392862141594230515010338956291 A[12,1] = -BigInt(2866556991825663971778295329101033887534912787724034363)//BigInt(868226711619262703011213925016143612030669233795338240) A[12,4] = -BigInt(16957088714171468676387054358954754000)//BigInt(143690415119654683326368228101570221) A[12,5] = -BigInt(4583493974484572912949314673356033540575)//BigInt(451957703655250747157313034270335135744) A[12,6] = BigInt(2346305388553404258656258473446184419154740172519949575)//BigInt(256726716407895402892744978301151486254183185289662464) A[12,7] = BigInt(1657121559319846802171283690913610698586256573484808662625)//BigInt(13431480411255146477259155104956093505361644432088109056) A[12,8] = BigInt(345685379554677052215495825476969226377187500)//BigInt(74771167436930077221667203179551347546362089) A[12,9] = -BigInt(3205890962717072542791434312152727534008102774023210240571361570757249056167015230160352087048674542196011)//BigInt(947569549683965814783015124451273604984657747127257615372449205973192657306017239103491074738324033259120) A[12,10] = BigInt(40279545832706233433100438588458933210937500)//BigInt(8896460842799482846916972126377338947215101) A[12,11] = -BigInt(6122933601070769591613093993993358877250)//BigInt(1050517001510235513198246721302027675953) A[13,1] = -BigInt(618675905535482500672800859344538410358660153899637)//BigInt(203544282118214047100119475340667684874292102389760) A[13,4] = -BigInt(4411194916804718600478400319122931000)//BigInt(40373053902469967450761491269633019) A[13,5] = -BigInt(16734711409449292534539422531728520225)//BigInt(1801243715290088669307203927210237952) A[13,6] = BigInt(135137519757054679098042184152749677761254751865630525)//BigInt(16029587794486289597771326361911895112703716593983488) A[13,7] = BigInt(38937568367409876012548551903492196137929710431584875)//BigInt(340956454090191606099548798001469306974758443147264) A[13,8] = -BigInt(6748865855011993037732355335815350667265625)//BigInt(7002880395717424621213565406715087764770357) A[13,9] = -BigInt(1756005520307450928195422767042525091954178296002788308926563193523662404739779789732685671)//BigInt(348767814578469983605688098046186480904607278021030540735333862087061574934154942830062320) A[13,10] = BigInt(53381024589235611084013897674181629296875)//BigInt(8959357584795694524874969598508592944141) α[1] = 44901867737754616851973//1014046409980231013380680 α[6] = 791638675191615279648100000//2235604725089973126411512319 α[7] = 3847749490868980348119500000//15517045062138271618141237517 α[8] = -13734512432397741476562500000//875132892924995907746928783 α[9] = BigInt(12274765470313196878428812037740635050319234276006986398294443554969616342274215316330684448207141)//BigInt(489345147493715517650385834143510934888829280686609654482896526796523353052166757299452852166040) α[10] = -9798363684577739445312500000//308722986341456031822630699 α[11] = 282035543183190840068750//12295407629873040425991 α[12] = -306814272936976936753//1299331183183744997286 αEEst[1] = 10835401739407019406577//244521829356935137978320 αEEst[6] = 13908189778321895491375000//39221135527894265375640567 αEEst[7] = 73487947527027243487625000//296504045773342769773399443 αEEst[8] = 68293140641257649609375000//15353208647806945749946119 αEEst[9] = BigInt(22060647948996678611017711379974578860522018208949721559448560203338437626022142776381)//BigInt(1111542009262325874512959185795727215759010577565736079641376621381577236680929558640) αEEst[10] = -547971229495642458203125000//23237214025700991642563601 αEEst[13] = -28735456870978964189//79783493704265043693 A = map(T,A) α = map(T,α) c = map(T,c) αEEst = map(T,αEEst) return(ExplicitRKTableau(A,c,α,8,adaptiveorder=7,αEEst=αEEst)) end """ Papakostas's Order 6 On Phase-Fitted modified Runge-Kutta Pairs of order 6(5), by Ch. Tsitouras and I. Th. Famelis, International Conference of Numerical Analysis and Applied Mathematics, Crete, (2006) """ function constructPapakostas6(T::Type = Float64) A = zeros(T,9,9) c = zeros(T,9) α = zeros(T,9) αEEst = zeros(T,9) c[2]=17//183 c[3]=12//83 c[4]=18//83 c[5]=71//125 c[6]=42//59 c[7]=199//200 c[8]=1 c[9]=1 A[2,1]=17//183 A[3,1]=3756//117113 A[3,2]=13176//117113 A[4,1]=9//166 A[4,2]=0 A[4,3]=27//166 A[5,1]=207751751//316406250 A[5,2]=0 A[5,3]=-526769377//210937500 A[5,4]=1524242129//632812500 A[6,1]=-4970082682619223281//2887511529739311186 A[6,2]=0 A[6,3]=97919278033879057//13556392158400522 A[6,4]=-407131674007930877068//74078904949579652469 A[6,5]=1237601855204268750000//1753200750473385108433 A[7,1]=176597685527535385020980411//42773485015591331328000000 A[7,2]=0 A[7,3]=-6793162515552646891859//401628967282547712000 A[7,4]=12704926019361287204873446554247//886659402653054716778496000000 A[7,5]=-50728836334509259632278125//32657591718008685915971584 A[7,6]=51536223982796190703//51293749413888000000 A[8,1]=299033520572337573523//66918720793812357519 A[8,2]=0 A[8,3]=-16550269823961899//902146153892364 A[8,4]=49920346343238033627496282//3215735869387500624775563 A[8,5]=-1686432488955761721093750//978844996793357447730403 A[8,6]=161901609084039//149698803705724 A[8,7]=-305146137600000//54760341991955873 A[9,1]=24503//381483 A[9,2]=0 A[9,3]=0 A[9,4]=1366847103121//4106349847584 A[9,5]=20339599609375//75933913767768 A[9,6]=35031290651//194765546144 A[9,7]=16620160000000//11001207123543 A[9,8]=-14933//11016 α[1]=24503//381483 α[2]=0 α[3]=0 α[4]=1366847103121//4106349847584 α[5]=20339599609375//75933913767768 α[6]=35031290651//194765546144 α[7]=16620160000000//11001207123543 α[8]=-14933//11016 αEEst[1]=61010485298317//979331468960880 αEEst[2]=0 αEEst[3]=0 αEEst[4]=320207313882553286621//941222813406992395200 αEEst[5]=6845867841119140625//29008216787127405534 αEEst[6]=124109197949158875473//562495660250110816320 αEEst[7]=19339714537078400000//16810691577722216811 αEEst[8]=-211029377951//210416202900 αEEst[9]=-1//150 A = map(T,A) α = map(T,α) c = map(T,c) αEEst = map(T,αEEst) return(ExplicitRKTableau(A,c,α,6,adaptiveorder=5,αEEst=αEEst,fsal=true)) end """ Lawson's Order 6 An Order 6 Runge-Kutta Process with an Extended Region of Stability, by J. D. Lawson, Siam Journal on Numerical Analysis, Vol. 4, No. 4 (Dec. 1967) pages 620-625. """ function constructLawson6(T::Type = Float64) A = zeros(T,7,7) c = zeros(T,7) α = zeros(T,7) c[2]=26//105-2//315*51^(1//2) c[3]=13//35-1//105*51^(1//2) c[4]=7//8 c[5]=1//2 c[6]=1//8 c[7]=1 A[2,1]=26//105-2//315*51^(1//2) A[3,1]=13//140-1//420*51^(1//2) A[3,2]=39//140-1//140*51^(1//2) A[4,1]=917//1024+133//2048*51^(1//2) A[4,2]=-3339//1024-567//2048*51^(1//2) A[4,3]=1659//512+217//1024*51^(1//2) A[5,1]=-653//5684-1313//45472*51^(1//2) A[5,2]=477//464+81//928*51^(1//2) A[5,3]=-2478015//5433904-21365//339619*51^(1//2) A[5,4]=14568//339619+1528//339619*51^(1//2) A[6,1]=9515//50176+1443//100352*51^(1//2) A[6,2]=-477//1024-81//2048*51^(1//2) A[6,3]=25166643//29980160-2073611//59960320*51^(1//2) A[6,4]=111//46844-275//46844*51^(1//2) A[6,5]=-141//320+21//320*51^(1//2) A[7,1]=-790//49-2243//3528*51^(1//2) A[7,2]=159//4+27//8*51^(1//2) A[7,3]=-44723003//702660-258368//175665*51^(1//2) A[7,4]=327392//316197-2528//105399*51^(1//2) A[7,5]=3158//135-56//45*51^(1//2) A[7,6]=448//27 α[1]=1//70 α[2]=0 α[3]=0 α[4]=256//945 α[5]=58//135 α[6]=256//945 α[7]=1//70 A = map(T,A) α = map(T,α) c = map(T,c) return(ExplicitRKTableau(A,c,α,6)) end """ Tsitouras-Papakostas's Order 6 Cheap Error Estimation for Runge-Kutta methods, by Ch. Tsitouras and S.N. Papakostas, Siam Journal on Scientific Computing, Vol. 20, Issue 6, Nov 1999. """ function constructTsitourasPapakostas6(T::Type = Float64) A = zeros(T,8,8) c = zeros(T,8) α = zeros(T,8) αEEst = zeros(T,8) c[2]=4//27 c[3]=2//9 c[4]=3//7 c[5]=11//16 c[6]=10//13 c[7]=1 c[8]=1 A[2,1]=4//27 A[3,1]=1//18 A[3,2]=1//6 A[4,1]=66//343 A[4,2]=-729//1372 A[4,3]=1053//1372 A[5,1]=13339//49152 A[5,2]=-4617//16384 A[5,3]=5427//53248 A[5,4]=95207//159744 A[6,1]=-6935//57122 A[6,2]=23085//48334 A[6,3]=33363360//273642941 A[6,4]=972160//118442467 A[6,5]=172687360//610434253 A[7,1]=611//1891 A[7,2]=-4617//7564 A[7,3]=6041007//13176488 A[7,4]=12708836//22100117 A[7,5]=-35840000//62461621 A[7,6]=6597591//7972456 A[8,1]=-11107621//14976720 A[8,2]=193023//332816 A[8,3]=1515288591//724706840 A[8,4]=-801736138//352888965 A[8,5]=814013312//606245145 A[8,6]=0 A[8,7]=0 α[1]=131//1800 α[2]=0 α[3]=1121931//3902080 α[4]=319333//1682928 α[5]=262144//2477325 α[6]=4084223//15177600 α[7]=1891//25200 α[8]=0 αEEst[1]=4093//60720 αEEst[2]=0 αEEst[3]=16273467//51284480 αEEst[4]=693889//5376020 αEEst[5]=146456576//626763225 αEEst[6]=15737111//93089280 αEEst[7]=1//25 αEEst[8]=1//23 A = map(T,A) α = map(T,α) c = map(T,c) αEEst = map(T,αEEst) return(ExplicitRKTableau(A,c,α,6)) end """ DormandLockyerMcCorriganPrince Order 6 Global Error Estimation Global Error estimation with Runge-Kutta triples, by J.R.Dormand, M.A.Lockyer, N.E.McCorrigan and P.J.Prince, Computers and Mathematics with Applications, 18 (1989) pages 835-846. """ function constructDormandLockyerMcCorriganPrince6(T::Type = Float64) A = zeros(T,9,9) c = zeros(T,9) α = zeros(T,9) αEEst = zeros(T,9) c[2]=4//39 c[3]=2//13 c[4]=3//13 c[5]=13021//22659 c[6]=39//67 c[7]=86//87 c[8]=1 c[9]=1 A[2,1]=4//39 A[3,1]=1//26 A[3,2]=3//26 A[4,1]=3//52 A[4,2]=0 A[4,3]=9//52 A[5,1]=1406640413621//2478209482476 A[5,2]=0 A[5,3]=-1755653396555//826069827492 A[5,4]=1321105868272//619552370619 A[6,1]=1463083752990765495750771//2783511518115684243554356 A[6,2]=0 A[6,3]=-417501847634533557363//213770948323146013636 A[6,4]=414679390177938970209021//208212903666744217281464 A[6,5]=48490458547529962724706855//2711140218644676453221942744 A[7,1]=-1146771707244809451668952985850//1178428995610817474751161698881 A[7,2]=0 A[7,3]=883524649813655720289029//257840992692019416968811 A[7,4]=-550972740958654450507278066587//325473716439106878117398000544 A[7,5]=-968282586950392419883943203143069455//32828460835559176127341032228568032 A[7,6]=3230428272165469542719//108684219530393291772 A[8,1]=-368234904360842614256649493//316816365518493517722015828 A[8,2]=0 A[8,3]=19247613365107236707//4836252322132133628 A[8,4]=-3627331815384766429266963559//1852940251413900055822878936 A[8,5]=-BigInt(1113629335962635330712822690622675431)//BigInt(31397585101055611361934971740947112) A[8,6]=7515696221383336//210977455127283 A[8,7]=-2466239729887929744//169565652664150899277 A[9,1]=265211783//3930519060 A[9,2]=0 A[9,3]=0 A[9,4]=53198489747//147124055808 A[9,5]=-165257255035734106911939//60068315060067285425920 A[9,6]=1687862952891371//536423949472320 A[9,7]=397538864947251//474823057340620 A[9,8]=-7142267//10794560 α[1]=265211783//3930519060 α[2]=0 α[3]=0 α[4]=53198489747//147124055808 α[5]=-165257255035734106911939//60068315060067285425920 α[6]=1687862952891371//536423949472320 α[7]=397538864947251//474823057340620 α[8]=-7142267//10794560 αEEst[1]=61909109615135759874805//909927606074377568473224 αEEst[2]=0 αEEst[3]=0 αEEst[4]=191026171670373376608531013//532182573462887006324068800 αEEst[5]=-16134045022121298883366692111903//6240618480477179290700503465280 αEEst[6]=123307988643675607258315033//41372548832224521099482580 αEEst[7]=1648598887654728061640361417//1901782716664438059984219160 αEEst[8]=-17//25 αEEst[9]=-259237562821839//28937895739220050 A = map(T,A) α = map(T,α) c = map(T,c) αEEst = map(T,αEEst) return(ExplicitRKTableau(A,c,α,6,fsal=true)) end """ constructTanakaKasugaYamashitaYazaki Order 6 D On the Optimization of Some Eight-stage Sixth-order Explicit Runge-Kutta Method, by M. Tanaka, K. Kasuga, S. Yamashita and H. Yazaki, Journal of the Information Processing Society of Japan, Vol. 34, No. 1 (1993), pages 62 to 74. """ function constructTanakaKasugaYamashitaYazaki6D(T::Type = Float64) A = zeros(T,9,9) c = zeros(T,9) α = zeros(T,9) αEEst = zeros(T,9) c[2]=1//250 c[3]=61//500 c[4]=333//1000 c[5]=56//125 c[6]=221//250 c[7]=31//40 c[8]=1 c[9]=1 A[2,1]=1//250 A[3,1]=-3477//2000 A[3,2]=3721//2000 A[4,1]=5893//960 A[4,2]=-155333//24000 A[4,3]=2//3 A[5,1]=24108461262911872104361127416//10623418096530694059554537475 A[5,2]=-5917919395479054157557671728//2612315925376400178578984625 A[5,3]=232280630406488701857197728//1274810171583683287146544497 A[5,4]=5471535080448445690971520//20898527403011201428631877 A[6,1]=BigInt(6437510603218083862126582571557692575111484748570978316432703829)//BigInt(462655382163643505993146289902892129798145865902401535270912000) A[6,2]=-25083125032238929237211514053//1741543950250933452385989750 A[6,3]=BigInt(191471592827838703414373042295007003190385413366181230722958611)//BigInt(201998644533947923598793335502244867385083328952030670310246400) A[6,4]=-BigInt(723497411149192638041504597211363403675892243966786723532767)//BigInt(584075332429951001052530262062711448991365353250616692288000) A[6,5]=BigInt(20368395873479894700136894919415569811519)//BigInt(12245317457885350045378058014027636736000) A[7,1]=-BigInt(1053235524590544505552847720336831796873415344431325525436387723231012409077771)//BigInt(151882847900737440255262782278159042659114250070600779935950570906179610869760) A[7,2]=110556263584594415870053017499108315//15711834237341169934278611772340752 A[7,3]=BigInt(27883202411953328597139497795565290131372048742640205053329023156822546337651)//BigInt(66313136270946971682878125476803367732416846682610519097035561761715919388672) A[7,4]=-BigInt(134178686375307065766767719823759096451245561636180411142603830822435042239)//BigInt(191743203036279335802347994313396039656741862830639889775698722784806666240) A[7,5]=BigInt(1192365386732594028658711096219658222748871123673792181)//BigInt(1339984907884349782764828286663427050229770700990709760) A[7,6]=1//16 A[8,1]=BigInt(139472860369418680318405579593202554799958275392817771935632674919979204348213061279)//BigInt(37743426679794053215512832350173350810071958850079214921681182582585691688393293824) A[8,2]=-764277739442303538051509575419150250//241488900344982122797870353007980909 A[8,3]=-BigInt(7371553444464925806572282862981502319484039203525663191183225508630714910425342675)//BigInt(4846779318912629207559289658412491740053883371241789941360845137312080681624033792) A[8,4]=BigInt(8801132642676112554927251590711440081546416618547191842516081677278335900180085)//BigInt(2071689842640920660980285733506294870283516476694921779219912685665930091210112) A[8,5]=-BigInt(298368164571373898500032163301238477118159212135073380804555)//BigInt(92828694411466377636989261345454987067336292167034303922176) A[8,6]=-240786613447142598518650//1482036690339563123748203 A[8,7]=52860332724128242560//47245264125077723987 A[9,1]=-BigInt(70335440697678472884389500304621917516220835017456192585613672496070708987720514364185979355621)//BigInt(5219173052726261361794285880851433852375067076864459834761850509124293962243718819226202095616) A[9,2]=BigInt(55057602528492150957512473843952756334951250)//BigInt(3825542672812430870769605520946554801303339) A[9,3]=-BigInt(77980782561050861485751338404689118536953083418300092316918550174770045337324006382574405615)//BigInt(39424363920896283457695961995433007195554371610645094834381872512385376722954438065398121984) A[9,4]=BigInt(2177533338045948502564188609357006384381697842038103984795932632439956259386338338574157919)//BigInt(387582376874948376545286159615164080595940845312336117877987331598679530792396290947668352) A[9,5]=-BigInt(1429465993644112145770349737782620061868013755109155817982734692969685809)//BigInt(295236772314145005435682097407317274403116518334589230312060489793912832) A[9,6]=-BigInt(162882286796273095967517045376244881875035)//BigInt(1226056358355165076488310595889719561316528) A[9,7]=61//43 A[9,8]=0 α[1]=4783097999//163657290888 α[2]=0 α[3]=1043945712500000//4811958427553991 α[4]=12731480000000//212031396478881 α[5]=3420427578125//10842916254408 α[6]=22394286718750//335690566308279 α[7]=601452300800//2318247773037 α[8]=11555307017//221477311314 α[9]=0 αEEst[1]=401359714631498171030059910//15267054426382006481740339287 αEEst[2]=0 αEEst[3]=1872374247198383241346435120000//8180262298260999309783933205623 αEEst[4]=215533585653762668869792000//12486209626895034583833718599 αEEst[5]=10643952265127448732528470125//29318848171314267144678228486 αEEst[6]=49992752929730967217476426250//454671173972183983205823101943 αEEst[7]=41077106106719584492498868480//192232681858829974979332840131 αEEst[8]=0 αEEst[9]=13//318 A = map(T,A) α = map(T,α) c = map(T,c) αEEst = map(T,αEEst) return(ExplicitRKTableau(A,c,α,6)) end """ constructTanakaKasugaYamashitaYazaki Order 6 C On the Optimization of Some Eight-stage Sixth-order Explicit Runge-Kutta Method, by M. Tanaka, K. Kasuga, S. Yamashita and H. Yazaki, Journal of the Information Processing Society of Japan, Vol. 34, No. 1 (1993), pages 62 to 74. """ function constructTanakaKasugaYamashitaYazaki6C(T::Type = Float64) A = zeros(T,9,9) c = zeros(T,9) α = zeros(T,9) αEEst = zeros(T,9) c[2]=1//100 c[3]=11//100 c[4]=33//100 c[5]=43//100 c[6]=177//200 c[7]=39//50 c[8]=1 c[9]=1 A[2,1]=1//100 A[3,1]=-99//200 A[3,2]=121//200 A[4,1]=931//600 A[4,2]=-1133//600 A[4,3]=2//3 A[5,1]=91958765679448742111//109106230010372040600 A[5,2]=-2894503641421476941//3306249394253698200 A[5,3]=42536885251557644//181843716683953401 A[5,4]=124767011204926940//545531150051860203 A[6,1]=BigInt(2432582077139652291040741961658278931048241123)//BigInt(1042517831971586211298945274530293397770600000) A[6,2]=-5086533956627126527//2404545014002689600 A[6,3]=-BigInt(28233375209799471901757749652001556425049843)//BigInt(193956805948202085823059585959124353073600000) A[6,4]=-BigInt(64143717120482360099075816554051641763250417)//BigInt(80815335811750869092941494149635147114000000) A[6,5]=368355687225024998241998861697//229321716687851233749224000000 A[7,1]=-BigInt(218146655544801488428299362829574333746140908535614198091)//BigInt(129671440168534479253493745415914920591854154832474480000) A[7,2]=3225815231796122073236882//1899026582543293621912725 A[7,3]=BigInt(90474621818961969641326264897298400503244755498598833923)//BigInt(192999352808981550516827900153919881811131765332055040000) A[7,4]=-BigInt(1469845402576002022391874744301978646600063888157276277)//BigInt(1827645386448688925348749054487877668666020505038400000) A[7,5]=BigInt(236537380548621166783390004339783981292483)//BigInt(228189687334919622052120277661453593600000) A[7,6]=1//16 A[8,1]=BigInt(20745045233020258374773011834908499926526821480199726607317)//BigInt(6703021612773885250214340203318686755448684469761576087200) A[8,2]=-10062163770239468036471200//3356073465031554356885103 A[8,3]=-BigInt(6600722659737487671039676673294259975246135762925896417)//BigInt(9474444736357442876691573864537554989039417965607475200) A[8,4]=BigInt(80225477702329113082476652949959330669319741271791060069)//BigInt(21801989308095252074205042131464260060005478841312656000) A[8,5]=-BigInt(572101553192975384084889663939632620168841743)//BigInt(192662470105149826757328803363450656933632000) A[8,6]=-223351084874296955//945158103640587564 A[8,7]=3931692514354//3491636633667 A[9,1]=BigInt(1042961260483134998083167734826788831236189026066652724982180836527)//BigInt(166222076505494822625791597828801737913859767695104851064439232400) A[9,2]=-4754359005740226898713768141325//838855890397374389475867956178 A[9,3]=-BigInt(42750676060439966787211838095564987064848978722439603628501619959)//BigInt(27488922212795009426487499382540855882395413779035427566212172800) A[9,4]=BigInt(205686747816311239126907254218859168711623965851131881271189228377)//BigInt(77312593723485964011996092013396157169237101253537140029971736000) A[9,5]=-BigInt(75946282944430327533335124997772458772050053465158929)//BigInt(62109575359095533697158937720595685916206114892672000) A[9,6]=-264375495891732597330199815785//1220640854393687869579464225198 A[9,7]=8//11 A[9,8]=0 α[1]=92725517//4525454934 α[2]=0 α[3]=203217125//939422946 α[4]=16794200//510230259 α[5]=7552075//21858018 α[6]=22967296000//562167528741 α[7]=229749950//798566769 α[8]=67671311//1203893922 α[9]=0 αEEst[1]=1773713828809813620047//76546041056100425133324 αEEst[2]=0 αEEst[3]=576936961199857272475//2856614455209710336064 αEEst[4]=28907772101756401885//257621374646831423844 αEEst[5]=26154985271235615605//103780652247492444224 αEEst[6]=-12285231659070809284400//206713297399385556110331 αEEst[7]=237440851439576891930//613972674142959478347 αEEst[8]=0 αEEst[9]=1//12 A = map(T,A) α = map(T,α) c = map(T,c) αEEst = map(T,αEEst) return(ExplicitRKTableau(A,c,α,6)) end """ constructTanakaKasugaYamashitaYazaki Order 6 B On the Optimization of Some Eight-stage Sixth-order Explicit Runge-Kutta Method, by M. Tanaka, K. Kasuga, S. Yamashita and H. Yazaki, Journal of the Information Processing Society of Japan, Vol. 34, No. 1 (1993), pages 62 to 74. """ function constructTanakaKasugaYamashitaYazaki6B(T::Type = Float64) A = zeros(T,9,9) c = zeros(T,9) α = zeros(T,9) αEEst = zeros(T,9) c[2]=1//5 c[3]=3//20 c[4]=2//5 c[5]=1//2 c[6]=3//4 c[7]=4//5 c[8]=1 c[9]=1 A[2,1]=1//5 A[3,1]=3//32 A[3,2]=9//160 A[4,1]=-71//400 A[4,2]=-53//400 A[4,3]=71//100 A[5,1]=20861387//246441600 A[5,2]=-2674671//27382400 A[5,3]=14952973//61610400 A[5,4]=555163//2053680 A[6,1]=197611764884981873//575894510914896000 A[6,2]=16376643//136912000 A[6,3]=-281745778473052331//1007815394101068000 A[6,4]=-761016212927233//4799120924290800 A[6,5]=7119424000//9814726677 A[7,1]=2141650356947829667//15684127070697870750 A[7,2]=124460793//710231000 A[7,3]=729393184893351907//54894444747442547625 A[7,4]=-857270676309077671//4182433885519432200 A[7,5]=496013770812769//855353429900550 A[7,6]=1//10 A[8,1]=59331555925548335981//3757538915366421931200 A[8,2]=-807212549//1636098400 A[8,3]=3149950445245624075843//6575693101891238379600 A[8,4]=213748931907436689007//250502594357761462080 A[8,5]=-229203323900597//10246103543907504 A[8,6]=-1306535//1043952 A[8,7]=35275//24856 A[9,1]=-82982860685611836065717//78861151462250930572800 A[9,2]=-4258977307//34337529600 A[9,3]=257897426390213500944749//138007015058939128502400 A[9,4]=7684590445667080762721//5257410097483395371520 A[9,5]=-418846778208596923//215039562243682176 A[9,6]=14813735//606173568 A[9,7]=10//13 A[9,8]=0 α[1]=253//6048 α[2]=0 α[3]=71200//292383 α[4]=725//7056 α[5]=124//441 α[6]=-160//1323 α[7]=10375//26208 α[8]=239//4284 α[9]=0 αEEst[1]=1061069//26544672 αEEst[2]=0 αEEst[3]=17566690//75486411 αEEst[4]=3286405//10322928 αEEst[5]=-148157//1935549 αEEst[6]=3649510//5806647 αEEst[7]=-27404525//115026912 αEEst[8]=0 αEEst[9]=2//21 A = map(T,A) α = map(T,α) c = map(T,c) αEEst = map(T,αEEst) return(ExplicitRKTableau(A,c,α,6)) end """ TanakaKasugaYamashitaYazaki Order 6 A On the Optimization of Some Eight-stage Sixth-order Explicit Runge-Kutta Method, by M. Tanaka, K. Kasuga, S. Yamashita and H. Yazaki, Journal of the Information Processing Society of Japan, Vol. 34, No. 1 (1993), pages 62 to 74. """ function constructTanakaKasugaYamashitaYazaki6A(T::Type = Float64) A = zeros(T,9,9) c = zeros(T,9) α = zeros(T,9) αEEst = zeros(T,9) c[2]=1//100 c[3]=63//500 c[4]=63//200 c[5]=91//200 c[6]=21//25 c[7]=121//200 c[8]=1 c[9]=1 A[2,1]=1//100 A[3,1]=-3339//5000 A[3,2]=3969//5000 A[4,1]=7409//2400 A[4,2]=-2751//800 A[4,3]=2//3 A[5,1]=879267900536971128212809//854114648586066093506400 A[5,2]=-32049516354411396813791//31633875873558003463200 A[5,3]=415928927703664756394//3202929932197747850649 A[5,4]=4945815476584113561248//16014649660988739253245 A[6,1]=BigInt(264364422150553063350374034500534568610052754809099)//BigInt(76594645373513572971811613560122656201571540854400) A[6,2]=-4942738251715985006237//988558621048687608225 A[6,3]=BigInt(857070668647981630049588712158864150151776268553)//BigInt(235675831918503301451728041723454326774066279552) A[6,4]=-BigInt(575165070977427883238172569267537197276341339157)//BigInt(157117221279002200967818694482302884516044186368) A[6,5]=41539505110086641686706761875//17218030314852652247321445632 A[7,1]=-BigInt(7384332843467384251748880311988170834105043085902030057004642579)//BigInt(1315398593371357506652159044711447383311024734178303038737049600) A[7,2]=5656750329476775131209405217353//831860122007390417396028410400 A[7,3]=-BigInt(333346623364947547771242996577246247659950898456372585977245839)//BigInt(190226873502934777885081461850578544663440500019631516371204096) A[7,4]=BigInt(7100883704616862706920092769411499000150383998766128086388523)//BigInt(4497089208107205150947552289611785925849657210865993294827520) A[7,5]=-BigInt(2179578925723796720158217363641743790096295)//BigInt(4632538471751007430155293544377007681933312) A[7,6]=1//16 A[8,1]=-BigInt(2032945134330261992294105795025746391376544519999264654886089003925834569)//BigInt(370204608795106233716533306294462150635128163745469969879178202980349056) A[8,2]=104400636802444703223540033778325//10372562449415222120508569935887 A[8,3]=-BigInt(1282736937871158931724790946085686543914309745805108288599439251128389375)//BigInt(121675640652937013878860597173704343215741424447811808281967661119415424) A[8,4]=BigInt(69591790067148620077849942355176681820069270105493261193720164528698455)//BigInt(5177686836295192079951514773349120987903890402034545033275219622102784) A[8,5]=-BigInt(277247257204675144642737406657510117972977064771225085)//BigInt(32594437558796549389422402191407749206480247247281408) A[8,6]=206057847729988755370340//362497215363232474360143 A[8,7]=482015452915471040//328927231085842477 A[9,1]=BigInt(149325572050522001715919965533803802996087500406626797668220209959218571)//BigInt(32753013690352528622466162111781793695504165546622151035212132158026496) A[9,2]=-BigInt(2802049259729837394660870725638222600)//BigInt(886389794903952023701108666481909127) A[9,3]=-BigInt(247046029472637555138263024224503808349122551711211979490561742618517516481875)//BigInt(62386925133959095591686249139496495395591773904073928005189977758511970759424) A[9,4]=BigInt(17088511387117619026306232328319322341432000611501090216270103426420066685195)//BigInt(2654762771657833854965372303808361506195394634215911830008084159936679606784) A[9,5]=-BigInt(1905127968294488657830239165534292922200553202371049568415)//BigInt(428517757803501533934438346448102436732262910333700649472) A[9,6]=BigInt(348519474109306281943020669241827485)//BigInt(1086125473870224085200713153977775988) A[9,7]=5//4 A[9,8]=0 α[1]=703450322//19272872619 α[2]=0 α[3]=1197916718750000//6147559103622993 α[4]=20671732000//146035199457 α[5]=5917108000//36937869969 α[6]=23017645375//95720436504 α[7]=23374664000//131057876103 α[8]=4163735081//86609369112 α[9]=0 αEEst[1]=572720177479183548778421//9495052764959466387724734 αEEst[2]=0 αEEst[3]=4729369086248342935597812500//45049042428818825645161957701 αEEst[4]=1262015363899500143599090//3413510704758155271620049 αEEst[5]=-106945493251106516170//843200958241401898187 αEEst[6]=711860694955815806929355//3483263158495231241410497 αEEst[7]=1763218383368811180347180//5312516873137942812595713 αEEst[8]=0 αEEst[9]=1//18 A = map(T,A) α = map(T,α) c = map(T,c) αEEst = map(T,αEEst) return(ExplicitRKTableau(A,c,α,6)) end """ Mikkawy-Eisa Order 6 A general four-parameter non-FSAL embedded Runge–Kutta algorithm of orders 6 and 4 in seven stages, by M.E.A. El-Mikkawy and M.M.M. Eisa, Applied Mathematics and Computation, Vol. 143, No. 2, (2003) pages 259 to 267. """ function constructMikkawyEisa(T::Type = Float64) A = zeros(T,7,7) c = zeros(T,7) α = zeros(T,7) αEEst = zeros(T,7) c[2]=3//19 c[3]=9//38 c[4]=342//683 c[5]=22//27 c[6]=9//11 c[7]=1 A[2,1]=3//19 A[3,1]=9//152 A[3,2]=27//152 A[4,1]=94474764//318611987 A[4,2]=-310753854//318611987 A[4,3]=375818328//318611987 A[5,1]=-1073687692//12631821129 A[5,2]=1421618//911979 A[5,3]=-797223265256//505937677851 A[5,4]=8812260835312//9612815879169 A[6,1]=-140767614//1662675883 A[6,2]=302157//185009 A[6,3]=-6054413848056//3590051357651 A[6,4]=1284501562654185//1332472080292672 A[6,5]=-34222143195//4324664670656 A[7,1]=666881867//2248862211 A[7,2]=-19133//162361 A[7,3]=337642670416//1442156395089699 A[7,4]=388258038731717539//1050707717099750592 A[7,5]=11900400291315//1115308365632 A[7,6]=-1730758084946//169374832839 α[1]=456023//6094440 α[2]=0 α[3]=85078761640//257589787311 α[4]=5499272526534791//27584466984914400 α[5]=4653828837//944530400 α[6]=-962601827//208639800 α[7]=162361//1977800 αEEst[1]=851969693//10815599520 αEEst[2]=0 αEEst[3]=6675679897963//21163704130305 αEEst[4]=14982723272271221//67174249160747520 αEEst[5]=86506623776547//20114719398400 αEEst[6]=-4 αEEst[7]=162361//1977800 A = map(T,A) α = map(T,α) c = map(T,c) αEEst = map(T,αEEst) return(ExplicitRKTableau(A,c,α,6)) end """ Chummund's First Order 6 method A three-dimensional family of seven-step Runge-Kutta methods of order 6, by G. M. Chammud (Hammud), Numerical Methods and programming, 2001, Vol.2, 2001, pages 159-166 (Advanced Computing Scientific journal published by the Research Computing Center of the Lomonosov Moscow State Univeristy) """ function constructChummund6(T::Type = Float64) A = zeros(T,7,7) c = zeros(T,7) α = zeros(T,7) c[2]=4//7 c[3]=5//7 c[4]=6//7 c[5]=1//2-1//10*5^(1//2) c[6]=1//2+1//10*5^(1//2) c[7]=1 A[2,1]=4//7 A[3,1]=115//112 A[3,2]=-5//16 A[4,1]=589//630 A[4,2]=5//18 A[4,3]=-16//45 A[5,1]=-29//6000*5^(1//2)+229//1200 A[5,2]=-187//1200*5^(1//2)+119//240 A[5,3]=34//375*5^(1//2)-14//75 A[5,4]=-3//100*5^(1//2) A[6,1]=71//2400-587//12000*5^(1//2) A[6,2]=187//480-391//2400*5^(1//2) A[6,3]=-38//75+26//375*5^(1//2) A[6,4]=27//80-3//400*5^(1//2) A[6,5]=1//4*5^(1//2)+1//4 A[7,1]=43//160*5^(1//2)-49//480 A[7,2]=51//32*5^(1//2)-425//96 A[7,3]=-4//5*5^(1//2)+52//15 A[7,4]=3//16*5^(1//2)-27//16 A[7,5]=-3//4*5^(1//2)+5//4 A[7,6]=5//2-1//2*5^(1//2) α[1]=1//12 α[2]=0 α[3]=0 α[4]=0 α[5]=5//12 α[6]=5//12 α[7]=1//12 A = map(T,A) α = map(T,α) c = map(T,c) return(ExplicitRKTableau(A,c,α,6)) end """ Chummund's Second Order 6 method A three-dimensional family of seven-step Runge-Kutta methods of order 6, by G. M. Chammud (Hammud), Numerical Methods and programming, 2001, Vol.2, 2001, pages 159-166 (Advanced Computing Scientific journal published by the Research Computing Center of the Lomonosov Moscow State Univeristy) """ function constructChummund62(T::Type = Float64) A = zeros(T,7,7) c = zeros(T,7) α = zeros(T,7) c[2]=750557//18870600 c[3]=748997//1685240 c[4]=624713//833636 c[5]=1//2-1//10*5^(1//2) c[6]=1//2+1//10*5^(1//2) c[7]=1 A[2,1]=750557//18870600 A[3,1]=-BigInt(2618524936181374161531835574563010037)//BigInt(1243120420984397996179713098690223560)+3789435780183611636743//41406596067467160927808051177*826321815619^(1//2) A[3,2]=79275599164011507825658766043059867//31078010524609949904492827467255589-3789435780183611636743//41406596067467160927808051177*826321815619^(1//2) A[4,1]=-BigInt(268026226098849050460940570252841432600176823648465011650077)//BigInt(5150236239747870879984989921782016384957547807596571966124)+BigInt(75921270833932035355544618340621749422299940322072982)//BigInt(1287559059936967719996247480445504096239386951899142991531)*826321815619^(1//2) A[4,2]=BigInt(31518353641767830515869136922708314461169361966520350023632789875)//BigInt(553066166067338591675176701039997652878686435650554094276321374)-BigInt(35636316203172669679112815645100268396443783685272961209975)//BigInt(553066166067338591675176701039997652878686435650554094276321374)*826321815619^(1//2) A[4,3]=-BigInt(2316654837809440895270062146979974449795559590002488266576771951)//BigInt(551916642155010749548578353874329493980040828668519599378437454)+BigInt(3018332469640261586051209426912089268215645947971953376787)//BigInt(551916642155010749548578353874329493980040828668519599378437454)*826321815619^(1//2) A[5,1]=7557928766176693537//21071504819547814620-3808101371666611447//21071504819547814620*5^(1//2)+1830930386299//21071504819547814620*826321815619^(1//2)+2746507673963//21071504819547814620*826321815619^(1//2)*5^(1//2) A[5,2]=-306676166043910279990585069725//1347799564492996272083268634084+281193211701441857849531780175//1347799564492996272083268634084*5^(1//2)-135831974232490143067425//1347799564492996272083268634084*826321815619^(1//2)-203756277349889111562225//1347799564492996272083268634084*826321815619^(1//2)*5^(1//2) A[5,3]=59599519641651017073742851703//154850028368562600922782726012-19070452237413934132405792717//154850028368562600922782726012*5^(1//2)+3249940179998629473139//154850028368562600922782726012*826321815619^(1//2)+4875109240133213731043//154850028368562600922782726012*826321815619^(1//2)*5^(1//2) A[5,4]=-179555217756675442542605916547//11203233050287498632877256911020-10652462513464648264279700899//2240646610057499726575451382204*5^(1//2)-79525200332352343539019//11203233050287498632877256911020*826321815619^(1//2)-119292669246564200876603//11203233050287498632877256911020*826321815619^(1//2)*5^(1//2) A[6,1]=-6195448976102809771//10535752409773907310-1420753767230413123//10535752409773907310*5^(1//2)+2689284093484//5267876204886953655*826321815619^(1//2)+57223580479//5267876204886953655*826321815619^(1//2)*5^(1//2) A[6,2]=583482891804126952797345387075//673899782246498136041634317042+55992621594804402807018397650//336949891123249068020817158521*5^(1//2)-254249689295131189339575//336949891123249068020817158521*826321815619^(1//2)-49153620024967407924975//673899782246498136041634317042*826321815619^(1//2)*5^(1//2) A[6,3]=19364321623328820251751318839//77425014184281300461391363006+1040656060288496357918934418//38712507092140650230695681503*5^(1//2)-11049041979964718221319//38712507092140650230695681503*826321815619^(1//2)-19382660452189380569035//77425014184281300461391363006*826321815619^(1//2)*5^(1//2) A[6,4]=BigInt(115027957037506741565069794992944458832461)//BigInt(237309179139566721503969545043331027555510)-BigInt(19371908299625437579435144814464589438224)//BigInt(118654589569783360751984772521665513777755)*5^(1//2)+BigInt(64263508455127995016232148201271346)//BigInt(118654589569783360751984772521665513777755)*826321815619^(1//2)-BigInt(41163600684913040546114995026601907)//BigInt(237309179139566721503969545043331027555510)*826321815619^(1//2)*5^(1//2) A[6,5]=43435632261//211822050005*5^(1//2)-43432990301//84728820002+205839//423644100010*826321815619^(1//2)*5^(1//2)-514//42364410001*826321815619^(1//2) A[7,1]=9047270149938488929//4214300963909562924+2216536302042479231//1404766987969854308*5^(1//2)-991800665293//1404766987969854308*826321815619^(1//2)*5^(1//2)-12588066760235//4214300963909562924*826321815619^(1//2) A[7,2]=-4301448087821718128020528522125//1347799564492996272083268634084-2525818490403297345388026853875//1347799564492996272083268634084*5^(1//2)+1510317586999119637060875//1347799564492996272083268634084*826321815619^(1//2)*5^(1//2)+5764153657065074502128625//1347799564492996272083268634084*826321815619^(1//2) A[7,3]=-491640814441543287886227446905//154850028368562600922782726012+752920605871714580844952275//1564141700692551524472552788*5^(1//2)+1711626851729573101365//1564141700692551524472552788*826321815619^(1//2)*5^(1//2)+204731138699301217060685//154850028368562600922782726012*826321815619^(1//2) A[7,4]=-BigInt(222449163212150849227576515026323419478375)//BigInt(94923671655826688601587818017332411022204)+BigInt(805495551976820677052288591080735143509)//BigInt(958824966220471602036240586033660717396)*5^(1//2)-BigInt(253684995628220543918448388547756365)//BigInt(94923671655826688601587818017332411022204)*826321815619^(1//2)+BigInt(882636009291931420795306294477883)//BigInt(958824966220471602036240586033660717396)*826321815619^(1//2)*5^(1//2) A[7,5]=-44506854521//84728820002*5^(1//2)+214493500755//42364410001-205839//84728820002*826321815619^(1//2)*5^(1//2)+2570//42364410001*826321815619^(1//2) A[7,6]=5//2-1//2*5^(1//2) α[1]=1//12 α[2]=0 α[3]=0 α[4]=0 α[5]=5//12 α[6]=5//12 α[7]=1//12 A = map(T,A) α = map(T,α) c = map(T,c) return(ExplicitRKTableau(A,c,α,6)) end """ Anton Hutas Second Order 6 method Une amélioration de la méthode de Runge-Kutta-Nyström pour la résolution numérique des équations différentielles du premièr ordre, by Anton Huta, Acta Fac. Nat. Univ. Comenian Math., Vol. 1, pages 201-224 (1956). """ function constructHuta62(T::Type = Float64) A = zeros(T,8,8) c = zeros(T,8) α = zeros(T,8) c[2]=1//9 c[3]=1//6