# Partial Derivative after Numerical Integration of a Complicated Expression with Singularity at Zero

I want to generate Fig no. 2,3,4 from https://arxiv.org/pdf/2212.08237

with the help of the following equations. In Fig 3 t_(opt) is the time Q_T reaches its max value Q_T(opt) for each T.

Here is my failed attempt

   Needs["NumericalCalculus"]
p = 0.13; \[Alpha] = 0.77; \[Mu] = 0.5;
\[CapitalGamma][t_?NumericQ, temp_?NumericQ] :=
NIntegrate[
p Exp[-0.5 \[Mu]^2 k^2]/(
Sqrt[k^4 + 2 \[Alpha] k^2] (k^2 + 2 \[Alpha]))
Sin[(t Sqrt[k^4 + 2 \[Alpha] k^2])/(2 temp)]^2 Coth[Sqrt[
k^4 + 2 \[Alpha] k^2]/(2 temp)], {k, 0, \[Infinity]},
Method -> "LocalAdaptive", PrecisionGoal -> 12];
data = ListDensityPlot[
Table[temp^2*(ND[\[CapitalGamma][t, temp], temp, temp])^2/(
Exp[2 \[CapitalGamma][t, temp]] - 1), {t, 1, 25, 2}, {temp,
0.1*10^-9, 10^-9, 0.1*10^-9}], InterpolationOrder -> 2]


Edit

Needs["NumericalDifferentialEquationAnalysis"];
Needs["NumericalCalculus"]
p = 0.13; \[Alpha] = 0.77; \[Mu] = 0.5;
points = gg[[All, 1]];
weight = gg[[All, 2]];
a\[CapitalGamma][t_?NumericQ, temp_?NumericQ, k_?NumericQ] :=
p Exp[-\[Mu]^2 k^2/2]/(
Sqrt[k^4 + 2 \[Alpha] k^2] (k^2 + 2 \[Alpha]))
Sin[(t Sqrt[k^4 + 2 \[Alpha] k^2])/(2 temp)]^2 Coth[Sqrt[
k^4 + 2 \[Alpha] k^2]/(2 temp)];

data = ParallelTable[{t, temp,
temp^2 (Evaluate[
weight.Table[
Evaluate[
Derivative[0, 1, 0][a\[CapitalGamma]][t, temp, k]], {k,
points}]])^2/(
Exp[2 weight.Table[
Evaluate[a\[CapitalGamma][t, temp, k]], {k, points}]] -
1)}, {t, 1, 25}, {temp, 0.1*10^-9, 1.0*10^-9, 0.1*10^-9}];


1. Do you realize that you have about 1.97506*10^10 oscillations per unit change in k for smallish k (gets worse as k increases)? I think the Levin rule is the only hope here.

2. And that the amplitude rapidly decays below what is computable by Mathematica? You'll probably be better off integrating k from 0 to 20. (Some of the computations of the integrand were taking 25+ GB with 8+ GB in active RAM in the kernel.

So see if this seems reasonable:

\[CapitalGamma][t_?NumericQ, temp_?NumericQ] :=
NIntegrate[
p  Exp[-1/2  \[Mu]^2  k^2]/(Sqrt[
k^4 + 2  \[Alpha]  k^2]  (k^2 +
2  \[Alpha]))  Sin[(t  Sqrt[
k^4 + 2  \[Alpha]  k^2])/(2  temp)]^2  Coth[
Sqrt[k^4 + 2  \[Alpha]  k^2]/(2  temp)], {k, 0, 20},
Method -> {"LevinRule",
Method -> {"GaussKronrodRule", "Points" -> 21}}
, MinRecursion -> 2, MaxRecursion -> 40,
WorkingPrecision -> MachinePrecision, PrecisionGoal -> 10];

start = SessionTime[];
Block[{p = 13/100, \[Alpha] = 77/100, \[Mu] = 5/10},
Table[With[{gamma = \[CapitalGamma][t, temp]},
temp^2*((\[CapitalGamma][t,
temp (1 + Sqrt@Sqrt[$$MachineEpsilon])] - 2 gamma + \[CapitalGamma][t, temp (1 + Sqrt@Sqrt[$$MachineEpsilon])])/(temp^2*
Sqrt[\$MachineEpsilon]))^2/(Exp[2  gamma] - 1)
], {t, 1, 25; 3, 2}, {temp, 1/10*10^-9, 10^-9; 2/10*10^-9,
1/10*10^-9}]
]
stop = SessionTime[];

NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.


NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.

< repeated>

(*  {{6.98765*10^24, 1.85143*10^24}, {4.22587*10^24, 1.11481*10^24}}  *)

stop - start

(*  26.240706  <-- long time for just 4 values *)


Using a higher WorkingPrecision made things go way too slow for me. Maybe it would improve things, if you have the time.

• A higher PrecisionGoal seems to be the problem. With PrecisionGoal -> 6 it works with a sufficiently high MaxRecursion. Commented Jul 1 at 19:14

Maybe this can get you going, not sure it is all correct, the graph looks a bit different, but maybe with more points and fixing ColorFunction and PlotLegends and so on you can replicate the figure from the paper. All runs in less than 15 seconds or so.

{p, a, m} = Rationalize[{0.13, 0.77, 0.5}]; G[(t_)?NumericQ, temp_] :=
NIntegrate[(p*Exp[-0.5*m^2*k^2]*Sin[(t*Sqrt[k^4 + 2*a*k^2])/(2*temp)]^2*
Coth[Sqrt[k^4 + 2*a*k^2]/(2*temp)])/(Sqrt[k^4 + 2*a*k^2]*(k^2 + 2*a)), {k, 0, Infinity},
MaxRecursion -> 44, PrecisionGoal -> 6];
Do[g[t] = (Interpolation[#1, Method -> "Spline", InterpolationOrder -> 3] & )[
ParallelTable[{temp, G[t, temp]}, {temp, 0.1/10^9, 10^(-9), 0.1/10^9}]];
gd[t] = Derivative[1][g[t]], {t, 1, 25}];
r = ParallelTable[temp^2*(gd[t][temp]^2/(Exp[2*G[t, temp]] - 1)),
{temp, 0.1/10^9, 10^(-9), 0.1/10^9}, {t, 1, 25}];
ListDensityPlot[r, ColorFunction -> "TemperatureMap", PlotLegends -> Automatic]
`
• Thank you for the answer @RolfMertig. But the Q_T value is coming different than the paper. Also how can I find max Q_T for each T and t? Commented Jul 2 at 7:24
• Anyway I tried this approach of integration in "Edit" from @AlexTrounev mathematica.stackexchange.com/questions/284427/… but I do not know what is going wrong! If anybody can look into it it will be very helpful. Commented Jul 2 at 7:24
• y-axis in the plot of Q_T is also coming different ! @RolfMertig. Commented Jul 3 at 6:56