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

Question

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.

Answer 1

(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.

Answer 2

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.