(* You should be able to copy and paste this text file directly into Mathematica.  These notes will be interpreted as commented code. *)


(* Problem 4ab ~ Code for exploring itineraries for a quadratic map *)


Clear[qc, x, a]
	(* Remove previous definitions, for safety *)


qcexact[x_] := x^2 - 15/4
	(* The quadratic function of interest, as an exact function *)

intI := -pplus <= x <= pplus
intI0 := -pplus <= x <= -Sqrt[-pplus + 15/4]
intI1 := Sqrt[-pplus + 15/4] <= x <= pplus
	(* Computes the intervals I, I0, and I1.  For numerical approximations of these values use N[] around them. *)

IntervalLocation[xval_] := 
 Which[intI0 /. x -> xval, 0, intI1 /. x -> xval, 1, True, "Unknown"]
	(* Helper function that determines which interval I0 or I1 a given input value x lies in.  If x lies in neither interval the function returns "Unknown" *)

ComputedItinerary[xval_] := 
 Map[IntervalLocation, NestList[qcexact, xval, 10]]
	(* Computes the 0th through 10th digits of the itinerary of a given x value under the quadratic map.  Note that it may be necessary to specify a high precision of the input number xval in order to get accurate results.  If more digits are desired, simply increase the number 10. *)

x1 := -1.50000`30
x2 := -2.2320508`30
x3 := 1.58367`30
	(* The three values from part 1b, specified to 30-digit precision *)

ComputedItinerary[x1]
	(* This will compute the first ten digits of the itinerary of the point x1 *)





(* Problem 4cd ~ Code for computing itinerary intervals for a quadratic map *)


qc[x_] := N[x^2 - 15/4]
	(* The quadratic function of interest, with approximated arithmetic because we will be computing lengthy inverse iterates and using exact arithmetic will quickly become unreasonable. *)

pplus := Max[x /. Solve[qc[x] == x, x]]
	(* Computes the larger fixed point *)

intI := -pplus <= x <= pplus
	(* Computes the starting interval *)

invimage0[interval_] := 
  Quiet[Reduce[ReplaceAll[interval, x -> qc[x]] && x < 0, x]];
invimage1[interval_] := 
  Quiet[Reduce[ReplaceAll[interval, x -> qc[x]] && x > 0, x]];
	(* Helper functions that compute the inverse image of an interval on the I0-side and the I1-side *)

IntervalFromItinerary[interval_, itinerarylist_] := 
 Which[
	itinerarylist == {}, 
		interval, 		(* if list is empty, just return the interval *)
	itinerarylist[[-1]] == 0,   
		IntervalFromItinerary[invimage0[interval], Drop[itinerarylist, -1]], (* if list ends in a 0, compute the inverse image on the 0 side *)
   	itinerarylist[[-1]] == 1, 
  		IntervalFromItinerary[invimage1[interval], Drop[itinerarylist, -1]] (* if list ends in a 1, compute the inverse image on the 1 side *)
]
(* Computes the successive inverse images of a starting interval so that the inverse images follow the list of 0s and 1s in the itinerarylist *)
(* Note that because we are computing this interval using successive inverse images, we use the term at the end of the itinerary list rather than the beginning to decide where to go next. *)

ItineraryInterval[itinerarylist_] := 
 IntervalFromItinerary[intI, itinerarylist]
	(* Computes the interval of points whose label is itinerarylist *)

ItineraryInterval[{0}]
	(* This is the interval whose itinerary is (0): the interval I_0 *)

ItineraryInterval[{0,0}]
	(* This is the interval whose itinerary is (00): the interval I_00 *)




(* Problem 4efg ~ Code for Newton's method on iterates *)


qc[x_] := N[x^2 - 15/4]
	(* The quadratic function of interest, with approximated arithmetic because we will be computing lengthy iterates and using exact arithmetic will quickly become unreasonable. *)

f[x_] := Nest[qc, x, 2] - x
	(* The second iterate of f, minus x, to search for 2-cycles *)
        (* To search for cycles of a different length, simply replace 2 with the desired length *)

nif[x_] := x - f[x] / f'[x]
	(* The Newton iterating function of f *)

NestList[nif, 2.49, 10]
	(* Iterates the Newton function to search for a point of order dividing 2 *)
        (* In this case the starting value 2.49 is very close to a fixed point of qc, so it will find the fixed point x=2.5 *)
	(* Replace the starting value with something else as suggested in the question to find the 2-cycle *)




(* Problem 6 ~ Code for trying to compute cycles directly *)


q[x_] := x^2 - 6
r1[x_] := Nest[q, x, 8]
list1 := NSolve[r1[x] == x, x]

(* This asks Mathematica to numerically compute the 256 points of order dividing 8 for q *)

r1[x] - x /. list1

(* This evaluates r1[x] - x where x is replaced, successively, by the values in list1 *)
(* If the computations are accurate, the results should be close to zero *)


r2[x_] := (Nest[q, x, 8] - x)/(Nest[q, x, 4] - x)
list2 = NSolve[r2[x] == 0, x]


(* This asks Mathematica to numerically compute the 256 points of order exactly 8 for q *)
(* This computation may take some time! *)

r1[x] - x /. list2

(* This evaluates r1[x] - x where x is replaced, successively, by the values in list2 *)
(* If the computations are accurate, the results should be close to zero *)