BYOM function derivatives.m (the model in ODEs)
Syntax: dX = derivatives(t,X,par,c)
- t is the time point, provided by the ODE solver
- X is a vector with the previous value of the states
- par is the parameter structure
- c is the external concentration (or scenario number)
Time t and scenario name c are handed over as single numbers by call_deri.m (you do not have to use them in this function). Output dX (as vector) provides the differentials for each state at t.
function dX = derivatives(t,X,par,c)
global glo % allow for global parameters in structure glo (handy for switches)
The state variables enter this function in the vector X. Here, we give them a more handy name.
Cw = X(1); % state 1 is the external concentration at previous time point Ci = X(2); % state 2 is the internal concentration at previous time point
The parameters enter this function in the structure par. The names in the structure are the same as those defined in the byom script file. The 1 between parentheses is needed as each parameter has 5 associated values.
kd = par.kd(1); % degradation rate constant, d-1 ke = par.ke(1); % elimination rate constant, d-1 Piw = par.Piw(1); % bioconcentration factor, L/kg Ct = par.Ct(1); % threshold external concentration that stops degradation, mg/L % % Optionally, include the threshold concentration as a global (see script) % Ct = glo.Ct; % threshold external concentration that stops degradation, mg/L
Calculate the derivatives
This is the actual model, specified as a system of two ODEs:
if Cw > Ct % if we are above the critical internal concentration ... dCw = -kd * Cw; % let the external concentration degrade else % otherwise ... dCw = 0; % make the change in external concentration zero end dCi = ke * (Piw*Cw-Ci); % first order bioconcentration dX = [dCw;dCi]; % collect both derivatives in one vector dX