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:
adult:=t->95*t*(t-70)/4056 + 5025/169:
* Death:=proc(p) if p < 22 then child(p) else
adult(p)
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});