Mathematica: ускорение численного интегрирования

[illuminates]

Столкнулся с бедой, что NIntegrate и Plot3D работают крайне медленно. Один график строится полдня. Приведу минимальный пример, в котором отображен основной алгоритм расчёта.

Сначала я численно решаю уравнение:

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[{m = 5.5, M = 300, Nc = 3, 
         c = ScientificForm[-44687.3983417778], 
         b = ScientificForm[161593.81818181818], 
         k1 = ScientificForm[16.485010961790245], 
         k2 = ScientificForm[-13.131344420001051], ma, mp, 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)])], 
          ma = Sqrt[-2*(M^2 - 2*(3*k1 + k2)*
                     (Sqrt[(c + M^2 + 2*m5^2)/(2*(k1 + k2))] + 
                          m*(b/(2*(c + M^2 + 2*m5^2))))^2 - c + 
                   2*m5^2)], mp = Sqrt[(2*b*m)/
                (Sqrt[(c + M^2 + 2*m5^2)/(2*(k1 + k2))] + 
                   m*(b/(2*(c + M^2 + 2*m5^2))))]}; 
        sols = gensol]]; 
qpC12 = Compile[{{q, _Complex}, {m5, _Complex}, 
         {\[Theta], _Complex}}, Evaluate[qp /. sols[[2]]], 
       RuntimeOptions -> "EvaluateSymbolically" -> False]; 
qp32 = Re[qpC12[q, m5, \[Theta]]]; 

Далее я провожу численное интегрирование:

  qqq[q_?NumericQ, m5_?NumericQ] := 
  NIntegrate[( 
     q^2 qp32^2 Sin[\[Theta]]^2 )/((2 (2 \[Pi])^2)  Sqrt[-(2 q Cos[\
\[Theta]] qp32) + qp32^2 + q^2]) , {\[Theta], 0, Pi/2}];

И строю график:

 Plot3D[qqq[q, m5], {q, 0.8, 150}, {m5, 0, 150}] 

Данный пример занимает 80 секунд. Читал советы, что ускорить расчёт можно через первообразную, но у меня функция нескольких параметров и я не сильно понимаю как это сделать для неё.
Пробовал так:

 f1[q_, m5_, \[Theta]_] := (q^2 qp32[q, 
      m5, \[Theta]]^2 Sin[\[Theta]]^2)/((2 (2 \[Pi])^2) Sqrt[-(2 q \
Cos[\[Theta]] qp32[q, m5, \[Theta]]) + qp32[q, m5, \[Theta]]^2 + q^2])

f2sol = ParametricNDSolve[{D[f2[q, m5, \[Theta]], \[Theta]] == 
    f1[q, m5, \[Theta]], D[f2[q, m5, \[Theta]], \[Theta]] == 0}, 
  f2, {\[Theta], 0, Pi/2}, {q, m5}, 
  Method -> {"ParametricSensitivity" -> None}]
Plot3D[f2sol[q, m5, \[Theta]], {q, 0.8, 150}, {m5, 0, 150}]

Но такой алгоритм не получилось завести. Для справки, для одной переменной, очень легко завести такой алгоритм. Например, простейший пример:

 f[x_] := Sin[x^2]
iF = NDSolveValue[{F2'[x] == f[x], F2[0] == 0}, F2, {x, 0, 5}]
Plot[iF[x], {x, 0, 5}]

Буду благодарен любым советам и предложениям! Возможно стоит попробовать сделать всё на питоне?

[Markiyan Hirnyk]

При построении графика получаю сообщение об ошибке

NIntegrate::inumr: The integrand (0.00832331 Re[(0.202667 Cos[\[Theta]] (<<1>>))/(2.11694*10^-7+<<5>>)-<<1>> <<1>>+1/2 Sqrt[<<1>>]]^2 Sin[\[Theta]]^2)/(\[Sqrt](0.657182 -1.62134 <<1>> Re[0.202667 Cos[<<1>>] Power[<<2>>] Plus[<<11>>]-<<1>>+1/2 Power[<<2>>]]+Re[0.202667 Cos[<<1>>] Power[<<2>>] Plus[<<11>>]-1/2 <<1>>+1/2 Power[<<2>>]]^2)) has evaluated to non-numerical values for all sampling points in the region with boundaries {{0,1.5708}}.

Задумайтесь над этим.

[illuminates]

Markiyan Hirnyk
к сожалению, пока не думается на счёт этого, можно подсказку?

[arseniiv]

Решил-таки глянуть, что тут. Итак, зачем, к примеру, вы использовали ScientificForm при присвоении значений переменным блока? Это формат вывода, на значения он никак не влияет.

illuminates в сообщении #1273535 писал(а):

к сожалению, пока не думается на счёт этого, можно подсказку?

Вычислите просто qqq[1, 1], к примеру. Сообщение буквально говорит, что подынтегральное выражение (integrand) вычислолось во что-то нечисловое на всех точках некоторого промежутка. Если посмотреть на выдаваемое невычисленное выражение

NIntegrate[(1^2 qp32^2 Sin[\[Theta]]^2)/((2 (2 \[Pi])^2) Sqrt[-(2 Cos[\[Theta]] qp32) + qp32^2 + 1^2]), {\[Theta], 0, \[Pi]/2}]

видно, что вы напортачили с определениями, и qp32 не было вычислено вовремя.

-- Сб дек 09, 2017 21:59:58 --

А, нет, не так. Само qp32 вычисляется — но в громадное символьное выражение, зависящее от переменных, и притом ещё и весьма долго. Вы бы эти параметры выше определили, а то ваш минимальный пример на самом деле не пример, а хотелось бы, если вам интересны какие-то соображения, иметь возможность запускать то же, что и у вас.

[illuminates]

arseniiv
Извините, я хотел чтоб на форуме код максимально коротко и поэтому использовал ScientificForm, но действительно, оказалось что если скопировать предоставленный мной код, он просто не будет работать. Давайте ещё раз. Мой рабочий код это:

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[{m = 5.5, M = 300, Nc = 3, c = \!\(\*
TagBox[
InterpretationBox[
RowBox[{"\"\<-4.46874\>\"", "*", 
SuperscriptBox["10", "\"\<4\>\""]}],
-44687.3983417778,
AutoDelete->True],
ScientificForm]\), b = \!\(\*
TagBox[
InterpretationBox[
RowBox[{"\"\<1.61594\>\"", "*", 
SuperscriptBox["10", "\"\<5\>\""]}],
161593.81818181818`,
AutoDelete->True],
ScientificForm]\), k1 = \!\(\*
TagBox[
InterpretationBox[
RowBox[{"\"\<1.6485\>\"", "*", 
SuperscriptBox["10", "\"\<1\>\""]}],
16.485010961790245`,
AutoDelete->True],
ScientificForm]\), k2 = \!\(\*
TagBox[
InterpretationBox[
RowBox[{"\"\<-1.31313\>\"", "*", 
SuperscriptBox["10", "\"\<1\>\""]}],
-13.131344420001051`,
AutoDelete->True],
ScientificForm]\), ma, mp, 
    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)])), 
    ma = \[Sqrt](-2 (M^2 - 
          2 (3 k1 + k2) (Sqrt[(c + M^2 + 2 m5^2)/(2 (k1 + k2))] + 
             m b/(2 (c + M^2 + 2 m5^2)))^2 - c + 2 m5^2)), 
    mp = Sqrt[
     2 b m (Sqrt[(c + M^2 + 2 m5^2)/(2 (k1 + k2))] + 
        m b/(2 (c + M^2 + 2 m5^2)))^-1]};
   sols = gensol]];
qpC12 = Compile[{{q, _Complex}, {m5, _Complex}, {\[Theta], _Complex}},
    Evaluate[qp /. sols[[2]]], 
   RuntimeOptions -> "EvaluateSymbolically" -> False] ;
qp32 = Re@qpC12[q, m5, \[Theta]];

далее

qqq[q_?NumericQ, m5_?NumericQ] := 
 NIntegrate[(q^2 qp32^2 Sin[\[Theta]]^2)/((2 (2 \[Pi])^2) Sqrt[-(2 q \
Cos[\[Theta]] qp32) + qp32^2 + q^2]), {\[Theta], 0, Pi/2}];

И далее

Plot3D[qqq[q,m5],{q,0.8,150},{m5,0,150}]

Постройка графика занимает 80 секунд. И мне хотелось бы понять, как можно ускорить код.

[Vince Diesel]

Опция PlotPoints -> 10 в Plot3D уполовинивает время (возможно, за счет некоторого огрубления графика).

[illuminates]

Vince Diesel
Спасибо удалось ускорить расчет графика до 6 часов с 12-ти. Конечно тоже плохо, ибо графиков нужно пару десятков...

[Vince Diesel]

То ли у меня компьютер поновее, то ли версия математики - график из примера строится за 17 секунд вместо 80.

Если у процессора несколько ядер, можно распараллелить вычисления, математика позволяет. Например, командой ParallelTable. Кроме того, функция из примера сильно (неограниченно?) растет в одном из углов прямоугольника. Возможно, программа берет много точек (чтобы определить масштаб по вертикали?). Если этот угол не очень то нужен, как вариант, вырезать его из области опцией RegionFunction и посмотреть, не сократится ли время расчета. Еще, если на компьютере установлен компилятор Си, в Compile можно выставить опцию CompilationTarget -> "C". Иногда это ускоряет вычисления, но далеко не всегда.