# How to solve and plot large equation with parameters?

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)

• Mathematica's inability to find a closed-form solution for qp doesn't necessarily mean you did something wrong. Maybe there isn't one. You might need to substitute numerical values for q and m5 before solving for qp. – Simon Woods Apr 9 '17 at 18:50
• Solve[expr, qp] returns in a few of seconds for me, if I leave the parameters as undefined symbols. – Michael E2 Apr 9 '17 at 19:05
• @MichaelE2 Yes. But unfortunately it is too large answer for ploting this. It is possible create table (qp,p) for this? (m5 - constant and let = 100) – illuminates Apr 9 '17 at 19:34
• @user2975438 See my answer. It seemed faster to show you than to try to describe it. Note the reason it seems not to plot for you is probably because of complex values. – Michael E2 Apr 9 '17 at 20:17

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 = ExperimentalOptimizeExpression[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.

• Thank you very much. May I ask you about add condition? It is imposible to put q^2 + qp^2 - q qp Cos[[Theta]]>0? – illuminates Apr 9 '17 at 21:42
• @illuminates Do you mean ConditionalExpression? – Michael E2 Apr 9 '17 at 22:01
• Yes. But if I write Plot3D[ConditionalExpression[Re@qpC[q, m5], q^2 + qp^2 - q qp Cos[\[Theta]] > 0], {q, 0, 10000}, {m5, 0, 1000}, AxesLabel -> Automatic] that it doesn't work. – illuminates Apr 10 '17 at 18:51
• @illuminates Isn't the condition equivalent to Im[qp] == 0? Both your condition and the condition Im[qp] < 10^-10 give the same plots as the Chop one above. – Michael E2 Apr 11 '17 at 2:24
• Yes, you are right. I have another question. Is it possible to create functions where is qpC? For example if I want to create function: 'z=qpC^2+qpC^3+m5' and after that to plot it: Plot3D[Re@z[q, m5], {q, 0, 10000}, {m5, 0, 1000}, AxesLabel -> Automatic]`, for some reason it doesn't plot . – illuminates Apr 11 '17 at 10:44