World  Population Growth

Part 1. Background

Enter the following line to load the needed Mathematica packages.

<< Graphics`Graphics`
<< Graphics`PlotField`

Answer the questions in step 1 here.

Here are the world population data as best as we know them.  The data columns are given separately as Years and Pops (populations in millions), and the columns are then threaded together into WorldData.  Enter the following commands:

Years = {1000, 1650, 1750, 1800, 1850, 1900, 1910, 1920, 1930, 1940, 
          1950, 1955, 1960, 1965, 1970, 1975, 1980, 
          1985};
Pops = {200, 545, 728, 906, 1171, 1608, 1750, 1834, 2070, 2295, 2517, 
          2780, 3005, 3345, 3707, 4086, 4454, 
          4850};
WorldData  =  Thread[{Years,  Pops}]

Enter the next command to plot the data

PopPlot  =  
        ListPlot[WorldData, 
          PlotStyle -> {RGBColor[0, 0, 1], 
              PointSize[0.015]}];

Answer to step 2:

By changing ListPlot to LogListPlot in the preceding command, you can construct a semilog plot of the data.  Do that, and then answer step 3.

Answer to step 3:

LogListPlot[WorldData, 
        PlotStyle -> {RGBColor[1, 0, 0], 
            PointSize[0.015]}];

Part 2. The Coalition Model

The von Foerster paper argues that the differential equation modeling growth of world population P as a function of time t might have the form

	(dP/dt) = k (P^( r + 1)),
	
where r and k are positive constants.

Explain why

	((P(i + 1) - P(i - 1))/(t(i + 1) - t(i - 1)))
	
is a good estimate of (dP/dt) at time (Subscript[t, i]).

Answer:

The following commands will construct the symmetric difference quotients as a list SDQ.

n  =  Length[Pops]
SDQ = 
      Table[(Pops[[i + 1]] - 
              Pops[[i - 1]])/(Years[[i + 1]] - 
              Years[[i - 1]]),  {i,  2,  n - 1}]

The following commands construct a log-log plot of SDQ versus population.  Note that there is no SDQ for  i = 1  or  i = n,  so we use a shortened list of population data.

Pops2  =  
        Delete[Pops,  {{1},  {n}}];
Pops2SDQData  =  
        Thread[{Pops2,  SDQ}];
LogLogListPlot[Pops2SDQData, 
        PlotStyle -> {RGBColor[0.2, 0, 0.8], 
            PointSize[0.015]}];

Is it possible that (dP/dt) is a power function of P?

Answer:

Use the built-in least squares procedure to find the line that best fits the population-slope data.  Note that we have to take log10 of the data entries to match the numbers in the plot.

LogLogData  =  
        Thread[{Log[10, Pops2],  Log[10, SDQ]}];
Clear[x];
Fit[LogLogData,  {1,  x},  x]

Define the constants k and r.

r  =  ?;
k  =  ?;

The following commands draw a slope field for our model differential equation and a sample solution.  (Click the lower-right corner of the image and drag to enlarge the plot.)

		(dP/dt) = k (P^( r + 1))

slope  =  
        PlotVectorField[{1,  k * P ^(1 + r)},  {t,  1000,  
            2000},  {P,  0,  8000},  ScaleFunction -> (1 &),  
          HeadLength -> 0,  HeadWidth -> 0,  
          AspectRatio ->  1,  Axes -> True,  
          DisplayFunction -> Identity];
diffEq  =  P'[t] ==  k*P[t]^(1 + r);
solution  =  
        NDSolve[{diffEq,  P[1940] ==  Pops[[10]]},  
          P,  {t,  1000,  2000}] ;
integralcurve  =  
        Plot[Evaluate[P[t] /. solution],  {t,  1000,  2000},  
          DisplayFunction -> Identity];
Show[slope,  integralcurve,  
        DisplayFunction -> $DisplayFunction];

The next command adds the original data plot to the slope field plot.

Show[PopPlot, slope,  integralcurve];

What do you think about the Coalition Model as a description of the historical data?

Answer:

Part 3. Solving the Coalition Model

We may use separation of variables to write the differential equation

		(dP/dt) = k (P^( r + 1))
in the form
		(dP/(P^(r + 1))) = k dt

and then integrate both sides to get

		\(TraditionalForm\`∫\(1\/P\^\((1 + r)\)\) dP\) = ∫ k dt
	
Carry out the integrations, and solve for P as a function of t.  (We first purge the values of k and r so you can do the calculation symbolically.)

Clear[k,  r];

Solution:

Show that your solution can be written in the form  

			P = 
        (1/((( r k (T - t) ))^( 1/ r)))      

for some constant T.  Specifically, how is T related to your constant of integration C?

Answer:

Part 4. When is Doomsday?

The following commands restore the values of k and r you computed in Part 2.  (If your values were different, change these commands accordingly.)

r  = 1.190456141;
k  = 10^(-5.957305457);

The next group of commands constructs a log-log plot of P versus (T - t) for your choice of Doomsday time T.  To help you find the best choice, a line is drawn between the first and last points.

T  =  ?;

TimeLeft  = Table[T - Years[[i]], {i, n}];
DoomsData  =  Thread[{TimeLeft, Pops}];
llplot = 
        LogLogListPlot[DoomsData,  
          AxesLabel -> {"T-t",  "P"},  
          Prolog -> {PointSize[0.015], RGBColor[0, 0, 1]}, 
          DisplayFunction -> Identity];
line  =  
        ListPlot[{{Log[10,  TimeLeft[[1]]],  
              Log[10,  Pops[[1]]]},  {Log[10,  TimeLeft[[n]]], 
              Log[10,  Pops[[n]]]}},  PlotJoined -> True,  
          PlotStyle ->  RGBColor[1, 0, 0], 
          DisplayFunction -> Identity];
Show[llplot, line,  
        DisplayFunction -> $DisplayFunction];

Enter the model function, and plot it together with the data.

P[t_] :=  
        1/(r*k*(T - t))^(1/r);
ModelPlot  =  
        Plot[P[t], {t, 1000, 2050},  PlotRange -> {0,  8000},  
          DisplayFunction -> Identity];
Show[ModelPlot, PopPlot, 
        DisplayFunction -> $DisplayFunction];

Experiment with adjustments in k and/or r to see if you can get a better fit.

k  = ?;
r  = ?;
ModelPlot  =  
        Plot[P[t], {t, 1000, 2050},  PlotRange -> {0,  8000},  
          DisplayFunction -> Identity];
Show[ModelPlot, PopPlot,  
        DisplayFunction -> $DisplayFunction];

Calculate your "predictions" for 1990 and 1995.

Answer:

Calculate predictions for 2000, 2010, and 2020.

Answer:

How do your predictions compare with Census Bureau estimates and projections?

Answer:

Part 5. Summary