Logistic Growth Model

Part 1. Background

Enter the following to load the necessary package.

<< Graphics`PlotField`

Part 2. Equilibria

Answer in the cells below the questions in Part 2.

Answer to question 1:

Answer to question 2:

Answer to question 3:

Answer to question 4:

Answer to question 5:

Answer to question 6:

Answer to question 7:

Answer to question 8:

Part 3. Inflection Points and Concavity

The following command defines a "derivative of P" function DP as the right-hand side of the logistic differential equation.  Calculate the second derivative of P by differentiating DP with respect to t.  The answer should contain a Chain Rule factor of ((dP(t))/dt) although you may have to factor the answer to see that.

DP[t_] :=  r*P[t]*(1 - P[t]/K);

What values of P make the second derivative have the value 0?

Answer:

Answers to the rest of Part 3:

Part 4. Symbolic Solutions

After separating variables, the logistic equation has the form

		(dP/(P (1 - (P/K)))) = r dt.

Integrate both sides, and solve for P using the input cell below.

The following commands define the logistic differential equation, LDE, and use the built-in differential equation solver to solve it.

LDE  =  P'[t] == DP[t];
DSolve[LDE, P[t], t]

Answer to question 4:

Answers to questions 5, 6, and 7:

The next group of commands will set values of r, K, and P[0], draw the direction field, and superimpose your solution-- which you have to fill in where you see ?.  Vary the values of  r and K to confirm visually that your solution formula is correct.

r  =  0.5;
K  =  1;
P0  =  0.1;
dP  =  DP[t] /. P[t] -> P;
dirfld  =  
        PlotVectorField[{1,  dP}, {t, 0, 10}, {P, 0, 1.5}, 
          ScaleFunction -> (1 &), HeadWidth -> 0,  
          DisplayFunction -> Identity];
solution   =  ?;
solplot  =  
        Plot[solution, {t, 0, 10},  PlotStyle -> RGBColor[0, 0, 1],  
          AxesLabel -> {"t", "P"}, 
          DisplayFunction -> Identity];
Show[solplot,  dirfld, 
        PlotLabel -> "Logistic Growth Model", 
        DisplayFunction -> $DisplayFunction];

Part 5. Fitting a Logistic Model to Data, I

Here is the population data from 1790 to 1940. We don't need the dates for this calculation, but we show them to help keep track of the year for each population.  The first command, incidentally, suppresses the warning that Mathematica would have generated about the similarity between our variable date and the built-in function Date.

Off[General::spell1];

date  =  {1790, 1800, 1810, 1820, 1830, 1840, 1850, 1860, 1870, 
          1880, 1890, 1900, 1910, 1920, 1930, 1940};
pop  =  {3.929, 5.308, 7.240, 9.638, 12.866, 17.069, 23.192, 
          31.443, 38.558, 50.156, 62.948, 75.996, 91.972, 105.711, 122.775, 
          131.669};

Next we construct a sequence that approximates the ratio of dP/dt to P at each relevant time using symmetric difference quotients.  The following calculations construct the sequence of data points to plot, where derpop represents the ratio of derivative to population, and newpop represents the shortened list of populations.

derpop  =  
        Table[(pop[[i + 2]] - pop[[i]])/(20*
                pop[[i + 1]]), {i, 1, Length[pop] - 2}];
newpop  =  Delete[pop,  {{1},  {Length[pop]}}];
popdata  =  Thread[{newpop, derpop}]
dataplot  =  
        ListPlot[popdata,  
          PlotStyle -> {RGBColor[0, 0, 1], PointSize[0.015]}, 
          AxesLabel -> {"P", "ratio"}];

Enter commands to estimate the slope and intercept of a line that will fit the plotted points.  Then use the slope and intercept to estimate r and K.  Be sure to define r and K before going on to the next step.

r  =  ?;
K   =  ?;

Check your work by plotting (r/(1 - (P/K))) against P on top of your ratio data:

fitplot  =  
        Plot[r*(1 - P/K), {P, 0, 130},  
          DisplayFunction -> Identity];
Show[dataplot, fitplot, 
        DisplayFunction -> $DisplayFunction];

Answer to question 5:

The following commands extend the logistic test to data through 1990. Enter the commands, and then answer question 6.

date  =  {1790, 1800, 1810, 1820, 1830, 1840, 1850, 1860, 1870, 
          1880, 1890, 1900, 1910, 1920, 1930, 1940, 1950, 1960, 1970, 1980, 
          1990};
pop  =  {3.929, 5.308, 7.240, 9.638, 12.866, 17.069, 23.192, 
          31.443, 38.558, 50.156, 62.948, 75.996, 91.972, 105.711, 122.775, 
          131.669, 150.697, 179.323, 203.185, 226.546, 248.710};
derpop  =  
        Table[(pop[[i + 2]] - pop[[i]])/(20*
                pop[[i + 1]]), {i, Length[pop] - 2}];
newpop  =  Delete[pop,  {{1}, {Length[pop]}}];
popdata  =  Thread[{newpop, derpop}]
dataplot  =  
        ListPlot[popdata, 
          PlotStyle -> {RGBColor[0, 0, 1], 
              PointSize[0.015]}];

Part 6. Fitting a Logistic Model to Data, II

The following commands will plot the solution,

		((K (Subscript[P, 0]))/((Subscript[P, 0]) + (K - (Subscript[P, 0])) (e^(-rt)))),

of the logistic differential equation, where r and K have the values you assigned in Part 5, and (Subscript[P, 0]) is the population at t = 0 in 1790, that is, P0 = 3.929.

P0  =  3.929;
P[t_] :=  K*P0/(P0 + (K - P0)*Exp[-r*t]);
curveplot  =  
        Plot[P[t], {t, 0, 150}, 
          AxesLabel -> {"t", "P"}];

Now we redefine the time scale to start at 1790, and we reset the population data to end at 1940. Then we add a plot of these data to your solution curve.

newtime  =  {0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 
          130, 140, 150};
pop  =  {3.929, 5.308, 7.240, 9.638, 12.866, 17.069, 23.192, 
          31.443, 38.558, 50.156, 62.948, 75.996, 91.972, 105.711, 122.775, 
          131.669};
popdata  =  Thread[{newtime, pop}]
dataplot  =  
        ListPlot[popdata, 
          PlotStyle -> {RGBColor[0, 0, 1], PointSize[0.015]},  
          DisplayFunction -> Identity];
Show[curveplot, dataplot,  
        DisplayFunction -> $DisplayFunction];

Next we construct a plot of ln(P/(K - P)) versus t.

lnquot  =  
        Table[Log[pop[[i]]/(K - pop[[i]])],  {i, 1, 
            Length[pop]}];
lnquotdata  =  Thread[{newtime, lnquot}]
lnquotplot  =  
        ListPlot[lnquotdata, 
          PlotStyle -> {RGBColor[0, 0, 1], 
              PointSize[0.015]}];

If necessary, enter a new definition of K and re-enter the commands to construct lnquotplot.  Continue until the plot is reasonably straight.  Then calculate r from your last plot.

If you have not already done so, enter your value for r here.  Then redefine the solution curve, and plot it with the data.

r  =  ?;
P[t_] :=  K*P0/(P0 + (K - P0)*Exp[-r*t]);
curveplot  =  
        Plot[P[t], {t, 0, 150}, 
          DisplayFunction -> Identity];
Show[curveplot, dataplot, 
        DisplayFunction -> $DisplayFunction];

If necessary, adjust K and r until you get a good fit.

Part 7. Summary

Answers to the questions in Part 7: