diff --git a/src/runge_kutta.jl b/src/runge_kutta.jl index a080602de..5712148ca 100644 --- a/src/runge_kutta.jl +++ b/src/runge_kutta.jl @@ -20,7 +20,8 @@ immutable TableauRKExplicit{Name, S, T} <: Tableau{Name, S, T} @assert istril(a) @assert S==length(c)==size(a,1)==size(a,2)==size(b,2) @assert size(b,1)==length(order) - @assert norm(sum(a,2)-c'',Inf)<1e-10 # consistency. + aa = T<:Rational ? convert(Matrix{Rational{BigInt}}, a) : a # to avoid overflow + @assert norm(sum(aa,2)-c'',Inf)<1e-10 # consistency. new(order,a,b,c) end end @@ -100,7 +101,7 @@ const bt_rk23 = TableauRKExplicit(:bogacki_shampine,(2,3), Rational{Int64}, 2/9 1/3 4/9 0], [7/24 1/4 1/3 1/8 2/9 1/3 4/9 0], - [0, 1//2, 3//4, 1] + [0, 1//2, 3//4, 1] ) # Fehlberg https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta%E2%80%93Fehlberg_method @@ -151,6 +152,28 @@ const bt_feh78 = TableauRKExplicit(:feh78, (7,8), Rational{Int64}, [0, 2//27, 1//9, 1//6 , 5//12, 1//2 , 5//6 , 1//6 , 2//3 , 1//3 , 1 , 0, 1] ) +# Dormand-Prince Order 7(8) coefficients. +# High order embedded Runge-Kutta formulae. Prince & Dormand 1981 +# http://dx.doi.org/10.1016/0771-050X(81)90010-3 +# Copied from DifferentialEquations.jl +const bt_dopri78 = TableauRKExplicit(:dopri78, (7,8), Rational{Int64}, + [0 0 0 0 0 0 0 0 0 0 0 0 0 + 1//18 0 0 0 0 0 0 0 0 0 0 0 0 + 1//48 1//16 0 0 0 0 0 0 0 0 0 0 0 + 1//32 0 3//32 0 0 0 0 0 0 0 0 0 0 + 5//16 0 -75//64 75//64 0 0 0 0 0 0 0 0 0 + 3//80 0 0 3//16 3//20 0 0 0 0 0 0 0 0 + 29443841//614563906 0 0 77736538//692538347 -28693883//1125000000 23124283//1800000000 0 0 0 0 0 0 0 + 16016141//946692911 0 0 61564180//158732637 22789713//633445777 545815736//2771057229 -180193667//1043307555 0 0 0 0 0 0 + 39632708//573591083 0 0 -433636366//683701615 -421739975//2616292301 100302831//723423059 790204164//839813087 800635310//3783071287 0 0 0 0 0 + 246121993//1340847787 0 0 -37695042795//15268766246 -309121744//1061227803 -12992083//490766935 6005943493//2108947869 393006217//1396673457 123872331//1001029789 0 0 0 0 + -1028468189//846180014 0 0 8478235783//508512852 1311729495//1432422823 -10304129995//1701304382 -48777925059//3047939560 15336726248//1032824649 -45442868181//3398467696 3065993473//597172653 0 0 0 + 185892177//718116043 0 0 -3185094517//667107341 -477755414//1098053517 -703635378//230739211 5731566787//1027545527 5232866602//850066563 -4093664535//808688257 3962137247//1805957418 65686358//487910083 0 0 + 403863854//491063109 0 0 -5068492393//434740067 -411421997//543043805 652783627//914296604 11173962825//925320556 -13158990841//6184727034 3936647629//1978049680 -160528059//685178525 248638103//1413531060 0 0], + [14005451//335480064 0 0 0 0 -59238493//1068277825 181606767//758867731 561292985//797845732 -1041891430//1371343529 760417239//1151165299 118820643//751138087 -528747749//2220607170 1//4 + 13451932//455176623 0 0 0 0 -808719846//976000145 1757004468//5645159321 656045339//265891186 -3867574721//1518517206 465885868//322736535 53011238//667516719 2//45 0], + [0, 1//18, 1//12, 1//8, 5//16, 3//8, 59//400, 93//200, 5490023248//9719169821, 13//20, 1201146811//1299019798, 1, 1] + ) ################################ # Fixed step Runge-Kutta methods @@ -173,7 +196,7 @@ function oderk_fixed{N,S}(fn, y0::AbstractVector, tspan, # TODO: instead of AbstractVector use a Holy-trait # Needed interface: - # On components: + # On components: # On y0 container: length, deepcopy, similar, setindex! # On time container: getindex, convert. length @@ -213,6 +236,7 @@ ode45_dp(fn, y0, tspan; kwargs...) = oderk_adapt(fn, y0, tspan, bt_dopri5; kwarg # Use Dormand-Prince version of ode45 by default const ode45 = ode45_dp ode78(fn, y0, tspan; kwargs...) = oderk_adapt(fn, y0, tspan, bt_feh78; kwargs...) +ode78_dp(fn, y0, tspan; kwargs...) = oderk_adapt(fn, y0, tspan, bt_dopri78; kwargs...) function oderk_adapt(fn, y0, tspan, btab::TableauRKExplicit; kwords...) # For y0 which don't support indexing. @@ -233,7 +257,7 @@ function oderk_adapt{N,S}(fn, y0::AbstractVector, tspan, btab_::TableauRKExplici # - note that the type of the components might change! # On y0 container: length, similar, setindex! # On time container: getindex, convert, length - + # For y0 which support indexing. Currently y0<:AbstractVector but # that could be relaxed with a Holy-trait. !isadaptive(btab_) && error("Can only use this solver with an adaptive RK Butcher table") diff --git a/test/runtests.jl b/test/runtests.jl index 9c88cd281..47bf34e38 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,7 +1,7 @@ using ODE using Base.Test -tol = 1e-2 +tol = 0.013 solvers = [ ## Non-stiff @@ -18,6 +18,7 @@ solvers = [ ODE.ode45_dp, ODE.ode45_fe, ODE.ode78, + ODE.ode78_dp, ## Stiff # fixed-step