Something is making Julia think I am using imaginary numbers?

New to Julia
1

ERROR: DomainError with -123.74760920723358: Exponentiation yielding a complex result requires a complex argument. Replace x^y with (x+0im)^y, Complex(x)^y, or similar

I figure I should start a new topic since I am running into a different error now. I hope by creating a new topic with the error in the heading it will help someone in the future.

EDIT: Based on a discussion with a colleague, the growth is too high too quickly and Julia cannot contain such large numbers? I can only plot when the time step is reduced to 1 second, but I don’t know what is causing such large growth in the formulas. I should not be generating any imaginary/complex numbers

Continuing from here my code is now

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 + xI*yTI*BT*(BI^h/(B0^h + BI^h))-q*BT*E
end r K xI yIB h B0 eIB xT yTI eTI q E
u0 = [155.0,107.0,93.0]
tspan = (0.0,10000.0)
p = (1.1, 450, 0.18, 10, 1.2, 80, 0.66, 0.06, 10, 0.85, 0.01, 23)
# r K xI yIB h B0 eIB xT yTI eTI q E
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))

And the error is

ERROR: LoadError: DomainError with -1.534065831079686:
Exponentiation yielding a complex result requires a complex argument.
Replace x^y with (x+0im)^y, Complex(x)^y, or similar.
Stacktrace:
 [1] throw_exp_domainerror(::Float64) at .\math.jl:35
 [2] ^ at .\math.jl:782 [inlined]
 [3] macro expansion at C:\owncloud\FoodWebsAnthropocene\Modelling\Julia\Julia foodchain code:4 [inlined]
 [4] (::getfield(Main, Symbol("##23#24")))(::Array{Float64,1}, ::Array{Float64,1}, ::Tuple{Float64,Int64,Float64,Int64,Float64,Int64,Float64,Float64,Int64,Float64,Float64,Int64}, ::Float64) at C:\Users\Helga\.julia\packages\ParameterizedFunctions\ozDxQ\src\ode_def_opts.jl:249
 [5] FoodChain at C:\Users\Helga\.julia\packages\ParameterizedFunctions\ozDxQ\src\maketype.jl:82 [inlined]
 [6] perform_step!(::OrdinaryDiffEq.ODEIntegrator{CompositeAlgorithm{Tuple{Tsit5,Rosenbrock23{0,false,LinSolveFactorize{typeof(LinearAlgebra.lu!)},DataType}},AutoSwitch{Tsit5,Rosenbrock23{0,false,LinSolveFactorize{typeof(LinearAlgebra.lu!)},DataType},Rational{Int64},Float64}},Array{Float64,1},Float64,Tuple{Float64,Int64,Float64,Int64,Float64,Int64,Float64,Float64,Int64,Float64,Float64,Int64},Float64,Float64,Float64,Array{Array{Float64,1},1},OrdinaryDiffEq.ODECompositeSolution{Float64,2,Array{Array{Float64,1},1},Nothing,Nothing,Array{Float64,1},Array{Array{Array{Float64,1},1},1},ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Tuple{Float64,Int64,Float64,Int64,Float64,Int64,Float64,Float64,Int64,Float64,Float64,Int64},FoodChain{getfield(Main, Symbol("##23#24")),Nothing,Nothing,Nothing,Nothing,Nothing,Any,Any},Nothing,DiffEqBase.StandardODEProblem},CompositeAlgorithm{Tuple{Tsit5,Rosenbrock23{0,false,LinSolveFactorize{typeof(LinearAlgebra.lu!)},DataType}},AutoSwitch{Tsit5,Rosenbrock23{0,false,LinSolveFactorize{typeof(LinearAlgebra.lu!)},DataType},Rational{Int64},Float64}},OrdinaryDiffEq.CompositeInterpolationData{FoodChain{getfield(Main, Symbol("##23#24")),Nothing,Nothing,Nothing,Nothing,Nothing,Any,Any},Array{Array{Float64,1},1},Array{Float64,1},Array{Array{Array{Float64,1},1},1},OrdinaryDiffEq.CompositeCache{Tuple{OrdinaryDiffEq.Tsit5Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}},OrdinaryDiffEq.Rosenbrock23Cache{Array{Float64,1},Array{Float64,1},Array{Float64,2},Array{Float64,2},OrdinaryDiffEq.Rosenbrock23ConstantCache{Float64,typeof(identity),typeof(identity)},DiffEqDiffTools.TimeGradientWrapper{FoodChain{getfield(Main, Symbol("##23#24")),Nothing,Nothing,Nothing,Nothing,Nothing,Any,Any},Array{Float64,1},Tuple{Float64,Int64,Float64,Int64,Float64,Int64,Float64,Float64,Int64,Float64,Float64,Int64}},DiffEqDiffTools.UJacobianWrapper{FoodChain{getfield(Main, Symbol("##23#24")),Nothing,Nothing,Nothing,Nothing,Nothing,Any,Any},Float64,Tuple{Float64,Int64,Float64,Int64,Float64,Int64,Float64,Float64,Int64,Float64,Float64,Int64}},LinSolveFactorize{typeof(LinearAlgebra.lu!)},DiffEqDiffTools.JacobianCache{Array{Float64,1},Array{Float64,1},Array{Float64,1},Val{:central},Float64,Val{true}},DiffEqDiffTools.GradientCache{Nothing,Array{Float64,1},Array{Float64,1},Val{:central},Float64,Val{true}}}},AutoSwitch{Tsit5,Rosenbrock23{0,false,LinSolveFactorize{typeof(LinearAlgebra.lu!)},DataType},Rational{Int64},Float64}}}},FoodChain{getfield(Main, Symbol("##23#24")),Nothing,Nothing,Nothing,Nothing,Nothing,Any,Any},OrdinaryDiffEq.CompositeCache{Tuple{OrdinaryDiffEq.Tsit5Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}},OrdinaryDiffEq.Rosenbrock23Cache{Array{Float64,1},Array{Float64,1},Array{Float64,2},Array{Float64,2},OrdinaryDiffEq.Rosenbrock23ConstantCache{Float64,typeof(identity),typeof(identity)},DiffEqDiffTools.TimeGradientWrapper{FoodChain{getfield(Main, Symbol("##23#24")),Nothing,Nothing,Nothing,Nothing,Nothing,Any,Any},Array{Float64,1},Tuple{Float64,Int64,Float64,Int64,Float64,Int64,Float64,Float64,Int64,Float64,Float64,Int64}},DiffEqDiffTools.UJacobianWrapper{FoodChain{getfield(Main, Symbol("##23#24")),Nothing,Nothing,Nothing,Nothing,Nothing,Any,Any},Float64,Tuple{Float64,Int64,Float64,Int64,Float64,Int64,Float64,Float64,Int64,Float64,Float64,Int64}},LinSolveFactorize{typeof(LinearAlgebra.lu!)},DiffEqDiffTools.JacobianCache{Array{Float64,1},Array{Float64,1},Array{Float64,1},Val{:central},Float64,Val{true}},DiffEqDiffTools.GradientCache{Nothing,Array{Float64,1},Array{Float64,1},Val{:central},Float64,Val{true}}}},AutoSwitch{Tsit5,Rosenbrock23{0,false,LinSolveFactorize{typeof(LinearAlgebra.lu!)},DataType},Rational{Int64},Float64}},OrdinaryDiffEq.DEOptions{Float64,Float64,Float64,Float64,typeof(DiffEqBase.ODE_DEFAULT_NORM),typeof(LinearAlgebra.opnorm),CallbackSet{Tuple{},Tuple{}},typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN),typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE),typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK),DataStructures.BinaryHeap{Float64,DataStructures.LessThan},DataStructures.BinaryHeap{Float64,DataStructures.LessThan},Nothing,Nothing,Int64,Array{Float64,1},Float64,Array{Float64,1}},Array{Float64,1},Float64}, ::OrdinaryDiffEq.Tsit5Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}}, ::Bool) at C:\Users\Helga\.julia\packages\OrdinaryDiffEq\Ih6Ud\src\perform_step\low_order_rk_perform_step.jl:614
 [7] solve!(::OrdinaryDiffEq.ODEIntegrator{CompositeAlgorithm{Tuple{Tsit5,Rosenbrock23{0,false,LinSolveFactorize{typeof(LinearAlgebra.lu!)},DataType}},AutoSwitch{Tsit5,Rosenbrock23{0,false,LinSolveFactorize{typeof(LinearAlgebra.lu!)},DataType},Rational{Int64},Float64}},Array{Float64,1},Float64,Tuple{Float64,Int64,Float64,Int64,Float64,Int64,Float64,Float64,Int64,Float64,Float64,Int64},Float64,Float64,Float64,Array{Array{Float64,1},1},OrdinaryDiffEq.ODECompositeSolution{Float64,2,Array{Array{Float64,1},1},Nothing,Nothing,Array{Float64,1},Array{Array{Array{Float64,1},1},1},ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Tuple{Float64,Int64,Float64,Int64,Float64,Int64,Float64,Float64,Int64,Float64,Float64,Int64},FoodChain{getfield(Main, Symbol("##23#24")),Nothing,Nothing,Nothing,Nothing,Nothing,Any,Any},Nothing,DiffEqBase.StandardODEProblem},CompositeAlgorithm{Tuple{Tsit5,Rosenbrock23{0,false,LinSolveFactorize{typeof(LinearAlgebra.lu!)},DataType}},AutoSwitch{Tsit5,Rosenbrock23{0,false,LinSolveFactorize{typeof(LinearAlgebra.lu!)},DataType},Rational{Int64},Float64}},OrdinaryDiffEq.CompositeInterpolationData{FoodChain{getfield(Main, Symbol("##23#24")),Nothing,Nothing,Nothing,Nothing,Nothing,Any,Any},Array{Array{Float64,1},1},Array{Float64,1},Array{Array{Array{Float64,1},1},1},OrdinaryDiffEq.CompositeCache{Tuple{OrdinaryDiffEq.Tsit5Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}},OrdinaryDiffEq.Rosenbrock23Cache{Array{Float64,1},Array{Float64,1},Array{Float64,2},Array{Float64,2},OrdinaryDiffEq.Rosenbrock23ConstantCache{Float64,typeof(identity),typeof(identity)},DiffEqDiffTools.TimeGradientWrapper{FoodChain{getfield(Main, Symbol("##23#24")),Nothing,Nothing,Nothing,Nothing,Nothing,Any,Any},Array{Float64,1},Tuple{Float64,Int64,Float64,Int64,Float64,Int64,Float64,Float64,Int64,Float64,Float64,Int64}},DiffEqDiffTools.UJacobianWrapper{FoodChain{getfield(Main, Symbol("##23#24")),Nothing,Nothing,Nothing,Nothing,Nothing,Any,Any},Float64,Tuple{Float64,Int64,Float64,Int64,Float64,Int64,Float64,Float64,Int64,Float64,Float64,Int64}},LinSolveFactorize{typeof(LinearAlgebra.lu!)},DiffEqDiffTools.JacobianCache{Array{Float64,1},Array{Float64,1},Array{Float64,1},Val{:central},Float64,Val{true}},DiffEqDiffTools.GradientCache{Nothing,Array{Float64,1},Array{Float64,1},Val{:central},Float64,Val{true}}}},AutoSwitch{Tsit5,Rosenbrock23{0,false,LinSolveFactorize{typeof(LinearAlgebra.lu!)},DataType},Rational{Int64},Float64}}}},FoodChain{getfield(Main, Symbol("##23#24")),Nothing,Nothing,Nothing,Nothing,Nothing,Any,Any},OrdinaryDiffEq.CompositeCache{Tuple{OrdinaryDiffEq.Tsit5Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}},OrdinaryDiffEq.Rosenbrock23Cache{Array{Float64,1},Array{Float64,1},Array{Float64,2},Array{Float64,2},OrdinaryDiffEq.Rosenbrock23ConstantCache{Float64,typeof(identity),typeof(identity)},DiffEqDiffTools.TimeGradientWrapper{FoodChain{getfield(Main, Symbol("##23#24")),Nothing,Nothing,Nothing,Nothing,Nothing,Any,Any},Array{Float64,1},Tuple{Float64,Int64,Float64,Int64,Float64,Int64,Float64,Float64,Int64,Float64,Float64,Int64}},DiffEqDiffTools.UJacobianWrapper{FoodChain{getfield(Main, Symbol("##23#24")),Nothing,Nothing,Nothing,Nothing,Nothing,Any,Any},Float64,Tuple{Float64,Int64,Float64,Int64,Float64,Int64,Float64,Float64,Int64,Float64,Float64,Int64}},LinSolveFactorize{typeof(LinearAlgebra.lu!)},DiffEqDiffTools.JacobianCache{Array{Float64,1},Array{Float64,1},Array{Float64,1},Val{:central},Float64,Val{true}},DiffEqDiffTools.GradientCache{Nothing,Array{Float64,1},Array{Float64,1},Val{:central},Float64,Val{true}}}},AutoSwitch{Tsit5,Rosenbrock23{0,false,LinSolveFactorize{typeof(LinearAlgebra.lu!)},DataType},Rational{Int64},Float64}},OrdinaryDiffEq.DEOptions{Float64,Float64,Float64,Float64,typeof(DiffEqBase.ODE_DEFAULT_NORM),typeof(LinearAlgebra.opnorm),CallbackSet{Tuple{},Tuple{}},typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN),typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE),typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK),DataStructures.BinaryHeap{Float64,DataStructures.LessThan},DataStructures.BinaryHeap{Float64,DataStructures.LessThan},Nothing,Nothing,Int64,Array{Float64,1},Float64,Array{Float64,1}},Array{Float64,1},Float64}) at C:\Users\Helga\.julia\packages\OrdinaryDiffEq\Ih6Ud\src\perform_step\composite_perform_step.jl:43
 [8] #__solve#202(::Base.Iterators.Pairs{Symbol,Real,Tuple{Symbol,Symbol},NamedTuple{(:default_set, :saveat),Tuple{Bool,Float64}}}, ::Function, ::ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Tuple{Float64,Int64,Float64,Int64,Float64,Int64,Float64,Float64,Int64,Float64,Float64,Int64},FoodChain{getfield(Main, Symbol("##23#24")),Nothing,Nothing,Nothing,Nothing,Nothing,Any,Any},Nothing,DiffEqBase.StandardODEProblem}, ::CompositeAlgorithm{Tuple{Tsit5,Rosenbrock23{0,false,LinSolveFactorize{typeof(LinearAlgebra.lu!)},DataType}},AutoSwitch{Tsit5,Rosenbrock23{0,false,LinSolveFactorize{typeof(LinearAlgebra.lu!)},DataType},Rational{Int64},Float64}}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at C:\Users\Helga\.julia\packages\OrdinaryDiffEq\Ih6Ud\src\solve.jl:7
 [9] (::getfield(DiffEqBase, Symbol("#kw##__solve")))(::NamedTuple{(:default_set, :saveat),Tuple{Bool,Float64}}, ::typeof(DiffEqBase.__solve), ::ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Tuple{Float64,Int64,Float64,Int64,Float64,Int64,Float64,Float64,Int64,Float64,Float64,Int64},FoodChain{getfield(Main, Symbol("##23#24")),Nothing,Nothing,Nothing,Nothing,Nothing,Any,Any},Nothing,DiffEqBase.StandardODEProblem}, ::CompositeAlgorithm{Tuple{Tsit5,Rosenbrock23{0,false,LinSolveFactorize{typeof(LinearAlgebra.lu!)},DataType}},AutoSwitch{Tsit5,Rosenbrock23{0,false,LinSolveFactorize{typeof(LinearAlgebra.lu!)},DataType},Rational{Int64},Float64}}) at .\none:0
 [10] #__solve#2(::Bool, ::Base.Iterators.Pairs{Symbol,Float64,Tuple{Symbol},NamedTuple{(:saveat,),Tuple{Float64}}}, ::Function, ::ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Tuple{Float64,Int64,Float64,Int64,Float64,Int64,Float64,Float64,Int64,Float64,Float64,Int64},FoodChain{getfield(Main, Symbol("##23#24")),Nothing,Nothing,Nothing,Nothing,Nothing,Any,Any},Nothing,DiffEqBase.StandardODEProblem}, ::Nothing) at C:\Users\Helga\.julia\packages\DifferentialEquations\88DSk\src\default_solve.jl:15
 [11] #__solve at .\none:0 [inlined]
 [12] #__solve#1 at C:\Users\Helga\.julia\packages\DifferentialEquations\88DSk\src\default_solve.jl:5 [inlined]
 [13] #__solve at .\none:0 [inlined]
 [14] #solve#442 at C:\Users\Helga\.julia\packages\DiffEqBase\nW6r3\src\solve.jl:41 [inlined]
 [15] (::getfield(DiffEqBase, Symbol("#kw##solve")))(::NamedTuple{(:saveat,),Tuple{Float64}}, ::typeof(solve), ::ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Tuple{Float64,Int64,Float64,Int64,Float64,Int64,Float64,Float64,Int64,Float64,Float64,Int64},FoodChain{getfield(Main, Symbol("##23#24")),Nothing,Nothing,Nothing,Nothing,Nothing,Any,Any},Nothing,DiffEqBase.StandardODEProblem}) at .\none:0
 [16] top-level scope at none:0
in expression starting at C:\owncloud\FoodWebsAnthropocene\Modelling\Julia\Julia foodchain code:12
2

When you look up the error by clicking on the blue line number, you will see

    if isnan(z) & !isnan(x+y)
        throw_exp_domainerror(x)
    end

in math.jl online 782. Apparently your simulation is producing not-a-numbers.

1 Like
3

That source line is deceiving and due to a mismatch in LLVM/Julia semantics. His simulation isn’t generating NaNs — it’s generating negative numbers. Exponentiating negative numbers by a fractional exponent is an error in Julia but LLVM produces NaNs.

In other words, the issue is this:

julia> (-1.534065831079686)^1.2
ERROR: DomainError with -1.534065831079686:
Exponentiation yielding a complex result requires a complex argument.
Replace x^y with (x+0im)^y, Complex(x)^y, or similar.
Stacktrace:
 [1] throw_exp_domainerror(::Float64) at ./math.jl:35
 [2] ^(::Float64, ::Float64) at ./math.jl:769
 [3] top-level scope at none:0
4

??? How does the error get misattributed? Are the line numbers in general in error?

5

No, the error is in the right place. I’m using master, so my line numbers are slightly different.

The point is that the simulation isn’t the thing directly generating NaNs. The exponentiation is. That isnan check occurs after computing the exponentiation. The exponentiation intrinsic is generating NaNs because you cannot represent the result of a negative number exponentiated by a fractional exponent in a floating point number. Julia chooses to throw an error in such a case instead of returning a NaN.

6

Thank you for your reply. I am unable to click on the blue line number - I am using Atom and the stacktrace is all grey, but when I copy code into this forum colours do appear. Is there a way to have colour in a code in Atom? It would make it a lot easier to read, and maybe it would make things clickable?

I am sorry if these are basic questions but I am just learning Julia.

7

Thank you for your reply. Do you have any ideas why it is generating a negative number? I am confused how it can, because there are no negative values for those parameters.

8

Petr is likely using a Jupyter notebook, which linkifies the line numbers in stack traces. It can be very handy when you’re debugging your own files, but in a #user:first-steps thread it’s just a distraction to dig deep into Julia’s internals.

I’m not a differential equations expert, but you’re solving for BB and BI and have both BB^h and BI^h in your equations. If either ever goes below zero, that’ll error.

9

A post was split to a new topic: Is exp implemented with an LLVM intrinsic?

10

When experiencing such behavior, it is a good idea to check your simulations for a shorter time periode. By playing around with your code, it works with tspan = (0,13.87), with result (I used the pyplot() backend):
image
… but it crashes with tspan = (0,13.88).

Based on the plot, I would guess that BI becomes negative between time 13.87 and 13.88. Thus, either your model is unstable, or the model is extremely sensitive to discretization errors.

In this case, I would have checked for typos in the model.

2 Likes
11

Thank you for your reply. Do you know of a way for Julia to state the tspan at which it will crash, or is it a guessing game?

12

I did trial and error… That is not too difficult: start with a very small number to get a first idea, and then increase until it crashes.

Your model appears to be similar to the following one, A Three-Species Food Chain System with Two Types of Functional Responses

This paper contains some stability proofs - some numerical values. [I cannot vouch for the quality of this paper – but it looks ok. I found several papers in Science Direct, but since I am outside of my job network, I cannot download those for free…)

13

Thank you, I will try to stabilize my model.