kill(all)$

maxima:population models,numeric,part1,version 1.0
peter.vlasschaert@gmail.com,27/01/2022

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

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];

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);

→ steady state model:equation

what to do:
----------------------------------------------------
1e step : p1, C[n+1]→C[n]; p2,S[n+1]→S[n]
2e step : constraint : C[n]+S[n]=1
-----------------------------------------------------

p3a:reverse(C[n]+S[n]=1);/*constraint */
p3b:C[n]=rhs(p1);
p3c:S[n] =rhs(p2);

p4:listofvars(p3a);

augcoefmatrix ([-p3a,-p3b,-p3c], p4);

steady state : result

p5:solve([p3a,p3b,p3c],p4)[1];

total population ; urban (=(C)ity) and suburban ((S)uburbs)

ptot: 1000 ; /*thousand people is total population*/

p5a:makelist(lhs(p5[i])=ptot*rhs(p5[i]),i,1,length(p4));

how to find : n for steady state (n = numbero of years).

initial population : (initial population : @ 
some time,use index:0)

C[0]:700;
S[0]:300;

first code : simple

for n:0 thru 40 do(
    C[n+1]:0.1*S[n]+0.85*C[n],
    S[n+1]:0.9*S[n]+0.15*C[n],
    C[n]:C[n+1],
    S[n]:S[n+1]
)$

rem : other

r1:[C[29],S[29]];

r2:[C[30],S[30]];

r3:[C[31],S[31]];

r4:[C[32],S[32]];

rem : only condition: C[0]+S[0] = ptot , see : above.
example : ptot = 1000 

steady state

"C[0]+S[0] = total population (every combination allowed)"
!  check yourself : "see:totalpopulation = ptot"
check : C[0]=500, S[0]=500 ,same result for steadystate.

first method : guess
------------------------

find steady state: 
1e) take very large : n  
      (example : n=40)
2e)  value before approximated value after it 
      (C[n+1]≅ C[n],S[n+1]≅S[n])

numerical values : numer,fpprec,fpprintprec.
-----------------------------------------------------

numer:true$

fpprintprec : 4;

fpprec:4;

for i:1 thru 32 do(
    print(i,C[i],S[i])
)
;

other way : list

lr1:makelist([i,C[i],S[i]],i,1,32);

analyse data :lr1

lr1[30];

[lr1[30][1],lr1[30][2],lr1[30][3]];

result :steady state : C[0]=700,S[0] = 300, n≥ 30 years

second method : use p5a
-----------------------------

first code :
for n:0 thru 40 do(
    C[n+1]:0.1*S[n]+0.85*C[n],
    S[n+1]:0.9*S[n]+0.15*C[n],
    C[n]:C[n+1],
    S[n]:S[n+1]
)

p5b:p5a;

p5c:[C[n]=400.0,S[n]=600.0];

second code : doesn't work here.

pv1a:(rhs(p5c[1]));

pv1b:(lr1[30][2]);

1e) part : numeric values .

is(pv1a=pv1b);

2e) part: fixed values,numeric.

is(400.0=400.0);

rem: how to get break out "loop,for".
-------------------------------------------

→ use : block statement

cp:for i:1 step 1 thru 7 do
block(if(i≥6) then 
  exit
    else
 ( n:i,
  print(i,"is less than 6"))
)$
print(n+1)$

rem : store i value in "n"→ n:i,see:code above.

→ use:without block statement

cp1: for i:1 step 1 thru 7 
      do if(i≥6) then return(i) else print(i,"is less than 6");

-----------------------------------------------------------------------