using ModelingToolkit
@variables t x(t) # independent and dependent variables
@parameters τ # parameters
D = Differential(t) # define an operator for the differentiation w.r.t. time
# your first ODE, consisting of a single equation, the equality indicated by ~
@named fol = ODESystem([ D(x) ~ (1 - x)/τ + 0.5*(t>5)])
prob = ODEProblem(fol, [x => 0.0], (0.0,10.0), [τ => 3.0])
prob1 = SteadyStateProblem(prob)
sol = solve(prob1, SSRootfind())
expect result should x=1,but result is x=2.5。
my understand is when solve steady state problem, only consider t=0，when t=0, steady state equation should be 0 ~ (1 - x)/τ, because 0.5*(t>5) = 0 when t = 0 , but final result not as my understand, so how to deal with this.
Steady state is reached when \frac{\mathrm{d}x}{\mathrm{d}t} = 0, which in your case is satisfied when \frac{1-x}{3} + 0.5 * (t>5) = 0.
Since steady states are associated with long time behavior of a dynamical system (that is, when t\to \infty), I suspect the algorithm is solving the equation
\frac{1-x}{3} + 0.5 = 0,
which has x = 2.5 as solution.
I want to get the correct intial value of x ,in this problem x=1 make the \frac{dx}{dt}=0, then simulation the system from the steady state as \frac{dx}{dt}=0, then add some disturbance to system.
Just to be clear: there is no “correct initial value” per se – the initial value is what it is. From your questions, it appears that what you really want is an initial value so that the state of the system does not change (is “steady”) initially, until the system is “disturbed”. The concept of “steady state” is normally related to the value that is approached when time goes to infinity, assuming the state then approaches a steady solution (some systems instead are unstable with value going to +/- infinity, or become periodic, sometimes referred to as a stationary solution).
With your given system, the steady state is correctly found by Juia (as others have pointed out). Two ways to find what you want (which is not the steady state solution of your listed model):
Find the initial value which is consistent with the initial constraint du/dt = 0, as suggested by others, or
Remove the drive term in your model (t>5) before asking for steady state solution.
Both the above problems can be solved in two ways:
By solving (in general): a set of implicit, nonlinear equations which in principle can have multiple solutions, or
Simulating the system for “a long time” until (hopefully) a steady state is reached. In this second case, the steady state that you find might depend on which initial solution you choose.
help?> SteadyStateProblem
search: SteadyStateProblem SteadyStateProblemExpr
Defines an Defines a steady state ODE problem. Documentation Page:
https://diffeq.sciml.ai/stable/types/steadystatetypes/
Mathematical Specification of a Steady State Problem
======================================================
To define an Steady State Problem, you simply need to give the function f which defines the ODE:
\frac{du}{dt} = f(u,p,t)
and an initial guess u_0 of where f(u,p,t)=0. f should be specified as f(u,p,t) (or in-place as f(du,u,p,t)), and u₀
should be an AbstractArray (or number) whose geometry matches the desired geometry of u. Note that we are not
limited to numbers or vectors for u₀; one is allowed to provide u₀ as arbitrary matrices / higher dimension tensors
as well.
Note that for the steady-state to be defined, we must have that f is autonomous, that is f is independent of t. But
the form which matches the standard ODE solver should still be used. The steady state solvers interpret the f by
fixing t=0.
solve result is not same as document description, “The steady state solvers interpret the f by fixing t=0.”
It seems like the documentation is not the same as the implementation for that as seen here
julia> function f(x,p,t)
@show t
(1 .- x) .* p .+ (0.5*(t>5))
end
f (generic function with 1 method)
julia> prob = ODEProblem(f, 0.0, (0,10), 3.0)
ODEProblem with uType Float64 and tType Int64. In-place: false
timespan: (0, 10)
u0: 0.0
julia> prob1 = SteadyStateProblem(prob)
SteadyStateProblem with uType Float64. In-place: false
u0: 0.0
julia> solve(prob1, SSRootfind())
t = Inf
t = Inf
t = Inf
t = Inf
t = Inf
...
and I have to say I find the implementation more logical than the documentation. For me steady state has always meant that the derivative of the states are and will remain zero for all future t, and seems like that is the definition from wikipedia as well.
But there should for certain be a fix to make docs and implementation be the same no matter what is chosen.
Could add option to calculate the initial steady state in SteadyStateProblem function to implement the document description “The steady state solvers interpret the f by fixing t=0.”
f(x) = 0 can have many solutions, at most one of which is the steady state solution of dx/dt = f(x), $x(0) = x_0 The slow way to get this right is to solve the initial value problem with a high-accuracy integrator. Pseudo-transient continuation is designed for this kind of problem.
My new book explains this in several places. Section 1.7 is a good place to start.
Yes, this is common in drug studies, where the subjects were at a steady state at time = 0 and then perturbed towards a new steady state due to a parameter or forcing function change (typically) starting at time = 0. Of course, time = 0 is really time = infinity with the original parameters and forcing function. However, when the forcing function term is time dependent, then passing time = infinity may be problematic. It’s not a problem if one has a separate function to handle the SS case as I showed in my example in: f does not fix t=0 · Issue #44 · SciML/SteadyStateDiffEq.jl · GitHub. I think it would be worthwhile to have an interface that allowed fixing time to a particular value while solving for the steady state, although my solution shown in issue 44 is an easy fix. But I think my solution would not work in MTK?
I won’t even go into the common case that at time = 0, not only is the system at steady state but the concentrations in a few but not all compartments are known. So it also becomes an identifiability and optimization problem on top of a steady state solve.
That’s not how it’s defined in drug studies. In drug studies, yes you start with steady state but that is defined as the time = infinity result of the no dose model. Then the time is reset to tad = 0 (a different definition of time, time after dose), and doses are applied to the system there.
For this kind of case, I am not sure why one wouldn’t use NonlinearProblem or SteadyStateProblem on the autonomous no-dose system and use events for the time dependence? Indeed, you wouldn’t want to handle that discontinuity in the ODE definition anyways for dosing, you’d want to use a callback to get better convergence properties.
because it is nice and convenient to use SteadyStateProblem to solve nonlinear problem at t = 0, the ODEProblem can directly convert to SteadyStateProblem, but could not directly convert to NonlinearProblem, then use NonlinearSolve to solve.
, it’s not a steady state. Is this for trying to do initialization or something?