16
$\begingroup$

This is an update to my question Rigorous numerical integration from three years ago.

Is there now a package for rigorous numerical integration that uses interval arithmetic and has access to a well-developed library of special functions?

By "well-developed", I mean something that, at the very least, includes the error function $\text{erf}$, the incomplete gamma function $\Gamma(s,x)$ and the like. In other words, I would like to integrate expressions involving these functions. VNODE-LP doesn't seem to satisfy this requirement, as it relies on older interval-arithmetic libraries lacking such basic special functions. (Besides, perhaps it is just me, but VNODE-LP actually seems harder to install than three years ago; I am not sure it is being adequately maintained.)

PS. An ideal library would have both arbitrary-precision and double-precision modes, the latter one fast. That isn't really an issue for what I am doing now, but it sounds like a basic desideratum, just as, say, an intuitive interface written in C++ (as opposed to plain C) would be.

Improve this question
$\endgroup$
4
  • 2
    $\begingroup$ You might have more luck at scicomp.stackexchange.com. $\endgroup$
    – Federico Poloni
    Commented Aug 28, 2016 at 17:34
  • $\begingroup$ I'm trying that too. $\endgroup$
    – H A Helfgott
    Commented Aug 28, 2016 at 21:09
  • $\begingroup$ the Maple tool intpakX will not do? at least the precision can be set arbitrarily. $\endgroup$
    – Carlo Beenakker
    Commented Aug 29, 2016 at 6:36
  • $\begingroup$ Does intpakX have integration? (And does it have a good library of special functions?) $\endgroup$
    – H A Helfgott
    Commented Aug 29, 2016 at 12:04

2 Answers 2

Reset to default
21
+50
$\begingroup$

I develop Arb, an arbitrary-precision interval library with special functions support. There is some code included for integration of complex analytic functions using Taylor series (documentation), which in principle should scale nicely to high precision. An example program is available here.

The last test case in the example program is $\tfrac{1}{\pi} \int_0^{\pi} \cos(10 t - (20+10i) \sin(t)) dt$ ($= J_{10}(20+10i)$). You can change $\cos$ to $\operatorname{erf}$ (_acb_poly_cos_series to _acb_hypgeom_erf_series) and check that it works. In fact, the new integrand peaks with a magnitude around $10^{38}$ while the result has magnitude around $10^{15}$, so this is a nice test case for handling severe cancellation. With 20d target accuracy and 60d working precision, the program gives

(-1.55...e+24 + 5.83...e+23j)  +/-  (1.27e+103, 1.27e+103j)

where the resulting midpoint is inaccurate but correctly enclosed, and with 80d target accuracy and 240d working precision, it gives

(868463919173069.624... + 85377491518329.173...j)  +/-  (1.9e-65, 2.37e-66j).

This takes 1.5 seconds.

I would not necessarily recommend using this code for anything serious, for a few reasons. First, the interface is very low level and inconvenient to work with. Some kind of high-level wrapper would make sense. (Note that Arb is a Sage package, but no Sage wrapper for this particular code has been written yet.)

Further, the code needs to bound the integrand on a surrounding contour to get error bounds for the series truncations (via the Cauchy integral formula). This is done semi-automatically by evaluating the integrand on "thick" input intervals, subdividing if necessary. This can be slow and might fail to converge. For low precision integration, computing the original integral directly by naive subdivision might be better.

The current version of the code also requires that the user input a radius that is known a priori to isolate the integration path from any singularities. Singularities on the integration path itself are obviously not supported at all.

A more robust way to do integration is with Taylor models, i.e. truncated Taylor series together with bounds on the approximation error that are propagated automatically upon function composition. This is more efficient, and works natively with real analytic functions, or even non-analytic functions like $|x|$. Taylor models could presumably also be extended to Puiseux-type series to support algebraic and logarithmic endpoint singularities. Unfortunately, I'm not aware of any work so far on implementing Taylor models for higher special functions. It's on my list of things to work on...

Some people have also worked on rigorous error bounds for the double exponential integration formula, but I'm not aware of any general-purpose implementations yet.

While this doesn't quite answer your question, you might have some success combining Arb's special functions with other packages for integration. The main caveat is that Arb doesn't always compute very tight enclosures for special functions, which is hard to do in general. For example, the $\operatorname{erf}(x)$ implementation currently computes reasonably tight enclosures, but the $\Gamma(s,x)$ implementation does not (see this). This doesn't matter much for algorithms that only sample the integrand at discrete points (it's enough to increase the precision by $p$ bits to circumvent $p$ bits of cancellation), but it can be an obstacle for algorithms that need enclosures on "thick" input intervals (subdividing requires $2^p$ function evaluations).

Improve this answer
$\endgroup$
11
  • $\begingroup$ Thanks - I'll try this. Still, I do not quite get why the following simpler strategy (which requires no use of complex analysis) couldn't be automated. $\endgroup$
    – H A Helfgott
    Commented Aug 31, 2016 at 19:51
  • $\begingroup$ Let us work with the simplest kind of numerical integration - higher-order versions should work in the same way. We want to approximate the integral by a sum of type $\sum_{n=1}^N f(a_i)$, where $a_i = a + i (b-a)/N$. The error term is $\leq ((b-a)/N)^2 \sum_{n=1}^N \max_{a_{i-1}\leq x\leq a_i} |f'(x)|$. Now, we can compute the maximum within the sum simply by using interval arithmetic to bound the image of [a_{i-1},a_i] under $f'$. All we need, then, is interval arithmetic and symbolic differentiation. $\endgroup$
    – H A Helfgott
    Commented Aug 31, 2016 at 20:04
  • $\begingroup$ Of course I can do this by differentiating $f$ by hand (or in SAGE) and then typing in the resulting expression (very hairy, in what I am doing) - however, this is error-prone, and the procedure ought to be automated. think?What do you $\endgroup$
    – H A Helfgott
    Commented Aug 31, 2016 at 20:05
  • $\begingroup$ You're right, and that's indeed easy to implement. The drawback is that for $(d+1)$th order error bounds, you need to bound $|f^{(d)}|$ on intervals, which can be difficult when $d$ is large. With the complex method, you only need to bound $|f|$ itself on intervals; the bounds can then be reused for any $d$. (Taylor methods avoid evaluation on intervals completely.) All of this only really matters if you need so high precision that say $d = 2$ or $d = 4$ would require too large $N$; for modest precision, it would make sense most of the time to use a low order method and increase $N$. $\endgroup$
    – Fredrik Johansson
    Commented Aug 31, 2016 at 23:33
  • $\begingroup$ In Arb, I'm fairly systematically implementing special functions with support for truncated power series input and output (= automatic differentiation), so you can get arbitrary derivatives easily. Again, this isn't yet exposed in Sage, though it might be of some use to note that I have a separate Python interface that does wrap the power series conveniently. For example, you can do ctx.cap=3; x=arb_series([0,1]); (1+x*x).sqrt().erf() to produce $a_0 + a_1 x + a_2 x^2 + O(x^3) = \operatorname{erf}(\sqrt{1 + x^2})$. $\endgroup$
    – Fredrik Johansson
    Commented Aug 31, 2016 at 23:43
1
$\begingroup$

From what I read in the description of INTLAB (a Matlab/Octave "Interval Laboratory"), it has some of the desired capabilities, including gamma, erf and erfc functions, integration of univariate functions, multiple precision interval arithmetic, ...

Improve this answer
$\endgroup$

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .