Best way to solve stupidly complicated parametric non linear equation with reduce/solve etc

Question

I have a parametric nonlinear function which is literally a nightmare. I know roots exist, they are real and the two parameters, p,e are both positive. What I was expecting from mathematica was to get a solution (in form of a root, no closed form) but even letting the program working all night I surrendered. I can't understand if it's me framing the problem not in an efficient way, it's my pc which needs serious upgrades or the problem that is simply too hard for methods like Reduce or Solve. If the case is the latter, I guess I'm doomed... Any hint towards the other two? Thanks for the help.

My attempts and the equation:

f[x_]:=(1/(8 (p^2+x^2)^3))p^2 (-2 p^6+p^5 (4-8 x)+2 p^3 (3-8 x) x^2-6 p x^4+p^4 (80000+2 x-9 x^2)+2 p^2 x (40000+60000 x+x^2-5 x^3)-3 x^3 (-80000+40000 x+x^3)+(4 Sqrt[10] e x (p^2+x^2)^2 (2 p x^3+p^4 (-2+3 x)+2 p^3 x (-3+4 x)+x^2 (-80000+40000 x+x^3)+2 p^2 (40000-60000 x-x^2+2 x^3)))/Sqrt[-e p^2 (-1+x) x^2 (p^2+x^2)^2 (-40000+p^2+2 p x+x^2)])
Reduce[f[x]==0 && x>=0 &&p>=0 && e>=0,x,Reals] (*stuck running*)
Solve[f[x]==0 && x>=0 &&p>=0 && e>=0,x,Reals] (*stuck running*)

Answer 1

If you are satisfied with an approximate answer, then you could try using NDSolveValue. Your function:

f[x_] := (1/(8 (p^2+x^2)^3))p^2 (-2 p^6+p^5 (4-8 x)+2 p^3 (3-8 x) x^2-6 p x^4+p^4 (80000+2 x-9 x^2)+2 p^2 x (40000+60000 x+x^2-5 x^3)-3 x^3 (-80000+40000 x+x^3)+(4 Sqrt[10] e x (p^2+x^2)^2 (2 p x^3+p^4 (-2+3 x)+2 p^3 x (-3+4 x)+x^2 (-80000+40000 x+x^3)+2 p^2 (40000-60000 x-x^2+2 x^3)))/Sqrt[-e p^2 (-1+x) x^2 (p^2+x^2)^2 (-40000+p^2+2 p x+x^2)])

To use NDSolveValue, we need to know a boundary condition. For example, here's the value of x when p is 1:

x1 = x /. Block[{p=1}, First @ Solve[f[x] == 0, x]]

Root[256006399839996 + 1023948800640 e + (255942399200020 - 3071999998080 e) #1 + (511955203840004 + 2304217598880 e) #1^2 + (1279846402079976 - 2048025605760 e) #1^3 + (-960279984200093 + 1024166395840 e) #1^4 + (192129617159615 + 4095897593600 e) #1^5 + (-384151995680463 - 512486397600 e) #1^6 + (-3455678391520375 + 2047846414080 e) #1^7 + (2880427177039798 - 332788480 e) #1^8 + (-576225624399594 - 1024102385280 e) #1^9 + (43199520578 + 256064008800 e) #1^10 + (-14402879738 - 25598080 e) #1^11 + (-359829 + 12802240 e) #1^12 + 360087 #1^13 + (9 + 160 e) #1^14 + 9 #1^15 &, 1]

Now, we can use NDSolveValue:

sol = NDSolveValue[
    {
    D[f[x[p,e]]==0, p], x[1, e] == x1},
    x,
    {p,.1,100},
    {e,.1,10000},
    MaxStepFraction->.0005,
    PrecisionGoal->10
]; //AbsoluteTiming

{19.2292, Null}

Check some random samples:

Block[{p = 50, e = 200}, f[sol[p, e]]]
Block[{p = 10, e = 2000}, f[sol[p, e]]]

6.42413*10^-9

8.0893*10^-9

Visualization:

Plot3D[sol[p,e], {p,.1,100}, {e,.1,10000}]

Answer 2

FWIW, here is half a solution: Take the numerator, rationalize it so that it is a polynomial, and find the roots. What is left to do is to select the ones which are positive when p and e are positive. This step takes along time (if it can be done at all), except when specific numeric values are given to p and e.

num = Simplify[ff, x >= 0 && p >= 0 && e >= 0] // Together // 
      Numerator // Simplify[#, x >= 0 && p >= 0 && e >= 0] & // 
    FactorList // #[[-1, 1]] & // Simplify

(*
p Sqrt[-e (-1 + x) (-40000 + p^2 + 2 p x + x^2)] (2 p^6 + 6 p x^4 + 
    p^5 (-4 + 8 x) + 2 p^3 x^2 (-3 + 8 x) + 
    p^4 (-80000 - 2 x + 9 x^2) + 3 x^3 (-80000 + 40000 x + x^3) + 
    2 p^2 x (-40000 - 60000 x - x^2 + 5 x^3)) - 
 4 Sqrt[10]
   e (2 p x^5 + 4 p^3 x^3 (-1 + 2 x) + p^6 (-2 + 3 x) + 
    2 p^5 x (-3 + 4 x) + p^2 x^3 (-80000 - 2 x + 5 x^2) + 
    x^4 (-80000 + 40000 x + x^3) + 
    p^4 (80000 - 120000 x - 4 x^2 + 7 x^3))
*)

Check there are two terms (the first obviously containing the radical):

Length@num

(*  2  *)

rat = num*MapAt[-# &, num, 1] // Expand // Simplify

(*
e (p^2 (-1 + x) (-40000 + p^2 + 2 p x + x^2) (2 p^6 + 6 p x^4 + 
      p^5 (-4 + 8 x) + 2 p^3 x^2 (-3 + 8 x) + 
      p^4 (-80000 - 2 x + 9 x^2) + 3 x^3 (-80000 + 40000 x + x^3) + 
      2 p^2 x (-40000 - 60000 x - x^2 + 5 x^3))^2 + 
   160 e (2 p x^5 + 4 p^3 x^3 (-1 + 2 x) + p^6 (-2 + 3 x) + 
      2 p^5 x (-3 + 4 x) + p^2 x^3 (-80000 - 2 x + 5 x^2) + 
      x^4 (-80000 + 40000 x + x^3) + 
      p^4 (80000 - 120000 x - 4 x^2 + 7 x^3))^2)
*)

Solve[rat == 0, x]

(* <15 Root objects> *)

These roots contain extraneous solutions, and one seems to need to substitution numerical values for the parameters to work with them. If such an approach is useful, then perhaps substituting the parameters in f[x] and dealing with the resulting equation f[x] == 0 might be better.