scipy.special.

roots_chebys#

scipy.special.roots_chebys(n, mu=False)[source]#

Gauss-Chebyshev (second kind) quadrature.

Compute the sample points and weights for Gauss-Chebyshev quadrature. The sample points are the roots of the nth degree Chebyshev polynomial of the second kind, \(S_n(x)\). These sample points and weights correctly integrate polynomials of degree \(2n - 1\) or less over the interval \([-2, 2]\) with weight function \(w(x) = \sqrt{1 - (x/2)^2}\). See 22.2.7 in [AS] for more details.

Parameters:
nint

quadrature order

mubool, optional

If True, return the sum of the weights, optional.

Returns:
xndarray

Sample points

wndarray

Weights

mufloat

Sum of the weights

References

[AS]

Milton Abramowitz and Irene A. Stegun, eds. Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables. New York: Dover, 1972.

Examples

>>> from scipy.special import roots_chebys
>>> roots, weights = roots_chebys(5)
>>> roots
array([-1.73205081e+00, -1.00000000e+00,  1.22464680e-16,  1.00000000e+00,
        1.73205081e+00])
>>> weights
array([0.26179939, 0.78539816, 1.04719755, 0.78539816, 0.26179939])

Verify that the values in roots are roots of the Chebyshev polynomial of first kind \(S_5(x)\).

>>> from scipy.special import eval_chebys
>>> eval_chebys(5, roots)
array([-1.33226763e-15, -1.77635684e-15,  3.67394040e-16, -8.88178420e-16,
        1.33226763e-15])

The values of \(S_5(x)\) evaluated at the roots are indeed very close to zero.

Verify that the sum of the weights equals the integral from -2 to 2 of \(\sqrt{1-(x/2)^2}\) which evaluates to \(\pi\). There are two ways to obtain the sum of weights, both resulting in \(\pi\) within numerical precision.

>>> sum(weights)
np.float64(3.141592653589793)
>>> roots, weights, sum_of_weights = roots_chebys(5, mu=True)
>>> sum_of_weights
3.141592653589793

Roots and weights obtained from the Chebyshev polynomial of the second kind \(S_n(x)\) are used in Gauss-Chebyshev quadrature where the integral from -2 to 2 of \(f(x)\sqrt{1-(x/2)^2}\) is evaluated. Roots and weights for order \(n\) are expected to yield the exact result for polynomials \(f(x)\) of a maximal order of \(2n-1\).

>>> f = lambda x: x**4
>>> weights @ f(roots)
np.float64(6.283185307179585)

The exact result is \(2\pi\).

>>> from math import pi
>>> 2*pi
6.283185307179586

In general, Gauss-Chebyshev quadrature will only yield an approximate value of the integral. Consider the integral from -2 to 2 of \(\cos(x)\sqrt{1-(x/2)^2}\).

>>> import numpy as np
>>> weights @ np.cos(roots)
np.float64(1.8118352216010754)

Check against the result of scipy.integrate.quad.

>>> from scipy.integrate import quad
>>> result, abserror = quad(lambda x: np.cos(x) * np.sqrt(1-(x/2)**2), -2, 2)
>>> result
1.8118344191919165

The latter result has an estimated absolute error of

>>> abserror
1.5300099409643053e-08