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")
