Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
expint gives unexpected result

Subject: expint gives unexpected result

From: Peter Mairhofer

Date: 14 Aug, 2014 22:17:01

Message: 1 of 6

I compare expint() with numerical calculation.

The real part matches reasonably well but the imaginary part has a -2pi
jump *after* zero.

At least for real values I would not expect this.

Furthermore, [1] states that Ei is ambiguous but E1 is not.

Where does this discrepancy come from and how to fix it?

Thanks
Peter

PS: Here is a small piece of code that compares the function with real
arguments:

Ts = 1e-4;
t = [ -10:Ts:-Ts Ts:Ts:1 ];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');



[1] http://en.wikipedia.org/wiki/Exponential_integral

Subject: expint gives unexpected result

From: Roger Stafford

Date: 15 Aug, 2014 05:45:11

Message: 2 of 6

Peter Mairhofer <63832452@gmx.net> wrote in message <lsjcgs$jh1$1@news.albasani.net>...
> I compare expint() with numerical calculation.
> .......
> Ts = 1e-4;
> t = [ -10:Ts:-Ts Ts:Ts:1 ];t=t(:);
> .... plot(t, real([Ts*cumtrapz(exp(-t)./t) , (expint(t(1))-expint(t))]));
- - - - - - - - - - - -
  In your cumtrapz integration you have integrated right across a singularity in the integrand. This will give you the wrong answer. To avoid the singularity you need to follow a semi-circular path in the complex plane around the singularity and that will produce a very different result. It has to do with the theory of residues of complex functions, and this integrand does indeed have a non-zero residue at z = 0. That is why the expint function has a branch cut along the negative real axis.

Roger Stafford

Subject: expint gives unexpected result

From: Peter Mairhofer

Date: 15 Aug, 2014 16:47:06

Message: 3 of 6

Hi Roger,

On 2014-08-14 22:45, Roger Stafford wrote:
> Peter Mairhofer <63832452@gmx.net> wrote in message
> <lsjcgs$jh1$1@news.albasani.net>...
>> I compare expint() with numerical calculation.
>> .......
>> Ts = 1e-4;
>> t = [ -10:Ts:-Ts Ts:Ts:1 ];t=t(:);
>> .... plot(t, real([Ts*cumtrapz(exp(-t)./t) , (expint(t(1))-expint(t))]));
> - - - - - - - - - - - -
> In your cumtrapz integration you have integrated right across a
> singularity in the integrand. This will give you the wrong answer.
> To
> avoid the singularity you need to follow a semi-circular path in the
> complex plane around the singularity and that will produce a very
> different result. It has to do with the theory of residues of complex
> functions, and this integrand does indeed have a non-zero residue at z =
> 0.

Thanks. But are you sure about this? I read somewhere on the internet
that this singularity is a removeable singularity and hence I can
integrate over it.

Furthermore, this happens even when everything is real (as with the t
vector above). So I think complex analysis theory should not even be
required but Cauchy Principal value instead. In this case contributions
from either side cancel as they should.

Finally, this integral is part of a well defined problem where I know
the solution. And the result of the *numerical* integration is *correct*
and not the other way round.

So I think I do want to integrate over the singularity.

Can it be that it is the other way round and expint purposely does not
integrate over the singularity?

In this case, the documentation on expint is confusing (at least for me)
because only the integral definition is given. A remark that expint does
not integrate over the singularity would be helpful.

> That is why the expint function has a branch cut along the negative
> real axis.

Unfortunately my complex analysis course is over 10 years ago and I
vagely remember it. Why do we need a branch cut here? Isn't a branch cut
only used for a closed path integration? For example this

http://i.stack.imgur.com/8VNw5.png

would be a brunch cut along the negative imaginary axis?

If not, would you mind telling me over which z E1(z) as defined by
MATLAB integrates for

z = 1, ..., inf
z = -10, ..., inf
z = 1j, ..., j*inf
z = -1j, ..., j*inf


Thanks so much!
Peter

Subject: expint gives unexpected result

From: Nasser M. Abbasi

Date: 15 Aug, 2014 19:01:07

Message: 4 of 6

On 8/15/2014 11:47 AM, Peter Mairhofer wrote:

>
> Thanks. But are you sure about this? I read somewhere on the internet
> that this singularity is a removeable singularity and hence I can
> integrate over it.
>

Not following the discussion, just jumping in :)

In analytical integration, some systems provide an option
to bypass removable singularity in integration. For example,

Integrate[1/x, {x, -1, 1}] does not work, since there is
a pole at x=0. But using Cauchy definition one can integrate it

Integrate[1/x, {x, -1, 1}, PrincipalValue -> True]
(* 0 *)

?PrincipalValue

"PrincipalValue is an option for Integrate that specifies whether
the Cauchy principal value should be found for a definite
integral."

You can read more about Cauchy principal value for integration here

http://en.wikipedia.org/wiki/Cauchy_principal_value

(do not know if this applies to what you are doing, since
I did not follow the previous discussion you had, and do not
know if Matlab int supports this option either)

--Nasser

Subject: expint gives unexpected result

From: Roger Stafford

Date: 15 Aug, 2014 23:47:11

Message: 5 of 6

Peter Mairhofer <63832452@gmx.net> wrote in message <lsldia$bcv$1@news.albasani.net>...
> Thanks. But are you sure about this? I read somewhere on the internet
> that this singularity is a removeable singularity and hence I can
> integrate over it.
>
> Furthermore, this happens even when everything is real (as with the t
> vector above). So I think complex analysis theory should not even be
> required but Cauchy Principal value instead. In this case contributions
> from either side cancel as they should.
>
> Finally, this integral is part of a well defined problem where I know
> the solution. And the result of the *numerical* integration is *correct*
> and not the other way round.
>
> So I think I do want to integrate over the singularity.
>
> Can it be that it is the other way round and expint purposely does not
> integrate over the singularity?
>
> In this case, the documentation on expint is confusing (at least for me)
> because only the integral definition is given. A remark that expint does
> not integrate over the singularity would be helpful.
- - - - - - - - - -
  Peter, the discrepancy between the expint result and cumtrapz result using the Cauchy principal value concept is due to the difference between the Ei and E1 functions. The matlab expint function computes the E1 function whose analytic continuation defines it over the entire complex plane. (In fact the theoretical continuation has infinitely many branches but expint gives only the principal branch - hence the branch cut on the negative real axis.) The Ei function is really only valid along the real axis due to the singularity at x = 0. It cannot be analytically continued into the complex plane.

   However, if you take only the real part of the output from expint, you should still get the correct answer for your integration. That is, the integral from a to b for exp(-t)/t for real a and b is given by

 real(expint(a)-exp(b))

even of a is negative and b is positive, which would represent integration across the singularity (courtesy of Cauchy.)

  Any differences you observed in the real parts were undoubtedly due to the inaccuracies of cumtrapz, combined with the fact that you were comparatively far away from either side of zero at -1e-4 and +1e-4. Even if you use one of the matlab 'quad' functions on two intervals much more closely-spaced on each side of zero, you can expect larger than usual errors due to the large values of 1/t encountered, which makes integration difficult.

  The Mathworks documentation for expint does state that for real positive x

 E1(-x) = -Ei(x) - pi*i

and from this you can deduce that real(expint(a)-exp(b)) will give you what you are seeking without resorting to an integration procedure.

Roger Stafford

Subject: expint gives unexpected result

From: Roger Stafford

Date: 16 Aug, 2014 04:33:08

Message: 6 of 6

"Roger Stafford" wrote in message <lsm65v$ra3$1@newscl01ah.mathworks.com>...
> real(expint(a)-exp(b))
> ......
> and from this you can deduce that real(expint(a)-exp(b)) will give you what you are
> ......
- - - - - - -
  I had typographical errors in two places in my last article. In each place real(expint(a)-exp(b)) should be replaced by real(expint(a)-expint(b))

Roger Stafford

Tags for this Thread

No tags are associated with this thread.

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us