diff --git a/src/ODE.jl b/src/ODE.jl index e8ccc55c1..d530451a9 100644 --- a/src/ODE.jl +++ b/src/ODE.jl @@ -120,8 +120,8 @@ function make_consistent_types(fn, y0, tspan, btab::Tableau) # On container: eltype, promote_type # On time container: eltype - Ty = typeof(y0) Eyf = typeof(y0[1]/(tspan[end]-tspan[1])) + Ty = typeof(similar(y0, Eyf, 1)) Et = eltype(tspan) @assert Et<:Real @@ -265,7 +265,8 @@ function ode23s(F, y0, tspan; reltol = 1.0e-5, abstol = 1.0e-8, h = tdir * min(abs(h), maxstep) y = y0 - tout = Array(typeof(t), 1) + # below promotion should maybe be replaced by using make_consistent_types + tout = Array(promote_type(typeof(t), Float64), 1) tout[1] = t # first output time yout = Array(typeof(y0), 1) yout[1] = deepcopy(y) # first output solution diff --git a/test/interface-tests.jl b/test/interface-tests.jl index dbb811dcc..c134aae03 100644 --- a/test/interface-tests.jl +++ b/test/interface-tests.jl @@ -45,19 +45,19 @@ ODE.isoutofdomain(y::CompSol) = any(isnan, vcat(y.rho[:], y.x, y.p)) ################################################################################ - + # define RHSs of differential equations # delta, V and g are parameters function rhs(t, y, delta, V, g) H = [[-delta/2 V]; [V delta/2]] - + rho_dot = -im*H*y.rho + im*y.rho*H x_dot = y.p p_dot = -y.x - + return CompSol( rho_dot, x_dot, p_dot) end - + # inital conditons rho0 = zeros(2,2); rho0[1,1]=1.; @@ -75,6 +75,9 @@ for solver in solvers end t,y2 = solver((t,y)->rhs(t, y, delta0, V0, g0), y0, linspace(0., endt, 500)) @test norm(y1[end]-y2[end])<0.1 + + # test that typeof(tspan)==Vector{Int} does not throw: + t,y2 = solver((t,y)->rhs(t, y, delta0, V0, g0), y0, [0,1]) end println("ok.") diff --git a/test/runtests.jl b/test/runtests.jl index 9c88cd281..3f761bdf4 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -40,6 +40,14 @@ for solver in solvers t,y=solver((t,y)->2t, 0., [0:.001:1;]) @test maximum(abs(y-t.^2)) < tol + # test typeof(tspan)==Vector{Int} does not throw + t,y=solver((t,y)->2t, 0., [0,1]) + if !(solver in [ODE.ode4s_s, ODE.ode4s_kr, ODE.ode23s]) + # test typeof(y0)==Vector{Int} does not throw + t,y=solver((t,y)->[2t], [0], [0,1]) + # test typeof(y0)==Int does not throw + t,y=solver((t,y)->2t, 0, [0,1]) + end # dy # -- = y ==> y = y0*e.^t