4
$\begingroup$

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.

$\endgroup$
  • 2
    $\begingroup$ 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. $\endgroup$ – Simon Woods Apr 9 '17 at 18:50
  • $\begingroup$ 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
  • $\begingroup$ @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) $\endgroup$ – illuminates Apr 9 '17 at 19:34
  • $\begingroup$ @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. $\endgroup$ – Michael E2 Apr 9 '17 at 20:17
6
$\begingroup$

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]

Mathematica graphics

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

Mathematica graphics

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

Mathematica graphics

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]

Mathematica graphics


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

Mathematica graphics

The higher setting for MaxRecursion bumps the plotting time up to 14.5 sec.

$\endgroup$
  • 1
    $\begingroup$ 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? $\endgroup$ – illuminates Apr 9 '17 at 21:42
  • $\begingroup$ @illuminates Do you mean ConditionalExpression? $\endgroup$ – Michael E2 Apr 9 '17 at 22:01
  • $\begingroup$ 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. $\endgroup$ – illuminates Apr 10 '17 at 18:51
  • $\begingroup$ @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. $\endgroup$ – Michael E2 Apr 11 '17 at 2:24
  • $\begingroup$ 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 . $\endgroup$ – illuminates Apr 11 '17 at 10:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.