Zombies in biological ODE : why is my solver not sticking to zero?

Specific Domains Numerics
, , ,
23

No need for the zombie voodoo. I’ve already written down what’s happening here for your colleague. See the issue:

Basically, the issue is that LU-factorization based solving, i.e. A\b, is not guaranteed to preserve zeros. This means that there can be a true zero, but A\b can give 2e-16 or something close to eps(Float64) like that. When that happens, you get a ki = 2e-16 where it should be zero, and that gives a positive (or negative) derivative to a “dead” term.

There are many solutions to it. For one, if you use an explicit integrator then LU-factorization is never used so you’re good. Or if you change to a Krylov method like @isaacsas mentioned then you’re fine. Or you can use the trick from the issue to use mutation in the f function:

    for i in (num_nutrient+1):(index_vessel-1)
        if iszero(p[1][i])
            bioS[i] = 0.0 #safety against https://github.com/SciML/OrdinaryDiffEq.jl/issues/1651
        end
    end

In almost all cases, mutating u is a bad idea. In this specific case, it’s actually a good way to improve numerical accuracy.

5 Likes