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

- Peter Machado on Why fry rice before boiling?
- haakon.io on Why fry rice before boiling?
- Joshua Engel on Why fry rice before boiling?
- Jon Church on Why fry rice before boiling?
- Lex on Does Google Analytics track 404 page responses as valid page views?

Recent Questions

- How can I transform graph image into a tikzpicture LaTeX code?
- How Do I Get The Ifruit App Off Of Gta 5 / Grand Theft Auto 5
- Iv’e designed a space elevator using a series of lasers. do you know anybody i could submit the designs too that could manufacture the concept and put it to use
- Need help finding a book. Female OP protagonist, magic
- Why is the WWF pending games (“Your turn”) area replaced w/ a column of “Bonus & Reward”gift boxes?

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