# How to use DifferentialEquations with complicated equations

**URL:** <https://discourse.julialang.org/t/how-to-use-differentialequations-with-complicated-equations/20740>\
**Category:** New to Julia\
**Tags:** diffeq\
**Created:** [February 13, 2019, 10:34am UTC](https://discourse.julialang.org/t/how-to-use-differentialequations-with-complicated-equations/20740 "2019-02-13T10:34:28Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![HelgavonLichtenstein](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/helgavonlichtenstein/32/5005_2.png) [@HelgavonLichtenstein](https://discourse.julialang.org/u/HelgavonLichtenstein)\
**Post date:** [February 13, 2019, 10:34am UTC](https://discourse.julialang.org/t/how-to-use-differentialequations-with-complicated-equations/20740/1 "2019-02-13T10:34:28Z")

</div>

Hello,

I understand the basic structure of DifferentialEquations but what if the equations you would like to use do not fit neatly like

```julia
function lorenz(du,u,p,t)
 du[1] = 10.0*(u[2]-u[1])
 du[2] = u[1]*(28.0-u[3]) - u[2]
 du[3] = u[1]*u[2] - (8/3)*u[3]
end

```

For example, the equations I need to use (1,2) could be written to fit within DifferentialEquations structure, but I would be unable to list the parameters (3,4 ect) neatly since they refer to other variables that (I assume) would need to be included within the DifferentialEquations structure. Maybe they need to listed as global constants and can be referred to without being in the DifferentialEquations structure?

![image](https://global.discourse-cdn.com/julialang/original/3X/c/9/c9d4892358a4c43f701e6abda6575eb593037920.png)

Does anyone know how this could be done?

```julia
using Distributions
#Constants
const Ropt = 100
const γ= 2
const k = 8.6173324e-5 #Boltzman eV/K
const T0 = 273.15 #K
const Temp = 293.15 # example of 20 C degrees
const q = 1.2 #Hill exponent
const NumSpecies = 50
const B = fill(430.0,50)

const m = sort(rand(Uniform(10^2, 10^12), NumSpecies))

#L is the success curve of consumers aka feeding efficiency, a probability matrix
L = @. (m/(m'*Ropt) * exp(1-m/(m'*Ropt)))^γ # The @. macro here broadcasts the following scalar expression across `m` and `m'`
# Weak links (Lij <p) are removed
L[L .<= 0.01] .= 0 # . lets you do things element by element in an array
L

#create a matrix of links:
#take the link probabilities of the previous matrix and determine if L[i,j] > p
FoodWeb = zero(L)
FoodWeb[L .> rand.(Uniform(0, 1))] .= 1 # Sets F[i,j] = 0 wherever L[i,j] < F[i,j]
FoodWeb

function ParameterMatrix(m, I, b, c, E, FoodWeb = FoodWeb)
    #FoodWeb (to match food web generated), m (bodymass), I (intercept), b (body scaling species i)
  #c (body mass scaling Pred), E (Activation energy)
    x = @. (exp(I)*(m^b)*(m'^c)*exp(E*((T0 - Temp)/(k*T0*Temp)))).*FoodWeb

end

#attack rate a (m^2/s) - corresponds to formula 4
a = ParameterMatrix(m, -13.1, 0.25, -0.8, -0.38)

#handling time HT (s) (TH in Binzer paper)
HT = ParameterMatrix(m,9.66, -0.45, 0.47, 0.26)

function ParameterVector(m, I, b, E, FoodWeb = FoodWeb)
    x = @.(exp(I)*(m^b)*exp(E*((T0 - Temp)/(k*T0*Temp)))).*FoodWeb
end
#Carrying capacity K (g/m2)
K = ParameterVector(m, 10, 0.28, 0.71)
#growth rate r (1/s)
r = ParameterVector(m, -15.68, -0.25, -0.84)
#metabolism x ([1/s])
x = ParameterVector(m, -16.54, -0.31, -0.69)
#Feeding Function F

function Feeding(a, B, q, HT) #corresponds to formula 3
    x = (a.*(B').^q)/(1+ sum((HT[:,j:j]).*(a[:,j:j]).*B.^q))
end

dBasal = function Basal(r,B,K,Feeding) #corresponds to formula 1
  for i in 1:NumSpecies #part of the equation corresponding to their growth
    dBasal[i] = r[i:i,:].*B'*(1-B'./K[:,j:j])
    for j in 1:NumSpecies # part of the equation corresponding to the sum (losses by predation)
      dBasal[j] = dBasal[j] - Feeding(j,i, a, B, q, HT).*B'
      end
    end
end

dConsumer = function Consumer(e,Feeding,B,x) #corresponds to formula 2
  for i in 1:NumSpecies #part of the equation corresponding to their growth
    dConsumer[i] = dConsumer[i:i,:] + e.*Feeding(i, j, a, B[j], q, HT).*B[i:i,:]
    for j in 1:NumSpecies # part of the equation corresponding to the sum (losses by predation)
      dConsumer[j] = dConsumer[:,j:j] - Feeding(j,i, a, B, q, HT).*B'
      end
    end
    dConsumer[i] - x.*B[i:i,:]
end

```

---

<div class="post-metadata">

**Author:** ![mauro3](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mauro3/32/292_2.png) [@mauro3](https://discourse.julialang.org/u/mauro3)\
**Post date:** [February 13, 2019, 10:54am UTC](https://discourse.julialang.org/t/how-to-use-differentialequations-with-complicated-equations/20740/2 "2019-02-13T10:54:35Z")

</div>

I’m not sure where the confusion is. The function passed to diffeq can do anything it likes, including calculating intermediate things like `F_ij`. You can pass in a `p`, which can be any type you like, to hold your parameters. Maybe looking at [GitHub - SciML/SciMLTutorials.jl: Tutorials for doing scientific machine learning (SciML) and high-performance differential equation solving with open source software.](https://github.com/JuliaDiffEq/DiffEqTutorials.jl) could help you?

Also, your equations 1&2 are both for `dBi/dt`, this means you can eliminate `dBi/dt` and end up with a system of algebraic equations which you’d solve with, e.g., [GitHub - JuliaNLSolvers/NLsolve.jl: Julia solvers for systems of nonlinear equations and mixed complementarity problems](https://github.com/JuliaNLSolvers/NLsolve.jl).

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [February 13, 2019, 2:38pm UTC](https://discourse.julialang.org/t/how-to-use-differentialequations-with-complicated-equations/20740/3 "2019-02-13T14:38:33Z")

</div>

> [@HelgavonLichtenstein](#):
>
> For example, the equations I need to use (1,2) could be written to fit within DifferentialEquations structure, but I would be unable to list the parameters (3,4 ect) neatly since they refer to other variables

You can always pass additional parameters to a functional interface by using lexical scoping. (This seems to be a FAQ for libraries that do integration, optimization, root-finding, etcetera.)

For example, the [QuadGK package](https://github.com/JuliaMath/QuadGK.jl) package lets you integrate a function of the form `f(x)` via `quadgk(x, a, b)`. But suppose you have a function `g(x, α, β)` that depends on two additional parameters `α,β` — does that mean you can’t use it with QuadGK, unless you rewrite it to use global variables? No! You would just call `quadgk(x -> g(x, α, β), a, b)`, constructing an anonymous function `x -> ...` that [“captures” the values](https://en.wikipedia.org/wiki/Closure_(computer_programming)) of α, β from e.g. local variables.

Because of this, I think DifferentialEquations didn’t really need an explicit `p` argument `f(x, p, t)` in its basic ODE interface, since parameters could have been passed using closures. My guess is that @ChrisRackauckas defined it this way to make it easier to support parameter estimation and adjoint/sensitivity analysis.

If what you really mean is that you have ODEs coupled to other systems of non-ODE equations that are expressed implicitly, this is known as a [differential algebraic equation (DAE)](https://en.wikipedia.org/wiki/Differential-algebraic_system_of_equations) and is also [supported by DifferentialEquations.jl](http://docs.juliadiffeq.org/latest/tutorials/dae_example.html).

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [February 13, 2019, 2:46pm UTC](https://discourse.julialang.org/t/how-to-use-differentialequations-with-complicated-equations/20740/4 "2019-02-13T14:46:24Z")

</div>

> [@stevengj](#):
>
> Because of this, I think DifferentialEquations didn’t really need an explicit `p` argument `f(x, p, t)` in its basic ODE interface, since parameters could have been passed using closures. My guess is that @ChrisRackauckas defined it this way to make it easier to support parameter estimation and adjoint/sensitivity analysis.

Yeah, that’s precisely it. We put it in there in order to formalize parameter handling for estimation and sensitivities, since in that case you need it. Also, a lot of people had this question for how to add parameters in, so it just seemed like a good idea to make parameters explicit (we’ve had a lot fewer questions after that, and most models do use parameters). But in reality, you can just use closures just fine as long as you aren’t estimating/calculating sensitivities.

---

<div class="post-metadata">

**Author:** ![HelgavonLichtenstein](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/helgavonlichtenstein/32/5005_2.png) [@HelgavonLichtenstein](https://discourse.julialang.org/u/HelgavonLichtenstein)\
**Post date:** [February 20, 2019, 2:12pm UTC](https://discourse.julialang.org/t/how-to-use-differentialequations-with-complicated-equations/20740/5 "2019-02-20T14:12:45Z")

</div>

Hello again,  
I have tried for several days to understand your responses but I feel I am not closer to understanding what is meant by your advice.

In the hope that I can calculate differential equations with more complicated equations, my code snippet looks like this (full code is at the end).

```julia
using DifferentialEquations
using ParameterizedFunctions # I had to install this package because I had an error that @ode_def was undefined?
using Plots
f = @ode_def_bare BiomassDensity begin
  dBasal = function Basal(r,B,K,Feeding)
    for i in 1:NumSpecies #part of the equation corresponding to their growth
      dBasal[i] = r[i:i,:].*B'*(1-B'./K[:,j:j])
      for j in 1:NumSpecies # part of the equation corresponding to the sum (losses by predation)
        dBasal[j] = dBasal[j] - Feeding(j,i, a, B, q, HT).*B'
        end
      end
  end
  dConsumer = function Consumer(e,Feeding,B,x)
    for i in 1:NumSpecies #part of the equation corresponding to their growth
      dConsumer[i] = dConsumer[i:i,:] + e.*Feeding(i, j, a, B[j], q, HT).*B[i:i,:]
      for j in 1:NumSpecies # part of the equation corresponding to the sum (losses by predation)
        dConsumer[j] = dConsumer[:,j:j] - Feeding(j,i, a, B, q, HT).*B'
        end
      end
      dConsumer[i] - x.*B[i:i,:]
  end
end

u0 = [450.0,450.0]
tspan = (0.0,1000.0)
prob = ODEProblem(f,u0,tspan)
sol=solve(prob,Rosenbrock23())

plot(sol,xlabel = "Time" ,ylabel = "Density", title = "Biomass Density", lw=0.5, layout = (2,1))

```

With this I am hoping that because `Feeding` and the different matrices/vectors of `a, B,r,K,HT` ect. are written before my attempt at differential equations, Julia will be able to call them and I won’t have to write them out within the differential equation function.

Anyway, when I run this code I get an error of `ERROR: syntax: invalid function name "internal_var___u[1]"` which possibly could be from `dBasal/dConsumer` being too long a name?

I also get the error `ERROR: UndefVarError: f not defined` when Julia reaches `prob = ODEProblem(f,u0,tspan)`, which I don’t understand as its defined at the beginning with `f = @ode_def_bare BiomassDensity begin`.

So my whole code is then (I feel weird posting this all the time but I understand it is needed for people to understand and help me - if there is some code posting etiquette please let me know):

```julia
using Distributions
using Plots
using DifferentialEquations
#Constants
const Ropt = 100
const γ= 2
const k = 8.6173324e-5 #Boltzman eV/K
const T0 = 273.15 #K
const Temp = 293.15 # example of 20 C degrees
const q = 1.2 #Hill exponent
const NumSpecies = 50
const B = fill(430.0,50)

const m = sort(rand(Uniform(10^2, 10^12), NumSpecies))

#L is the success curve of consumers aka feeding efficiency, a probability matrix
L = @. (m/(m'*Ropt) * exp(1-m/(m'*Ropt)))^γ # The @. macro here broadcasts the following scalar expression across `m` and `m'`
# Weak links (Lij <p) are removed
L[L .<= 0.01] .= 0 # . lets you do things element by element in an array
L

#create a matrix of links:
#take the link probabilities of the previous matrix and determine if L[i,j] > p
FoodWeb = zero(L)
FoodWeb[L .> rand.(Uniform(0, 1))] .= 1 # Sets F[i,j] = 0 wherever L[i,j] < F[i,j]
FoodWeb

function ParameterMatrix(m, I, b, c, E, FoodWeb = FoodWeb)
    #FoodWeb (to match food web generated), m (bodymass), I (intercept), b (body scaling species i)
  #c (body mass scaling Pred), E (Activation energy)
    x = @. (exp(I)*(m^b)*(m'^c)*exp(E*((T0 - Temp)/(k*T0*Temp)))).*FoodWeb
end

#attack rate a (m^2/s)
a = ParameterMatrix(m, -13.1, 0.25, -0.8, -0.38)

#handling time HT (s) (TH in Binzer paper)
HT = ParameterMatrix(m,9.66, -0.45, 0.47, 0.26)

function ParameterVector(m, I, b, E, FoodWeb = FoodWeb)
    x = @.(exp(I)*(m^b)*exp(E*((T0 - Temp)/(k*T0*Temp)))).*FoodWeb
end
#Carrying capacity K (g/m2)
K = ParameterVector(m, 10, 0.28, 0.71)
#growth rate r (1/s)
r = ParameterVector(m, -15.68, -0.25, -0.84)
#metabolism x ([1/s])
x = ParameterVector(m, -16.54, -0.31, -0.69)
#Feeding Function F

function Feeding(a, B, q, HT)
    x = (a.*(B').^q)/(1+ sum((HT[:,j:j]).*(a[:,j:j]).*B.^q))
end

f = @ode_def_bare BiomassDensity begin
  dBasal = function Basal(r,B,K,Feeding)
    for i in 1:NumSpecies #part of the equation corresponding to their growth
      dBasal[i] = r[i:i,:].*B'*(1-B'./K[:,j:j])
      for j in 1:NumSpecies # part of the equation corresponding to the sum (losses by predation)
        dBasal[j] = dBasal[j] - Feeding(j,i, a, B, q, HT).*B'
        end
      end
  end
  dConsumer = function Consumer(e,Feeding,B,x)
    for i in 1:NumSpecies #part of the equation corresponding to their growth
      dConsumer[i] = dConsumer[i:i,:] + e.*Feeding(i, j, a, B[j], q, HT).*B[i:i,:]
      for j in 1:NumSpecies # part of the equation corresponding to the sum (losses by predation)
        dConsumer[j] = dConsumer[:,j:j] - Feeding(j,i, a, B, q, HT).*B'
        end
      end
      dConsumer[i] - x.*B[i:i,:]
  end
end

u0 = [450.0,450.0]
tspan = (0.0,1000.0)
prob = ODEProblem(f,u0,tspan)
sol=solve(prob,Rosenbrock23())

plot(sol,xlabel = "Time" ,ylabel = "Density", title = "Biomass Density", lw=0.5, layout = (2,1))

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [February 20, 2019, 2:27pm UTC](https://discourse.julialang.org/t/how-to-use-differentialequations-with-complicated-equations/20740/6 "2019-02-20T14:27:38Z")

</div>

Don’t use the macro here. Just define a function. This part of the tutorial shows how to do this:

[http://docs.juliadiffeq.org/latest/tutorials/ode\_example.html#Example-2:-Solving-Systems-of-Equations-1](http://docs.juliadiffeq.org/latest/tutorials/ode_example.html#Example-2:-Solving-Systems-of-Equations-1)

Or maybe some of the examples here would be good to follow:

[https://nbviewer.jupyter.org/github/JuliaDiffEq/DiffEqTutorials.jl/blob/master/Introduction/OptimizingDiffEqCode.ipynb](https://nbviewer.jupyter.org/github/JuliaDiffEq/DiffEqTutorials.jl/blob/master/Introduction/OptimizingDiffEqCode.ipynb)

---

<div class="post-metadata">

**Author:** ![HelgavonLichtenstein](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/helgavonlichtenstein/32/5005_2.png) [@HelgavonLichtenstein](https://discourse.julialang.org/u/HelgavonLichtenstein)\
**Post date:** [February 20, 2019, 3:08pm UTC](https://discourse.julialang.org/t/how-to-use-differentialequations-with-complicated-equations/20740/7 "2019-02-20T15:08:31Z")

</div>

Thank you for your reply.

I’ve changed the code snippet, but I get a method error at `du[1] = function Basal(r,B,K,Feeding)`. I am still unclear on whether I can have a function within a differential equation function.

```julia
ERROR: MethodError: Cannot `convert` an object of type getfield(Main, Symbol("#Basal#39")){Array{Float64,1}} to an object of type Float64
Closest candidates are:
  convert(::Type{Float64}, ::SymEngine.BasicType{Val{:RealDouble}}) at C:\Users\Helga\.julia\packages\SymEngine\6KyFJ\src\numerics.jl:35
  convert(::Type{Float64}, ::SymEngine.Basic) at C:\Users\Helga\.julia\packages\SymEngine\6KyFJ\src\numerics.jl:160
  convert(::Type{T<:Real}, ::SymEngine.Basic) where T<:Real at C:\Users\Helga\.julia\packages\SymEngine\6KyFJ\src\numerics.jl:165

```

```julia
function BiomassDensity(du,u,p,t)
  du[1] = function Basal(r,B,K,Feeding)
    for i in 1:NumSpecies #part of the equation corresponding to their growth
      u[1][i] = @. r[i:i,:]*B'*(1-B'/K[:,j:j])
      for j in 1:NumSpecies # part of the equation corresponding to the sum (losses by predation)
        u[1][j] = @. u[1][j] - Feeding(j,i, a, B, q, HT)*B'
        end
      end
  end
  du[2] = function Consumer(e,Feeding,B,x)
    for i in 1:NumSpecies #part of the equation corresponding to their growth
      u[2][i] = @. u[2][i:i,:] + e*Feeding(i, j, a, B[j], q, HT)*B[i:i,:]
      for j in 1:NumSpecies # part of the equation corresponding to the sum (losses by predation)
        u[2][j] = @. u[2][:,j:j] - Feeding(j,i, a, B, q, HT)*B'
        end
      end
      u[2][i] = @. u[2][i] - x*B[i:i,:]
  end
end

u0 = [450.0,450.0]
tspan = (0.0,1000.0)
prob = ODEProblem(BiomassDensity,u0,tspan)
sol=solve(prob,Rosenbrock23())

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [February 20, 2019, 3:09pm UTC](https://discourse.julialang.org/t/how-to-use-differentialequations-with-complicated-equations/20740/8 "2019-02-20T15:09:51Z")

</div>

I have no idea what you’re doing with those inner anonymous functions. You may want to follow the tutorials to see the proper syntax.

---

<div class="post-metadata">

**Author:** ![ElOceanografo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/eloceanografo/32/624_2.png) [@ElOceanografo](https://discourse.julialang.org/u/ElOceanografo)\
**Post date:** [February 20, 2019, 9:12pm UTC](https://discourse.julialang.org/t/how-to-use-differentialequations-with-complicated-equations/20740/9 "2019-02-20T21:12:07Z")

</div>

It seems like a lot of your confusion is related to defining and calling functions, not from DifferentialEquations. In your last example, you are assigning a function object (a.k.a. an anonymous function) to the elements of the vector `du`, which is not what you want. I’m guessing what you want to do is to _call_ those functions (`Basal` and `Consumer`) and assign the values they return to `du[1]` and `du[2]`. I think you should work through the functions section in the manual and make sure you have a good grasp of defining, calling, and composing functions before trying to make them work in DifferentialEquations.

When you come back to this problem, here’s some code that may help you get started, based on [this Lotka-Volterra system](https://en.wikipedia.org/wiki/Competitive_Lotka%E2%80%93Volterra_equations#4-dimensional_example) from Wikipedia. Notice how I have a helper function `interaction` which is defined outside, but then called within, the main simulation function:

```julia
using DifferentialEquations
using Plots
using LinearAlgebra

function interaction(B, i, p)
    return 1 - dot(B, p[:A][i, :]) / p[:K][i]
end

function lotka_volterra!(dB, B, par, t)
    for i in 1:length(B)
        dB[i] = par[:r][i] * B[i] * interaction(B, i, par)
    end
end

A = [1 1.09 1.52 0;
     0 1 0.44 1.36;
     2.33 0 1 0.47;
     1.21 0.51 0.35 1]
r = [1, 0.72, 1.53, 1.27]
K = ones(4)

params = Dict(:A => A, :r => r, :K => K)
B0 = [0.1, 0.5, 0.05, 0.4]
tspan = (0.0, 1000.0)
prob = ODEProblem(lotka_volterra!, B0, tspan, params)
sol = solve(prob)
plot(sol)

```

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [February 20, 2019, 11:33pm UTC](https://discourse.julialang.org/t/how-to-use-differentialequations-with-complicated-equations/20740/10 "2019-02-20T23:33:04Z")

</div>

Perhaps it is somewhat difficult to give you a straight answer to your questions because it is not totally clear what you try to do/what your models do. At least not to me. You present 4 equations. Equations (1) and (2) seem to be balance equations for some quantity B\_i.

- Is index i in B\_i the species type? E.g., if you have two species, could B\_1 be “hare” and B\_2 be “fox”? In your simulation code, you indicate 50 species. Does that mean B\_1, …, B\_{50}?
- Are Eqs. (1) and (2) two possible choices of models, or will they occur simultaneously?
- Are r\_i and K\_i in Eq. (1) _constants_? What about e\_i and x\_i in Eq. (2)?
- In your Eqs. (1) and (2), when you write \sum\_i F\_{ij}B\_i, do you mean that, or do you mean \sum\_j F\_{ij} B\_j. Likewise, in Eq. (2), do you mean \sum\_j e\_j F\_{ij}B\_i or do you mean \sum\_j e\_j F\_{ij} B\_j?

Your Eq. (3) seems to be some sort of “Michaelis-Menten” saturation function (or: rather a Hill model). Your code indicates that q is a Hill exponent.

- Is H a constant? Independent of i/j?
- Is T\_i some “temperature”?
- What is a\_i in the numerator of the expression for F\_{ij}? A constant? Or related to a\_{ij} in Eq. (4)?
- What is B in the numerator of the expression for F\_{ij}? Any relationship to B\_i?
- What is the sum index in the numerator of the expression for F\_{ij}? Do you mean \sum\_i?

In your Eq. (4), this looks like a “reaction expression” with an Arrhenius temperature dependence. I assume E is “activation energy” – not per mol, but per individual, and that E/k has the unit of temperature. The term \frac{T\_0 - T}{T T\_0} = \frac{1}{T} - \frac{1}{T\_0} supports the impression that this is an Arrhenius temperature dependence. Still:

- Is d a “reaction constant”, or is it the differential operator so that dm\_i is the differential of m\_i? Probably the first, but using d always creates uncertainty.

OK – my questions reflect uncertainty in the problem that you consider. I assume this is some on-going research you do, so perhaps you don’t want to reveal every detail.

Now, looking at your code, I assume that `Feeding` returns F\_{ij} – although I don’t understand why you introduce identifier `x` within that function (is it related to x\_i in Eq. (1)).

Also, I don’t understand what functions `Basal` and `Consumer` are meant to do, related to your Eqs. (1)-(2), and what is `FoodWeb`, etc.

OK – it would be easier to help if the problem was described somewhat clearer. It will still be possible to help without understanding your model, but it is necessary to understand what you want to achieve.

---

<div class="post-metadata">

**Author:** ![HelgavonLichtenstein](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/helgavonlichtenstein/32/5005_2.png) [@HelgavonLichtenstein](https://discourse.julialang.org/u/HelgavonLichtenstein)\
**Post date:** [May 9, 2019, 9:59am UTC](https://discourse.julialang.org/t/how-to-use-differentialequations-with-complicated-equations/20740/11 "2019-05-09T09:59:31Z")

</div>

Hello again, I re-did many things with the code and chose a different paper to base my formulas on. I want to thank everyone who commented so far - this is a learning process and hopefully I did not cause too many headaches.

Based on @ElOceanografo 's example, a snippet of code now looks like:

```julia
using Distributions, DifferentialEquations, Random
const D = 0.25

S = zeros(2) .- 1
while (S[1] <= 0.0 || S[2] <= 0.0)
    S[1] = rand(Normal(10,2)) #uniform from (0, 1] i.e. 0 is excluded
    S[2] = rand(Normal(10,2))
end

const N1 = rand(Uniform(S[1]/2,S[1]))
const N2 = rand(Uniform(S[2]/2,S[2]))

NumPlantSpecies = 8
NumAnimalSpecies = 12

K = rand(Uniform(0.1,0.2), 2, NumPlantSpecies)

const v1 = 1
const v2 = 0.5

r = @. BM[1:NumPlantSpecies] ^ (-0.25)

#initial biomass of a Plant species with random values uniformly distributed between 0 (ex.) and 10 (in.)
BioP = rand(Uniform(0,10), NumPlantSpecies)
#initial biomass of an Animalspecies with random values uniformly distributed between 0 (ex.) and 10 (in.)
BioA = rand(Uniform(0,10), NumAnimalSpecies)
BioS = vcat(BioP,BioA)

CalcG = function(K, N1, N2, i)
    return (min(N1/(K[1,i] + N1), N2/(K[2,i] + N2)))
end

function dNdT(dN, D, S, N1, N2, K, v1, v2, r, BioS)
    dN[1] = D*(S[1] - N1)
    dN[2] = D*(S[2] - N2)
    for i = 1:NumPlantSpecies
        dN[1] = dN[1] - v1 * r[i] * CalcG(K, N1, N2, i) * BioS[i]
        dN[2] = dN[2] - v2 * r[i] * CalcG(K, N1, N2, i) * BioS[i]
    end
end

u0 = [BioS]
tspan = (0.0,1000.0)
p = (D, S, N1, N2, K, v1, v2, r, BioS)
# K r yIB yTI eIB eTI B0 h xI xT q E
prob = ODEProblem(dNdT, u0, tspan, p)
sol=solve(prob,Rosenbrock23())

```

I _think_ this is the right set up now. I am generating a strange error though, but it only appears when I run this part of the code. `ERROR: MethodError: no method matching zero(::Type{Array{Float64,1}})`

Could someone please confirm that this is the right setup for DifferentialEquations (i.e. does it _look_ like it should work or am I still not understanding something?)? Does anyone have an idea for the error?

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [May 9, 2019, 10:33am UTC](https://discourse.julialang.org/t/how-to-use-differentialequations-with-complicated-equations/20740/12 "2019-05-09T10:33:31Z")

</div>

Having a self-contained example would make it much easier to reproduce the error and help you.

> [@Please read: make it easier to help you](https://discourse.julialang.org/t/psa-make-it-easier-to-help-you/14757):
>
> Welcome to the Julia Discourse! We are enthusiastic about helping Julia programmers, both beginner and experienced. This public service announcement (PSA) outlines best practices when asking for help. Following these points makes it easier for us to help you and more likely you’ll get a prompt, useful answer. Keywords are highlighted to make it easier to refer to specific points. Choose a descriptive title that captures the key part of your question, eg “plots with multiple axes” instead of …

---

<div class="post-metadata">

**Author:** ![HelgavonLichtenstein](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/helgavonlichtenstein/32/5005_2.png) [@HelgavonLichtenstein](https://discourse.julialang.org/u/HelgavonLichtenstein)\
**Post date:** [May 9, 2019, 12:46pm UTC](https://discourse.julialang.org/t/how-to-use-differentialequations-with-complicated-equations/20740/13 "2019-05-09T12:46:17Z")

</div>

Thank you for your reply, I have edited my request creating what I hope is a minimal working example.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [May 9, 2019, 1:23pm UTC](https://discourse.julialang.org/t/how-to-use-differentialequations-with-complicated-equations/20740/14 "2019-05-09T13:23:17Z")

</div>

> [@HelgavonLichtenstein](#):
>
> `for i = 1:NumPlantSpecies`

[Don’t use global variables](https://docs.julialang.org/en/latest/manual/performance-tips/#Avoid-global-variables-1) in performance-critical code (e.g. the rhs of an ODE). Either make this variable `const`, pass it in as one of the parameters, or (best) simply do `for i in eachindex(r)` or `for r in r` since `r` already knows its length. You are also computing `r[i] * CalcG(K, N1, N2, i) * BioS[i]` twice each loop iteration where you could be computing it only once. You also might consider using StaticArrays.jl for things like your ODE state, since that should be much more efficient than 2-component arrays.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 9, 2019, 1:25pm UTC](https://discourse.julialang.org/t/how-to-use-differentialequations-with-complicated-equations/20740/15 "2019-05-09T13:25:30Z")

</div>

> [@HelgavonLichtenstein](#):
>
> dNdT(dN, D, S, N1, N2, K, v1, v2, r, BioS)

This is not an acceptable function signature for diffeq. You may want to check out this tutorial for more information on defining a differential equation and optimizing the code: [http://juliadiffeq.org/DiffEqTutorials.jl/html/introduction/optimizing\_diffeq\_code.html](http://juliadiffeq.org/DiffEqTutorials.jl/html/introduction/optimizing_diffeq_code.html)

---

<div class="post-metadata">

**Author:** ![ElOceanografo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/eloceanografo/32/624_2.png) [@ElOceanografo](https://discourse.julialang.org/u/ElOceanografo)\
**Post date:** [May 10, 2019, 1:06am UTC](https://discourse.julialang.org/t/how-to-use-differentialequations-with-complicated-equations/20740/17 "2019-05-10T01:06:34Z")

</div>

Let’s forget optimizing for a second and focus on getting it to run. There a a number of issues with this code, and I’ll second @ChrisRackauckas in encouraging you to work through some of the tutorials to get more familiar with how DifferentialEquations works. A few specific suggestions, though:

1. Your notation and variable names are confusing (at least to me). For instance, your ODE function is called `dNdT`, which implies that the letter `N` represents biomass or population. However, you define your biomass vector as `BioS` above it, and then create a new variable pointing to the same vector called u0 below it. Additionally, you have two global constants called `N1` and `N2`. I have a hunch some of your issues may be related to, or worsened by, unclear variable names.
2. You are confusing **state variables** and **parameters**. State variables are the quanitities you are simulating; they change over the course of the model run. For DifferentialEquations they have to be in a single vector (this would be `BioS`, `u`, or whatever you decide to call it). Paramters are the model settings and constants which do not change during the simulation (`D`, `K`, `N1`, etc.). Here you have put some of both in the tuple `p = (D, S, N1, N2, K, v1, v2, r, BioS)`.
  - Side note: because DifferentialEquations takes the parameters as an argument to `ODEProblem`, you don’t actually need to define them as global constants.

3. Your main function `dNdt` doesn’t make sense to me. I think `dN` is supposed to be `NumPlantSpecies + NumAnimalSpecies` elements long…but you only update the first two elements, and the loop only goes over `1:NumPlantSpecies`.
4. `K` and `r` are only defined for the number of plant species. Is that what you want?
5. `r` is based on a variable called `BM` which is not in this code snippet. Paste your code into a freshly started Julia REPL to make sure it runs as-is. Without a MWE, it is really difficult for us to figure out what your code is trying to do, or what isn’t working.
6. You still haven’t explained exactly what you are trying to model. Is there a paper or reference where the equations you quoted above came from? Knowing that may help us to get you started.

Again, I’d really encourage you to work through a few tutorials. Make sure you understand what the equations mean before you try to transcribe them into code. Make sure your function `dNdT` works and produces the expected results when you call it on its own, before you define an `ODEProblem`. Start simple, make your model more complex in incremental steps, and only add new features or dynamics once you’ve got the simpler version working.

---

<div class="post-metadata">

**Author:** ![nilshg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nilshg/32/2283_2.png) [@nilshg](https://discourse.julialang.org/u/nilshg)\
**Post date:** [May 10, 2019, 9:33am UTC](https://discourse.julialang.org/t/how-to-use-differentialequations-with-complicated-equations/20740/18 "2019-05-10T09:33:20Z")

</div>

> [@HelgavonLichtenstein](#):
>
> ERROR: MethodError: no method matching zero(::Type{Array{Float64,1}})

I don’t think the code you posted actually contains the line that is causing your error. It appears your error comes from a confusion of the functions `zero` and `zeros` (I noted something similar when I commented on your related issue last night [here](https://discourse.julialang.org/t/very-slow-array-creation-and-benchmark-iterate-error/24020/5))

Consider:

```julia
help?> zero
search: zero zeros iszero set_zero_subnormals get_zero_subnormals count_zeros

  zero(x)

  Get the additive identity element for the type of x (x can also specify the
  type itself).

  Examples
  ≡≡≡≡≡≡≡≡≡≡

  julia> zero(1)
  0

  julia> zero(big"2.0")
  0.0

  julia> zero(rand(2,2))
  2×2 Array{Float64,2}:
   0.0 0.0
   0.0 0.0

```

Which is different from:

```julia
help?> zeros
search: zeros count_zeros set_zero_subnormals get_zero_subnormals leading_zeros

  zeros([T=Float64,] dims...)

  Create an Array, with element type T, of all zeros with size specified by
  dims. See also fill, ones.

  Examples
  ≡≡≡≡≡≡≡≡≡≡

  julia> zeros(1)
  1-element Array{Float64,1}:
   0.0

  julia> zeros(Int8, 2, 3)
  2×3 Array{Int8,2}:
   0 0 0
   0 0 0

```

Your use case seems largely to be the creation of arrays filled with zeros, for which `zeros` is the appropriate function. Your error likely comes from calling:

```julia
julia> zero(Array{Float64,1})
ERROR: MethodError: no method matching zero(::Type{Array{Float64,1}})

```

i.e. calling the `zero` function with a Type, rather than in instance/object of a type.

More generally I think others have already pointed out that many of your issues and questions fundamentally boil down to relatively basic misunderstandings of usage and behaviour of variables, functions, values, types etc. in Julia. I think it might be helpful for you to work through an introductory resource that clarifies some of these questions and hopefully helps you to avoid issues you seem to keep running into; I always found Ben Lauwen’s [Think Julia](https://benlauwens.github.io/ThinkJulia.jl/latest/book.html) an excellent resource for this tasks and I’d recommend you take a couple of days to work through it step-by-step!

---

<div class="post-metadata">

**Author:** ![HelgavonLichtenstein](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/helgavonlichtenstein/32/5005_2.png) [@HelgavonLichtenstein](https://discourse.julialang.org/u/HelgavonLichtenstein)\
**Post date:** [May 29, 2019, 12:41pm UTC](https://discourse.julialang.org/t/how-to-use-differentialequations-with-complicated-equations/20740/19 "2019-05-29T12:41:12Z")

</div>

Hello Everyone, I am back with a functioning code that generates differential equations. Wooo! 😌

Apologies if any of this seems obvious, but the things I learned along the way (in case anyone else struggling is reading this):

Making something a constant (`const`) even if the value changes, is okay (and recommended). The const only keeps the type the same e.g. `zxcv = 12.45` can change to `zxcv = 15.23` but it cannot change to `zxcv = 12` (this is a different data type).

`zeros(size(qwer))` is used when it is an **array**  
`zeros(asdf)` is used when it is a _tuple_

`for i in length(qwer)` is faster when it is an **array**  
`for i in 1:asdf` is faster when it is a _tuple_

So if you have formulas like [this](https://www.nature.com/articles/ncomms12718);

![image](https://global.discourse-cdn.com/julialang/original/3X/1/1/1118354ef16fe8300bf436d7fd8353a9e2aed059.png)

that feed into differential formulas like this;

![image](https://global.discourse-cdn.com/julialang/original/3X/d/f/df4a448f4cf45473d3349e4ede982428c35a7f97.png)

First you need to make a function that calculates **Fij**. I had difficulty with making a sum function ( ![image](https://global.discourse-cdn.com/julialang/original/3X/d/a/dafc07b8401c896b91285f1705677b85adc61f6d.png) ) because I thought it would be

```julia
sumF = 0.0
sumF = sumF + b[(pred), k]*bioS[k]^(1+q)

```

where `pred` is either `i` or `j` depending on the formula (because sometimes the predator becomes the prey) and `bioS` is the biomass of a species,

The feeding function then looks like

```julia
function FF(pred, prey, bioS, bm, b, h, ω, q)
    sumF = 0.0

    for k in 1:bm
        if foodWeb[(pred), k] != 0.0
            sumF = sumF + b[(pred), k]*bioS[k]^(1+q)
        end
    end
    return (((ω[pred]*b[pred,prey]*bioS[prey]*bioS[prey]^(1+q))/(1+c*bioS[pred]+ω[pred]*h[pred]*sumF))/bm[pred])
end

```

This then gets put into another function (_exactly_ how EIOceanografo showed [here](https://discourse.julialang.org/t/how-to-use-differentialequations-with-complicated-equations/20740/9), so apologies for being slow on the uptake)

The dAdT equation (second formula) then looks something like:

```julia
function dAdT(dA, foodWeb, e, bioS, X, num_plant, num_animal, num_nutrient, bm, b, h, ω, q)
    dA = zeros(num_animal)
    for i in 1:num_animal
        into = 0.0
        out = 0.0
        for j in num_nutrient+1:arr_size
            if foodWeb[i,j-num_nutrient] !=0.0
                into = into+e[j-num_nutrient]*FF(pred, prey, bioS, bm, b, h, ω, q) #here you put all the parameters needed to run FF
            end
            if foodWeb[j-num_nutrient, i+num_plant] !=0.0
                out = out+FF(pred, prey, bioS, bm, b, h, ω, q)*bioS[j] #here you put all the parameters needed to run FF
            end
        end
        into = into*bioS[num_nutrient+num_plant+i]
        met = X[num_plant+i]*bioS[num_nutrient+num_plant+i]
        dA[i] = into-out-met
    end
    return dA
end

```

And then you need to create the function to run the differential equation, which looks something like

```julia
function dBdT(dB, bioS, p, t) #the change in biomass (bioS) is what we are interested in

    dA = dB

    dB = dAdT(dA, foodWeb, e, bioS, X, num_plant, num_animal, num_nutrient, bm, b, h, ω, q) #put all the parameters needed to produce dAdT here
    return dB
end

tspan = (0.0,500.0)
p = (foodWeb, num_plant, num_animal, num_nutrient, bm, b, h, ω, q) #do not put the variable that is being changed (e.g. bioS) in here, only parameters

prob = ODEProblem(dBdT,bioS,tspan,p)
sol=solve(prob)

```

The code snippets I post won’t work by themselves, and they are also edited because I have two more differential functions that go into dBdT, but I hope this is enough for whoever is also trying to figure out nested functions.

My next task is how to put a callback into this 🤣 But I will make a new topic for it.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 29, 2019, 1:13pm UTC](https://discourse.julialang.org/t/how-to-use-differentialequations-with-complicated-equations/20740/20 "2019-05-29T13:13:53Z")

</div>

> [@HelgavonLichtenstein](#):
>
> function dBdT(dB, bioS, p, t) #the change in biomass (bioS) is what we are interested in dA = dB dB = dAdT(dA, foodWeb, e, bioS, X, num\_plant, num\_animal, num\_nutrient, bm, b, h, ω, q) #put all the parameters needed to produce dAdT here return dB end

should just be

```julia
function dAdT!(dA, bioS, p, t)
    foodWeb, num_plant, num_animal, num_nutrient, bm, b, h, ω, q = p
    for i in 1:num_animal
        into = 0.0
        out = 0.0
        for j in num_nutrient+1:arr_size
            if foodWeb[i,j-num_nutrient] !=0.0
                into = into+e[j-num_nutrient]*FF(pred, prey, bioS, bm, b, h, ω, q) #here you put all the parameters needed to run FF
            end
            if foodWeb[j-num_nutrient, i+num_plant] !=0.0
                out = out+FF(pred, prey, bioS, bm, b, h, ω, q)*bioS[j] #here you put all the parameters needed to run FF
            end
        end
        into = into*bioS[num_nutrient+num_plant+i]
        met = X[num_plant+i]*bioS[num_nutrient+num_plant+i]
        dA[i] = into-out-met
    end
    return nothing
end

```

as your ODE.

---

<div class="post-metadata">

**Author:** ![nilshg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nilshg/32/2283_2.png) [@nilshg](https://discourse.julialang.org/u/nilshg)\
**Post date:** [May 29, 2019, 2:36pm UTC](https://discourse.julialang.org/t/how-to-use-differentialequations-with-complicated-equations/20740/21 "2019-05-29T14:36:48Z")

</div>

Hmm, I wonder about this:

> `zeros(size(qwer))` is used when it is an **array**  
> `zeros(asdf)` is used when it is a _tuple_

when _what_ is an array or tuple? If you’re trying to generate an array of zeros the same length of an existing array or tuple, then `zeros(length(x))` is how you do it in both cases:

```julia
julia> a = [1,2,3]
3-element Array{Int64,1}:
 1
 2
 3

julia> t = (1,2,3)
(1, 2, 3)

julia> zeros(length(a))
3-element Array{Float64,1}:
 0.0
 0.0
 0.0

julia> zeros(length(t)) == zeros(length(a))
true

```

Just calling `zeros()` on a tuple though gives something different:

```julia
julia> zeros(t)
1×2×3 Array{Float64,3}:
[:, :, 1] =
 0.0 0.0

[:, :, 2] =
 0.0 0.0

[:, :, 3] =
 0.0 0.0

```

Similarly, I’m not sure about

> `for i in length(qwer)` is faster when it is an **array**  
> `for i in 1:asdf` is faster when it is a _tuple_

which seem to be different things - the first one is not a range, but only a single integer:

```julia
julia> for i in length(a)
           println(i)
       end
3

```

while the second one doesn’t work if `asdf` is a tuple (although it isn’t clear from your post what `asdf` actually is)