#OK to post homework
#Julian Herman, 10/11/21, Assignment 11

#1)
f = k*x*(1 - x);
SFPe(f, x);
                   [        [k - 1        ]]
                   [[0, k], [-----, -k + 2]]
                   [        [  k          ]]

#The SFPe(f,x) function outputs a list of lists, where each of the nested lists
#consists of a fixed point in the first position, followed by the value of the
#derivative of f evaluated at that fixed point.

#The fixed point is a STABLE fixed point if the absolute value of the derivative
#of f evaluated at the fixed point is less than 1.
#Solving for the first fixed point, x=0:
#x=0 is stable if |k|<1 ... -1<k<1

#Solving for the second fixed point, x=(k-1)/k:
#      |-k+2|<1
# -1<-k+2   -k+2<1
#   3>k       k>1
#       1<k<3
#THEREFORE, x=(k-1)/k is a STABLE fixed point for 1<k<3...
#For EACH value of k, there is only one stable fixed point... for k from 0 to 1
#it is 0, for k from 1 to 3, it is (k-1)/k...
#EXAMPLE: when k=2: x=(2-1)/2 = 1/2 is THE STABLE fixed point...
#Assuming k cannot be less than 1 (because then the population will go extinct,
#which means that x=0 is a STABLE fixed point, but this is irrelevant) AND that
#k must be an integer, the only STABLE fixed point for k<3 is x=1/2.




#2)i)
Digits := 2;

f := 3.1*x*(1 - x);
convert(Orb(f, x, 0.5, 1000, 1020), set);
                          {0.53, 0.78}

f := 3.2*x*(1 - x);
convert(Orb(f, x, 0.5, 1000, 1020), set);
                          {0.51, 0.80}

f := 3.3*x*(1 - x);
convert(Orb(f, x, 0.5, 1000, 1020), set);
                          {0.50, 0.82}

f := 3.4*x*(1 - x);
convert(Orb(f, x, 0.5, 1000, 1020), set);
                          {0.44, 0.85}

f := 3.5*x*(1 - x);
convert(Orb(f, x, 0.5, 1000, 1020), set);
                    {0.38, 0.46, 0.84, 0.88}

f := 3.49*x*(1 - x);
convert(Orb(f, x, 0.5, 1000, 1020), set);
                    {0.38, 0.46, 0.84, 0.88}

f := 3.48*x*(1 - x);
convert(Orb(f, x, 0.5, 1000, 1020), set);
                    {0.38, 0.46, 0.84, 0.88}

f := 3.47*x*(1 - x);
convert(Orb(f, x, 0.5, 1000, 1020), set);
                    {0.38, 0.46, 0.84, 0.88}

f := 3.46*x*(1 - x);
convert(Orb(f, x, 0.5, 1000, 1020), set);
                    {0.38, 0.46, 0.84, 0.88}

f := 3.45*x*(1 - x);
convert(Orb(f, x, 0.5, 1000, 1020), set);
                          {0.44, 0.85}

f := 3.451*x*(1 - x);
convert(Orb(f, x, 0.5, 1000, 1020), set);
                    {0.38, 0.46, 0.84, 0.88}

#THE SECOND PERIOD-DOUBLING EVENT / SECOND BIFURCATION POINT OCCURS AT ABOUT k=3.45

#ii)
f := Comp(k*x*(1 - x), x);
a := SFPe(f, x)[3][2];
solve(abs(a) = 1, k);
                      (1/2)       (1/2)
                 1 - 6     , 1 + 6     , -1, 3

#k=1+sqrt(6) is the EXACT second bifurcation point.

Digits := 20:

evalf(1 + sqrt(6));
                     3.4494897427831780982




#3)
#Model by Varley, Gradwell, and Hassell (1973):

#NOTE: In this model, using N in place of x

f := N^(-b)*lambda*N/alpha;
                            (-b)
                           N     lambda N
                      f := --------------
                               alpha

a := SFPe(f, N);
     [
     [
     [
     [
     [
     [
a := [[0, 0],
     [

  [                                        (-b)                ]]
  [                     /   /    /alpha \\\                    ]]
  [                     |   |  ln|------|||                    ]]
  [   /    /alpha \\    |   |    \lambda/||                    ]]
  [   |  ln|------||    |exp|- ----------||     lambda (-1 + b)]]
  [   |    \lambda/|    \   \      b     //                    ]]
  [exp|- ----------|, - ---------------------------------------]]
  [   \      b     /                     alpha                 ]]


solve(abs(a[2][2]) < 1, b);
Warning, solutions may have been lost
                     RealRange(1, Open(2))

#N=0 is a fixed point according to SFPe(), but I don't think it exists because
#you must divide by 0, which is undefined...

#N=exp(-ln(alpha/lambda)/b) is a stable fixed point: if and only if 1 <= b < 2
#OR (the same): 0<b<2.


#The second model:
f := N*exp(r*(1 - N/K));
                               /  /    N\\
                     f := N exp|r |1 - -||
                               \  \    K//

a := SFPe(f, N);
           a := [[0, exp(r)], [K, exp(0) - r exp(0)]]

#N=0 is a fixed point. It is stable for -infinity < r < 0... Since r cannot be
#negative (it is defined as a POSITIVE constant), N=0 is not stable on any
#defined range.

#N=K is a fixed point. It is stable when |1-r| < 1 OR 0<r<2.
#-1<1-r... 2>r
#1-r<1... r>0


#The third model proposed by Hassell (1975):
#I'm not sure about this one because SFPe() returns _Z and RootOf:
f := lambda*N(N*a + 1)^(-b);
                                         (-b)
                   f := lambda N(N a + 1)

c := SFPe(f, N);
         [
         [[      /                    (-b)     \    /
    c := [[RootOf\-lambda _Z(_Z a + 1)     + _Z/, - \lambda
         [

            /                    (-b)     \/        /
      RootOf\-lambda _Z(_Z a + 1)     + _Z/\a RootOf\
                        (-b)     \    \
    -lambda _Z(_Z a + 1)     + _Z/ + 1/^(-b) b

       /      /                    (-b)     \\/        /
      D\RootOf\-lambda _Z(_Z a + 1)     + _Z//\a RootOf\
                        (-b)     \    \  \//
    -lambda _Z(_Z a + 1)     + _Z/ + 1/ a/ \

            /                    (-b)     \/        /
      RootOf\-lambda _Z(_Z a + 1)     + _Z/\a RootOf\
                        (-b)     \    \\]
    -lambda _Z(_Z a + 1)     + _Z/ + 1//],

      [                (-b)            ]]
      [     lambda 0(1)     b D(0)(1) a]]
      [0, - ---------------------------]]
      [                0(1)            ]]




#4)
f := (z[1] + z[2])*1/(z[1] + z[2]);
                             f := 1

convert(Orbk(2, z, f, [1.1, 5.3], 10000, 10020), set);
                              {1}
#a,b=1,1                 EQ. point=1
#explicit: (1+a)/(1+b) = 2/2 = 1

f := (z[1] + z[2])/(2*z[1] + z[2]);
convert(Orbk(2, z, f, [1.1, 5.3], 10000, 10020), set);
                         {0.667, 0.670}
#a,b=1,2                 EQ. point = about 0.67
#explicit: (1+a)/(1+b) = 2/3 = 0.67

f := (z[1] + z[2])/(3*z[1] + z[2]);
convert(Orbk(2, z, f, [1.1, 5.3], 10000, 10020), set);
                         {0.498, 0.503}
#a,b=1,3                 EQ. point = about 0.5
#explicit: (1+a)/(1+b) = 2/4 = 0.5

f := (z[1] + z[2])/(4*z[1] + z[2]);
convert(Orbk(2, z, f, [1.1, 5.3], 10000, 10020), set);
                            {0.400}
#a,b=1,4                 EQ. point = about 0.4
#explicit: (1+a)/(1+b) = 2/5 = 0.4
#and so on for the below...

f := (z[1] + 2*z[2])/(z[1] + z[2]);
convert(Orbk(2, z, f, [1.1, 5.3], 10000, 10020), set);
                             {1.50}

f := (z[1] + 2*z[2])/(2*z[1] + z[2]);
convert(Orbk(2, z, f, [1.1, 5.3], 10000, 10020), set);
                             {1.00}

f := (z[1] + 2*z[2])/(3*z[1] + z[2]);
convert(Orbk(2, z, f, [1.1, 5.3], 10000, 10020), set);
                     {0.748, 0.753, 0.756}
#a,b=2,3            EQ. point = about 0.75
#explicit: (1+a)/(1+b) = 3/4 = 0.75

f := (z[1] + 2*z[2])/(4*z[1] + z[2]);
convert(Orbk(2, z, f, [1.1, 5.3], 10000, 10020), set);
                         {0.591, 0.609}

f := (z[1] + 3*z[2])/(z[1] + z[2]);
convert(Orbk(2, z, f, [1.1, 5.3], 10000, 10020), set);
                             {2.00}

f := (z[1] + 3*z[2])/(2*z[1] + z[2]);
convert(Orbk(2, z, f, [1.1, 5.3], 10000, 10020), set);
                          {1.32, 1.35}
#a,b=3,2                 EQ. point = about 1.33
#explicit: (1+a)/(1+b) = 4/3 = 1.33


f := (z[1] + 3*z[2])/(3*z[1] + z[2]);
convert(Orbk(2, z, f, [1.1, 5.3], 10000, 10020), set);
                         {0.808, 1.24}
#a,b=3,3               EQ. point = about 1? (taking average)
#explicit: (1+a)/(1+b) = 4/4 = 1

f := (z[1] + 3*z[2])/(4*z[1] + z[2]);
convert(Orbk(2, z, f, [1.1, 5.3], 10000, 10020), set);
                         {0.418, 1.60}
#a,b=3,4               EQ. point = about 1? (taking average)
#explicit: (1+a)/(1+b) = 4/5 = 0.8

f := (z[1] + 4*z[2])/(z[1] + z[2]);
convert(Orbk(2, z, f, [1.1, 5.3], 10000, 10020), set);
                             {2.50}

f := (z[1] + 4*z[2])/(2*z[1] + z[2]);
convert(Orbk(2, z, f, [1.1, 5.3], 10000, 10020), set);
                          {1.62, 1.71}

f := (z[1] + 4*z[2])/(3*z[1] + z[2]);
convert(Orbk(2, z, f, [1.1, 5.3], 10000, 10020), set);
                         {0.628, 2.39}
#a,b=4,3               EQ. point = ?
#explicit: (1+a)/(1+b) = 5/4 = 1.25

f := (z[1] + 4*z[2])/(4*z[1] + z[2]);
convert(Orbk(2, z, f, [1.1, 5.3], 10000, 10020), set);
                   {0.377, 0.378, 2.64, 2.65}
#a,b=4,4               EQ. point = ?
#explicit: (1+a)/(1+b) = 5/5 = 1