Tuesday, June 15, 2010

Assorted special functions update

Over the week-and-a-half since the last blog update, I've gotten a bunch more work done on special functions in mpmath.

Function plots in the documentation

I have started adding graphics of special functions to the documentation. So far, I've done most of the Bessel-type functions and orthogonal polynomials. Many more to come!

Example screenshot (see the Bessel functions page):

New inhomogeneous Bessel functions

There exists a large number of lesser-known special functions which are essentially variations of Bessel functions. These include functions which solve the generalized (inhomogeneous) Bessel differential equation

with some specific right-hand side g(z). New additions to mpmath in this category are the Anger function (angerj()), Weber function (webere()), and Lommel functions (lommels1(), lommels2()). See commits here and here.

More information about Anger-Weber functions and Lommel functions can be found in the DLMF.

In the near future, I will probably further improve the implementations of the main Bessel functions. The Bessel functions are mostly implemented as generic hypergeometric functions, but the standard cases can be tuned a great deal with special-purpose code.

Airy functions and related functions

The Airy functions Ai and Bi have been present for quite some time in mpmath. In a recent commit, I have rewritten them for improved rigor and better performance at high precision. There are also some new features, such as the ability to evaluate derivatives or iterated integrals of arbitrary order.


>>> from mpmath import *
>>> mp.dps = 25; mp.pretty = True
>>> airyai(1.5, derivative=5)
>>> diff(airyai, 1.5, 5)
>>> airyai(1.5, derivative=100)
>>> airybi(0, derivative=1000)
>>> airybi(0, derivative=1001)
>>> airybi(0, derivative=1002)


>>> airyai(5, derivative=-1)
>>> quad(airyai, [0,5])
>>> airyai(-100000, derivative=-1)

Also, functions for computing the zeros of Ai and Bi (and the first derivatives) have been added:

>>> airyaizero(1)
>>> airyaizero(2)
>>> airybizero(1)
>>> airybizero(1, derivative=1)
>>> airybizero(1, derivative=1, complex=True)
(0.2149470745374305676088329 + 1.100600143302797880647194j)
>>> airybizero(10000)
>>> airybizero(10000, complex=True)
(652.3059222438076432024695 + 1129.846189716375208308414j)

I have also implemented two new functions related to Airy functions: the Scorer functions Gi and Hi. These are available as scorergi() and scorerhi() respectively.

Here are two plots of the Gi-function, which can also be seen in the documentation:

Interval gamma functions

The interval arithmetic context now implements gamma, rgamma (reciprocal gamma function), factorial as well as loggamma for real as well as complex arguments (commit). For example:

>>> iv.dps = 10
>>> iv.gamma('50.3')
[1.96282982095908e+63, 1.96282982457481e+63]
>>> iv.gamma(iv.mpc('2.7','5.9'))
([0.00269836072064322, 0.00269836072271801] +
[0.0120124287790304, 0.0120124287810768]*j)

As a "practical" example, consider evaluating the Riemann-Siegel theta function which involves computing the difference of two log-gamma functions. For input with a large real part, the imaginary part in the result suffers from massive cancellation and may end up with the wrong sign:

>>> mp.dps = 15
>>> mp.siegeltheta(10**50 + 0.25j)
(5.61456887916465e+51 - 0.143091235731175j)
>>> mp.dps = 10; nprint(mp.siegeltheta(10**50 + 0.25j).imag)
>>> mp.dps = 100; nprint(mp.siegeltheta(10**50 + 0.25j).imag)

With interval arithmetic, the sign uncertainty is reflected in the output:

>>> iv.dps = 15
>>> iv.siegeltheta(10**50 + 0.25j)
([5.6145688791646467648e+51, 5.6145688791646474294e+51] +
[-5.0706024009129187319e+30, 5.070602400912917606e+30]*j)
>>> iv.dps = 50
>>> iv.siegeltheta(10**50 + 0.25j)
([5614568879164647368060513633451316140100495086670736.0, 5614568879164647368060
513633451316140100495086670744.0] +
8814838558203, 14.16147419395632497823207158108086766104408814838560341]*j)

As another example, consider evaluating the gamma function of a huge argument. The digits in the answer may be "wrong" because the input is converted from decimal to binary, and the gamma function is sensitive to the input being perturbed:

>>> mp.dps = 15
>>> mp.gamma('123456789012345.1')
>>> mp.dps = 30
>>> mp.gamma('123456789012345.1')
>>> mp.dps = 60
>>> mp.gamma('123456789012345.1')

With interval arithmetic, the uncertainty in the input is propagated correctly:

>>> iv.dps = 15
>>> iv.nprint(iv.gamma('123456789012345.1'), mode='diff')
[6.11545e+1686076589184486, 1.01533e+1686076589184487]
>>> iv.dps = 30
>>> iv.nprint(iv.gamma('123456789012345.1'), 20, mode='diff')
7.4903201854034[185631, 219359]e+1686076589184486
>>> iv.dps = 60
>>> iv.nprint(iv.gamma('123456789012345.1'), 50, mode='diff')
7.49032018540342058679709881225047421518964527[60898, 87505]e+1686076589184486

Rewritten Lambert W function

Lastly, the Lambert W function has received a much-needed rewrite (commit) mainly to improve evaluation very close to the branch cut along the negative axis and particularly near the branch point at -1/e for the k = -1, 0, 1 branches.

With the previous implementation, results were frequently inaccurate or ended up on the wrong branch in this region. Here are some hard cases that now work perfectly:

>>> mp.dps = 1000
>>> x = -1/e + mpf('1e-900')
>>> y = -1/e - mpf('1e-900')
>>> z = -1/e + mpf('1e-900')*1j
>>> w = -1/e - mpf('1e-900')*1j
>>> mp.dps = 25
>>> lambertw(x,0); lambertw(y,0); lambertw(z,0); lambertw(w,0)
(-1.0 + 2.331643981597124203363536e-450j)
(-1.0 + 1.648721270700128146848651e-450j)
(-1.0 - 1.648721270700128146848651e-450j)
>>> lambertw(x,1); lambertw(y,1); lambertw(z,1); lambertw(w,1)
(-3.088843015613043855957087 + 7.461489285654254556906117j)
(-3.088843015613043855957087 + 7.461489285654254556906117j)
(-3.088843015613043855957087 + 7.461489285654254556906117j)
(-1.0 + 1.648721270700128146848651e-450j)
>>> lambertw(x,-1); lambertw(y,-1); lambertw(z,-1); lambertw(w,-1)
(-1.0 - 2.331643981597124203363536e-450j)
(-1.0 - 1.648721270700128146848651e-450j)
(-3.088843015613043855957087 - 7.461489285654254556906117j)

To finish this post, I present the following Mathematica bug:

remote1:frejohl:[~]$ math
Mathematica 7.0 for Linux x86 (32-bit)
Copyright 1988-2008 Wolfram Research, Inc.

In[1]:= Im[LambertW[0,-1/E-10^(-900)]]

Out[1]= 0


Aaron Meurer said...

Here's what Maple 12 gives:

> evalf(Im(LambertW(0, -1/exp(1)-10^(-900))));

Fredrik Johansson said...

The Maple result seems fine, because its evalf function isn't intended to do more than blindly evaluate the expression with floating-point arithmetic.

But it's a bug by Mathematica to claim symbolically that the imaginary part is zero (due to faulty error control in the numerics, most likely).

Kris said...

I think the mpmath documentation, with the new color figures, will become a nice resource. All the recent improvements you have been adding are great!