home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / e / e032 / 3.ddi / FILES / GRAPHICS.PAK / IMPLICIT.M < prev    next >
Encoding:
Text File  |  1992-07-29  |  8.0 KB  |  231 lines
  1.  
  2. (* :Title: ImplicitPlot *)
  3.  
  4. (* :Authors: Jerry B. Keiper, Wolfram Research, Inc.,
  5.         contour plot method: Theo Gray, Jerry Glynn, Dan Grayson *)
  6.  
  7. (* :Summary:  Plot[ ] requires a function.  Many simple graphs (e.g.,
  8.     circles, ellipses, etc.) are not functions.  ImplicitPlot[ ] allows
  9.     the user to easily sketch figures defined by equations. *)
  10.  
  11. (* :Context: Graphics`ImplicitPlot` *)
  12.  
  13. (* :Package Version: 2.1 *)
  14.  
  15. (* :History:
  16.     V2.0 by Jerry B. Keiper, April 1991
  17.     V2.1 Modifications by John M. Novak
  18. *)
  19.  
  20. (* :Keywords: solution set, graphics *)
  21.  
  22. (* :Sources: The contour plot alternate method is from:
  23.     Gray, Theodore and Glynn, Jerry, Exploring Mathematics
  24.         with Mathematica, (Addison-Wesley, 1991) *)
  25.  
  26. (* :Warning: *)
  27.  
  28. (* :Mathematica Version: 2.0 *)
  29.  
  30. (* :Limitation:
  31.     ImplicitPlot[ ] relies on Solve[ ] for much of the work.  If Solve[ ]
  32.     fails, no plot can be made.
  33.  
  34.     Subscripted variables (e.g., x[1], x[2]) cannot be used.
  35. *)
  36.  
  37. BeginPackage["Graphics`ImplicitPlot`","Utilities`FilterOptions`"]
  38.  
  39. ImplicitPlot::usage =
  40.     "ImplicitPlot[eqn, {x, a, b}] draws a graph of the set of points
  41.     that satisfy the equation eqn.  The variable x is associated with
  42.     the horizontal axis and ranges from a to b.  The remaining
  43.     variable in the equation is associated with the vertical axis.
  44.     ImplicitPlot[eqn, {x, a, x1, x2, ..., b}] allows the user to specify
  45.     values of x where special care must be exercised.
  46.     ImplicitPlot[{eqn1,eqn2,..},{x,a,b}] allows more than one equation
  47.     to be plotted, with PlotStyles set as in the Plot function.
  48.     ImplicitPlot[eqn,{x,a,b},{y,a,b}] uses a contour plot method of
  49.     generating the plot.  This form does not allow specification
  50.     of intermediate points..."
  51.  
  52. Options[ImplicitPlot] =
  53.     {AspectRatio -> Automatic, Axes -> Automatic, AxesLabel -> None,
  54.     AxesOrigin -> Automatic, AxesStyle -> Automatic,
  55.     Background -> Automatic, ColorOutput -> Automatic,
  56.     DefaultColor -> Automatic, Epilog -> {}, Frame -> False,
  57.     FrameLabel -> None, FrameStyle -> Automatic, FrameTicks -> Automatic,
  58.     GridLines -> None, PlotLabel -> None, PlotPoints -> 25,
  59.     PlotRange -> Automatic, PlotRegion -> Automatic,
  60.     PlotStyle -> Automatic, Prolog -> {}, RotateLabel -> True,
  61.     Ticks -> Automatic, DefaultFont :> $DefaultFont,
  62.     DisplayFunction :> $DisplayFunction}
  63.  
  64. Begin["`Private`"]
  65.  
  66. ImplicitPlot[eqns:{__Equal}, xr:{_,_,_},yr:{_,_,_},opts___] :=
  67.     Module[{ps,df},
  68.     {ps,df} = {PlotStyle,DisplayFunction}/.{opts}/.Options[ImplicitPlot];
  69.     ps = cyclestyles[ps,Length[eqns]];
  70.     gr = MapThread[ImplicitPlot[#1,xr,yr,
  71.         ContourStyle->#2,DisplayFunction->Identity,opts]&, {eqns,ps}];
  72.     gr = Select[gr,Head[#] === ContourGraphics &];
  73.     Show[gr,FilterOptions[ContourGraphics,opts,
  74.         Sequence @@ Options[ImplicitPlot]],DisplayFunction->df]/;
  75.             gr =!= {}
  76.     ]
  77.         
  78. ImplicitPlot[eqns:{__Equal}, {x_,a_,m___,b_}, opts___] :=
  79.     Module[{ps, df, gr, ln},
  80.     {ps,df} = {PlotStyle,DisplayFunction}/.{opts}/.Options[ImplicitPlot];
  81.     ps = cyclestyles[ps,Length[eqns]];
  82.     gr = MapThread[makegr[#1, {x,a,m,b},
  83.         PlotStyle->#2,DisplayFunction->Identity,opts]&, {eqns,ps}];
  84.     gr = Select[gr, (# =!= $Failed)&];
  85.     Show[Graphics[gr], FilterOptions[Graphics, opts,
  86.         Sequence @@ Options[ImplicitPlot]],DisplayFunction->df]/;
  87.             gr=!={}
  88.     ]
  89.  
  90. ImplicitPlot[lhs_ == rhs_, xr:{_,_,_},yr:{_,_,_},opts___] :=
  91.     With[{ps = PlotStyle/.{opts}/.Options[ImplicitPlot],
  92.             copts = FilterOptions[ContourPlot,opts,
  93.                 Sequence @@ Options[ImplicitPlot]]},
  94.     ContourPlot[lhs - rhs,xr,yr,
  95.         copts,
  96.         ContourStyle->ps,
  97.         Contours->{0},
  98.         ContourLines->True,
  99.         ContourShading->False,
  100.         ContourSmoothing->4]
  101.     ]
  102.  
  103. ImplicitPlot[eqn_Equal, {x_,a_,m___,b_}, opts___] :=
  104.     Module[{ps, df, gr},
  105.     {ps,df} = {PlotStyle,DisplayFunction}/.{opts}/.Options[ImplicitPlot];
  106.     gr = makegr[eqn, {x,a,m,b}, PlotStyle->ps, opts];
  107.     Show[Graphics[gr], FilterOptions[Graphics, opts,
  108.         Sequence @@ Options[ImplicitPlot]], DisplayFunction->df]/;
  109.             gr=!=$Failed
  110.     ]
  111.  
  112. cyclestyles[ps_,ln_] :=
  113.     Module[{style = ps},
  114.         If[Head[ps] =!= List,
  115.             style = {ps},
  116.             If[Length[ps] == 0, 
  117.                 style = {{}}]
  118.         ];
  119.         While[Length[style] < ln, style = Join[style,style]];
  120.         Take[style,ln]
  121.     ]
  122.  
  123. ImplicitPlot::var =
  124. "Equation `1` does not have a single variable other than `2`."
  125.  
  126. findy[f_, x_] :=
  127.     Module[{nf},
  128.     nf = Select[Union[Cases[f,
  129.             (_Symbol | _[(_?NumberQ)...]),
  130.                 Infinity]],
  131.         (!(NumberQ[N[#]] || #===x))&];
  132.     If[Length[nf] == 1,
  133.         nf[[1]],
  134.         (* else *)
  135.         Message[ImplicitPlot::var, eqn, x];
  136.         $Failed
  137.         ]
  138.     ]
  139.  
  140. ImplicitPlot::epfail = "Equation `1` could not be solved for points to plot."
  141.  
  142. makegr[eqn_Equal, {x_, a_, m___, b_}, opts___] :=
  143.     Module[{f = eqn[[1]] - eqn[[2]], ranges, plots, ar, y},
  144.     If[(y = findy[eqn, x]) === $Failed, Return[$Failed]];
  145.     ranges = Solve[f == 0 && D[f, y] == 0, {x, y}];
  146.     If[ListQ[ranges] && Length[ranges] > 0, ranges = N[x /. ranges]];
  147.     If[!VectorQ[ranges, NumberQ],
  148.         Message[ImplicitPlot::epfail, eqn];
  149.         Return[$Failed]];
  150.     ranges = Select[Chop[ranges], FreeQ[#, Complex]&];
  151.     ranges = Sort[Select[ranges, (a < # < b)&]];
  152.     ranges = Union[Sort[Join[ranges, N[{a, m, b}]]]];
  153.     ar = N[b-a]/10^8;
  154.     ranges = Transpose[{Drop[ranges+ar, -1], Drop[ranges-ar, 1]}];
  155.     (* ranges is now a (sorted) list of disjoint intervals with small
  156.         gaps between them where singularities probably exist. *)
  157.     plots = Map[rangeplot[f, x, y, #, opts]&, ranges]
  158.     ];
  159.  
  160. distx[{x_, y_List}] :=
  161.     Module[{yy = Sort[Select[Chop[y], FreeQ[#, Complex]&]]},
  162.         Transpose[{Table[x, {Length[yy]}], yy}]
  163.     ]
  164.  
  165. evenup[{}] = {};
  166.  
  167. evenup[xys_] :=
  168.     Module[{ll = Length /@ xys},
  169.         If[Max[ll] == Min[ll], Return[xys]];
  170.         If[Length[ll] < 3, {}, Drop[Drop[xys, 1], -1]]
  171.     ]
  172.  
  173. rangeplot[f_, x_, y_, {a_, b_}, opts___] :=
  174.     Module[{pp, ps, j, multipoints, mdpt, len},
  175.     {pp, ps} = {PlotPoints - 1, PlotStyle} /. {opts} /.
  176.                         Options[ImplicitPlot];
  177.     If[ps === Automatic,ps = {}];
  178.     mdpt = (a+b)/2;
  179.     len = (b-a)/2;
  180.     multipoints = N[{#, y /. Solve[f==0 /. x -> #, y]}& /@
  181.             Table[N[mdpt + len Cos[j Pi/pp]], {j, pp, 0, -1}]];
  182.     multipoints = evenup[Select[distx /@ multipoints, (Length[#] > 0)&]];
  183.     If[Length[Dimensions[multipoints]] =!= 3,
  184.         multipoints = {}];
  185.     (* connect the dots to form the various curves *)
  186.     If[Length[multipoints] > 0,
  187.         Flatten[{ps,Line[#]}]& /@ Transpose[multipoints, {2,1,3}],
  188.         (* else *)
  189.         {}]
  190.     ];
  191.  
  192. Protect[ImplicitPlot];
  193.  
  194. End[] (* "`Private`" *)
  195.  
  196. EndPackage[]  (* "Graphics`ImplicitPlot`" *)
  197.  
  198. (* :Tests: *)
  199. (* :Examples:
  200.  
  201. ImplicitPlot[x^2 + 2 y^2 == 3, {x, -2, 2}] (* ellipse *)
  202. ImplicitPlot[(x^2 + y^2)^2 == (x^2 - y^2), {x, -2, 2}] (*lemniscate *)
  203. ImplicitPlot[(x^2 + y^2)^2 == 2 x y, {x, -2, 2}] (* lemniscate *)
  204. ImplicitPlot[x^3 + y^3 == 3 x y, {x, -3, 3}] (* folium of Descarte *)
  205. ImplicitPlot[x^2 + y^2 == x y + 3, {x, -3, 3}] (* ellipse *)
  206. ImplicitPlot[x^2 + y^2 == 3 x y + 3, {x, -10, 10},
  207.     PlotRange -> {{-10,10},{-10,10}}] (* hyperpola *)
  208. ImplicitPlot[{(x^2 + y^2)^2 == (x^2 - y^2),
  209.     (x^2 + y^2)^2 == 2 x y}, {x,-2,2},
  210.     PlotStyle->{GrayLevel[0],Hue[0]}] (* combined plots *)
  211. ImplicitPlot[{(x^2 + y^2)^2 == (x^2 - y^2),
  212.     (x^2 + z^2)^2 == 2 x z}, {x,-2,2},
  213.     PlotStyle->{GrayLevel[0],Dashing[{.01}]}] (* combined plots *)
  214. ImplicitPlot[{a == b, x^2 + 2 y^2 == 3}, {x, -1, 1}] (* one bad plot *)
  215. ImplicitPlot[x^2 + y^2 == Pi, {x, -2, 2}] (* OK eqn with 3 symbols *)
  216. ImplicitPlot[Sin[x] == Cos[y], {x, 1.5, Pi/2, 1.7}]
  217. (* contour method *)
  218. ImplicitPlot[Sin[2 x] + Cos[3 y] == 1,{x,-2 Pi,2 Pi},{y,-2 Pi,2 Pi}]
  219. ImplicitPlot[x^2 + x y + y^2 == 1,{x,-2Pi,2Pi},{y,-2Pi,2Pi}]
  220. ImplicitPlot[x^3 + x y + y^2 == 1,{x,-2Pi,2Pi},{y,-2Pi,2Pi}]
  221. ImplicitPlot[x^3 - x^2 == y^2 - y,{x,-1,2},{y,-1,2}]
  222. (* failure cases *)
  223. ImplicitPlot[a == b, {x, -1, 1}] (* bad plot *)
  224. ImplicitPlot[x^y == y^x, {x, -1, 1}] (* bad plot *)
  225. ImplicitPlot[{a == b, c == d}, {x, -1, 1}] (* bad plots *)
  226. ImplicitPlot[x^2 + y^2 == z, {x, -2, 2}] (* bad eqn with 3 vars *)
  227. ImplicitPlot[Sin[x] == y, {x, -3, 3}] (* Solve fails... *)
  228. ImplicitPlot[Sin[x] == Cos[y], {x, -5, 5}]
  229. ImplicitPlot[x^y == y^x, {x, -3, 3}]
  230. *)
  231.