Solving an ODE (NPZD)
This script shows you how to solve a system of ODEs (the NPZD model in this case). It uses the built in Matlab function ode45 which uses an explicit Runga-Kutta (4,5) method. It has an adaptive time step so time is not equally spaced. Please note that this script defines functions at the end, which is only supported by MATLAB 2016b or later.
We begin by defining all the parameters in the NPZD model. Since we want to use them inside of the function that defines the ODE we define them to be "global variables" (this means that they can be used inside of a function which normally has a separate scope).
global kn mu mp kp k Im alpha mz beta mu = 2; kn = 1; mp = 0.1; mz = 0.2; alpha = 0.7; beta = 0.3; kp = 0.7; Im = 1.5; k = 0.1;
Warning: The value of local variables may have been changed to match the globals. Future versions of MATLAB will require that you declare a variable to be global before you use that variable. Warning: The value of local variables may have been changed to match the globals. Future versions of MATLAB will require that you declare a variable to be global before you use that variable.
We define the ICs in the order: N; P; Z; D.
V0 = [1.4; 0.2; 0.2; 0.2];
Running the solver
tic % To use the ode45 function you need to provide: the function that defines % the DE, a vector indicating the initial time and final time, and the ICs % of your variable. Please see the npzd function at the end for more about % the DE. [t, V] = ode45(@func_npzd,[0 60],V0); toc plot(t,V)
You can try changing the various parameters and see how the system changes. Some interesting parameters to change are the ICs, k and mu.
ode45 is the standard ODE solver and will work for most cases. If you have a stiff ODE system, this method will not work and you will have to use another method such as ode23s. To see more about the different built in Matlab ODE solvers check out: https://www.mathworks.com/help/matlab/ordinary-differential-equations.html
NPZD is the differential equation for modeling Nutrient, Plankton, Zooplankton and Detritus.
function dvdt = func_npzd(t,V)
All variables are contained in the vector V and their derivatives are returned in dvdt. Requires that the parameters have been defined globally.
Explicitly declare names
We define each variable explicitly to make it easier to write the equations
N=V(1); P=V(2); Z=V(3); D=V(4);
Call the globally defined variables to have them defined in our scope.
global kn mu mp kp k Im mz alpha beta
The set of DEs in the form dY/dt = f(Y);
dPdt = (N/(kn + N))*mu*P - mp*P-(P/(kp+P))*Im*Z; dZdt = alpha*(P/(kp+P))*Im*Z-mz*Z; dDdt = mp*P + (1-alpha-beta)*(P/(kp+P))*Im*Z-k*D; dNdt = -dPdt-dZdt-dDdt;
Pass the derivatives back
dvdt = [dNdt; dPdt; dZdt; dDdt];
Elapsed time is 8.619955 seconds.