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.
HankelH1[]
is built-in? $\endgroup$PrecisionGoal
andWorkingPrecision
. You can play with them a bit. On the other hand you can get the both integrals the same assuming an appropriateMethod
, e.g. adding this option to your definitionsMethod -> "NewtonCotesRule"
provides the same correct values. $\endgroup$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 frommyHankelH1[]
compares favorably to the built-in function. (As a matter of fact, usingMethod -> "DoubleExponential"
on your two functions seems to work nicely.) $\endgroup$"GaussKronrodRule"
and"ClenshawCurtisRule"
work faster than"NewtonCotesRule"
. $\endgroup$