Gauss-Legendre Integration for Vector-Valued Functions#

Routine to perform Gauss-Legendre integration for vector-functions.

EXAMPLES:

We verify that \(\int_0^1 n x^{n-1} \, dx = 1\) for \(n=1, \dots, 4\):

sage: from sage.numerical.gauss_legendre import integrate_vector
sage: prec = 100
sage: K = RealField(prec)
sage: N = 4
sage: V = VectorSpace(K, N)
sage: f =  lambda x: V([(n+1)*x^n for n in range(N)])
sage: I = integrate_vector(f, prec)
sage: max([c.abs() for c in I-V(N*[1])])
0.00000000000000000000000000000

>>> from sage.all import *
>>> from sage.numerical.gauss_legendre import integrate_vector
>>> prec = Integer(100)
>>> K = RealField(prec)
>>> N = Integer(4)
>>> V = VectorSpace(K, N)
>>> f =  lambda x: V([(n+Integer(1))*x**n for n in range(N)])
>>> I = integrate_vector(f, prec)
>>> max([c.abs() for c in I-V(N*[Integer(1)])])
0.00000000000000000000000000000

AUTHORS:

  • Nils Bruin (2017-06-06): initial version

  • Linden Disney-Hogg (2021-06-17): documentation and integrate_vector method changes

Note

The code here is directly based on mpmath (see http://mpmath.org), but has a highly optimized routine to compute the nodes.

sage.numerical.gauss_legendre.estimate_error(results, prec, epsilon)[source]#

Routine to estimate the error in a list of quadrature approximations.

The method used is based on Borwein, Bailey, and Girgensohn. As mentioned in mpmath: Although not very conservative, this method seems to be very robust in practice.

The routine takes a list of vector results and, under assumption that these vectors approximate a given vector approximately quadratically, gives an estimate of the maximum norm of the error in the last approximation.

INPUT:

  • results – list. List of approximations to estimate the error from. Should be at least length 2.

  • prec – integer. Binary precision at which computations are happening.

  • epsilon – multiprecision float. Default error estimate in case of insufficient data.

OUTPUT:

An estimate of the error.

EXAMPLES:

sage: from sage.numerical.gauss_legendre import estimate_error
sage: prec = 200
sage: K = RealField(prec)
sage: V = VectorSpace(K, 2)
sage: a = V([1, -1])
sage: b = V([1, 1/2])
sage: L = [a + 2^(-2^i)*b for i in [0..5]]
sage: estimate_error(L, prec, K(2^(-prec)))
2.328235...e-10

>>> from sage.all import *
>>> from sage.numerical.gauss_legendre import estimate_error
>>> prec = Integer(200)
>>> K = RealField(prec)
>>> V = VectorSpace(K, Integer(2))
>>> a = V([Integer(1), -Integer(1)])
>>> b = V([Integer(1), Integer(1)/Integer(2)])
>>> L = [a + Integer(2)**(-Integer(2)**i)*b for i in (ellipsis_range(Integer(0),Ellipsis,Integer(5)))]
>>> estimate_error(L, prec, K(Integer(2)**(-prec)))
2.328235...e-10
sage.numerical.gauss_legendre.integrate_vector(f, prec, epsilon=None)[source]#

Integrate a one-argument vector-valued function numerically using Gauss-Legendre.

This function uses the Gauss-Legendre quadrature scheme to approximate the integral \(\int_0^1 f(t) \, dt\).

INPUT:

  • f – callable. Vector-valued integrand.

  • prec – integer. Binary precision to be used.

  • epsilon – multiprecision float (default: \(2^{(-\text{prec}+3)}\)). Target error bound.

OUTPUT:

Vector approximating value of the integral.

EXAMPLES:

sage: from sage.numerical.gauss_legendre import integrate_vector
sage: prec = 200
sage: K = RealField(prec)
sage: V = VectorSpace(K, 2)
sage: epsilon = K(2^(-prec + 4))
sage: f = lambda t:V((1 + t^2, 1/(1 + t^2)))
sage: I = integrate_vector(f, prec, epsilon=epsilon)
sage: J = V((4/3, pi/4))                                                        # needs sage.symbolic
sage: max(c.abs() for c in (I - J)) < epsilon                                   # needs sage.symbolic
True

>>> from sage.all import *
>>> from sage.numerical.gauss_legendre import integrate_vector
>>> prec = Integer(200)
>>> K = RealField(prec)
>>> V = VectorSpace(K, Integer(2))
>>> epsilon = K(Integer(2)**(-prec + Integer(4)))
>>> f = lambda t:V((Integer(1) + t**Integer(2), Integer(1)/(Integer(1) + t**Integer(2))))
>>> I = integrate_vector(f, prec, epsilon=epsilon)
>>> J = V((Integer(4)/Integer(3), pi/Integer(4)))                                                        # needs sage.symbolic
>>> max(c.abs() for c in (I - J)) < epsilon                                   # needs sage.symbolic
True

We can also use complex-valued integrands:

sage: prec = 200
sage: Kreal = RealField(prec)
sage: K = ComplexField(prec)
sage: V = VectorSpace(K, 2)
sage: epsilon = Kreal(2^(-prec + 4))
sage: f = lambda t: V((t, K(exp(2*pi*t*K.0))))                                  # needs sage.symbolic
sage: I = integrate_vector(f, prec, epsilon=epsilon)                            # needs sage.symbolic
sage: J = V((1/2, 0))
sage: max(c.abs() for c in (I - J)) < epsilon                                   # needs sage.symbolic
True

>>> from sage.all import *
>>> prec = Integer(200)
>>> Kreal = RealField(prec)
>>> K = ComplexField(prec)
>>> V = VectorSpace(K, Integer(2))
>>> epsilon = Kreal(Integer(2)**(-prec + Integer(4)))
>>> f = lambda t: V((t, K(exp(Integer(2)*pi*t*K.gen(0)))))                                  # needs sage.symbolic
>>> I = integrate_vector(f, prec, epsilon=epsilon)                            # needs sage.symbolic
>>> J = V((Integer(1)/Integer(2), Integer(0)))
>>> max(c.abs() for c in (I - J)) < epsilon                                   # needs sage.symbolic
True
sage.numerical.gauss_legendre.integrate_vector_N(f, prec, N=3)[source]#

Integrate a one-argument vector-valued function numerically using Gauss-Legendre, setting the number of nodes.

This function uses the Gauss-Legendre quadrature scheme to approximate the integral \(\int_0^1 f(t) \, dt\). It is different from integrate_vector by using a specific number of nodes rather than targeting a specified error bound on the result.

INPUT:

  • f – callable. Vector-valued integrand.

  • prec – integer. Binary precision to be used.

  • N – integer (default: 3). Number of nodes to use.

OUTPUT:

Vector approximating value of the integral.

EXAMPLES:

sage: from sage.numerical.gauss_legendre import integrate_vector_N
sage: prec = 100
sage: K = RealField(prec)
sage: V = VectorSpace(K,1)
sage: f = lambda t: V([t])
sage: integrate_vector_N(f, prec, 4)
(0.50000000000000000000000000000)

>>> from sage.all import *
>>> from sage.numerical.gauss_legendre import integrate_vector_N
>>> prec = Integer(100)
>>> K = RealField(prec)
>>> V = VectorSpace(K,Integer(1))
>>> f = lambda t: V([t])
>>> integrate_vector_N(f, prec, Integer(4))
(0.50000000000000000000000000000)

Note

The nodes and weights are calculated in the real field with prec bits of precision. If the vector space in which f takes values is over a field which is incompatible with this field (e.g. a finite field) then a TypeError occurs.

sage.numerical.gauss_legendre.nodes(prec)[source]#

Compute the integration nodes and weights for the Gauss-Legendre quadrature scheme, caching the output

Works by calling nodes_uncached.

INPUT:

  • degree – integer. The number of nodes. Must be 3 or even.

  • prec – integer (minimal value 53). Binary precision with which the nodes and weights are computed.

OUTPUT:

A list of (node, weight) pairs.

EXAMPLES:

The nodes for the Gauss-Legendre scheme are roots of Legendre polynomials. The weights can be computed by a straightforward formula (note that evaluating a derivative of a Legendre polynomial isn’t particularly numerically stable, so the results from this routine are actually more accurate than what the values the closed formula produces):

sage: from sage.numerical.gauss_legendre import nodes
sage: L1 = nodes(24, 53)
sage: P = RR['x'](sage.functions.orthogonal_polys.legendre_P(24, x))            # needs sage.symbolic
sage: Pdif = P.diff()                                                           # needs sage.symbolic
sage: L2 = [((r + 1)/2, 1/(1 - r^2)/Pdif(r)^2)                                  # needs sage.symbolic
....:        for r, _ in RR['x'](P).roots()]
sage: all((a[0] - b[0]).abs() < 1e-15 and (a[1] - b[1]).abs() < 1e-9            # needs sage.symbolic
....:      for a, b in zip(L1, L2))
True

>>> from sage.all import *
>>> from sage.numerical.gauss_legendre import nodes
>>> L1 = nodes(Integer(24), Integer(53))
>>> P = RR['x'](sage.functions.orthogonal_polys.legendre_P(Integer(24), x))            # needs sage.symbolic
>>> Pdif = P.diff()                                                           # needs sage.symbolic
>>> L2 = [((r + Integer(1))/Integer(2), Integer(1)/(Integer(1) - r**Integer(2))/Pdif(r)**Integer(2))                                  # needs sage.symbolic
...        for r, _ in RR['x'](P).roots()]
>>> all((a[Integer(0)] - b[Integer(0)]).abs() < RealNumber('1e-15') and (a[Integer(1)] - b[Integer(1)]).abs() < RealNumber('1e-9')            # needs sage.symbolic
...      for a, b in zip(L1, L2))
True
sage.numerical.gauss_legendre.nodes_uncached(degree, prec)[source]#

Compute the integration nodes and weights for the Gauss-Legendre quadrature scheme

We use the recurrence relations for Legendre polynomials to compute their values. This is a version of the algorithm that in [Neu2018] is called the REC algorithm.

INPUT:

  • degree – integer. The number of nodes. Must be 3 or even.

  • prec – integer (minimal value 53). Binary precision with which the nodes and weights are computed.

OUTPUT:

A list of (node, weight) pairs.

EXAMPLES:

The nodes for the Gauss-Legendre scheme are roots of Legendre polynomials. The weights can be computed by a straightforward formula (note that evaluating a derivative of a Legendre polynomial isn’t particularly numerically stable, so the results from this routine are actually more accurate than what the values the closed formula produces):

sage: from sage.numerical.gauss_legendre import nodes_uncached
sage: L1 = nodes_uncached(24, 53)
sage: P = RR['x'](sage.functions.orthogonal_polys.legendre_P(24, x))            # needs sage.symbolic
sage: Pdif = P.diff()                                                           # needs sage.symbolic
sage: L2 = [((r + 1)/2, 1/(1 - r^2)/Pdif(r)^2)                                  # needs sage.symbolic
....:        for r, _ in RR['x'](P).roots()]
sage: all((a[0] - b[0]).abs() < 1e-15 and (a[1] - b[1]).abs() < 1e-9            # needs sage.symbolic
....:      for a, b in zip(L1, L2))
True

>>> from sage.all import *
>>> from sage.numerical.gauss_legendre import nodes_uncached
>>> L1 = nodes_uncached(Integer(24), Integer(53))
>>> P = RR['x'](sage.functions.orthogonal_polys.legendre_P(Integer(24), x))            # needs sage.symbolic
>>> Pdif = P.diff()                                                           # needs sage.symbolic
>>> L2 = [((r + Integer(1))/Integer(2), Integer(1)/(Integer(1) - r**Integer(2))/Pdif(r)**Integer(2))                                  # needs sage.symbolic
...        for r, _ in RR['x'](P).roots()]
>>> all((a[Integer(0)] - b[Integer(0)]).abs() < RealNumber('1e-15') and (a[Integer(1)] - b[Integer(1)]).abs() < RealNumber('1e-9')            # needs sage.symbolic
...      for a, b in zip(L1, L2))
True

Todo

It may be worth testing if using the FLINT/Arb algorithm for finding the nodes and weights in src/acb_calc/integrate_gl_auto_deg.c has better performance.