Changing the righthand side in a callback

Specific Domains Modelling & Simulations
1

Hi everyone,
I am wondering how to correctly change the right hand side inside a DifferentialEquations callback.
If I don’t reinit the system runs through but gives wrong results (compared to just running new differential equations). If I reinit the integrator I get an infinite loop.
Any pointers?

mutable struct ode_switch
    original_ode
    fault_ode
    fault::Bool
end

function (pgs::ode_switch)(u, t, integrator)
    println("time $t") # To make sure the call back is really called
    if ! pgs.fault
        integrator.f = rhs(pgs.fault_ode)
        integrator.u = find_valid_ic(integrator.f, integrator.u) # This takes care of the DAE/Massmatrix aspect
        # reinit!(integrator, t0=t) # If I use reinit, I get an infinite loop with the callback being called over and over.
        # step!(integrator) # I tried to do a step in order for that not to happen but it didn't help
        pgs.fault = true
    else
        integrator.f = rhs(pgs.original_ode)
        integrator.u = find_valid_ic(integrator.f, integrator.u) # This takes care of the DAE/Massmatrix aspect
        # reinit!(integrator, t0=t) # If I use reinit, I get an infinite loop with the callback being called over and over.
        # step!(integrator) # I tried to do a step in order for that not to happen but it didn't help
        pgs.fault = false
    end
end

timespan = (0., 10.)
faulttime = (1., 2.)
pgs = ode_switch(func1,func2)
cb = FunctionCallingCallback(pgs; funcat=faultime, func_start = false)
prob = ODEProblem(func1, x0, timespan, callback=cb)
solve(prob, Rodas4(autodiff=false))
2

I would just use a DiscerteCallback to change a parameter that is the boolean in p, and then check for that in f. By making it integrator-local it won’t have re-use issues.

3

Thanks for the quick reply, I saw that as a design in the other answers. In this case it’s a bit of an awkward fit because I typically will get an ODEFunction and a fault function that changes the ODEFunction from the user. So it’s a bit of effort to build a wrapper that switches between them, thus I thought I’d try it this way.

4

For example, the different functions (func1, func2) could have different mass matrices, that could not be handled by a change of parameters, correct?

5

That probably could be, since it can be a DiffEqOperator with an update_func, but that is not well documented…

6

In that case I think we will stick to calling solve several times with initial conditions taken from the last simulation.

Is there a documented way to combine several solutions into one? (or the next solve on top of the existing solution object)?

7

There’s no documented way to do that. Instead of doing that, why not use the integrator interface?

http://docs.juliadiffeq.org/latest/basics/integrator.html