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
]