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

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*)

• You can find numerical x,p,e with NMinimize[{Abs[f[x, p, e]]^2, x >= 0, p >= 0, e >= 0}, {x, p, e}] . This gives {1.75528*10^-21, {x -> 1.99839, p -> 0.010061, e -> 0.0544189}} Aug 16, 2020 at 12:38
• ^ note with the above you need to redefine f[x_]:=... as f[x_, p_, e_]:=... . I don't think it's likely you'll get a nice symbolic solution out of Reduce or Solve, or a general solution for x for arbitrary p,e because it's a difficult mixed radical high order polynomial equation. Aug 16, 2020 at 12:56
• @flinty - Due to the complicated form of the function, I recommend using arbitrary-precision rather than machine precision for the minimization. arg = NMinimize[{Abs[f[x, p, e]]^2, x > 0, p > 0, e > 0}, {x, p, e}, WorkingPrecision -> 75][[2]] // N evaluates to {x -> 5.07773, p -> 0.00319756, e -> 0.00363597} and the function value f[x, p, e] /. arg evaluates to 5.81432*10^-19 Aug 16, 2020 at 14:23
• f = f[x] // Simplify; Solve[f == 0 && e >= 0 && p >= 0, Reals, Method -> Reduce] so it is not easy to solve x by e and p Aug 16, 2020 at 14:47
• Would you be interested in an InterpolatingFunction approximation to the solution? If so, what domain of interest would p and e have? Aug 16, 2020 at 16:52

## 2 Answers

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}]


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.