Analytic expression for a Complex Hilbert Transform in Mathematica

Question

I need to solve the following integral equations for a problem I'm working on -

$\displaystyle \frac{-i}{2 \pi}$ $\int_{-a}^{a} \mathrm{dt}\,\, \frac{e^{i k t}}{t + i \tau}$ and $\displaystyle \frac{-i}{2 \pi}$ $\int_{a}^{\infty} \mathrm{dt}\,\, \frac{e^{i k t}}{t + i \tau}$

where $\tau, k \in \mathbb{R}$.

Mathematica can do these numerically and I get sensible plots as a function of $\tau$ for given $a$ and $k$, where by sensible I mean that the integrals go to 0 as $\tau \to \pm \infty$.

However, I'd like to find an analytic expression for these, even if it is in terms of exponential integrals, etc. When I try to do this, Mathematica spits out an answer that agrees with the numerical integration at small $\tau$, but at large $\tau$ the analytic expressions break down.

For concreteness, suppose I take $\frac{-i}{2 \pi}$ $\int_{a}^{\infty} \mathrm{dt}\,\, \frac{e^{i k t}}{t + i \tau}$ with k = 1 and a = 1.

Doing this integral numerically and plotting the real part vs $\tau$,

Plot[Re[-I/(2 π)
    NIntegrate[E^(I  t)/(
    t + I τ), {t, 1, ∞}]], {τ, -200, 200}, 
 PlotRange -> All]

But now, suppose I do this integral analytically.

Integrate[-I/(2 π) E^(I t)/(t + I τ), {t, 1, ∞}]

This gives an output

    (E^τ (π + 2 I CosIntegral[1 + I τ] + 2 I SinhIntegral[I - τ]))/(4 π)

and when I plot the Real part of this, the result looks like -

I can't really understand why this breaks down and would like to have analytic expressions for both those integrals for an values of $\tau, k$ and $a$.

Answer 1

This is due to the sum of two very large numbers (coming from CosIntegral and SinhIntegral) being carried out without sufficient machine precision used to represent them. You can fix it giving an appropriate value of WorkingPrecision as an option to plot.

You can see quite clearly that the problem comes from this by plotting the two functions (the one coming from the numerical integration and the other from the numerical evaluation of the symbolic integration) at varying values of WorkingPrecision. The lesser WorkingPrecision is the sooner the problem arises:

GraphicsRow@Table[
   Plot[
    {
     Re[E^\[Tau] (\[Pi] + 2 I CosIntegral[1 + I \[Tau]] + 
         2 I SinhIntegral[I - \[Tau]])]/(4 \[Pi]),
     Re[-I/(2 \[Pi]) NIntegrate[
        E^(I t)/(t + I \[Tau]), {t, 1, \[Infinity]}]]
     },
    {\[Tau], 0, 100},
    PlotRange -> All,
    WorkingPrecision -> wp,
    ImageSize -> Medium
    ],
   {wp, {20, 40, 80}}
   ]~Monitor~{wp}

EDIT1: How to obtain a specific value with wanted precision

To obtain a correct particular value you can try enforcing the required precision of the result with the second argument of N. In this case even small values of this parameter seem to work, probably because Mathematica automatically uses the required precision during the computation to obtain the requested one at the end:

f[\[Tau]_] := 
 Re[E^\[Tau] (\[Pi] + 2 I CosIntegral[1 + I \[Tau]] + 
     2 I SinhIntegral[I - \[Tau]])]/(4 \[Pi])
N[f[100]]
(* Out=6.72029*10^42 *)
N[f[100], 2]
(* Out=0.0013 *)

Indeed, you can use this approach to obtain a more reliable plot. It may be necessary to modify the value of $MaxExtraPrecision to avoid the 50 digits internal default limit of Mathematica:

f[\[Tau]_] := 
 Re[E^\[Tau] (\[Pi] + 2 I CosIntegral[1 + I \[Tau]] + 
     2 I SinhIntegral[I - \[Tau]])]/(4 \[Pi])
Block[{$MaxExtraPrecision = 1000},
 Show[
  ListLinePlot[
   Table[
    {\[Tau], N[f[\[Tau]], 4]},
    {\[Tau], 0, 200}
    ],
   PlotStyle -> {Thick, Green},
   PlotRange -> All
   ],
  Plot[
   Re[-I/(2 \[Pi]) NIntegrate[
      E^(I t)/(t + I \[Tau]), {t, 1, \[Infinity]}]],
   {\[Tau], 0, 200},
   PlotRange -> All,
   ImageSize -> Large,
   PlotStyle -> {Red, Dashed}
   ]
  ]
 ]

EDIT 2: How to reliably evaluate the function at non-integer values of $k$

If more in general we are interested in evaluating the function at other values of $k$ other than 1, we must be careful in the way we give the value of $k$. If we use a value like k=-1.2 in the definition of $f$, Mathematica will not be able to use its arbitrary precision engine, as well explained for example in this answer. A way around this is to Rationalize the value of $k$ before the evaluation of $N$. Here is a working example of the correct evaluation of $f(\tau,k=-1.2)$ for $\tau=1,...,40$:

f[\[Tau]_, k_] := 
 f[\[Tau], k] = 
  Re@Integrate[-I/(2 \[Pi]) E^(I k t)/(t + I \[Tau]), {t, 
     1, \[Infinity]}, 
    Assumptions -> {\[Tau] \[Element] Reals, k \[Element] Reals}]
ListLinePlot[
  Table[
   {\[Tau], N[f[\[Tau], Rationalize[-1.2]], {Infinity, 4}]},
   {\[Tau], 0, 40}
   ],
  PlotRange -> All,
  ImageSize -> Large
  ]~Monitor~{\[Tau]}

Note that this may take a while to evaluate because the larger $\tau$ gets the more digits Mathematica has to use to correctly compute the result. The numerical integration is definitely faster in these cases.

You can also use this approach with Plot, if you are very patient. A way to get the feel of the hardness of such a computation without having to wait a very long time to get the complete result is for example using the dynamicPlot function given in this answer, which allows to see the plot while it's being drawn. To do this evaluate the function dynamicPlot in the linked answer and then use the following code:

f[\[Tau]_, k_] := 
 f[\[Tau], k] = 
  Re@Integrate[-I/(2 \[Pi]) E^(I k t)/(t + I \[Tau]), {t, 
     1, \[Infinity]}, 
    Assumptions -> {\[Tau] \[Element] Reals, k \[Element] Reals}]
g[\[Tau]_] := N[f[\[Tau], Rationalize[-1.2]], {Infinity, 4}]
dynamicPlot[
 g,
 {x, 0, 40},
 PlotRange -> All,
 ImageSize -> Large,
 WorkingPrecision -> Infinity
 ]

Answer 2

So upon thinking about this more and playing with different assumptions when trying to integrate the above function(s), what I found was the following -

$\int_{a}^{\infty} dt \frac{e^{i k t}}{t + i \tau} = e^{k\tau}\,\Gamma\left(k(-ia + \tau) \right)$, where $\Gamma$ is the incomplete Gamma function.

Now, suppose I try doing this integral with Mathematica, with the assumption that $a>0$. I get

E^(k \[Tau]) (Gamma[0, k (-I a + \[Tau])] - Log[-I k] - Log[a + I \[Tau]] + Log[k (-I a + \[Tau])])

Now for large values of $\tau$, if you want to plot this function, you have to follow @glance's answer given above. On examining the answer, however, you see that the Log part of this answer is 0, leaving one with just the $\Gamma$-function which Mathematica has no problems evaluating for any $\tau$ or $k$!

The integral $\int_{-a}^a dt \frac{e^{ikt}}{t + i \tau}$ is a little more complicated, but in Mathematica, it evaluates to (with assumptions $\tau \in \text{Reals}, a >0$)

E^(k \[Tau]) (-Gamma[0, k (-I a + \[Tau])] + Gamma[0, k (I a + \[Tau])] - Log[a - I \[Tau]] + Log[a + I \[Tau]] + Log[-I \[Tau]] - Log[I \[Tau]] - Log[k (-I a + \[Tau])] + Log[k (I a + \[Tau])])

Again, looking at the Log part of this answer, one finds that it is

0 if $k \tau > 0$ ; $2 \pi i$ if $k>0$ and $\tau < 0$ ; $-2 \pi i$ if $k < 0$ and $\tau > 0$.

Then, this function can be built with the $\Gamma$-functions and a piece-wise part, which matches the numerically integrated function and is evaluated by Mathematica for any $\tau$ and $k$ without any issues whatsoever.