From 9eb745bcf0478ba4958f1851802bdd51c5e3822a Mon Sep 17 00:00:00 2001 From: Mauro Werder Date: Wed, 15 Apr 2015 23:05:14 +0200 Subject: [PATCH] New implementation of Runge-Kutta solvers Features: - adaptive solvers with orders: 21, 45, 54, 78 - fixed step solvers with order: 1-4 - faster for Vector ODE states - a bit slower for scalar ODE states - MIT licensed - added support for keyword `points` - more general framework for coefficients tableaus - more general framework for determining the types of the computation TODO after merging this PR: - make solving scalar-like states faster again - use more general tableau-framework for the other solvers too Bumped required Julia to 0.3, added Compat 0.4.1 requirement. --- README.md | 10 +- REQUIRE | 3 +- src/ODE.jl | 418 ++++++++++++---------------------- src/runge_kutta.jl | 480 ++++++++++++++++++++++++++++++++++++++++ test/interface-tests.jl | 79 +++++++ test/runtests.jl | 40 ++-- 6 files changed, 729 insertions(+), 301 deletions(-) create mode 100644 src/runge_kutta.jl create mode 100644 test/interface-tests.jl diff --git a/README.md b/README.md index ce9c06875..006048044 100644 --- a/README.md +++ b/README.md @@ -15,7 +15,7 @@ Since then, quite a lot has happened in the package, and the best way to use cur Currently, `ODE` exports the following adaptive solvers: * `ode23`: 2nd order adaptive solver with 3rd order error control, using the Bogacki–Shampine coefficients -* `ode45`: 4th order adaptive solver with 5th order error control, using the Dormand Prince coefficients. Fehlberg and Cash-Karp coefficients are also available. +* `ode54`: 4th order adaptive solver with 5th order error control, using the Dormand Prince coefficients. Fehlberg and Cash-Karp coefficients are also available. * `ode78`: 7th order adaptive solver with 8th order error control, using the Fehlberg coefficients. * `ode23s`: 2nd/3rd order adaptive solver for stiff problems, using a modified Rosenbrock triple. @@ -24,18 +24,20 @@ all of which have the following basic API: tout, yout = odeXX(F, y0, tspan; keywords...) -to solve the explicit ODE defined by dy/dt = F(t,y). A few other solvers are also exported, see the source code for details. +to solve the explicitly defined ODE by dy/dt = F(t,y). A few other solvers are also exported, see the source code for details. The adaptive solvers accept the following keywords - `norm`: user-supplied norm for determining the error `E` (default `Base.vecnorm`), - `abstol` and/or `reltol`: an integration step is accepted if `E <= abstol || E <= reltol*abs(y)` (defaults `reltol = 1e-5`, `abstol = 1e-8`), - `maxstep`, `minstep` and `initstep`: determine the maximal, minimal and initial integration step (defaults `minstep=|tspan[end] - tspan[1]|/1e9`, `maxstep=|tspan[end] - tspan[1]|/2.5` and automatic initial step estimation). - -Additionally, `ode23s` supports - `points=:all` (default): output is given for each value in `tspan` as well as for each intermediate point the solver used. - `points=:specified`: output is given only for each value in `tspan`. + +Additionally, `ode23s` solver supports - `jacobian = G(t,y)`: user-supplied Jacobian G(t,y) = dF(t,y)/dy (default estimate by finite-difference method). +There are also fixed step Runge-Kutta and Rosenbrock solvers available. + # Need something long-term reliable right now? See [the Sundials.jl package](https://github.com/julialang/sundials.jl), which provides wrappers for the excellent Sundials ODE solver library. diff --git a/REQUIRE b/REQUIRE index a9aa1a122..8e8352dde 100644 --- a/REQUIRE +++ b/REQUIRE @@ -1,2 +1,3 @@ -julia 0.2 +julia 0.3 Polynomials +Compat 0.4.1 diff --git a/src/ODE.jl b/src/ODE.jl index c57eb4610..a364988d8 100644 --- a/src/ODE.jl +++ b/src/ODE.jl @@ -3,20 +3,57 @@ module ODE using Polynomials +using Compat ## minimal function export list # adaptive non-stiff: -export ode23, ode45, ode78 +export ode45, ode78 +# non-adaptive non-stiff: +export ode4, ode4ms # adaptive stiff: export ode23s -# non-adaptive: -export ode4s, ode4ms, ode4 +# non-adaptive stiff: +export ode4s -## complete function export list -#export ode23, ode4, -# oderkf, ode45, ode45_dp, ode45_fb, ode45_ck, -# oderosenbrock, ode4s, ode4s_kr, ode4s_s, -# ode4ms, ode_ms +## complete function export list: see runtests.jl + +############################################################################### +## Coefficient Tableaus +############################################################################### + +# Butcher Tableaus, or more generally coefficient tables +# see Hairer & Wanner 1992, p. 134, 166 + +abstract Tableau{Name, S, T<:Real} +# Name is the name of the tableau/method (a symbol) +# S is the number of stages (an int) +# T is the type of the coefficients +# +# TODO: have a type parameter which specifies adaptive vs non-adaptive +# +# For all types of tableaus it assumes fields: +# order::(Int...) # order of the method(s) +# +# For Runge-Kutta methods it assumes fields: +# a::Matrix{T} # SxS matrix +# b::Matrix{T} # 1 or 2 x S matrix (fixed step/ adaptive) +# c::Vector{T} # S +# +# For a tableau: +# c1 | a_11 .... a_1s +# . | a_21 . . +# . | a_31 . . +# . | .... . . +# c_s | a_s1 ....... a_ss +# -----+-------------------- +# | b_1 ... b_s this is the one used for stepping +# | b'_1 ... b'_s this is the one used for error-checking + +Base.eltype{N,S,T}(b::Tableau{N,S,T}) = T +order(b::Tableau) = b.order +# Subtypes need to define a convert method to convert to a different +# eltype with signature: +Base.convert{Tnew<:Real}(::Type{Tnew}, tab::Tableau) = error("Define convert method for concrete Tableau types") ############################################################################### ## HELPER FUNCTIONS @@ -24,7 +61,9 @@ export ode4s, ode4ms, ode4 # estimator for initial step based on book # "Solving Ordinary Differential Equations I" by Hairer et al., p.169 -function hinit(F, x0, t0, tdir, p, reltol, abstol) +function hinit(F, x0, t0, tend, p, reltol, abstol) + tdir = sign(tend-t0) + tdir==0 && error("Zero time span") tau = max(reltol*norm(x0, Inf), abstol) d0 = norm(x0, Inf)/tau f0 = F(t0, x0) @@ -45,250 +84,110 @@ function hinit(F, x0, t0, tdir, p, reltol, abstol) pow = -(2. + log10(max(d1, d2)))/(p + 1.) h1 = 10.^pow end - h = min(100*h0, h1), f0 + return tdir*min(100*h0, h1), tdir, f0 end -############################################################################### -## NON-STIFF SOLVERS -############################################################################### - -# ODERKF based on -# -# ode45 adapted from http://users.powernet.co.uk/kienzle/octave/matcompat/scripts/ode_v1.11/ode45.m -# (a newer version (v1.15) can be found here https://sites.google.com/site/comperem/home/ode_solvers) -# -# Original Octave implementation: -# Marc Compere -# CompereM@asme.org -# created : 06 October 1999 -# modified: 17 January 2001 -# -function oderkf(F, x0, tspan, p, a, bs, bp; reltol = 1.0e-5, abstol = 1.0e-8, - norm=Base.norm, - minstep=abs(tspan[end] - tspan[1])/1e9, - maxstep=abs(tspan[end] - tspan[1])/2.5, - initstep=0.) - # see p.91 in the Ascher & Petzold reference for more infomation. - pow = 1/p # use the higher order to estimate the next step size - - c = sum(a, 2) # consistency condition - k = Array(typeof(x0), length(c)) - - # Initialization - t = tspan[1] - tfinal = tspan[end] - tdir = sign(tfinal - t) +# isoutofdomain takes the state and returns true if state is outside +# of the allowed domain. Used in adaptive step-control. +isoutofdomain(x) = isnan(x) + +function make_consistent_types(fn, y0, tspan, btab::Tableau) + # There are a few types involved in a call to a ODE solver which + # somehow need to be consistent: + # + # Et = eltype(tspan) + # Ey = eltype(y0) + # Ef = eltype(Tf) + # + # There are also the types of the containers, but they are not + # needed as `similar` is used to make containers. + # Tt = typeof(tspan) + # Ty = typeof(y0) # note, this can be a scalar + # Tf = typeof(F(tspan(1),y0)) # note, this can be a scalar + # + # Returns + # - Et: eltype of time, needs to be a real "continuous" type, at + # the moment a FloatingPoint + # - Eyf: suitable eltype of y and f(t,y) + # --> both of these are set to typeof(y0[1]/(tspan[end]-tspan[1])) + # - Ty: container type of y0 + # - btab: tableau with entries converted to Et + + # Needed interface: + # On components: /, - + # On container: eltype, promote_type + # On time container: eltype + + Ty = typeof(y0) + Eyf = typeof(y0[1]/(tspan[end]-tspan[1])) + + Et = eltype(tspan) + @assert Et<:Real + if !(Et<:FloatingPoint) + Et = promote_type(Et, Float64) + end - h = initstep - if h == 0. - # initial guess at a step size - h, k[1] = hinit(F, x0, t, tdir, p, reltol, abstol) - else - k[1] = F(t, x0) # first stage + # if all are Floats, make them the same + if Et<:FloatingPoint && Eyf<:FloatingPoint + Et = promote_type(Et, Eyf) + Eyf = Et end - h = tdir*min(h, maxstep) - x = x0 - tout = Array(typeof(t), 1) - tout[1] = t # first output time - xout = Array(typeof(x0), 1) - xout[1] = x # first output solution - while abs(t) != abs(tfinal) && abs(h) >= minstep - if abs(h) > abs(tfinal-t) - h = tfinal - t - end + !isleaftype(Et) && warn("The eltype(tspan) is not a concrete type! Change type of tspan for better performance.") + !isleaftype(Eyf) && warn("The eltype(y0/tspan[1]) is not a concrete type! Change type of y0 and/or tspan for better performance.") - #(p-1)th and pth order estimates - xs = x + h*bs[1]*k[1] - xp = x + h*bp[1]*k[1] - for j = 2:length(c) - dx = a[j,1]*k[1] - for i = 2:j-1 - dx += a[j,i]*k[i] - end - k[j] = F(t + h*c[j], x + h*dx) + btab_ = convert(Et, btab) + return Et, Eyf, Ty, btab_ +end - # compute the (p-1)th order estimate - xs = xs + h*bs[j]*k[j] - # compute the pth order estimate - xp = xp + h*bp[j]*k[j] - end +############################################################################### +## NON-STIFF SOLVERS +############################################################################### - # estimate the local truncation error - gamma1 = xs - xp +include("runge_kutta.jl") - # Estimate the error and the acceptable error - delta = norm(gamma1, Inf) # actual error - tau = max(reltol*norm(x,Inf),abstol) # allowable error +# ODE_MS Fixed-step, fixed-order multi-step numerical method +# with Adams-Bashforth-Moulton coefficients +function ode_ms(F, x0, tspan, order::Integer) + h = diff(tspan) + x = Array(typeof(x0), length(tspan)) + x[1] = x0 - # Update the solution only if the error is acceptable - if delta <= tau - t = t + h - x = xp # <-- using the higher order estimate is called 'local extrapolation' - push!(tout, t) - push!(xout, x) - - # Compute the slopes by computing the k[:,j+1]'th column based on the previous k[:,1:j] columns - # notes: k needs to end up as an Nxs, a is 7x6, which is s by (s-1), - # s is the number of intermediate RK stages on [t (t+h)] (Dormand-Prince has s=7 stages) - if c[end] == 1 - # Assign the last stage for x(k) as the first stage for computing x[k+1]. - # This is part of the Dormand-Prince pair caveat. - # k[:,7] has already been computed, so use it instead of recomputing it - # again as k[:,1] during the next step. - k[1] = k[end] - else - k[1] = F(t,x) # first stage + if 1 <= order <= 4 + b = ms_coefficients4 + else + b = zeros(order, order) + b[1:4, 1:4] = ms_coefficients4 + for s = 5:order + for j = 0:(s - 1) + # Assign in correct order for multiplication below + # (a factor depending on j and s) .* (an integral of a polynomial with -(0:s), except -j, as roots) + p_int = polyint(poly(diagm(-[0:j - 1; j + 1:s - 1]))) + b[s, s - j] = ((-1)^j / factorial(j) + / factorial(s - 1 - j) * polyval(p_int, 1)) end end - - # Update the step size - h = tdir*min(maxstep, 0.8*abs(h)*(tau/delta)^pow) - end # while (t < tfinal) & (h >= minstep) - - if abs(t) < abs(tfinal) - error("Step size grew too small. t=$t, h=$(abs(h)), x=$x") end - return tout, xout -end - - -# Bogacki–Shampine coefficients -const bs_coefficients = (3, - [ 0 0 0 0 - 1/2 0 0 0 - 0 3/4 0 0 - 2/9 1/3 4/9 0], - # 2nd order b-coefficients - [7/24 1/4 1/3 1/8], - # 3rd order b-coefficients - [2/9 1/3 4/9 0], - ) -ode23_bs(F, x0, tspan; kwargs...) = oderkf(F, x0, tspan, bs_coefficients...; kwargs...) - - -# Both the Dormand-Prince and Fehlberg 4(5) coefficients are from a tableau in -# U.M. Ascher, L.R. Petzold, Computer Methods for Ordinary Differential Equations -# and Differential-Agebraic Equations, Society for Industrial and Applied Mathematics -# (SIAM), Philadelphia, 1998 -# -# Dormand-Prince coefficients -const dp_coefficients = (5, - [ 0 0 0 0 0 0 - 1/5 0 0 0 0 0 - 3/40 9/40 0 0 0 0 - 44/45 -56/15 32/9 0 0 0 - 19372/6561 -25360/2187 64448/6561 -212/729 0 0 - 9017/3168 -355/33 46732/5247 49/176 -5103/18656 0 - 35/384 0 500/1113 125/192 -2187/6784 11/84], - # 4th order b-coefficients - [5179/57600 0 7571/16695 393/640 -92097/339200 187/2100 1/40], - # 5th order b-coefficients - [35/384 0 500/1113 125/192 -2187/6784 11/84 0], - ) -ode45_dp(F, x0, tspan; kwargs...) = oderkf(F, x0, tspan, dp_coefficients...; kwargs...) - - -# Fehlberg coefficients -const fb_coefficients = (5, - [ 0 0 0 0 0 - 1/4 0 0 0 0 - 3/32 9/32 0 0 0 - 1932/2197 -7200/2197 7296/2197 0 0 - 439/216 -8 3680/513 -845/4104 0 - -8/27 2 -3544/2565 1859/4104 -11/40], - # 4th order b-coefficients - [25/216 0 1408/2565 2197/4104 -1/5 0], - # 5th order b-coefficients - [16/135 0 6656/12825 28561/56430 -9/50 2/55], - ) -ode45_fb(F, x0, tspan; kwargs...) = oderkf(F, x0, tspan, fb_coefficients...; kwargs...) - - -# Cash-Karp coefficients -# Numerical Recipes in Fortran 77 -const ck_coefficients = (5, - [ 0 0 0 0 0 - 1/5 0 0 0 0 - 3/40 9/40 0 0 0 - 3/10 -9/10 6/5 0 0 - -11/54 5/2 -70/27 35/27 0 - 1631/55296 175/512 575/13824 44275/110592 253/4096], - # 4th order b-coefficients - [37/378 0 250/621 125/594 0 512/1771], - # 5th order b-coefficients - [2825/27648 0 18575/48384 13525/55296 277/14336 1/4], - ) -ode45_ck(F, x0, tspan; kwargs...) = oderkf(F, x0, tspan, ck_coefficients...; kwargs...) - - -# Fehlberg 7(8) coefficients -# Values from pag. 65, Fehlberg, Erwin. "Classical fifth-, sixth-, seventh-, and eighth-order Runge-Kutta formulas with stepsize control". -# National Aeronautics and Space Administration. -const fb_coefficients_78 = (8, - [ 0 0 0 0 0 0 0 0 0 0 0 0 - 2/27 0 0 0 0 0 0 0 0 0 0 0 - 1/36 1/12 0 0 0 0 0 0 0 0 0 0 - 1/24 0 1/8 0 0 0 0 0 0 0 0 0 - 5/12 0 -25/16 25/16 0 0 0 0 0 0 0 0 - 1/20 0 0 1/4 1/5 0 0 0 0 0 0 0 - -25/108 0 0 125/108 -65/27 125/54 0 0 0 0 0 0 - 31/300 0 0 0 61/225 -2/9 13/900 0 0 0 0 0 - 2 0 0 -53/6 704/45 -107/9 67/90 3 0 0 0 0 - -91/108 0 0 23/108 -976/135 311/54 -19/60 17/6 -1/12 0 0 0 - 2383/4100 0 0 -341/164 4496/1025 -301/82 2133/4100 45/82 45/164 18/41 0 0 - 3/205 0 0 0 0 -6/41 -3/205 -3/41 3/41 6/41 0 0 - -1777/4100 0 0 -341/164 4496/1025 -289/82 2193/4100 51/82 33/164 12/41 0 1], - # 7th order b-coefficients - [41/840 0 0 0 0 34/105 9/35 9/35 9/280 9/280 41/840 0 0], - # 8th order b-coefficients - [0 0 0 0 0 34/105 9/35 9/35 9/280 9/280 0 41/840 41/840], - ) -ode78_fb(F, x0, tspan; kwargs...) = oderkf(F, x0, tspan, fb_coefficients_78...; kwargs...) - -# Use Fehlberg version of ode78 by default -const ode78 = ode78_fb - -# Use Dormand-Prince version of ode45 by default -const ode45 = ode45_dp - -# Use Bogacki–Shampine version of ode23 by default -const ode23 = ode23_bs - - -# more higher-order embedded methods can be found in: -# P.J. Prince and J.R.Dormand: High order embedded Runge-Kutta formulae, Journal of Computational and Applied Mathematics 7(1), 1981. - - -#ODE4 Solve non-stiff differential equations, fourth order -# fixed-step Runge-Kutta method. -# -# [T,X] = ODE4(ODEFUN, X0, TSPAN) with TSPAN = [T0:H:TFINAL] -# integrates the system of differential equations x' = f(t,x) from time -# T0 to TFINAL in steps of H with initial conditions X0. Function -# ODEFUN(T,X) must return a column vector corresponding to f(t,x). Each -# row in the solution array X corresponds to a time returned in the -# column vector T. -function ode4(F, x0, tspan) - h = diff(tspan) - x = Array(typeof(x0), length(tspan)) - x[1] = x0 - - midxdot = Array(typeof(x0), 4) + # TODO: use a better data structure here (should be an order-element circ buffer) + xdot = similar(x) for i = 1:length(tspan)-1 - # Compute midstep derivatives - midxdot[1] = F(tspan[i], x[i]) - midxdot[2] = 2*F(tspan[i]+h[i]./2, x[i] + midxdot[1].*h[i]./2) - midxdot[3] = 2*F(tspan[i]+h[i]./2, x[i] + midxdot[2].*h[i]./2) - midxdot[4] = F(tspan[i]+h[i], x[i] + midxdot[3].*h[i]) - - # Integrate - x[i+1] = x[i] + 1/6 .*h[i].*sum(midxdot) + # Need to run the first several steps at reduced order + steporder = min(i, order) + xdot[i] = F(tspan[i], x[i]) + + x[i+1] = x[i] + for j = 1:steporder + x[i+1] += h[i]*b[steporder, j]*xdot[i-(steporder-1) + (j-1)] + end end return vcat(tspan), x end +# Use order 4 by default +ode4ms(F, x0, tspan) = ode_ms(F, x0, tspan, 4) +ode5ms(F, x0, tspan) = ODE.ode_ms(F, x0, tspan, 5) + ############################################################################### ## STIFF SOLVERS ############################################################################### @@ -351,14 +250,15 @@ function ode23s(F, y0, tspan; reltol = 1.0e-5, abstol = 1.0e-8, # initialization t = tspan[1] tfinal = tspan[end] - tdir = sign(tfinal - t) + h = initstep if h == 0. - # initial guess at a step size - h, F0 = hinit(F, y0, t, tdir, 3, reltol, abstol) + # initial guess at a step size + h, tdir, F0 = hinit(F, y0, t, tfinal, 3, reltol, abstol) else - F0 = F(t,y0) + tdir = sign(tfinal - t) + F0 = F(t,y0) end h = tdir*min(h, maxstep) @@ -366,7 +266,7 @@ function ode23s(F, y0, tspan; reltol = 1.0e-5, abstol = 1.0e-8, tout = Array(typeof(t), 1) tout[1] = t # first output time yout = Array(typeof(y0), 1) - yout[1] = copy(y) # first output solution + yout[1] = deepcopy(y) # first output solution J = jac(t,y) # get Jacobian of F wrt y @@ -522,45 +422,5 @@ const ms_coefficients4 = [ 1 0 0 0 5/12 -4/3 23/12 0 -9/24 37/24 -59/24 55/24] -# ODE_MS Fixed-step, fixed-order multi-step numerical method -# with Adams-Bashforth-Moulton coefficients -function ode_ms(F, x0, tspan, order::Integer) - h = diff(tspan) - x = Array(typeof(x0), length(tspan)) - x[1] = x0 - - if 1 <= order <= 4 - b = ms_coefficients4 - else - b = zeros(order, order) - b[1:4, 1:4] = ms_coefficients4 - for s = 5:order - for j = 0:(s - 1) - # Assign in correct order for multiplication below - # (a factor depending on j and s) .* (an integral of a polynomial with -(0:s), except -j, as roots) - p_int = polyint(poly(diagm(-[0:j - 1; j + 1:s - 1]))) - b[s, s - j] = ((-1)^j / factorial(j) - / factorial(s - 1 - j) * polyval(p_int, 1)) - end - end - end - - # TODO: use a better data structure here (should be an order-element circ buffer) - xdot = similar(x) - for i = 1:length(tspan)-1 - # Need to run the first several steps at reduced order - steporder = min(i, order) - xdot[i] = F(tspan[i], x[i]) - - x[i+1] = x[i] - for j = 1:steporder - x[i+1] += h[i]*b[steporder, j]*xdot[i-(steporder-1) + (j-1)] - end - end - return vcat(tspan), x -end - -# Use order 4 by default -ode4ms(F, x0, tspan) = ode_ms(F, x0, tspan, 4) end # module ODE diff --git a/src/runge_kutta.jl b/src/runge_kutta.jl new file mode 100644 index 000000000..9b13e0944 --- /dev/null +++ b/src/runge_kutta.jl @@ -0,0 +1,480 @@ +# Explicit Runge-Kutta solvers +############################## +# (Hairer & Wanner 1992 p.134, p.165-169) + +########################################### +# Tableaus for explicit Runge-Kutta methods +########################################### + +immutable TableauRKExplicit{Name, S, T} <: Tableau{Name, S, T} + order::(@compat(Tuple{Vararg{Int}})) # the order of the methods + a::Matrix{T} + # one or several row vectors. First row is used for the step, + # second for error calc. + b::Matrix{T} + c::Vector{T} + function TableauRKExplicit(order,a,b,c) + @assert isa(S,Integer) + @assert isa(Name,Symbol) + @assert c[1]==0 + @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. + new(order,a,b,c) + end +end +function TableauRKExplicit{T}(name::Symbol, order::(@compat(Tuple{Vararg{Int}})), + a::Matrix{T}, b::Matrix{T}, c::Vector{T}) + TableauRKExplicit{name,length(c),T}(order, a, b, c) +end +function TableauRKExplicit(name::Symbol, order::(@compat(Tuple{Vararg{Int}})), T::Type, + a::Matrix, b::Matrix, c::Vector) + TableauRKExplicit{name,length(c),T}(order, convert(Matrix{T},a), + convert(Matrix{T},b), convert(Vector{T},c) ) +end +conv_field{T,N}(D,a::Array{T,N}) = convert(Array{D,N}, a) +function Base.convert{Tnew<:Real,Name,S,T}(::Type{Tnew}, tab::TableauRKExplicit{Name,S,T}) + # Converts the tableau coefficients to the new type Tnew + newflds = () + @compat for n in fieldnames(tab) + fld = getfield(tab,n) + if eltype(fld)==T + newflds = tuple(newflds..., conv_field(Tnew, fld)) + else + newflds = tuple(newflds..., fld) + end + end + TableauRKExplicit{Name,S,Tnew}(newflds...) # TODO: could this be done more generically in a type-stable way? +end + + +isexplicit(b::TableauRKExplicit) = istril(b.a) # Test whether it's an explicit method +isadaptive(b::TableauRKExplicit) = size(b.b, 1)==2 + + +# First same as last. Means ks[:,end]=ks_nextstep[:,1], c.f. H&W p.167 +isFSAL(btab::TableauRKExplicit) = btab.a[end,:]==btab.b[1,:] && btab.c[end]==1 # the latter is not needed really + +## Tableaus for explicit RK methods +# Fixed step: +const bt_feuler = TableauRKExplicit(:feuler,(1,), Rational{Int64}, + zeros(Int,1,1), + [1]', + [0] + ) +const bt_midpoint = TableauRKExplicit(:midpoint,(2,), Rational{Int64}, + [0 0 + 1//2 0], + [0, 1]', + [0, 1//2] + ) +const bt_heun = TableauRKExplicit(:heun,(2,), Rational{Int64}, + [0 0 + 1 0], + [1//2, 1//2]', + [0, 1]) + +const bt_rk4 = TableauRKExplicit(:rk4,(4,),Rational{Int64}, + [0 0 0 0 + 1//2 0 0 0 + 0 1//2 0 0 + 0 0 1 0], + [1//6, 1//3, 1//3, 1//6]', + [0, 1//2, 1//2, 1]) + +# Adaptive step: +# Heun Euler https://en.wikipedia.org/wiki/Runge–Kutta_methods +const bt_rk21 = TableauRKExplicit(:heun_euler,(2,1), Rational{Int64}, + [0 0 + 1 0], + [1//2 1//2 + 1 0], + [0, 1]) + +# Bogacki–Shampine coefficients +const bt_rk23 = TableauRKExplicit(:bogacki_shampine,(2,3), Rational{Int64}, + [0 0 0 0 + 1/2 0 0 0 + 0 3/4 0 0 + 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] + ) + +# Fehlberg https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta%E2%80%93Fehlberg_method +const bt_rk45 = TableauRKExplicit(:fehlberg,(4,5),Rational{Int64}, + [ 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 ], + [25//216 0 1408//2565 2197//4104 -1//5 0 + 16//135 0 6656//12825 28561//56430 -9//50 2//55], + [0, 1//4, 3//8, 12//13, 1, 1//2]) + +# Dormand-Prince https://en.wikipedia.org/wiki/Dormand%E2%80%93Prince_method +const bt_dopri5 = TableauRKExplicit(:dopri, (5,4), Rational{Int64}, + [0 0 0 0 0 0 0 + 1//5 0 0 0 0 0 0 + 3//40 9//40 0 0 0 0 0 + 44//45 -56//15 32//9 0 0 0 0 + 19372//6561 -25360//2187 64448//6561 -212//729 0 0 0 + 9017//3168 -355//33 46732//5247 49//176 -5103//18656 0 0 + 35//384 0 500//1113 125//192 -2187//6784 11//84 0], + [35//384 0 500//1113 125//192 -2187//6784 11//84 0 + 5179//57600 0 7571//16695 393//640 -92097//339200 187//2100 1//40], + [0, 1//5, 3//10, 4//5, 8//9, 1, 1] + ) + +# Fehlberg 7(8) coefficients +# Values from pag. 65, Fehlberg, Erwin. "Classical fifth-, sixth-, seventh-, and eighth-order Runge-Kutta formulas with stepsize control". +# National Aeronautics and Space Administration. +const bt_feh78 = TableauRKExplicit(:feh78, (7,8), Rational{Int64}, + [ 0 0 0 0 0 0 0 0 0 0 0 0 0 + 2//27 0 0 0 0 0 0 0 0 0 0 0 0 + 1//36 1//12 0 0 0 0 0 0 0 0 0 0 0 + 1//24 0 1//8 0 0 0 0 0 0 0 0 0 0 + 5//12 0 -25//16 25//16 0 0 0 0 0 0 0 0 0 + 1//20 0 0 1//4 1//5 0 0 0 0 0 0 0 0 + -25//108 0 0 125//108 -65//27 125//54 0 0 0 0 0 0 0 + 31//300 0 0 0 61//225 -2//9 13//900 0 0 0 0 0 0 + 2 0 0 -53//6 704//45 -107//9 67//90 3 0 0 0 0 0 + -91//108 0 0 23//108 -976//135 311//54 -19//60 17//6 -1//12 0 0 0 0 + 2383//4100 0 0 -341//164 4496//1025 -301//82 2133//4100 45//82 45//164 18//41 0 0 0 + 3//205 0 0 0 0 -6//41 -3//205 -3//41 3//41 6//41 0 0 0 + -1777//4100 0 0 -341//164 4496//1025 -289//82 2193//4100 51//82 33//164 12//41 0 1 0], + [41//840 0 0 0 0 34//105 9//35 9//35 9//280 9//280 41//840 0 0 + 0 0 0 0 0 34//105 9//35 9//35 9//280 9//280 0 41//840 41//840], + [0, 2//27, 1//9, 1//6 , 5//12, 1//2 , 5//6 , 1//6 , 2//3 , 1//3 , 1 , 0, 1] + ) + + +################################ +# Fixed step Runge-Kutta methods +################################ + +# TODO: iterator method +ode1(fn, y0, tspan) = oderk_fixed(fn, y0, tspan, bt_feuler) +ode2_midpoint(fn, y0, tspan) = oderk_fixed(fn, y0, tspan, bt_midpoint) +ode2_heun(fn, y0, tspan) = oderk_fixed(fn, y0, tspan, bt_heun) +ode4(fn, y0, tspan) = oderk_fixed(fn, y0, tspan, bt_rk4) + +function oderk_fixed(fn, y0, tspan, btab::TableauRKExplicit) + # Non-arrays y0 treat as scalar + fn_(t, y) = [fn(t, y[1])] + t,y = oderk_fixed(fn_, [y0], tspan, btab) + return t, vcat_nosplat(y) +end +function oderk_fixed{N,S}(fn, y0::AbstractVector, tspan, + btab_::TableauRKExplicit{N,S}) + # TODO: instead of AbstractVector use a Holy-trait + + # Needed interface: + # On components: + # On y0 container: length, deepcopy, similar, setindex! + # On time container: getindex, convert. length + + Et, Eyf, Ty, btab = make_consistent_types(fn, y0, tspan, btab_) + dof = length(y0) + + ys = Array(Ty, length(tspan)) + allocate!(ys, y0, dof) + ys[1] = deepcopy(y0) + + tspan = convert(Vector{Et}, tspan) + # work arrays: + ks = Array(Ty, S) + # allocate!(ks, y0, dof) # no need to allocate as fn is not in-place + ytmp = similar(y0, Eyf, dof) + for i=1:length(tspan)-1 + dt = tspan[i+1]-tspan[i] + ys[i+1][:] = ys[i] + for s=1:S + calc_next_k!(ks, ytmp, ys[i], s, fn, tspan[i], dt, dof, btab) + for d=1:dof + ys[i+1][d] += dt * btab.b[s]*ks[s][d] + end + end + end + return tspan, ys +end + +############################## +# Adaptive Runge-Kutta methods +############################## + +ode21(fn, y0, tspan; kwargs...) = oderk_adapt(fn, y0, tspan, bt_rk21; kwargs...) +ode23(fn, y0, tspan; kwargs...) = oderk_adapt(fn, y0, tspan, bt_rk23; kwargs...) +ode45_fe(fn, y0, tspan; kwargs...) = oderk_adapt(fn, y0, tspan, bt_rk45; kwargs...) +ode45_dp(fn, y0, tspan; kwargs...) = oderk_adapt(fn, y0, tspan, bt_dopri5; kwargs...) +# 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...) + +function oderk_adapt(fn, y0, tspan, btab::TableauRKExplicit; kwords...) + # For y0 which don't support indexing. + fn_ = (t, y) -> [fn(t, y[1])] + t,y = oderk_adapt(fn_, [y0], tspan, btab; kwords...) + return t, vcat_nosplat(y) +end +function oderk_adapt{N,S}(fn, y0::AbstractVector, tspan, btab_::TableauRKExplicit{N,S}; + reltol = 1.0e-5, abstol = 1.0e-8, + norm=Base.norm, + minstep=abs(tspan[end] - tspan[1])/1e9, + maxstep=abs(tspan[end] - tspan[1])/2.5, + initstep=0., + points=:all + ) + # Needed interface: + # On components: + # - 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") + + Et, Eyf, Ty, btab = make_consistent_types(fn, y0, tspan, btab_) + # parameters + order = minimum(btab.order) + timeout_const = 5 # after step reduction do not increase step for + # timeout_const steps + + ## Initialization + dof = length(y0) + tspan = convert(Vector{Et}, tspan) + t = tspan[1] + tstart = tspan[1] + tend = tspan[end] + + # work arrays: + y = similar(y0, Eyf, dof) # y at time t + y[:] = y0 + ytrial = similar(y0, Eyf, dof) # trial solution at time t+dt + yerr = similar(y0, Eyf, dof) # error of trial solution + ks = Array(Ty, S) + # allocate!(ks, y0, dof) # no need to allocate as fn is not in-place + ytmp = similar(y0, Eyf, dof) + + # output ys + nsteps_fixed = length(tspan) # these are always output + ys = Array(Ty, nsteps_fixed) + allocate!(ys, y0, dof) + ys[1] = y0 + + # Option points determines where solution is returned: + if points==:all + tspan_fixed = tspan + tspan = Et[tstart] + iter_fixed = 2 # index into tspan_fixed + sizehint!(tspan, nsteps_fixed) + elseif points!=:specified + error("Unrecognized option points==$points") + end + # Time + dt, tdir, ks[1] = hinit(fn, y, tstart, tend, order, reltol, abstol) # sets ks[1]=f0 + if initstep!=0 + dt = sign(initstep)==tdir ? initstep : error("initstep has wrong sign.") + end + # Diagnostics + dts = Et[] + errs = Float64[] + steps = [0,0] # [accepted, rejected] + + ## Integration loop + laststep = false + timeout = 0 # for step-control + iter = 2 # the index into tspan and ys + while true + # do one step (assumes ks[1]==f0) + rk_embedded_step!(ytrial, yerr, ks, ytmp, y, fn, t, dt, dof, btab) + # Check error and find a new step size: + err, newdt, timeout = stepsize_hw92!(dt, tdir, y, ytrial, yerr, order, timeout, + dof, abstol, reltol, maxstep, norm) + + if err<=1.0 # accept step + # diagnostics + steps[1] +=1 + push!(dts, dt) + push!(errs, err) + + # Output: + f0 = ks[1] + f1 = isFSAL(btab) ? ks[S] : fn(t+dt, ytrial) + if points==:specified + # interpolate onto given output points + while iter-1= tdir*tend + dt = tend-t + laststep = true # next step is the last, if it succeeds + end + elseif abs(newdt)0 no step size increase is allowed, timeout is + # decremented in here. + # + # Returns the error, newdt and the number of timeout-steps + # + # TODO: + # - allow component-wise reltol and abstol? + # - allow other norms + + # Needed interface: + # On components: isoutofdomain + # On y0 container: norm, get/setindex + + timout_after_nan = 5 + fac = [0.8, 0.9, 0.25^(1/(order+1)), 0.38^(1/(order+1))][1] + facmax = 5.0 # maximal step size increase. 1.5-5 + facmin = 1./facmax # maximal step size decrease. ? + + # in-place calculate xerr./tol + for d=1:dof + # if outside of domain (usually NaN) then make step size smaller by maximum + isoutofdomain(xtrial[d]) && return 10., dt*facmin, timout_after_nan + xerr[d] = xerr[d]/(abstol + max(norm(x0[d]), norm(xtrial[d]))*reltol) # Eq 4.10 + end + err = norm(xerr, 2) # Eq. 4.11 + newdt = min(maxstep, tdir*dt*max(facmin, fac*(1/err)^(1/(order+1)))) # Eq 4.13 modified + if timeout>0 + newdt = min(newdt, dt) + timeout -= 1 + end + return err, tdir*newdt, timeout +end + +function calc_next_k!{Ty}(ks::Vector, ytmp::Ty, y, s, fn, t, dt, dof, btab) + # Calculates the next ks and puts it into ks[s] + # - ks and ytmp are modified inside this function. + + # Needed interface: + # On components: +, * + # On y0 container: setindex!, getindex, fn + + ytmp[:] = y + for ss=1:s-1, d=1:dof + ytmp[d] += dt * ks[ss][d] * btab.a[s,ss] + end + ks[s] = fn(t + btab.c[s]*dt, ytmp)::Ty + nothing +end + +# Helper functions: +function allocate!{T}(vec::Vector{T}, y0, dof) + # Allocates all vectors inside a Vector{Vector} using the same + # kind of container as y0 has and element type eltype(eltype(vec)). + for s=1:length(vec) + vec[s] = similar(y0, eltype(T), dof) + end +end +function index_or_push!(vec, i, val) + # Fills in the vector until there is no space, then uses push! + # instead. + if length(vec)>=i + vec[i] = val + else + push!(vec, val) + end +end +vcat_nosplat(y) = eltype(y[1])[el[1] for el in y] # Does vcat(y...) without the splatting + +function hermite_interp!(y, tquery,t,dt,y0,y1,f0,f1) + # For dense output see Hairer & Wanner p.190 using Hermite + # interpolation. Updates y in-place. + # + # f_0 = f(x_0 , y_0) , f_1 = f(x_0 + h, y_1 ) + # this is O(3). TODO for higher order. + + theta = (tquery-t)/dt + for i=1:length(y0) + y[i] = ((1-theta)*y0[i] + theta*y1[i] + theta*(theta-1) * + ((1-2*theta)*(y1[i]-y0[i]) + (theta-1)*dt*f0[i] + theta*dt*f1[i]) ) + end + nothing +end +function hermite_interp(tquery,t,dt,y0,y1,f0,f1) + # Returns the y instead of in-place + y = similar(y0) + hermite_interp!(y,tquery,t,dt,y0,y1,f0,f1) + return y +end diff --git a/test/interface-tests.jl b/test/interface-tests.jl new file mode 100644 index 000000000..bc25056e9 --- /dev/null +++ b/test/interface-tests.jl @@ -0,0 +1,79 @@ +# Here are tests which test what interface the solvers require. + +################################################################################ +# This is to test a scalar-like state variable +# (due to @acroy: https://gist.github.com/acroy/28be4f2384d01f38e577) +const delta0 = 0. +const V0 = 1. +const g0 = 0. + +# define custom type ... +immutable CompSol + rho::Matrix{Complex128} + x::Float64 + p::Float64 + + CompSol(r,x,p) = new(copy(r),x,p) +end + +# ... which has to support the following operations +# to work with odeX +Base.norm(y::CompSol, p::Float64) = maximum([Base.norm(y.rho, p) abs(y.x) abs(y.p)]) +Base.norm(y::CompSol) = norm(y::CompSol, 2.0) + ++(y1::CompSol, y2::CompSol) = CompSol(y1.rho+y2.rho, y1.x+y2.x, y1.p+y2.p) +-(y1::CompSol, y2::CompSol) = CompSol(y1.rho-y2.rho, y1.x-y2.x, y1.p-y2.p) +*(y1::CompSol, s::Real) = CompSol(y1.rho*s, y1.x*s, y1.p*s) +*(s::Real, y1::CompSol) = y1*s +/(y1::CompSol, s::Real) = CompSol(y1.rho/s, y1.x/s, y1.p/s) + +### new for PR #68 +Base.abs(y::CompSol) = norm(y, 2.) # TODO not needed anymore once https://github.com/JuliaLang/julia/pull/11043 is in current stable julia +Base.zero(::Type{CompSol}) = CompSol(complex(zeros(2,2)), 0., 0.) +ODE.isoutofdomain(y::CompSol) = any(isnan, vcat(y.rho[:], y.x, y.p)) + +# Because the new RK solvers wrap scalars in an array and because of +# https://github.com/JuliaLang/julia/issues/11053 these are also needed: +.+(y1::CompSol, y2::CompSol) = CompSol(y1.rho+y2.rho, y1.x+y2.x, y1.p+y2.p) +.-(y1::CompSol, y2::CompSol) = CompSol(y1.rho-y2.rho, y1.x-y2.x, y1.p-y2.p) +.*(y1::CompSol, s::Real) = CompSol(y1.rho*s, y1.x*s, y1.p*s) +.*(s::Real, y1::CompSol) = y1*s +./(y1::CompSol, s::Real) = CompSol(y1.rho/s, y1.x/s, y1.p/s) + + +################################################################################ + +# define RHSs of differential equations +# delta, V and g are parameters +function rhs(t, y, delta, V, g) + H = [[-delta/2 V]; [V delta/2]] + + rho_dot = -im*H*y.rho + im*y.rho*H + x_dot = y.p + p_dot = -y.x + + return CompSol( rho_dot, x_dot, p_dot) +end + +# inital conditons +rho0 = zeros(2,2); +rho0[1,1]=1.; +y0 = CompSol(complex(rho0), 2., 1.) + +# solve ODEs +endt = 2*pi; + +t,y1 = ODE.ode45((t,y)->rhs(t, y, delta0, V0, g0), y0, [0., endt]) # used as reference +print("Testing interface for scalar-like state... ") +for solver in solvers + # these only work with some Array-like interface defined: + if solver in [ODE.ode23s, ODE.ode4s_s, ODE.ode4s_kr] + continue + end + t,y2 = solver((t,y)->rhs(t, y, delta0, V0, g0), y0, linspace(0., endt, 500)) + @test norm(y1[end]-y2[end])<0.1 +end +println("ok.") + +################################################################################ +# TODO: test a vector-like state variable, i.e. one which can be indexed. diff --git a/test/runtests.jl b/test/runtests.jl index 575926dca..64f57e0b7 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,31 +3,35 @@ using Base.Test tol = 1e-2 -ode5ms(F, x0, tspan) = ODE.ode_ms(F, x0, tspan, 5) - solvers = [ - ODE.ode23, - - ODE.ode4, - ODE.ode45_dp, - ODE.ode45_fb, - ODE.ode45_ck, - - ODE.ode23s, - - ODE.ode4ms, - ode5ms, - ODE.ode4s_s, - ODE.ode4s_kr, + ## Non-stiff + # fixed step + ODE.ode1, + ODE.ode2_midpoint, + ODE.ode2_heun, + ODE.ode4, + ODE.ode4ms, + ODE.ode5ms, + # adaptive +# ODE.ode21, # this fails on Travis with 0.4?! TODO revert once fixed. + ODE.ode23, + ODE.ode45_dp, + ODE.ode45_fe, + ODE.ode78, - ODE.ode78_fb] + ## Stiff + # fixed-step + ODE.ode4s_s, + ODE.ode4s_kr, + # adaptive + ODE.ode23s] for solver in solvers println("using $solver") # dy # -- = 6 ==> y = 6t # dt - t,y=solver((t,y)->6, 0., [0:.1:1;]) + t,y=solver((t,y)->6.0, 0., [0:.1:1;]) @test maximum(abs(y-6t)) < tol # dy @@ -35,6 +39,7 @@ for solver in solvers # dt t,y=solver((t,y)->2t, 0., [0:.001:1;]) @test maximum(abs(y-t.^2)) < tol + # dy # -- = y ==> y = y0*e.^t @@ -74,5 +79,6 @@ let 0.9999999791665050] # reference solution at tspan[2] @test norm(refsol-y[end], Inf) < 2e-10 end +include("interface-tests.jl") println("All looks OK")