# My code works in R but generates an error in DifferentialEquations.jl
**URL:** https://discourse.julialang.org/t/my-code-works-in-r-but-generates-an-error-in-differentialequations-jl/17778
**Category:** New to Julia
**Created:** [November 20, 2018, 3:45pm UTC](https://discourse.julialang.org/t/my-code-works-in-r-but-generates-an-error-in-differentialequations-jl/17778 "2018-11-20T15:45:17Z")
**Posts on this page:** 7
**Page:** 1
### Author:  [@HelgavonLichtenstein](https://discourse.julialang.org/u/HelgavonLichtenstein)
#### Post date: [November 20, 2018, 3:45pm UTC](https://discourse.julialang.org/t/my-code-works-in-r-but-generates-an-error-in-differentialequations-jl/17778/1 "2018-11-20T15:45:17Z")
Hello again, I am still trying to determine why my model’s population goes below 0 which generates an error of [imaginary numbers](https://discourse.julialang.org/t/something-is-making-julia-think-i-am-using-imaginary-numbers/17341). 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)
```julia
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.
```julia
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")
```
---
### Author:  [@ElOceanografo](https://discourse.julialang.org/u/ElOceanografo)
#### Post date: [November 20, 2018, 4:48pm UTC](https://discourse.julialang.org/t/my-code-works-in-r-but-generates-an-error-in-differentialequations-jl/17778/3 "2018-11-20T16:48:57Z")
The system of equations appears to be [stiff](https://en.wikipedia.org/wiki/Stiff_equation) 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:
```julia
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](http://docs.juliadiffeq.org/latest/solvers/ode_solve.html) and [step sizes](http://docs.juliadiffeq.org/latest/basics/common_solver_opts.html#Stepsize-Control-1) if the defaults don’t work.
Edit: typo
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [November 20, 2018, 6:00pm UTC](https://discourse.julialang.org/t/my-code-works-in-r-but-generates-an-error-in-differentialequations-jl/17778/4 "2018-11-20T18:00:40Z")
> [@ElOceanografo](#):
>
> `solve` in `DifferentialEquations` doesn’t automatically detect stiffness, but you can warn it ahead of time by providing the optional `alg_hints` argument:
`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](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.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [November 20, 2018, 7:16pm UTC](https://discourse.julialang.org/t/my-code-works-in-r-but-generates-an-error-in-differentialequations-jl/17778/5 "2018-11-20T19:16:18Z")
See the development discussion in [https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/issues/566](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.
---
### Author:  [@BLI](https://discourse.julialang.org/u/BLI)
#### Post date: [November 20, 2018, 7:25pm UTC](https://discourse.julialang.org/t/my-code-works-in-r-but-generates-an-error-in-differentialequations-jl/17778/6 "2018-11-20T19:25:51Z")
In the mean time…
```julia
sol=solve(prob,Rosenbrock23(),saveat=0.05)
```
works and leads to:
---
### Author:  [@ElOceanografo](https://discourse.julialang.org/u/ElOceanografo)
#### Post date: [November 20, 2018, 7:41pm UTC](https://discourse.julialang.org/t/my-code-works-in-r-but-generates-an-error-in-differentialequations-jl/17778/7 "2018-11-20T19:41:35Z")
Ah, thanks for the clarification.
---
### Author:  [@BLI](https://discourse.julialang.org/u/BLI)
#### Post date: [November 20, 2018, 7:44pm UTC](https://discourse.julialang.org/t/my-code-works-in-r-but-generates-an-error-in-differentialequations-jl/17778/8 "2018-11-20T19:44:55Z")
… in addition to `Rosenbrock23()`, the following methods work with r=0.3:
- `BS3()`
- `Rodas4()`
- `QNDF()`
- `CVODE_BDF()`
- `Trapezoid()`
- `TRBDF2()`
… in the “Translation from MATLAB/Python/R” part of [http://docs.juliadiffeq.org/stable/solvers/ode\_solve.html](http://docs.juliadiffeq.org/stable/solvers/ode_solve.html)