# “Beautifying” a ContourPlot

The following code to draw a ContourPlot creates a thick line (in some places only) along the diagonal. Since in some places the line is thinner, I wonder if there is anything I can do to reduce the thickness throughout?

ContourPlot[((256 ((1 - p1) p1)^(7/
2) Sqrt[((1 - p2) p2)/((1 - p1) p1)])/(4 (1 - p1) p1 -
4 (1 - p2) p2)^2 - (512 (1 - p1) p1 ((1 - p2) p2)^(5/2))/(4 (1 - p1) p1 +
4 (1 - p2) p2)^2 + (256 ((1 - p1) p1)^(5/2) (1 - p2) p2 (-2 +
Sqrt[((1 - p2) p2)/((1 - p1) p1)]))/(4 (1 - p1) p1 - 4 (1 - p2) p2)^2 +
16 (1 - p1)^2 p1^2 (3 Sqrt[
1/((1 - p2) p2)] - (32 ((1 - p2) p2)^(3/2))/(4 (1 - p1) p1 +
4 (1 - p2) p2)^2) - (16 ((1 - p1) p1)^(3/
2) (4 (1 - p1) p1 Sqrt[((1 - p2) p2)/((1 - p1) p1)] +
4 (1 - p2) p2 (-3 + 2 Sqrt[((1 - p2) p2)/((1 - p1) p1)])))/(4 (1 -
p1) p1 - 4 (1 - p2) p2) + (2 (-384 (1 - p1)^3 p1^3 +
768 (1 - p1)^2 p1^2 (1 - p2) p2 +
128 (1 - p2)^3 p2^3 Sqrt[((1 - p2) p2)/((1 - p1) p1)] -
128 (1 - p1) p1 (1 - p2)^2 p2^2 (1 +
3 Sqrt[((1 - p2) p2)/((1 - p1) p1)])))/(Sqrt[
1/((1 - p1) p1)] Sqrt[((1 - p2) p2)/((1 - p1) p1)] (4 (1 - p1) p1 -
4 (1 - p2) p2)^2))/(8 Sqrt[(1 - p2) p2]), {p1, 0.5, 1}, {p2, 0.5, 1},
PlotPoints -> 100, MaxRecursion -> 2, PlotLegends -> Automatic,
FrameLabel -> {"\!$$\*SubscriptBox[\(p$$, $$1$$]\)",
"\!$$\*SubscriptBox[\(p$$, $$2$$]\)"}]


• The 'line' in the plot is caused by the function having the value of 1/0 when p1 == p2. Why would you want to hide that? – Rohit Namjoshi Jan 10 at 21:28
• I don't want to hide this fact. But the way the ContourPlot shows this fact is crude in my view. Its not homogeneous along the 45 degrees. Like I said, the line is thick some places and thin at others. – user120911 Jan 10 at 21:40
• Try removing the PlotPoints and MaxRecursion options. – Rohit Namjoshi Jan 10 at 21:43
• That makes things much worse. – user120911 Jan 10 at 21:59
• PlotPoints -> 200, MaxRecursion -> 0 gives a fairly consistent thickness for me. – egwene sedai Jan 10 at 22:18

Let fun be the OP's function. If you simplify the expression, the removable discontinuity along p1 == p2 cancels out. Here's one way to help with the algebra: rewrite the function in terms of u1^2 == p1(1 - p1) and u2^2 == p2(1 - p2). Then the function becomes a rational function, which is much easier to simplify. Back substitution brings the simplified expression back to an equivalent function of p1 and p2.

Simplify[fun /. {
1 - p1 -> 1 - p1,     (* trick to replace just p1 in p1(1-p1) *)
p1 -> u1^2/(1 - p1),
1 - p2 -> 1 - p2,     (* trick to replace just p2 in p2(1-p2) *)
p2 -> u2^2/(1 - p2)},
u1 > 0 && u2 > 0]
fun2 = % /. {u1 -> Sqrt[p1 (1 - p1)], u2 -> Sqrt[p2 (1 - p2)]};
(*
(2 u1 u2 (u1^4 + 2 u1^3 u2 - 2 u1^2 u2^2 + 2 u1 u2^3 + u2^4)) /
((u1 + u2)^2 (u1^2 + u2^2))
*)

ContourPlot[fun2, {p1, 0.5, 1}, {p2, 0.5, 1},
FrameLabel -> {HoldForm@Subscript[p, 1], HoldForm@Subscript[p, 2]}]


As you've written the function, I agree with @RohitNamjoshi : "Why would you want to hide that?"

But if you're willing (or better yet "able") to use the limit as p2 -> p1 as the result for when p1 = p2 and also use the Exclusions->None option, then the line disappears. (Or you could add in a nice uniform line afterwards.)

FullSimplify[Limit[((256 ((1 - p1) p1)^(7/2) Sqrt[((1 - p2) p2)/((1 - p1) p1)])/(4 (1 - p1) p1 -
4 (1 - p2) p2)^2 - (512 (1 - p1) p1 ((1 - p2) p2)^(5/2))/(4 (1 - p1) p1 +
4 (1 - p2) p2)^2 + (256 ((1 - p1) p1)^(5/2) (1 - p2) p2 (-2 +
Sqrt[((1 - p2) p2)/((1 - p1) p1)]))/(4 (1 - p1) p1 - 4 (1 - p2) p2)^2 +
16 (1 - p1)^2 p1^2 (3 Sqrt[1/((1 - p2) p2)] - (32 ((1 - p2) p2)^(3/2))/(4 (1 - p1) p1 +
4 (1 - p2) p2)^2) - (16 ((1 - p1) p1)^(3/2) (4 (1 - p1) p1 Sqrt[((1 - p2) p2)/((1 - p1) p1)] +
4 (1 - p2) p2 (-3 + 2 Sqrt[((1 - p2) p2)/((1 - p1) p1)])))/(4 (1 - p1) p1 -
4 (1 - p2) p2) + (2 (-384 (1 - p1)^3 p1^3 + 768 (1 - p1)^2 p1^2 (1 - p2) p2 +
128 (1 - p2)^3 p2^3 Sqrt[((1 - p2) p2)/((1 - p1) p1)] - 128 (1 - p1) p1 (1 - p2)^2 p2^2 (1 + 3 Sqrt[((1 - p2) p2)/((1 - p1) p1)])))/(Sqrt[
1/((1 - p1) p1)] Sqrt[((1 - p2) p2)/((1 - p1) p1)] (4 (1 -
p1) p1 - 4 (1 - p2) p2)^2))/(8 Sqrt[(1 - p2) p2]), p2 -> p1], Assumptions -> 0 < p1 < 1]

(* -(-1 + p1) p1 *)


So we write the function as

f[p1_, p2_] := If[p1 == p2, p1 (1 - p1),
((256 ((1 - p1) p1)^(7/2) Sqrt[((1 - p2) p2)/((1 - p1) p1)])/(4 (1 - p1) p1 -
4 (1 - p2) p2)^2 - (512 (1 - p1) p1 ((1 - p2) p2)^(5/2))/(4 (1 - p1) p1 +
4 (1 - p2) p2)^2 + (256 ((1 - p1) p1)^(5/2) (1 - p2) p2 (-2 +
Sqrt[((1 - p2) p2)/((1 - p1) p1)]))/(4 (1 - p1) p1 - 4 (1 - p2) p2)^2 +
16 (1 - p1)^2 p1^2 (3 Sqrt[1/((1 - p2) p2)] - (32 ((1 - p2) p2)^(3/
2))/(4 (1 - p1) p1 + 4 (1 - p2) p2)^2) - (16 ((1 - p1) p1)^(3/
2) (4 (1 - p1) p1 Sqrt[((1 - p2) p2)/((1 - p1) p1)] + 4 (1 - p2) p2 (-3 +
2 Sqrt[((1 - p2) p2)/((1 - p1) p1)])))/(4 (1 - p1) p1 -
4 (1 - p2) p2) + (2 (-384 (1 - p1)^3 p1^3 + 768 (1 - p1)^2 p1^2 (1 - p2) p2 +
128 (1 - p2)^3 p2^3 Sqrt[((1 - p2) p2)/((1 - p1) p1)] -
128 (1 - p1) p1 (1 - p2)^2 p2^2 (1 + 3 Sqrt[((1 - p2) p2)/((1 - p1) p1)])))/(Sqrt[
1/((1 - p1) p1)] Sqrt[((1 - p2) p2)/((1 - p1) p1)] (4 (1 -
p1) p1 - 4 (1 - p2) p2)^2))/(8 Sqrt[(1 - p2) p2])]


With the contour plot given by the following:

ContourPlot[f[p1, p2], {p1, 0.5, 1}, {p2, 0.5, 1}, PlotPoints -> 100, PlotLegends -> Automatic,
FrameLabel -> {"\!$$\*SubscriptBox[\(p$$, $$1$$]\)", "\!$$\*SubscriptBox[\(p$$, $$2$$]\)"},
Exclusions -> None]


Clear["Global*"]

expr[p1_, p2_] = ((256 ((1 - p1) p1)^(7/
2) Sqrt[((1 - p2) p2)/((1 - p1) p1)])/(4 (1 - p1) p1 -
4 (1 - p2) p2)^2 - (512 (1 -
p1) p1 ((1 - p2) p2)^(5/2))/(4 (1 - p1) p1 +
4 (1 - p2) p2)^2 + (256 ((1 - p1) p1)^(5/2) (1 - p2) p2 (-2 +
Sqrt[((1 - p2) p2)/((1 - p1) p1)]))/(4 (1 - p1) p1 -
4 (1 - p2) p2)^2 +
16 (1 - p1)^2 p1^2 (3 Sqrt[
1/((1 - p2) p2)] - (32 ((1 - p2) p2)^(3/2))/(4 (1 - p1) p1 +
4 (1 - p2) p2)^2) - (16 ((1 - p1) p1)^(3/
2) (4 (1 - p1) p1 Sqrt[((1 - p2) p2)/((1 - p1) p1)] +
4 (1 - p2) p2 (-3 + 2 Sqrt[((1 - p2) p2)/((1 - p1) p1)])))/(4 (1 -
p1) p1 -
4 (1 - p2) p2) + (2 (-384 (1 - p1)^3 p1^3 +
768 (1 - p1)^2 p1^2 (1 - p2) p2 +
128 (1 - p2)^3 p2^3 Sqrt[((1 - p2) p2)/((1 - p1) p1)] -
128 (1 - p1) p1 (1 - p2)^2 p2^2 (1 +
3 Sqrt[((1 - p2) p2)/((1 - p1) p1)])))/(Sqrt[
1/((1 - p1) p1)] Sqrt[((1 - p2) p2)/((1 - p1) p1)] (4 (1 - p1) p1 -
4 (1 - p2) p2)^2))/(8 Sqrt[(1 - p2) p2]);


The function domain is

fd = FunctionDomain[expr[p1, p2], {p1, p2}]

(* 0 < p1 < 1 && 0 < p2 < 1 && p1 - p2 != 0 &&
p1 + p2 != 1 && -p1 + p1^2 - p2 + p2^2 != 0 *)


Start by simplifying the expression

expr2[p1_, p2_] = expr[p1, p2] // Simplify[#, fd] &;

LeafCount /@ {expr[p1, p2], expr2[p1, p2]}

(* {489, 357} *)


For the case p2 == p1

Limit[expr2[p1, p2], p2 -> p1]

(* -(-1 + p1) p1 *)

expr2[p_, p_] = p (1 - p);

ContourPlot[
expr2[p1, p2], {p1, 1/2, 1}, {p2, 1/2, 1},
PlotPoints -> 100,
MaxRecursion -> 2,
PlotLegends -> Automatic,
FrameLabel -> {"\!$$\*SubscriptBox[\(p$$, $$1$$]\)",
"\!$$\*SubscriptBox[\(p$$, $$2$$]\)"}] //
Quiet
`