My code works in R but generates an error in DifferentialEquations.jl


#1

Hello again, I am still trying to determine why my model’s population goes below 0 which generates an error of imaginary numbers. I am new to Julia, using the latest version on Atom (1.0 I think) and working on Windows10.

I am trying to recreate a three species food chain but keep running into a Exponentiation yielding a complex result requires a complex argument. Replace x^y with (x+0im)^y, Complex(x)^y, or similar error in Julia.

But, my code works in R so I am stumped as what to do.

Julia code (that works)

using DifferentialEquations
f = @ode_def_bare FoodChain begin
dBB = r*(1- BB/K)*BB -xI*yIB*BI*((BB^h/(B0^h + BB^h))/eIB)
dBI = -xI*BI + xI*yIB*BI*(BB^h/(B0^h + BB^h))-xT *yTI*BT*((BI^h/(B0^h + BI^h))/eTI)
dBT = -xT*BT + xT*yTI*BT*(BI^h/(B0^h + BI^h))
end r K xI yIB h B0 eIB xT yTI eTI q E
u0 = [155.0,107.0,93.0]
tspan = (0.0,1000.0)
p = (1.1, 700, 0.18, 10, 1.2, 80, 0.66, 0.06, 10, 0.85)
# r K xI yIB h B0 eIB xT yTI eTI
prob = ODEProblem(f,u0,tspan,p)
sol=solve(prob, saveat=0.05)
 
using Plots
 
plot(sol,xlabel = "Time" ,ylabel = "Density", lw=0.5, layout = (3,1))

If I change r from 1.1 to 0.3 (first parameter), it will generate the previously mentioned error BUT it will work in R. In R the populations do not go below zero, which I thought is the reason for triggering the Julia error.

I am completely discombobulated.

library("deSolve")
library("EMD") # this package provides a function that returns local maxima and minima of a time series, not only the global extrema

model.rma <-function(t, x, parms)
{
  BB<-x[1]
  BI<-x[2]
  BT<-x[3]
  
  dBB <- r*(1- BB/K)*BB -xI*yIB*BI*((BB^h/(B0^h + BB^h))/eIB)
  dBI <- -xI*BI + xI*yIB*BI*(BB^h/(B0^h + BB^h))-xT *yTI*BT*((BI^h/(B0^h + BI^h))/eTI)
  dBT <- -xT*BT + xT*yTI*BT*(BI^h/(B0^h + BI^h))
  
  return(list(c(dBB,dBI,dBT)))
}

tmax <- 1000
BB.0 <- 155.0   # initial values
BI.0 <- 107.0
BT.0 <- 93.0

r <- 0.3
K <- 700
xI <- 0.18      # metabolic rate of species 2
yIB <- 10     # maximum ingestion rate (relative to metabolic rate) of species 2
h <- 1.2     #hill exponent
B0 <- 80 # half-saturation density
eIB <- 0.66          # assimilation efficiency
xT <- 0.06      # metabolic rate of species 3
yTI <- 10         # maximum ingestion rate (relative to metabolic rate) of species 3
eTI <- 0.85          # assimilation efficiency

#parms 
xstart <-c(BB=BB.0, BI=BI.0, BT=BT.0)

times <-seq(0, tmax, 0.05)

out1 <-as.data.frame(lsoda(xstart, times, model.rma, parms))

par(mfrow=c(3,1))
par(las=1)
par(pty="m") 

# Check if some values are lower than 0

# for basal species
index.lower.than.zero = which(out1$BB <= 10e-4)
out1$time[index.lower.than.zero]

# for intermediate species
index.lower.than.zero = which(out1$BI <= 10e-4)
out1$time[index.lower.than.zero]

# for Tertiary species
index.lower.than.zero = which(out1$BT <= 10e-4)
out1$time[index.lower.than.zero]

# plot time series
plot(
  c(out1$time,out1$time,out1$time),
  c(out1$BB,out1$BI,out1$BT),
  xlab="time",
  ylab="prey density",
  type="n"
)
lines(out1$time,out1$BB,col="darkgreen")

plot(
  c(out1$time,out1$time,out1$time),
  c(out1$BB,out1$BI,out1$BT),
  xlab="time",
  ylab="consumer density",
  type="n"
)
lines(out1$time,out1$BI,col="orange")

plot(
  c(out1$time,out1$time,out1$time),
  c(out1$BB,out1$BI,out1$BT),
  xlab="time",
  ylab="top predator density",
  type="n"
)
lines(out1$time,out1$BT,col="red")

#3

The system of equations appears to be stiff when r == 0.3. In a nutshell, this means that the rates of change in the ODE system are fast enough that the solver can overshoot the true solution if it’s using the wrong algorithm or the step sizes are too big. In your case, one of the population densities goes negative, which is obviously impossible, and leads to the error.

If you look at the help file for the lsoda function in deSolve in R, it says that it’s calling the Fortran library of the same name, and, critically, switches automatically between stiff and non-stiff methods. solve in DifferentialEquations doesn’t automatically detect stiffness doesn’t seem to detect that stiffness in this case, but you can warn it ahead of time by providing the optional alg_hints argument:

sol = solve(prob, alg_hints=[:stiff], saveat=0.05)

The above change makes the model work for me.

More generally, the issue of stiffness/numerical stability is always present when solving ODEs/PDEs numerically. If your solution is producing unrealistic values (e.g. negative, infinite, etc. depending on the context) a first step should always be checking to make sure a) that your parameter values are correct and b) the timestep and/or integration algorithm are appropriate. One of the nice things about DifferentialEquations is that it makes it easy to try different solvers and step sizes if the defaults don’t work.

Edit: typo


#4

solve does automatically detect stiffness, since in this case it defaults to AutoTsit5

http://docs.juliadiffeq.org/latest/solvers/ode_solve.html#Pre-Built-Stiffness-Detecting-and-Auto-Switching-Algorithms-1

Which should automatically switch to Rosenbrock23() when the stiffness is detected mid run. It seems to be a little too strict to switch here. I’ll work with @YingboMa to get this tweaked, fixed, and add it to our tests.


#5

See the development discussion in https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/issues/566 .

Basically, LSODA’s stiffness detection worked here and ours didn’t with the defaults because they use a much lower default tolerance on the integrator, reltol=1e-6,abstol=1e-6 instead of reltol=1e-3,abstol=1e-6. If you pass in reltol=1e-3 (default abstol=1e-6) it works.

reltol=1e-3 is a default that’s around in most places like MATLAB because it means you get fast results to “plotting accuracy” (about 3 decimal places, so beyond the human eye) as fast as possible. Then anyone can choose to make it more accurate. But there are a few cases where an integrator will slip up with that default. I’m not entirely sure we should change the default here unless this comes up often (in which case, we may just want to special case the stiffness detection algorithms), so for now my recommendation is to just try lower reltols if this comes up and we’ll see if we can make something a bit more robust without effecting default speed too much.


#6

In the mean time…

sol=solve(prob,Rosenbrock23(),saveat=0.05)

works and leads to:


#7

Ah, thanks for the clarification.


#8

… in addition to Rosenbrock23(), the following methods work with r=0.3: