Steady state solutions#

Solving steady state solutions for an ODE system is to find a combination of state variables such that their derivatives are all zeroes. Their are two ways:

  • Running the ODE solver until a steady state is reached (DynamicSS() and DifferentialEquations.jl)

  • Using a root-finding algorithm to find a steady state (SSRootfind() and NonlinearSolve.jl)

Defining a steady state problem#

using DifferentialEquations

model(u, p, t) = 1 - p * u

p = 1.0
u0 = 0.0
prob = SteadyStateProblem(model, u0, p)
alg = DynamicSS(Tsit5())
DynamicSS{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Float64}(Tsit5(; stage_limiter! = trivial_limiter!, step_limiter! = trivial_limiter!, thread = static(false),), Inf)

Solve the problem. The result should be close to 1.0.

sol = solve(prob, alg) ##
retcode: Success
u: 0.9999994687788004

Modelingtoolkit solving nonlinear systems#

Use NonlinearSolve.jl and NonlinearSystem().

using ModelingToolkit
using NonlinearSolve

@variables x y z
@parameters σ ρ β

eqs = [
    0 ~ σ * (y - x),
    0 ~ x * (ρ - z) - y,
    0 ~ x * y - β * z
]

@named ns = NonlinearSystem(eqs, [x, y, z], [σ, ρ, β])

guess = [x => 1.0, y => 0.0, z => 0.0]
ps = [σ => 10.0, ρ => 26.0, β => 8 / 3]
prob = NonlinearProblem(ns, guess, ps)
sol = solve(prob, NewtonRaphson()) ## The results should be all zeroes
retcode: Success
u: 3-element Vector{Float64}:
 0.0
 0.0
 0.0

Another nonlinear system example with structural_simplify().

@parameters t
@variables u1(t) u2(t) u3(t) u4(t) u5(t)

eqs = [
    0 ~ u1 - sin(u5)
    0 ~ u2 - cos(u1)
    0 ~ u3 - hypot(u1, u2)
    0 ~ u4 - hypot(u2, u3)
    0 ~ u5 - hypot(u4, u1)
]

@named sys = NonlinearSystem(eqs, [u1, u2, u3, u4, u5], [])
\[\begin{split} \begin{align} 0 =& \mathrm{u1}\left( t \right) - \sin\left( \mathrm{u5}\left( t \right) \right) \\ 0 =& \mathrm{u2}\left( t \right) - \cos\left( \mathrm{u1}\left( t \right) \right) \\ 0 =& \mathrm{u3}\left( t \right) - \mathrm{hypot}\left( \mathrm{u1}\left( t \right), \mathrm{u2}\left( t \right) \right) \\ 0 =& \mathrm{u4}\left( t \right) - \mathrm{hypot}\left( \mathrm{u2}\left( t \right), \mathrm{u3}\left( t \right) \right) \\ 0 =& \mathrm{u5}\left( t \right) - \mathrm{hypot}\left( \mathrm{u4}\left( t \right), \mathrm{u1}\left( t \right) \right) \end{align} \end{split}\]

You can simplify the problem using structural_simplify(). There will be only one state variable left. The solve can solve the problem faster.

simple_sys = structural_simplify(sys)
prob = NonlinearProblem(simple_sys, [u5 => 0.0])
sol = solve(prob, NewtonRaphson())
retcode: Success
u: 1-element Vector{Float64}:
 1.6069926947050055

The answer should be 1.6 and 1.0.

@show sol[u5] sol[u1];
sol[u5] = 1.6069926947050055
sol[u1] = 
0.9993449829954304

This notebook was generated using Literate.jl.