I have a recurring problem in my ODE-based models. In most case the outputs of the ODE system have to be strictly positive, but because the precision of the solver is not perfect, outputs that should be 0 are sometimes estimated as small negative values. This leads to problems when I use these outputs in further computations, e.g. ratios may become negative or greater than 1, sqrt transformations don’t work, etc.
Until now I used a silly way around, checking for negative values and replacing them by a small positive value:
y = integrate_ode_rk45(ode, init, t0, theta, xrs, xis, 1.0E-7, 1.0E-7, 1.0E3);
for(i in 1:S) for(j in 1:C) y[i,j] = y[i,j]<0 ? 1e-10 : y[i,j];
I found that forcing y to be positive leads to errors and interrupts sampling.
Is there a more principled way to deal with this issue?
the ODE solver controls “0” within the limits of the absolute error only - there are no more guarantees than that. I would suggest you to add something slightly larger than absolute error to your solution (to all states)… this is a bit dirty, but it works just fine. So in your case just add 1E-6 to all outputs and then you will be fine. An alternative is to integrate log(y) instead of y by dividing your ODE through with y - please search the forum, this problem comes up once in a while.
well… ODE is not so often in general, but I have given this answer probably now the 3rd time… but no worries. Probably everyone uses different terms for this.