Take the 2-minute tour ×

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):

admissible region for Hankel function contour

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):

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

share|improve this question
I presume you're aware that HankelH1[] is built-in? –  Ｊ. Ｍ. Jun 7 '13 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. –  Matthew Jun 7 '13 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. –  Artes Jun 7 '13 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.) –  Ｊ. Ｍ. Jun 7 '13 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". –  Ｊ. Ｍ. Jun 7 '13 at 9:31

1 Answer 1

up vote 2 down vote accepted

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

integrand for Hankel function contour integral

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.

share|improve this answer

Your Answer


