Skip to content
This repository was archived by the owner on Sep 9, 2024. It is now read-only.

Added DOPRI78 tableau from DifferentialEquations.jl #101

Closed
mauro3 wants to merge 1 commit into from
Closed

Added DOPRI78 tableau from DifferentialEquations.jl #101

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
32 changes: 28 additions & 4 deletions src/runge_kutta.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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

Expand Down Expand Up @@ -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.
Expand All @@ -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")
Expand Down
3 changes: 2 additions & 1 deletion test/runtests.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
using ODE
using Base.Test

tol = 1e-2
tol = 0.013

solvers = [
## Non-stiff
Expand All @@ -18,6 +18,7 @@ solvers = [
ODE.ode45_dp,
ODE.ode45_fe,
ODE.ode78,
ODE.ode78_dp,

## Stiff
# fixed-step
Expand Down