kill(all)$

maxima:population models,numeric,part2,version 1.0
peter.vlasschaert@gmail.com,08/04/2022

ref:elementary linear algebra 
      auth: C.H.Edwards,David E.penny
      Prentice Hall

ref: see website
maxima:population models,numeric,part1,version 1.0

packages and option
--------------------------

ratprint:false$

problem :discrete odes "representation".
==============================
odes :(o)rdinary (d)ifferential (e)quations (s)ystem

probblem : move people suburbs to city , city to suburbs.
Given total population 1000000 people,15% people of the city moves to the suburbs
and 10% people moves from suburbs to the city.
--------------------------------------------------------------------------
? how many people after n years are in city are in the suburbs
? What "n" value for steady state .
---------------------------------------------------------------------------

solve : problem
------------------

15% people : city       → suburbs : S[n+1]  : part of total "0.15*C[n]"
10% people :suburbs  → city       : C[n+1] : part of total "0.10* S[n]"

model population : interaction between(city,suburbs)
---------------------------------------------------------------

p1:C[n+1]=(1-0.15)*C[n]+0.10*S[n];
p2:S[n+1]=0.15*C[n]+(1-0.1)*S[n];

p3:[p1,p2];

p3a:sort(listofvars(p3));

p3b1:makelist(p3a[i],i,1,2);
p3b2:makelist(p3a[i],i,2,3);

p31a:(-1)*coefmatrix(p3,p3b1);

augcoefmatrix : {p1,p2}

pa1:augcoefmatrix ([-p1,-p2], [C[n], S[n]]);

pa2:sort(listofvars(p1));
pa3:sort(listofvars(p2));

→ pa,pb  = Transition Matrix of the population model.

pa:(-1)*coefmatrix([p1,p2],[C[n],S[n]]);

pb:float(pa);

intial condition:

first step : calculation, eigenvalues.
------------------------------------------

characteristic equation : find λ

p4:λ*ident(length(pb));

use : pa 

p4a:determinant(pa-p4)=0;

p4a1:expand(p4a);

p4a2:sort(solve(%,listofvars(%)));

second step : calculation, eigenvectors.
---------------------------------------------

p5:transpose([x,y]);

p5a:transpose(makelist(0,i,1,2));

simp:false;

p5b1:(pa-p4).p5=p5a;

simp:true;

p5b2:p5b1;

a) first : eigenvalue.
------------------------

p5b1a:subst([p4a2[1]],p5b1);

p5b1a1:flatten(args(lhs(p5b1a)));
p5b1a2:flatten(args(rhs(p5b1a)));

p5b2aa1:map("=",p5b1a1,p5b1a2);

p5b2aa2:listofvars(p5b2aa1[1]);

p5b3:flatten(solve(p5b2aa1,p5b2aa2));

p5b4a:subst([%r1=1],p5b3);

p5b4b:makelist(rhs(p5b4a[ii]),ii,1,2);

b) second : eigenvalue.
---------------------------

use: p5b2

p6a1:subst([p4a2[2]],p5b1);

p6b1:flatten(args(lhs(p6a1)));
p6b2:flatten(args(rhs(p6a1)));

p6b3:map("=",p6b1,p6b2);

p6b4a:flatten(solve(p6b3,p5b2aa2));
p6b4b:subst([%r2=3],p6b4a);
p6b4c:makelist(rhs(p6b4b[ii]),ii,1,2);

third step : eigenvector "matrix" 

p7a:[p5b4b,p6b4c];

p7b:transpose(apply(matrix,p7a));

P:p7b;

fourth step :inverse(P) = P^(-1)

PP:invert(P);

fifth step : Diagonal matrix:D[i,i],i=1,2

D: ident( 2);

first eigenvalue :D[1,1]

D[1,1]:rhs(p4a2[1]);

second eigenvalue :D[2,2]

D[2,2]:rhs(p4a2[2]);

D:D;

simp:false;

check :A = P.D.P^(-1);

A:P.D.PP;

simp:true;

A:A;

check : against A ? pa

p8:is(A=pa);

set : B =A to find A^k , integer. number of years
----------------------------------------------------------

ref :https://maxima.sourceforge.io/docs/manual/de/maxima_70.html

first : step : load packages
-------------------------------

load("solve_rec")$

second : step : reccurency
-------------------------------

rec:V[k+1]=B*V[k];/*general model: B=matrix */

p10a:solve_rec(rec,V[k]);

initial condition : k =0

p10b:reverse(subst([k=0],p10a));/*V[0]=vectopr*/

V[0]:transpose([700000,300000]);

final : result V[k], k=0,....,m
m ≥ 30, steadystate: see : part 1, "above"

p10c:subst([p10b],p10a);/*V[k]= vector*/

see before : A^k , k = integer (see:check,above)
----------------------------------
rem :A^2 =P.D.P^(-1).P.D.P^(-1)=P.D.I.D.P^(-1)
         A^2 =P.D.D.P^(-1) =P.D^2.P^(-1) 
         -------------------------------------------------
        A^k  =P.D^k.P^(-1)

simp:false;

construct : C^k, C=A=B
-----------------------------

p11a:C^k = P.D^k.PP;

simp:true;

p11aa:'limit((3/4)^k,k,100)=float(limit((3/4)^k,k,100));

p11bb:factor(D^k);

result : very high values,see p11aa
estimate : book . k ≥ 40 : infinity ,reach steady state.

p11cc:subst([p11bb[1,1]=0],p11bb);

C[k]= P.p11cc.PP;

steadysate : population "1 miljoen people"

P.p11cc.PP.V[0];

see also : part 1

evolution problem:plot "dynamics package."
keywords : power matrix,eigenvalue,eigenvector,reccurency.