# Numerical contour integrations in the complex plane - contour deformation gives different answer for analytic kernel

I am trying to do a contour integration in Mathematica numerically. In particular, I'm checking the identity:

$$H_m^{(1)}(z) =\frac{i^{-m}}{\pi}\int_{-\pi/2 + i \infty}^{\pi/2 - i \infty} \exp[i m \beta + i z \cos\beta]\mathrm d\beta$$

In Stratton's Electromagnetic Theory (I am a physicist), he claims (p. 368) that you can freely deform the contour since the kernel is an analytic function, assuming the contour stays within a certain region. Here's a diagram (from here):

This is all fine, and I have no problem with this; however, when I try to evaluate the contour integration in Mathematica with code that looks like:

Hmz[m_, z_] := Exp[-I m Pi/2]/Pi NIntegrate[Exp[I m B + I z Cos[B]],
{B, -Pi/2 + I Infinity, -Pi/2 + I Pi/2, Pi/2 - I Pi/2, Pi/2 - I Infinity}]


or alternatively (along a different, but apparently valid, contour):

Hmzother[m_, z_] :=  Exp[-I m Pi/2]/Pi NIntegrate[Exp[I m B + I z Cos[B]],
{B,  -Pi/2 + I Infinity, -Pi/2, Pi/2, Pi/2 - I Infinity}]


I get different answers - at least for the imaginary part (particularly for small $z$). For example, evaluating for $m = 3$ and $z = 1$ gives (not actual lines of code):

Hmz[3, 1] == 0.0195634 - 2.72402 I
Hmzother[3, 1] == 0.0195634 - 0.248805 I


For reference, the true value should be HankelH1[3,1] == 0.0195634 - 5.82152 I. The only difference between the definition of Hmz and Hmzother here are the contours, yet they should give the same answer. Ideally, I'd like to use the contour {B, I Infinity, 0, Pi, -I Infinity}, but this gives an even worse disagreement.

Does anybody know how to get Mathematica to do the contour integral correctly in the complex plane? Thanks all.

• I presume you're aware that HankelH1[] is built-in? Commented Jun 7, 2013 at 8:35
• I am indeed. This integral is a merely a stepping stone to an angular spectrum representation of cylindrical functions, as well as a test case to check I am doing things the right way. Commented Jun 7, 2013 at 8:39
• @Matthew The imaginary part of your integral strongly depends on PrecisionGoal and WorkingPrecision. You can play with them a bit. On the other hand you can get the both integrals the same assuming an appropriate Method, e.g. adding this option to your definitions Method -> "NewtonCotesRule" provides the same correct values. Commented Jun 7, 2013 at 8:48
• Personally, I prefer using double-exponential quadrature for situations like these. Here's an example, showing the use of a different contour: myHankelH1[m_?NumericQ, z_?NumericQ] := Exp[-I m Pi/2] NIntegrate[Exp[I m (ArcTan[u] - I u) + I z Cos[ArcTan[u] - I u]] (1/(1 + u^2) - I), {u, -5, 5}, Method -> "DoubleExponential"]/Pi. The result from myHankelH1[] compares favorably to the built-in function. (As a matter of fact, using Method -> "DoubleExponential" on your two functions seems to work nicely.) Commented Jun 7, 2013 at 9:07
• As I said, double-exponential quadrature works nicely even with your two example contours. If you're looking out for speed, then "GaussKronrodRule" and "ClenshawCurtisRule" work faster than "NewtonCotesRule". Commented Jun 7, 2013 at 9:31

(This is more of a plausibility argument than a rigorous answer.)

I had previously mentioned that one would want to avoid contours that go too near to the imaginary axis, since the resulting integrals will tend to be ill-conditioned. If my words were not sufficiently convincing, then maybe a few pictures might help you see what I'm seeing:

With[{m = 3, z = 1},
Table[Plot3D[With[{β = u + I v}, f[Exp[I m β + I z Cos[β]]]],
{u, -π, π}, {v, -6, 6},
BoundaryStyle -> None, BoxRatios -> {1, 3/2, 3/2}, ClippingStyle -> None,
Mesh -> False, PlotLabel -> f, PlotRange -> {-10, 10},
PlotStyle -> Directive[Pink, Specularity[White, 50], Opacity[0.8]]],
{f, {Re, Im}}] // GraphicsRow]


The pictures given above are plots of the real and imaginary parts of $\exp(i m \beta + i z \cos\beta)$ for $m=3$ and $z=1$, within the region $-\pi < \Re \beta < \pi$. As you might notice, the integrand varies rather wildly on the imaginary axis; it stands to reason that a contour that lies too near to it will give function values that vary wildly as well, making numerical computations of the contour integral unstable. You'll thus want a contour that tries to keep away from the "forest". In particular, contours that lie on $\Re \beta=-\pi/2$ for $\Im \beta > 0$ and $\Re \beta=\pi/2$ for $\Im \beta < 0$ are at a comfortable distance from the "forest", which makes them suitable for quadrature.

This should really be a comment, but I don't have enough reputation to do that, so I post it as an answer instead.

I took a close look at your first integration path, they all lay within region where the integrand is actually well-behaved: smooth change and only one or two oscillations, there shouldn't be any significant error in your Hmz function. So I changed the $\infty$ of integration bound to some number large enough for the integrand to vanish. The result with $\infty$ replaced by $10$ is:

$$0.0195634 - 5.82152 i$$

So the error should be caused by Mathematica's way of handling infinite integration bound, not the integrand.