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

#1.)
f := (2 + x)/(4 + x);
                                2 + x
                           f := -----
                                4 + x

a := Orb(f, x, 1.0, 1, 1000);
a := [0.6000000000, 0.5652173913, 0.5619047619, 0.5615866388,

  0.5615560641, 0.5615531253, 0.5615528428, 0.5615528157,

  0.5615528131, 0.5615528128, [...981 terms...], 0.5615528128,

  0.5615528128, 0.5615528128, 0.5615528128, 0.5615528128,

  0.5615528128, 0.5615528128, 0.5615528128, 0.5615528128,

  0.5615528128]

#There is a steady state at x=0.5615528128
# x = (2+x) / (4+x)
# x^2+4x-x-2=0
# x^2+3x-2=0
# using fsolve yields the expected fixed points:
fsolve(x^2 + 3*x - 2 = 0);
                   -3.561552813, 0.5615528128

# x=0.5615528128 is a stable fixed point.
SFP(f, x);
                         [0.561552813]

#abs(f'(0.5615528128)) < 1 ?
abs(subs(x = 0.5615528128, diff(f, x)));
                         0.0961179680
#x=0.5615528128 IS STABLE!

#f'(x)= ((4+x)-(2+x)) / (4+x)^2
#abs(f'(-3.561552813)) < 1 ?
abs(subs(x = -3.561552813, diff(f, x)));
                          10.40388204
#x=-3.561552813 IS NOT STABLE!




#2.)
Orb(x*(1 - x), x, 0.5, 1, 1000);
 [0.25, 0.1875, 0.15234375, 0.1291351318, 0.1124592495,

   0.09981216670, 0.08984969808, 0.08177672983, 0.07508929629,

   0.06945089387, [...981 terms...], 0.0009993245315,

   0.0009983258820, 0.0009973292274, 0.0009963345618,

   0.0009953418792, 0.0009943511737, 0.0009933624394,

   0.0009923756705, 0.0009913908610, 0.0009904080051]
#When k=1, the orbit approaches 0 as more terms are computed.


Orb(2*x*(1 - x), x, 0.5, 1, 1000);
[0.50, 0.5000, 0.50000000, 0.5000000000, 0.5000000000,

  0.5000000000, 0.5000000000, 0.5000000000, 0.5000000000,

  0.5000000000, [...981 terms...], 0.5000000000, 0.5000000000,

  0.5000000000, 0.5000000000, 0.5000000000, 0.5000000000,

  0.5000000000, 0.5000000000, 0.5000000000, 0.5000000000]
#When k=2, the orbit converges to the fixed point 0.5.


Orb(2.5*x*(1 - x), x, 0.5, 1, 1000);
[0.625, 0.5859375, 0.6065368652, 0.5966247410, 0.6016591485,

  0.5991635438, 0.6004164790, 0.5997913268, 0.6001042278,

  0.5999478590, [...981 terms...], 0.6000000000, 0.6000000000,

  0.6000000000, 0.6000000000, 0.6000000000, 0.6000000000,

  0.6000000000, 0.6000000000, 0.6000000000, 0.6000000000]
#When k=2.5, the orbit converges to the fixed point 0.6.


Orb(3.1*x*(1 - x), x, 0.5, 1, 1000);
[0.775, 0.5405625, 0.7698995192, 0.5491781734, 0.7675026726,

  0.5531711926, 0.7662357553, 0.5552674201, 0.7655310881,

  0.5564290480, [...981 terms...], 0.5580141245, 0.7645665203,

  0.5580141245, 0.7645665203, 0.5580141245, 0.7645665203,

  0.5580141245, 0.7645665203, 0.5580141245, 0.7645665203]
#When k=3.1, the orbit oscillates between two values: 0.5580141245, 0.7645665203


Orb(3.5*x*(1 - x), x, 0.5, 1, 1000);
[0.875, 0.3828125, 0.8269348143, 0.5008976952, 0.8749971794,

  0.3828199039, 0.8269408878, 0.5008837956, 0.8749972662,

  0.3828196760, [...981 terms...], 0.5008842111, 0.8749972637,

  0.3828196827, 0.8269407062, 0.5008842111, 0.8749972637,

  0.3828196827, 0.8269407062, 0.5008842111, 0.8749972637]
#When k=3.5, the orbit oscillates between four values: 0.3828196827, 0.8269407062, 0.5008842111, 0.8749972637




#3.)
a := rsolve({a(0) = 0.5, a(1) = 0.7, a(n) = (2*a(n - 1) + a(n - 2))/(a(n - 1) + 7*a(n - 2))}, a(n), 'makeproc');
[seq(a(i), i = 0 .. 1000)];
[0.5, 0.7, 0.4523809524, 0.2998220641, 0.3034843505,

  0.3774773620, 0.4230595506, 0.3991635773, 0.3634451637,

  0.3566181271, [...981 terms...], 0.3750000000, 0.3749999999,

  0.3750000000, 0.3750000001, 0.3750000000, 0.3749999999,

  0.3750000000, 0.3750000001, 0.3750000000, 0.3749999999]
#It oscillates between the three values (hitting 0.3750000000 twice) with a period of 4:
0.3750000000, 0.3750000001, 0.3750000000, 0.3749999999
#This could be because of rounding in Maple, but I think not. IF it is, the steady state
#would be 0.3750000000.