How to use DifferentialEquations with complicated equations


#1

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

#2

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

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

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.


#5

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

#6

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:


#7

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

#8

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.


#9

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)

#10

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.