# An intergral which can be evaluated by NIntegrate only with symbolic preprocessing

I am calculating the integral from $s=-\infty$ to $s+\infty$ of the following function

f[s_, t_] := Exp[-I*s*t]/(s + I*10^-6)


which is essentially the Fourier representation of Heaviside step function, modulo $2 \pi \mathrm{i}$ factor. I define the following two functions

fAnalytic[s_, t_] := f[s, t]
fNumerical[s_?NumericQ, t_?NumericQ] := f[s, t]


and, after setting $t=0.1$, NIntegrate is able to grecefully integrate the first one, but fails on the second one, as per the following screenshot Clearly, then, in the first case numerical integration succeeds because NIntegrate can do some sort of symbolic preprocessing before crunching numbers.

I am interested, however, in evaluating the integral in a fully numerical way. The motivation is that I am dealing with a function with with a similar behaviour, whose form I know only numerically.

I have tried performing variable changes, for instance mapping the infinite integration region to the interval $[0,1]$, but it does not change the situations, fAnalytic is always calculated correctly, while the numerical integration for fNumerical always fails.

Edit: as requested in a comment I tried many different integrations methods, but no one gets close to the real result:

iMethods = {"Trapezoidal", "AdaptiveMonteCarlo", "GlobalAdaptive","LocalAdaptive", "DoubleExponential", "AdaptiveMonteCarlo","AdaptiveQuasiMonteCarlo", "DuffyCoordinates", "Oscillatory","BooleRule", "ClenshawCurtisRule", "GaussBerntsenEspelidRule","GaussKronrodRule", "LobattoKronrodRule", "LobattoPeanoRule","MultiPanelRule", "NewtonCotesRule", "PattersonRule","SimpsonThreeEightsRule", "TrapezoidalRule"};

NIntegrate[fNumeric[s, 0.1], {s, -Infinity, Infinity},Method -> #] & /@ iMethods

• You could try different integration rules and compare their results and whether warning messages appear or not. – Marius Ladegård Meyer Apr 28 '17 at 9:32
• @MariusLadegårdMeyer: many thanks for you answer. I tried many different integration methods, and I edited my question to reflect that. However, that doesn't seem to help, it looks like no method gets even close to the real answer. – zakk Apr 28 '17 at 9:45
• You're going to have to do some symbolic processing, imo. Either by the human user or by Mathematica. – Michael E2 Apr 28 '17 at 10:19
• Michael E2: my problem is that I have a similar integral, in which the integrand function f is effectively a black box, hindering any kind of symbolic processing... – zakk Apr 28 '17 at 11:16
• Yes, I understood that. But oscillatory strategies exist because such integrals are difficult and perhaps cannot be computed via generic methods. You say your integral is similar. That suggests perhaps something is known about the form inside the black box. The period of the oscillatory kernel, perhaps. The truncation error of Ryhor's approach can be seen to be on average to be around 2/(t m) at most, where the integration interval is {s, -m, m}. – Michael E2 Apr 28 '17 at 13:18

Is this sufficient?

First the antiderivative via NDSolve (takes a long time, 103 sec.). The substitution maps $(-\pi/2, \pi/2)$ to $(-\infty, \infty)$. The errors are just warnings that integration is stopping at (or just before, really) the singularities of tangent.

sol = NDSolve[{y'[t] == fNumerical[Tan[t], 0.1] Sec[t]^2, y == 0}, y, {t, -Pi/2, Pi/2}]


NDSolve::ndsz: At t == -1.5708, step size is effectively zero; singularity or stiff system suspected.

NDSolve::ndsz: At t == 1.5707962897386536, step size is effectively zero; singularity or stiff system suspected.

(*  {{y -> InterpolatingFunction[{{-1.5708, 1.5708}}, <>]}}  *)


Find difference at endpoints to get integral:

y@y["Domain"] /. sol // Flatten // Differences // First
(*  0. - 6.28319 I  *)


This gives the range of sampling (any numerical techinique is going to have finite sampling):

Flatten[y["Domain"] /. sol] // Tan
(*  {-2.6986*10^7, 2.6986*10^7}  *)

Length@Flatten[y["Grid"] /. sol]  (* number of _steps_ *)
(*  3815939  *)


I think NDSolve, with its local error estimation, is probably the best one can do for an oscillatory integral whose oscillatory nature cannot be analyzed.

• Thanks for you answer. I am afraid this is not sufficient. Even if SymbolicProcessing is set to 0, Mathematica still does some analytical computation before crunching numbers... If I use your method to integrate fAnalytic I get the same answer as you get, on the other hand, if I use it to integrate fNumerical I get -167.721 - 211.129 I. – zakk Apr 28 '17 at 10:13
• @zakk I changed my answer. – Michael E2 Apr 28 '17 at 12:30
• @zakk You're welcome. :) – Michael E2 Apr 28 '17 at 13:19

You have a very sharp behavior near s=0, so I think that you have to divide integration range to take into account for this pecularities.

NIntegrate[Im[fNumerical[s, 0.1]], {s, -10000, -0.0001}] +
NIntegrate[Im[fNumerical[s, 0.1]], {s, 0.0001, +10000}] +
NIntegrate[Im[fNumerical[s, 0.1]], {s, -0.0001, +0.0001}]


-6.282058269

• This works, because Infinity is replaced by a large number, not because a separate integration is performed near s == 0. Compare, for instance, NIntegrate[fNumerical[s, 0.1], {s, -10000, 10000}, MaxRecursion -> 100]`. – bbgodfrey Apr 28 '17 at 10:46
• Thanks @Ryhor for you solution. However, I would like something that is working without a cutoff, unlike the 10000 cutoff you have used. I might consider using a cutoff if nothing else works, though... Thanks again! – zakk Apr 28 '17 at 13:13