From 3a200a430190f368a82f414fd866ad986bd75133 Mon Sep 17 00:00:00 2001 From: Andrei Berceanu Date: Wed, 25 Jun 2014 19:59:39 +0200 Subject: [PATCH 1/3] added min,max,init-step keywords and output at specified points --- src/ODE.jl | 37 +++++++++++++++++++++++++++++++------ 1 file changed, 31 insertions(+), 6 deletions(-) diff --git a/src/ODE.jl b/src/ODE.jl index 724ac1bb3..08355a7d7 100644 --- a/src/ODE.jl +++ b/src/ODE.jl @@ -122,6 +122,25 @@ function ode23(F, y0, tspan; reltol = 1.e-5, abstol = 1.e-8) end # ode23 +# helper functions +maxeps(x::FloatingPoint, y::FloatingPoint) = max(eps(abs(x)), eps(abs(y))) + +# `isapprox` with configurable tolerance +function approxeq(x::FloatingPoint, y::FloatingPoint; + rtol::FloatingPoint=cbrt(maxeps(x,y)), atol::FloatingPoint=sqrt(maxeps(x,y))) + abs(x-y) <= atol + rtol*max(abs(x), abs(y)) +end + +# an extension of the `in` statement for floating point values +function approxin(c::FloatingPoint, span::AbstractVector{Float64}; atol::FloatingPoint=.1) +#for some strange reason AbstractVector{FloatingPoint} does not work here! + truth = map(elem -> approxeq(c, elem; atol=atol), span) + for elem in truth + elem && return true + end + return false +end + # ode45 adapted from http://users.powernet.co.uk/kienzle/octave/matcompat/scripts/ode_v1.11/ode45.m @@ -181,7 +200,11 @@ end # ode23 # 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) +function oderkf(F, x0, tspan, p, a, bs, bp; reltol = 1.0e-5, abstol = 1.0e-8, + initstep = sign(tspan[end] - tspan[1])*abs(tspan[end] - tspan[1])/100, + minstep = abs(tspan[end] - tspan[1])/1e9, + maxstep = abs(tspan[end] - tspan[1])/2.5, + points = :all) # see p.91 in the Ascher & Petzold reference for more infomation. pow = 1/p # use the higher order to estimate the next step size @@ -191,9 +214,9 @@ function oderkf(F, x0, tspan, p, a, bs, bp; reltol = 1.0e-5, abstol = 1.0e-8) t = tspan[1] tfinal = tspan[end] tdir = sign(tfinal - t) - hmax = abs(tfinal - t)/2.5 - hmin = abs(tfinal - t)/1e9 - h = tdir*abs(tfinal - t)/100 # initial guess at a step size + hmax = maxstep + hmin = minstep + h = initstep x = x0 tout = t # first output time xout = Array(typeof(x0), 1) @@ -234,8 +257,10 @@ function oderkf(F, x0, tspan, p, a, bs, bp; reltol = 1.0e-5, abstol = 1.0e-8) if delta <= tau t = t + h x = xp # <-- using the higher order estimate is called 'local extrapolation' - tout = [tout; t] - push!(xout, x) + if points == :all || approxin(t, tspan; atol=.02) + tout = [tout; t] + push!(xout, x) + end # 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), From 8684b86c5a49535d297b375d4ab9daebefe93a6d Mon Sep 17 00:00:00 2001 From: Andrei Berceanu Date: Wed, 25 Jun 2014 20:35:07 +0200 Subject: [PATCH 2/3] using the standard Base.isapprox --- src/ODE.jl | 15 +++------------ 1 file changed, 3 insertions(+), 12 deletions(-) diff --git a/src/ODE.jl b/src/ODE.jl index 08355a7d7..5f4a65428 100644 --- a/src/ODE.jl +++ b/src/ODE.jl @@ -122,19 +122,11 @@ function ode23(F, y0, tspan; reltol = 1.e-5, abstol = 1.e-8) end # ode23 -# helper functions -maxeps(x::FloatingPoint, y::FloatingPoint) = max(eps(abs(x)), eps(abs(y))) - -# `isapprox` with configurable tolerance -function approxeq(x::FloatingPoint, y::FloatingPoint; - rtol::FloatingPoint=cbrt(maxeps(x,y)), atol::FloatingPoint=sqrt(maxeps(x,y))) - abs(x-y) <= atol + rtol*max(abs(x), abs(y)) -end +# helper functions # an extension of the `in` statement for floating point values -function approxin(c::FloatingPoint, span::AbstractVector{Float64}; atol::FloatingPoint=.1) -#for some strange reason AbstractVector{FloatingPoint} does not work here! - truth = map(elem -> approxeq(c, elem; atol=atol), span) +function approxin{T<:FloatingPoint}(c::FloatingPoint, span::AbstractVector{T}; atol::FloatingPoint=.1) + truth = map(elem -> isapprox(c, elem; atol=atol), span) for elem in truth elem && return true end @@ -142,7 +134,6 @@ function approxin(c::FloatingPoint, span::AbstractVector{Float64}; atol::Floatin end - # 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) # From 75485186f73d52359981aef7a27ed7bd94d2cb35 Mon Sep 17 00:00:00 2001 From: Andrei Berceanu Date: Fri, 27 Jun 2014 17:40:57 +0200 Subject: [PATCH 3/3] implemented algo for initial step selection --- src/ODE.jl | 32 ++++++++++++++++++++++++++++++-- 1 file changed, 30 insertions(+), 2 deletions(-) diff --git a/src/ODE.jl b/src/ODE.jl index 5f4a65428..95c1c777b 100644 --- a/src/ODE.jl +++ b/src/ODE.jl @@ -191,14 +191,42 @@ end # CompereM@asme.org # created : 06 October 1999 # modified: 17 January 2001 + +# estimator for initial step based on book +# "Solving Ordinary Differential Equations I" by Hairer et al., p.169 +function hinit(F, x0, t0, p, reltol, abstol) + tau = max(reltol*norm(x0, Inf), abstol) + d0 = norm(x0, Inf)/tau + f0 = F(t0, x0) + d1 = norm(f0, Inf)/tau + if d0 < 1e-5 || d1 < 1e-5 + h0 = 1e-6 + else + h0 = 1e-2d0/d1 + end + # perform Euler step + x1 = x0 + h0*f0 + f1 = F(t0 + h0, x1) + # estimate second derivative + d2 = norm(f1 - f0, Inf)/(tau*h0) + if max(d1, d2) <= 1e-15 + h1 = max(1e-6, 1e-3h0) + else + pow = -(2. + log10(max(d1, d2)))/(p+1.) + h1 = 10.^pow + end + h = min(100.0h0, h1) +end + + function oderkf(F, x0, tspan, p, a, bs, bp; reltol = 1.0e-5, abstol = 1.0e-8, - initstep = sign(tspan[end] - tspan[1])*abs(tspan[end] - tspan[1])/100, + initstep = hinit(F, x0, tspan[1], p, reltol, abstol), minstep = abs(tspan[end] - tspan[1])/1e9, maxstep = abs(tspan[end] - tspan[1])/2.5, points = :all) # see p.91 in the Ascher & Petzold reference for more infomation. pow = 1/p # use the higher order to estimate the next step size - + @show initstep c = sum(a, 2) # consistency condition # Initialization