I'm trying to compute the following. Let $$f(t) := \exp\left(\frac{1}{t^2 - 1}\right)1_{-1 < t < 1}.$$ This function looks like (it's the standard example of the bump function on Wikipedia).


Next, consider the Fourier transform of $f$ given by $$\widehat{f}(\xi) = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}f(t)e^{-i t\omega}\, dt.$$ My question is, how does one numerically compute the integral $$\int_{-\infty}^{\infty}|\widehat{f}(\omega)|\, d\omega?$$

I'm doing this via the commands

   f[t_] = Piecewise[{{Exp[1/(t^2 - 1)], -1 < t < 1}}]
   NIntegrate[Abs[NFourierTransform[f[t], t, ω]], {ω, -Infinity,Infinity}]

This outputs

  NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in ω near {ω} = {12.4765}. NIntegrate obtained 1.1329578878244777` and 0.000052199655872313335` for the integral and error estimates.

Is there a way to get rid of this issue?

  • $\begingroup$ The warning goes away if you add the option PrecisionGoal -> 4, because the numbers in the message imply four digits of precision. If that is not accurate enough, you can set MaxRecursion -> 20 higher. (I believe the problem is that it is an oscillatory integral, but the function is numerical, perhaps limiting NIntegrate from using a more sophisticated strategy.) $\endgroup$ – Michael E2 Aug 4 '16 at 20:27

Three methods:

All simplify the problem by using the even symmetry of the transform to change the interval to {0, Infinity}, which saves some work for NIntegrate[]. (NIntegrate[] would do this automatically, if it could determine the integrand was symmetric.)

1. PrecisionGoal -> 4

2 * NIntegrate[
   Abs[NFourierTransform[f[t], t, ω]], {ω, 0, Infinity},
    PrecisionGoal -> 4] // AbsoluteTiming
(*  {28.6487, 1.13297} *)

2. MaxRecursion -> 20

2 * NIntegrate[
   Abs[NFourierTransform[f[t], t, ω]], {ω, 0, Infinity},
    MaxRecursion -> 20] // AbsoluteTiming
(*  {86.8169, 1.13296}  *)

3. Feed zeros to NIntegrate[] as singular points.

The derivative is undefined at the zeros of f[t], so they are "weak" singularities. Again, if NIntegrate[] could find them, it could adjust itself, but f[t] essentially behaves like a numerical black box.

{poszeros} = 
     NDSolve[{x'[ω] == Sin[ω], x[0] == 0., 
       WhenEvent[Re@NFourierTransform[f[t], t, ω] == 0, 
        Sow[ω]]}, x, {ω, 0, 200}, 
      Method -> "Extrapolation", MaxStepSize -> 3]; // AbsoluteTiming
(*  {29.1027, Null}  *)

2 * NIntegrate[
   Abs[NFourierTransform[f[t], t, ω]], {ω, 0, 
    Sequence @@ Sort[poszeros], Infinity}] // AbsoluteTiming
(*  {28.8482, 1.13296}  *)
| improve this answer | |

One possibility is to use NDSolve to create a FourierTransform function, and then use NDSolve again to integrate the function. To do this, first note that the function is even, so that the FourierTransform can be converted to a FourierCosTransform. This is useful because the end point is problematic for NDSolve. So, here is a Fourier transform interpolating function:

ft = NDSolveValue[
    {D[FT[t, ω], t] == Sqrt[2/π] Piecewise[{{Exp[1/(t^2-1)], t<1}}] Cos[t ω], FT[0, ω]==0}, 
    FT[1, ω], 
    {t, 0, 1}, 
    {ω, -500, 500}, 
]; //AbsoluteTiming

{0.550217, Null}

I used a cutoff of $\left| \omega \right| <500$ for the frequency domain. Now, we could use NIntegrate to compute the $L^1$ norm, but this is rather slow. It is much faster to use NDSolveValue again:

NDSolveValue[{g'[ω] == 2 Abs @ ft, g[0] == 0}, g[500], {ω, 0, 500}] //AbsoluteTiming

{0.291507, 1.13296}

in good agreement with @MichaelE2's answer.

| improve this answer | |

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.