Get Started

using FindSteadyStates
using DifferentialEquations

Create Differential Equation System

The standard way to create DE system is to create a derivative function in the form of f(du,u,p,t) with preallocated derivatives (du), initial variables ('u'), parameters (p) and time ('t')

function bistable_ode!(du, u, p ,t)
	s1, s2 = u
	K1, K2, k1, k2, k3, k4, n1 , n2  = p
	du[1] = k1 / (1 + (s2/K2)^n1) - k3*s1
	du[2] = k2/  (1 + (s1/K1)^n2) - k4*s2
end

NameTuple may cause error

One may want to use NameTuple like LabelledArrays to create Tuple array and apply to the ode function. However, this may cause getindex undefined error for Jacobian Calculation via ModelingToolkit, but still safe with solving steady states with NameTuple as returns.

Create DE Problem for Solving

de = DEsteady(func = bistable_ode!,
         p =  [1.,1.,20.,20.,5.,5.,4.,4.],
         u0 = [3.,1.],
         method = SSRootfind()
         )

Method for solving steady state

The DEsteady.method will further be used for DifferentialEquation.solve. Both dynamic solver (Dynamicss which reaches steady states by solving DE) and root finding (SSRootfind) can get the stable point. However, only SSRootfind can find the saddle point if the initial values is near by.

The system is now well organized and ready to be solved.

sol = solve(de)

u: 2-element Vector{Float64}:
 3.999999765270408
 0.015564205973794263

Exploring the initial states

To find all the steady states, one needs to sample the initial state by grid or random search. To begin with, the ParameterGrid (grid search) and ParameterRandom (random search) are useful generator for iterating and returning the parameter set.

    param_rand = ParameterRandom(
        methods= [
            Uniform(0.,100.),
            Log_uniform(0.001,100.)
        ],
        len= 100
    )

100-element ParameterRandom:
 [24.35987673478841, 2.346518234169877]
 [12.815682704126697, 1.9670007373318323]
 [39.978536251948825, 99.07942522523452]
 [31.482679067327844, 95.79460114896938]
 [65.14362211440334, 0.018610679106675933]
 [6.403263704023998, 0.004743911569750526]
 [33.6421007178003, 2.0333940766390586]
 [58.28925049633143, 1.114099465349585]
 [24.35938669614075, 0.02999796439283893]
 [8.37601852204144, 0.0012776428013927738]
 ⋮
 [17.385926332034995, 0.0017145924864299248]
 [32.098389458310095, 0.8833154996316388]
 [0.30900110714018325, 0.021759874707872405]
 [56.06934214848218, 0.11007644363961876]
 [83.50494016483592, 72.93141530958678]
 [73.80165362520121, 7.455134494126275]
 [29.591636491270368, 0.26929879719681454]
 [27.776464750956855, 2.8572021084758985]
 [58.649895555775025, 0.0045185366363317166]

     param_grid = ParameterGrid([
            [0.,100.,100],
            [0.,100.,100]
        ])

10000-element ParameterGrid:
 [0.0, 0.0]
 [1.0101010101010102, 0.0]
 [2.0202020202020203, 0.0]
 [3.0303030303030303, 0.0]
 [4.040404040404041, 0.0]
 [5.050505050505051, 0.0]
 [6.0606060606060606, 0.0]
 [7.070707070707071, 0.0]
 [8.080808080808081, 0.0]
 [9.090909090909092, 0.0]
 ⋮
 [91.91919191919193, 100.0]
 [92.92929292929294, 100.0]
 [93.93939393939395, 100.0]
 [94.94949494949496, 100.0]
 [95.95959595959597, 100.0]
 [96.96969696969697, 100.0]
 [97.97979797979798, 100.0]
 [98.98989898989899, 100.0]
 [100.0, 100.0]

Solve solutions with multi-threading

Random Search

    sols_rand = solve(de, param_rand)

EnsembleSolution Solution of length 100 with uType:
SciMLBase.NonlinearSolution{Float64, 1, Vector{Float64}, Vector{Float64}, SteadyStateProblem{Vector{Float64}, true, Vector{Float64}, ODEFunction{true, typeof(Main.bistable_ode!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, SSRootfind{SteadyStateDiffEq.var"#2#4"}, Nothing, Nothing}

Grid Search

    sols_grid = solve(de, param_grid)

EnsembleSolution Solution of length 10000 with uType:
SciMLBase.NonlinearSolution{Float64, 1, Vector{Float64}, Vector{Float64}, SteadyStateProblem{Vector{Float64}, true, Vector{Float64}, ODEFunction{true, typeof(Main.bistable_ode!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, SSRootfind{SteadyStateDiffEq.var"#2#4"}, Nothing, Nothing}