From d278c38dda4774726ed260d1cb7d7675d30dca6b Mon Sep 17 00:00:00 2001 From: Mauro Werder Date: Wed, 25 May 2016 17:57:25 +0200 Subject: [PATCH 1/2] Added test for using eltype(tspan)==Int Also fixed ode23s to allow for it. --- src/ODE.jl | 3 ++- test/interface-tests.jl | 11 +++++++---- test/runtests.jl | 2 ++ 3 files changed, 11 insertions(+), 5 deletions(-) diff --git a/src/ODE.jl b/src/ODE.jl index e8ccc55c1..1a6312240 100644 --- a/src/ODE.jl +++ b/src/ODE.jl @@ -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..8c062cb0d 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -40,6 +40,8 @@ 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]) # dy # -- = y ==> y = y0*e.^t From 5999a8c43706002ba1a3f05b889c3f92e3029cd6 Mon Sep 17 00:00:00 2001 From: Mauro Werder Date: Mon, 20 Jun 2016 09:30:38 +0200 Subject: [PATCH 2/2] Also added test for using an Int as y0 and fixed runge-kutta solvers to cope with it. --- src/ODE.jl | 2 +- test/runtests.jl | 6 ++++++ 2 files changed, 7 insertions(+), 1 deletion(-) diff --git a/src/ODE.jl b/src/ODE.jl index 1a6312240..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 diff --git a/test/runtests.jl b/test/runtests.jl index 8c062cb0d..3f761bdf4 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -42,6 +42,12 @@ for solver in solvers # 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