Re(2): RE: ExpIntegralEi

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

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/

**Follow-Ups**:**Re: Re(2): RE: ExpIntegralEi***From:*Carl Woll <carlw@u.washington.edu>

**Re: Solving difficult integral**

**Re: [Q] Implementing identities as rules**

**V4.0 Compatibility (Financial Derivatives book)**

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