Interesting discrepencies between integrate functions

Question

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

Int=-24.1155.

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

Re[Int']=2Int,

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?

Thanks!

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

Integrate[(
 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

Integrate[(
 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?

Thanks.

Answer 1

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)*)

Answer 2

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;

With

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 π) //
  N

You can also FullSimplify the integrand first and it still gives you the same answer.