I am currently trying to compute some integrals for a physics problem. I have a number of real parameters which are

w0 = 10^(10);
d = 3.18*10^(-26);
ra = 1.33*10^(-7);
eps = 8.85*10^{-12};
hbar = 1.05*10^{-34};
c = 3*10^8;
R = 15 ra;.

I also define a function

g[w0_] := w0^3 d^2/(3 Pi eps hbar c^3).

Now, I tried to perform the following integration numerically with NIntegrate

Int=(3 g[w0]/(4 Pi w0)) NIntegrate[(wk/(w0 - wk)) (Sin[wk R/c]/(wk R/c)+  
      Cos[wk R/c]/(wk R/c)^2 - 
      Sin[wk R/c]/(wk R/c)^3) - (wk/(w0 + wk)) (Sin[
        wk R/c]/(wk R/c) + Cos[wk R/c]/(wk R/c)^2 - 
      Sin[wk R/c]/(wk R/c)^3), {wk, 0, Infinity}, MaxRecursion -> 100]

I'm told the numerical integration is converging too slowly, but I still get an answer which is


At this point I don't really know what to make of this answer, so I try using Integrate instead of NIntegrate. I perform the integral

Int'=(3 g[w0]/(4 Pi w0)) Integrate[(wk/(w0 - wk)) (Sin[wk R/c]/(wk R/c) + 
      Cos[wk R/c]/(wk R/c)^2 - 
      Sin[wk R/c]/(wk R/c)^3) - (wk/(w0 + wk)) (Sin[
        wk R/c]/(wk R/c) + Cos[wk R/c]/(wk R/c)^2 - 
      Sin[wk R/c]/(wk R/c)^3), {wk, 0, Infinity}, 
  PrincipalValue -> True]

which gives me the answer

Int'=-48.231 + 0.00320736 I.

The first interesting thing is that


the numerical integral is giving me exactly half the answer of Integrate in the real part. Things get even more interesting by noting that one can actually compute this integral analytically with residues, which gives

Int=-(d^2 w0^3/(4 Pi hbar c^3 eps)) (Cos[w0 R/c]/(w0 R/c) - Sin[w0 R/c]/(w0 R/c)^2
     +(1 - Cos[w0 R/c])/(w0 R/c)^3),

which for the parameters I give at the start of the post gives me the answer

Int = -24.1155.

So, NIntegrate and Integrate disagree by exactly a factor of 1/2 and moreoever my exact calculation agrees exactly with NIntegrate. It seems like Integrate is getting something wrong? Could this be to do with taking the principal value? Another odd thing is that when you split the principal part integral up into the sum of two principal part integrals Integrate doesn't give you the same answer, so the principal part integration is apparently not linear? Can anybody explain what is happening here?


Later amendment:

Yet further odd behaviour is found if you Simplify the integrand. By fiddling around with the integrand one can pin down discrepencies in the results of Integrate based on how the integrand's denominator is expressed. Specifically

 2 c (c R wk Cos[(R wk)/c] + (-c^2 + R^2 wk^2) Sin[(R wk)/c]))/(
 R^3 wk (w0 - wk) (w0 + wk)), {wk, 0, Infinity}, 
 PrincipalValue -> True]

gives exactly 2/3 the result of

 2 c (c R wk Cos[(R wk)/c] + (-c^2 + R^2 wk^2) Sin[(R wk)/c]))/(
 R^3 wk (w0^2 - wk^2)), {wk, 0, Infinity}, 
 PrincipalValue -> True]

yet the only difference is that I've multiplied out the brackets in the denominator of the first integrand to get the denominator in the second integrand, that is, I used

(w0 - wk) (w0 + wk) = w0^2-wk^2.

There is clearly something wrong with Integrate for this to be happening?! If you Integrate the expression which results from taking the difference of the two integrands above then you get the message

Infinity::indet: Indeterminate expression 0. \[Infinity] encountered.

along with the answer {0.}.

If you take the difference of the integrands and then Simplify the result you get 0 obviously and then the integral evaluates to 0 identically. But evaluating the integrals separately and then taking the difference gives you a non-zero result because mathematica seems to think that the two integrals differ by a factor of 2/3, even though the integrands are identically equal! It seems that Integrate and PrincipalValue are somewhat unreliable for these problems? Does anybody know if this is to do with the PrincipalValue?


  • $\begingroup$ If I FullSimplify your integrand, Integrate and NIntegrate give 0 and -1.18039*10^14, respectively. Not sure that helps much, except maybe to indicate that the results are not all that robust? I would imagine it has something to do with the disparity in the orders of magnitude of your parameters $\endgroup$ – aardvark2012 Jun 27 '17 at 12:48
  • $\begingroup$ Did you try the Method->"ExtrapolatingOscillatory" or Method->"DoubleExponentialOscillatory" option for NIntegrate? The integrands look like they could be hard to integrate with the default strategy because of the oscillating integrand. For more information, NIntegrate Integration Strategies has some good advice how to get the most out of NIntegrate. $\endgroup$ – Thies Heidecke Jun 27 '17 at 12:56
  • $\begingroup$ Thanks for these replys. The integrand is highly oscillatory, and there is indeed some disparity in orders of magnitude of the parameters. Ordinarily I wouldn't put any stock in the results of the integrations, but it seems odd that my pen and paper calculation would give exactly what NIntegrate gives and that this is exactly half what integrate gives! I'll look into it further. $\endgroup$ – Adam Jun 27 '17 at 13:58
  • $\begingroup$ When I FullSimplify the integrand NIntegrate gives me the same answer which is -24.1155, but Integrate now gives me -72.3465 - 725278i whose real part is now exactly three times -24.1155 ?? $\endgroup$ – Adam Jun 27 '17 at 14:13
  • 3
    $\begingroup$ Have you tried Integrate but using exact rather than approximate numeric values for the parameters? When given approx values Integrate might try to discard small regions around singularities and this can have bad effects if the decision of what is "small" is not actually small enough. There could be other issues as well, since Integrate uses various symbolic methods that do not always play nice with approximate values. $\endgroup$ – Daniel Lichtblau Jun 27 '17 at 15:43

Following Daniel Lichtbalu's suggestion, note that the integral can be, and should be, performed analytically. When possible, this is much more robust than numerical integration, and the obtained results are more useful.

To do so, you should first change variables to non dimensional ones. It makes life much easier and is a good tip in general. Here the natural units for w are c/R. So we should define ω=w R/c:

integrand = (wk/(w0 - wk)) (Sin[wk R/c]/(wk R/c) + 
      Cos[wk R/c]/(wk R/c)^2 - 
      Sin[wk R/c]/(wk R/c)^3) - (wk/(w0 + wk)) (Sin[wk R/c]/(wk R/c) +
       Cos[wk R/c]/(wk R/c)^2 - Sin[wk R/c]/(wk R/c)^3);

dimensionlessintegrand = 
 Simplify[integrand R/c /. {wk -> ω c/R, w0 -> ω0 c/R}]
(*result: -((2 R (ω Cos[ω] + (-1 + ω^2) Sin[ω]))/(
 c ω (ω^2 - ω0^2)))*)

Here it is in LaTeX form, where its clear that this is a nice function:

$$\frac{2 R \left(\left(\omega ^2-1\right) \sin (\omega )+\omega \cos (\omega )\right)}{c \omega \left(\omega ^2-\text{$\omega $0}^2\right)}$$

Integrate can handle this:

Integrate[dimensionlessintegrand, {ω, 0, Infinity}, Assumptions -> ω0>0, PrincipalValue->True]

(*result: (π R (-1 +Cos[ω0] - ω0^2 Cos[ω0] + ω0 Sin[ω0]))/(c ω0^2)*)
  • $\begingroup$ Thanks. Very helpful. $\endgroup$ – Adam Jun 27 '17 at 17:23

Oscillatory integrals are always interesting. If as per @Daniel Lichtblau you replace your constants

w0 = 10^(10);
d = 3.18*10^(-26);
ra = 1.33*10^(-7);
eps = 8.85*10^{-12};  (*   <----- should change { } to ( )   *)
hbar = 1.05*10^{-34}; (*   <----- should change { } to ( )   *)
c = 3*10^8;
R = 15 ra;


w0 = 10^(10);
d = 318*10^(-28);
ra = 133*10^(-9);
eps = 885*10^(-14);
hbar = 105*10^(-36);
c = 3*10^8;
R = 15 ra;

Then both NIntegrate and Integrate give the same answer,

(* -24.1155  *)

But the Integrate also gives an unreduced answer

(224720000 (3999999982311 Cos[133/2000000] + 
    2000000 (-2000000 + 133 Sin[133/2000000])))/(26234255187 π) //

You can also FullSimplify the integrand first and it still gives you the same 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.