Ensuring certain outputs in an ODE problem reach steady state

Specific Domains Modelling & Simulations
1

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!

2

You can use the TerminateSteadyState callback with a custom test function implementing your termination condition: Steady State Callbacks · DiffEqCallbacks.jl

Use it as follows:

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

3 Likes
3

Thank you for your answer! I want to make sure that my_test_function is correct:

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?

4

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

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?

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):

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.)

5

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.

1 Like