age_structured.html

MTH 455 A DISCRETE LINEAR AGE-STRUCTURED MODEL
JR 2/1/2005 age_structured.mws
See Haberman, Section 35, p138.

> # SET UP THE MATRIX A
restart:
with(LinearAlgebra):
M:=3;
b:=[0,1.7,0.5,0];       # birth rates
d:=[0.2,0.1,0.6,0.9]; # death rates
A:=[seq([seq(0,col=0..M)],row=0..M)]: #fill A with zeros
for col from 0 to M-1 do
  A[1,col+1]:=b[col+1];
  A[col+2,col+1]:=1-d[col+1];
end do:
evalm(A);

> # SOME USEFUL PROCEDURES

# MUTLIPLY VECTOR N BY MATRIX A
apply_A:=(A,N)->convert(evalm(A&*N),list):

# NORMALIZE A VECTOR
normaliz:=proc(N)
  local tot,i,Nlocal;
  Nlocal:=N;
  tot:=0;
  for i from 1 to nops(N) do tot:=tot+abs(N[i]); end do:
  unassign('i'):
  for i from 1 to nops(N) do Nlocal[i]:=N[i]/tot; end do:
#   for i from 1 to nops(N) do Nlocal[i]:=N[i]/N[1]; end do:
  return Nlocal;
end proc:

# MAKE A HISTOGRAM OF A VECTOR
my_histogram:=proc(N)
  local n,i,r;
  n:=nops(N);
  for i from 0 to n-1 do
    r[i]:=plottools[rectangle]([i,N[i+1]],[i+1,0],color=pink);
  end do:
  return plots[display](seq(r[i],i=0..n-1));
end proc:

> N:=[1,0,0,0]; # WHAT HAPPENS IF WE START WITH JUST BABIES?
for t from 1 to 30 do
N:=apply_A(A,N);
end do;

> # SAME EXPERIMENT, BUT MAKE AN ANIMATED HISTOGRAM THIS TIME
N:=[1,0,0,0];
histo[0]:=my_histogram(convert(N,list)):
for t from 1 to 30 do
N:=apply_A(A,N);
histo[t]:=my_histogram(convert(N,list)):
end do: unassign('t'):
plots[display]( seq( histo[t], t=0..30 ), insequence=true );

> # SAME EXPERIMENT, BUT THIS TIME NORMALIZE POPULATION
# IN ORDER TO SEE BETTER WHAT'S HAPPENING TO THE AGE-DISTRIBUTION
N:=[1,0,0,0];
histo[0]:=my_histogram(convert(N,list)):
for t from 1 to 30 do
  N:=apply_A(A,N);
  N:=normaliz(N);
  print(N);
  histo[t]:=my_histogram(convert(N,list)):
end do: unassign('t'):
plots[display]( seq( histo[t], t=0..30 ), insequence=true );

THE RELEVANCE OF THE EIGENDATA OF MATRIX A

> Am:=convert(A,Matrix); # "Eigenvectors" needs input of type "Matrix"
#Eigenvalues(Am);
(vals,vecs):=Eigenvectors(Am):
eigvals:=vals;
eigvecs:=evalm(vecs);

> C1:=Column(vecs,1); # THE DOMINANT EIGENVECTOR
# COMPARE WITH SIMULATION (SCALE SO FIRST ELEMENT IS 1)
C1/C1[1];
N/N[1]; # N STORES THE POPULATION AT THE END OF OUR SIMULATION

OBSERVE THAT THE AGE DISTRIBUTION OF THE POPULATION IS SETTLING INTO THAT
GIVEN BY THE DOMINANT EIGENVECTOR OF MATRIX A.

> ?Matrix

> ?Vector

> ?Column