# Using DifferentialEquations.jl for a simple food-chain

Hello, I am a complete Julia beginner and am working through Example 2: Solving Systems of Equations and Defining Parameterized Functions, trying to adapt my food-chain formulas to Julia. I have also watched Intro to solving differential equations in Julia.

I am using Julia 1.0 on Windows 10, and have updated all packages and have DifferentialEquations (also updated) installed.

The start equations and parameters are:

``````   du[1] = u1*(1-u1)-(1/e)*y2*x2*u1^(1-q)/(u012^(1+q)+u1^(q+1))*u2
du[2] = y2*x2*u1^(1+q)/(u012^(1+q)+u1^(1+q))*u2-(1/e)*y3*x3*u2^(1+q)/(u023^(1+q)+u2^(1+q))*u3-x2*u2
du[3] = y3*x3*u2^(1+q)/(u023^(1+q)+u2^(1+q))*u3-x3*u3
``````

Initial values = u1= 0.01, u2 = 0.01, u3 = 0.01
e = 1
y2 = 2
x2 = 3
q = 4
u012 = 5
y3 = 6
x3 = 7
u023 = 8
t-max = 10000
t-interval = 0.05

Based on the tutorial my code currently looks like this:

using DifferentialEquations
function Bio(du,u,p,t)
du[1] = u[1](1-u[1])-(1/[1])[2]*[3]*u[1]^(1-[4])/([5]^(1+[4])u[1]^([4]+1))u[2]
du[2] = [2]
[3]u[1]^(1+[4])/([5]^(1+[4])+u[1]^(1+[4]))u[2]-(1/[1])[6][7]*u[2]^(1+[4])/([8]^(1+[4])+u[2]^(1+[4]))*u[3]-[3]u[2]
du[3] = [6]
[7]*u[2]^(1+[4])/([8]^(1+[4])+u[2]^(1+[4]))*u[3]-[7]*u[3]
end
u0 = [0.01,0.01,0.01]
p=(1,2.009,0.4,5,0.08,0,0.16129,0.5)
tspan = (0.01, 10000)
prob=ODEProblem(Bio,u0,tspan,p)
sol=solve(prob)

I generate an error once sol=solve(prob) runs.

ERROR: MethodError: no method matching -(::Int32, ::Array{Int32,1})
Closest candidates are:
-(::Union{Int128, Int16, Int32, Int64, Int8, UInt128, UInt16, UInt32, UInt64, UInt8}) at int.jl:51
-(::T<:Union{Int128, Int16, Int32, Int64, Int8, UInt128, UInt16, UInt32, UInt64, UInt8}, ::T<:Union{Int128, Int16, Int32, Int64, Int8, UInt128, UInt16, UInt32, UInt64, UInt8}) where T<:Union{Int128, Int16, Int32, Int64, Int8, UInt128, UInt16, UInt32, UInt64, UInt8} at int.jl:52
-(::Union{Int16, Int32, Int8}, ::BigInt) at gmp.jl:458
â€¦
Stacktrace:
[1] Bio(::Array{Float64,1}, ::Array{Float64,1}, ::Tuple{Int32,Float64,Float64,Int32,Float64,Int32,Float64,Float64}, ::Float64) at .\REPL[12]:2
[2] (::ODEFunction{true,typeof(Bio),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing})(::Array{Float64,1}, ::Array{Float64,1}, ::Vararg{Any,N} where N) at C:\Users\Helga.julia\packages\DiffEqBase\mcoTV\src\diffeqfunction.jl:106
[3] initialize!..

â€¦ followed by 5 pages of text.

What am I doing wrong? And how may I incorporate time intervals of 0.05? Any assistant/advice would be greatly appreciated.

In Julia 1.0, you need to use `.-` (i.e. broadcasting `-` by adding a dot) to subtract a scalar and an array. Julia 0.6 would let you do this with `-`, 0.7 added a deprecation warning, and 1.0 throws an error. Looks like the DifferentialEquations example problems havenâ€™t been updated yet. For a deeper dive on Juliaâ€™s dot-broadcasting syntax, check out this blog post.

Also, I just noticed your code has a bunch of numbers in bracketsâ€“not sure what your intention is, but each of these is going to evaluate to a single-element vector, which is probably not what you want?

`````` du[2] = [2][3]u[1]^(1+[4])/([5]^(1+[4])+u[1]^(1+[4]))u[2]-(1/[1])[6][7]*u[2]^(1+[4])/([8]^(1+[4])+u[2]^(1+[4]))*u[3]-[3]u[2]
``````

(Also, make sure you use triple backticks (```) to quote your code blocks so they format correctly.)

This error:

``````ERROR: MethodError: no method matching -(::Int32, ::Array{Int32,1})
``````

means that you should be â€śbroadcastingâ€ť one of operations, or there is a typo where you expect a scalar but are currently getting a vector. It seems to me most likely that you have a typo in your code, if you can re-quote it in backticks it may be more clear. Based on how your code currently reads, are the `[1]` and so on supposed to be `p[1]`??

I thought adjustable variables need to be replaced with `[1]` and to be specified in p, so `p=(1,2.009...)` would correspond with `[1][2]...` ect.

``````e =1				e=[1]
y2 = 2.009			y2=[2]
x2 = 0.4			x2=[3]
q = 0				q=[4]
u012 = 0.16129			u012=[5]
y3 =5				y3=[6]
x3 = 0.08			x3=[7]
u023 = 0.5			u023=[8]
``````

My code block is (no change from previous)

``````using DifferentialEquations
function Bio(du,u,p,t)
du[1] = u[1]*(1-u[1])-(1/[1])*[2]*[3]*u[1]^(1-[4])/([5]^(1+[4])u[1]^([4]+1))*u[2]
du[2] = [2]*[3]*u[1]^(1+[4])/([5]^(1+[4])+u[1]^(1+[4]))*u[2]-(1/[1])*[6]*[7]*u[2]^(1+[4])/([8]^(1+[4])+u[2]^(1+[4]))*u[3]-[3]*u[2]
du[3] = [6]*[7]*u[2]^(1+[4])/([8]^(1+[4])+u[2]^(1+[4]))*u[3]-[7]*u[3]
end
u0 = [0.01,0.01,0.01]
p=(1,2.009,0.4,5,0.08,0,0.16129,0.5)
tspan = (0.01, 10000)
prob=ODEProblem(Bio,u0,tspan,p)
sol=solve(prob)
``````

@ElOceanografo So I should replace each `-` with `.-` whenever I am subtracting a substituted number `[2]` such as

`du[1] = u[1]*(1-u[1])-(1/[1])*[2]*[3]*u[1]^(1.-[4])/([5]^(1+[4])u[1]^([4]+1))*u[2]` ?

@platawiec So I should place `p[#]` in front of each number that is an adjustable variable like

``````du[1] = u[1]*(1-u[1])-(1/p[1])*p[2]*p[3]*u[1]^(1-p[4])/(p[5]^(1+p[4])u[1]^(p[4]+1))*u[2]
``````

I am sorry if these are very basic questions, but I am still a very basic person. Thank you for your patience.

Yes, this is correct. Have you tried running the code with this change?

Just like you write `u[1]` to get the first entry of `u`, you write `p[1]` to get the first entry of `p` (in this case `p[1] = 1`, etc.)

I have put the pâ€™s in and it now generates:

``````julia> function Bio(du,u,p,t)
du[1] = u[1]*(1-u[1])-(1/p[1])*p[2]*p[3]*u[1]^(1-p[4])/p([5]^(1+p[4])u[1]^(p[4]+1))*u[2]
du[2] = p[2]*p[3]*u[1]^(1+p[4])/(p[5]^(1+p[4])+u[1]^(1+p[4]))*u[2]-(1/p[1])*p[6]*p[7]*u[2]^(1+p[4])/(p[8]^(1+p[4])+u[2]^(1+p[4]))*u[3]-p[3]*u[2]
du[3] = p[6]*p[7]*u[2]^(1+p[4])/(p[8]^(1+p[4])+u[2]^(1+p[4]))*u[3]-p[7]*u[3]
end
Bio (generic function with 1 method)

julia> u0 = [0.01,0.01,0.01]
3-element Array{Float64,1}:
0.01
0.01
0.01

julia> p=(1,2.009,0.4,5,0.08,0,0.16129,0.5)
(1, 2.009, 0.4, 5, 0.08, 0, 0.16129, 0.5)

julia> tspan = (0.01, 10000)
(0.01, 10000)

julia> prob=ODEProblem(Bio,u0,tspan,p)
ODEProblem with uType Array{Float64,1} and tType Float64. In-place: true
timespan: (0.01, 10000.0)
u0: [0.01, 0.01, 0.01]

julia> sol=solve(prob)
ERROR: MethodError: no method matching ^(::Array{Int32,1}, ::Float64)
Closest candidates are:
^(::Missing, ::Number) at missing.jl:93
^(::Float64, ::Float64) at math.jl:767
^(::Irrational{:â„Ż}, ::Number) at mathconstants.jl:91
...

``````

What is the next step?

You have one typo in there. The error message is giving you a hint where to look for it, so read it carefully.

I am joining in this late so maybe it does more harm than good, but this probably wouldâ€™ve been a good use-case for the `@ode_def` macro.

``````f = @ode_def FoodChain begin
du1 = u1*(1-u1)-(1/e)*y2*x2*u1^(1-q)/(u012^(1+q)+u1^(q+1))*u2
du2 = y2*x2*u1^(1+q)/(u012^(1+q)+u1^(1+q))*u2-(1/e)*y3*x3*u2^(1+q)/(u023^(1+q)+u2^(1+q))*u3-x2*u2
du3 = y3*x3*u2^(1+q)/(u023^(1+q)+u2^(1+q))*u3-x3*u3
end e y2 x2 q u012 y3 x3 u023
u0 = [0.01,0.01,0.01]
tspan = (0.0,10000.0)
p = (1,2,3,4,5,6,7,8)
prob = ODEProblem(f,u0,tspan,p)
solve(prob,saveat=0.05)
``````

This might be a silly question, but does the error message mean that there is an error somewhere in the formula where there is a `^(`?

Yes, somewhere you took the exponential of a vector of ints, likely `p^`. If you shared the error message it would be easier to find out though.

Oh my.

I found the error of

``````p([5]^(1
``````

instead of `(p[5]^(1`

@ChrisRackauckas the error message is

``````ERROR: MethodError: no method matching ^(::Array{Int32,1}, ::Float64)
Closest candidates are:
^(::Missing, ::Number) at missing.jl:93
^(::Float64, ::Float64) at math.jl:767
^(::Irrational{:â„Ż}, ::Number) at mathconstants.jl:91
...
``````

Unless you mean the accompanying 5 pages of text?

`retcode: Success` (Which I must say I really enjoy this about Julia, that and `Everything is okay` when installing)

I will try using `@ode_def` from now on to make things simpler.

Thank you everyone for your help!

1 Like

@ChrisRackauckas Sorry to bother you again, but when I run the `@ode_def` macro I get:

``````julia> f = @ode_def FoodChain begin
du1 = u1*(1-u1)-(1/e)*y2*x2*u1^(1+q)/(u012^(1+q)+u1^(q+1))*u2
du2 = y2*x2*u1^(1+q)/(u012^(1+q)+u1^(1+q))*u2-(1/e)*y3*x3*u2^(1+q)/(u023^(1+q)+u2^(1+q))*u3-x2*u2
du3 = y3*x3*u2^(1+q)/(u023^(1+q)+u2^(1+q))*u3-x3*u3
end e y2 x2 q u012 y3 x3 u023
â”Ś Warning: Failed to build the parameter derivatives.
â”” @ ParameterizedFunctions C:\Users\Helga\.julia\packages\ParameterizedFunctions\47E7u\src\ode_def_opts.jl:226
(::FoodChain{getfield(Main, Symbol("##3#8")),getfield(Main, Symbol("##4#9")),getfield(Main, Symbol("##5#10")),getfield(Main, Symbol("##6#11")),getfield(Main, Symbol("##7#12")),Nothing,Expr,Expr}) (generic function with 2 methods)
``````

What is the reason?

It just means that the Jacobian of the parameters doesnâ€™t seem to be buildable, not sure why it isnâ€™t, but that shouldnâ€™t hurt your use case? You can use `@ode_def_bare` to turn off all extra computations.