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?
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
) because I thought it would be
But I will make a new topic for it.