# Multi-dimensional integral in the complex plane with poles and essential singularity

I've passed the last week searching a way to numerically integrate this multi-dimensional integral in the complex plane at the poles and avoiding the singularity at z=0: $$\oint_{C}\oint_{C\ auound\ the\ poles\ only}\frac{z_{1}z_{2}}{z_{2}^{-1}p+(1-p)z_{1}-1}e^{\frac{1}{z_{1}}+z_{1}+\frac{1}{z_{2}}+z_{2}} \, \mathrm{d}z_{1}\mathrm{d}z_{2}$$ where C is a circle of radius 1 centred at 0 and $p\in]0,1[$ and I used this topic Finding residues of multi-dimensional complex functions to compute this integral.
The fact is that when I give the whole integral to Mathematica, I don't have the contribution of the poles because the contour cannot be deformed due to the singularity. So I'm trying to find a way to exclude the singularity of the contour, maybe with a small contour around z=0 using piecewise or something similar (see the figure here http://img834.imageshack.us/img834/6499/v2uq.png where red=singularity and green=poles) to get only the contribution of the poles.

Mathematica code

NContourIntegrate[f_, par : (z_ -> g_), {t_, a_, b_}] :=
NIntegrate[Evaluate[D[g, t] (f /. par) /. t -> t1], {t1, a, b}]

Clear[Pinfz1];
Pinfz1[p_?NumericQ, z2_?NumericQ] :=
NContourIntegrate[1/(2*I*Pi)^2*z1*z2/(p + (1 - p)*z1*z2-z2)*E^(1/z1 + z1 + 1/z2 + z2),
z1 -> Exp[I t],
{t, 0, 2*Pi}]

Clear[Pinfz];
Pinfz[p_?NumericQ] := NContourIntegrate[Pinfz1[p, z2], z2 -> Exp[I t], {t, 0, 2*Pi }]

-
Where does this computation arise? I'm wondering if a double contour integration is even what's actually wanted. – Daniel Lichtblau Jun 18 '13 at 16:03
This come from a physical situation (particle diffusion on lattice), so there is actually 2 integrals – Jprog Jun 18 '13 at 22:54

I'm not familiar with the physics so I cannot say whether integrating over this set of 2 real dimension, in C^2, is what is wanted. I think the code below will cover the product space of the contours that are requested.

ii = z1*z2/(p/z2 + (1 - p) z1 - 1)*Exp[1/z1 + z1 + 1/z2 + z2];
i1 = (ii /. {z1 -> Exp[I*t1], z2 -> Exp[I*t2]})*I*Exp[I*t1]*I*
Exp[I*t2];
f[pval_] :=
Quiet[NIntegrate[
Evaluate[i1 /. p -> pval], {t1, 0, 2*Pi}, {t2, 0, 2*Pi}]]
f[1/2]
f[1/3]

(* Out[114]= 89.3519154317 - 5.40968238685*10^-8 I

Out[115]= 62.3532282953 + 0.000124769719229 I *)


Also:

f[.1]

(* Out[116]= 34.4365223939 + 0.000414467508089 I *)

f[.9]

(* Out[117]= 170.924622109 + 0.00264992865898 I *)


Offhand I cannot say much about how accurate these might be. When I did them at higher working precision I got consistent results, but that does not prove anything.

-