# 8.2: Code for two species

$$\newcommand{\vecs}{\overset { \rightharpoonup} {\mathbf{#1}} }$$ $$\newcommand{\vecd}{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash {#1}}}$$$$\newcommand{\id}{\mathrm{id}}$$ $$\newcommand{\Span}{\mathrm{span}}$$ $$\newcommand{\kernel}{\mathrm{null}\,}$$ $$\newcommand{\range}{\mathrm{range}\,}$$ $$\newcommand{\RealPart}{\mathrm{Re}}$$ $$\newcommand{\ImaginaryPart}{\mathrm{Im}}$$ $$\newcommand{\Argument}{\mathrm{Arg}}$$ $$\newcommand{\norm}{\| #1 \|}$$ $$\newcommand{\inner}{\langle #1, #2 \rangle}$$ $$\newcommand{\Span}{\mathrm{span}}$$ $$\newcommand{\id}{\mathrm{id}}$$ $$\newcommand{\Span}{\mathrm{span}}$$ $$\newcommand{\kernel}{\mathrm{null}\,}$$ $$\newcommand{\range}{\mathrm{range}\,}$$ $$\newcommand{\RealPart}{\mathrm{Re}}$$ $$\newcommand{\ImaginaryPart}{\mathrm{Im}}$$ $$\newcommand{\Argument}{\mathrm{Arg}}$$ $$\newcommand{\norm}{\| #1 \|}$$ $$\newcommand{\inner}{\langle #1, #2 \rangle}$$ $$\newcommand{\Span}{\mathrm{span}}$$

Below is computer code for two-species dynamics—a logical expansion of the code for one-species dynamics you saw earlier. The code here produced the graph in Figure 8.1.2 and, with other values for Ni, ri, si,j, also produced the graphs in Figures 8.1.4 and 8.1.6.

N1=.01; N2=.01;

r1=0.5; r2=0.8;

s11=-0.08; s12=-0.03; s21=-0.09; s22=-0.06;

t=0; dt=.0001; tmax=75; step=0;

print(c(t, N1, N2));

while(t<tmax)

{ dN1=(r1+s11*N1+s12*N2)*N1*dt;

dN2=(r2+s21*N1+s22*N2)*N2*dt;

N1=N1+dN1; if(N1<0) N1=0;

N2=N2+dN2; if(N2<0) N2=0;

t=t+dt; step=step+1;

if(step==1000)

{ print(c(t, N1, N2)); step=0; }

This code uses a time step dt of 1/10,000 but displays only every 1000th time step, using the variable named step. This is a reliable way to display output periodically, because step takes on only integer values 0, 1, 2.... Watching variable t for the same purpose may not be reliable; t has rounding errors in its fractional portion, which typically is represented only approximately by the computer.

By the way, here is a crucial requirement for coding multi-species dynamics: you must update all of the Ni together. You might be tempted to shorten the code by writing the following, getting rid of the dN1 and dN2 variables.

N1=N1+(r1+s11*N1+s12*N2)*N1*dt; if(N1<0) N1=0;

N2=N2+(r2+s21*N1+s22*N2)*N2*dt; if(N2<0) N2=0;

This shortened code, however, contains a serious bug that could be difficult to detect. (A “bug” is the computerese term for any mistake in computer code, typically one that goes undetected.) The calculation of N2 on the second line uses the new value of N1, not the present value. This will generate subtly incorrect results that could be costly—if, for example, you were using the code to project the course of an epidemic.

Careful code reviews with experienced colleagues are one way to help avoid such bugs. If the problem is important enough, having the code written independently by two or more people not communicating with each other, except by accepting the same specifications and comparing the final results, is a way to approach correctness. This is like independent replication of scientific experiments. Other methods of ensuring that code is correct are addressed later.

This page titled 8.2: Code for two species is shared under a CC BY-NC 4.0 license and was authored, remixed, and/or curated by Clarence Lehman, Shelby Loberg, & Adam Clark (University of Minnesota Libraries Publishing) via source content that was edited to the style and standards of the LibreTexts platform; a detailed edit history is available upon request.