numerical integration with arbitrary precision

[burcin]

From the sage-devel:

On Feb 20, 2010, at 12:40 PM, John H Palmieri wrote:
...
> I was curious about this, so I tried specifying the number of digits:
>
> sage: h = integral(sin(x)/x^2, (x, 1, pi/2)); h
> integrate(sin(x)/x^2, x, 1, 1/2*pi)
> sage: h.n()
> 0.33944794097891573
> sage: h.n(digits=14)
> 0.33944794097891573
> sage: h.n(digits=600)
> 0.33944794097891573
> sage: h.n(digits=600) == h.n(digits=14)
> True
> sage: h.n(prec=50) == h.n(prec=1000)
> True
>
> Is there an inherit limit in Sage on the accuracy of numerical
> integrals?  

The _evalf_ function defined on line 179 of sage/symbolic/integration/integral.py calls the gsl numerical_integral() function and ignores the precision.

We should raise a NotImplementedError for high precision, or find a way to do arbitrary precision numerical integration.

CC: @sagetrac-maldun @fredrik-johansson @kcrisman @sagetrac-mariah @sagetrac-bober @eviatarbach @mforets

Component: symbolics

Keywords: numerics, integration, sd32

Work Issues: add more arbitrary precision tests

Author: Stefan Reiterer

Reviewer: Paul Zimmermann

Issue created by migration from https://trac.sagemath.org/ticket/8321

[zimmermann6]

comment:1

Note that PARI/GP can do (arbitrary precision) numerical integration:

sage: %gp
gp: \p100
   realprecision = 115 significant digits (100 digits displayed)
gp: intnum(x=1,Pi/2,sin(x)/x^2)
0.3394479409789156796919271718652186179944769882691752577491475103140509215988087842006743201412102786

I don't know why it does not work from Sage:

sage: gp.intnum(x=1,2,sin(x)/x^2)
------------------------------------------------------------
   File "<ipython console>", line 1
SyntaxError: non-keyword arg after keyword arg (<ipython console>, line 1)

[aghitza]

comment:2

There is an example-doctest in the file interfaces/maxima.py, for the method nintegral. It says:

        Note that GP also does numerical integration, and can do so to very
        high precision very quickly::
        
            sage: gp('intnum(x=0,1,exp(-sqrt(x)))')            
            0.5284822353142307136179049194             # 32-bit
            0.52848223531423071361790491935415653021   # 64-bit
            sage: _ = gp.set_precision(80)
            sage: gp('intnum(x=0,1,exp(-sqrt(x)))')
            0.52848223531423071361790491935415653021675547587292866196865279321015401702040079

[mwhansen]

comment:3

mpmath also supports it and can handle python functions.

[sagetrac-maldun]

comment:4

I think I have the solution for this trac with mpmath:

def _evalf_(self, f, x, a, b, parent = None): 
        """
        Returns numerical approximation of the integral
        
        EXAMPLES::

            sage: from sage.symbolic.integration.integral import definite_integral
            sage: h = definite_integral(sin(x)/x^2, x, 1, 2); h
            integrate(sin(x)/x^2, x, 1, 2)
            sage: h.n() # indirect doctest
            0.4723991772689525...

        TESTS:

        Check if #3863 is fixed::

            sage: integrate(x^2.7 * e^(-2.4*x), x, 0, 3).n()
            0.154572952320790

        #Check if #8321 is fixed:

            sage: d = definite_integral(sin(x)/x^2, x, 1, 2)
            sage: d.n(77)
            0.4723991772689525199904
        """
        #from sage.gsl.integration import numerical_integral
        # The gsl routine returns a tuple, which also contains the error.
        # We only return the result.
        #return numerical_integral(f, a, b)[0]

        #Lets try mpmath instead:

        import sage.libs.mpmath.all as mpmath

        try:
            precision = parent.prec()
            mpmath.mp.prec = precision
        except AttributeError:
            precision = mpmath.mp.prec

        mp_f = lambda z: \
               f(x = mpmath.mpmath_to_sage(z,precision))

        return mpmath.call(mpmath.quad,mp_f,[a,b])

The tests just run fine:

sage: sage: sage: from sage.symbolic.integration.integral import definite_integral
sage: sage: h = definite_integral(sin(x)/x^2, x, 1, 2); h
integrate(sin(x)/x^2, x, 1, 2)
sage: h.n()
0.472399177268953
sage: h.n(77)
0.4723991772689525199904
sage: h.n(100)
0.47239917726895251999041133798

greez maldun

[sagetrac-maldun]

Attachment: trac_8321_numeric_int_mpmath.patch.gz

Numerical evaluation of symbolic integration with arbitrary precision with help of mpmath