I have to perform a numerical integral on $[0,1]$ with a complicated integrand (possibly with integrable singularities at $0$ and/or $1$) and I am looking for a platform (sagemath, maple, python, mathematica, pari, matlab,…) that allows to do this in a satisfactory way. Ideally, though, I would like to go on working with sagemath.
The (real) integrand I am considering is a linear combination of functions of the form
$$\mathrm{Im}\left(\log\left(\frac{1-\zeta^a(1-t)^{1/N}}{\zeta^bt^{1/N}}\right) \cdot\frac{\zeta^bt^{1/N-1}-\zeta^a(1-t)^{1/N-1}}{\zeta^bt^{1/N}+\zeta^a(1-t)^{1/N}} \right),$$
where $t$ is the variable, $N,a,b \in \mathbb{N}$ are fixed and $\zeta=e^{2\pi i /N}$.
Concretely, I would like to be able to:
- Perform the integral with arbitrary precision, i.e., I will set the number of digits for the result and the program guarantees that the result is exact up to this digit. Probably this is just hopeless: one can always have pathological functions and since all numerical integrals are (to my knowledge) based on pointwise evaluation on some grid, one can never obtain a maximal bound on the precision, am I right?
- Alternatively and more humbly, I would also appreciate a way to set the number of digits of the result (even if they are not guaranteed to be all exact) and to get an error estimation (even if it is not a maximal bound).
Sagemath's "numerical_integral" deals with integrable singularities (can perform this integral) and provides an estimation of the error, although it is too big in some cases: that is why I need a way to reduce it. In principle it can also adapt the number of digits, but apparently this last feature does not work with the numerical integration: even changing the number of digits, the result is always the same. Also, one can in principle set the tolerance for the absolute and/or relative error of the numerical integration, but also this does not seem to change the result.
Any help/hint is appreciated, in particular on how to implement (part of) my goal with sagemath, but also for other platforms or for theoretical results.
Best Answer
mpmath
(an arbitrary precision mathematical library for Python 2.5 or higher, including Python 3.x) sounds exactly like what you'd be interested in. The numerical integration documentation is here. I used it to good effect in this MSE answer - not quite the same application, but I was computing an oscillatory integral on the interval $[0,\infty)$. The result was accurate to within $\sim10^{-50}$ after $\sim33$ seconds of calculation.I haven't personally used Sagemath, however, at least part of
mpmath
is documented on Sagemath.org, see this link.Edit: Yes, if you pass the option "
error=True
" to thequad
function, then you will receive an estimate of the integral & an estimate of the error. Below I've posted an example, from this MSE question, computed in Python 2.7.12. The only change I've made is to format the long decimal numbers in a more readable format. Unfortunately I don't know the details of the error estimation routine. Also, you can see in this instance that the estimate is not very precise.