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

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.

  • $\begingroup$ I presume you're aware that HankelH1[] is built-in? $\endgroup$ Commented Jun 7, 2013 at 8:35
  • $\begingroup$ 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. $\endgroup$
    – George
    Commented Jun 7, 2013 at 8:39
  • $\begingroup$ @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. $\endgroup$
    – Artes
    Commented Jun 7, 2013 at 8:48
  • 1
    $\begingroup$ 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.) $\endgroup$ Commented Jun 7, 2013 at 9:07
  • 1
    $\begingroup$ 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". $\endgroup$ Commented Jun 7, 2013 at 9:31

2 Answers 2


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


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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.