Re: Complex path integral wrong
- To: mathgroup at smc.vnet.net
- Subject: [mg131373] Re: Complex path integral wrong
- From: "Dr. Wolfgang Hintze" <weh at snafu.de>
- Date: Tue, 2 Jul 2013 00:45:02 -0400 (EDT)
- Delivered-to: l-mathgroup@mail-archive0.wolfram.com
- Delivered-to: l-mathgroup@wolfram.com
- Delivered-to: mathgroup-outx@smc.vnet.net
- Delivered-to: mathgroup-newsendx@smc.vnet.net
- References: <kqomj6$p76$1@smc.vnet.net> <kqrj4e$4bn$1@smc.vnet.net>
Am Montag, 1. Juli 2013 11:45:50 UTC+2 schrieb Kevin J. McCann:
> A followup on my earlier reply.
>
>
>
> If I compare the analytical and numerical results from Mathematica, the
>
> discrepancy occurs on the first leg of the integration
>
> (1+I -> 1+I R). It is clear that the numerical integration is correct.
>
>
>
> Could this be some kind of branch cut issue with the hypergeometric
>
> function evaluation?
>
>
>
> Kevin
>
>
>
> On 6/30/2013 3:26 AM, Dr. Wolfgang Hintze wrote:
>
> > I suspect this is a bug
>
> > In[361]:= $Version
>
> > Out[361]= "8.0 for Microsoft Windows (64-bit) (October 7, 2011)"
>
> >
>
> > The follwing path integral comes out wrong:
>
> >
>
> > R = 3 \[Pi] ;
>
> > Integrate[Exp[I s]/(
>
> > Exp[s] - 1 ), {s, 1 + I, 1 + I R, -1 + I R, -1 + I, 1 + I}] // FullSimplify
>
> >
>
> > Out[351]= 0
>
> >
>
> > It should have the value
>
> >
>
> > In[356]:= (2 \[Pi] I) Residue[Exp[I s]/(Exp[s] - 1 ), {s, 2 \[Pi] I}]
>
> >
>
> > Out[356]= (2 \[Pi] I) E^(-2 \[Pi])
>
> >
>
> > Without applying FullSimplify the result of the integration is
>
> >
>
> > In[357]:= R = 3*Pi;
>
> > Integrate[
>
> > Exp[I*s]/(Exp[s] - 1), {s, 1 + I, 1 + I*R, -1 + I*R, -1 + I, 1 + I}]
>
> >
>
> > Out[358]=
>
> > I*E^((-1 - I) - 3*Pi)*((-E)*Hypergeometric2F1[I, 1, 1 + I, -(1/E)] +
>
> > E^(3*Pi)*Hypergeometric2F1[I, 1, 1 + I, E^(-1 + I)]) +
>
> > I*E^(-I - 3*Pi)*(Hypergeometric2F1[I, 1, 1 + I, -(1/E)] -
>
> > E^(2*I)*Hypergeometric2F1[I, 1, 1 + I, -E]) +
>
> > I*E^I*(Hypergeometric2F1[I, 1, 1 + I, -E]/E^(3*Pi) -
>
> > Hypergeometric2F1[I, 1, 1 + I, E^(1 + I)]/E) +
>
> > I*E^(-1 - I)*(-Hypergeometric2F1[I, 1, 1 + I, E^(-1 + I)] +
>
> > E^(2*I)*Hypergeometric2F1[I, 1, 1 + I, E^(1 + I)])
>
> >
>
> > which, numerically, is
>
> >
>
> > In[359]:= N[%]
>
> >
>
> > Out[359]= -2.7755575615628914*^-17 + 2.7755575615628914*^-17*I
>
> >
>
> > i.e. zero.
>
> >
>
> > On simpler functions like 1, s and s^2 (instead of Exp[I s]) it works out fine, but not so with e.g. Sin[s] in which case we get 0 again (instead of Sinh[2 \[Pi]]).
>
> >
>
> > The integration topic seems to be full of pitfalls in Mathematica...
>
> >
>
> > Best regards,
>
> > Wolfgang
>
> >
Kevin,
I think you are right, the branch cut might be the reason.
Indeed there seem to be some deeper difficulties with the hypergeometric function 2F1(a,b,c;z) (and/or with branch cuts) in Mathematica.
Because this function appears frequently as a result of indefinite integration which Mathematica transforms into definite integrals by taking limits, it deserves a somewhat closer look.
The function 2F1[a,b,c,z] as a function of z (if it doesn't reduce to a polynomial in z) has a branch cut in the complex z-plane ranging from +1 to +oo.
Let us investigate the behaviour on both side of the cut with an example.
We have
In[280]:= i0 = Simplify[TrigToExp[Integrate[Sin[s]/(Exp[s] - 1), s]]]
Out[280]= 1/2 E^(-I s) (Hypergeometric2F1[-I, 1, 1 - I, E^s] +
E^(2 I s) Hypergeometric2F1[I, 1, 1 + I, E^s])
Now, let's take the first part only
In[296]:= f = Hypergeometric2F1[-I, 1, 1 - I, E^s];
We now compute the values of f on both side of the cut with different methods
1) replacement with a numerical values close to -0 and +0
In[298]:= f /. {s -> 2 + I*y} /. y -> -0.001
Out[298]= -0.07641223288902461 - 0.0652022010883159*I
In[299]:= f /. {s -> 2 + I*y} /. y -> Plus[0.001]
Out[299]= -2.6886865517869416 + 5.642501876662279*I
2) Limit to the same numerical values as in 1)
This leads to the same results as in 1)
In[300]:= Limit[f /. {s -> 2 + I*y}, y -> -0.001]
Out[300]= -0.07641223288902478 - 0.06520220108831606*I
In[301]:= Limit[f /. {s -> 2 + I*y}, y -> Plus[0.001]]
Out[301]= -2.688686551786946 + 5.642501876662289*I
3) Limit to exactly y=0 making the distinction using Direction
In[305]:= Limit[f /. {s -> 2 + I*y}, y -> 0, Direction -> -1];
N[%]
Out[306]= -0.07649229821994431 - 0.06513712192483992*I
In[307]:= Limit[f /. {s -> 2 + I*y}, y -> 0, Direction -> Plus[1]];
N[%]
Out[308]= -0.07649229821994431 - 0.06513712192483992*I
Obviously, the result here is independent of Direction, and it is (approximately) equal to the approximation from the negative side (y<0).
I would consider this indifference with respect to Direction a bug.
And this might well be a cause of some of the difficulties in Integrate.
I'm not sure because I don't know which limiting procedure is used by Mathematica in order to create definite integrals out of indefinite ones.
Maybe others who know more can clarify the situation.
PS: By the way, you can easily check that Mathematica behaves correctly with other functions with a similar branch cut, e.g. the function g[z_] := Sqrt[1 - z] at z=2+Iy, so it seems that there is a particular difficulty with the hypergeometric function. Which is especially unpleasant in view of the importance of this type of function in Mathematica's integration algorithm.
Best regards,
Wolfgang