TransWikia.com

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

Mathematica Asked by j.foobles on May 9, 2021

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)

One Answer

(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

Answered by Anton Antonov on May 9, 2021

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP