9
$\begingroup$

Suppose I have the following ContourPlot3D:

ContourPlot3D[
p1 p2^6 (C p1 - (-1 + C) p2)^3 ((-1 + C) p1 - C p2)^2 
== 
(-1 + p1) (-1 + p2)^6 (-1 + C p1 - (-1 + C) p2)^3 (1 + (-1 + C) p1 - C p2)^2,
{p1, 0, 0.3}, {p2, 0.5, 1}, {C, 0, 1}, 
 Lighting -> ({"Directional", White, #} & /@ Tuples[{-1, 1}, 3]), 
 Mesh -> None, ContourStyle -> {Directive[Red, Opacity[0.5]]}, 
 AxesLabel -> {"\!\(\*SubscriptBox[\(p\), \(1\)]\)", 
 "\!\(\*SubscriptBox[\(p\), \(2\)]\)", "C"}]

which looks like:

enter image description here

My question is: how to project this transparent red surface onto the $p_1,p_2$ plane and then convert the projection to a standard 2d plot (like what we will get by using Plot[])

Additional requirement: It would be nice if the resulting 2D plot is clean, because I may need to plot other curves on top of the 2D plot.

$\endgroup$
  • $\begingroup$ RegionPlot[]? $\endgroup$ – J. M. will be back soon Jun 1 '13 at 6:40
  • $\begingroup$ I think RegionPlot[] works only for flat plane. $\endgroup$ – wdg Jun 1 '13 at 6:50
9
$\begingroup$

First let's save the plot in cp:

cp = ContourPlot3D[
  p1 p2^6 (c0 p1 - (-1 + c0) p2)^3 ((-1 + c0) p1 - c0 p2)^2 ==
   (-1 + p1) (-1 + p2)^6 (-1 + p1 - (-1 + c0) p2)^3 (1 + (-1 + c0) p1 - c0 p2)^2,
  {p1, 0, 0.3}, {p2, 0.5, 1}, {c0, 0, 1},
  Lighting -> ({"Directional", White, #} & /@ Tuples[{-1, 1}, 3]), 
  Mesh -> None, ContourStyle -> {Directive[Red, Opacity[0.5]]}, 
  AxesLabel -> {"\!\(\*SubscriptBox[\(p\), \(1\)]\)", 
    "\!\(\*SubscriptBox[\(p\), \(2\)]\)", "C"}];

[By the way, you should avoid using capital C, which is used by Mathematica for constant parameters.]

The easy way

Project the coordinates of the polygons in cp onto the p1/p2 plane:

Graphics[
  Cases[cp, GraphicsComplex[pts_, g_, rest___] :>
   GraphicsComplex[
     Most /@ pts,  (* the first two coordinates are p1, p2 *)
     {Directive[Lighter[Red, 0.5]], Cases[g, _Polygon, Infinity]}
    ],
   Infinity],
  Method -> {"TransparentPolygonMesh" -> True}]

Envelope Envelope with axes

I replaced the semi-opacity with a semi-light red. This should be fine if you keep the graphics in the background. Otherwise, it's harder.

Graphics styling. To get axes on the graph and other features, add Axes -> True and other options:

Graphics[
  Cases[cp, GraphicsComplex[pts_, g_, rest___] :>
   GraphicsComplex[
     Most /@ pts,  (* the first two coordinates are p1, p2 *)
     {Directive[Lighter[Red, 0.5]], Cases[g, _Polygon, Infinity]}
    ],
   Infinity],
  Method -> {"TransparentPolygonMesh" -> True},
  Axes -> True]

(See above.) Others to consider are AxesOrigin -> {0, 0}, using Frame -> True instead of Axes, and setting the ratio height/width with AspectRatio -> 1/GoldenRatio.

Edit - The hard way -- Envelopes

First, call the equation eqn and change C to c0:

eqn = (p1 p2^6 (C p1 - (-1 + C) p2)^3 ((-1 + C) p1 - C p2)^2 ==
    (-1 + p1) (-1 + p2)^6 (-1 + C p1 - (-1 + C) p2)^3 (1 + (-1 + C) p1 - C p2)^2) /. C -> c0
p1 p2^6 (c0 p1 - (-1 + c0) p2)^3 ((-1 + c0) p1 - c0 p2)^2 ==
 (-1 + p1) (-1 + p2)^6 (-1 + c0 p1 - (-1 + c0) p2)^3 (1 + (-1 + c0) p1 - c0 p2)^2

The boundaries of the projection will be either the edges of the 3D plot or the envelope of the surface. The edges of the 3D graph are where c is 0 or 1. There are many solutions which we and we have to extract the ones corresponding to the plot. Luckily the only ones in the plot region are the boundaries and they can be written in the form p2 = func[p1]; hence, the boundaries can picked out by inspecting the value of p2 at some intermediate value of p1, say p1 = 0.15.

sols = Extract[#, Position[p2 /. # /. p1 -> 0.15, p_Real /; 0.5 <= p <= 1]] &@
  Solve[eqn && (0 == c0 || c0 == 1), {c0, p2}]
{{c0 -> 1, p2 -> (-1 + p1 + Sqrt[p1 - p1^2])/(-1 + 2 p1)},
 {c0 -> 0, p2 -> 1 - p1 ...omitted...}}

Here we have to set up the differential equation for the envelope, which will be the points on the surface where the partial derivative with respect to c0 is zero (envelope == 0); the condition "on the surface" is modeled by the integral curve being orthogonal to the gradient (surface == 0). One of the boundary values of the differential equation is {p1, p2, c0} == {0, 1, 0.} (I've change the value of c0 three times because I've been fooled by the Mathematica output, the numerics being so tricky). The DE is degenerate there, so we have to find a good initial point for NDSolve. As with sols there are many solutions and we have to extract the initial point we want. There is only one envelope line in the plot region, so finding it presents no trouble. The solutions found will be of the form p2 = fp2[p1], c0 = fc0[p1], where fp2, fc0 will be the InterpolatingFunction returned by NDSolve.

fnsrule = Thread[{p2, c0} -> {fp2[p1], fc0[p1]}];
envelope = D[Subtract @@ eqn, {c0}] /. fnsrule;
surface = (D[Subtract @@ eqn, {{p1, p2, c0}}] /. fnsrule) . {1, fp2'[p1], fc0'[p1]};
initsol = NSolve[D[Subtract @@ eqn, {c0}] == 0 && eqn == 0 /. p1 -> 1/100, {p2, c0}];
initpt = First @ Extract[initsol, 
   Position[{p2, c0} /. initsol, {p_Real, c_Real} /; 0.5 <= p <= 1 && 0 <= c <= 1]]
{p2 -> 0.751159, c0 -> 0.267272}

Expect warning messages. The choice of p1 -> 1/100 was through trial and error. (It's far enough away from p1 == 0 to stay out of serious trouble and give a good solution. Getting much closer tended to call for increasing WorkingPrecision, but increasing WorkingPrecision tended to cause the kernel to crash.) As you can see below, the DE can be solved for almost all the requested domain.

bdy = NDSolve[{envelope == 0, surface == 0,
        fp2[0.01] == p2 /. initpt, fc0[0.01] == c0 /. initpt},
       {fp2, fc0}, {p1, 0, 0.3}]
{{fp2 -> InterpolatingFunction[{{1.91025*10^-9, 0.3}}, <>],
  fc0 -> InterpolatingFunction[{{1.91025*10^-9, 0.3}}, <>]}}

However, the problem still has one more tricky pitfall lurking: The envelope ("green") appears to be the lower boundary of the projection and to be tangent at about p1 == 0.2 to the red line, which corresponds to the boundary of the ContourPlot where c0 == 0. In fact, the envelope passes out of the plot region into a region where c0 < 0, and the red line becomes the lower boundary, which can be seen by inspecting the value of c0 at p1 = 0.3.

boundariesPlot = Plot[Evaluate[Append[p2 /. sols, fp2[p1] /. First@bdy]], {p1, 0, 0.3}, 
  PlotStyle -> {Blue, Red, Darker@Green}]

Three potential boundary lines

fc0[0.3] /. First@bdy
-0.488191

So we need to find where the envelope crosses c0 == 0:

tanPt = p1 /. FindRoot[fc0[p1] /. bdy, {p1, 0.2}]
0.196357

Then we can plot the projection of the region: The first solution in sols (where c0 == 1) is the upper boundary. The lower is equal either fp2 or the second solution in sols depending on p1; you can use Piecewise or an If statement.

pl = With[{fp2 = fp2 /. First@bdy},
      Plot[
       {p2 /. First@sols,
        Piecewise[
          {{fp2[p1], fp2["Domain"][[1, 1]] <= p1 <= tanPt},
           {p2 /. Last@sols, True}}]
         },
       {p1, 0, 0.3}, Filling -> {1 -> {2}}, 
       FillingStyle -> Directive[Red, Opacity[0.5]], PlotStyle -> None, 
       PlotRange -> {0.5, 1}]]

Plot of projection

You can make a Polygon of the projection from a plot, if you omit PlotStyle -> None. It may be more convenient to combine it with other graphics. Extract the two graphs with Cases and stitch the points together, reversing the direction of one of the lines.

envPts = Join[First[#], Reverse@Last[#]] &@Cases[
    With[{fp2 = fp2 /. First@bdy},
     Plot[Evaluate[
       {p2 /. First@sols,
        Piecewise[
         {{fp2[p1], fp2["Domain"][[1, 1]] <= p1 <= tanPt},
          {p2 /. Last@sols, True}}]
        }],
      {p1, 0, 0.3}]], Line[pts_] :> pts, Infinity];
envPlot = Graphics[{Opacity[0.5], Red, Polygon[envPts]},
   Axes -> True, PlotRange -> {.5, 1}, AspectRatio -> 1/GoldenRatio];

Show[envPlot, boundariesPlot]

Comparison of projection and boundaries

Here we can compare the easy and hard ways. The easy way is recolored Blue and so differences will show up as Blue or light Red. Close inspection shows that the easy way has greater granularity typical of surface plots.

Show[proj /. _RGBColor -> Blue, envPlot,
 AspectRatio -> 1/GoldenRatio, PlotRange -> {0.5, 1}]
Show[proj /. _RGBColor -> Blue, envPlot,
 AspectRatio -> 1, PlotRange -> {{0, 0.02}, {.7, 1}}]

Comparison of projections Close up

While some steps might be automated, it seems that in any particular case, you're going to have to inspect what's going on and make an intelligent choice. The "easy way" above is fairly straightforward.

Edit 2 - 3rd way - MeshFunctions

We can use MeshFunctions to find the envelope. It can be a good way to solve for a curve on a surface.

This way still has its hangups. The graph gets extremely flat near p1 == 0 and it's hard to compute anything. Still, we can do well enough if we increase the precision and subdivide the plot region as we approach p1 == 0. The following goes down to p1 = 3*10^(-10) an order of magnitude at a time.

Table[
   mymesh[i] = With[{prec = 200},
     ContourPlot3D[Evaluate @ eqn,
      {p1, 3*10^(-i - 1), 3*10^(-i)}, {p2, SetPrecision[0.5, prec], 1}, {c0, 0, 1},
      PlotPoints -> 50, MaxRecursion -> 3, Mesh -> {{0}}, 
      BoundaryStyle -> None,  ContourStyle -> None,
      MeshFunctions -> {Function[{p1, p2, c0},
         Evaluate[CubeRoot@D[Subtract @@ eqn, {c0}]]]}, 
      WorkingPrecision -> prec]
     ], {i, 9}]; // Timing
{4.208199, Null}

The curve is still a ways from reaching p2 == 1, but it doesn't matter for the sake of plotting. We are concerned with p2 as a function of p1, and we don't need p1 much below 10^-4 or 10^-5 -- we can just connect a straight line to the endpoint. Here the segments of the envelope are colored differently:

With[{i1 = 1, i2 = 9}, 
 Graphics3D[
  MapThread[Cases[#1, 
     GraphicsComplex[pts_, g_, rest___] :> GraphicsComplex[pts, {#2, Thick,
        Select[Cases[g, _Line, Infinity], 
         Last @ pts[[ #[[1, 1]] ]] < 0.5 &&
           (First@pts[[ #[[1, 1]] ]] > 0.1 || Last@pts[[ #[[1, 1]] ]] > 0.001) &]
        }]] &,
   {Table[mymesh[i], {i, i1, i2}], Table[Hue[i/9], {i, i1, i2}]}],
  PlotRange -> {{0, 0.2}, {0.6, 1}, {0, 0.3}}, BoxRatios -> {1, 1, 1}, 
  Axes -> True, 
  AxesLabel -> {"\!\(\*SubscriptBox[\(p\), \(1\)]\)",
                "\!\(\*SubscriptBox[\(p\), \(2\)]\)", "C"}]
 ]

Envelope on surface

Collect the points down to p1 == 3 * 10^-5:

envelopePoints = Append[With[{i1 = 1, i2 = 4},
    Flatten[
     Map[Cases[#1, GraphicsComplex[pts_, g_, rest___] :>
         pts[[
           Select[Cases[g, _Line, Infinity], 
             Last@pts[[#[[1, 1]]]] < 0.5 &&
               (First@pts[[#[[1, 1]]]] > 0.1 ||
                  Last@pts[[#[[1, 1]]]] > 0.001) &][[1, 1]]
           ]] ] &,
      Table[mymesh[i], {i, i1, i2}]],
     2]
    ], {0., 1., 0.}];

Then join them with the other parts of the boundary and compare with the other methods.

envMeshPlot = 
 Graphics[{Opacity[0.5], Red, 
   Polygon @ Join[Table[{p1, p2 /. Last@sols}, {p1, 0.3, 0.2, -0.01}], 
     Most /@ envelopePoints, 
     Table[{p1, p2 /. First@sols} /. p1 -> p^2, {p, 0, Sqrt[0.3], Sqrt[0.3]/100}]]}]

Show[proj /. _RGBColor -> Blue, envMeshPlot,
 AspectRatio -> 1, PlotRange -> {{0, 0.005}, {.7, 1}}, Axes -> False]
Show[envMeshPlot, envPlot /. _RGBColor -> Blue,
 AspectRatio -> 1, PlotRange -> {{0, 0.005}, {.7, 1}}, Axes -> False]

Envelope and comparisions

The first method is easiest; the drawback is that the overlapping polygons make Opacity give an nonuniform appearance.

$\endgroup$
  • $\begingroup$ your easy way already solve my question, but I would like to see roughly what the hard way is if you could show me. $\endgroup$ – wdg Jun 3 '13 at 10:12
  • $\begingroup$ by the way, how can I add the axes and the axes labels to your 2d plot? $\endgroup$ – wdg Jun 3 '13 at 10:18
  • $\begingroup$ @wdg I responded to both comments in the updated answer. $\endgroup$ – Michael E2 Jun 3 '13 at 17:29
  • $\begingroup$ very grateful for your in-depth answers! $\endgroup$ – wdg Jun 5 '13 at 11:41
9
$\begingroup$

Add ViewPoint -> {0, 0, \[Infinity]} to your ContourPlot3D command and to get the projection onto the $p_1$ $p_2$ plane

enter image description here

$\endgroup$
  • $\begingroup$ It would be nice if you can clean up the region, because I need to plot other curves on top of that region. $\endgroup$ – wdg Jun 1 '13 at 7:39
  • $\begingroup$ What do you mean by "clean up"? If you change the options like Lighting and etc you can change the appearance. Or the Opacity you have specified makes it have the different shades. $\endgroup$ – bill s Jun 1 '13 at 7:46
  • $\begingroup$ How about Lighting->None and Opacity[0.1]? $\endgroup$ – bill s Jun 1 '13 at 8:02
  • $\begingroup$ The color of the region is still not uniform... $\endgroup$ – wdg Jun 1 '13 at 11:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.