scipy.special.

roots_gegenbauer#

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

Gauss-Gegenbauer quadrature.

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

Parameters:

nint: quadrature order.
alphafloat: alpha must be > -0.5.
mubool, optional: If True, return the sum of the weights, optional.

Returns:

xndarray: Sample points
wndarray: Weights
mufloat: Sum of the weights

See also

scipy.integrate.fixed_quad

References

[AS]

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

Examples

Special cases of Gauss-Gegenbauer quadrature are the Gauss-Chebyshev first kind quadrature (\(\alpha=0\)) and Gauss-Chebyshev second kind quadrature (\(\alpha=1\)). Therefore, roots and weights obtained from roots_gegenbauer should agree with those obtained from roots_chebyt and roots_chebyu for the appropriate values of \(\alpha\).

>>> from scipy.special import roots_chebyt, roots_chebyu, roots_gegenbauer
>>> roots_gegenbauer(5, 0)  
(array([-0.95105652, -0.58778525,  0.        ,  0.58778525,  0.95105652]),
 array([0.62831853, 0.62831853, 0.62831853, 0.62831853, 0.62831853]))
>>> roots_chebyt(5)  
(array([-0.95105652, -0.58778525,  0.        ,  0.58778525,  0.95105652]),
 array([0.62831853, 0.62831853, 0.62831853, 0.62831853, 0.62831853]))

>>> roots_gegenbauer(5, 1)  
(array([-0.8660254, -0.5      ,  0.       ,  0.5      ,  0.8660254]),
 array([0.13089969, 0.39269908, 0.52359878, 0.39269908, 0.13089969]))
>>> roots_chebyu(5)  
(array([-8.66025404e-01, -5.00000000e-01,  6.12323400e-17,  5.00000000e-01,
         8.66025404e-01]),
 array([0.13089969, 0.39269908, 0.52359878, 0.39269908, 0.13089969]))

The sum of weights should equal the integral from -1 to 1 of \((1-x^2)^{\alpha-1/2}\) which evaluates to \(\sqrt{\pi}\Gamma(\alpha+1/2)/\Gamma(\alpha+1)\).

>>> alpha = 0.7
>>> roots, weights, sum_of_weights = roots_gegenbauer(5, alpha, mu=True)
>>> sum(weights)
np.float64(1.7910437497388674)
>>> sum_of_weights
np.float64(1.7910437497388672)
>>> from math import gamma, pi, sqrt
>>> sqrt(pi) * gamma(alpha+0.5) / gamma(alpha+1)
1.7910437497388667

Roots and weights obtained from the Gegenbauer polynomial \(C^\alpha_n(x)\) are used in Gauss-Gegenbauer quadrature where the integral from -1 to 1 of \(f(x)(1-x^2)^{\alpha-1/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(0.29265420747367116)

This result is indeed very close to the exact value of \(3\sqrt{\pi}\Gamma(\alpha+1/2)/4\Gamma(\alpha+3)\).

>>> (3*sqrt(pi)/4) * gamma(alpha+0.5) / gamma(alpha+3)
0.2926542074736711

In general, Gauss-Gegenbauer quadrature will only yield an approximate value of the integral. Consider the integral from -1 to 1 of \(\cos(x)(1-x^2)^{\alpha-1/2}\) which evaluates to \(2^\alpha\sqrt{\pi}\Gamma(\alpha+1/2)J_\alpha(1)\) where \(J_\alpha\) is the Bessel function of first kind and order \(\alpha\).

>>> import numpy as np
>>> weights @ np.cos(roots)
np.float64(1.5395778712347201)
>>> from scipy.special import jv
>>> 2**alpha * sqrt(pi) * gamma(alpha+0.5) *jv(alpha, 1)
np.float64(1.5395778706293377)