Remove ability to store multiple independent variables

Not currently supported.

Author: HarrisonGrodin

README.md

@@ -46,7 +46,7 @@ Each operation builds an `Operation` type, and thus `eqs` is an array of
 analyzed by other programs. We can turn this into a `DiffEqSystem` via:
 
 ```julia
-de = DiffEqSystem(eqs,[t],[x,y,z],[σ,ρ,β])
+de = DiffEqSystem(eqs,t,[x,y,z],[σ,ρ,β])
 de = DiffEqSystem(eqs)
 ```
 

src/equations.jl

@@ -13,7 +13,7 @@ Base.:~(lhs::Number    , rhs::Expression) = Equation(lhs, rhs)
 
 
 _is_dependent(x::Variable) = x.subtype === :Unknown && !isempty(x.dependents)
-_is_parameter(ivs) = x -> x.subtype === :Parameter && x ∉ ivs
+_is_parameter(iv) = x -> x.subtype === :Parameter && x ≠ iv
 _subtype(subtype::Symbol) = x -> x.subtype === subtype
 
 function extract_elements(eqs, predicates)

src/systems/diffeqs/diffeqsystem.jl

@@ -16,26 +16,28 @@ get_args(eq::DiffEq) = Expression[eq.var, eq.rhs]
 
 struct DiffEqSystem <: AbstractSystem
     eqs::Vector{DiffEq}
-    ivs::Vector{Variable}
+    iv::Variable
     dvs::Vector{Variable}
     ps::Vector{Variable}
     jac::RefValue{Matrix{Expression}}
-    function DiffEqSystem(eqs, ivs, dvs, ps)
+    function DiffEqSystem(eqs, iv, dvs, ps)
         jac = RefValue(Matrix{Expression}(undef, 0, 0))
-        new(eqs, ivs, dvs, ps, jac)
+        new(eqs, iv, dvs, ps, jac)
     end
 end
 
 function DiffEqSystem(eqs)
     dvs, = extract_elements(eqs, [_is_dependent])
     ivs = unique(vcat((dv.dependents for dv ∈ dvs)...))
-    ps, = extract_elements(eqs, [_is_parameter(ivs)])
-    DiffEqSystem(eqs, ivs, dvs, ps)
+    length(ivs) == 1 || throw(ArgumentError("one independent variable currently supported"))
+    iv = first(ivs)
+    ps, = extract_elements(eqs, [_is_parameter(iv)])
+    DiffEqSystem(eqs, iv, dvs, ps)
 end
 
-function DiffEqSystem(eqs, ivs)
-    dvs, ps = extract_elements(eqs, [_is_dependent, _is_parameter(ivs)])
-    DiffEqSystem(eqs, ivs, dvs, ps)
+function DiffEqSystem(eqs, iv)
+    dvs, ps = extract_elements(eqs, [_is_dependent, _is_parameter(iv)])
+    DiffEqSystem(eqs, iv, dvs, ps)
 end
 
 
@@ -44,12 +46,12 @@ function generate_ode_function(sys::DiffEqSystem; version::FunctionVersion = Arr
     param_exprs = [:($(sys.ps[i].name) = p[$i]) for i in eachindex(sys.ps)]
     sys_exprs = build_equals_expr.(sys.eqs)
     if version === ArrayFunction
-        dvar_exprs = [:(du[$i] = $(Symbol("$(sys.dvs[i].name)_$(sys.ivs[1].name)"))) for i in eachindex(sys.dvs)]
+        dvar_exprs = [:(du[$i] = $(Symbol("$(sys.dvs[i].name)_$(sys.iv.name)"))) for i in eachindex(sys.dvs)]
         exprs = vcat(var_exprs,param_exprs,sys_exprs,dvar_exprs)
         block = expr_arr_to_block(exprs)
         :((du,u,p,t)->$(toexpr(block)))
     elseif version === SArrayFunction
-        dvar_exprs = [:($(Symbol("$(sys.dvs[i].name)_$(sys.ivs[1].name)"))) for i in eachindex(sys.dvs)]
+        dvar_exprs = [:($(Symbol("$(sys.dvs[i].name)_$(sys.iv.name)"))) for i in eachindex(sys.dvs)]
         svector_expr = quote
             E = eltype(tuple($(dvar_exprs...)))
             T = StaticArrays.similar_type(typeof(u), E)

src/systems/diffeqs/first_order_transform.jl

@@ -12,9 +12,8 @@ end
 
 function ode_order_lowering(sys::DiffEqSystem; kwargs...)
     eqs = sys.eqs
-    ivs = sys.ivs
     eqs_lowered = ode_order_lowering(eqs; kwargs...)
-    DiffEqSystem(eqs_lowered, ivs)
+    DiffEqSystem(eqs_lowered, sys.iv)
 end
 ode_order_lowering(eqs; naming_scheme = "_") = ode_order_lowering!(deepcopy(eqs), naming_scheme)
 function ode_order_lowering!(eqs, naming_scheme)

test/system_construction.jl

@@ -11,7 +11,7 @@ using Test
 eqs = [D(x) ~ σ*(y-x),
        D(y) ~ x*(ρ-z)-y,
        D(z) ~ x*y - β*z]
-de = DiffEqSystem(eqs,[t],[x,y,z],[σ,ρ,β])
+de = DiffEqSystem(eqs,t,[x,y,z],[σ,ρ,β])
 ModelingToolkit.generate_ode_function(de)
 ModelingToolkit.generate_ode_function(de;version=ModelingToolkit.SArrayFunction)
 jac_expr = ModelingToolkit.generate_ode_jacobian(de)
@@ -20,10 +20,11 @@ f = ODEFunction(de)
 ModelingToolkit.generate_ode_iW(de)
 
 # Differential equation with automatic extraction of variables
-de2 = DiffEqSystem(eqs, [t])
+de2 = DiffEqSystem(eqs, t)
 
 function test_vars_extraction(de, de2)
-    for el in (:ivs, :dvs, :ps)
+    @test de.iv == de2.iv
+    for el in (:dvs, :ps)
         names2 = sort(collect(var.name for var in getfield(de2,el)))
         names = sort(collect(var.name for var in getfield(de,el)))
         @test names2 == names
@@ -37,7 +38,7 @@ test_vars_extraction(de, de2)
 @Unknown u(t) u_tt(t) u_t(t) x_t(t)
 eqs = [D3(u) ~ 2(D2(u)) + D(u) + D(x) + 1
        D2(x) ~ D(x) + 2]
-de = DiffEqSystem(eqs, [t])
+de = DiffEqSystem(eqs, t)
 de1 = ode_order_lowering(de)
 lowered_eqs = [D(u_tt) ~ 2u_tt + u_t + x_t + 1
                D(x_t)  ~ x_t + 2
@@ -51,7 +52,7 @@ a = y - x
 eqs = [D(x) ~ σ*a,
        D(y) ~ x*(ρ-z)-y,
        D(z) ~ x*y - β*z]
-de = DiffEqSystem(eqs,[t],[x,y,z],[σ,ρ,β])
+de = DiffEqSystem(eqs,t,[x,y,z],[σ,ρ,β])
 ModelingToolkit.generate_ode_function(de)
 jac = ModelingToolkit.calculate_jacobian(de)
 f = ODEFunction(de)
@@ -75,8 +76,8 @@ ModelingToolkit.generate_nlsys_function(ns)
 _x = y / C
 eqs = [D(x) ~ -A*x,
        D(y) ~ A*x - B*_x]
-de = DiffEqSystem(eqs,[t],[x,y],[A,B,C])
-test_vars_extraction(de, DiffEqSystem(eqs,[t]))
+de = DiffEqSystem(eqs,t,[x,y],[A,B,C])
+test_vars_extraction(de, DiffEqSystem(eqs,t))
 test_vars_extraction(de, DiffEqSystem(eqs))
 @test eval(ModelingToolkit.generate_ode_function(de))([0.0,0.0],[1.0,2.0],[1,2,3],0.0) ≈ -1/3