4
$\begingroup$

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

NIntegrate works on fAnalytic but fails on fNumerical.

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
Improve this question
$\endgroup$
5
  • $\begingroup$ You could try different integration rules and compare their results and whether warning messages appear or not. $\endgroup$
    – Marius Ladegård Meyer
    Apr 28, 2017 at 9:32
  • $\begingroup$ @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. $\endgroup$
    – zakk
    Apr 28, 2017 at 9:45
  • $\begingroup$ You're going to have to do some symbolic processing, imo. Either by the human user or by Mathematica. $\endgroup$
    – Michael E2
    Apr 28, 2017 at 10:19
  • $\begingroup$ 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... $\endgroup$
    – zakk
    Apr 28, 2017 at 11:16
  • 1
    $\begingroup$ 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}. $\endgroup$
    – Michael E2
    Apr 28, 2017 at 13:18

2 Answers 2

Reset to default
5
$\begingroup$

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] == 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.

Improve this answer
$\endgroup$
3
  • $\begingroup$ 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. $\endgroup$
    – zakk
    Apr 28, 2017 at 10:13
  • $\begingroup$ @zakk I changed my answer. $\endgroup$
    – Michael E2
    Apr 28, 2017 at 12:30
  • $\begingroup$ @zakk You're welcome. :) $\endgroup$
    – Michael E2
    Apr 28, 2017 at 13:19
2
$\begingroup$

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

Improve this answer
$\endgroup$
2
  • $\begingroup$ 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]. $\endgroup$
    – bbgodfrey
    Apr 28, 2017 at 10:46
  • $\begingroup$ 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! $\endgroup$
    – zakk
    Apr 28, 2017 at 13:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.