Introduction

## PART I: FIRST ORDER EQUATIONS Section 3 A Changing Population Heterogeneous in Age James V. Herod*

Web page maintained by Evans M. Harrell, II, harrell@math.gatech.edu.

In this model we keep up with the changes in time of a population that is not homogeneous in age. We suppose that the function P(a,t) is the population density of a species changing in age as time evolves. That is,

P(y,t) [[Delta]]y

represents an approximation of the number of the species at time t in the age bracket [y, y + [[Delta]]y].

This suggests that the number of teenagers in this species at time t would be

I(13,20, ) P(y,t) dy.

We also keep track of the death rate for individuals at time t and age a with D(a,t). With this notation, the number of people that die in the age bracket [52,54] and in the time bracket from t = 92 to t = 96 would be

I(92,96, ) I(52,54, ) D(y,t) P(y,t) dy dt.

For example, the structure of D(y,t) for humans, and as a function of y, initially decreases slightly and then rises as y increases, rising more sharply after y = 75. On the other hand, one would expect that if D were zero at some time t and s > 0, then

P(y,t+s) = P( y-s,t).

If D is not zero,

P(a,t+[[Delta]]t) - P(a-[[Delta]]t,t) ~

- D(a-[[Delta]]t,t) P(a-[[Delta]]t,t) [[Delta]]t.

Adding and subtracting P(a,t) to the left side, and dividing by [[Delta]]t gives

F(P(a,t+[[Delta]]t) - P(a,t),[[Delta]]t) + F(P(a-[[Delta]]t,t) - P(a,t),-[[Delta]]t) ~ - D(a-[[Delta]]t,t) P(a-[[Delta]]t,t).

This leads to the partial differential equation

F([[partialdiff]] ,[[partialdiff]]t)P(y,t) + F([[partialdiff]] ,[[partialdiff]]y)P(y,t) = - D(y,t) P(y,t). (3.1)

The solution for the equation (1) in case y > t is

P(y,t) = P( y-t,0) exp( - I(0,t, )D( y-t+u,u) du ). (3.2)

A verification for this solution can be made with MAPLE. The syntax for the verification is a derivation of the solution follows:

* P:=(y,t)->h(y-t)*exp(-int(De(y-t+u,u),u=0..t));

* diff(P(y,t),t)+diff(P(y,t),y)+De(y,t)*P(y,t);

* simplify(");

We derive a solution for equation (3.2). The solution for the equation (3.2) is a function in the upper-right-quarter plane t > 0, y > 0. The function has continuous partial derivatives connected by the equation (3.2). Suppose that Y(t) = t + [[eta]], where [[eta]] > 0. This defines a family of lines in the right-quarter plane.

Suppose that we had a solution S of (3.2) and that z(t) and h(t) are defined by

z(t) = P(Y(t),t), h(t) = D(Y(t),t)

Then z[[minute]](t) = F([[partialdiff]] ,[[partialdiff]]t)P(Y(t),t) + F([[partialdiff]] ,[[partialdiff]]y)P(Y(t),t) = - h(t) z(t) (3.3)

with

z(0) = P(Y(0),0) = P([[eta]],0).

Thus, z(t) satisfies an ordinary differential equation whose graph lies above the characteristic line Y(t) = t + [[eta]]. The surface that would be the graphs of all such z's as [[eta]] changes would be a surface over the upper-right-quarter plane that could be thought of as a solution surface for (3.2).

The solution of (3.3) is

z(t) = z(0) exp(- I(0,t, ) h(s) ds)

= P([[eta]],0) exp(- I(0,t, ) D(Y(s),s) ds)

= P([[eta]],0) exp(- I(0,t, ) D(s+[[eta]],s) ds)

We have only to unravel the connection of [[eta]], t, and y. If a point {t,y} is chosen in the quarter plane, then the [[eta]] from which the characteristic curve originates that goes through this point {[[eta]],0} is [[eta]] = y - t. Hence, we have that

P(y,t) = P(Y(t),t) = z(t) = P(y-t,0) exp(- I(0,t, ) D(s+y-t,s) ds).

This is the solution (3.2).

Equation (3.2) shows how the already born population progresses in time -- already born because y > t. If y < t, then the population depends on the unborn population at t = 0. If we suppose the size of unborn population, and we suppose that D(a,t) = 0 for a < 0, then equation (3.2) can still be interpreted:

for 0 < y < t

I(0,t, )D( y-t+u,u) du = I(0,t-y, )D( y-t+u,u) du + I(t-y,t, )D( y-t+u,u) du.

Because, if 0 < y < t-y, then y - t < y - t + u < 0 and the first integral is zero. By a change of variables, the second integral can be transformed.

I(0,t-y, ) 0 du + I(0,y, )D( z, z+t-y) dz.

Thus, if 0 < y < t and P(a,0) is known for negative a's,

P(y,t) = P( y-t,0) exp( - I(0,y, )D( z, z+t-y) dz ). (3.4)

A discrete model.

Assume an initial population P1(n), n = 1, 2, ... 100, and a function D that changes with age, but not in t. Also, suppose that

Pn(1) = ISU(i=21,30, )Pn-1(i)/10

and

Pn(i+1) = Pn-1(i) - D(i) Pn-1(i).

To see whether this is a growing population, one could compute

Tot= ISU(1,100, )Pn(i)

for each 1 <= n <= 10.

Here is a way to construct a discrete model of an aging population.

We first make a hypothetical function D. The definition changes at age 22.

* a:=-25/5016: b:= 25/152: c:=-2797/2508:

child:=t->a*t^3+b*t^2+ c*t+3:

* Death:=proc(p) if p < 22 then child(p) else

fi end;

* plot(Death,0..100);

We take a hypothetical population. What follows initializes the first population for the first year. We take the birth to be a function of the population only between ages 21 and 30. Also, it will be of interest if this hypothetical model is increasing or decreasing in size.

* P1:=t->(100-t)*(25+t):plot(P1(t),t=0..100);sum('P1(p)','p'=1..100)

Advance the 2nd year of population

* P2:=proc(p) if p <= 1 then birth1 else

if 1 < p and p < 100 then P1(p-1)*(1-Death(p-1)/1000) else

if p = 100 then 0

fi fi fi end;

sum('P2(p)','p'=1..100);

birth2:=sum('S2(p)','p'=21..30)/10;

* plot(P2,1..100);

Advance the 3rd year of population

* P3:=proc(p) if p <= 1 then birth2 else

if 1 < p and p < 100 then P2(p-1)*(1-Death(p-1)/1000) else

if p = 100 then 0

fi fi fi end;

sum('P3(p)','p'=1..100);

birth3:=sum('P3(p)','p'=21..30)/10;

* plot(P3,1..100);

Advance the 4th year of population

* P4:=proc(p) if p <= 1 then birth3 else

if 1 < p and p < 100 then P3(p-1)*(1-Death(p-1)/1000) else

if p = 100 then 0

fi fi fi end;

sum('P4(p)','p'=1..100);

birth4:=sum('P4(p)','p'=21..30)/10;

* plot(P4,1..100);

Special Case with Birth Rate and other simplifying assumptions: Suppose that P(y,0) is known and that the birth rate is proportional to the total population at time t. That is, there is a number k so that the number of births is

P(0,t) = kI(0,*, ) P(y,t) dy.

Suppose also that D(y,t) = c is constant in t and y. From equation (3.2) with y > t we have how the population changes

P(y,t) = P(y-t,0) e-ct.

Equation (3.4) is another matter. We must predict unborn populations from these simplifying assumptions; we need to determine P(y-t,0) for t > y. Define

u(t) = P(0,t).

Then, P(y,t) = P( y-t,0) = P(0,t-y) = u(t-y).

The simplifying assumptions predict the level of the unborn population. We should be able to determine u(t-y) of all y > 0, t > 0 and, consequently, determine P(y,t).

Here's how to find u:

u(t) = k I(0, *, )P(y,t) dy = k I(0,t, )P(y,t) dy + k I(t,*, ) P(y,t) dy.

With D [[equivalence]] c, and using, first, equation (3) and then (2), these last sum to

k I(0,t, ) P( y-t,0) e-cy dy + k I(t,*, )P( y-t,0) e-ct dy

= k I(0,t, ) u(t-y) e-cy dy + k I(t,*, )P( y-t,0) e-ct dy.

Now, use the change of variable z = y-t for the second integral to get

k I(0,t, ) u(t-y) e-cy dy + k e-ct I(0,*, )P(z,0) dz.

Again, a change of variables of z = t - y changes this last to

k e-ct I(0,t, ) u(z) ecz dz + k e-ct I(0,*, )P(y,0) dy.

In summary,

ect u(t) = k B( I(0,t, ) u(z) ecz dz + I(0,*, )P(y,0) dy).

This is an integral equation in u that must be solved. This may can be done by taking the derivative of both sides. The solution is, for 0 < y < t,

P(y,t) = exp( (k-c)(t-y)) k I(0,*, )P(y,0) dy. (3.3')

Information on death rates in the United States

In the US there is a high infant mortality. Here is an understanding of that

mortality table as data in MAPLE. The units are deaths per 1000 population.

* yrs:=([1, 3, 9, 19, 29, 39, 49, 59, 69, 79, 89]);

DR:=([11.2, .6, .3, 1.5, 1.9, 2.9, 6.5, 16.5, 37, 83.5, 181.9]);

* pts:=[seq([yrs[i],DR[i]], i=1..9)];

* plot(pts,style=POINT);

Curiously, the growth seems exponential after nine years. Here's an attempt to fit the data. First, take the logarithm of the data. Then, plot logarithm of the data; if it looks like a line, the death rate can be approximated with

an exponential growth.

* lnpts:=[seq([yrs[i],ln(DR[i]) ], i=1..9)];

* plot(lnpts,style=POINT);

We find a least squares linear fit to this "near-linear" data.

* with(linalg):

* leastsqrs({a*yrs[1]+b=ln(DR[1]),

a*yrs[2]+b=ln(DR[2]),

a*yrs[3]+b=ln(DR[3]),

a*yrs[4]+b=ln(DR[4]),

a*yrs[5]+b=ln(DR[5]),

a*yrs[6]+b=ln(DR[6]),

a*yrs[7]+b=ln(DR[7]),

a*yrs[8]+b=ln(DR[8]),

a*yrs[9]+b=ln(DR[9])},{a,b});

* assign(");

Now we can make an approximation of the "death-rate" function and compare its graph with the data.

* with(plots):

dpl:=plot(pts,style=POINT):

fpl:=plot(exp(a*x+b),x=9..90):

display({dpl,fpl});