From 5c72d0869b5ae30aab1c5db81abb05c64a344d6e Mon Sep 17 00:00:00 2001 From: acroy Date: Wed, 26 Feb 2014 13:08:32 +0100 Subject: [PATCH] Duck typing and consistent return types (Vector{typeof(y0)}) for all solvers; see #7. Tests adjusted accordingly. --- src/ODE.jl | 105 ++++++++++++++++++++++++++------------------------ test/tests.jl | 9 +++-- 2 files changed, 60 insertions(+), 54 deletions(-) diff --git a/src/ODE.jl b/src/ODE.jl index a57ba88bc..1f088c96a 100644 --- a/src/ODE.jl +++ b/src/ODE.jl @@ -43,8 +43,7 @@ export ode23, ode4, ode45, ode4s, ode4ms # Initialize variables. # Adapted from Cleve Moler's textbook # http://www.mathworks.com/moler/ncm/ode23tx.m -function ode23{T}(F::Function, tspan::AbstractVector, y0::AbstractVector{T}; reltol = 1.e-5, abstol = 1.e-8) - +function ode23(F, tspan, y0; reltol = 1.e-5, abstol = 1.e-8) if reltol == 0 warn("setting reltol = 0 gives a step size of zero") end @@ -58,7 +57,8 @@ function ode23{T}(F::Function, tspan::AbstractVector, y0::AbstractVector{T}; rel y = y0 tout = t - yout = y.' + yout = Array(typeof(y0),1) + yout[1] = y tlen = length(t) @@ -101,7 +101,7 @@ function ode23{T}(F::Function, tspan::AbstractVector, y0::AbstractVector{T}; rel t = tnew y = ynew tout = [tout; t] - yout = [yout; y.'] + push!(yout, y) s1 = s4 # Reuse final function value to start new step end @@ -118,7 +118,7 @@ function ode23{T}(F::Function, tspan::AbstractVector, y0::AbstractVector{T}; rel end # while (t != tfinal) - return (tout, yout) + return tout, yout end # ode23 @@ -180,9 +180,7 @@ end # ode23 # CompereM@asme.org # created : 06 October 1999 # modified: 17 January 2001 - -function oderkf{T}(F::Function, tspan::AbstractVector, x0::AbstractVector{T}, a, b4, b5; reltol = 1.0e-5, abstol = 1.0e-8) - +function oderkf(F, tspan, x0, a, b4, b5; reltol = 1.0e-5, abstol = 1.0e-8) # see p.91 in the Ascher & Petzold reference for more infomation. pow = 1/5 @@ -196,10 +194,11 @@ function oderkf{T}(F::Function, tspan::AbstractVector, x0::AbstractVector{T}, a, h = (tfinal - t)/100 # initial guess at a step size x = x0 tout = t # first output time - xout = x.' # first output solution + xout = Array(typeof(x0), 1) + xout[1] = x # first output solution - k = zeros(eltype(x), length(c), length(x)) - k[1,:] = F(t,x) # first stage + k = Array(typeof(x0), length(c)) + k[1] = F(t,x) # first stage while t < tfinal && h >= hmin if t + h > tfinal @@ -207,14 +206,14 @@ function oderkf{T}(F::Function, tspan::AbstractVector, x0::AbstractVector{T}, a, end for j = 2:length(c) - k[j,:] = F(t + h.*c[j], x + h.*(a[j,1:j-1]*k[1:j-1,:]).') + k[j] = F(t + h.*c[j], x + h.*(a[j,1:j-1]*k[1:j-1])[1]) end # compute the 4th order estimate - x4 = x + h.*(b4*k).' + x4 = x + h.*(b4*k)[1] # compute the 5th order estimate - x5 = x + h.*(b5*k).' + x5 = x + h.*(b5*k)[1] # estimate the local truncation error gamma1 = x5 - x4 @@ -228,7 +227,7 @@ function oderkf{T}(F::Function, tspan::AbstractVector, x0::AbstractVector{T}, a, t = t + h x = x5 # <-- using the higher order estimate is called 'local extrapolation' tout = [tout; t] - xout = [xout; x.'] + push!(xout, x) # 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), @@ -238,9 +237,9 @@ function oderkf{T}(F::Function, tspan::AbstractVector, x0::AbstractVector{T}, a, # 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[end,:] + k[1] = k[end] else - k[1,:] = F(t,x) # first stage + k[1] = F(t,x) # first stage end end @@ -252,7 +251,7 @@ function oderkf{T}(F::Function, tspan::AbstractVector, x0::AbstractVector{T}, a, println("Step size grew too small. t=", t, ", h=", h, ", x=", x) end - return (tout, xout) + return tout, xout end # Both the Dormand-Prince and Fehlberg 4(5) coefficients are from a tableau in @@ -316,59 +315,65 @@ const ode45 = ode45_dp # ODEFUN(T,X) must return a column vector corresponding to f(t,x). Each # row in the solution array X corresponds to a time returned in the # column vector T. -function ode4{T}(F::Function, tspan::AbstractVector, x0::AbstractVector{T}) +function ode4(F, tspan, x0) h = diff(tspan) - x = Array(T, length(tspan), length(x0)) - x[1,:] = x0 + x = Array(typeof(x0), length(tspan)) + x[1] = x0 - midxdot = Array(T, 4, length(x0)) + midxdot = Array(typeof(x0), 4) for i = 1:length(tspan)-1 # Compute midstep derivatives - midxdot[1,:] = F(tspan[i], x[i,:]') - midxdot[2,:] = F(tspan[i]+h[i]./2, x[i,:]' + midxdot[1,:]'.*h[i]./2) - midxdot[3,:] = F(tspan[i]+h[i]./2, x[i,:]' + midxdot[2,:]'.*h[i]./2) - midxdot[4,:] = F(tspan[i]+h[i], x[i,:]' + midxdot[3,:]'.*h[i]) + midxdot[1] = F(tspan[i], x[i]) + midxdot[2] = F(tspan[i]+h[i]./2, x[i] + midxdot[1].*h[i]./2) + midxdot[3] = F(tspan[i]+h[i]./2, x[i] + midxdot[2].*h[i]./2) + midxdot[4] = F(tspan[i]+h[i], x[i] + midxdot[3].*h[i]) # Integrate - x[i+1,:] = x[i,:] + 1./6.*h[i].*[1 2 2 1]*midxdot + x[i+1] = x[i] + 1./6.*(h[i].*[1 2 2 1]*midxdot)[1] end - return (tspan, x) + return tspan, x end #ODEROSENBROCK Solve stiff differential equations, Rosenbrock method # with provided coefficients. -function oderosenbrock{T}(F::Function, G::Function, tspan::AbstractVector, x0::AbstractVector{T}, gamma, a, b, c) +function oderosenbrock(F, G, tspan, x0, gamma, a, b, c) h = diff(tspan) - x = Array(T, length(tspan), length(x0)) - x[1,:] = x0 + x = Array(typeof(x0), length(tspan)) + x[1] = x0 solstep = 1 while tspan[solstep] < maximum(tspan) ts = tspan[solstep] hs = h[solstep] - xs = reshape(x[solstep,:], size(x0)) + xs = x[solstep] dFdx = G(ts, xs) - jac = eye(size(dFdx,1))./gamma./hs-dFdx + # FIXME + if size(dFdx,1) == 1 + jac = 1/gamma/hs - dFdx[1] + else + jac = eye(dFdx)./gamma./hs - dFdx + end - g = zeros(size(a,1), length(x0)) - g[1,:] = jac \ F(ts + b[1].*hs, xs) + g = Array(typeof(x0), size(a,1)) + g[1] = (jac \ F(ts + b[1].*hs, xs)) for i = 2:size(a,1) - g[i,:] = jac \ (F(ts + b[i].*hs, xs + (a[i,1:i-1]*g[1:i-1,:]).') + (c[i,1:i-1]*g[1:i-1,:]).'./hs) + g[i] = (jac \ (F(ts + b[i].*hs, xs + (a[i,1:i-1]*g[1:i-1])[1]) + (c[i,1:i-1]*g[1:i-1])[1]./hs)) end - - x[solstep+1,:] = x[solstep,:] + b*g + x[solstep+1] = x[solstep] + (b*g)[1] solstep += 1 end - return (tspan, x) + return tspan, x end -function oderosenbrock{T}(F::Function, tspan::AbstractVector, x0::AbstractVector{T}, gamma, a, b, c) +function oderosenbrock(F, tspan, x0, gamma, a, b, c) # Crude forward finite differences estimator as fallback - function jacobian(F::Function, t::Number, x::AbstractVector) + # FIXME: This doesn't really work if x is anything but a Vector or a scalar + function jacobian(F, t, x) ftx = F(t, x) - dFdx = zeros(length(x), length(x)) - for j = 1:length(x) - dx = zeros(size(x)) + lx = max(length(x),1) + dFdx = zeros(eltype(x), lx, lx) + for j = 1:lx + dx = zeros(eltype(x), lx) # The 100 below is heuristic dx[j] = (x[j]+(x[j]==0))./100 dFdx[:,j] = (F(t,x+dx)-ftx)./dx[j] @@ -412,10 +417,10 @@ ode4s_s(F, G, tspan, x0) = oderosenbrock(F, G, tspan, x0, s4_coefficients...) const ode4s = ode4s_s # ODE_MS Fixed-step, fixed-order multi-step numerical method with Adams-Bashforth-Moulton coefficients -function ode_ms{T}(F::Function, tspan::AbstractVector, x0::AbstractVector{T}, order::Integer) +function ode_ms(F, tspan, x0, order::Integer) h = diff(tspan) - x = zeros(T, length(tspan), length(x0)) - x[1,:] = x0 + x = Array(typeof(x0), length(tspan)) + x[1] = x0 if 1 <= order <= 4 b = [ 1 0 0 0 @@ -438,10 +443,10 @@ function ode_ms{T}(F::Function, tspan::AbstractVector, x0::AbstractVector{T}, or 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,:] + b[steporder,1:steporder]*xdot[i-(steporder-1):i,:].*h[i] + xdot[i] = F(tspan[i], x[i]) + x[i+1] = x[i] + (b[steporder,1:steporder]*xdot[i-(steporder-1):i])[1].*h[i] end - return (tspan, x) + return tspan, x end # Use order 4 by default diff --git a/test/tests.jl b/test/tests.jl index 87df8384e..7fe72ed0a 100644 --- a/test/tests.jl +++ b/test/tests.jl @@ -19,19 +19,19 @@ for solver in solvers # dy # -- = 6 ==> y = 6t # dt - t,y=solver((t,y)->6, [0:.1:1], [0.]) + t,y=solver((t,y)->6, [0:.1:1], 0.) @test maximum(abs(y-6t)) < tol # dy # -- = 2t ==> y = t.^2 # dt - t,y=solver((t,y)->2t, [0:.001:1], [0.]) + t,y=solver((t,y)->2t, [0:.001:1], 0.) @test maximum(abs(y-t.^2)) < tol # dy # -- = y ==> y = y0*e.^t # dt - t,y=solver((t,y)->y, [0:.001:1], [1.]) + t,y=solver((t,y)->y, [0:.001:1], 1.) @test maximum(abs(y-e.^t)) < tol # dv dw @@ -40,7 +40,8 @@ for solver in solvers # # y = [v, w] t,y=solver((t,y)->[-y[2], y[1]], [0:.001:2*pi], [1., 2.]) - @test maximum(abs(y-[cos(t)-2*sin(t) 2*cos(t)+sin(t)])) < tol + 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 println("All looks OK")