Спасибо за предложение, алгоритм удалось организовать но почему-то два разных способа решения приводят к двум разным ответам.
Первый способ такой:
Код:
eqn = 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;
With[{gensol = Solve[eqn , qp]},
Block[{\[Theta] = Pi/6, ma = 980, mp = 139,
j},(*subs vals when gensol is evaluated*)
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)]));
sols = gensol]];
qpC13 = Compile[{{q, _Complex}, {m5, _Complex}},
Evaluate[qp /. sols[[3]]],
RuntimeOptions -> "EvaluateSymbolically" -> False] ;
Plot3D[Re@qpC13[q, m5], {q, 0, 10000}, {m5, 0, 2000},
AxesLabel -> Automatic]
Картинка которая получается:
Второй способ немного иной:
Код:
eqn = j -
Sqrt[q^2 + qp^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;
With[{gensol = Solve[eqn, qp]},
Block[{\[Theta] = Pi/6, ma = 980, mp = 139, j},
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)]));
sols = gensol]];
opt = Experimental`OptimizeExpression[qp /. sols[[3]]];
Block[{\[Theta] = Pi/6, f},
f[q0_, m50_] :=
Block[{q = SetPrecision[q0, 40], m5 = SetPrecision[m50, 40], qp},
qp = First@opt;
Re@qp /; Im@qp == 0];
Plot3D[f[q, m5], {q, 0, 10000}, {m5, 0, 1000}, MaxRecursion -> 3,
AxesLabel -> Automatic]]
Картинка получается другая:
Не подскажите, почему так происходит? Уже лысину прочесал - не пойму в чём дело.