1
$\begingroup$

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<y<\infty$, 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}].
$\endgroup$
4
  • 2
    $\begingroup$ Please include also Mathematica code (in input form) next to the equations in your question. $\endgroup$
    – Domen
    May 24, 2023 at 17:57
  • $\begingroup$ 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 $\endgroup$
    – Michael E2
    May 24, 2023 at 19:14
  • $\begingroup$ @MichaelE2 Thank you! I'll do that. $\endgroup$
    – HadamardN2
    May 24, 2023 at 19:15
  • $\begingroup$ 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.) $\endgroup$
    – Michael E2
    May 24, 2023 at 19:18

1 Answer 1

1
$\begingroup$

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.

$\endgroup$
4
  • $\begingroup$ 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 $\endgroup$
    – HadamardN2
    May 24, 2023 at 20:04
  • $\begingroup$ 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 $\endgroup$
    – HadamardN2
    May 24, 2023 at 20:08
  • 1
    $\begingroup$ @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. $\endgroup$
    – Michael E2
    May 24, 2023 at 20:37
  • $\begingroup$ Thank you very much! $\endgroup$
    – HadamardN2
    May 25, 2023 at 2:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.