MathGroup Archive 1999

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re(2): Re(2): RE: ExpIntegralEi

  • To: mathgroup at smc.vnet.net
  • Subject: [mg18755] Re(2): [mg18681] Re(2): [mg18491] RE: [mg18463] ExpIntegralEi
  • From: Andrzej Kozlowski <andrzej at tuins.ac.jp>
  • Date: Tue, 20 Jul 1999 01:33:21 -0400
  • Sender: owner-wri-mathgroup at wolfram.com

Actually, the original example was omewhat different since the
singularity that appears in 

 beta c2        c3 + beta c5
E        (-1 + E            )
 -----------------------------
         c3 + beta c5

is removable and not, as I incorrectly wrote, a simple pole as in your
example. Nevertheless you are  right: the source of the problem is the
same. It could be avoided if Mathematica used a less simple minded
algorithm when integrating in the complex plane analytic functions with
branch cuts. It seems that it should always be possible to avoid crossing
the branch-cut, for example by using a path consisting of a vertical and
horizontal line segments going around it. I am not an expert but my
impression is that this should not be too difficult to implement so
perhaps we might get it in the next version?

Andrzej

On Fri, Jul 16, 1999, Carl Woll <carlw at u.washington.edu> wrote:

>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
>
>
>
>
>


Andrzej Kozlowski
Toyama International University
JAPAN
http://sigma.tuins.ac.jp/
http://eri2.tuins.ac.jp/



  • Prev by Date: Compile problem: Mathematica 4
  • Next by Date: RE: Re: Rationalizing the denominator (better solution!)
  • Previous by thread: Re: Compile problem: Mathematica 4
  • Next by thread: Limit of an Integer Function