0
$\begingroup$

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;
r = 0.1; (*radial coordinate*)

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): enter image description here

(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:

enter image description here

(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)

enter image description here

and r =0.8

enter image description here

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)

$\endgroup$
0
1
$\begingroup$

(Extended comment, not an answer.)

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

enter image description here

$\endgroup$

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.