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

Add keywords minstep, maxstep and initstep to ode23s and oderkf #59

Merged
merged 3 commits into from
Mar 28, 2015
Merged

Add keywords minstep, maxstep and initstep to ode23s and oderkf #59

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
c04447a
Add keywords maxstep, minstep and initstep to ode23s and oderkf. New …
acroy Mar 24, 2015
42f0fd5
Fix initialization of array tout.
acroy Mar 24, 2015
a1e3d17
Fix indentation, remove hmax and hmin, make hinit return inital rhs. …
acroy Mar 24, 2015
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
103 changes: 78 additions & 25 deletions src/ODE.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,41 @@ export ode4s, ode4ms, ode4
# oderosenbrock, ode4s, ode4s_kr, ode4s_s,
# ode4ms, ode_ms

###############################################################################
## HELPER FUNCTIONS
###############################################################################

# 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
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The indentation in this if statement and the next makes it more difficult to understand the control flow here, IMO.

h0 = 1e-6
else
h0 = 0.01*(d0/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
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I haven't studied the algorithm, but is 1e-15 here really intended to be equivalent to eps()? If so, we might want to explicitly say eps(<some quantity>) or at least eps(typeof(<some quantity>)) to stay generic (what if the user wants to solve stuff with Float32, or BigFloat?).

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In the book they use those values, but don't say anything about the FP precision. But I guess your suggestion makes sense.

h1 = max(1e-6, 1e-3*h0)
else
pow = -(2. + log10(max(d1, d2)))/(p + 1.)
h1 = 10.^pow
end
h = min(100*h0, h1), f0
end

###############################################################################
## NON-STIFF SOLVERS
###############################################################################

# NOTE: in the near future ode23 will be replaced by ode23_bs

#ODE23 Solve non-stiff differential equations.
#
Expand Down Expand Up @@ -61,14 +96,12 @@ function ode23(F, y0, tspan; reltol = 1.e-5, abstol = 1.e-8)
hmax = abs(0.1*(tfinal-t))
y = y0

tout = [t]
tout = Array(typeof(t), 1)
tout[1] = t # first output time
yout = Array(typeof(y0),1)
yout[1] = y

tlen = length(t)
yout[1] = y # first output solution

# Compute initial step size.

s1 = F(t, y)
r = norm(s1./max(abs(y), threshold), Inf) + realmin() # TODO: fix type bug in max()
h = tdir*0.8*reltol^(1/3)/r
Expand Down Expand Up @@ -186,28 +219,37 @@ end # ode23
# [email protected]
# 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,
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)
hmax = abs(tfinal - t)/2.5
hmin = abs(tfinal - t)/1e9
h = tdir*abs(tfinal - t)/100 # initial guess at a step size

h = initstep
if h == 0.
# initial guess at a step size
h, k[1] = hinit(F, x0, t, p, reltol, abstol)
else
k[1] = F(t, x0) # first stage
end
h = tdir*min(h, maxstep)
x = x0
tout = [t] # first output time
tout = Array(typeof(t), 1)
tout[1] = t # first output time
xout = Array(typeof(x0), 1)
xout[1] = x # first output solution

k = Array(typeof(x0), length(c))
k[1] = F(t,x) # first stage

while abs(t) != abs(tfinal) && abs(h) >= hmin
while abs(t) != abs(tfinal) && abs(h) >= minstep
if abs(h) > abs(tfinal-t)
h = tfinal - t
end
Expand Down Expand Up @@ -257,8 +299,8 @@ function oderkf(F, x0, tspan, p, a, bs, bp; reltol = 1.0e-5, abstol = 1.0e-8)
end

# Update the step size
h = min(hmax, 0.8*h*(tau/delta)^pow)
end # while (t < tfinal) & (h >= hmin)
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")
Expand All @@ -267,6 +309,7 @@ function oderkf(F, x0, tspan, p, a, bs, bp; reltol = 1.0e-5, abstol = 1.0e-8)
return tout, xout
end


# Bogacki–Shampine coefficients
const bs_coefficients = (3,
[ 0 0 0 0
Expand Down Expand Up @@ -302,6 +345,7 @@ const dp_coefficients = (5,
)
ode45_dp(F, x0, tspan; kwargs...) = oderkf(F, x0, tspan, dp_coefficients...; kwargs...)


# Fehlberg coefficients
const fb_coefficients = (5,
[ 0 0 0 0 0
Expand All @@ -317,6 +361,7 @@ const fb_coefficients = (5,
)
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,
Expand Down Expand Up @@ -364,7 +409,7 @@ const ode78 = ode78_fb
# Use Dormand Prince version of ode45 by default
const ode45 = ode45_dp

# some higher-order embedded methods can be found in:
# 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.


Expand Down Expand Up @@ -427,7 +472,7 @@ function fdjacobian(F, x, t)
end

# ODE23S Solve stiff systems based on a modified Rosenbrock triple
# (also used by MATLABS ODE23s); see Sec. 4.1 in
# (also used by MATLAB's ODE23s); see Sec. 4.1 in
#
# [SR97] L.F. Shampine and M.W. Reichelt: "The MATLAB ODE Suite," SIAM Journal on Scientific Computing, Vol. 18, 1997, pp. 1–22
#
Expand All @@ -436,7 +481,11 @@ end
function ode23s(F, y0, tspan; reltol = 1.0e-5, abstol = 1.0e-8,
jacobian=nothing,
points=:all,
norm=Base.norm)
norm=Base.norm,
minstep=abs(tspan[end] - tspan[1])/1e18,
maxstep=abs(tspan[end] - tspan[1])/2.5,
initstep=0.)


# select method for computing the Jacobian
if typeof(jacobian) == Function
Expand All @@ -456,21 +505,25 @@ function ode23s(F, y0, tspan; reltol = 1.0e-5, abstol = 1.0e-8,
tfinal = tspan[end]
tdir = sign(tfinal - t)

hmax = abs(tfinal - t)/10
hmin = abs(tfinal - t)/1e18
h = tdir*abs(tfinal - t)/100 # initial step size
h = initstep
if h == 0.
# initial guess at a step size
h, F0 = hinit(F, y0, t, 3, reltol, abstol)
else
F0 = F(t,y0)
end
h = tdir*min(h, maxstep)

y = y0
tout = Array(typeof(t), 1)
tout[1] = t # first output time
yout = Array(typeof(y0), 1)
yout[1] = copy(y) # first output solution

F0 = F(t,y)

J = jac(t,y) # get Jacobian of F wrt y

while abs(t) < abs(tfinal) && hmin < abs(h)
while abs(t) < abs(tfinal) && minstep < abs(h)
if abs(t-tfinal) < abs(h)
h = tfinal - t
end
Expand Down Expand Up @@ -531,7 +584,7 @@ function ode23s(F, y0, tspan; reltol = 1.0e-5, abstol = 1.0e-8,
end

# update of the step size
h = tdir*min( hmax, abs(h)*0.8*(delta/err)^(1/3) )
h = tdir*min( maxstep, abs(h)*0.8*(delta/err)^(1/3) )
end

return tout, yout
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
Expand Up @@ -63,7 +63,8 @@ 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-8, reltol=1e-8,
maxstep=1e11/10, minstep=1e11/1e18)

refsol = [0.2083340149701255e-07,
0.8333360770334713e-13,
Expand Down