# a problem in numerical integration in Mathematica

I am trying to numerically compute an integral of the following type in Mathematica. I can show that this integral can be computed analytically. Still, I want to make sure Mathematica can also give me the same result using NIntegrate since I will compute some complicated integral of the same kind numerically later. Here is the integral,

$$\int_0^{3} d\omega e^{-i\omega t}\omega\Big(\int_0 ^{\infty}d\omega' \frac{e^{i\omega' t}}{\omega'-\omega-i\epsilon}-\int_0 ^{\infty}d\omega' \frac{e^{-i\omega' t}}{\omega'+\omega+i\epsilon}\Big),$$

where $$0\leq\epsilon\ll 1$$ after the computation of the integral. One can change the variable of integration in the second integral (inside the parenthesis) such that $$\omega' \to -\omega'$$ and then combine the two integrals with the result $$\Big(\int_0 ^{\infty}d\omega' \frac{e^{i\omega' t}}{\omega'-\omega-i\epsilon}-\int_0 ^{\infty}d\omega' \frac{e^{-i\omega' t}}{\omega'+\omega+i\epsilon}\Big)=\int_{-\infty} ^{\infty}d\omega' \frac{e^{i\omega' t}}{\omega'-\omega-i\epsilon}.$$

Then, we can use complex integration to solve this integral. The answer is $$2\pi i e^{i\omega t}$$. Plugging this result into the original double integral, we find $$\int_0^{3} d\omega e^{-i\omega t}\omega (2\pi i e^{-i\omega t})=2\pi i\int_0^3 d\omega \omega .$$

So, we can see that this integral is $$t$$ independent. However, if we use the NIntegrate and compute the original integral numerically, we get a t-dependent result. I am going to make this clear with an example. We know that $$2\pi i\int_0^3 d\omega \omega =0. + 28.2743 I$$. Next, we try numerical integration. The original integral can be written as $$\int_0^{3} d\omega e^{-i\omega t}\omega\Big(\int_0 ^{\infty}d\omega' \frac{e^{i\omega' t }}{\omega'-\omega-i\epsilon}-\int_0 ^{\infty}d\omega' \frac{e^{-i\omega' t}}{\omega'+\omega+i\epsilon}\Big)$$. The first integral can be written as a principal value integral plus a Dirac delta part: $$\int_0^{3} d\omega e^{-i\omega t}\omega \Bigg(PV(\int_0 ^{\infty}d\omega' \frac{e^{i\omega' t }}{\omega'-\omega})+I\pi e^{i\omega t}-\int_0 ^{\infty}d\omega' \frac{e^{-i\omega' t}}{\omega'+\omega+i\epsilon}\Bigg)$$

We can simply set $$\epsilon=0$$ in the last integral since it doesn't have any singularity. I computed these integrals using NIntegrate and for t=1, the result was 9.22151244253655*^-13 + 28.274333882309378 I. For t=10, the answer is 1.83864*10^-7 + 28.2743 I. The answers are t dependent and also the real part is not equal to the correct answer. As you can see, the correct answer is a pure imaginary number, but the numerical answer has a small real part. I have included my code here. I will appreciate any advice or suggestion. Thanks!

I1[t] is the last integral (in the last equation). I have performed the integral over y first. That gave me an exponential integral.

I1[t_]:=NIntegrate[Exp[-I*x*t](E^(I t x) ExpIntegralEi[-I t(x+Infinity)]-E^(I t x) ExpIntegralEi[-I t (x + 0)])*(x), {x, 0, 3}]


The contribution of the second term ($$i\pi e^{i\omega t}$$) to the integral has been computed by

I2 = NIntegrate[I*Pi*x, {x, 0, 3}]


This part is t independent and is equal to 0. + 14.1372 I. The principal value integral is computed by

I3[ts_]:=NIntegrate[NIntegrate[Exp[-I*x*ts]*(Exp[I*y*ts]/(y - x))*(x), {y, 0, x, 3},PrincipalValue -> True], {x, 0, 3}]


Since we $$3, there is no singularly in the integral so there is no need to use the principal value. For that piece, I used

I4[ts_]:=NIntegrate[Exp[-i*x*t] (x)*(E^(I t x) ExpIntegralEi[I t (-x + Infinity)] -E^(I t x) ExpIntegralEi[I t (-x + 3)]), {x, 0, 3}].

• Please include also Mathematica code (in input form) next to the equations in your question. May 24 at 17:57
• You can format inline code and code blocks by selecting the code and clicking the {} button above the edit window. The edit window help button ? is useful for learning how to format your questions and answers. You may also find the meta Q&A, How to copy code from Mathematica so it looks good on this site, helpful May 24 at 19:14
• @MichaelE2 Thank you! I'll do that. May 24 at 19:15
• I fixed it already in this case. (You had indented just one of the lines of code. You need to indent all 4 space. But there's an easy way to do that, the {} button.) May 24 at 19:18

Here's a straightforward way:

i10[t_, w_?NumericQ] := NIntegrate[
Exp[I v t]/(v - w), {v, -Infinity, w, Infinity},
Method -> "PrincipalValue",
WorkingPrecision -> MachinePrecision];
i1[t_] := NIntegrate[E^(-I t w) w*i10[t, w], {w, 0, 3},
WorkingPrecision -> MachinePrecision]

i1[1/2]

(*  1.69022*10^-8 + 14.1372 I  *)


Notes: (1) I think your 28.27 is off by a factor of 2. (2) I don't think ExpIntegralEi represents your integral.

You can check the interior integral this way:

Assuming[t > 0 && w > 0,
Sqrt[2 Pi] FourierTransform[1/(v - w), v, t]
]

(*  I E^(I t w) \[Pi]  *)

i2[t_] :=
NIntegrate[E^(-I t w) w*(I E^(I t w) \[Pi]), {w, 0, 3},
WorkingPrecision -> MachinePrecision]

i2[1/2]

(*  0. + 14.1372 I  *)


And i2[t] may be computed exactly, of course.

• Thank you for your answer! I have included all codes now. Also, that extra 14.1372 *I is coming from the Dirac delta part. @MichaelE2 May 24 at 20:04
• The problem is that the real part of the numerical answer is not equal to 0. Also, it is changing with t. I am wondering how I can get rid of that. That seems to be a numerical error. Thank you! @MichaelE2 May 24 at 20:08
• @HadamardN2 "the real part of the numerical answer is not equal to 0": Depends on what you mean by equal. It's equal up to round-off error. Set WorkingPrecision -> 24 (in both places) and Chop[] the result. May 24 at 20:37
• Thank you very much! May 25 at 2:39