It is feasible to solve the general equation with symbolic parameters and then substitute their numeric values.
eqn = j - Sqrt[q^2 + qp^2 - q qp Cos[θ]] -
√(qp^2 + 1/2 (16 m5^2 + ma^2 + mp^2 -
Sqrt[(-(16 m5^2) - ma^2 - mp^2)^2 -
4 (ma^2 mp^2 - 16 m5^2 qp^2)])) == 0;
With[{gensol = Solve[eqn, qp]},
Block[{θ = Pi/6, ma = 980, mp = 139, j}, (* subs vals when gensol is evaluated *)
j = √(q^2 + 1/2 (16 m5^2 + ma^2 + mp^2 +
Sqrt[(-(16 m5^2) - ma^2 - mp^2)^2 -
4 (ma^2 mp^2 - 16 m5^2 q^2)]));
sols = gensol
]];
Compiling the long expression is a good thing to do, if machine precision is sufficient. Note the function is complex in some parts. Even where the solution appears to be real, the imaginary parts can be great enough that Plot3D
won't draw them. One can use Re
, Im
or Chop
to plot the parts as desired. Chop
will result in a plot of the real part only where the imaginary part is less than 10^-10
in magnitude.
There are four solutions. Below we plot the 4th.
qpC = Compile[{{q, _Complex}, {m5, _Complex}},
Evaluate[qp /. sols[[4]]], (* index = 1, 2, 3, or 4 *)
RuntimeOptions -> "EvaluateSymbolically" -> False]
Plot3D[Re@qpC[q, m5],
{q, 0, 10000}, {m5, 0, 1000}, AxesLabel -> Automatic]

Plot3D[Im@qpC[q, m5],
{q, 0, 10000}, {m5, 0, 1000}, AxesLabel -> Automatic]

Plot3D[Chop@qpC[q, m5],
{q, 0, 10000}, {m5, 0, 1000}, AxesLabel -> Automatic]

One can plot the expression without compiling, but it takes 104 sec. instead of just 0.14 sec.
Plot3D[Evaluate[Chop@qp /. sols[[4]]],
{q, 0, 10000}, {m5, 0, 1000}, AxesLabel -> Automatic]

Update
This is pretty nice. Optimizing the expression saves a lot of time and memory! And it's fast, much faster than Simplify
, which ran longer than I would wait.
opt = Experimental`OptimizeExpression[qp /. sols[[4]]]; // AbsoluteTiming
(* {0.015315, Null} *)
ByteCount /@ {qp /. sols[[4]], opt}
(* {2858680, 25912} *)
Plotting it is not as fast as plotting the compiled version (3.7 sec. vs. 0.14 sec), but much faster than qp /. sols[[4]]
(176 sec.).
Another advantage is one can evaluate with arbitrary precision. Higher precision makes a nicer plot, too.
Here's an example that implements the OP's desired constraint via a pattern restriction on a function. One problem is that the function in the constraint is sometimes complex-valued. This causes an error in the inequality. Hends I added a condition Im[qp] == 0
. I would have used WorkingPrecision
instead of SetPrecision
; but somehow machine precision numbers got passed and numerical errors resulted that gave messages.
Block[{θ = Pi/6, f},
f[q0_, m50_] :=
Block[{q = SetPrecision[q0, 40], m5 = SetPrecision[m50, 40], qp},
qp = First@opt;
Re@qp /; Im@qp == 0 && q^2 + Re@qp^2 - q Re@qp Cos[θ] > 0];
Plot3D[f[q, m5], {q, 0, 10000}, {m5, 0, 1000},
MaxRecursion -> 3, AxesLabel -> Automatic]
]

The higher setting for MaxRecursion
bumps the plotting time up to 14.5 sec.
Solve[expr, qp]
returns in a few of seconds for me, if I leave the parameters as undefined symbols. $\endgroup$ – Michael E2 Apr 9 '17 at 19:05