#Maple code for Lecture 5

#DMB:="Dynamical Models in Biology" by S.P. Ellner and J. Gukenheimer

Help5:=proc(): print(`RecToSeq(INI,REC,N), GrowthC(INI,REC,K) , GrowthCe(REC), LeslieM(SUR,FER) `):end:



#RecToSeq(INI,REC,N): Inputs two lists of numbers, INI and REC (of the same length, let's call it k) and a positive integer N larger than their length
#outputs the list of the first N members of the sequence satisfying the linear recurrence with constants coefficients or order k
#f(n)=REC[1]*f(n-1)+...+REC[k]*f(n-k)

RecToSeq:=proc(INI,REC,N) local i,k,L,newguy:
if not (type(INI,list) and type(REC,list) and nops(INI)=nops(REC) and type(N,integer) and N>=nops(INI)) then
 print(`bad innput`):
 RETURN(FAIL):
fi:

k:=nops(INI):

L:=INI:

while nops(L)<N do
 newguy:=add(REC[i]*L[-i],i=1..k):
 L:=[op(L),newguy]:
od:
L:
end:


#GrowthC(INI,REC,K): The estimate of the growth constant of the sequence given by the inital conditions and recurrence pair [INI,REC] using K terms
GrowthC:=proc(INI,REC,K) local L,a,b:
L:=RecToSeq(INI,REC,K):

a:=L[-1]/L[-2]:
b:=L[-2]/L[-3]:

if abs(a-b)<1/10^(Digits+3) then
RETURN(evalf(a)):
else
print(` make `, K, `bigger `):
 RETURN(FAIL):
fi:
end:



#GrowthCe(REC): The  EXACT  growth constant of the recurrence REC using the characteristic equation
GrowthCe:=proc(REC) local x,i:
evalf([solve(1-add(REC[i]/x^i,i=1..nops(REC)))])[1]:

end:


#LeslieMod(SUR,FER):In a population with A age-groups and survival vector 
#(following the notation in the book DMB, p.33
#SUR=[p[0], ...p[A-1]] where p[i] is the probability of somebody of age i will still be alive the next year and
#FER=[f[0],..., f[A]] where f[i] is the fertility factor (of course f[0]=0 in real life, and also in real life the
#younger ages can't have babies, so SUR has A components and FER has A+1 components
#Outputs the recurrence vector REC that enables the computation of the population growth. TRY:
#LelieM([9/10,9/10],[0,1/2,1]]);
LeslieM:=proc(SUR,FER) local i,L,A:

 if not (type(SUR,list) and type(FER,list) and nops(SUR)+1=nops(FER)) then
  print(`bad input`):
 RETURN(FAIL):
 fi:

A:=nops(SUR):

L[0]:=1:

for i from 1 to A do
 L[i]:=L[i-1]*SUR[i]:
od:

[seq(FER[i]*L[i],i=1..A)]:

end: