Rigorous numerical integration

October 3, 2013

As promised in the last post, I have started a complex analysis module in Arb. It is now possible to evaluate integrals

$$\int_a^b f(t) dt$$

where $f$ is a given analytic function and $a$ and $b$ are (finite) complex numbers. With rigorous error bounds, and to arbitrary precision!

As a consequence, it suddenly becomes possible to rigorously evaluate various special functions with appropriate integral representations (such as: elliptic integrals, the error function, exponential integrals, Bessel functions). This is a fairly dumb way to evaluate any of those special functions, but might do fine if the arguments are not too large and you just want a couple of values!

The algorithm consists of approximating the integrand by high-order Taylor polynomials, integrating the polynomials, and using the Cauchy integral formula to bound the truncation error. The user must provide an inner radius (used for the Taylor series expansions) and an outer radius (used for the Cauchy bound). The outer radius must be small enough that any singularities of the integrand are avoided (the user must therefore have some knowledge about the integrand, though it would be possible to write a function that finds an admissible outer radius automatically).

There is some art involved in choosing the best combination of parameters (depending on the specific integral): a smaller inner radius gives faster convergence but requires more evaluation points; a larger outer radius gives faster convergence at higher precision but also adds a lot of overhead.

In general, this algorithm is not very efficient, because computing many Taylor series coefficients is typically more expensive than doing many pointwise evaluations. It is conceptually simple, though, and there are special cases where it can actually perform very well.

A rather severe limitation is that the integrand really must be analytic, even at the endpoints, though it should be possible to extend the algorithm in the future to support explicitly factored-out endpoint singularities at least of the simplest type, e.g. allowing $f(t) = \int_a^b (t-a)^{\alpha} (b-t)^{\beta} g(t)$ with $\alpha, \beta, g$ given as input.

I have a work-in-progress demo program that evaluates three integrals.

The first is an elementary integral:

$$I_1 = \int_0^{100} \sin t \, dt$$

The second defines an incomplete elliptic integral of the first kind:

$$I_2 = F(\phi,m) = \int_0^{\phi} \frac{1}{\sqrt{1-m \sin^2 t}} \, dt$$

The third defines a Bessel function:

$$I_3 = J_n(z) = \frac{1}{\pi} \int_0^{\pi} \cos(tn-z \sin t) \, dt$$

We evaluate $I_2$ with $m = 1/2, \phi = 6+6i$. The integrand has singularities (at the points $(n+1/2) \pi \pm \log(1+\sqrt 2) i, n \in \mathbb{Z}$), so we must choose a path that gives the intended branch cut and a radius bound that separates the path from the singularities. We use the two-segment path $0 \to 6 \to 6+6i$ (this gives the branch cut for F that is consistent with the ellipf function in mpmath).

We evaluate $I_3$ with $n = 10$, $z = 20 + 10i$.

For reference, the values are approximately (using mpmath):

>>> mp.dps = 30
>>> print cos(0) - cos(100)
0.137681127712316065898061486049
>>> print ellipf(6+6j, 0.5)
(7.41433976790976840857250785206 + 1.84734298077474928677131567616j)
>>> print besselj(10, 20+10j)
(596.964215758266079893402499222 + 326.082687277832010489600057171j)

Here is the output of the program, with 30 digits:

(0.137681127712316065898061486048 + 0j)  +/-  (9.92e-29, 7.14e-29j)
cpu/wall(s): 1.27 1.394

(7.41433976790976840857250785205 + 1.84734298077474928677131567616j)  +/-  (2.49e-26, 2.49e-26j)
cpu/wall(s): 0.27 0.269

(596.96421575826607989340249922 + 326.08268727783201048960005717j)  +/-  (3.1e-27, 1.62e-27j)
cpu/wall(s): 2.45 2.454

100 digits:

(0.1376811277123160658980614860491574644899159914644891707198378873072789119490733758969048943157227149 + 0j)  +/-  (6.38e-100, 6.38e-100j)
cpu/wall(s): 1.42 1.545

(7.414339767909768408572507852061324751189822032351426492644763432404197856221946818841218954992362325 + 1.847342980774749286771315676163486111543446141032818463594864121106831066453673386671719988277396978j)  +/-  (9.48e-100, 2.62e-100j)
cpu/wall(s): 0.73 0.731

(596.9642157582660798934024992220192684521348281360896022586915678041303218514604715239969732511226485 + 326.0826872778320104896000571707543269024938241093598189516330860745933923186652586820087014046850108j)  +/-  (2.24e-97, 1.17e-97j)
cpu/wall(s): 3.33 3.332

300 digits:

(0.137681127712316065898061486049157464489915991464489170719837887307278911949073375896904894315722714932864392444837669518894471931980661458906553793051115068984106183459664059666773393595959288592869686373065385439151640649879054463782064508146529471957980879841224514023584138013318423417984137632747 + 0j)  +/-  (3.21e-300, 3.21e-300j)
cpu/wall(s): 1.79 1.926

(7.4143397679097684085725078520613247511898220323514264926447634324041978562219468188412189549923623251566956438639953295695770605410270673115560512440033820571777902211147943411850962243637143164586102261799674430882651622752124627843210626527003790788555082128824014775933331176228142376021223472975 + 1.84734298077474928677131567616348611154344614103281846359486412110683106645367338667171998827739697796984424819881452682759690239505182352096597727525387956138651930750604129983788765757287394003899415339522809137192108561098920269805739597828700917063531775146524096353778246948853516104000684186535j)  +/-  (1.21e-299, 3.13e-300j)
cpu/wall(s): 4.82 4.826

(596.964215758266079893402499222019268452134828136089602258691567804130321851460471523996973251122648617462219670949426902892296098782688463795925609029007129880820159114690096234183689878456781035491966023711148660594815823674343542010947756638434097480007609768484852174367316584740218820032971668111 + 326.082687277832010489600057170754326902493824109359818951633086074593392318665258682008701404685010831106763353048630147084498489960527131607737608751912171247853797898416224954765015116463421805124531090166990515228922943510982793932883627473390959458940870940801400761943735359287238906056404067229j)  +/-  (2.93e-297, 1.53e-297j)
cpu/wall(s): 13.07 13.074

Here is how long it takes to evaluate the same integrals using the mpmath quad function:

          30 digits    100 digits     300 digits
I1        0.0047       0.016          0.22
I2        0.022        0.19           1.53
I3        0.031        0.14           0.65

Indeed, mpmath.quad is at least an order of magnitude faster, which is unsurprising since it uses a sample-based quadrature algorithm (tanh-sinh a.k.a. double-exponential quadrature), with only heuristic estimation of the error. Now the challenge is to try to match (or beat) the speed of mpmath with an algorithm that actually includes rigorous error bounds (for example, Gaussian or tanh-sinh quadrature with explicit error bounds based on the endpoint derivatives)!

fredrikj.net | Blog index | RSS feed | Follow me on Mastodon | Become a sponsor