>

BDH1-1.mws

Modelling Via Differential Equations

Section 1.1 - Blanchard, Devaney, Hall

11 January 2000

updated to Release 5.1 by Doug Meade

18 January 1998

Doug Meade (modified by M^2)

>

Introduction

Getting Started

Every Maple worksheet should begin by re-initializing the Maple "kernel" and loading the additional packages that we are most likely to use.

> restart;

> with( plots ):

> with( DEtools ):

>

Unlimited Population Growth

One model for unlimited population growth is represented by the following first-order ordinary differential equation.

> MODEL := diff( P(t), t ) = k * P(t);

Question

Identify all independent variables, dependent variables, and parameters in the problem.

Equilibrium Solution

(See the text for a definition of equilibrium solution.)

> EQUILeqn := rhs( MODEL ) = 0;

> EQUILsol := solve( EQUILeqn, { P(t) } );

Qualitative Analysis

> InitCond := P(0) = p[0];

The text provides a method by which we can determine that solutions to this model are increasing if > 0 and > 0 and decreasing if < 0 and > 0. This type of analysis can be done for other equations but is often quite difficult.

>

Analytic Solution

We will learn specific techniques for solving special classes of differential equations. The Maple command for finding explicit solutions to differential equations is dsolve . The arguments are the initial value problem (a set containing the differential equation and the initial condition) and the function for which we are trying to solve.

> SOLN := dsolve( { MODEL, InitCond }, P(t) );

Warning

While Maple is good at solving many differential equations, there are times when the results produced by dsolve are quite different from what we obtain by hand. If you encounter this, you should not automatically assume Maple is correct. Instead, I recommend that you use Maple to test which (if either) of the solutions actually satisfies the differential equation (and initial condition).

>

The U.S. Population

The computations used to determine values for the parameters k and from the census data are easily reproduced in Maple. The idea is that two data points can be used to construct a system of two equations for the two unknown constants. While it is possible to exactly duplicate the steps shown in the text to determine k and , it is simpler to ask Maple to solve the system.

>

> EQ1790 := subs( P(t)= 3.9, t= 0, SOLN );

> EQ1800 := subs( P(t)= 5.3, t=10, SOLN );

> PARAMETERS := solve( { EQ1790, EQ1800 }, { k, p[0] } );

Thus, this model predicts the U.S. population is given by the following function.

> solP := eval( rhs(SOLN), PARAMETERS );

>

Here is the actual data from the text.

> data:= [1790,3.9], [1800,5.3], [1810,7.2], [1820,9.6], [1830,12], [1840,17], [1850,23], [1860,31], [1870,38], [1880, 50], [1890,62], [1900,75], [1910,91], [1920,105], [1930,122], [1940,131], [1950,151], [1960,179], [1970,203], [1980,226], [1990,249], [2000,`NA`], [2010,'NA'], [2020,'NA'], [2030,'NA'], [2040,'NA'], [2050,'NA'];

> shifted_data:= seq( [ 10*(k-1), data[k,2] ], k = 1 .. nops([data]) );

>

To evaluate the reasonableness of this prediction we will often plot the function. However, since the text provides a table with census values from 1790 through 2050, let's produce predicted populations for each decade (and a few decades into the next millenium). Notice that t in the formula is not "date" -- the year 1790 corresponds to the value t = 0. We produce a second list, called shifted_data, that uses t = 0 as the starting time.

> soln_points:= NULL:

> for t from 0 to 260 by 10 do

> printf( `%4d %10.4f %10.4a\n`, t+1790, solP, shifted_data[t/10+1,2] );

> soln_points:= soln_points, [t, solP]:

> od:

1790 3.9000 3.9

1800 5.3000 5.3

1810 7.2026 7.2

1820 9.7881 9.6

1830 13.3018 12

1840 18.0768 17

1850 24.5659 23

1860 33.3844 31

1870 45.3685 38

1880 61.6547 50

1890 83.7871 62

1900 113.8645 75

1910 154.7390 91

1920 210.2863 105

1930 285.7737 122

1940 388.3592 131

1950 527.7701 151

1960 717.2261 179

1970 974.6918 203

1980 1324.5812 226

1990 1800.0719 249

2000 2446.2516 NA

2010 3324.3932 NA

2020 4517.7651 NA

2030 6139.5269 NA

2040 8343.4597 NA

2050 11338.5478 NA

> t:= 't':

> plot( [ [soln_points], [shifted_data] ], color=[red, blue], style = POINT );

>

Question

How well does the plotted solution (or at least the points we've sampled at 10 year intervals) fit the real word data?. Can you explain why they are not exactly the same?

Question

Determine the values of and when the populations in 1790 and 1890 are used as data for the system of equations. Repeat the calculations using the populations in 1790 and 1990 as data for the system of equations. How do the different choices of data change the quality of the prediction?

>

Logistic Population Growth

The logistic population model is motivated and derived in the text (pp. 8-9). The resulting differential equation is

> MODEL := diff( P(t), t ) = k * ( 1 - P(t)/N ) * P(t);

with the initial condition

> InitCond := P(0) = p[0];

>

Equilibrium Solution

The equilibrium solutions are found in the same manner as in the first example: determine all constant populations (i.e., for which the derivative, dP/dt, and hence the RHS of the model, is zero).

> EQUILeqn := rhs( MODEL ) = 0;

> EQUILsol := solve( EQUILeqn, { P(t) } );

Note that there are now two equilibria. In addition to the zero solution, the carrying capacity of the population, , is an equilibrium solution.

>

Qualitative Analysis

The text gives a good description of the qualitative analysis of the logistic model. Assume > 0. The basic idea is that (for equations of this type - we'll be more specific shortly) the solutions are either always increasing or always decreasing in regions between the equilibria. In this case, for 0 < P < N, the RHS of the model is positive > 0, and so the solutions are increasing (to the carrying capacity). By the same reasoning, for P > N or P < 0, < 0 and solutions are decreasing (to the carrying capacity).

Questions

Analytic Solution

We will learn specific techniques for solving special classes of differential equations. The Maple command for finding explicit solutions to differential equations is dsolve . The arguments are the initial value problem (a set containing the differential equation and the initial condition) and the function for which we are trying to solve.

> SOLNlog := dsolve( { MODEL, InitCond }, P(t) );

>

The U.S. Population

To conclude this investigation of the models, lets determine the values of the parameters that fit the U.S. population in 1790, 1900, and 1990. (Why are three years used for this fit?)

> EQ1790 := eval( SOLNlog, [P(t)=3.9, t=0] );

> EQ1900 := eval( SOLNlog, [P(t)=75, t=110] );

> EQ1990 := eval( SOLNlog, [P(t)=249, t=200] );

>

The solution(s) to this system of three equations and three unknowns can be found using the solve command.

> solve( { EQ1790, EQ1900, EQ1990 }, { k, N, p[0] } );

Note that while there are many solutions, all but one of them involves complex numbers. Since we know our answer should be real-valued (WHY?) we can assume that these are the parameter values that we should use for this problem. There is a fancy way to have Maple select all real-valued solutions but it's simpler to simply copy and paste the portion of the output that interests us.

> PARAMETERS := {p[0] = 3.900000000, k = .2934545722e-1, N = 302.9459184};

(optional) Fancy method of extracting real-valued solution from multiple solutions of a system of equations

> PARAMETERS := op( remove( has, {%}, I ) );

>

The logistic model for the U.S. population based on the 1790, 1900, and 1990 populations is

> solPlog := eval( rhs( SOLNlog ), PARAMETERS );

> soln_points:= NULL:

A table of predicted populations for 1790 - 2050 is easily created as before

> for t from 0 to 260 by 10 do

> printf( `%4d %10.4f %10.4f %10.4a\n`, t+1790, solP, solPlog, shifted_data[t/10+1,2] );

> soln_points:= soln_points, [t, solPlog]:

> od:

1790 3.9000 3.9000 3.9

1800 5.3000 5.2072 5.3

1810 7.2026 6.9425 7.2

1820 9.7881 9.2380 9.6

1830 13.3018 12.2612 12

1840 18.0768 16.2190 17

1850 24.5659 21.3605 23

1860 33.3844 27.9729 31

1870 45.3685 36.3678 38

1880 61.6547 46.8528 50

1890 83.7871 59.6840 62

1900 113.8645 75.0000 75

1910 154.7390 92.7479 91

1920 210.2863 112.6205 105

1930 285.7737 134.0360 122

1940 388.3592 156.1821 131

1950 527.7701 178.1285 151

1960 717.2261 198.9777 179

1970 974.6918 218.0050 203

1980 1324.5812 234.7437 226

1990 1800.0719 249.0000 249

2000 2446.2516 260.8112 NA

2010 3324.3932 270.3747 NA

2020 4517.7651 277.9754 NA

2030 6139.5269 283.9271 NA

2040 8343.4597 288.5338 NA

2050 11338.5478 292.0674 NA

> t := 't':

>

Note that the logistic model does a much better job of following the actual U.S. populations (as reported in Table 1.1 on p. 7). A final comparison of the two solutions can be performed by looking at the plot of the solutions to each model. Here we want to overlay two different kinds of plots (not just two curves in one kind of plot). The way to do this is to give the plots names, end the commands with colons, and then ``display'' the two together

> plotA:= plot( [ solP, solPlog ], t=0..300, 0..400, color=[red, green], title=`Unlimited and Logistic Growth of the US Population (in millions)` ):

> plotB:= plot( [ [shifted_data] ], color = blue, style = POINT):

> display([plotA, plotB]);

Note how one solution appears to level off around 300 (million) and the other appears to continue to grow forever. Which curve corresponds with which model?

>

> restart:

>