diff --git a/src/algorithms.jl b/src/algorithms.jl deleted file mode 100644 index 54e21e547..000000000 --- a/src/algorithms.jl +++ /dev/null @@ -1,306 +0,0 @@ -################################ -# Fixed step Runge-Kutta methods -################################ - -type Step - t; y; dy; dt -end - -# 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])/1e18, - 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 - islaststep = abs(t+dt-tend)<=eps(tend) ? true : 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 - islaststep = 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) -function hermite_interp!(y,t,step0::Step,step1::Step) - # 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. - - y0, y1 = step0.y, step1.y - dy0, dy1 = step0.dy, step1.dy - - dt = step1.t-step0.t - theta = (t-step0.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*dy0[i] + theta*dt*dy1[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/src/dense.jl b/src/dense.jl index afdd98460..aa255474b 100644 --- a/src/dense.jl +++ b/src/dense.jl @@ -30,8 +30,8 @@ function DenseProblem(args...; tspan = [Inf], points = :all, method = bt_feuler, end -# create an instance of the zero-th step for RKProblem -function Step(problem :: AbstractProblem) +# create an instance of the zero-th step for a problem +function Step(problem) t0 = problem.t0 y0 = problem.y0 dy0 = problem.F(t0,y0) diff --git a/src/iterators.jl b/src/iterators.jl index 2528916dd..1fd6de08f 100644 --- a/src/iterators.jl +++ b/src/iterators.jl @@ -1,313 +1,4 @@ -#################### -# Iterator methods # -#################### +import Base: start, next, done, call -# common structures and functions - -const timeout_const = 5 - -type TempArrays - y; ks; yerr -end - -type State - t; dt; y; tmp :: TempArrays - timeout -end - -abstract AbstractProblem - - -##################### -# Fixed step method # -##################### - -# type FixedState <: AbstractState -# t; dt; y; tmp :: TempArrays -# end - - -immutable FixedProblem <: AbstractProblem - F - method - y0 - t0 - dt0 - tstop -end - - -# # generates the iterator from TableauRKExplicit -# function call(tab::TableauRKExplicit, fn, y0, t0, dt0; tstop = Inf) -# return FixedProblem(fn, tab, y0, t0, dt0, tstop) -# end - -# function start(problem :: FixedProblem) -# t0 = problem.t0 -# dt0 = problem.dt0 -# y0 = problem.y0 -# tmp = TempArrays(similar(y0), Array(typeof(y0), S(problem.btab)), Void()) -# return FixedState(t0, dt0, y0, tmp) -# end - - -function next(prob :: FixedProblem, state :: State) - dof = length(state.y) - for s=1:S(prob.method) - calc_next_k!(state.tmp, s, state, prob) - for d=1:dof - state.y[d] += state.dt * prob.method.b[s]*state.tmp.ks[s][d] - end - end - state.t += state.dt - return ((state.t,state.y), state) -end - - -# function done(prob :: FixedProblem, state :: FixedState) -# return state.t >= prob.tstop # || abs(state.dt) < prob.minstep -# end - - -######################## -# Adaptive step method # -######################## - -immutable AdaptiveProblem <: AbstractProblem - F - method - y0 - t0 - dt0 - tstop - reltol - abstol - minstep - maxstep -end - - -# function AdaptiveProblem(fn, y0, t0; tstop = Inf, method = bt_feuler, reltol = 1e-5, abstol = 1e-5, minstep = 1e-10, maxstep = 1/minstep) -# return AdaptiveProblem{:adaptive}(fn, method, y0, t0, tstop, reltol, abstol, minstep, maxstep) -# end - - -function solver(F, y0, t0; - tstop = Inf, - method = bt_feuler, - reltol = 1e-5, - abstol = 1e-5, - minstep = 1e-10, - maxstep = 1/minstep, - dt0 = hinit(F, y0, t0, tstop, method, reltol, abstol)) - if isadaptive(method) - return AdaptiveProblem(F, method, y0, t0, dt0, tstop, reltol, abstol, minstep, maxstep) - else - return FixedProblem(F, method, y0, t0, dt0, tstop) - end -end - - -function start(problem :: AbstractProblem) - t0, dt0, y0 = problem.t0, problem.dt0, problem.y0 - - tmp = TempArrays(similar(y0), Array(typeof(y0), S(problem.method)), similar(y0)) - tmp.ks[1] = problem.F(t0,y0) # we assume that ks[1] is already initialized - - timeout = 0 # for step control - return State(t0,dt0,y0,tmp,timeout) -end - - -function done(prob :: AbstractProblem, state :: State) - return state.t >= prob.tstop || state.dt < prob.minstep -end - - -function next(prob :: AdaptiveProblem, state :: State) - - # the initial values - dt = state.dt # dt is the previous stepisze, it is - # modified inside the loop - timeout = state.timeout - - # for aesthetical reasons we extract the temporary componen - tmp = state.tmp - - # The while loop continues until we either find a stepsize which - # leads to a small enough error or the stepsize reaches - # prob.minstep - - while true - - # do one step (assumes ks[1]==f0), changes only tmp - err, newdt, timeout = rk_trial_step!(tmp, state, prob, dt, timeout) - - if abs(newdt) < prob.minstep # minimum step size reached, break - println("Warning: dt < minstep. Stopping.") - # passing the newdt to state will result in done() - state.dt = newdt - break - end - - if err > 1 # error is too large, repeat the step with smaller dt - # redo step with smaller dt and reset the timeout - dt = newdt - timeout = timeout_const - else - # step is accepted - - # preload ks[1] for the next step - if isFSAL(prob.method) - tmp.ks[1] = tmp.ks[S(prob.method)] - else - tmp.ks[1] = prob.F(state.t+dt, state.tmp.y) - end - - # Swap bindings of y and ytrial, avoids one copy - state.y, state.tmp.y = state.tmp.y, state.y - - # Update state with the data from the step we have just - # made: - state.t += dt - state.dt = newdt - state.timeout = timeout - break - end - end - return ((state.t,state.y),state) -end - - -function rk_trial_step!(tmp, state, prob, dt, timeout) - - # tmp.y and tmp.yerr and tmp.ks are updated after this step - rk_embedded_step!(tmp, state, prob, dt) - - # changes tmp.yerr (via in place update) - err, newdt, timeout = stepsize_hw92!(tmp, state, prob, dt, timeout) - - return err, newdt, timeout -end - - -function rk_embedded_step!(tmp :: TempArrays, state :: State, prob :: AdaptiveProblem, dt) - # Does one embedded R-K step updating ytrial, yerr and ks. - # - # Assumes that ks[:,1] is already calculated! - # - # Modifies tmp.y and tmp.yerr only - - y = state.y - dof = length(y) - b = prob.method.b - - tmp.y[:] = 0 - tmp.yerr[:] = 0 - - for d=1:dof - tmp.y[d] += b[1,1]*tmp.ks[1][d] - tmp.yerr[d] += b[2,1]*tmp.ks[1][d] - end - - for s=2:S(prob.method) - calc_next_k!(state.tmp, s, state, prob) - for d=1:dof - tmp.y[d] += b[1,s]*tmp.ks[s][d] - tmp.yerr[d] += b[2,s]*tmp.ks[s][d] - end - end - - for d=1:dof - tmp.yerr[d] = dt * (tmp.y[d]-tmp.yerr[d]) - tmp.y[d] = y[d] + dt * tmp.y[d] - end -end - - -function stepsize_hw92!(tmp, state, prob, dt, timeout) - # Estimates the error and a new step size following Hairer & - # Wanner 1992, p167 (with some modifications) - # - # If timeout>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 - - order = minimum(prob.method.order) - 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. ? - dof = length(state.y) - - # in-place calculate yerr./tol - for d=1:dof - - # if outside of domain (usually NaN) then make step size smaller by maximum - if isoutofdomain(tmp.y[d]) - return 10., dt*facmin, timout_after_nan - end - - tmp.yerr[d] = tmp.yerr[d]/(prob.abstol + max(norm(prob.y0[d]), norm(tmp.y[d]))*prob.reltol) # Eq 4.10 - end - - err = norm(tmp.yerr, 2) # Eq. 4.11 - newdt = min(prob.maxstep, 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, newdt, timeout -end - - -function hinit(F, y0, t0, tstop, method, reltol, abstol) - # Returns first step size - order = minimum(method.order) - tau = max(reltol*norm(y0, Inf), abstol) - d0 = norm(y0, Inf)/tau - f0 = F(t0, y0) - d1 = norm(f0, Inf)/tau - if d0 < 1e-5 || d1 < 1e-5 - h0 = 1e-6 - else - h0 = 0.01*(d0/d1) - end - # perform Euler step - y1 = y0 + h0*f0 - f1 = F(t0 + h0, y1) - # estimate second derivative - d2 = norm(f1 - f0, Inf)/(tau*h0) - if max(d1, d2) <= 1e-15 - h1 = max(1e-6, 1e-3*h0) - else - pow = -(2 + log10(max(d1, d2)))/(order+1) - h1 = 10^pow - end - return min(100*h0, h1, tstop-t0) -end - - -# For clarity we pass the TempArrays part of the state separately, -# this is the only part of state that can be changed here -function calc_next_k!(tmp :: TempArrays, i, state :: State, prob :: AbstractProblem) - dof = length(state.y) - t, dt, a, c = state.t, state.dt, prob.method.a, prob.method.c - - tmp.y[:] = state.y - for j=1:i-1 - for d=1:dof - tmp.y[d] += dt * tmp.ks[j][d] * a[i,j] - end - end - tmp.ks[i] = prob.F(t + c[i]*dt, tmp.y) - - nothing -end +include("rk.jl") +include("dense.jl") diff --git a/src/rk.jl b/src/rk.jl new file mode 100644 index 000000000..a75c04c44 --- /dev/null +++ b/src/rk.jl @@ -0,0 +1,290 @@ +# This file contains the implementation of explicit Runkge-Kutta +# solver from (Hairer & Wanner 1992 p.134, p.165-169). + +# include the Butcher tableaus. +include("tableaus.jl") + +#################### +# Iterator methods # +#################### + +# common structures and functions + +type TempArrays + y; ks; yerr +end + +type State + t; dt; y; tmp :: TempArrays + timeout +end + + +immutable Problem{MethodType} + F + method + y0 + t0 + dt0 + tstop + reltol + abstol + minstep + maxstep +end + + +function solver(F, y0, t0; + tstop = Inf, + method = bt_feuler, + reltol = 1e-5, + abstol = 1e-5, + minstep = 1e-10, + maxstep = 1/minstep, + dt0 = hinit(F, y0, t0, tstop, method, reltol, abstol)) + + if isadaptive(method) + methodtype = :adaptive + else + methodtype = :fixed + end + + return Problem{:fixed}(F, method, y0, t0, dt0, tstop, reltol, abstol, minstep, maxstep) + +end + + +function call(btab :: TableauRKExplicit, args...; opt_args...) + return solver(args...; opt_args..., method = btab) +end + + +function start(problem :: Problem) + t0, dt0, y0 = problem.t0, problem.dt0, problem.y0 + + tmp = TempArrays(similar(y0), Array(typeof(y0), S(problem.method)), similar(y0)) + tmp.ks[1] = problem.F(t0,y0) # we assume that ks[1] is already initialized + + timeout = 0 # for step control + return State(t0,dt0,y0,tmp,timeout) +end + + +function done(prob :: Problem, state :: State) + return state.t >= prob.tstop || state.dt < prob.minstep +end + +##################### +# Fixed step method # +##################### + +function next(prob :: Problem{:fixed}, state :: State) + dof = length(state.y) + for s=1:S(prob.method) + calc_next_k!(state.tmp, s, state, prob) + for d=1:dof + state.y[d] += state.dt * prob.method.b[s]*state.tmp.ks[s][d] + end + end + state.t += state.dt + return ((state.t,state.y), state) +end + +######################## +# Adaptive step method # +######################## + +function next(prob :: Problem{:adaptive}, state :: State) + + const timeout_const = 5 + + # the initial values + dt = state.dt # dt is the previous stepisze, it is + # modified inside the loop + timeout = state.timeout + + # for aesthetical reasons we extract the temporary componen + tmp = state.tmp + + # The while loop continues until we either find a stepsize which + # leads to a small enough error or the stepsize reaches + # prob.minstep + + while true + + # do one step (assumes ks[1]==f0), changes only tmp + err, newdt, timeout = rk_trial_step!(tmp, state, prob, dt, timeout) + + if abs(newdt) < prob.minstep # minimum step size reached, break + println("Warning: dt < minstep. Stopping.") + # passing the newdt to state will result in done() + state.dt = newdt + break + end + + if err > 1 # error is too large, repeat the step with smaller dt + # redo step with smaller dt and reset the timeout + dt = newdt + timeout = timeout_const + else + # step is accepted + + # preload ks[1] for the next step + if isFSAL(prob.method) + tmp.ks[1] = tmp.ks[S(prob.method)] + else + tmp.ks[1] = prob.F(state.t+dt, state.tmp.y) + end + + # Swap bindings of y and ytrial, avoids one copy + state.y, state.tmp.y = state.tmp.y, state.y + + # Update state with the data from the step we have just + # made: + state.t += dt + state.dt = newdt + state.timeout = timeout + break + end + end + return ((state.t,state.y),state) +end + + +########################## +# Lower level algorithms # +########################## + + +function rk_trial_step!(tmp, state, prob, dt, timeout) + + # tmp.y and tmp.yerr and tmp.ks are updated after this step + rk_embedded_step!(tmp, state, prob, dt) + + # changes tmp.yerr (via in place update) + err, newdt, timeout = stepsize_hw92!(tmp, state, prob, dt, timeout) + + return err, newdt, timeout +end + + +function rk_embedded_step!(tmp :: TempArrays, state :: State, prob :: Problem, dt) + # Does one embedded R-K step updating ytrial, yerr and ks. + # + # Assumes that ks[:,1] is already calculated! + # + # Modifies tmp.y and tmp.yerr only + + y = state.y + dof = length(y) + b = prob.method.b + + tmp.y[:] = 0 + tmp.yerr[:] = 0 + + for d=1:dof + tmp.y[d] += b[1,1]*tmp.ks[1][d] + tmp.yerr[d] += b[2,1]*tmp.ks[1][d] + end + + for s=2:S(prob.method) + calc_next_k!(state.tmp, s, state, prob) + for d=1:dof + tmp.y[d] += b[1,s]*tmp.ks[s][d] + tmp.yerr[d] += b[2,s]*tmp.ks[s][d] + end + end + + for d=1:dof + tmp.yerr[d] = dt * (tmp.y[d]-tmp.yerr[d]) + tmp.y[d] = y[d] + dt * tmp.y[d] + end +end + + +function stepsize_hw92!(tmp, state, prob, dt, timeout) + # Estimates the error and a new step size following Hairer & + # Wanner 1992, p167 (with some modifications) + # + # If timeout>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 + + order = minimum(prob.method.order) + 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. ? + dof = length(state.y) + + # in-place calculate yerr./tol + for d=1:dof + + # if outside of domain (usually NaN) then make step size smaller by maximum + if isoutofdomain(tmp.y[d]) + return 10., dt*facmin, timout_after_nan + end + + tmp.yerr[d] = tmp.yerr[d]/(prob.abstol + max(norm(prob.y0[d]), norm(tmp.y[d]))*prob.reltol) # Eq 4.10 + end + + err = norm(tmp.yerr, 2) # Eq. 4.11 + newdt = min(prob.maxstep, 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, newdt, timeout +end + + +function hinit(F, y0, t0, tstop, method, reltol, abstol) + # Returns first step size + order = minimum(method.order) + tau = max(reltol*norm(y0, Inf), abstol) + d0 = norm(y0, Inf)/tau + f0 = F(t0, y0) + d1 = norm(f0, Inf)/tau + if d0 < 1e-5 || d1 < 1e-5 + h0 = 1e-6 + else + h0 = 0.01*(d0/d1) + end + # perform Euler step + y1 = y0 + h0*f0 + f1 = F(t0 + h0, y1) + # estimate second derivative + d2 = norm(f1 - f0, Inf)/(tau*h0) + if max(d1, d2) <= 1e-15 + h1 = max(1e-6, 1e-3*h0) + else + pow = -(2 + log10(max(d1, d2)))/(order+1) + h1 = 10^pow + end + return min(100*h0, h1, tstop-t0) +end + + +# For clarity we pass the TempArrays part of the state separately, +# this is the only part of state that can be changed here +function calc_next_k!(tmp :: TempArrays, i, state :: State, prob :: Problem) + dof = length(state.y) + t, dt, a, c = state.t, state.dt, prob.method.a, prob.method.c + + tmp.y[:] = state.y + for j=1:i-1 + for d=1:dof + tmp.y[d] += dt * tmp.ks[j][d] * a[i,j] + end + end + tmp.ks[i] = prob.F(t + c[i]*dt, tmp.y) + + nothing +end diff --git a/src/runge_kutta.jl b/src/runge_kutta.jl deleted file mode 100644 index 4a3095cd5..000000000 --- a/src/runge_kutta.jl +++ /dev/null @@ -1,16 +0,0 @@ -# Explicit Runge-Kutta solvers -############################## -# (Hairer & Wanner 1992 p.134, p.165-169) - -import Base: start, next, done, call - -include("tableaus.jl") - -# include("algorithms.jl") - -include("iterators.jl") - -# include("iterators.jl") -# include("fixed.jl") -# include("variable.jl") -include("dense.jl")