# NIntegrate blowing up/behaving weirdly at the turning point of the integrand

I'm performing (what should be) a straightforward numerical integration (Fourier transform) of a function with no poles / singularities (at least in a particular parameter regime):

<< FourierSeries

ClearAll["Global*"]

l = 1;
ϵ = 10^-4;

p = 0;
P = 80;

Data1 = ConstantArray[Null, {P, 2}];

While[p < P, p++;

w = -10 + 20 p/P;

σ = (
2 (-2 l^2 r^2 Sin[s/(2 Sqrt[l^-2 - r^2])]^2 + (1 - l^2 r^2) 2 Sinh[
s/(2 Sqrt[l^-2 - r^2]) - I ϵ]^2))/l^2;

RF = NInverseFourierTransform[-1/(4 π^2)/σ, s,
w, Method -> {"LocalAdaptive", "SymbolicProcessing" -> 0},
MinRecursion -> 6] // Chop;

Data1[[p, 1]] = w;
Data1[[p, 2]] = RF;

];

ListLinePlot[Data1, PlotRange -> All,
LabelStyle -> {FontSize -> 16, FontFamily -> "CMU Serif", Black},
PlotLabel -> StringForm["r = .", {r // N}]] // Print;


When the parameter r is small, the numerics seem to work fine/as expected, producing plots like this (r = 0.1):

(x-axis is the energy, y-axis is the transition rate of a two-level system - this shows a roughly thermal Planck spectrum).

When increasing r beyond some critical point, the numerics seem to blow up, giving a weird answer. For example, r = 0.8 yields:

(see especially the magnitude of the y-axis).

I plotted the integrand for these two values below. The numerics seem to blow up when the function bifurcates from having one turning point to two turning points:

r=0.1 (integrand plotted as a function of s)

and r =0.8

Why does NIntegrate not like this seemingly inconspicuous change in behaviour of the integrand? Any help would be greatly appreciated! (If there's an analytic solution, even better :P)

Maybe this behavior is expected?

If "plain" Method->"GlobalAdaptive" is used the messages indicate multiple problems encountered by NIntegrate. Try to do your investigations using rationals for the parameters and carefully selecting precision and accuracy goals.

Here is example code:

AbsoluteTiming[
Table[
Block[{},
l = 1;
\[Epsilon] = 10^-4;
p = 0;
P = 80;

Data1 = ConstantArray[Null, {P, 2}];

While[p < P,
p++;
w = -10 + 20 p/P;
\[Sigma] = (2 (-2 l^2 r^2 Sin[s/(2 Sqrt[l^-2 - r^2])]^2 + (1 -
l^2 r^2) 2 Sinh[
s/(2 Sqrt[l^-2 - r^2]) - I \[Epsilon]]^2))/l^2;
RF = NInverseFourierTransform[-1/(4 \[Pi]^2) E^(-I w s)/\[Sigma],
s, w, Method -> {"GlobalAdaptive", "SymbolicProcessing" -> 0,
"SingularityHandler" -> None, "MaxErrorIncreases" -> 100000},
MinRecursion -> 6, MaxRecursion -> 100, WorkingPrecision -> 30,
PrecisionGoal -> 4, AccuracyGoal -> 4];
Data1[[p, 1]] = w;
Data1[[p, 2]] = RF;
];
ListLinePlot[Data1, PlotRange -> All,
LabelStyle -> {FontSize -> 16, FontFamily -> "CMU Serif", Black},
PlotLabel -> StringForm["r = .", {r // N}]]
],
{r, Range[7, 8, 1/10]/10}]
]