Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
30 changes: 30 additions & 0 deletions src/systems/nonlinear/nonlinearsystem.jl
Original file line number Diff line number Diff line change
Expand Up @@ -232,6 +232,32 @@ function NonlinearSystem(eqs; kwargs...)
return NonlinearSystem(eqs, collect(allunknowns), collect(new_ps); kwargs...)
end

"""
$(TYPEDSIGNATURES)

Convert an `ODESystem` to a `NonlinearSystem` solving for its steady state (where derivatives are zero).
Any differential variable `D(x) ~ f(...)` will be turned into `0 ~ f(...)`. The returned system is not
simplified. If the input system is `complete`d, then so will the returned system.
"""
function NonlinearSystem(sys::AbstractODESystem)
eqs = equations(sys)
obs = observed(sys)
subrules = Dict(D(x) => 0.0 for x in unknowns(sys))
eqs = map(eqs) do eq
fast_substitute(eq, subrules)
end

nsys = NonlinearSystem(eqs, unknowns(sys), [parameters(sys); get_iv(sys)];
parameter_dependencies = parameter_dependencies(sys),
defaults = merge(defaults(sys), Dict(get_iv(sys) => Inf)), guesses = guesses(sys),
initialization_eqs = initialization_equations(sys), name = nameof(sys),
observed = obs)
if iscomplete(sys)
nsys = complete(nsys; split = is_split(sys))
end
return nsys
end

function calculate_jacobian(sys::NonlinearSystem; sparse = false, simplify = false)
cache = get_jac(sys)[]
if cache isa Tuple && cache[2] == (sparse, simplify)
Expand Down Expand Up @@ -529,6 +555,10 @@ function DiffEqBase.NonlinearProblem{iip}(sys::NonlinearSystem, u0map,
return remake(NonlinearProblem{iip}(f, u0, p, pt; filter_kwargs(kwargs)...))
end

function DiffEqBase.NonlinearProblem(sys::AbstractODESystem, args...; kwargs...)
NonlinearProblem(NonlinearSystem(sys), args...; kwargs...)
end

"""
```julia
DiffEqBase.NonlinearLeastSquaresProblem{iip}(sys::NonlinearSystem, u0map,
Expand Down
50 changes: 50 additions & 0 deletions test/nonlinearsystem.jl
Original file line number Diff line number Diff line change
Expand Up @@ -390,3 +390,53 @@ end
@test !any(isequal(p[1]), parameters(sys))
@test is_parameter(sys, p)
end

@testset "Can convert from `ODESystem`" begin
@variables x(t) y(t)
@parameters p q r
@named sys = ODESystem([D(x) ~ p * x^3 + q, 0 ~ -y + q * x - r], t;
defaults = [x => 1.0, p => missing], guesses = [p => 1.0],
initialization_eqs = [p^3 + q^3 ~ 4r], parameter_dependencies = [r ~ 3p])
nlsys = NonlinearSystem(sys)
defs = defaults(nlsys)
@test length(defs) == 3
@test defs[x] == 1.0
@test defs[p] === missing
@test isinf(defs[t])
@test length(guesses(nlsys)) == 1
@test guesses(nlsys)[p] == 1.0
@test length(initialization_equations(nlsys)) == 1
@test length(parameter_dependencies(nlsys)) == 1
@test length(equations(nlsys)) == 2
@test all(iszero, [eq.lhs for eq in equations(nlsys)])
@test nameof(nlsys) == nameof(sys)
@test !ModelingToolkit.iscomplete(nlsys)

sys1 = complete(sys; split = false)
nlsys = NonlinearSystem(sys1)
@test ModelingToolkit.iscomplete(nlsys)
@test !ModelingToolkit.is_split(nlsys)

sys2 = complete(sys)
nlsys = NonlinearSystem(sys2)
@test ModelingToolkit.iscomplete(nlsys)
@test ModelingToolkit.is_split(nlsys)

sys3 = structural_simplify(sys)
nlsys = NonlinearSystem(sys3)
@test length(equations(nlsys)) == length(ModelingToolkit.observed(nlsys)) == 1

prob = NonlinearProblem(sys3, [q => 2.0])
@test prob.f.initialization_data.initializeprobmap === nothing
sol = solve(prob)
@test SciMLBase.successful_retcode(sol)
@test sol.ps[p^3 + q^3]≈sol.ps[4r] atol=1e-10

@testset "Differential inside expression also substituted" begin
@named sys = ODESystem([0 ~ y * D(x) + x^2 - p, 0 ~ x * D(y) + y * p], t)
nlsys = NonlinearSystem(sys)
vs = ModelingToolkit.vars(equations(nlsys))
@test !in(D(x), vs)
@test !in(D(y), vs)
end
end
Loading