How to use DifferentialEquations with complicated equations

Hello,

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

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

Does anyone know how this could be done?

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
1 Like

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 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., https://github.com/JuliaNLSolvers/NLsolve.jl.

3 Likes

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 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 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) and is also supported by DifferentialEquations.jl.

4 Likes

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.

4 Likes

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

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

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

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

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

2 Likes

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.

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

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.

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 from Wikipedia. Notice how I have a helper function interaction which is defined outside, but then called within, the main simulation function:

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)
4 Likes

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.

4 Likes

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:

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?

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

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

Don’t use global variables 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.

1 Like

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

3 Likes

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.

4 Likes

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)

Consider:

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:

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> 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 an excellent resource for this tasks and I’d recommend you take a couple of days to work through it step-by-step!

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

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;

image

that feed into differential formulas like this;

image

First you need to make a function that calculates Fij. I had difficulty with making a sum function (image ) because I thought it would be

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,

but instead it is

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

(But I don’t understand why)

The feeding function then looks like

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, so apologies for being slow on the uptake)

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

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

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 :rofl: But I will make a new topic for it.

1 Like

should just be

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.

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

1 Like