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
):
(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)
(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}]
]
Answered by Anton Antonov on May 9, 2021
Get help from others!
Recent Answers
Recent Questions
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP