The following mass matrix DAE solvers fail to solve when used with callbacks that modify the current state: Rodas4, Rodas42, Rodas4P, Rodas4P2, Rodas5.
MWE:
# original problem
using DifferentialEquations
function f(du, u, p, t)
du[1] = -p[1]*u[1] + p[2]*u[2]*u[3]
du[2] = p[1]*u[1] - p[2]*u[2]*u[3] - p[3]*u[2]*u[2]
du[3] = u[1] + u[2] + u[3] - 1.
end
M = [1 0 0; 0 1 0; 0 0 0]
p = [0.04, 10^4, 3e7]
u0 = [1.,0.,0.]
tspan = (0., 1e6)
prob = ODEProblem(ODEFunction(f, mass_matrix = M), u0, tspan, p)
sol = solve(prob, Rodas5()) # solves fine
# problem with callback that resets the state to the initial conditions
tstops = [0.1]
affect!(integrator) = integrator.u .= u0
cb = PresetTimeCallback(tstops, affect!)
cprob = ODEProblem(ODEFunction(f, mass_matrix = M), u0, tspan, p; callback = cb)
csol = solve(cprob, Rodas5()) # aborts due to instability right after callback (at t = 0.1)
Note that while this example uses Rodas5 as the solver, the same result is obtained when using any of the aforementioned solvers.
I tested the same problem using the following solvers and they do not have the same issue: Ros3P, Rodas3, RosShamp4, Veldd4, Velds4, GRK4T, GRK4A, Ros4LStab
The following mass matrix DAE solvers fail to solve when used with callbacks that modify the current state:
Rodas4
,Rodas42
,Rodas4P
,Rodas4P2
,Rodas5
.MWE:
Note that while this example uses
Rodas5
as the solver, the same result is obtained when using any of the aforementioned solvers.I tested the same problem using the following solvers and they do not have the same issue:
Ros3P
,Rodas3
,RosShamp4
,Veldd4
,Velds4
,GRK4T
,GRK4A
,Ros4LStab