Mathematica and hypergeometric functions

Recently I did some numerical experiments in Mathematica involving the hypergeometric function _2F_1. The results were clearly wrong (a positive-definite matrix having negative eigenvalues, for example), so I spent a couple of hours checking the code. Finally, I discovered a bug, but it was not in my program!

This post is not about the hypergeometric function, it is enough to know that it accepts four arguments, usually typed as _2F_1(a,b;c;x). For some parameters a, b, c it reduces to elementary functions, for example _2F_1(\frac{1}{2},\frac{1}{2};\frac{3}{2};-x) = \mathrm{arsinh}(\sqrt{x}) / \sqrt{x}. Of course, Wolfram Research products are aware of that (link). A sligth modification of the parameters of the hypergeometric function does not change the value of the function significantly, it is smooth in all four arguments. But according to Wolfram Alpha, _2F_1(\frac{1}{2},\frac{1}{2}-\frac{1}{1000};\frac{3}{2}-\frac{1}{1000};-x) is not even continuous! Click on an image to see the plot at the Wolfram Alpha site.

There is clearly something wrong with the way Wolfram Alpha (and its kernel, Mathematica) computes _2F_1. This issue is related to finite precision calculations: forcing Mathematica to work with higher precision helps a lot, but it does not resolve the problem completely. Numerical instability is nothing spectacular, it is a common issue in approximate computations. Here, however, the error cannot be explained simply by accumulation of numerical errors, the result is just completely wrong: zero instead of somethign close to 0.8. This kind of bug is quite unexpected in such an advanced application.

Let us play around this bug for a while. Wolfram Alpha asked about _2F_1(\frac{1}{2},\frac{1}{2}-\frac{1}{1000};\frac{3}{2}-\frac{1}{1000};-2), answers correctly (link). The numerical approximation (N[…]) is “result: 0.81072238…”, which is fine; the next line, however, reads “number name: zero” (link). Things get even worse if we request five digits only: “decimal approximation: 0.810722, result: 0.0” (link).

I was in fact quite lucky to find this issue: there are not too many parameters for which _2F_1 is not computed accurately. In the contour plot of _2F_1(\frac{1}{2},\frac{1}{2}+p;\frac{3}{2}+q;-2) they are all concentrated near the line p = q, when p is close to 1. And this is exactly the case I was interested in.

Nearly a week ago, I submitted this bug to Wolfram Support Center. Unfortunately, apart from an automated confirmation letter (claiming that they typically answer within three business days), I have received no reply.

The above bug might be caused by Mathematica using a wrong hypergeometric identity (perhaps this one?) in numerical evaluation. If this is the case, it will be rather easy to fix. Until then, one can simply use another hypergeometric relation (like this one) explicitly. For an example, see the last two (theoretically identical) plots on the right.

Let me end this post with a few words about what I was trying to do when I discovered this bug. Unfortunately, this will be a little bit technical for those not used to non-local operators and jump processes.

Let D be an open set in \mathbf{R}^d. Consider a Brownian motion X_t (with variance \mathrm{Var} \, X_t = 2t) starting at a fixed point x \in D, and let \tau_D denote the (random) first time when X_t hits the complement of D. The mean amount of time that X_t spent in a Borel set A \subseteq D before \tau_D can be expressed as \int_A G_D(x, y) dy for some function G_D(x, y). It can be proved that G_D(x, y) is the Green function for -\Delta_D, the Laplace operator in D with zero Dirichlet boundary conditions. In fact, this is probably the easiest way to construct the Green function for a general open set.

Replace in the above paragraph the Brownian motion X_t with the isotropic \alpha-stable Lévy process (here \alpha \in (0, 2)), and G_D(x, y) will become the Green function for the fractional power of the Laplace operator (-\Delta)^{\alpha/2}, with zero exterior condition (this is no longer a local operator, so instead of boundary conditions, one must give exterior condition). This operator is usually denoted as ((-\Delta)^{\alpha/2})_D, and it is important to note that this is something else than (-\Delta_D)^{\alpha/2}. Fractional powers of -\Delta are inverse operators to Riesz potentials, and they are important examples of pseudo-differential operators. It is not my point to give a detailed introduction to probabilistic potential theory and Riesz potential theory herem and there is much literature on these topics, inlcuding classical textbooks by Landkoff, Bliedtner and Hansen, and Blumenthal and Getoor.

If D is a ball, then the formula for the Green function for -\Delta_D is relatively easy to find using Kelvin transform. Suprisingly, similar method works for ((-\Delta)^{\alpha/2})_D, due to the calculation of M. Riesz in 1938, in this case the formula for the Green function for a ball is:

\displaystyle \begin{aligned} G_D(x, y) & = \frac{\Gamma(\frac{d}{2}) |x - y|^{\alpha - d}}{2^\alpha \pi^{d/2} (\Gamma(\frac{\alpha}{2}))^2} \int_0^{\frac{(1 - x^2)(1 - y^2)}{|x - y|^2}} \frac{s^{\alpha/2 - 1}}{(1 + s)^{d/2}} \, ds \\ & = \frac{\Gamma(\frac{d}{2}) (1 - x^2)^{\alpha / 2} (1 - y^2)^{\alpha / 2}}{2^\alpha \pi^{d/2} \Gamma(\frac{\alpha}{2}) \Gamma(1 + \frac{\alpha}{2}) |x - y|^{d}} \, {_2 F_1}\Bigl(\frac{\alpha}{2}, \, \frac{d}{2}; \, 1 + \frac{\alpha}{2}; \, -\frac{(1 - x^2)(1 - y^2)}{|x - y|^2}\Bigr) . \end{aligned}

While I was working on numerical bounds for the eigenvalues of ((-\Delta)^{\alpha/2})_D for D = (-1, 1) (that is, a one-dimensional ball), I needed to compute an array of values of G_D(x, y). Numerical approximation to the integral worked fine, but it was very slow. For that reason, I tried the formula with the hypergeometric function. It worked much faster, but the results were different (and plainly wrong). While I was checking the code, I plotted the graph of G_D(x, y) for \alpha = 0.99. The correct picture should be like the one below (prepared using a better hypergeometric identity).

Instead, Mathematica produced the following image.

The erroneous plot can be reproduced in lower resolution using Wolfram Alpha (link).

I will update this post as soon as Wolfram contacts me.

Update, Mar 11, 2011:

I have just received a copy of Mathematica 8 and istalled it on my computer. The first thing I checked was, clearly, whether Wolfram fixed the hypergeometric bug. And they did! A huge surprise to me, as they did not respond to any of my two bug submissions. And the bug is still present on Wolfram Alpha! Does Wolfram Alpha run an older version of Mathematica then?

Update, Mar 23, 2011:

I have just received a kind email from Wolfram Research confirming that the issue has been resolved in Mathematica 8.

This entry was posted in bugs, numerical math and tagged , , , . Bookmark the permalink.

2 Responses to Mathematica and hypergeometric functions

  1. juanmarqz says:

    by this way, is it posible to receive a free-copy of Mathematica8? he.he

Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s