# Ensuring certain outputs in an ODE problem reach steady state

**URL:** <https://discourse.julialang.org/t/ensuring-certain-outputs-in-an-ode-problem-reach-steady-state/116610>\
**Category:** Modelling & Simulations\
**Created:** [July 4, 2024, 12:59am UTC](https://discourse.julialang.org/t/ensuring-certain-outputs-in-an-ode-problem-reach-steady-state/116610 "2024-07-04T00:59:54Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![Nick1](https://avatars.discourse-cdn.com/v4/letter/n/47e85d/32.png) [@Nick1](https://discourse.julialang.org/u/Nick1)\
**Post date:** [July 4, 2024, 12:59am UTC](https://discourse.julialang.org/t/ensuring-certain-outputs-in-an-ode-problem-reach-steady-state/116610/1 "2024-07-04T00:59:55Z")

</div>

I am interested in ensuring that certain outputs in an ODE problem reach a steady state except for a few because they follow an oscillatory behavior. I am familiar with the `SteadyStateProblem` but in my case does not work because of the oscillations of certain variables. Is there any way around it? Thank you!

---

<div class="post-metadata">

**Author:** ![danielwe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielwe/32/35657_2.png) [@danielwe](https://discourse.julialang.org/u/danielwe)\
**Post date:** [July 4, 2024, 3:10am UTC](https://discourse.julialang.org/t/ensuring-certain-outputs-in-an-ode-problem-reach-steady-state/116610/2 "2024-07-04T03:10:31Z")

</div>

You can use the `TerminateSteadyState` callback with a custom `test` function implementing your termination condition: [Steady State Callbacks · DiffEqCallbacks.jl](https://docs.sciml.ai/DiffEqCallbacks/stable/steady_state/)

Use it as follows:

```julia
using OrdinaryDiffEq, DiffEqCallbacks

# define my_test_function, abstol, reltol

prob = ODEProblem(...)
callback = TerminateSteadyState(abstol, reltol, my_test_function)
sol = solve(prob, alg; callback)

```

For more details about callback usage, see [Event Handling and Callback Functions · DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/features/callback_functions/#Using-Callbacks)

---

<div class="post-metadata">

**Author:** ![Nick1](https://avatars.discourse-cdn.com/v4/letter/n/47e85d/32.png) [@Nick1](https://discourse.julialang.org/u/Nick1)\
**Post date:** [July 4, 2024, 10:52pm UTC](https://discourse.julialang.org/t/ensuring-certain-outputs-in-an-ode-problem-reach-steady-state/116610/3 "2024-07-04T22:52:51Z")

</div>

Thank you for your answer! I want to make sure that my\_test\_function is correct:

```julia
abstol = 1e-8;
reltol = 1e-6;

function test_function(integrator, abstol, reltol, min_t)

    testval = SciMLBase.get_du(integrator); 
    testval = testval - integrator.u;

    (abs(testval[68] + testval[73]) < abstol) | (abs(testval[68] + testval[73]) < reltol*abs(integrator.u[68] + integrator.u[73]));
    (abs(testval[76]) < abstol) | (abs(testval[76]) < reltol*abs(integrator.u[76]));
    (abs(testval[34] + testval[65]) < abstol) | (abs(testval[34] + testval[65]) < reltol*abs(integrator.u[34] + integrator.u[65]));
    (abs(testval[3] + testval[4] + testval[54]) < abstol) | (abs(testval[3] + testval[4] + testval[54]) < reltol*abs(integrator.u[3] + integrator.u[4] + integrator.u[54]));
    (abs(testval[35]) < abstol) | (abs(testval[35]) < reltol*abs(integrator.u[35]));
    (abs(testval[28] + testval[30]) < abstol) | (abs(testval[28] + testval[30]) < reltol*abs(integrator.u[28] + integrator.u[30]));

end

cbT = TerminateSteadyState(abstol, reltol, test_function);

```

Essentially I changed the `allDerivPass` function which is the default. Do you think that is correct?

---

<div class="post-metadata">

**Author:** ![danielwe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielwe/32/35657_2.png) [@danielwe](https://discourse.julialang.org/u/danielwe)\
**Post date:** [July 5, 2024, 5:46am UTC](https://discourse.julialang.org/t/ensuring-certain-outputs-in-an-ode-problem-reach-steady-state/116610/4 "2024-07-05T05:46:04Z")

</div>

`allDerivPass` is a great starting point, but I’m not sure your modification is correct.

> [@Nick1](#):
>
> ```julia
> testval = testval - integrator.u;
> 
> ```

This line should only be included if your problem is a discrete problem, i.e., if you’re studying a discrete dynamical system and using `DiscreteProblem` instead of `ODEProblem` in your code. Is this your use case?

> [@Nick1](#):
>
> ```julia
> (abs(testval[68] + testval[73]) < abstol) | (abs(testval[68] + testval[73]) < reltol*abs(integrator.u[68] + integrator.u[73]));
> (abs(testval[76]) < abstol) | (abs(testval[76]) < reltol*abs(integrator.u[76]));
> (abs(testval[34] + testval[65]) < abstol) | (abs(testval[34] + testval[65]) < reltol*abs(integrator.u[34] + integrator.u[65]));
> (abs(testval[3] + testval[4] + testval[54]) < abstol) | (abs(testval[3] + testval[4] + testval[54]) < reltol*abs(integrator.u[3] + integrator.u[4] + integrator.u[54]));
> (abs(testval[35]) < abstol) | (abs(testval[35]) < reltol*abs(integrator.u[35]));
> (abs(testval[28] + testval[30]) < abstol) | (abs(testval[28] + testval[30]) < reltol*abs(integrator.u[28] + integrator.u[30]));
> 
> ```

Only the last of these lines will be the return value from the function and determine when your integration terminates. The other lines are discarded. I assume you want to take the union of these conditions and terminate if all of them are true. In that case, I would write your function something like this (assuming you’re working with an ODEProblem, not a DiscreteProblem):

```julia
function checktol(testval, u, indices, abstol, reltol)
    val = abs(sum(testval[i] for i in indices))
    scale = abs(sum(u[i] for i in indices))
    return (val < abstol) | (val < reltol * scale)
end

function test_function(integrator, abstol, reltol, min_t)
    if (min_t !== nothing) && (integrator.t < min_t)
        return false
    end
    testval = SciMLBase.get_du(integrator)
    return all(
        checktol(testval, integrator.u, indices, abstol, reltol)
        for indices in [[68, 73], [76], [34], [3, 4, 54], [28, 30]]
    )
end

```

I’ve simplified `allDerivPass` somewhat in that you can only pass scalars as `abstol` and `reltol` here. If you want to pass arrays to specify different tolerances for different components, you need to make a small modification to the iteration. You can look at `allDerivPass` to figure it out (they use `Iterators.cycle` to allow both scalars and arrays).

Note that semicolons are not needed in Julia.

(Also note that I had pasted the wrong link at the end of my previous post. See the updated version.)

---

<div class="post-metadata">

**Author:** ![Nick1](https://avatars.discourse-cdn.com/v4/letter/n/47e85d/32.png) [@Nick1](https://discourse.julialang.org/u/Nick1)\
**Post date:** [July 7, 2024, 11:31pm UTC](https://discourse.julialang.org/t/ensuring-certain-outputs-in-an-ode-problem-reach-steady-state/116610/5 "2024-07-07T23:31:19Z")

</div>

Thank you very much for your response! The code works smoothly. I also compared its output with the original way I was doing with the `SteadyStateProblem` and both give almost identical results but it is much faster now.