ode.jl

#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).
#
#   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]
#      y_0 = [1 0]'
#      F = (t, y) -> [0 1; -1 0]*y
#      ode23(F, tspan, y_0)
#
#   See also ODE23.

# Initialize variables.
# Adapted from Cleve Moler's textbook
# http://www.mathworks.com/moler/ncm/ode23tx.m

function ode23(F::Function, tspan::Vector, y_0)

    rtol = 1.e-5
    atol = 1.e-8

    t0 = tspan[1]
    tfinal = tspan[2]
    tdir = sign(tfinal - t0)
    threshold = atol / rtol
    hmax = abs(0.1*(tfinal-t0))
    t = t0
    y = y_0[:]

    tout = t
    yout = y.'

    tlen = length(t)

    # Compute initial step size.

    s1 = F(t, y)
    r = norm(s1./max(abs(y), threshold), Inf) + realmin()
    h = tdir*0.8*rtol^(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 <= rtol
            t = tnew
            y = ynew
            tout = [tout; t]
            yout = [yout; y.']
            s1 = s4   # Reuse final function value to start new step
        end
            
        # Compute a new step size.
        
        h = h*min(5, 0.8*(rtol/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

# 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.
#
# Both the Dormand-Prince and Fehlberg 4(5) coefficients are from a tableu 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
#
# 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.
#
# Marc Compere
# CompereM@asme.org
# created : 06 October 1999
# modified: 17 January 2001

# Dormand Prince
function ode45_dp(F, tspan, x0)
    tol = 1.0e-5
    
    # see p.91 in the Ascher & Petzold reference for more infomation.
    pow = 1/6; 

    a_ = zeros(7,7)
    b4_ = zeros(7)
    b5_ = zeros(7)
    c_ = zeros(7)

    a_[1,1]=0;
    a_[2,1]=1/5;
    a_[3,1]=3/40;       a_[3,2]=9/40;
    a_[4,1]=44/45;      a_[4,2]=-56/15;      a_[4,3]=32/9;
    a_[5,1]=19372/6561; a_[5,2]=-25360/2187; a_[5,3]=64448/6561; a_[5,4]=-212/729;
    a_[6,1]=9017/3168;  a_[6,2]=-355/33;     a_[6,3]=46732/5247; a_[6,4]=49/176;   a_[6,5]=-5103/18656;
    a_[7,1]=35/384;     a_[7,2]=0;           a_[7,3]=500/1113;   a_[7,4]=125/192;  a_[7,5]=-2187/6784;  a_[7,6]=11/84;
    # 4th order b-coefficients
    b4_[1]=5179/57600; b4_[2]=0; b4_[3]=7571/16695; b4_[4]=393/640; b4_[5]=-92097/339200; b4_[6]=187/2100; b4_[7]=1/40;
    # 5th order b-coefficients
    b5_[1]=35/384; b5_[2]=0; b5_[3]=500/1113; b5_[4]=125/192; b5_[5]=-2187/6784; b5_[6]=11/84; b5_[7]=0;
    for i=1:7
        c_[i]=sum(a_[i,:])
    end

    # Initialization
    t0 = tspan[1]
    tfinal = tspan[2]
    t = t0
    hmax = (tfinal - t)/2.5
    hmin = (tfinal - t)/1e9
    h = (tfinal - t)/100  # initial guess at a step size
    x = x0[:]             # this always creates a column vector, x
    tout = t              # first output time
    xout = x.'            # first output solution
        
    # k_ needs to be initialized as an Nx7 matrix where N=number of rows in x
    # (just for speed so octave doesn't need to allocate more memory at each stage value)
    k_ = zeros(length(x),7)
    
    # Compute the first stage prior to the main loop.  This is part of the Dormand-Prince pair caveat
    # Normally, during the main loop the last stage for x[k] is the first stage for computing x[k+1].
    # So, the very first integration step requires 7 function evaluations, then each subsequent step
    # 6 function evaluations because the first stage is simply assigned from the last step's last stage.
    # note: you can see this by the last element in c_ is 1.0, thus t+c_[7]*h = t+h, ergo, the next step.
    
    k_[:,1] = F(t,x) # first stage
    
    # The main loop using Dormand-Prince 4(5) pair
    while (t < tfinal) & (h >= hmin)
        if t + h > tfinal; h = tfinal - t; 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),
        #        s is the number of intermediate RK stages on [t (t+h)] (Dormand-Prince has s=7 stages)

        for j = 1:6
            k_[:,j+1] = F(t+c_[j+1]*h, x + reshape(h*k_[:,1:j]*a_[j+1,1:j]', length(x)))
        end
                
	# compute the 4th order estimate
        x4=x + h* (k_*b4_) # k_ is Nxs (or Nx7) and b4_ is a 7x1
                    
        # compute the 5th order estimate
        x5=x + h*(k_*b5_) # k_ is Nxs (or Nx7) and b5_ is a 7x1
                
        # estimate the local truncation error
        gamma1 = x5 - x4
                
        # Estimate the error and the acceptable error
        delta = norm(gamma1, Inf)       # actual error
        tau = tol*max(norm(x,Inf), 1.0) # allowable error
                    
	# 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'
            tout = [tout; t]
            xout = [xout; x.']
        end

	# Update the step size
        if (delta==0.0); delta = 1e-16; end
        h = min(hmax, 0.8*h*(tau/delta)^pow)
                    
        # 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
	# again as k_[:,1] during the next step.
        k_[:,1]=k_[:,7]

    end # while (t < tfinal) & (h >= hmin)

    if (t < tfinal)
      println("Step size grew too small. t=", t, ", h=", h, ", x=", x)
    end

    return (tout, xout)
    
end # ode45_dp

# Fehlberg
function ode45_fb(F, tspan, x0)
    tol = 1.0e-5

    # see p.91 in the Ascher & Petzold reference for more infomation.
    pow = 1/6; 

    a_ = zeros(6,6)
    b4_ = zeros(5)
    b5_ = zeros(6)
    c_ = zeros(6)

    a_[1,1]=0;
    a_[2,1]=1/4;
    a_[3,1]=3/32;      a_[3,2]=9/32;
    a_[4,1]=1932/2197; a_[4,2]=-7200/2197; a_[4,3]=7296/2197;
    a_[5,1]=439/216;   a_[5,2]=-8;         a_[5,3]=3680/513;   a_[5,4]=-845/4104;
    a_[6,1]=-8/27;     a_[6,2]=2;          a_[6,3]=-3544/2565; a_[6,4]=1859/4104; a_[6,5]=-11/40;
    # 4th order b-coefficients (guaranteed to be a column vector)
    b4_[1]=25/216; b4_[2]=0; b4_[3]=1408/2565; b4_[4]=2197/4104; b4_[5]=-1/5;
    # 5th order b-coefficients (also guaranteed to be a column vector)
    b5_[1]=16/135; b5_[2]=0; b5_[3]=6656/12825; b5_[4]=28561/56430; b5_[5]=-9/50; b5_[6]=2/55;
    for i=1:6
        c_[i]=sum(a_[i,:])
    end

    # Initialization
    t0 = tspan[1]
    tfinal = tspan[2]
    t = t0
    hmax = (tfinal - t)/2.5
    hmin = (tfinal - t)/1e9
    h = (tfinal - t)/100  # initial guess at a step size
    x = x0[:]             # this always creates a column vector, x
    tout = t              # first output time
    xout = x.'            # first output solution
        
    # k_ needs to be initialized as an Nx6 matrix where N=number of rows in x
    # (just for speed so octave doesn't need to allocate more memory at each stage value)
    k_ = zeros(length(x),6)
    
    while (t < tfinal) & (h >= hmin)
        if t + h > tfinal; h = tfinal - t; 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 Nx6, a_ is 6x5, which is s by (s-1),  (RK-Fehlberg has s=6 stages)
        #        s is the number of intermediate RK stages on [t (t+h)]
        
        k_[:,1]=F(t,x) # first stage
        
        for j = 1:5
            k_[:,j+1] = F(t+c_[j+1]*h, x + reshape(h*k_[:,1:j]*a_[j+1,1:j]', length(x)))
        end
        
        # compute the 4th order estimate
        x4=x + h* (k_[:,1:5]*b4_) # k_[:,1:5] is an Nx5 and b4_ is a 5x1
        
	# compute the 5th order estimate
        x5=x + h*(k_*b5_) # k_ is the same as k_[:,1:6] and is an Nx6 and b5_ is a 6x1
        
        # estimate the local truncation error
        gamma1 = x5 - x4
        
        # Estimate the error and the acceptable error
        delta = norm(gamma1, Inf)       # actual error
        tau = tol*max(norm(x, Inf),1.0) # allowable error
        
	# 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'
            tout = [tout; t]
            xout = [xout; x.']
        end 

        # Update the step size
        if (delta==0.0); delta = 1e-16; end
        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)
    end

    return (tout, xout)
    
end # ode45_fb

# Use Dormand Prince version of ode45 by default
const ode45 = ode45_dp