How to solve and plot large equation with parameters?

Question

I have a very large expression:

j - Sqrt[q^2 + qp^2 - 
   2 q qp Cos[\[Theta]]] - \[Sqrt](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

where

\[Theta] = Pi/6; ma = 980; mp = 139;
 j = \[Sqrt](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)])) 

And $qp$ is real and positive. I want to find qp as function $q>0$ and $m5>0$ (if it is imposible then function only $q>0$ and $m5=constant>0$). I atemted to solve this the next way:

Solve[ j - Sqrt[
     q^2 + qp^2 - 
      2 q qp Cos[\[Theta]]] - \[Sqrt](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 && qp > 0 , qp, 
  Reals];

Mathematica hangs after two hour calculation. Where I am wrong? After colculation I want to plot function form: $z=q^2+m5^2$ (in reality this is more complicated)

---

Link to Wolfram Community version.

Answer 1

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.