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


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)#

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

sage.numerical.gauss_legendre.integrate_vector(f, prec, epsilon=None)#

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


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

sage.numerical.gauss_legendre.integrate_vector_N(f, prec, N=3)#

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)


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)#

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

sage.numerical.gauss_legendre.nodes_uncached(degree, prec)#

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


Todo

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