diff --git a/src/ODE.jl b/src/ODE.jl index e8ccc55c1..93f5745fb 100644 --- a/src/ODE.jl +++ b/src/ODE.jl @@ -62,14 +62,14 @@ Base.convert{Tnew<:Real}(::Type{Tnew}, tab::Tableau) = error("Define convert met # estimator for initial step based on book # "Solving Ordinary Differential Equations I" by Hairer et al., p.169 -function hinit(F, x0, t0, tend, p, reltol, abstol) +function hinit(F, x0, t0, tend, p, reltol, abstol, norm) # Returns first step, direction of integration and F evaluated at t0 tdir = sign(tend-t0) tdir==0 && error("Zero time span") - tau = max(reltol*norm(x0, Inf), abstol) - d0 = norm(x0, Inf)/tau + tau = max(reltol.*abs(x0), abstol) + d0 = norm(x0./tau) f0 = F(t0, x0) - d1 = norm(f0, Inf)/tau + d1 = norm(f0./tau) if d0 < 1e-5 || d1 < 1e-5 h0 = 1e-6 else @@ -79,7 +79,7 @@ function hinit(F, x0, t0, tend, p, reltol, abstol) x1 = x0 + tdir*h0*f0 f1 = F(t0 + tdir*h0, x1) # estimate second derivative - d2 = norm(f1 - f0, Inf)/(tau*h0) + d2 = norm((f1 - f0)./tau)/h0 if max(d1, d2) <= 1e-15 h1 = max(1e-6, 1e-3*h0) else @@ -220,6 +220,40 @@ function fdjacobian(F, x, t) return dFdx end +## Function to get the Scaled error given the error estimate +# The following 3 functions are for the stiff solver +# If both reltol and abstol are scalars then restore previous behaviour +function getScaledError(y,ynew,errEstimate,h,reltol::Number,abstol::Number,norm) + return (abs(h)/6)*norm(errEstimate)/max(reltol*max(norm(y),norm(ynew)), abstol) # scaled error estimate +end +# If atleast one of them is a vector then perform component-wise computations +function getScaledError(y,ynew,errEstimate,h,reltol::Number,abstol::Vector,norm) + return (abs(h)/6)*norm((errEstimate)./max(reltol.*max(abs(y),abs(ynew)),abstol)) # scaled error estimate +end +function getScaledError(y,ynew,errEstimate,h,reltol::Vector,abstol,norm) + return (abs(h)/6)*norm((errEstimate)./max(reltol.*max(abs(y),abs(ynew)),abstol)) # scaled error estimate +end +# The following 2 functions are for the non-stiff solvers since it requires dof check +# As per the old code, the error estimate is rewritten with the scaled error +# Estimates the error and a new step size following Hairer & +# Wanner 1992, p167 (with some modifications) +function getScaledError!(y,ynew,errEstimate,reltol::Number,abstol::Number,norm,dof) + for d=1:dof + # if outside of domain (usually NaN) then make step size smaller by maximum + isoutofdomain(ynew[d]) && return true + errEstimate[d] = errEstimate[d]/(abstol + max(abs(y[d]), abs(ynew[d]))*reltol) # Eq 4.10 + end + return false +end +function getScaledError!(y,ynew,errEstimate,reltol::Vector,abstol::Vector,norm,dof) + for d=1:dof + # if outside of domain (usually NaN) then break and return + isoutofdomain(ynew[d]) && return true + errEstimate[d] = errEstimate[d]/(abstol[d] + max(abs(y[d]), abs(ynew[d]))*reltol[d]) # Eq 4.10 + end + return false +end + # ODE23S Solve stiff systems based on a modified Rosenbrock triple # (also used by MATLAB's ODE23s); see Sec. 4.1 in # @@ -251,13 +285,17 @@ function ode23s(F, y0, tspan; reltol = 1.0e-5, abstol = 1.0e-8, # initialization t = tspan[1] + + # Component-wise Tolerance check similar to the non-stiff solvers + @assert length(abstol) == 1 || length(abstol) == length(y0) "Dimension of Absolute tolerance does not match the dimension of the problem" + @assert length(reltol) == 1 || length(reltol) == length(y0) "Dimension of Absolute tolerance does not match the dimension of the problem" tfinal = tspan[end] h = initstep if h == 0. # initial guess at a step size - h, tdir, F0 = hinit(F, y0, t, tfinal, 3, reltol, abstol) + h, tdir, F0 = hinit(F, y0, t, tfinal, 3, reltol, abstol, norm) else tdir = sign(tfinal - t) F0 = F(t,y0) @@ -300,11 +338,19 @@ function ode23s(F, y0, tspan; reltol = 1.0e-5, abstol = 1.0e-8, F2 = F(t + h, ynew) k3 = W\(F2 - e32*(k2 - F1) - 2*(k1 - F0) + T ) - err = (abs(h)/6)*norm(k1 - 2*k2 + k3) # error estimate - delta = max(reltol*max(norm(y),norm(ynew)), abstol) # allowable error + # If reltol and abstol are vectors + # restore old behvaiour + # else perform component-wise computations for error + errEstimate = k1 - 2*k2 + k3; + err = getScaledError(y,ynew,errEstimate,h,reltol,abstol,norm) + #if length(abstol) == 1 && length(reltol) == 1 + # err = (abs(h)/6)*norm(k1 - 2*k2 + k3)/max(reltol*max(norm(y),norm(ynew)), abstol) # scaled error estimate + #else + # err = (abs(h)/6)*norm((k1 - 2*k2 + k3)./max(reltol.*max(abs(y),abs(ynew)),abstol)) # scaled error estimate + #end # check if new solution is acceptable - if err <= delta + if err <= 1 if points==:specified || points==:all # only points in tspan are requested @@ -334,7 +380,7 @@ function ode23s(F, y0, tspan; reltol = 1.0e-5, abstol = 1.0e-8, end # update of the step size - h = tdir*min( maxstep, abs(h)*0.8*(delta/err)^(1/3) ) + h = tdir*min( maxstep, abs(h)*0.8*(1/err)^(1/3) ) end return tout, yout diff --git a/src/runge_kutta.jl b/src/runge_kutta.jl index a080602de..e6ffe2437 100644 --- a/src/runge_kutta.jl +++ b/src/runge_kutta.jl @@ -250,6 +250,10 @@ function oderk_adapt{N,S}(fn, y0::AbstractVector, tspan, btab_::TableauRKExplici t = tspan[1] tstart = tspan[1] tend = tspan[end] + + # Component-wise reltol and abstol as per Solving Ordinary Differential Equations I by Hairer and Wanner + @assert length(abstol) == 1 || length(abstol) == length(y0) "Dimension of Absolute tolerance does not match the dimension of the problem" + @assert length(reltol) == 1 || length(reltol) == length(y0) "Dimension of Relative tolerance does not match the dimension of the problem" # work arrays: y = similar(y0, Eyf, dof) # y at time t @@ -276,7 +280,7 @@ function oderk_adapt{N,S}(fn, y0::AbstractVector, tspan, btab_::TableauRKExplici error("Unrecognized option points==$points") end # Time - dt, tdir, ks[1] = hinit(fn, y, tstart, tend, order, reltol, abstol) # sets ks[1]=f0 + dt, tdir, ks[1] = hinit(fn, y, tstart, tend, order, reltol, abstol, norm) # sets ks[1]=f0 if initstep!=0 dt = sign(initstep)==tdir ? initstep : error("initstep has wrong sign.") end @@ -409,11 +413,9 @@ function stepsize_hw92!(dt, tdir, x0, xtrial, xerr, order, 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 + # Get the scaled error and if outside of domain (usually NaN) then make step size smaller by maximum + getScaledError!(x0,xtrial,xerr,reltol,abstol,norm,dof) && return 10., dt*facmin, timout_after_nan + 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 diff --git a/test/interface-tests.jl b/test/interface-tests.jl index dbb811dcc..816d6773c 100644 --- a/test/interface-tests.jl +++ b/test/interface-tests.jl @@ -45,19 +45,42 @@ ODE.isoutofdomain(y::CompSol) = any(isnan, vcat(y.rho[:], y.x, y.p)) ################################################################################ - +# TODO: test a vector-like state variable, i.e. one which can be indexed. +function getZeroVector(y::Array{CompSol,1}) + numElems = length(y) + tmp = Array{CompSol,1}(numElems) + for i = 1:numElems + tmp[i] = CompSol(complex(zeros(2,2)), 0., 0.) + end + return tmp +end +Base.zeros(y::Array{CompSol,1}) = getZeroVector(y) + +function getAbsVector(y::Array{CompSol,1}) + numElems = length(y) + absVec = Array{Float64}(numElems) + for i = 1:numElems + absVec[i] = norm(y[i], 2.) + end + return absVec +end +Base.abs(y::Array{CompSol,1}) = getAbsVector(y) + + +################################################################################ + # 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.; @@ -77,6 +100,3 @@ for solver in solvers @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 9c88cd281..141390dad 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -75,7 +75,7 @@ let ydot end t = [0., 1e11] - t,y = ode23s(f, [1.0, 0.0, 0.0], t; abstol=1e-8, reltol=1e-8, + t,y = ode23s(f, [1.0, 0.0, 0.0], t; abstol=1e-9*ones(3), reltol=1e-9*ones(3), maxstep=1e11/10, minstep=1e11/1e18) refsol = [0.2083340149701255e-07, @@ -83,6 +83,7 @@ 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")