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


(* Problem 1 ~ Code for exploring iterates of a quadratic function *)

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

qc[x_] := x^2 - 1.76
	(* The function of interest *)

NestList[qc, 0.0`, 700][[400 ;; 405]]
	(* This computes the 400th through 405th values in the critical orbit *)

qc[{0.0236, 0.0241}]
	(* Evaluates qc on the list of values given, which in this case are the endpoints of the first interval I1 given in part (b).  By default Mathematica will attempt to thread functions over a list, but be careful doing that in general. *)

Max[Abs[qc'[{0.0236, 0.0241}]]]
	(* Evaluates the derivative qc' on the list of values given and then computes the maximum of the absolute values of those quantities. *)



(* Problem 2 ~ Code for studying the quadratic family orbit diagrams *)

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

f[a_][x_] := x^2 + a;
	(* This is the function of interest.  Note the syntax for separating the parameter a from the function variable x *)

AsymptoticOrbit[a_] := Drop[NestList[f[a], 0.0`, 700], 200];
        (* This computes the 200th through 700th points in the orbit of the critical point 0 *)

OrbitPoints[a_] := {a, #} & /@ AsymptoticOrbit[a];
        (* This creates the list of points to plot at horizontal coordinate a *)

orbPicquad[lendpt_,rendpt_,numpts_] := 
 Graphics[{PointSize[0.002], Opacity[0.1], 
   Point[Flatten[Table[OrbitPoints[a], {a, lendpt, rendpt, (rendpt-lendpt)/numpts}], 
     1]]}, AspectRatio -> 1/2, Axes -> True, 
  AxesOrigin -> {Automatic, -2.5}, ImageSize -> 1000, TicksStyle -> 18]
	(* This draws the orbit diagram for the family where the parameter ranges from lendpt to rendpt, with a total number of numpts values drawn in the middle.  The larger numpts is, the more accurate the picture will be, but the longer it will take to render, as Mathematica has to compute 700*numpts values and then draw 500*numpts points.  Suggested value is 500 or 1000. *)

quadvalues[a_]:= NestList[f[a], 0.0`, 700][[500 ;; 520]]
	(* This computes the 500th through 520th values in the critical orbit *)


orbPicquad[-1.5,-1.2,500]
	(* This plots the orbit diagram on the range [-1.5,-1.2] with 500 intermediate points.  Mathematica on my computer reports taking 0.64 seconds to compute the points and then 10.82 seconds to draw the picture. *)







(* Problem 3 ~ Code to produce a general orbit diagram *)

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

f[a_][x_] := a Cos[x];
       (* This is the function of interest.  Note the syntax for separating the parameter a from the function variable x *)

AsymptoticOrbit[a_] := Drop[NestList[f[a], 0.0, 500], 200];
       (* This computes the 200th through 500th points in the orbit of the critical point 0 *)

OrbitPoints[a_] := {a, #} & /@ AsymptoticOrbit[a];
       (* This creates the list of points to plot at horizontal coordinate a *)

orbPic = Graphics[{PointSize[0.002], Opacity[0.1], Point[Flatten[Table[OrbitPoints[a], {a, 0, 8, 0.002}], 1]]}, 
    AspectRatio -> 1/2, Axes -> True, AxesOrigin -> {0, -8}]

       (* This will plot the points in each column in the range 0 <= a <= 8 with step size 0.002 and point size 0.002 *)
       (* Several other plot options are included to draw the coordinate axes and align them at the bottom of the plot*)

(* This code will take a nontrivial amount of time to run, on the order of 5-60 seconds, because Mathematica is computing 2 million points and then plotting 60% of them to make the picture *)



FindRoot[x Tan[x]==-1, {x,-3}]
	(* This will search numerically for a solution to x tan(x) = -1 with starting value x=-3.  Mathematica uses an efficient implementation of Newton's method to search for roots. *)






(* Here is code that can produce an orbit diagram when there are multiple critical points *)
(* This is not needed for the homework assignment, but is provided in case you want to explore some of the examples yourself *)
(* See here for additional discussion: http://mathematica.stackexchange.com/questions/13277/mathematica-code-for-bifurcation-diagram *)




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

f[a_][x_] := a*x - 1.0 x^3;
cps[a_] = x /. Quiet[Solve[D[f[a][x], x] == 0, x], Solve::ratnz]
       (* This generates a list of critical points of f, approximated *)
       (* Solve will require f to be a polynomial map for this to work properly *)

criticalOrbits[a_, cp_] := 
  Module[{try}, 
   If[Head[cp] === Real, 
    try = NestWhileList[f[a], cp, Abs[#] < 100 &, 1, 500];
    If[Abs[Last[try]] >= 100, try = {}, 
     try = Drop[{a, #} & /@ try, 100]], {}]];
       (* This computes the orbit of each critical point, with some graceful handling of the cases where the critical orbit goes to infinity or becomes non-real *)

points[k_] := 
 Partition[
  Flatten[Table[criticalOrbits[a, cps[a][[k]]], {a, -1, 4, 0.002}]], 2]
       (* This computes the list of points to plot at horizontal coordinate a *)
       (* Suggestion: increase the step size to make the code run faster in order to get a rough estimate of an appropriate range for a *)
       (* Then re-run the code with a smaller step size.  You may also want to adjust Pointsize[] in the code below *)
       (* In this case, the range [-1,4] includes points that will not give a bounded or real critical orbit *)

orbdiag = 
 Graphics[{Opacity[0.02], PointSize[0.002], 
   Table[{ColorData[1, k], Point[points[k]]}, {k, 1, 
     Length[cps[a]]}]}, Frame -> False, AspectRatio -> 1/2, 
  Axes -> True, AxesOrigin -> {0, -2}]

       (* This will plot all of the points to create the diagram, with each critical point shown in a different color *)
       (* In this case there are two critical orbits, shown in red and blue *)



(* Here is code that can export graphics to a file *)


Directory[]
       (* This will return the location of Mathematica's current output directory. *)
       (* If you wish to change the directory, use SetDirectory["C:/yourdirectorypathname"] *)

Export["orbit_diagram_1.png", orbPic, ImageResolution -> 100]
       (* This will export your orbit diagram in png format to the indicated filename in the output directory *)
       (* Mathematica supports many export filetypes.  If you are using a TeX compiler, the eps graphic format is well supported.  I typically use png. *)
       (* Different graphics formats have varying quality.  Consult wikipedia or the CS department for more information on graphics formats. *)

Export["orbit_diagram_1.png", orbPic, ImageResolution -> 400, ImageSize -> 900]
       (* Higher resolutions and larger output sizes will take correspondingly longer. This computation may take 60 seconds or longer. *)
       (* ImageResolution is in dpi, while ImageSize represents width in pixels *)