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

WIP: Generalize oderkf to eventually replace ode23 #16

Merged
tomasaschan merged 7 commits into from
May 6, 2014
Merged

WIP: Generalize oderkf to eventually replace ode23 #16

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
af2c6a6
Generalize oderkf to support orders != 5. Add ode23_bs based on oderkf.
acroy Mar 2, 2014
be5ca6e
Add performance comparison for adaptive solvers.
acroy Mar 2, 2014
dcc2661
Some cosmetics & new test added (B5).
acroy Mar 4, 2014
4e1c417
Devectorization in oderkf, oderosenbrock and ode_ms (fixes problems r…
acroy Mar 12, 2014
f1f0041
Cleaned up performance test.
acroy Mar 13, 2014
c7dc86f
Added backward propagation in oderkf (part of #2).
acroy Mar 14, 2014
17a11e2
Remove performance test for now.
acroy May 5, 2014
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
96 changes: 68 additions & 28 deletions src/ODE.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -181,18 +181,19 @@ end # ode23
# CompereM@asme.org
# created : 06 October 1999
# modified: 17 January 2001
function oderkf(F, x0, tspan, a, b4, b5; 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)
# see p.91 in the Ascher & Petzold reference for more infomation.
pow = 1/5
pow = 1/p # use the higher order to estimate the next step size

c = sum(a, 2)
c = sum(a, 2) # consistency condition

# Initialization
t = tspan[1]
tfinal = tspan[end]
hmax = (tfinal - t)/2.5
hmin = (tfinal - t)/1e9
h = (tfinal - t)/100 # initial guess at a step size
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
x = x0
tout = t # first output time
xout = Array(typeof(x0), 1)
Expand All @@ -201,23 +202,29 @@ function oderkf(F, x0, tspan, a, b4, b5; reltol = 1.0e-5, abstol = 1.0e-8)
k = Array(typeof(x0), length(c))
k[1] = F(t,x) # first stage

while t < tfinal && h >= hmin
if t + h > tfinal
while abs(t) != abs(tfinal) && abs(h) >= hmin
if abs(h) > abs(tfinal-t)
h = tfinal - t
end

#(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)
k[j] = F(t + h.*c[j], x + h.*(a[j,1:j-1]*k[1:j-1])[1])
end

# compute the 4th order estimate
x4 = x + h.*(b4*k)[1]
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)

# compute the 5th order estimate
x5 = x + h.*(b5*k)[1]
# 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

# estimate the local truncation error
gamma1 = x5 - x4
gamma1 = xs - xp

# Estimate the error and the acceptable error
delta = norm(gamma1, Inf) # actual error
Expand All @@ -226,14 +233,14 @@ function oderkf(F, x0, tspan, a, b4, b5; reltol = 1.0e-5, abstol = 1.0e-8)
# Update the solution only if the error is acceptable
if delta <= tau
t = t + h
x = x5 # <-- using the higher order estimate is called 'local extrapolation'
x = xp # <-- using the higher order estimate is called 'local extrapolation'
tout = [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 && t != tspan[1]
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
Expand All @@ -248,20 +255,35 @@ function oderkf(F, x0, tspan, a, b4, b5; reltol = 1.0e-5, abstol = 1.0e-8)
h = min(hmax, 0.8*h*(tau/delta)^pow)
end # while (t < tfinal) & (h >= hmin)

if t < tfinal
println("Step size grew too small. t=", t, ", h=", h, ", x=", x)
if abs(t) < abs(tfinal)
println("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, tspan, x0) = oderkf(F, tspan, x0, bs_coefficients...)


# 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 = ([ 0 0 0 0 0 0
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
Expand All @@ -276,7 +298,8 @@ const dp_coefficients = ([ 0 0 0 0 0
ode45_dp(F, x0, tspan; kwargs...) = oderkf(F, x0, tspan, dp_coefficients...; kwargs...)

# Fehlberg coefficients
const fb_coefficients = ([ 0 0 0 0 0
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
Expand All @@ -291,7 +314,8 @@ ode45_fb(F, x0, tspan; kwargs...) = oderkf(F, x0, tspan, fb_coefficients...; kwa

# Cash-Karp coefficients
# Numerical Recipes in Fortran 77
const ck_coefficients = ([ 0 0 0 0 0
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
Expand All @@ -307,6 +331,10 @@ ode45_ck(F, x0, tspan; kwargs...) = oderkf(F, x0, tspan, ck_coefficients...; kwa
# Use Dormand Prince version of ode45 by default
const ode45 = ode45_dp

# some 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.
#
Expand Down Expand Up @@ -383,15 +411,23 @@ function oderosenbrock(F, x0, tspan, gamma, a, b, c; jacobian=nothing)
if size(dFdx,1) == 1
jac = 1/gamma/hs - dFdx[1]
else
jac = eye(dFdx)./gamma./hs - dFdx
jac = eye(dFdx)/gamma/hs - dFdx
end

g = Array(typeof(x0), size(a,1))
g[1] = (jac \ F(ts + b[1].*hs, xs))
g[1] = (jac \ F(ts + b[1]*hs, xs))
x[solstep+1] = x[solstep] + b[1]*g[1]

for i = 2:size(a,1)
g[i] = (jac \ (F(ts + b[i].*hs, xs + (a[i,1:i-1]*g[1:i-1])[1]) + (c[i,1:i-1]*g[1:i-1])[1]./hs))
dx = zero(x0)
dF = zero(x0/hs)
for j = 1:i-1
dx += a[i,j]*g[j]
dF += c[i,j]*g[j]
end
g[i] = (jac \ (F(ts + b[i]*hs, xs + dx) + dF/hs))
x[solstep+1] += b[i]*g[i]
end
x[solstep+1] = x[solstep] + (b*g)[1]
solstep += 1
end
return [tspan], x
Expand Down Expand Up @@ -457,7 +493,11 @@ function ode_ms(F, x0, tspan, order::Integer)
# 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] + (b[steporder,1:steporder]*xdot[i-(steporder-1):i])[1].*h[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 [tspan], x
end
Expand Down
4 changes: 3 additions & 1 deletion test/runtests.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,11 +5,13 @@ tol = 1e-2

solvers = [
ODE.ode23,
ODE.ode23_bs,

ODE.ode4,
ODE.ode45_dp,
ODE.ode45_fb,
ODE.ode45_ck,

ODE.ode4ms,
ODE.ode4s_s,
ODE.ode4s_kr]
Expand Down