Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question
    
I presume you're aware that HankelH1[] is built-in? –  J. M. 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
1  
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.) –  J. M. Jun 7 '13 at 9:07
1  
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". –  J. M. Jun 7 '13 at 9:31
show 7 more comments

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
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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