Added tests

src/rosw.jl

@@ -7,7 +7,9 @@
 export ode_rosw, ode_rosw_fixed
 
 # TODO:
-# - AD Jacobian
+# - Jacobian coloring: https://github.com/mlubin/ReverseDiffSparse.jl
+#      http://wiki.cs.purdue.edu:9835/coloringpage/abstracts/acyclic-SISC.pdf
+# - AD Jacobian: https://github.com/mlubin/ReverseDiffSparse.jl
 # - fix fixed step solver
 
 # Rosenbrock-W methods are typically specified for autonomous DAE:
@@ -172,7 +174,9 @@ function oderosw_adapt{N,S}(fn, x0::AbstractVector, tspan, btab::TableauRosW{N,S
 
     # FIXME: add interpolation
     if length(tspan)>2
-        error("specified output times not supported yet")
+        warn("specified output times not supported yet")
+        tspan = [tspan[1], tspan[end]]
+        points=:all
     else
         points=:all
     end

src/runge_kutta.jl

@@ -464,7 +464,6 @@ function hermite_interp!(y, tquery,t,dt,y0,y1,f0,f1)
     #
     # f_0 = f(x_0 , y_0) , f_1 = f(x_0 + h, y_1 )
     # this is O(3). TODO for higher order.
-error("")
     theta = (tquery-t)/dt
     for i=1:length(y0)
         y[i] = ((1-theta)*y0[i] + theta*y1[i] + theta*(theta-1) *

test/implicit-eqs.jl

@@ -0,0 +1,78 @@
+using ODE
+using Base.Test
+
+tol = 1e-2
+
+solvers = [
+           ## Stiff
+           ODE.ode_rosw
+           ]
+
+for solver in solvers
+    println("using $solver")
+    # dy
+    # -- = 6 ==> y = 6t
+    # dt
+    f = (res,y,ydot)->(res[:]=6.0-ydot; nothing)
+    t,y=solver(f, [0.], [0.,1])
+    y = vcat(y...)
+    @test maximum(abs(y-6t)) < tol
+
+    # dy
+    # -- = 2t ==> y = t.^2
+    # dt
+    # t,y=solver((t,y)->2t, 0., [0:.001:1;])
+    # @test maximum(abs(y-t.^2)) < tol
+
+
+    # dy
+    # -- = y ==> y = y0*e.^t
+    # dt
+    f = (res,y,ydot)->(res[:]=y-ydot;nothing)
+    t,y=solver(f, [1.], [0,1])
+    y = vcat(y...)
+    @test maximum(abs(y-e.^t)) < tol
+
+    t,y=solver((res,y,ydot)->(res[:]=y-ydot), [1.], [1, 0])
+    y = vcat(y...)
+    @test maximum(abs(y-e.^(t-1))) < tol
+
+    # dv       dw
+    # -- = -w, -- = v ==> v = v0*cos(t) - w0*sin(t), w = w0*cos(t) + v0*sin(t)
+    # dt       dt
+    #
+    # y = [v, w]
+
+    function ff(res,y,ydot)
+        res[1] = -y[2]-ydot[1]
+        res[2] = y[1] -ydot[2]
+    end
+
+    t,y=solver(ff, [1., 2.], [0, 2*pi])
+    ys = hcat(y...).'   # convert Vector{Vector{Float}} to Matrix{Float}
+    @test maximum(abs(ys-[cos(t)-2*sin(t) 2*cos(t)+sin(t)])) < tol
+end
+
+# rober testcase from http://www.unige.ch/~hairer/testset/testset.html
+let
+    println("ROBER test case")
+    function f(t, y)
+        ydot = similar(y)
+        ydot[1] = -0.04*y[1] + 1.0e4*y[2]*y[3]
+        ydot[3] = 3.0e7*y[2]*y[2]
+        ydot[2] = -ydot[1] - ydot[3]
+        ydot
+    end
+    f!(res, y, ydot) = (res[:] = f(0, y) - ydot; nothing)
+
+    t = [0., 1e11]
+    t,y = ode_rosw(f!, [1.0, 0.0, 0.0], t; abstol=1e-8, reltol=1e-8,
+                                        maxstep=1e11/10, minstep=1e11/1e18)
+
+    refsol = [0.2083340149701255e-07,
+              0.8333360770334713e-13,
+              0.9999999791665050] # reference solution at tspan[2]
+    @test norm(refsol-y[end], Inf) < 3e-10
+end
+
+println("All looks OK")

test/runtests.jl

@@ -81,5 +81,6 @@ let
     @test norm(refsol-y[end], Inf) < 2e-10
 end
 include("interface-tests.jl")
+include("implicit-eqs.jl")
 
 println("All looks OK")