diff --git a/LICENSE.md b/LICENSE.md index 0ef0e5e59..c8c001a77 100644 --- a/LICENSE.md +++ b/LICENSE.md @@ -1,18 +1,7 @@ -ode23 -===== +oderkf +====== -The ode23 function is adapted from Chapter 7 of the book Numerical Computing -with Matlab by Cleve Moler. The original source is at: -http://www.mathworks.in/moler/ncm/ode23tx.m - -The code carries the following notice: - -> Copyright 2012 Cleve Moler and The MathWorks, Inc. - -ode45 -===== - -The ode45 function is adapted from this implementation for Octave, and +The oderkf function is adapted from this implementation for Octave, and is available under the GPL. http://users.powernet.co.uk/kienzle/octave/matcompat/scripts/ode_v1.11/ode45.m @@ -35,3 +24,9 @@ ode4, ode4s, ode4ms =================== These are implemented by Patrick O'Leary, and available under the MIT License. + +ode23s +====== + +The ode23s functions is implemented by Alexander Croy, and available under the +MIT License. diff --git a/README.md b/README.md index 24750b2ae..ce9c06875 100644 --- a/README.md +++ b/README.md @@ -16,16 +16,26 @@ Currently, `ODE` exports the following adaptive solvers: * `ode23`: 2nd order adaptive solver with 3rd order error control, using the Bogacki–Shampine coefficients * `ode45`: 4th order adaptive solver with 5th order error control, using the Dormand Prince coefficients. Fehlberg and Cash-Karp coefficients are also available. -* `ode78`: 7th order adaptive solver with 8th order error control, using the Fehlberg coefficients +* `ode78`: 7th order adaptive solver with 8th order error control, using the Fehlberg coefficients. -* `ode23s`: 2nd/3rd order adaptive solver for stiff problems, using a modified Rosenbrock triple +* `ode23s`: 2nd/3rd order adaptive solver for stiff problems, using a modified Rosenbrock triple. all of which have the following basic API: - tout, yout = odeXX(F, y0, tspan) + tout, yout = odeXX(F, y0, tspan; keywords...) to solve the explicit ODE defined by dy/dt = F(t,y). A few other solvers are also exported, see the source code for details. +The adaptive solvers accept the following keywords +- `norm`: user-supplied norm for determining the error `E` (default `Base.vecnorm`), +- `abstol` and/or `reltol`: an integration step is accepted if `E <= abstol || E <= reltol*abs(y)` (defaults `reltol = 1e-5`, `abstol = 1e-8`), +- `maxstep`, `minstep` and `initstep`: determine the maximal, minimal and initial integration step (defaults `minstep=|tspan[end] - tspan[1]|/1e9`, `maxstep=|tspan[end] - tspan[1]|/2.5` and automatic initial step estimation). + +Additionally, `ode23s` supports +- `points=:all` (default): output is given for each value in `tspan` as well as for each intermediate point the solver used. +- `points=:specified`: output is given only for each value in `tspan`. +- `jacobian = G(t,y)`: user-supplied Jacobian G(t,y) = dF(t,y)/dy (default estimate by finite-difference method). + # Need something long-term reliable right now? See [the Sundials.jl package](https://github.com/julialang/sundials.jl), which provides wrappers for the excellent Sundials ODE solver library. diff --git a/src/ODE.jl b/src/ODE.jl index 52a47323a..32c884175 100644 --- a/src/ODE.jl +++ b/src/ODE.jl @@ -52,173 +52,17 @@ end ## NON-STIFF SOLVERS ############################################################################### -# NOTE: in the near future ode23 will be replaced by ode23_bs - -#ODE23 Solve non-stiff differential equations. -# -# ODE23(F,TSPAN,Y0) with TSPAN = [T0 TFINAL] integrates the system -# of differential equations dy/dt = f(t,y) from t = T0 to t = TFINAL. -# The initial condition is y(T0) = Y0. -# -# The first argument, F, is a function handle or an anonymous function -# that defines f(t,y). This function must have two input arguments, -# t and y, and must return a column vector of the derivatives, dy/dt. -# -# With two output arguments, [T,Y] = ODE23(...) returns a column -# vector T and an array Y where Y(:,k) is the solution at T(k). +# ODERKF based on # -# More than four input arguments, ODE23(F,TSPAN,Y0,RTOL,P1,P2,...), -# are passed on to F, F(T,Y,P1,P2,...). -# -# ODE23 uses the Runge-Kutta (2,3) method of Bogacki and Shampine (BS23). -# -# Example -# tspan = [0, 2*pi] -# y0 = [1, 0] -# F = (t, y) -> [0 1; -1 0]*y -# ode23(F, tspan, y0) -# -# See also ODE23. - -# Initialize variables. -# Adapted from Cleve Moler's textbook -# http://www.mathworks.com/moler/ncm/ode23tx.m -function ode23(F, y0, tspan; reltol = 1.e-5, abstol = 1.e-8) - if reltol == 0 - warn("setting reltol = 0 gives a step size of zero") - end - - threshold = abstol / reltol - - t = tspan[1] - tfinal = tspan[end] - tdir = sign(tfinal - t) - hmax = abs(0.1*(tfinal-t)) - y = y0 - - tout = Array(typeof(t), 1) - tout[1] = t # first output time - yout = Array(typeof(y0),1) - 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 - - # The main loop. - - while t != tfinal - - hmin = 16*eps()*abs(t) - if abs(h) > hmax; h = tdir*hmax; end - if abs(h) < hmin; h = tdir*hmin; end - - # Stretch the step if t is close to tfinal. - - if 1.1*abs(h) >= abs(tfinal - t) - h = tfinal - t - end - - # Attempt a step. - - s2 = F(t+h/2, y+h/2*s1) - s3 = F(t+3*h/4, y+3*h/4*s2) - tnew = t + h - ynew = y + h*(2*s1 + 3*s2 + 4*s3)/9 - s4 = F(tnew, ynew) - - # Estimate the error. - - e = h*(-5*s1 + 6*s2 + 8*s3 - 9*s4)/72 - err = norm(e./max(max(abs(y), abs(ynew)), threshold), Inf) + realmin() - - # Accept the solution if the estimated error is less than the tolerance. - - if err <= reltol - t = tnew - y = ynew - push!(tout, t) - push!(yout, y) - s1 = s4 # Reuse final function value to start new step - end - - # Compute a new step size. - - h = h*min(5, 0.8*(reltol/err)^(1/3)) - - # Exit early if step size is too small. - - if abs(h) <= hmin - println("Step size ", h, " too small at t = ", t) - t = tfinal - end - - end # while (t != tfinal) - - return tout, yout - -end # ode23 - - - # 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) # -# ode45 (v1.11) integrates a system of ordinary differential equations using -# 4th & 5th order embedded formulas from Dormand & Prince or Fehlberg. -# -# The Fehlberg 4(5) pair is established and works well, however, the -# Dormand-Prince 4(5) pair minimizes the local truncation error in the -# 5th-order estimate which is what is used to step forward (local extrapolation.) -# Generally it produces more accurate results and costs roughly the same -# computationally. The Dormand-Prince pair is the default. -# -# This is a 4th-order accurate integrator therefore the local error normally -# expected would be O(h^5). However, because this particular implementation -# uses the 5th-order estimate for xout (i.e. local extrapolation) moving -# forward with the 5th-order estimate should yield errors on the order of O(h^6). -# -# The order of the RK method is the order of the local *truncation* error, d. -# The local truncation error is defined as the principle error term in the -# portion of the Taylor series expansion that gets dropped. This portion of -# the Taylor series exapansion is within the group of terms that gets multipled -# by h in the solution definition of the general RK method. Therefore, the -# order-p solution created by the RK method will be roughly accurate to -# within O(h^(p+1)). The difference between two different-order solutions is -# the definition of the "local error," l. This makes the local error, l, as -# large as the error in the lower order method, which is the truncation -# error, d, times h, resulting in O(h^(p+1)). -# Summary: For an RK method of order-p, -# - the local truncation error is O(h^p) -# - the local error (used for stepsize adjustment) is O(h^(p+1)) -# -# This requires 6 function evaluations per integration step. -# -# The error estimate formula and slopes are from -# Numerical Methods for Engineers, 2nd Ed., Chappra & Cannle, McGraw-Hill, 1985 -# -# Usage: -# (tout, xout) = ode45(F, tspan, x0) -# -# INPUT: -# F - User-defined function -# Call: xprime = F(t,x) -# t - Time (scalar). -# x - Solution column-vector. -# xprime - Returned derivative COLUMN-vector; xprime(i) = dx(i)/dt. -# tspan - [ tstart, tfinal ] -# x0 - Initial value COLUMN-vector. -# -# OUTPUT: -# tout - Returned integration time points (column-vector). -# xout - Returned solution, one solution column-vector per tout-value. -# # Original Octave implementation: # Marc Compere # 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, norm=Base.norm, minstep=abs(tspan[end] - tspan[1])/1e9, @@ -406,9 +250,13 @@ ode78_fb(F, x0, tspan; kwargs...) = oderkf(F, x0, tspan, fb_coefficients_78...; # Use Fehlberg version of ode78 by default const ode78 = ode78_fb -# Use Dormand Prince version of ode45 by default +# Use Dormand-Prince version of ode45 by default const ode45 = ode45_dp +# Use Bogacki–Shampine version of ode23 by default +const ode23 = ode23_bs + + # 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. @@ -508,7 +356,7 @@ function ode23s(F, y0, tspan; reltol = 1.0e-5, abstol = 1.0e-8, h = initstep if h == 0. # initial guess at a step size - h, F0 = hinit(F, y0, t, 3, reltol, abstol) + h, F0 = hinit(F, y0, t, 3, reltol, abstol) else F0 = F(t,y0) end diff --git a/test/runtests.jl b/test/runtests.jl index cb1bc2f75..ef743b88e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -7,7 +7,6 @@ ode5ms(F, x0, tspan) = ODE.ode_ms(F, x0, tspan, 5) solvers = [ ODE.ode23, - ODE.ode23_bs, ODE.ode4, ODE.ode45_dp,