2

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?

Michael E2
  • 235,386
  • 17
  • 334
  • 747
MOON
  • 3,864
  • 23
  • 49

1 Answers1

7

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[1000]
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

xzczd
  • 65,995
  • 9
  • 163
  • 468
  • 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 07 '15 at 13:42
  • @yashar See my edit. – xzczd Jan 07 '15 at 14:03
  • If I put a real number in int[1000.0], I get: 5.84348167853146910^-45 - 3.31235971546484310^26 I

    However, int[1000] returns: 5.843481678531469*10^-45 + 0.001000200240724069 I

     – MOON Jan 07 '15 at 14:29
  • @yashar Yeah, that's expected result. 1000 is a number with infinity precision, while 1000.0 is a number with MachinePrecison. There're many posts about this issue in this site, have a search :) – xzczd Jan 07 '15 at 14:35
  • Could you please give me some of those topics about this issue? I have no clue what search for. Moreover, Limit[int[x], x -> 1000.0] returns 5.843481678531469*10^-45, it omits the imaginary part. – MOON Jan 07 '15 at 15:24
  • 1
    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 07 '15 at 18:41
  • 1
    @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 07 '15 at 23:53
  • @yashar Just search MachinePrecision and WorkingPrecision, and you'll see lots of posts about how and why precision dramatically influence the result. I used to collect just 4 of them here. – xzczd Jan 08 '15 at 03:13
  • @MichaelE2 is this "quickly" concept documented. I can see no indication that N[expr] should differ (in this case wildly) from N[expr,$MachinePrecision] – george2079 Jan 08 '15 at 21:51
  • @george2079 "The main advantage of using the built-in floating-point capabilities of your computer is speed." From tutorial/MachinePrecisionNumbers; see also tutorial/ArbitraryPrecisionNumbers. Try timing N on expr = Sin[Range[10^5]]; the $MachinePrecision calculation takes about 25 times longer. Also compare the number of digits in FullForm and ByteCount of N[Sin[1]] with those of N[Sin[1], $MachinePrecision]. – Michael E2 Jan 08 '15 at 23:33