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

I've spent the last week searching for a way to numerically calculate the pole contributions to this multi-dimensional complex contour integral while avoiding the singularity at z=0: $$\oint_{C_1(0)}\oint_{C_p}\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}} \quad \text d z_1 \, \text d z_2$$ where $p \in (0,1)$, $C_1(0)$ is a circle of radius 1 centered at 0 and $C_p$ is a contour surrounding the poles only.

I used Finding residues of multi-dimensional complex functions to compute this integral. 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 to get only the contribution of the poles.

## 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= 89.3519154317 - 5.40968238685*10^-8 I

Out= 62.3532282953 + 0.000124769719229 I *)


Also:

f[.1]

(* Out= 34.4365223939 + 0.000414467508089 I *)

f[.9]

(* Out= 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.