Re: Re(2): RE: ExpIntegralEi
- To: mathgroup at smc.vnet.net
- Subject: [mg18750] Re: [mg18681] Re(2): [mg18491] RE: [mg18463] ExpIntegralEi
- From: Carl Woll <carlw at u.washington.edu>
- Date: Sat, 17 Jul 1999 02:36:59 -0400
- Organization: Physics Department, U of Washington
- References: <199907150546.BAA15932@smc.vnet.net.>
- Sender: owner-wri-mathgroup at wolfram.com
Andrzej,
I thought I would try to shed more light on what's going on here. Rather than
analyzing the somewhat complicated situation which initiated this thread, it's
easier to see what's going on with a simple example. Consider
Integrate[ 1/(t-1) , {t, -.1 I, I}]
As written, this integral is not well defined. When integrating over the real
line, we clearly want the integration contour to be a straight line. However,
when integrating over the complex plane, we need to define the contour of
integration. By default NIntegrate chooses a straight path, so let's make the
natural assumption that the integration contour is again straight, and then
revisit the issue.
The integrand has a simple pole at t=1, and the integral along the straight path
from -.1 I to I is well defined. On the other hand, the indefinite integral is
Log[t-1], and this function has a branch cut. Mathematica by default chooses the
branch cut of the Log function to run along the negative imaginary axis. With
this choice of branch cut, the integration contour passes through the branch
cut. Since Integrate works by first evaluating the indefinite integral, and then
plugging in the limits of integration, when Mathematica evaluates the integral
using Integrate it gets Log[t-I]-Log[t+.1 I]. When evaluated this answer is off
by 2 Pi I. This is of course directly related to the residue of the integrand at
the simple pole t=1, and would be correct if the contour of integration wasn't
straight, but wrapped around t=1.
Carl Woll
Dept of Physics
U of Washington
Andrzej Kozlowski wrote:
> I have had an interesting exchange with Hendrik van Hees concerning this
> problem (one of his messages to me is included below). I think he has
> produced the clearest explanation of the problem involved in symbolic
> integration in this case. However, I do not agree with his suggestion
> that the "symbolic answer" is as correct or incorrect as the numerical
> one. In my opinion, there is no ambiguity here: the numerical answer is
> right. As this seems to raise faily interesting issues relating to both
> mathematics and Mathematica I would like to try to explain this a little
> more precisely. ( However, I have to first say that I am treading on thin
> ice as this area is not my speciality and my knowledge does not go beyond
> what would be expected of any professional mathematician in an area in
> which he has never worked i.e. almost nil).
>
> An integral like the one considered here is simply an integral of a
> complex valued function over a measurable subset of R^2. This is just a
> special case of a more general concept of an integral of a function with
> values in a Banach space (called a Bochner integral) and evaluating it
> simply amounts to evaluating its real and imaginary parts separately and
> taking the resulting complex number. All that is required for such an
> integral to be well defined is that the real and imaginary parts of the
> function be integrable real functions. Such an integral can always quite
> anambigiously be evaluated simply by evaluating its real and imaginary
> parts separately. Each of these can be reduced two repeated comptation of
> ordinary integrals by Fubini's theorem:
>
> Thus with
>
> In[1]:=
> c2 = -0.05018627683354541 - 0.153047656745338 I;
> In[2]:=
> c3 = -0.7828709924214918 + 0.2780791279205129 I;
> In[3]:=
> c5 = -0.6758555487562639 - 0.04753624179417532 I;
>
> In[5]:=
> f[x_, y_] := Exp[beta*c2 + s*(c3 + beta*c5)]
> In[6]:=
> f1[x_, y_] := Re[f[x, y]]
> In[7]:=
> f2[x_, y_] := Im[f[x, y]]
>
> We get:
>
> In[8]:=
> NIntegrate[f[x, y], {s, 0, 1}, {beta, 0, 1}]
> Out[8]=
> 0.587252 + 0.0191685 I
>
> In[9]:=
> NIntegrate[f2[x, y], {s, 0, 1}, {beta, 0, 1}]
> Out[9]=
> 0.0191685
>
> In[8]:=
> NIntegrate[f1[x, y], {s, 0, 1}, {beta, 0, 1}]
> Out[8]=
> 0.587252
>
> All these computations are instantenous and there are no ambiguities or
> singularities involved, since the functions are continuous and thus
> integrable. There is no need for complex analysis here, both the
> integrals of f1 and f2 can be computed as ordinary Lebesgue integrals on R^2.
>
> However, ths situation is quite different when we try to use symbolic
> integration. Mathematica is unable to calclulate Integrate[f1[x, y], {s,
> 0, 1}, {beta, 0, 1}] or Integrate[f2[x, y], {s, 0, 1}, {beta, 0, 1}]:
>
> In[10]:=
> Integrate[f1[x, y], {s, 0, 1}, {beta, 0, 1}]
> Out[10]=
> Integrate[Re[Power[E, (-0.0501863 - 0.153048 I) beta +
>
> (-0.782871 + 0.278079 I -
>
> (0.675856 + 0.0475362 I) beta) s]], {s, 0, 1},
>
> {beta, 0, 1}]
>
> Instead of viewing the integral as essentially two independent real
> integrals Mathematica considers Exp as an analytic function in the
> complex plane and tries to express the integral in terms of special
> functions which are themselves defined as integrals in the complex plane.
> However, as these function are not themselves analytic but have branch
> cuts this introduces ambiguities in the answer. These ambiguities however
> are not genuine, they result merely from the need to express the answers
> in terms of already defined functions. One can see what happens by
> considering the Henrik's computations below:
>
> In[19]:=
> Clear[c2, c3, c5]
>
> In[20]:=
> Integrate[Exp[beta*c2 + s*(c3 + beta*c5)], {s, 0, 1}] // Simplify
> Out[20]=
> beta c2 c3 + beta c5
> E (-1 + E )
> -----------------------------
> c3 + beta c5
>
> This is function has a simple pole at c3=-beta*c5 . Numerical integration
> of such a function shoudl noto cause any problem as long as c2 and c5 are
> not both zero, but symbolic integration leads to new complications:
>
> In[21]:=
> Integrate[(E^(beta*c2)*(-1 + E^(c3 + beta*c5)))/
> (c3 + beta*c5), {beta, 0, 1}]
> Out[21]=
> c2 c3 c2 c3 + c3 c5
> ExpIntegralEi[-----] - ExpIntegralEi[-------------]
> c5 c5
> --------------------------------------------------- -
> (c2 c3)/c5
> E c5
>
> 1 c2 c3
> -- (ExpIntegralEi[c2 + -----] -
> c5 c5
>
> c2 c3 + c3 c5 (c2 c3)/c5
> ExpIntegralEi[c2 + c5 + -------------]) / E
> c5
> Now we have an answer, but in terms of a function which has branch cuts.
> This sort of answer is inherently ambiguous. With an unlucky choice of
> the parameters it will produce numerically incorrect answer. In this I
> agree with Henrik. However, this does no tchange the fact that the answer
> one gets from symbolic integration is just simply wrong! The original
> integral definitely is well defined and has a unique correct value.
>
>
> On Tue, Jul 13, 1999, Hendrik van Hees <h.vanhees at gsi.de> wrote:
>
> >Hello again,
> >
> >here is Mathematica's result:
> >
> >In[2]:= Integrate[Exp[beta*c2+s*(c3+beta*c5)],{s,0,1},{beta,0,1}]
> >
> > c2 c3 c2 c3
> > ExpIntegralEi[c2 + -----] - ExpIntegralEi[-----]
> > c5 c5
> >Out[2]= -(------------------------------------------------) -
> > (c2 c3)/c5
> > c5 E
> >
> > c3 (c2 + c5) c2 c3
> > ExpIntegralEi[------------] - ExpIntegralEi[c2 + c3 + ----- + c5]
> > c5 c5
> >> -----------------------------------------------------------------
> > (c2 c3)/c5
> > c5 E
> >
> >In[3]:= Simplify[%]
> >
> > c2 c3 c2 c3
> >Out[3]= -((ExpIntegralEi[c3 + -----] - ExpIntegralEi[-----] +
> > c5 c5
> >
> > c2 (c3 + c5) c2 c3
> >> ExpIntegralEi[------------] - ExpIntegralEi[c2 + c3 + ----- +
> c5]) /
> > c5 c5
> >
> > (c2 c3)/c5
> >> (c5 E ))
> >
> >This function is not unique because ExpIntegralEi has a branch cut along
> >the negative real axis. Thus it is an analytic function in the right
> >half plane of its argument. Thus if the combinations of the constants
> >c2,c3,c5 are on the other (left) half plane its value is not unique. One
> >must specify on which sheet of the Riemann surface you want to evaluate
> >your integral.
> >
> >The reason is simply seen. If you take the inner integral for its own it
> >results:
> >
> >In[4]:= Integrate[Exp[beta*c2+s*(c3+beta*c5)],{s,0,1}]//Simplify
> >
> > beta c2 c3 + beta c5
> > E (-1 + E )
> >Out[4]= -----------------------------
> > c3 + beta c5
> >
> >Thus if along the path of integration for beta (beta goes along the real
> >axis from 0 to 1 in the outer integral) the denominator becomes 0 there
> >is a ambiguity and this is the reason for the above mentioned branch
> >cut. Thus the values you get out from NIntegrate could eventually depend
> >on the gridding you take for integrating. I think one can check this
> >explicitly with the given definite values of the parameters. But these I
> >have to look up first.
> >
> >Greetings,
> >
> >--
> >Hendrik van Hees Phone: ++49 6159 71-2755
> >c/o GSI-Darmstadt SB3 3.162 Fax: ++49 6159 71-2990
> >Planckstr. 1 mailto:h.vanhees at gsi.de
> >D-64291 Darmstadt http://theory.gsi.de/~vanhees/vanhees.html
> >
>
> Andrzej Kozlowski
> Toyama International University
> JAPAN
> http://sigma.tuins.ac.jp/
> http://eri2.tuins.ac.jp/
--
Carl Woll
Dept of Physics
U of Washington
- References:
- Re(2): RE: ExpIntegralEi
- From: Andrzej Kozlowski <andrzej@tuins.ac.jp>
- Re(2): RE: ExpIntegralEi