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




(* Code for producing a bifurcation diagram for problem 3 *)

f[x_] := L x - x^3
fixedpointdiagram = ContourPlot[f[x] - x == 0, {L, -2, 3}, {x, -2, 2}]

(* This code will produce the bifurcation diagram for the fixed points of the family f[x] = L x - x^3 *)
(* The function ContourPlot[rel , {x,xmin,xmax}, {y,ymin,ymax}] will plot all points (x,y) satisfying the relation rel inside the region xmin <= x <= xmax, ymin <= y <= ymax *)

g[x_] := Simplify[((f[f[x]] - x)/(f[x] - x))]
twocyclediagram = ContourPlot[g[x] == 0, {L, -2, 3}, {x, -2, 2}, ContourStyle -> RGBColor[{0.882, 0.612, 0.141}]]

(* This code will produce the bifurcation diagram for the period-2 points of the family f[x] = L x - x^3 *)
(* The option ContourStyle -> RGBColor[{0.882, 0.612, 0.141}] is used to change the color of the contour *)

fulldiagram1 = Show[fixedpointdiagram, twocyclediagram]

(* The function Show[] will combine any number of graphics objects.  In our case, it will combine the two contour plots to make the full bifurcation diagram *)


fulldiagram2 = ContourPlot[{f[x] - x == 0, g[x] == 0}, {L, -2, 3}, {x, -2, 2}]

(* ContourPlot is also capable of plotting several contours at once: for example, ContourPlot[{rel1, rel2, rel3}, {x,xmin,xmax}, {y,ymin,ymax}] will plot the three relations rel1, rel2, rel3 *)

(* In this case, Mathematica automatically gives each contour a different color *)





(* Code for producing a bifurcation diagram for problem 2a *)

f[x_] := L Sin[x]

(* Defines the function of interest *)

fulldiagram2a = 
 ContourPlot[
 {f[f[x]] == x , f[x] == x},
  {L, -4, 2}, {x, -4, 4}, 
  PlotPoints -> 40, 
  ContourStyle -> { RGBColor[{0.882, 0.612, 0.141}], RGBColor[{0.369, 0.506, 0.710}]}
]

(* This plots the contours f[f[x]] = x and f[x] = x identifying period-2 and period-1 points.  Because it is not possible to give a simple algebraic expression for solving f[f[x]] = x that doesn't include the solutions of f[x] - x, we plot the period-2 point contour first, and then the fixed point contour afterwards on top of it.  Because of this, we change the colors of the contours using ContourStyle at the end to be consistent with the other examples. *)

(* The PlotPoints option increases the number of points used to compute the contours from the default setting 20, which smooths out the picture where the contours intersect.  Try removing this line and looking at the difference in the plot. *)




(* Here is code to compute Schwarzian derivatives *)

SchwarzianDerivative[f_][x_] := Derivative[3][f][x]/Derivative[1][f][x] - 3/2 (Derivative[2][f][x]/Derivative[1][f][x])^2
       (* This will define a new operator, SchwarzianDerivative.  To use it, apply it to a function and then apply that new "function" to an argument *)

g[x_] := Sin[2x]
SchwarzianDerivative[g][x]
       (* This should return the Schwarzian derivative of sin(2x), which is -4-6tan(2x)^2 *)




(* Here is code to produce the orbit diagram for problem 6 *)



Clear[f]
(* Removes any previous definition of f *)

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 *)





(* 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, a]

(* Clears the definition of f and a, 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 0.002 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 to make the diagram better *)
       (* 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 *)
       (* If the diagram is blurry, try adjusting PointSize[], and also changing the step size in the point list when computing the list of points *)



(* 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 *)