Mathematical Modeling of Social Phenomena

Dr. Courtney Brown

Assignment #6

Consider the following logistic model of population growth: dN/dt = rN(k - N), with parameters r and k both given positive values. Use this model to describe the spread of rumors following at least two scenarios. That is, describe the possible birth and death of rumors as described by this model under two situations which you feel are realistic, where the population N represents people who believe a given rumor. Hint: Note that parameter r is the initial exponential growth rate parameter that controls the growth of the logistic model from its initial conditions, and k is the carrying capacity (or limit) of the population N that is vulnerable to believing the rumor. Note that this limit will usually not be equal to the total population. Think about why this would be so. You can also re-write this by using N as a proportion of a population, in which case dN/dt = rN(1 - N/k), and 1 is now the population limit (for a very gullible population).

What is the equilibrium value for each of your scenarios, and where is the inflection point where positive feedback switches to negative feedback in the rumor growth process? Draw at least one phase-plane diagram and one over-time plot for each scenario without solving explicitly for N. That is, use a Runge-Kutta to do your definite integration. Hand in all of this together with some interpretive text.

The following R code get you started. To see how it works, cut and paste it into R and run it. Then play with the values of the parameters and initial values for N. The R code below uses an Euler (or first-order Runge-Kutta) to do the integration. You get a higher grade for this assignment if you successfully re-write this code to use a fourth-order Runge-Kutta. To do this, you will want to put the model into an R function, and then call it multiple times in the Runge-Kutta. If you do this, be sure to email me the R code so that I can check that it works (and thus give you the higher grade!). You can find an example (plus an explanation) of a fourth-order Runge-Kutta written using SAS in your text on differential equations.

r <- 0.1 # The growth rate parameter
k <- 100 # The limit of N. Note that N can exceed this limit temporarily due to exogenous forces.
h <- 0.1 # The step size for the Runge-Kutta
N <- 1 # The initial population that believes the rumor
time <- 0 # The beginning of time
mymodel <- 0 # Initializing the derivative model
iterations <- 1000 # The number of times we want to iterate the Runge-Kutta
for (i in 1:iterations) {
mymodel[i] <- r*N[i]*(1 - (N[i]/k))
if (i <= iterations-1) {
N[i+1] <- N[i] + (h*mymodel[i]) # This is the Euler. A fourth-order Runge-Kutta will be more accurate.
time[i+1] <- time[i] + h }
}

mydata <- as.data.frame(cbind(N, time, mymodel))
attach(mydata)
parameters <- cbind(r, k)
parameters
mydata
plot(time, N, type="l", ylab="Population N", xlab="Time", lwd=2)
title(main = list("Figure 1: Rumor believers (N) over time", cex=1.5,
col="black", font=1))

windows()

plot(N, mymodel, type="l", ylab="Model", xlab="N", lwd=2)
title(main = list("Figure 2:Phase Plot with Initial N Value Below Equilibrium", cex=1.5,
col="black", font=1))

* Below is a SAS program that does the same thing as the R program above;

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;