5
$\begingroup$

My last question on varying plot styles was perhaps too minimal of an example (and a duplicate to boot!). Here's a more realistic example that still gives me trouble.

I'd like to plot an isocline where one function (dg[x1,x2]) equals zero and style it differently depending on the sign of another function (dg2[x1,x2]). For example,

dg[x1_, x2_] := (2 (x1 + E^(4 (x1 - x2)^2) x1 - 2 x1^3 - 2 x2 + 2 x1^2 x2 + 2 E^(2 (x1 - x2)^2) (x1 - x2) (-1 + x2^2)))/((-1 + E^(4 (x1 - x2)^2)) (-1 + x1^2));
dg2[x1_, x2_] := 2/((-1 + E^(4 (x1 - x2)^2)) (-1 + x1^2)^2) (-1 + 7 x1^2 - 10 x1^4 + 8 x1^6 + E^(4 (x1 - x2)^2) (1 - 7 x1^2 + 2 x1^4) - 8 x1 x2 + 24 x1^3 x2 - 16 x1^5 x2 + 8 x2^2 - 16 x1^2 x2^2 + 8 x1^4 x2^2 - 8 E^(2 (x1 - x2)^2) (x1 - x2) (-1 + x2) (1 + x2) (x1^3 + x2 - x1^2 x2));

Here's what the isocline looks like and the regions that should be styled differently:

iso = ContourPlot[dg[x1, x2] == 0, {x1, -1, 1}, {x2, -1, 1},
  Exclusions -> {x1 == x2}, PlotPoints -> 50];
cp = ContourPlot[dg2[x1, x2], {x1, -1, 1}, {x2, -1, 1}, 
  PlotPoints -> 50, Contours -> {0}, ContourShading -> {White, LightRed}];
Show[cp, iso]

Mathematica graphics

My approach so far is to extract the contours, then ListLinePlot and color them using ColorFunction:

ExtractPlotPoints[plot_Graphics] := Cases[Normal@plot, Line[x_] :> x, \[Infinity]];

isocol = ListLinePlot[ExtractPlotPoints[iso],
  ColorFunction -> (If[dg2[#1, #2] > 0, Red, Black] &),
  ColorFunctionScaling -> False, Frame -> True, Axes -> False, 
  AspectRatio -> 1, PlotRange -> {{-1, 1}, {-1, 1}}];
Show[cp, isocol]

Mathematica graphics

This works great, but I'd like other style options besides varying the color (say, thick vs thin). This is where I'm stuck. I tried MeshFunctions but got nothing:

ListLinePlot[ExtractPlotPoints[iso], 
  MeshFunctions -> {dg2[#1, #2] &}, Mesh -> {{0}}, 
  MeshShading -> {Directive[Black, Thick], Directive[Red, Thin]},
  MeshStyle -> None, Frame -> True, Axes -> False, AspectRatio -> 1]

Mathematica graphics

dg2 is sketchy when x1==x2 due to the denominator equalling zero, but it has a well-defined limit.

Limit[dg2[x1, x2], {x2 -> x1}]
(* (-2 + 4 x1^2)/(-1 + x1^2) *)

Defining a better dg2 when x1==x2 helps, but isn't quite right:

dg2[x1_, x1_] = (-2 + 4 x1^2)/(-1 + x1^2);
isonew = ListLinePlot[ExtractPlotPoints[iso], 
  MeshFunctions -> {dg2[#1, #2] &}, Mesh -> {{0}}, 
  MeshShading -> {Directive[Black, Thick], Directive[Red, Thin]}, 
  MeshStyle -> None, Frame -> True, Axes -> False, AspectRatio -> 1];
Show[cp, isonew]

Mathematica graphics

Any thoughts? The functions dg and dg2 may differ from this example.

(The eventual application is isoclines of evolutionary dynamics as in Fig. 6a of this paper.)

$\endgroup$
7
$\begingroup$
cp2 = Quiet@ ContourPlot[{ConditionalExpression[dg[x1, x2], dg2[x1, x2] < 0] ==  0,
 ConditionalExpression[dg[x1, x2], dg2[x1, x2] > 0] ==  0}, 
{x1, -1, 1}, {x2, -1, 1}, PlotPoints -> 50,  Exclusions -> {x1 == x2}, 
ContourStyle -> {Directive[Black, Thick, Dashed], Directive[Red, Thin]}];

Show[cp, cp2]

enter image description here

Alternatively,

dg2b[x1_, x2_] :=  If[x1 == x2, (-2 + 4 x1^2)/(-1 + x1^2), 
   2/((-1 + E^(4 (x1 - x2)^2)) (-1 + x1^2)^2) (-1 + 7 x1^2 - 
   10 x1^4 + 8 x1^6 + E^(4 (x1 - x2)^2) (1 - 7 x1^2 + 2 x1^4) - 
   8 x1 x2 + 24 x1^3 x2 - 16 x1^5 x2 + 8 x2^2 - 16 x1^2 x2^2 + 
   8 x1^4 x2^2 - 8 E^(2 (x1 - x2)^2) (x1 - x2) (-1 + x2) (1 + x2) (x1^3 + x2 - x1^2 x2))];

isonewb = ListLinePlot[ExtractPlotPoints[iso], 
   MeshFunctions -> {dg2b[#1, #2] &}, Mesh -> {{0}}, 
   MeshShading -> {Directive[Black, Thick], Directive[Red, Thin]}, 
   MeshStyle -> None, Frame -> True, Axes -> False, AspectRatio -> 1];
Show[cp, isonewb]

enter image description here

$\endgroup$
  • $\begingroup$ Thanks, that's elegant. Unfortunately for me, it's significantly slower (4.33s vs 1.2s). Of course I can spare 3.13s in this example, but sometimes dg and dg2 are a lot more computationally expensive, so I'd like to save that 3.6X speedup if possible. But maybe it's an inevitable trade-off. $\endgroup$ – Chris K Feb 21 '18 at 4:28
  • $\begingroup$ Huh, when I run your edited version, I get the same as my last figure where the parts don't line up (M11.2 on MacOS10.13.3). ¯\_(ツ)_/¯ $\endgroup$ – Chris K Feb 22 '18 at 3:32
  • $\begingroup$ @ChrisK, version 9.0 Windows 10 here. $\endgroup$ – kglr Feb 22 '18 at 3:54
0
$\begingroup$

My colleague suggested the following variation, which is both fast and accurate (avoids potential problems I had with MeshFunctions). Idea: make ContourPlot, extract contours, make a separate ListLinePlot for each desired style using Opacity[0] in ColorFunction to blank out the undesired parts, then combine.

iso1 = ListLinePlot[ExtractPlotPoints[iso], ColorFunctionScaling -> False,
  ColorFunction -> (If[dg2[#1, #2] > 0, {Red}, {Opacity[0]}] &),
  PlotStyle -> Thin];
iso2 = ListLinePlot[ExtractPlotPoints[iso], ColorFunctionScaling -> False,
  ColorFunction -> (If[dg2[#1, #2] < 0, {Black}, {Opacity[0]}] &),
  PlotStyle -> Thick];
Show[cp, iso1, iso2, Frame -> True, Axes -> False, AspectRatio -> 1,
  PlotRange -> {{-1, 1}, {-1, 1}}]

Mathematica graphics

$\endgroup$

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.