# Integration of a high oscillatory function

I want to get the numerical result of the integration below:

Integrate[Exp[I* s*1000.0]*1/(1 + (10)^2*s^2), {s, 0, Infinity}]
(*5.8434816785318*10^-45 - 5.204153095728048*10^26 I*)


The result is wrong in my opinion because when the Exp[I*...] have a high frequency this tends to make the integral zero or finite, this is called wave rotating approximation in physics.

I tried:

 NIntegrate[Exp[I* s*1000.0]*1/(1 + (10)^2*s^2), {s, 0, Infinity},
MaxRecursion -> 1000]

(* lots of erros ...........
-0.000021305291177576493 + 0.0010002002551860643 I*)


Because the function 1/(1 + (10)^2*s^2) decays fast enough I think it isn't necessary to integrate up to infinity so I change it to a numerical number like 10^3, using bigger upper limits almost has no effect on the result so it's a good approximation and I consider this as the right answer:

NIntegrate[Exp[I* s*1000.0]*1/(1 + (10)^2*s^2), {s, 0, 10^3},
MaxRecursion -> 1000]
(*1.1026966286690337*10^-14 + 0.0010002002470217504 I*)


Questions:

Why does MMA give the wrong result in the first example by using Integrate?

I really couldn't get the errors in the second example and I am not familiar with either Integrate's options or that of NIntegrate. Is it possible to get the correct result when the upper bound of Integration is infinity?

I think this is merely a matter of precision:

Integrate[Exp[I*s*1000]*1/(1 + 10^2*s^2), {s, 0, Infinity}]
N[%, 16]

π/(20 E^100) + 1/20 I Sqrt[π] MeijerG[{{1/2}, {}}, {{1/2, 1/2}, {0}}, 2500]
0.*10^-45 + 0.0010002002407240688 I


If you want to use NIntegrate, then a quick search in the document shows that "LevinRule" is your friend:

NIntegrate[Exp[I*s*1000]*1/(1 + 10^2*s^2), {s, 0, Infinity}, Method -> "LevinRule"]

-5.3668*10^-18 + 0.0010002 I


A little different from the result generated by Integrate but it's acceptable. If you add a WorkingPrecision -> (* a number >= 16* ) then you'll get better result.

Regarding to the comment:

int[x_] = Integrate[Exp[I*s*x]*1/(1 + 10^2*s^2), {s, 0, Infinity}, Assumptions -> x ∈ Reals]
int
N[%, 16]

1/20 E^(-(Abs[x]/10)) π + 1/20 I Sqrt[π]
MeijerG[{{1/2}, {}}, {{1/2, 1/2}, {0}}, x^2/400] Sign[x]
π/(20 E^100) + 1/20 I Sqrt[π] MeijerG[{{1/2}, {}}, {{1/2, 1/2}, {0}}, 2500]
5.843481678531469*10^-45 + 0.001000200240724069 I


The assumption is necessary. I'm using v9.0.1

• If you change 1000in the exponential function to something else like x and define a function out of it, and then put x then the result is wrong. I think MMA solves the integrate wrong when we use Integrate. I don't know in this case how the precision matters. – MOON Jan 7 '15 at 13:42
• @yashar See my edit. – xzczd Jan 7 '15 at 14:03
• If I put a real number in int[1000.0], I get: 5.843481678531469*10^-45 - 3.312359715464843*10^26 I However, int returns: 5.843481678531469*10^-45 + 0.001000200240724069 I – MOON Jan 7 '15 at 14:29
• note that MiejerG eval breaks down dramatically if you do not request extended precision: N[MeijerG[{{1/2}, {}}, {{1/2, 1/2}, {0}}, 2500],16]-> 0.01128 ... N[MeijerG[{{1/2}, {}}, {{1/2, 1/2}, {0}}, 2500]]-> -3.7376*10^27 obviously an implementation bug, er.. limitation. – george2079 Jan 7 '15 at 18:41
• @yashar A bit of algebra leads to the exponential integral: Exp[I*s*1000]*Apart[1/(1 + (10)^2*s^2) /. s -> u/I] /. u -> I s. Expand and integrate each part. Use exact coefficients as @xzczd recommends above. Numerically it is better behaved. As for N[expr] vs N[expr, prec]: The first does machine precision quickly without precision tracking, and the second tracks the precision and ensures an accurate result. The first way, depending on the algorithm, can accumulate a lot of error. – Michael E2 Jan 7 '15 at 23:53