Open dlfivefifty opened 4 years ago
Hi,
If the first one is real and the second one is complex, you can just use complex types.
Your example works if you just change it to the following:
julia> solve(ODEProblem(((λ,v),_,t) -> [v,-Vp(λ)], [λ_0,v_0], (0.0, 10.0)); reltol=1E-6)
which is funny, because the part (λ,v) in the function definition works. I suspect it does (λ,v) = u internally (in Julia).
I am trying to think about this but I do not see anything wrong. However, I see your point. I didn't know about the support of broadcasted functions in tuples.
Finally, try:
julia> function f(du, (λ,v),_,t)
du[1] = v
du[2] = -Vp(λ)
end
julia> solve(ODEProblem(f, [λ_0,v_0], (0.0, 10.0)); reltol=1E-6)
Thanks!
Tuples don't have vector semantics which is what causes the issue. Here's the amount of type piracy that you'd need to do in order to do this:
using DifferentialEquations, Plots
V = x -> x^2/2 # Potential
Vp = x -> x # Force
λ_0 = 2.3 # initial location
v_0 = 1.2 # initial position
Base.zero(x::NTuple{N}) where {N} = ntuple(i->zero(x[i]),N)
Base.size(x::Tuple) = (length(x),)
Base.:*(x::Float64,y::Tuple) = x .* y
Base.:+(x::NTuple{N},y::NTuple{N}) where {N} = ntuple(i->x[i] + y[i],N)
λv = solve(ODEProblem(((λ,v),_,t) -> (v,-Vp(λ)), (λ_0,v_0), (0.0, 10.0)), Tsit5(); reltol=1E-6)
Honestly, I like those definitions so I'd add them to DiffEqBase, but I am sure that would cause a frenzy haha.
However, to support stiff ODE solvers... you're asking for linear algebra on tuples. Not impossible, but 🤷♂ . Do you have a proposed fix for:
u = (λ_0,v_0)
J = false .* u .* u'
? That would be required to support stiff ODEs. What is the tuple equivalent to a matrix? What factorization code should it call? That all seems to be begging for static arrays.
tuples do not support addition or multiplication but they do support broadcasting addition and multiplication, which is what the "time-stepping" algorithms should be using.
If you broadcast static arrays then you lose a ton of performance right now, so the out-of-place functions do not broadcast.
We can also get auto-switch algorithms to work, if they never switch to the implicit method:
using StaticArrays
OrdinaryDiffEq._vec(x::Tuple) = SVector(x)
Base.adjoint(x::Tuple) = adjoint(SVector(x))
Base.:-(x::NTuple{N},y::NTuple{N}) where {N} = ntuple(i->x[i] - y[i],N)
λv = solve(ODEProblem(((λ,v),_,t) -> (v,-Vp(λ)), (λ_0,v_0), (0.0, 10.0)); reltol=1E-6)
For future reference to myself, playing around:
using DifferentialEquations, Plots
V = x -> x^2/2 # Potential
Vp = x -> x # Force
λ_0 = 2.3 # initial location
v_0 = 1.2 # initial position
Base.zero(x::NTuple{N}) where {N} = ntuple(i->zero(x[i]),N)
Base.size(x::Tuple) = (length(x),)
Base.:*(x::Float64,y::Tuple) = x .* y
Base.:+(x::NTuple{N},y::NTuple{N}) where {N} = ntuple(i->x[i] + y[i],N)
λv = solve(ODEProblem(((λ,v),_,t) -> (v,-Vp(λ)), (λ_0,v_0), (0.0, 10.0)), Tsit5(); reltol=1E-6)
u = (λ_0,v_0)
J = false .* _vec(u) .* _vec(u)'
using StaticArrays
OrdinaryDiffEq._vec(x::Tuple) = SVector(x)
Base.adjoint(x::Tuple) = adjoint(SVector(x))
Base.:-(x::NTuple{N},y::NTuple{N}) where {N} = ntuple(i->x[i] - y[i],N)
using ForwardDiff
ForwardDiff.derivative(f,x::Tuple) = ForwardDiff.derivative(f,SVector(x))
λv = solve(ODEProblem(((λ,v),_,t) -> (v,-Vp(λ)), (λ_0,v_0), (0.0, 10.0)); reltol=1E-6)
λv = solve(ODEProblem(((λ,v),_,t) -> (v,-Vp(λ)), (λ_0,v_0), (0.0, 10.0)), TRBDF2(autodiff=false); reltol=1E-6)
If the first one is real and the second one is complex, you can just use complex types
That’s roughly 4x slower. Also that was just a model example.
A simpler approach might be to just create a VectorTuple
type that has vector semantics but allows different types for the entries. Then a Tuple
would auto convert to this type at the beginning. Linear algebra would require some work though... as my actual usage is just in lecture notes this is not urgent at all.
Using an automatic conversion to TupleVec
, I can get it to be good for non-stiff ODEs. It's not robust though, and fails at the stiff ODE solver in some deep ways that really need StaticArray types of properties.
using OrdinaryDiffEq, StaticArrays, BenchmarkTools
struct TupleVec{T}
tup::T
end
Base.:+(x::TupleVec,y::TupleVec) = TupleVec(map(+,x.tup,y.tup))
Base.:-(x::TupleVec,y::TupleVec) = TupleVec(map(-,x.tup,y.tup))
Base.:-(x::TupleVec) = TupleVec(map(-,x.tup))
Base.zero(x::TupleVec) = TupleVec(map(zero,x.tup))
Base.size(x::TupleVec) = (length(x.tup),)
Base.length(x::TupleVec) = length(x.tup)
Base.iterate(x::TupleVec,i) = iterate(x.tup,i)
Base.:*(x::Number,y::TupleVec) = TupleVec(x .* y.tup)
Base.:*(x::TupleVec,y::Number) = TupleVec(x.tup .* y)
Base.:/(x::TupleVec,y::Number) = TupleVec(x.tup ./ y)
Base.eltype(x::TupleVec) = promote_type(map(typeof,x.tup)...)
RecursiveArrayTools.recursive_unitless_bottom_eltype(x::TupleVec) = promote_type(map(x->typeof(one(x)),x.tup)...)
ArrayInterfaceCore.zeromatrix(x::TupleVec) = false .* SVector(x.tup) .* SVector(x.tup)'
Base.vec(x::TupleVec) = x
(f::SciMLBase.ODEFunction)(u::TupleVec,p,t) = TupleVec(f.f(u.tup,p,t))
DiffEqBase.function get_concrete_u0(prob, isadapt, t0, kwargs)
if eval_u0(prob.u0)
u0 = prob.u0(prob.p, t0)
elseif haskey(kwargs, :u0)
u0 = kwargs[:u0]
else
u0 = prob.u0
end
if u0 isa Tuple
u0 = SciMLBase.TupleVec(prob.u0)
end
isadapt && eltype(u0) <: Integer && (u0 = float.(u0))
_u0 = handle_distribution_u0(u0)
if isinplace(prob) && (_u0 isa Number || _u0 isa SArray || _u0 isa SciMLBase.TupleVec)
throw(IncompatibleInitialConditionError())
end
_u0
end
V(x) = x^2/2 # Potential
Vp(x) = x # Force
λ_0 = 2.3 # initial location
v_0 = 1.2 # initial position
function ftup((λ,v),_,t)
(v,-Vp(λ))
end
function fsa(u,_,t)
SA[u[2],-Vp(u[1])]
end
@benchmark solve(ODEProblem(ftup, (λ_0,v_0), (0.0, 10.0)), Tsit5(); reltol=1E-6)
@benchmark solve(ODEProblem(fsa, SA[λ_0,v_0], (0.0, 10.0)), Tsit5(); reltol=1E-6)
using ForwardDiff
ForwardDiff.derivative(f,x::SciMLBase.TupleVec) = ForwardDiff.jacobian((args...)->SVector(f(args...)),SVector(x.tup))
solve(ODEProblem(ftup, (λ_0,v_0), (0.0, 10.0)), TRBDF2(); reltol=1E-6)
For now, I added an informative error:
julia> solve(prob,Tsit5())
ERROR: Tuple type used as a state. Since a tuple does not have vector
properties, it will not work as a state type in equation solvers.
Instead, change your equation from using tuple constructors `()`
to static array constructors `SA[]`. For example, change:
function ftup((a,b),p,t)
return b,-a
end
u0 = (1.0,2.0)
tspan = (0.0,1.0)
ODEProblem(ftup,u0,tspan)
to:
using StaticArrays
function fsa(u,p,t)
SA[u[2],u[1]]
end
u0 = SA[1.0,2.0]
tspan = (0.0,1.0)
ODEProblem(ftup,u0,tspan)
This will be safer and fast for small ODEs. For more information, see:
https://diffeq.sciml.ai/stable/tutorials/faster_ode_example/#Further-Optimizations-of-Small-Non-Stiff-ODEs-with-StaticArrays
Stacktrace:
[1] get_concrete_u0(prob::ODEProblem{Tuple{Float64, Float32}, Tuple{Float64, Float64}, false, SciMLBase.NullParameters, ODEFunction{false, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, isadapt::Bool, t0::Float64, kwargs::Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol}, NamedTuple{(:u0, :p), Tuple{Tuple{Float64, Float32}, SciMLBase.NullParameters}}})
@ DiffEqBase C:\Users\accou\.julia\dev\DiffEqBase\src\solve.jl:907
[2] get_concrete_problem(prob::ODEProblem{Tuple{Float64, Float32}, Tuple{Float64, Float64}, false, SciMLBase.NullParameters, ODEFunction{false, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, isadapt::Bool; kwargs::Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol}, NamedTuple{(:u0, :p), Tuple{Tuple{Float64, Float32}, SciMLBase.NullParameters}}})
@ DiffEqBase C:\Users\accou\.julia\dev\DiffEqBase\src\solve.jl:801
[3] solve_up(prob::ODEProblem{Tuple{Float64, Float32}, Tuple{Float64, Float64}, false, SciMLBase.NullParameters, ODEFunction{false, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, sensealg::Nothing, u0::Tuple{Float64, Float32}, p::SciMLBase.NullParameters, args::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ DiffEqBase C:\Users\accou\.julia\dev\DiffEqBase\src\solve.jl:722
[4] solve_up(prob::ODEProblem{Tuple{Float64, Float32}, Tuple{Float64, Float64}, false, SciMLBase.NullParameters, ODEFunction{false, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, sensealg::Nothing, u0::Tuple{Float64, Float32}, p::SciMLBase.NullParameters, args::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False})
@ DiffEqBase C:\Users\accou\.julia\dev\DiffEqBase\src\solve.jl:713
[5] solve(prob::ODEProblem{Tuple{Float64, Float32}, Tuple{Float64, Float64}, false, SciMLBase.NullParameters, ODEFunction{false, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, args::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}; sensealg::Nothing, u0::Nothing, p::Nothing, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ DiffEqBase C:\Users\accou\.julia\dev\DiffEqBase\src\solve.jl:708
[6] solve(prob::ODEProblem{Tuple{Float64, Float32}, Tuple{Float64, Float64}, false, SciMLBase.NullParameters, ODEFunction{false, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, args::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False})
@ DiffEqBase C:\Users\accou\.julia\dev\DiffEqBase\src\solve.jl:703
[7] top-level scope
@ REPL[3]:1
Fairly regularly (in teaching) I need to turn a second order ODE to a first order one. I'm surprised tuples can't be used to do this, see below. I also don't understand the conceptional reason why: tuples do not support addition or multiplication but they do support broadcasting addition and multiplication, which is what the "time-stepping" algorithms should be using.
I realise StaticArrays gives a work around but this is less than ideal for pedagogical usages, and also practical reasons when one has mixed types (say first argument is complex and the second real).