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 = 10.0*(u-u) du = u*(28.0-u) - u du = u*u - (8/3)*u 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