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.

Contents

Initialization

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.

Notes

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 ODE

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

Equations

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;

Return

Pass the derivatives back

dvdt = [dNdt; dPdt; dZdt; dDdt];

end

Elapsed time is 8.619955 seconds.

Published with MATLAB® R2017a