#OK to post homework
#Julian Herman, 9/20/21, Assignment 5

#0)
#   a(n+3) = -6*a(n-1) - 5*a(n+1) substitute n=n-3
#   a(n) = 0*a(n-1) - 5*a(n-2) + 0*a(n-3) - 6*a(n-4) : canonical form

read `/Users/julianherman/Documents/Rutgers/Fall 2021/Dynamical Models In Biology/HW/M5.txt`;

RecToSeq([1, 2, 4, 11], [0, -5, 0, -6], 1000)[1000];
1818014589793496842119263353977165955901169251300081152017301791\

  62903000957919477420992513491076776799335003400559596244171485\

  81612767396466425154660618131176283941650552170945484194399749\

  3283513047867597347184546959401904109745684403540309




#1)
GrowthC([1, 1, 2, 3, 5, 8, 13, 21, 34, 55], [1, 1, 1, 1, 1, 1, 1, 1, 1, 1], 1000);
                          1.999018633

GrowthCe([1, 1, 1, 1, 1, 1, 1, 1, 1, 1]);
                          1.999018633



#2)
Sur := [seq(99/100, i = 0 .. 44)];

Fer := [seq(0, i = 0 .. 15), seq(1/2, i = 16 .. 30), seq(1/4, i = 31 .. 45)];

Rec := LeslieMod(Sur, Fer);

GrowthCe(Rec);
                          1.085998871

M := LeslieMat(Sur, Fer);
with(LinearAlgebra);
evalf(Eigenvalues(Matrix(M)))[1];
                          1.085998871



#3)
GrowthCe([0, 0, 0, 0.16, 0.41]);
                          0.8879729192
#the growth rate of salmons



#4)
PlantGseq := proc(alpha, beta, gamma, sigma, INI, K) option remember; RecToSeq(INI, [alpha*sigma*gamma, beta*sigma^2*(1 - alpha)*gamma], K); end proc;
PlantGseq(0.5, 0.25, 2, 0.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(0.6, 0.3, 2, 0.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 are correct!



#5)
PlantGseq2 := proc(alpha, beta, gamma, sigma) option remember; GrowthCe([alpha*sigma*gamma, beta*sigma^2*(1 - alpha)*gamma]); end proc;
PlantGseq2(0.5, 0.25, 2, 0.8);
                          0.9656854250 #eventual extinction

PlantGseq2(0.1, 0.2, 0.15, 0.8);
                          0.1375902732 #eventual extinction

PlantGseq2(0.05, 0.1, 0.5, 0.3);
                         0.06924093449 #eventual extinction

PlantGseq2(0.5, 0.4, 2, 0.95);
                          1.240914486 #explosion

PlantGseq2(0.6, 0.5, 4, 0.95);
                          2.561829807 #explosion

PlantGseq2(0.6, 0.8, 7, 0.95);
                          4.444821422 #explosion

PlantGseq2(0.5502, 0.30099, 1.9009, 0.799035);
                          1.000000231 #stability