#OK to post Homework 
#Jeton Hida, September 20, 2021, Assignment 5
> read "/Users/jeton/Desktop/Math 336/M5.txt"
> Help5()
           RecToSeq(INI,REC,N), GrowthC(INI,REC,K) , GrowthCe(REC)
           LeslieMod(SUR,FER): e.g. LeslieMod([9/10,9/10],[0,1,1]);
           LeslieMat(SUR,FER); e.g. LeslieMat([9/10,9/10],[0,1,1]);
> #Number 0 eq:={6*a(n-1)+a(n+3)+5*a(n+1)=0};

> subs(n=n-3,eq)
                     {6 a(n - 4) + a(n) + 5 a(n - 2) = 0}
> #Cannonical form would be a(n)=0*a(n-1)-5*a(n-2)+0*a(n-3)-6*a(n-4)
> print(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 

   nops(INI) <= N) then print(`bad innput`); RETURN(FAIL) end if; k := 

   nops(INI); L := INI; while nops(L) < N do newguy := add(REC[i]*L[-i], i = 

   1 .. k); L := [op(L), newguy] end do; L end proc
> RecToSeq([1,2,4,11],[0,-5,0,-6],1000)[1000];
181801458979349684211926335397716595590116925130008115201730179162903000957919\

  4774209925134910767767993350034005595962441714858161276739646642515466061813\

  1176283941650552170945484194399749328351304786759734718454695940190410974568\

  4403540309
> #Number 1
> print(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) end if 

    end proc
> GrowthC([0,1,1,2,3,5,8,13,21],[1,1,1,1,1,1,1,1,1],100)
                                 1.998029470
> print(GrowthCe)
proc (REC) local x, i; evalf([solve(1-add(REC[i]/x^i, i = 1 .. nops(REC)))])[1\

  ] end proc
> GrowthCe([1,1,1,1,1,1,1,1,1])
                                 1.998029470
> #Number 2
> print(LeslieMod)
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) end if; 

   A := nops(SUR); L[0] := 1; for i to A do L[i] := L[i-1]*SUR[i] end do; 

   [seq(FER[i+1]*L[i], i = 0 .. A)] end proc
> 
> LeslieMod([.99,.99,.99,.99,.99,.99,.99,.99,.99,.99,.99,.99,.99,.99,.99,.99,.99
> ,.99,.99,.99,.99,.99,.99,.99,.99,.99,.99,.99,.99],[.5,.5,.5,.5,.5,.5,.5,.5,.5,.5
> ,.5,.5,.5,.5,.5,.25,.25,.25,.25,.25,.25,.25,.25,.25,.25,.25,.25,.25,.25,.25])
  [0.5, 0.495, 0.49005, 0.4851495, 0.480298005, 0.4754950250, 0.4707400747, 

    0.4660326740, 0.4613723472, 0.4567586238, 0.4521910375, 0.4476691271, 

    0.4431924358, 0.4387605115, 0.4343729064, 0.2150145887, 0.2128644428, 

    0.2107357984, 0.2086284404, 0.2065421560, 0.2044767344, 0.2024319671, 

    0.2004076474, 0.1984035710, 0.1964195352, 0.1944553399, 0.1925107865, 

    0.1905856786, 0.1886798218, 0.1867930236]
> GrowthCe(%);
                                 1.489452967
> print(LeslieMat)
proc (SUR, FER) local i, A; if not (type(SUR, list) and type(FER, list) and 

   nops(SUR)+1 = nops(FER)) then print(`bad input`); RETURN(FAIL) end if; A 

   := nops(SUR); matrix([FER, seq([0$i-1, SUR[i], 0$A+1-i], i = 1 .. A)]) end 

   proc
> 
> LeslieMat([.99,.99,.99,.99,.99,.99,.99,.99,.99,.99,.99,.99,.99,.99,.99,.99,.99
> ,.99,.99,.99,.99,.99,.99,.99,.99,.99,.99,.99,.99],[.5,.5,.5,.5,.5,.5,.5,.5,.5,.5
> ,.5,.5,.5,.5,.5,.25,.25,.25,.25,.25,.25,.25,.25,.25,.25,.25,.25,.25,.25,.25])

> 
Loading LinearAlgebra
> with(LinearAlgebra):
> 
> Eigenvalues(LeslieMat([.99,.99,.99,.99,.99,.99,.99,.99,.99,.99,.99,.99,.99,.99
> ,.99,.99,.99,.99,.99,.99,.99,.99,.99,.99,.99,.99,.99,.99,.99],[.5,.5,.5,.5,.5,.5
> ,.5,.5,.5,.5,.5,.5,.5,.5,.5,.25,.25,.25,.25,.25,.25,.25,.25,.25,.25,.25,.25,.25,
> .25,.25]))

> abs(1.48945296740263 + 0.*I)
                                 1.489452967
> #Number 3
> #S(t)=.16*S(t-4)+.41*S(t-5) GrowthCe([0,0,0,.16,.41])
                                 0.8879729192
> #Number 4
> PlantGseq:=proc(alpha,beta,gamma,sigma,INI,K) local a, b, c, D,REC, newplant:
if not type(INI,list) then RETURN (FAIL): 
fi: 
a:=(alpha*sigma*gamma):
b:=(beta*sigma*sigma*(1-alpha)*gamma): 
c:=nops(INI): 
D:=INI: 
REC:=[a,b]: 
while nops(D)<K do: 
newplant:=add(REC[i]*D[-i],i=1..c): 
D:=[op(D),newplant]: 
od: 
D:
end:
> PlantGseq(.5,.25,2,.8,[100,80],21)
[100, 80, 80.00000, 76.8000000, 74.24000000, 71.68000000, 69.22240000, 

  66.84672000, 64.55296000, 62.33784320, 60.19874816, 58.13305344, 56.13824246, 

  54.21188252, 52.35162481, 50.55520105, 48.82042081, 47.14516882, 45.52740239, 

  43.96514892, 42.45650352]
> PlantGseq(.6,.3,2,.8,[100,96],21)
[100, 96, 107.5200, 117.964800, 129.7612800, 142.6902221, 156.9139458, 

  172.5546061, 189.7544040, 208.6686153, 229.4681472, 252.3409206, 277.4935912, 

  305.1534130, 335.5702921, 369.0190446, 405.8018797, 446.2511298, 490.7322533, 

  539.6471367, 593.4377253]
> #The two tables on p 18 are correct as according to the program
> #Number 5 
PlantGseq:=proc(alpha,beta,gamma,sigma) local a,b,REC,x:
a:=(alpha*sigma*gamma): 
b:=(beta*sigma*sigma*(1-alpha)*gamma): 
REC:=[a,b]:
evalf([solve(1-add(REC[i]/x^i,i=1..nops(REC)))])[1]: 
end:

> PlantGseq(.5,.25,2,.8)
                                 0.9656854250
> PlantGseq(.7,.33,1,.9)
                                 0.7385740785
> PlantGseq(.20,.70,2,.99)
                                 1.264262632
> PlantGseq(.5,.25,2,1)
                                 1.207106781
> PlantGseq(.5,1,1,1)
                                      1.
> PlantGseq(.25,1,1,1)
                                      1.
> PlantGseq(.5,1,1,.9)
                                 0.9000000000
> PlantGseq(1,.5,.5,.8)
                                 0.4000000000
> #Growth constants equal to 1 are stable, >1 are explosions, <1 are
> extinctions.