NULL;

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

#Conjecture 4: For any non-negative initial conditions, every positive solution of the below 3rd-order rational difference equation converges to a period-six solution of the form: [phi, psi, psi/phi, 1/phi, 1/psi, phi/psi];

;

NULL;

F := (1 + z[1])/(z[2] + z[3]);

NULL;

#The below visualizes how the periodic orbit is achieved within the first few iterations

Orbk(3, z, F, [0.5, 0.5, 0.9], 1, 50);

dataplot(%);

Orbk(3, z, F, [0.5, 0.5, 0.9], 1004, 1009);

#For different initial conditions, the orbit still converges to a period-six solution of the same form: [phi, psi, psi/phi, 1/phi, 1/psi, phi/psi];, but of different values from the above.

Orbk(3, z, F, [5.1, 0.1, 0.9], 1, 50);

dataplot(%);

rollerCoaster := Orbk(3, z, F, [5.1, 0.1, 0.9], 1005, 1010);

#Below we show that the third term in the above list 'rollerCoaster' is equal to the second term divided by first term and so on...

rollerCoaster[1] = rollerCoaster[6]/rollerCoaster[5];

rollerCoaster[2] = rollerCoaster[1]/rollerCoaster[6];

rollerCoaster[3] = rollerCoaster[2]/rollerCoaster[1];

rollerCoaster[4] = rollerCoaster[3]/rollerCoaster[2];

rollerCoaster[5] = rollerCoaster[4]/rollerCoaster[3];

rollerCoaster[6] = rollerCoaster[5]/rollerCoaster[4];

#Assuming that this 3rd order difference equation always converges to a period of 6, then if you convert the underlying transformation F := (1 + z[1])/(z[2] + z[3]); into a system of three unknowns and compose that system with itself six times, you would expect for that composed system to have a fixed point. This fixed point corresponds to each of the unique values in the period as they repeat themselves every 6 iterations of the original transformation F := (1 + z[1])/(z[2] + z[3]);after it has converged to said period.

#NOTE: The underlying transformation F itself CANNOT technically have a fixed point because it maps from R^3 to R... HOWEVER, Fcomp6 DOES have a fixed point because it maps from R^3 to R^3

Help(ToSys);

Fsys := ToSys(3, z, F);

Help(ComK);

Fcomp6 := ComK(Fsys, 6);

Help(FP);

FP(Fcomp6, [z[1], z[2], z[3]])[1];

#This is true for z[1] being any of the values and proceeding in sequential order: [phi, psi, psi/phi, 1/phi, 1/psi, phi/psi];

#For example: ifz[1] = phi, z[2] = psi, z[3] = z[2]/z[1] and z[2]/z[1] = psi/phi .. ();and so on...

#Another example: ifz[1] = 1/phi, z[2] = 1/psi, z[3] = z[2]/z[1] and z[2]/z[1] = phi/psi and phi/psi = phi/psi;

#Below we will demonstrate the part of the conjecture that states that every positive solution of the difference equation:x(n)=(1+x(n-1))/(x(n-2)+x(n-3)) , ;converges to a positive solution of the difference equation: y(n+1)=(y(n))/(y(n-1 )), ;which is the same as y(n) = y(n - 1)/y(n - 2);

orbOriginal := Orbk(3, z, F, [11.2, 1.3, 6.3], 1, 100);

dataplot(orbOriginal);

{op(Orbk(3, z, F, [11.2, 1.3, 6.3], 1000, 1100))};

Y := z[1]/z[2];

orbNew := Orbk(2, z, Y, [0.1540079945, 0.2224474660], 1, 100);

dataplot(orbNew);

{op(Orbk(2, z, Y, [0.1540079945, 0.2224474660], 1000, 1100))};

#Lets try starting Y's orbit with the same initial conditions as we did for F's.

orbNew := Orbk(2, z, Y, [11.2, 1.3], 1, 100);

dataplot(orbNew);

{op(Orbk(2, z, Y, [11.2, 1.3], 1000, 1100))};

#NOTE: Starting both with the same initial conditions does NOT result in achieving the same period six orbit! This is shown above.

#The above is verified by the fixed point of Fcomp6 being [z[1],z[2],(z[2])/(z[1])], ;therefore, for any two sequential values of the period six orbit obtained from F: [`x__1`,`x__2`],;if we use them as initial conditions for Y, then Y's ultimate orbit will automatically fall within the same period because the rule for computing the next value is the same:x(n - 1)/x(n - 2);.

#Below is a jacobian analysis to determine stability of the only positive fixed point of Fsys:

Fsys;

FP(Fsys);

NULL;

SFP(Fsys);

Help(JAC);

j := JAC(Fsys);

jSpecific := subs({z[1] = 1.0, z[2] = 1.0, z[3] = 1.0}, j);

Help(IsDisStable);

IsDisStable(jSpecific);

#We know from calculus that the equilibrium solution for Fsys [1.0,1.0,1.0] is a STABLE equilibrium solution as shown above.

Orbk(3, z, F, [1.0, 1.0, 1.0], 1, 100);

dataplot(%);

#The above verifies that [1.0,1.0,1.0] is an equilibrium solution. Lets check if it is STABLE numerically by starting the orbit at initial conditions very close to the equilibrium solution:[1.01,1.01,1.01]

Orbk(3, z, F, [1.01, 1.01, 1.01], 1, 1000);

dataplot(%);

{op(Orbk(3, z, F, [1.01, 1.01, 1.01], 1000, 1100))};

{op(Orbk(3, z, F, [1.01, 1.01, 1.01], 10000, 10100))};

#Even starting at initial conditions very close to the equilibrium solution, the third-order difference equation still converges to a period six solution that appears to never get any closer to the equilibrium solution [1.0,1.0,1.0] than it does after the first few iterations, but it is a stable oscillation as described by Keshet. Additionally, the arithmetic mean of the six values in the stable oscillation approach 1.0 with initial conditions close to 1.0

sumOfSix := ((((0.9872939141 + 0.9914995178) + 0.9957583404) + 1.004259728) + 1.008573360) + 1.012869608;

sumOfSix/6;

NULL;

#Now lets move on to a numeric based verification of the conjecture:

checkPeriodSixOrbit := proc(k, F) local ListOfLists, r, list, i; ListOfLists := []; r := evalf(rand(1.0 .. 10.0)); for i to k do ListOfLists := [op(ListOfLists), Orbk(3, z, F, [r(), r(), r()], 1000, 1005)]; end do; for list in ListOfLists do for i to 6 do if evalb(list[((i + 1) mod 6) + 1] - 0.0001 <= list[(i mod 6) + 1]/list[i] and list[(i mod 6) + 1]/list[i] <= list[((i + 1) mod 6) + 1] + 0.0001) = false then print(false); return ; end if; end do; end do; print(true); return ; end proc;

checkPeriodSixOrbit(1000, (1 + z[1])/(z[2] + z[3]));

#The above shows that for 100 orbits of F, each with different initial conditions (ranging from 1 to 10), the 1000th-1005th values of each orbit match the format of the period 6 solution proposed in conjecture four: that every third value is the ratio of the prior two.

NULL;

checkPeriodSixOrbit(100, (2 + z[1])/(z[2] + z[3]));

NULL;

checkPeriodSixOrbit(100, 0.5*(1 + z[1])/(z[2] + z[3]));

NULL;

checkPeriodSixOrbit(100, (0.5 + z[1])/(z[2] + z[3]));

NULL;

checkPeriodSixOrbit(100, (10 + z[1])/(z[2] + z[3]));

#The above shows that the resulting solution of period six: [phi, psi, psi/phi, 1/phi, 1/psi, phi/psi];appears to only hold true for the original difference equation.

#For the 3rd order difference equation: x(n)=0.5*(1+x(n-1))/(x(n-2)+x(n-3)), ;lets further explore!

Fsys2 := ToSys(3, z, 0.5*F);

{op(Orbk(3, z, 0.5*(1 + z[1])/(z[2] + z[3]), [0.1, 0.5, 0.8], 1000, 1012))};

#The above approaches a fixed point around x(n) = 0.640388;

Orbk(3, z, 0.5*(1 + z[1])/(z[2] + z[3]), [0.1, 0.5, 0.8], 1, 1000);

dataplot(%);

FP(Fsys2);

#The fixed point it approaches is x(n) = 0.6403882032;

SFP(Fsys2);

Orbk(3, z, 0.5*(1 + z[1])/(z[2] + z[3]), [0.6403882032, 0.6403882032, 0.6403882032], 1, 2000);

dataplot(%);

#The above shows the fixed point, but since it is not exact, maple is still plotting the minute oscillations. Note: Graphs don't reveal the entire truth all the time!

NULL;

Orbk(3, z, 0.5*(1 + z[1])/(z[2] + z[3]), [0.7403882032, 0.7403882032, 0.7403882032], 1, 2000);

dataplot(%);

#The above shows the plot with initial conditions varying ever so slightly from the fixed point x(n) = 0.6403882032;demonstrating that it is indeed stable!

NULL;

#Lets explore another variation: Changing the numerator parameter from 1.0 to 0.9

Orbk(3, z, (0.9 + z[1])/(z[2] + z[3]), [0.9, 0.1, 0.9], 1, 1000);

dataplot(%);

Fsys3 := ToSys(3, z, (0.9 + z[1])/(z[2] + z[3]));

FP(Fsys3);

SFP(Fsys3);

Orbk(3, z, (0.9 + z[1])/(z[2] + z[3]), [0.96, 0.96, 0.96], 1, 1000);

dataplot(%);

#The above shows that the fixed point x(n) = 0.9658910532;is not stable because starting the orbit with initial conditions very,very,very close to the fixed point [0.96,0.96,0.96] and the system diverges.

Orbk(3, z, (0.9 + z[1])/(z[2] + z[3]), [0.9658910532, 0.9658910532, 0.9658910532], 1, 1000);

dataplot(%);

#The above shows that starting the orbit with initial conditions as close as possible (with 10 digits of accuracy) to the fixed point and it stabilizes at the fixed point forever on...

NULL;

#Some other interesting plots below:

Orbk(3, z, (0 + z[1])/(z[2] + z[3]), [1.0, 1.0, 1.0], 1, 500);

dataplot(%);

jsys := ToSys(3, z, (0 + z[1])/(z[2] + z[3]));

jfp := FP(jsys);

j := JAC(jsys);

m := [];

for i in jfp do

print(`for fixed point`);

print(i);

print(`the associated jacobian matrix is`);

jm := subs({z[1] = i[1], z[2] = i[2], z[3] = i[3]}, j);

print(Matrix(jm));

print(`stable = true unstable = false`);

print(IsDisStable(jm));

print(`#########################`);

end do;

Orbk(3, z, (0.9 + z[1])/(z[2] + z[3]), [500.01, 10.0, 1.3], 1, 1000);

dataplot(%);

NULL;

Orbk(3, z, (1.0 + z[1])/(z[2]^1.01 + z[3]), [1.01, 1.5, 1.3], 1, 1000);

dataplot(%);

NULL;

Orbk(3, z, (1.0 + z[1])/(0.03 + z[2]^1.01 + z[3]), [1.01, 1.5, 1.3], 1, 1000);

dataplot(%);

NULL;

Orbk(3, z, (0.9 + z[1])/(0.03 + z[2]^1.1 + z[3]), [1.01, 10.5, 1.3], 1, 1000);

dataplot(%);

NULL;

Orbk(3, z, (0.9 + z[1])/(0.03 + z[2]^1.01 + z[3]), [1.01, 1.5, 1.2], 1, 1000);

dataplot(%);

{op(Orbk(3, z, (0.9 + z[1])/(0.03 + z[2]^1.01 + z[3]), [1.01, 1.5, 1.2], 10000, 10020))};

nops(%);