Solving steady state in julia

I was looking for some solvers for my example model to find steady state values for states. I found SSRootfind() and DynamicSS(). Upon reading I felt. It is much better to use DynamicSS as for SSRootfind() your initial values assumption for states should be really good.

Here maybe I’m wrong but for my example model as I change initial condition I find different steady state values (I’m using DynamicSS()). In my head I always thought it should be same irrespective of your initial values.

Can someone shed light on this? and also as I’m getting only single value here for states, does this reflect it has only one solution for steadystate?

function F_simple(x, p, t; γ = 1/18, R₀ = 3.0, σ = 1/5.2)
    s, e, i, r = x

    return [-γ*R₀*s*i;       # ds/dt = -γR₀si
             γ*R₀*s*i -  σ*e;# de/dt =  γR₀si -σe
             σ*e - γ*i;      # di/dt =         σe -γi
                   γ*i;      # dr/dt =             γi
            ]
end

i_0 =  0.7           
e_0 = 0.1         
r_0 = 0.2
x_0 = [s_0, e_0, i_0, r_0]  # initial condition

tspan = (0.0, 350.0)  # ≈ 350 days
problem = ODEProblem(F_simple, x_0, tspan)
prob = SteadyStateProblem(problem)
sol = solve(prob, DynamicSS(Rodas5P()))
julia> sol = solve(prob, DynamicSS(Rodas5P()))
u: 4-element Vector{Float64}:
 0.010268836444676122
 1.6803043675976122e-9
 1.365844689856922e-7
 0.98973102529055

Now if I change initial condition:

i_0 =  0.5           
e_0 = 0.4         
r_0 = 0.1

output:

julia> sol = solve(prob, DynamicSS(Rodas5P()))
u: 4-element Vector{Float64}:
 0.007545374053424994
 7.285142390091502e-10
 8.02329983116032e-8
 1.1924545449850648

A multistable dynamical system will have different steady states for different initial conditions. See e.g. Multistability in Physical and Living Systems: Characterization and Applications | SpringerLink .

Likely your dynamical system is multistable? It is non-linear so it is possible.

1 Like

Julia has a software package for analyzing multistable dynamical systems. I believe it is straght-forward to use, you may use your system like the example here:

https://juliadynamics.github.io/Attractors.jl/dev/examples/#Newton’s-fractal-(basins-of-2D-map)