Complex exponential integral: Mathematica and MATLAB give unexpected results

1.2k Views Asked by At

I currently compare analytical vs. numerical evaluation of the complex exponential integral and find mismatches: The imaginary part differs by $\pm \pi$ and the real part has a large error when integration does not include the singularity.

While Wikipedia and other sources like http://www.cs.utah.edu/~vpegorar/research/2011_JGT/paper.pdf state that $E_i(z)$ might be ambiguous, $E_1(z)$ should not be.

I use $E_1(z)$ as implemented as expint in MATLAB and ExpIntegralE in Mathematica - both results match. The definition is as follows:

$$ E_1(z) = \int_z^{\infty} \frac{e^{-v}}{v} \, dv $$

I then compare the two functions when $z \in \mathbb{R}$:

Ts = 1e-4;
t=-3:Ts:3;t=t(:);
figure, plot(t, real([Ts*cumtrapz(exp(-t)./t) , (expint(t(1))-expint(t))]));legend('cumtrapz','E_1');
figure, plot(t, imag([Ts*cumtrapz(exp(-t)./t) , (expint(t(1))-expint(t))]));legend('cumtrapz','E_1');

While the real part matches, the imaginary part has the phase jump of $-\pi$:

comparison of real part comparison of imaginary part

Similarly, I can do this when $z$ is imaginary:

tau = [ -10:Ts:-Ts Ts:Ts:7 ];
w = 1i*pi*tau;w = w(:);
figure; plot(imag(w), real([cumtrapz(exp(-w)./w)*Ts*1i*pi, (expint(w(1))-expint(w))]));
figure; plot(imag(w), imag([cumtrapz(exp(-w)./w)*Ts*1i*pi, (expint(w(1))-expint(w))]));

The phase shift is $+\pi$ (I only show the imaginary part):

comparison of imaginary part

Furthermore, I observe a slight divergence in the real part when I do not include the singularity (i.e., when I start integration at or above zero). This is probably because the contributions cancel when I include the the singularity but it bothers me that the difference is so huge! And what could be done about it if I just want to start the integration from zero?

tau = [ -n:Ts:0-Ts Ts:Ts:N-n-Ts ];
w = 1i*pi*gamma*tau;w = w(:);
figure; plot(imag(w), real([cumtrapz(exp(-w)./w)*Ts*1i*pi, (expint(w(1))-expint(w))]));
figure; plot(imag(w), imag([cumtrapz(exp(-w)./w)*Ts*1i*pi, (expint(w(1))-expint(w))]));

I only show the real part:

comparison of real part