Using DifferentialEquations.jl for a simple food-chain

New to Julia
,
1

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.

2

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

3

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]??

4

@ElOceanografo @platawiec

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.

5

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

6

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?

7

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

8

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

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 ^(?

10

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.

11

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
12

@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?

13

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.