Mathematical Modeling of Social Phenomena

Dr. Courtney Brown

Assignment #6

Consider the following models of population growth, all of which are based on the logistic form.

1. dN/dt = -aN + bN²

2. dN/dt = -aN - bN²

3. dN/dt = aN + bN²

4. dN/dt = aN - bN²

with parameters a and b both given positive values. Use these models to describe the spread of rumors in various settings. That is, describe the possible birth and death of rumors as modeled by each of the four statements above. Hint: You can re-write each of the equations as dN/dt = NF(N), where the function F(N) = a - bN, for example. Then you can write about the interpretation of the parameters a and b as they are found in F(N). Note also that if you write the model as dN/dt = rN(k - N), then a=rk and b=r, where the parameter r is the initial exponetial growth rate parameter that controls the growth of the logistic model from its initial conditions, and k is the carrying capacity of the population N. Thus, feel free to mess around with the parameters to give new interpretations to the model. Choose interesting values of a and b to produce plots for each of the four stories you will have to tell.

Calculate the equilibrium values for each of the equations, and draw at least one phase-plane diagram and one over-time plot for each of the equations without solving explicitly for N. That is, use a first-order Runge-Kutta (Euler) approximation to do your definite integration. Hand in all of this together with some interpretive text.

The following SAS code may help you. Cut and paste it into your SAS editor and run it. Then play with the values of the parameters and initial values for N.

GOPTIONS lfactor=10 hsize=6 in vsize=6 in horigin=1 in vorigin=3 in;
options nocenter;
**********************************************************;
* CLASS, NOTE THAT IF YOU BEGIN A LINE WITH AN ASTERISK *
* THEN YOU CAN PUT NOTES IN YOUR PROGRAM FILES. THIS IS
* LIKE A COMMENT CARD IN SPSS. HOWEVER, REMEMBER
* TO EVENTUALLY PUT A FINAL SEMICOLON AT THE END OF YOUR COMMENTS.;
***********************************************************;
* NOTE THAT I INDENT SOME STATEMENTS. THIS
* IS JUST FOR NEATNESS.;
***********************************************************;
* COPYRIGHT (c) Courtney Brown 2004, All Rights Reserved;
* Permission granted to use this file and computer code for any nonprofit and
* educational purposes, including classroom instruction.
* No further permission required.
* Please cite source as "From www.courtneybrown.com";
***********************************************************;
DATA;
A=0.3;B=0.6;H=0.5;N=0.01;time=0;
DO T=1 TO 100;
DERIV=N*(A-(B*N));
NNEXT=N + (H*(DERIV));
time=time+H;
OUTPUT;
N=NNEXT;
END;
symbol1 color=black v=NONE f=centb i=join;
symbol2 color=black f=centb v='.' height=2 interpol=R;
symbol3 color=black f=centb v='.' height=2;
PROC GPLOT;
axis1 color=black
value=(h=1.5 f=swissb c=black)
label=(h=1.3 a=90 r=0 f=swissb c=black 'N');
axis2 color=black
value=(h=1.5 f=swissb c=black)
label=(h=1.3 f=swissb c=black 'Time');
PLOT NNEXT*time=1/
vaxis=axis1 haxis=axis2 vminor=0 hminor=0;
TITLE 'Figure 1: Time Plot with Initial N Value Below Equilibrium';

PROC GPLOT;
axis1 color=black
value=(h=1.5 f=swissb c=black)
label=(h=1.3 a=90 r=0 f=swissb c=black 'dN/dt');
axis2 color=black
value=(h=1.5 f=swissb c=black)
label=(h=1.3 f=swissb c=black 'N');
PLOT DERIV*N=1/
vaxis=axis1 haxis=axis2 vminor=0 hminor=0;
TITLE 'Figure 2: Phase Plot with Initial N Value Below Equilibrium';

DATA;
A=0.3;B=0.6;H=0.5;N=0.95;time=0;
DO T=1 TO 100;
DERIV=N*(A-(B*N));
NNEXT=N + (H*(DERIV));
time=time+H;
OUTPUT;
N=NNEXT;
END;
PROC GPLOT;
axis1 color=black
value=(h=1.5 f=swissb c=black)
label=(h=1.3 a=90 r=0 f=swissb c=black 'N');
axis2 color=black
value=(h=1.5 f=swissb c=black)
label=(h=1.3 f=swissb c=black 'Time');
PLOT NNEXT*time=1/
vaxis=axis1 haxis=axis2 vminor=0 hminor=0;
TITLE 'Figure 3: Time Plot with Initial N Value Above Equilibrium';

PROC GPLOT;
axis1 color=black
value=(h=1.5 f=swissb c=black)
label=(h=1.3 a=90 r=0 f=swissb c=black 'dN/dt');
axis2 color=black
value=(h=1.5 f=swissb c=black)
label=(h=1.3 f=swissb c=black 'N');
PLOT DERIV*N=1/
vaxis=axis1 haxis=axis2 vminor=0 hminor=0;
TITLE 'Figure 4: Phase Plot with Initial N Value Above Equilibrium';
run;
quit;