1
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Bill Watts Dec 30 '18 at 22:39
  • $\begingroup$ @BillWatts isn't that what p1 -> p1[p2] does? (Perhaps an obvious newbie question). $\endgroup$ – user120911 Dec 30 '18 at 22:43
4
$\begingroup$
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

enter image description here

**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"}]
(* *)

enter image description here

$\endgroup$
  • 2
    $\begingroup$ 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} ! $\endgroup$ – Ulrich Neumann Dec 31 '18 at 8:47
  • $\begingroup$ 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? $\endgroup$ – user120911 Jan 1 '19 at 16:29
  • $\begingroup$ See comment by @UlrichNeumann (which you should have already received) above your comment. $\endgroup$ – Bob Hanlon Jan 1 '19 at 16:40
  • $\begingroup$ Oh, I missed that. Thank you. $\endgroup$ – user120911 Jan 1 '19 at 17:28
4
$\begingroup$

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]

enter image description here

$\endgroup$
  • $\begingroup$ You understand correctly. $\endgroup$ – user120911 Dec 30 '18 at 23:01
  • $\begingroup$ 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? $\endgroup$ – user120911 Dec 30 '18 at 23:14
  • $\begingroup$ Perhaps WorkingPrecision->.. might help. Alternatively you could evaluate the contourline using NDSolve! $\endgroup$ – Ulrich Neumann Dec 30 '18 at 23:21

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.