Population growth and the earth's human carrying capacity
by Prof. Matt Miller, Department of Mathematics, University of South Carolina email: miller@math.sc.edu
A Maple worksheet to examine the results presented in the paper by Cohen, J. E. ,1995,
Science, 269 : 341-348.
Pay close attention to the sections enclosed between stars ***like this*** ; these are problems that you will be expected to discuss!
Cohen says on p. 341 that the annual rate of increase of the earth's population was 0.04% per year from AD 1 to 1650, rose to a peak of 2.1% per year around 1965-1970, and is now 1.6%. (You might want to compare these assertions with your own analyses of the data from problem 5 on ws1!) What are the consequences of these differences in rates? Using the data from Fig 1. we can predict the earth's population if the growth rate was 0.04%, 1.6%, or 2.1% over the period from AD 1 to the present. We use an exponential growth model with r the Malthusian parameter and P the population size.
> restart: with (DEtools): with(plots):
Warning, the name changecoords has been redefined
> Popgrowth:=diff(P(t), t)=r * P(t); # This is our equation in symbolic form (r is a parameter).
We will run the model to the year 2000. The initial population in the (fake) year AD 0 was 0.252 billion.
> vars:=[ P(t) ]; domain:= 0 .. 2000; init:=[ P(0)= 0.252 ];
The substitute command is very handy for making changes in the values of a parameter. Here t is the independent variable, P is the dependent variable, and r is treated as constant for any given run of the model. But we are free to adjust r and to examine the effect of such a change on the model.
> diffeq:= subs(r = 0.0004, Popgrowth); # This is our equation with r = 0.0004.
> DEplot( diffeq, vars, domain, [init], stepsize=20, title=`World Population in Billions: r=0.0004` , arrows = LINE );
In this case we can also get a formula solution.
> formula:= dsolve( {diffeq, P(0) = 0.252}, P(t) );
> evalf( subs( t = 1650, formula) );
> evalf(subs( t = 2000, formula) );
Was a growth rate of 0.0004 (0.04%) fast enough to reach the current (year 2000 for convenience) population size of 2.6 billion? What if we used the current growth rate of 1.6% instead? Go back to where we substituted r = .0004 into the DE, and try r = .016 instead.
> # Use the mouse to skip to the next visible > if you like. To open the optional section, just hit Return now, or click on the little box; to close it, click on the little box.
Optional section for those comfortable with piecewise defined functions
> restart: with(plots): with(DEtools):
SInce none of the listed rates held for all of recent (AD) recorded history, let's try to see if we can use the figures Cohen gave us in a more realistic way. Recall that r = .04% up to 1650, then increases to r = 2.1% around 1965, and decreases to r = 1.6% in the present. Let's plot r as a function of t. Can you see how we can get a formula for r at the times in between as well? (HInt: what if you use style = LINE?)
> plot( [ [0, .0004], [1650, 0.0004], [1965, .021], [1995, .016] ], style = POINT, title = `r vs. t`);
Ok, did you figure it out? If not, take a look at the next line, and see if you can make sense of it.
> r:= piecewise( t < 1650, 0.0004, t < 1965, (.021-.0004)/(1965-1650) * (t - 1650) + .0004,(.016-.021)/(1995-1965) * (t - 1965) + .021);
> formula:= dsolve( {diff(P(t), t) = r * P(t), P(0) = 0.252}, P(t));
For reasons that I do not understand, sometimes the above command works, and sometimes it doesn't. If it works, it gives some spurious "undefined" answers, but for whatever it's worth, let's see what it gives for the population in 1995 and 2000.
> evalf( subs( t = 1995, formula)); evalf( subs( t = 2000, formula));
If that doesn't work we can do it piece by piece, using the final value for one segment to give the beginning value of the next one. First we do the computation symbolically.
> formula1:= dsolve( {diff(P(t), t) = 0.0004 * P(t), P(0) = p[0]}, P(t) );
> formula2:= dsolve({diff(P(t), t) = ( m[1]*t + b[1]) * P(t), P(1650)= p[1650]}, P(t));
> formula3:= dsolve({diff(P(t), t) = ( m[2] * t + b[2]) * P(t), P(1965)= p[1965]}, P(t));
Now we start filling in the parameter values, and the population values at the end of each segment, beginning of the next. Thanks to Doug Meade for helping out with the syntax for this section!
> p[0]:= 0.252;
> p[1650]:= evalf(subs( t = 1650, rhs(formula1)));
> m[1]:= (.021-.0004)/(1965-1650); b[1]:= (.021-.0004)/(1965-1650)*(- 1650) + .0004;
> m[2]:= (.016-.021)/(1995-1965); b[2]:= (.016-.021)/(1995-1965)*(-1965) + 0.021;
> p[1965]:= evalf(subs( t = 1965, rhs(formula2)));
> evalf(subs( t = 1995, formula3));
> evalf(subs( t = 2000, formula3));
Of course we can also get a graph by using DEplot (or by plotting our formula solution).
> Popgrowth:= diff(P(t), t) = r * P(t);
> DEplot(Popgrowth, [P(t)], 0 .. 2000, {[0, 0.252]}, stepsize=20, title=`World population in billions: r defined piecewise`);
> formula:= piecewise(t<1650, rhs(formula1), t<1965, rhs(formula2), rhs(formula3));
> plot( formula, t = 0 .. 2000, title=`World population in billions: r defined piecewise`);
How do these results compare to Cohen's Fig. 1, or with his reference in the text to the current, 1995, world population as 5.7 billion? The axis scale is in billions, as in Cohen's graph. Note, however, that his graph is on a LOG scale and your graphs are on a linear scale. Lets take advantage of our formula answer here.
> logplot( formula, t = 0 .. 2000, axes=FRAME, title=`World population in billions: r defined piecewise`);
Is it possible that separate growth rates, over successive time periods, can give population predictions consistent with the data of Cohen's Fig. 1? If we have not succeeded in doing this what do you think has gone wrong? In any event is this a plausible model? a useful one for predicting future growth?
>
Cohen mentions in the Fig. 1 legend that the population growth was "faster than exponential" between 1400 and 1970.
1. What does he mean by this?
2. What evidence does he have?
3. How could you rewrite the equation P' = rP to describe a population that has faster than exponential growth?
Cohen proposes "idealized mathematical models" for the race betwen human population growth and human carrying capacity. Suppose that it is possible to define human carrying capacity K(t) as a numerical quantity measured in numbers of individuals. Suppose also that P(t) is the total number of individuals in the population at time t. The logistic model then gives us the rate equation P' = rP(K - P). (NOTE: this is not the standard version of the logistic equation, which would have the right hand side divided by K; qualitatively, however, the results are the same.) In the simplest case K(t) is taken to be constant, reflecting the stability of the natural environment over a time span that includes many generations: let's say, to begin with, K1= 0.252789 billion. Again the constant r is the Malthusian parameter, in this simulation roughly 1.4%. Cohen's simulations run for 2500 years with a step size of 20 years, and again an initial population of 0.252 billion.
> Bounded_growth:= diff(P(t),t) = r * P *(K1 - P); vars:= [ P(t) ]; domain:= 0.. 2500; init:= [ P(0)= 0.252];
> DEplot(subs( {r = 0.0014829, K1 = 0.252789}, Bounded_growth), vars, domain, [init], arrows = LINE, stepsize=20);
How does this pattern of population growth compare to what you modelled above? What effect does the substitution of other constant values of r and K1 have on the solutions?
Cohen's intuition is that, at least for human populations, the population itself can affect the environment in such a way that the carrying capacity will change. Thus he goes on to investigate a more complicated model in which the carrying capacity K(t) is allowed to vary, since "every human being represents hands to work, and not just another mouth to feed" (George Bush). He assumes the rate of change of carrying capacity is proportional to the rate of change of population size, with the Condorcet parameter c as the constant of proportionality.
We define Eq4 and Eq5 to accomodate Cohen's time-varying carrying capcity K(t) rather than the fixed carrying capacity K1 that we used before. Do you see where Eq5 comes from?
> Eq4:= diff(P(t), t) = r*P(t)*(K(t)- P(t)); Eq5:= diff(K(t), t) = c* ( r*P(t)*(K(t) - P(t)) );
> vars:= [ P(t), K(t) ]; domain:= 0 .. 2500 ; init2:=[P(0)=0.252, K(0)=0.252789];
First try modelling the populaton with the Condorcet parameter c < 1, say, c = 0.5. This means that each additional person increases the carrying capacity by less than s/he consumes.
> rate_eqns:= subs( { c = 0.5, r = 0.0014829 }, [ Eq4, Eq5 ] );
> Kplot:= DEplot(rate_eqns, vars, domain, [init2], stepsize=20, scene=[t, K], arrows=NONE):
> Pplot:= DEplot(rate_eqns, vars, domain, [init2], stepsize=20, scene=[t, P], arrows=NONE):
> display( {Kplot, Pplot }, title = `Condorcet growth model (c=0.5)`);
***Can you identify which curve is which? (If not, try adding the option linecolor=BLUE to the DEplot command that produces Kplot.) How do they compare to Cohen's Fig. 4? We can now look at the case where each person increases resources by MORE than s/he consumes (c > 1): try using c = 2.5 in the lines above. Is the result consistent with Cohen's discussion at the end of page 343? Try running the model with c = 1.0 (each person adds to the carrying capacity just as much as s/he consumes). Does the population grow as Cohen describes in each case?***
For the mathematically intrepid, footnote 35 gives Cohen's reasoning. Can you put this in more intuitive terms: for instance, what does he mean by "virtual carrying capacity" and "virtual Malthusian parameter"?
***It is also interesting to consider the case of a negative Condorcet parameter. What would happen to the world's present population in 50 years, assuming a present carrying capacity of 18 billion, a 1.6% intrinsic growth rate, and a Condorcet value c = -0.5? (Go back and change all the relevant parameter values, including initial conditions and the domain.)***
We define Eq6 to accomodate Cohen's model of dilution, but not depletion, of resources; L is the Mill parameter.
> Eq4:= Eq4; Eq6:= diff(K(t),t) = (L/P(t))*(r*P(t)*(K(t) - P(t)));
> rate_equations:= subs( { L = 3.7 , r = 0.0014829 } , [ Eq4, Eq6 ] );
> Kplot:= DEplot(rate_equations, vars, domain, [init2], stepsize=20, scene=[t, K], arrows=NONE, linecolor=BLUE): Pplot:= DEplot(rate_equations, vars, domain, [init2], stepsize=20, scene=[t, P], arrows=NONE ):
> display({Kplot, Pplot}, title = `Cohen Fig. 4`);
***Cohen claims that the long term behavior of K(t) and P(t) is independent of the value of r. Is this true, and if so, how do you account for it? What effect does the Malthus parameter have on the solutions? (Recall you looked at similar questions in the context of the logistic model.) Read footnote 39. Why did Cohen chose 3.7 for the parameter L? What is he assuming about the shape of the relation between K(t) and time that led him to use this value? Do you believe he made a reasonable choice? Do you think Cohen's curve-fitting approach is likely to give a better estimate of the upper bound on human populations than the estimates in Fig 3? Why or why not?***
