Addition of Variable and Fixed Adam Moulton Method

We add two new solvers for ODE.jl, and relocate the old multistep solver
ode_ms, namely ode_am and ode113. ode_am is a fixed step method similiar
to ode_ms, but uses a PECE step. ode113 is a variable step size and
variable order method using the same underlying theory as ode_am.For more
documentation, refer to Hairer et al Volume 1 on Adam Methods.

src/ODE.jl

@@ -8,9 +8,9 @@ using Compat
 
 ## minimal function export list
 # adaptive non-stiff:
-export ode23, ode45, ode78
+export ode23, ode45, ode78, ode113
 # non-adaptive non-stiff:
-export ode4, ode4ms
+export ode4, ode4ms, ode4am
 # adaptive stiff:
 export ode23s
 # non-adaptive stiff:
@@ -93,7 +93,7 @@ end
 # of the allowed domain.  Used in adaptive step-control.
 isoutofdomain(x) = isnan(x)
 
-function make_consistent_types(fn, y0, tspan, btab::Tableau)
+function make_consistent_types(fn, y0, tspan)
     # There are a few types involved in a call to a ODE solver which
     # somehow need to be consistent:
     #
@@ -113,7 +113,6 @@ function make_consistent_types(fn, y0, tspan, btab::Tableau)
     # - Eyf: suitable eltype of y and f(t,y)
     #   --> both of these are set to typeof(y0[1]/(tspan[end]-tspan[1]))
     # - Ty: container type of y0
-    # - btab: tableau with entries converted to Et
 
     # Needed interface:
     # On components: /, -
@@ -138,6 +137,19 @@ function make_consistent_types(fn, y0, tspan, btab::Tableau)
     !isleaftype(Et) && warn("The eltype(tspan) is not a concrete type!  Change type of tspan for better performance.")
     !isleaftype(Eyf) && warn("The eltype(y0/tspan[1]) is not a concrete type!  Change type of y0 and/or tspan for better performance.")
 
+    return Et, Eyf, Ty
+end
+
+
+function make_consistent_types(fn, y0, tspan, btab::Tableau)
+    # Returns
+    # - Et: eltype of time, needs to be a real "continuous" type, at
+    #       the moment a AbstractFloat
+    # - Eyf: suitable eltype of y and f(t,y)
+    #   --> both of these are set to typeof(y0[1]/(tspan[end]-tspan[1]))
+    # - Ty: container type of y0
+    # - btab: tableau with entries converted to Et
+    Et, Eyf, Ty = make_consistent_types(fn,y0,tspan)
     btab_ = convert(Et, btab)
     return Et, Eyf, Ty, btab_
 end
@@ -147,48 +159,7 @@ end
 ###############################################################################
 
 include("runge_kutta.jl")
-
-# ODE_MS Fixed-step, fixed-order multi-step numerical method
-#   with Adams-Bashforth-Moulton coefficients
-function ode_ms(F, x0, tspan, order::Integer)
-    h = diff(tspan)
-    x = Array(typeof(x0), length(tspan))
-    x[1] = x0
-
-    if 1 <= order <= 4
-        b = ms_coefficients4
-    else
-        b = zeros(order, order)
-        b[1:4, 1:4] = ms_coefficients4
-        for s = 5:order
-            for j = 0:(s - 1)
-                # Assign in correct order for multiplication below
-                #  (a factor depending on j and s) .* (an integral of a polynomial with -(0:s), except -j, as roots)
-                p_int = polyint(poly(diagm(-[0:j - 1; j + 1:s - 1])))
-                b[s, s - j] = ((-1)^j / factorial(j)
-                               / factorial(s - 1 - j) * polyval(p_int, 1))
-            end
-        end
-    end
-
-    # TODO: use a better data structure here (should be an order-element circ buffer)
-    xdot = similar(x)
-    for i = 1:length(tspan)-1
-        # 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]
-        for j = 1:steporder
-            x[i+1] += h[i]*b[steporder, j]*xdot[i-(steporder-1) + (j-1)]
-        end
-    end
-    return vcat(tspan), x
-end
-
-# Use order 4 by default
-ode4ms(F, x0, tspan) = ode_ms(F, x0, tspan, 4)
-ode5ms(F, x0, tspan) = ODE.ode_ms(F, x0, tspan, 5)
+include("adams_methods.jl")
 
 ###############################################################################
 ## STIFF SOLVERS
@@ -419,10 +390,5 @@ ode4s_s(F, x0, tspan; jacobian=nothing) = oderosenbrock(F, x0, tspan, s4_coeffic
 # Use Shampine coefficients by default (matching Numerical Recipes)
 const ode4s = ode4s_s
 
-const ms_coefficients4 = [ 1      0      0     0
-                          -1/2    3/2    0     0
-                          5/12  -4/3  23/12 0
-                          -9/24   37/24 -59/24 55/24]
-
 
 end # module ODE