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))
sage: max(c.abs() for c in (I - J)) < epsilon
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))))
sage: I = integrate_vector(f, prec, epsilon=epsilon)
sage: J = V((1/2, 0))
sage: max(c.abs() for c in (I - J)) < epsilon
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 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(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
sage: L1 = nodes(24, 53)
sage: P = RR['x'](sage.functions.orthogonal_polys.legendre_P(24, x))
sage: Pdif = P.diff()
sage: L2 = [((r + 1)/2, 1/(1 - r^2)/Pdif(r)^2)
....:        for r, _ in RR['x'](P).roots()]
sage: all((a[0] - b[0]).abs() < 1e-15 and (a[1] - b[1]).abs() < 1e-9
....:      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.