2012
Oct
03

## Your computer algebra system is out to get you!

OK, maybe it's not what you think, the CAS's of the world are not going to rise up and enslave us :-P (Not yet, anyway.) But it is always a good idea to view their results with a healthy dose of skepticism, because like any computer program, they have bugs, and every once in a while one can actually mess up your calculation.

Here's my example:

$\int_0^{2\pi} \cos x\,e^{ik\cos(x - y)}\udc x$

Plug this into Mathematica 8.0 and it will tell you that the integral is equal to $2\pi i J_1(k)$, where $J_1$ is a Bessel function of the first kind:

In[1]  := Integrate[Cos[x] Exp[I k Cos[x - y]], {x, 0, 2 Pi}]
Out[1] := ConditionalExpression[2 I \[Pi] BesselJ[1, k], k \[Element] Reals]


But it's not! If you shift the integration variable $x\to x + y$, Mathematica tips its hand and tells you the real answer, $2\pi i J_1(k)\cos y$:

In[2]  := Integrate[Cos[x + y] Exp[I k Cos[x]], {x, 0, 2 Pi}]
Out[2] := 2 Pi I BesselJ[1, k] Cos[y]


You can also do this integral out by hand (useful because when Mathematica contradicts itself once, maybe you can't be sure which of its results to trust). One way is to use the relation

$e^{iz\cos\theta} = \sum_{n = -\infty}^{\infty}i^nJ_n(z)e^{in\theta}$

You can plug this with $z\to k, \theta\to x-y$ right into the integral, interchange the order of the summation and the integration, and get

$\sum_{n = -\infty}^{\infty}\int_0^{2\pi} i^n J_n(k) \cos x\, e^{in(x - y)}\udc x$

The cosine splits into two terms, $\cos x = \frac{1}{2}(e^{ix} + e^{-ix})$, so with a little algebraic rearrangement you get

$\sum_{n = -\infty}^{\infty} i^n J_n(k) \frac{e^{-iny}}{2} \int_0^{2\pi} \bigl(e^{i(n+1)x} + e^{i(n-1)x}\bigr)\udc x$

Now, $\int_0^{2\pi}e^{imx}\udc x$ is zero except when $m = 0$, when it is equal to $2\pi$. So the two terms in that integral are zero except in the $n=-1$ and $n=+1$ terms of the sum over $n$, respectively.

$-i J_{-1}(k) \frac{e^{iy}}{2}(2\pi) + i J_{1}(k) \frac{e^{-iy}}{2}(2\pi)$

And finally, there is a useful identity for remapping Bessel functions of negative integer order on to their positive-order counterparts, i.e. in this case $J_{-1}(k) = -J_{1}(k)$, so you're left with

$i J_{1}(k) \frac{e^{iy}}{2}(2\pi) + i J_{1}(k) \frac{e^{-iy}}{2}(2\pi) = 2\pi i J_{1}(k) \cos y$

Hooray!

blog comments powered by Disqus