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

[WIP] adding output at specified points #43

Closed
wants to merge 3 commits into from
+51 −7
Closed

[WIP] adding output at specified points #43

Changes from 1 commit
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
3a200a4
added min,max,init-step keywords and output at specified points
Jun 25, 2014
8684b86
using the standard Base.isapprox
Jun 25, 2014
7548518
implemented algo for initial step selection
Jun 27, 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
37 changes: 31 additions & 6 deletions src/ODE.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)))
Copy link
Contributor

Choose a reason for hiding this comment

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

Why do you have abs in eps(abs(...))? According to the documentation eps(x) is the difference between x and the next larger number.


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

Choose a reason for hiding this comment

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

I think you wanted the signature to look like

function approxin{T<:FloatingPoint}(c::FloatingPoint, span::AbstractVector{T}; atol::FloatingPoint=.1)

There are lots of discussion about this topic, but it boils down to AbstractVector{Float64} <: AbstractVector{FloatingPoint} not being true.

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
Expand Down Expand Up @@ -181,7 +200,11 @@ 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,
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,
Copy link
Contributor

Choose a reason for hiding this comment

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

It's a bit dangerous to use the difference between tspan[end] and tspan[1] here, since they don't really tell us anything about the problem - what if the user wants to integrate on [0,Inf] and break when a certain condition is met? What if the user wants to use this library for a really long-running task and doesn't care that the solver will take over 1e12 steps?

A better approach would be to base an estimate on some property of the ODE system itself, e.g. initstep = norm(F(t,x0)) / 100 for some well-behaved norm. Maybe we could even base it on stability properties of the method? There seems to be quite a lot of literature on the subject of stability regions for RK methods.

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

Expand All @@ -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)
Expand Down Expand Up @@ -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),
Expand Down