# A problem with FindRoot

I have the following

With[{p2 = 0.75}, FindRoot[1/(64 (p1 - p2)^2 (-(-1 + p2) p2)^(
3/2) (-1 + p1 + p2)^2) (1 - p1) p1 (-128 (-(-1 + p1) p1)^(3/2) (-1 + p2)^2 p2^2 + (384 (-1 + p1)^3 p1^3)/Sqrt[1/(p2 - p2^2)] + 128 Sqrt[-(-1 + p1) p1] (-1 + p2)^3 p2^3 (-3 + 2 Sqrt[((-1 + p2) p2)/((-1 + p1) p1)]) + 128 Sqrt[1/(p1 - p1^2)] (3 (-1 + p1) p1 (-1 + p2)^3 p2^3 - (-1 + p2)^4 p2^4 + (2 (-(-1 + p2) p2)^(7/2))/Sqrt[1/(p1 - p1^2)] + ((-1 + p2)^4 p2^4)/(((-1 + p2) p2)/((-1 + p1) p1))^(3/2) + 3 (-1 + p1)^4 p1^4 Sqrt[((-1 + p2) p2)/((-1 + p1) p1)] - 6 (-1 + p1)^4 p1^4 (((-1 + p2) p2)/((-1 + p1) p1))^(3/2)) + 4 (1 - p1) p1 (-128 (-(-1 + p2) p2)^(5/2) - 96/(1/(p2 - p2^2))^(5/2)) + 256 (-1 + p1)^2 p1^2 ((-(-1 + p2) p2)^(3/2) + 3/(1/(p2 - p2^2))^(3/2))) == 0, {p1, 0.5}]]


which produces the solution

{p1 -> 0.882142}


Now, I try to generate a series of such solution. I start by creating the equation

eqn = 1/(64 (p1 - p2)^2 (-(-1 + p2) p2)^(3/2) (-1 + p1 + p2)^2) (1 - p1) p1(-128 (-(-1 + p1) p1)^(3/2) (-1 + p2)^2 p2^2 + (384 (-1 + p1)^3 p1^3)/Sqrt[1/(p2 - p2^2)] + 128 Sqrt[-(-1 + p1) p1] (-1 + p2)^3 p2^3 (-3 + 2 Sqrt[((-1 + p2) p2)/((-1 + p1) p1)]) + 128 Sqrt[1/(p1 - p1^2)] (3 (-1 + p1) p1 (-1 + p2)^3 p2^3 - (-1 + p2)^4 p2^4 + (2 (-(-1 + p2) p2)^(7/2))/Sqrt[1/(p1 -p1^2)] + ((-1 + p2)^4 p2^4)/(((-1 + p2) p2)/((-1 + p1) p1))^(3/2) + 3 (-1 + p1)^4 p1^4 Sqrt[((-1 + p2) p2)/((-1 + p1) p1)] - 6 (-1 + p1)^4 p1^4 (((-1 + p2) p2)/((-1 + p1) p1))^(3/2)) +4 (1 - p1) p1 (-128 (-(-1 + p2) p2)^(5/2) - 96/(1/(p2 - p2^2))^(5/2)) + 256 (-1 + p1)^2 p1^2 ((-(-1 + p2) p2)^(3/2) + 3/(1/(p2 - p2^2))^(3/2))) == 0 /. p1 -> p1[p2]


and attempt to setup an initial condition so that I can then use NDSolveValue. However, something goes wrong. The following initial condition should output {p1 -> 0.882142} (or should it?). It seems to do nothing.

p11 = p1[0.75] /. FindRoot[eqn /. p2 -> .75, {p1[.75], .75}]


What am I doing wrong?

• I don't see where you have p1 defined as a function that would take a value. MMa naturally doesn't know what p1[.75] is. – Bill Watts Dec 30 '18 at 22:39
• @BillWatts isn't that what p1 -> p1[p2] does? (Perhaps an obvious newbie question). – user120911 Dec 30 '18 at 22:43

Clear["Global*"]

expr = ((1 -
p1) p1 (-128 ((1 - p1) p1)^(3/2) (-1 +
p2)^2 p2^2 + (384 (-1 + p1)^3 p1^3)/Sqrt[1/(p2 - p2^2)] +
128 Sqrt[(1 - p1) p1] (-1 + p2)^3 p2^3 (-3 +
2 Sqrt[((-1 + p2) p2)/((-1 + p1) p1)]) +
128 Sqrt[1/(p1 -
p1^2)] (3 (-1 +
p1) p1 (-1 + p2)^3 p2^3 - (-1 +
p2)^4 p2^4 + (2 ((1 - p2) p2)^(7/2))/
Sqrt[1/(p1 -
p1^2)] + ((-1 +
p2)^4 p2^4)/(((-1 + p2) p2)/((-1 + p1) p1))^(3/2) +
3 (-1 + p1)^4 p1^4 Sqrt[((-1 + p2) p2)/((-1 + p1) p1)] -
6 (-1 + p1)^4 p1^4 (((-1 + p2) p2)/((-1 + p1) p1))^(3/2)) +
4 (1 - p1) p1 (-128 ((1 - p2) p2)^(5/2) - 96/(1/(p2 - p2^2))^(5/2)) +
256 (-1 + p1)^2 p1^2 (((1 - p2) p2)^(3/2) +
3/(1/(p2 - p2^2))^(3/2))))/(64 (p1 - p2)^2 ((1 - p2) p2)^(3/
2) (-1 + p1 + p2)^2) // Simplify;


The domain over which expr is real is given by FunctionDomain

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

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


Simplifying expr within the specified domain

expr2 = expr // FullSimplify[#, fd] &;


The effect of the simplification is significant

LeafCount /@ {expr, expr2}

(* {335, 184} *)


Since you are interested in expr2 == 0, it is equivalent to Numerator@expr2 == 0

LeafCount[Numerator@expr2]

(* 166 *)


The simplified form plots more efficiently and precisely

ContourPlot[Numerator@expr2 == 0, {p1, 0, 1}, {p2, 0, 1},
Exclusions -> {p1 == p2, p1 + p2 == 1},
MaxRecursion -> 5,
FrameLabel -> {p1, p2}] // Quiet


**EDIT:**

Clear["Global*"]

expr[p1_, p2_] = ((1 -
p1) p1 (-128 (-(-1 + p1) p1)^(3/2) (-1 +
p2)^2 p2^2 + (384 (-1 + p1)^3 p1^3)/Sqrt[1/(p2 - p2^2)] +
128 Sqrt[-(-1 + p1) p1] (-1 + p2)^3 p2^3 (-3 +
2 Sqrt[((-1 + p2) p2)/((-1 + p1) p1)]) +
128 Sqrt[1/(p1 -
p1^2)] (3 (-1 +
p1) p1 (-1 + p2)^3 p2^3 - (-1 +
p2)^4 p2^4 + (2 (-(-1 + p2) p2)^(7/2))/
Sqrt[1/(p1 -
p1^2)] + ((-1 + p2)^4 p2^4)/(((-1 + p2) p2)/((-1 + p1) p1))^(3/
2) + 3 (-1 + p1)^4 p1^4 Sqrt[((-1 + p2) p2)/((-1 + p1) p1)] -
6 (-1 + p1)^4 p1^4 (((-1 + p2) p2)/((-1 + p1) p1))^(3/2)) +
4 (1 - p1) p1 (-128 (-(-1 + p2) p2)^(5/2) -
96/(1/(p2 - p2^2))^(5/2)) +
256 (-1 + p1)^2 p1^2 ((-(-1 + p2) p2)^(3/2) +
3/(1/(p2 - p2^2))^(3/2))))/(64 (p1 - p2)^2 (-(-1 + p2) p2)^(3/
2) (-1 + p1 + p2)^2);

expr2[p1_, p2_] =
FullSimplify[expr[p1, p2], FunctionDomain[expr[p1, p2], {p1, p2}]];

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

(* {341, 184} *)

Manipulate[
Plot[expr2[p1, p2], {p1, 0.9, 1}], {p2, 0.9001, 0.9999, 0.0001,
Appearance -> "Labeled"}]
(* *)


• Your tricky simplification without Exclusions and increased Workingprecision ContourPlot[Numerator@expr2 == 0, {p1, 0, 1}, {p2, 0, 1}, MaxRecursion -> 5, FrameLabel -> {p1, p2}, WorkingPrecision -> 50] // Quiet evaluates the contourlines including the points {0,0},{1,0],{0,1},{1,1} ! – Ulrich Neumann Dec 31 '18 at 8:47
• Bob Hanlon, as you will notice from my edit, there are solutions closer to p1=1 (or p1 = 0) than what the contour plot show. Can your solution show the needed level of precision? – user120911 Jan 1 '19 at 16:29
• See comment by @UlrichNeumann (which you should have already received) above your comment. – Bob Hanlon Jan 1 '19 at 16:40
• Oh, I missed that. Thank you. – user120911 Jan 1 '19 at 17:28

If I understand your question right you want to solve an equation depending on tw0 variables p1,p2.

expr[p1,p2]:=((1 - p1) p1 (-128 ((1 - p1) p1)^(3/2) (-1 + p2)^2 p2^2 + (
384 (-1 + p1)^3 p1^3)/Sqrt[1/(p2 - p2^2)] +
128 Sqrt[(1 - p1) p1] (-1 + p2)^3 p2^3 (-3 +
2 Sqrt[((-1 + p2) p2)/((-1 + p1) p1)]) +
128 Sqrt[1/(
p1 - p1^2)] (3 (-1 + p1) p1 (-1 + p2)^3 p2^3 - (-1 +
p2)^4 p2^4 + (2 ((1 - p2) p2)^(7/2))/Sqrt[1/(
p1 - p1^2)] + ((-1 +
p2)^4 p2^4)/(((-1 + p2) p2)/((-1 + p1) p1))^(3/2) +
3 (-1 + p1)^4 p1^4 Sqrt[((-1 + p2) p2)/((-1 + p1) p1)] -
6 (-1 + p1)^4 p1^4 (((-1 + p2) p2)/((-1 + p1) p1))^(3/2)) +
4 (1 - p1) p1 (-128 ((1 - p2) p2)^(5/2) - 96/(1/(p2 - p2^2))^(
5/2)) + 256 (-1 + p1)^2 p1^2 (((1 - p2) p2)^(3/2) +
3/(1/(p2 - p2^2))^(3/2))))/(64 (p1 - p2)^2 ((1 - p2) p2)^(
3/2) (-1 + p1 + p2)^2)


ContourPlot shows the possible solutions p1[p2] in the range 0<p1,p2<1

ContourPlot[Evaluate[expr[p1,p2] == 0], {p1, 0, 1}, {p2, 0, 1}, Exclusions -> None,MaxRecursion -> 5, FrameLabel -> {p1, p2},WorkingPrecison->50]


• You understand correctly. – user120911 Dec 30 '18 at 23:01
• I am a little surprised by the "gap" at the top-right / bottom-left. Is that due to the level of precision used? If so, how do I amend your code to increase precision? – user120911 Dec 30 '18 at 23:14
• Perhaps WorkingPrecision->.. might help. Alternatively you could evaluate the contourline using NDSolve! – Ulrich Neumann Dec 30 '18 at 23:21