-1

I'm currently developing a model to integrate microbiology and geochemistry that uses lsoda to solve a great number of differential equations. The model is far too large to be posted here because it's made of several modules, but I have some very weird happening.

These are my initial values enter image description here

I have initialised them as zero because I don't want any kind of microbe, just to check how the chemistry would change without microbes. However after 5 or 6 steps I start seeing values that are different from zero in some of my microbial variables: enter image description here.

I wonder if maybe lsoda is doing some kind of round and that's why I get these values, because otherwise I cannot explain where these values are popping out from. If this is the case, does anyone know how to stop this kind of round-ups?

Mario Toubes-Rodrigo
  • 11
  • 3
    Welcome to Stack Overflow. Please don’t use images of code or data as they cannot be used without a lot of unnecessary effort. It’s really helpful if your question is reproducible. [How to ask a good question](https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example) – Peter Sep 01 '21 at 13:47
  • 2
    I appreciate that your model is complex, but it's basically going to be impossible to answer this question without seeing more detail. Floating-point math is intrinsically imprecise; while (for example) multiplying by 0 should always give zero, and adding 0 should never change a value, most other computations will incur roundoff error (e.g. `sqrt(2)^2-2 !=0`, and other more subtle examples ...) It may be worth trying Euler integration (`ode(..., method="euler")` to minimize the chances that anything complicated is happening elsewhere. – Ben Bolker Sep 01 '21 at 15:02
  • 1
    Check if that behavior changes if you tighten the error tolerances. Is the absolute tolerance compatible with the scale of the components? Produce the result of the ODE function along these points, non-zero components where zero is expected may hint at some programming error. – Lutz Lehmann Sep 01 '21 at 18:07
  • 1
    Please provide enough code so others can better understand or reproduce the problem. – Community Sep 04 '21 at 22:20

1 Answers1

2

Short answer

No, lsoda does not create "fictional" values.

Detailed Answer

As @Ben Bolker stated, "floating point arithmetics is intrinsically imprecise". I want to add that all solvers including lsoda compute approximations, and that nonlinear models can magnify errors.

You can check your model, if all derivatives are exactly zero by calling the model function separately, without the solver, see the following example of a Lorenz system, where all states are set to zero:

library(deSolve)

Lorenz <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    dX <-  a * X + Y * Z
    dY <-  b * (Y - Z)
    dZ <- -X * Y + c * Y - Z
    list(c(dX, dY, dZ))
  })
}

parameters <- c(a = -8/3, b = -10, c = 28)
state      <- c(X = 0, Y = 0, Z = 0)
times      <- seq(0, 500, by = 1)

range(Lorenz(0, state, parameters))

that gives

[1] 0 0

Counter example

If we modify the model, so that one derivative is slightly different from zero due to rounding error, e.g. by using Ben's example:

Fancy  <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    dX <-  a * X + Y * Z + 27 * (sqrt(X + Y + Z + 2)^2 - 2)
    dY <-  b * (Y - Z)
    dZ <- -X * Y + c * Y - Z
    list(c(dX, dY, dZ))
  })
}
range(Fancy(0, state, parameters))

then it gives:

[1] 0.000000e+00 1.199041e-14

Depending on tolerances, such a model may then fluctuate or even blow up. Surprisingly, a smaller tolerance may sometimes be worse:

out1 <- ode(y = state, times = times, func = Lorenz, parms = parameters, atol=1e-6)
out2 <- ode(y = state, times = times, func = Fancy, parms = parameters, atol=1e-6)
out3 <- ode(y = state, times = times, func = Fancy, parms = parameters, atol=1e-8)

par(mfrow=c(3,1))
plot(out1[,1:2], type="l", main="Lorenz, derivatives = 0")
plot(out2[,1:2], type="l", main="Fancy, small error in derivatives, atol=1e-6")
plot(out3[,1:2], type="l", main="Fancy, small error in derivatives, atol=1e-8")

Conclusion

  • lsodaworks well within the limitations of floating point arithmetics.
  • models should be carefully designed and tested.
tpetzoldt
  • 5,338
  • 2
  • 12
  • 29