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_vectorby 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
precbits of precision. If the vector space in whichftakes values is over a field which is incompatible with this field (e.g. a finite field) then aTypeErroroccurs.
- 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.chas better performance.