`special.orthogonal`¶

A collection of functions to find the weights and abscissas for Gaussian Quadrature.

These calculations are done by finding the eigenvalues of a tridiagonal matrix whose entries are dependent on the coefficients in the recursion formula for the orthogonal polynomials with the corresponding weighting function over the interval.

Many recursion relations for orthogonal polynomials are given:

The recursion relation of interest is

where has a different normalization than .

The coefficients can be found as:

where

assume:

For the mathematical background, see [golub.welsch-1969-mathcomp] and [abramowitz.stegun-1965].

References¶

 [golub.welsch-1969-mathcomp] Golub, Gene H, and John H Welsch. 1969. Calculation of Gauss Quadrature Rules. Mathematics of Computation 23, 221-230+s1–s10.
 [abramowitz.stegun-1965] Abramowitz, Milton, and Irene A Stegun. (1965) Handbook of Mathematical Functions: with Formulas, Graphs, and Mathematical Tables. Gaithersburg, MD: National Bureau of Standards. http://www.math.sfu.ca/~cbm/aands/
 [townsend.trogdon.olver-2014] Townsend, A. and Trogdon, T. and Olver, S. (2014) Fast computation of Gauss quadrature nodes and weights on the whole real line. :arXiv:`1410.5286`.
 [townsend.trogdon.olver-2015] Townsend, A. and Trogdon, T. and Olver, S. (2015) Fast computation of Gauss quadrature nodes and weights on the whole real line. IMA Journal of Numerical Analysis :doi:`10.1093/imanum/drv002`.

Module Contents¶

Classes¶

 `orthopoly1d`(self,roots,weights=None,hn=1.0,kn=1.0,wfunc=None,limits=None,monic=False,eval_func=None)

Functions¶

 `_gen_roots_and_weights`(n,mu0,an_func,bn_func,f,df,symmetrize,mu) [x,w] = gen_roots_and_weights(n,an_func,sqrt_bn_func,mu) `roots_jacobi`(n,alpha,beta,mu=False) rGauss-Jacobi quadrature. `jacobi`(n,alpha,beta,monic=False) rJacobi polynomial. `roots_sh_jacobi`(n,p1,q1,mu=False) Gauss-Jacobi (shifted) quadrature. `sh_jacobi`(n,p,q,monic=False) rShifted Jacobi polynomial. `roots_genlaguerre`(n,alpha,mu=False) rGauss-generalized Laguerre quadrature. `genlaguerre`(n,alpha,monic=False) rGeneralized (associated) Laguerre polynomial. `roots_laguerre`(n,mu=False) rGauss-Laguerre quadrature. `laguerre`(n,monic=False) rLaguerre polynomial. `roots_hermite`(n,mu=False) rGauss-Hermite (physicst’s) quadrature. `_compute_tauk`(n,k,maxit=5) Helper function for Tricomi initial guesses `_initial_nodes_a`(n,k) rTricomi initial guesses `_initial_nodes_b`(n,k) rGatteschi initial guesses `_initial_nodes`(n) Initial guesses for the Hermite roots `_pbcf`(n,theta) rAsymptotic series expansion of parabolic cylinder function `_newton`(n,x_initial,maxit=5) Newton iteration for polishing the asymptotic approximation `_roots_hermite_asy`(n) rGauss-Hermite (physicst’s) quadrature for large n. `hermite`(n,monic=False) rPhysicist’s Hermite polynomial. `roots_hermitenorm`(n,mu=False) rGauss-Hermite (statistician’s) quadrature. `hermitenorm`(n,monic=False) rNormalized (probabilist’s) Hermite polynomial. `roots_gegenbauer`(n,alpha,mu=False) rGauss-Gegenbauer quadrature. `gegenbauer`(n,alpha,monic=False) rGegenbauer (ultraspherical) polynomial. `roots_chebyt`(n,mu=False) rGauss-Chebyshev (first kind) quadrature. `chebyt`(n,monic=False) rChebyshev polynomial of the first kind. `roots_chebyu`(n,mu=False) rGauss-Chebyshev (second kind) quadrature. `chebyu`(n,monic=False) rChebyshev polynomial of the second kind. `roots_chebyc`(n,mu=False) rGauss-Chebyshev (first kind) quadrature. `chebyc`(n,monic=False) rChebyshev polynomial of the first kind on . `roots_chebys`(n,mu=False) rGauss-Chebyshev (second kind) quadrature. `chebys`(n,monic=False) rChebyshev polynomial of the second kind on . `roots_sh_chebyt`(n,mu=False) rGauss-Chebyshev (first kind, shifted) quadrature. `sh_chebyt`(n,monic=False) rShifted Chebyshev polynomial of the first kind. `roots_sh_chebyu`(n,mu=False) rGauss-Chebyshev (second kind, shifted) quadrature. `sh_chebyu`(n,monic=False) rShifted Chebyshev polynomial of the second kind. `roots_legendre`(n,mu=False) rGauss-Legendre quadrature. `legendre`(n,monic=False) rLegendre polynomial. `roots_sh_legendre`(n,mu=False) rGauss-Legendre (shifted) quadrature. `sh_legendre`(n,monic=False) rShifted Legendre polynomial.
class `orthopoly1d`(roots, weights=None, hn=1.0, kn=1.0, wfunc=None, limits=None, monic=False, eval_func=None)
`__init__`(roots, weights=None, hn=1.0, kn=1.0, wfunc=None, limits=None, monic=False, eval_func=None)
`__call__`(v)
`_scale`(p)
`_gen_roots_and_weights`(n, mu0, an_func, bn_func, f, df, symmetrize, mu)

[x,w] = gen_roots_and_weights(n,an_func,sqrt_bn_func,mu)

Returns the roots (x) of an nth order orthogonal polynomial, and weights (w) to use in appropriate Gaussian quadrature with that orthogonal polynomial.

The polynomials have the recurrence relation
P_n+1(x) = (x - A_n) P_n(x) - B_n P_n-1(x)

an_func(n) should return A_n sqrt_bn_func(n) should return sqrt(B_n) mu ( = h_0 ) is the integral of the weight over the orthogonal

interval
`roots_jacobi`(n, alpha, beta, mu=False)

Computes the sample points and weights for Gauss-Jacobi quadrature. The sample points are the roots of the n-th degree Jacobi polynomial, . These sample points and weights correctly integrate polynomials of degree or less over the interval with weight function .

n : int
alpha : float
alpha must be > -1
beta : float
beta must be > 0
mu : bool, optional
If True, return the sum of the weights, optional.
x : ndarray
Sample points
w : ndarray
Weights
mu : float
Sum of the weights

`jacobi`(n, alpha, beta, monic=False)

rJacobi polynomial.

Defined to be the solution of

for ; is a polynomial of degree .

n : int
Degree of the polynomial.
alpha : float
Parameter, must be greater than -1.
beta : float
Parameter, must be greater than -1.
monic : bool, optional
If True, scale the leading coefficient to be 1. Default is False.
P : orthopoly1d
Jacobi polynomial.

For fixed , the polynomials are orthogonal over with weight function .

`roots_sh_jacobi`(n, p1, q1, mu=False)

Computes the sample points and weights for Gauss-Jacobi (shifted) quadrature. The sample points are the roots of the n-th degree shifted Jacobi polynomial, . These sample points and weights correctly integrate polynomials of degree or less over the interval with weight function

n : int
p1 : float
(p1 - q1) must be > -1
q1 : float
q1 must be > 0
mu : bool, optional
If True, return the sum of the weights, optional.
x : ndarray
Sample points
w : ndarray
Weights
mu : float
Sum of the weights

`sh_jacobi`(n, p, q, monic=False)

rShifted Jacobi polynomial.

Defined by

where is the nth Jacobi polynomial.

n : int
Degree of the polynomial.
p : float
Parameter, must have .
q : float
Parameter, must be greater than 0.
monic : bool, optional
If True, scale the leading coefficient to be 1. Default is False.
G : orthopoly1d
Shifted Jacobi polynomial.

For fixed , the polynomials are orthogonal over with weight function .

`roots_genlaguerre`(n, alpha, mu=False)

Computes the sample points and weights for Gauss-generalized Laguerre quadrature. The sample points are the roots of the n-th degree generalized Laguerre polynomial, . These sample points and weights correctly integrate polynomials of degree or less over the interval with weight function .

n : int
alpha : float
alpha must be > -1
mu : bool, optional
If True, return the sum of the weights, optional.
x : ndarray
Sample points
w : ndarray
Weights
mu : float
Sum of the weights

`genlaguerre`(n, alpha, monic=False)

rGeneralized (associated) Laguerre polynomial.

Defined to be the solution of

where ; is a polynomial of degree .

n : int
Degree of the polynomial.
alpha : float
Parameter, must be greater than -1.
monic : bool, optional
If True, scale the leading coefficient to be 1. Default is False.
L : orthopoly1d
Generalized Laguerre polynomial.

For fixed , the polynomials are orthogonal over with weight function .

The Laguerre polynomials are the special case where .

laguerre : Laguerre polynomial.

`roots_laguerre`(n, mu=False)

Computes the sample points and weights for Gauss-Laguerre quadrature. The sample points are the roots of the n-th degree Laguerre polynomial, . These sample points and weights correctly integrate polynomials of degree or less over the interval with weight function .

n : int
mu : bool, optional
If True, return the sum of the weights, optional.
x : ndarray
Sample points
w : ndarray
Weights
mu : float
Sum of the weights

`laguerre`(n, monic=False)

rLaguerre polynomial.

Defined to be the solution of

is a polynomial of degree .

n : int
Degree of the polynomial.
monic : bool, optional
If True, scale the leading coefficient to be 1. Default is False.
L : orthopoly1d
Laguerre Polynomial.

The polynomials are orthogonal over with weight function .

`roots_hermite`(n, mu=False)

Computes the sample points and weights for Gauss-Hermite quadrature. The sample points are the roots of the n-th degree Hermite polynomial, . These sample points and weights correctly integrate polynomials of degree or less over the interval with weight function .

n : int
mu : bool, optional
If True, return the sum of the weights, optional.
x : ndarray
Sample points
w : ndarray
Weights
mu : float
Sum of the weights

For small n up to 150 a modified version of the Golub-Welsch algorithm is used. Nodes are computed from the eigenvalue problem and improved by one step of a Newton iteration. The weights are computed from the well-known analytical formula.

For n larger than 150 an optimal asymptotic algorithm is applied which computes nodes and weights in a numerically stable manner. The algorithm has linear runtime making computation for very large n (several thousand or more) feasible.

 [townsend.trogdon.olver-2014] Townsend, A. and Trogdon, T. and Olver, S. (2014) Fast computation of Gauss quadrature nodes and weights on the whole real line. :arXiv:`1410.5286`.
 [townsend.trogdon.olver-2015] Townsend, A. and Trogdon, T. and Olver, S. (2015) Fast computation of Gauss quadrature nodes and weights on the whole real line. IMA Journal of Numerical Analysis :doi:`10.1093/imanum/drv002`.
`_compute_tauk`(n, k, maxit=5)

Helper function for Tricomi initial guesses

For details, see formula 3.1 in lemma 3.1 in the original paper.

n : int
k : ndarray of type int
Index of roots to compute
maxit : int
Number of Newton maxit performed, the default value of 5 is sufficient.
tauk : ndarray
Roots of equation 3.1

initial_nodes_a roots_hermite_asy

`_initial_nodes_a`(n, k)

rTricomi initial guesses

Computes an initial approximation to the square of the k-th (positive) root of the Hermite polynomial of order . The formula is the one from lemma 3.1 in the original paper. The guesses are accurate except in the region near .

n : int
k : ndarray of type int
Index of roots to compute
xksq : ndarray
Square of the approximate roots

initial_nodes roots_hermite_asy

`_initial_nodes_b`(n, k)

rGatteschi initial guesses

Computes an initial approximation to the square of the k-th (positive) root of the Hermite polynomial of order . The formula is the one from lemma 3.2 in the original paper. The guesses are accurate in the region just below .

n : int
k : ndarray of type int
Index of roots to compute
xksq : ndarray
Square of the approximate root

initial_nodes roots_hermite_asy

`_initial_nodes`(n)

Initial guesses for the Hermite roots

Computes an initial approximation to the non-negative roots of the Hermite polynomial of order . The Tricomi and Gatteschi initial guesses are used in the region where they are accurate.

n : int
xk : ndarray
Approximate roots

roots_hermite_asy

`_pbcf`(n, theta)

rAsymptotic series expansion of parabolic cylinder function

The implementation is based on sections 3.2 and 3.3 from the original paper. Compared to the published version this code adds one more term to the asymptotic series. The detailed formulas can be found at [parabolic-asymptotics]. The evaluation is done in a transformed variable where and .

n : int
theta : ndarray
Transformed position variable
U : ndarray
Value of the parabolic cylinder function .
Ud : ndarray
Value of the derivative of the parabolic cylinder function.

roots_hermite_asy

`_newton`(n, x_initial, maxit=5)

Newton iteration for polishing the asymptotic approximation to the zeros of the Hermite polynomials.

n : int
x_initial : ndarray
Initial guesses for the roots
maxit : int
Maximal number of Newton iterations. The default 5 is sufficient, usually only one or two steps are needed.
nodes : ndarray
weights : ndarray

roots_hermite_asy

`_roots_hermite_asy`(n)

rGauss-Hermite (physicst’s) quadrature for large n.

Computes the sample points and weights for Gauss-Hermite quadrature. The sample points are the roots of the n-th degree Hermite polynomial, . These sample points and weights correctly integrate polynomials of degree or less over the interval with weight function .

This method relies on asymptotic expansions which work best for n > 150. The algorithm has linear runtime making computation for very large n feasible.

n : int
nodes : ndarray
weights : ndarray

roots_hermite

 [townsend.trogdon.olver-2014] Townsend, A. and Trogdon, T. and Olver, S. (2014) Fast computation of Gauss quadrature nodes and weights on the whole real line. :arXiv:`1410.5286`.
 [townsend.trogdon.olver-2015] Townsend, A. and Trogdon, T. and Olver, S. (2015) Fast computation of Gauss quadrature nodes and weights on the whole real line. IMA Journal of Numerical Analysis :doi:`10.1093/imanum/drv002`.
`hermite`(n, monic=False)

rPhysicist’s Hermite polynomial.

Defined by

is a polynomial of degree .

n : int
Degree of the polynomial.
monic : bool, optional
If True, scale the leading coefficient to be 1. Default is False.
H : orthopoly1d
Hermite polynomial.

The polynomials are orthogonal over with weight function .

`roots_hermitenorm`(n, mu=False)

Computes the sample points and weights for Gauss-Hermite quadrature. The sample points are the roots of the n-th degree Hermite polynomial, . These sample points and weights correctly integrate polynomials of degree or less over the interval with weight function .

n : int
mu : bool, optional
If True, return the sum of the weights, optional.
x : ndarray
Sample points
w : ndarray
Weights
mu : float
Sum of the weights

For small n up to 150 a modified version of the Golub-Welsch algorithm is used. Nodes are computed from the eigenvalue problem and improved by one step of a Newton iteration. The weights are computed from the well-known analytical formula.

For n larger than 150 an optimal asymptotic algorithm is used which computes nodes and weights in a numerical stable manner. The algorithm has linear runtime making computation for very large n (several thousand or more) feasible.

`hermitenorm`(n, monic=False)

rNormalized (probabilist’s) Hermite polynomial.

Defined by

is a polynomial of degree .

n : int
Degree of the polynomial.
monic : bool, optional
If True, scale the leading coefficient to be 1. Default is False.
He : orthopoly1d
Hermite polynomial.

The polynomials are orthogonal over with weight function .

`roots_gegenbauer`(n, alpha, mu=False)

Computes the sample points and weights for Gauss-Gegenbauer quadrature. The sample points are the roots of the n-th degree Gegenbauer polynomial, . These sample points and weights correctly integrate polynomials of degree or less over the interval with weight function .

n : int
alpha : float
alpha must be > -0.5
mu : bool, optional
If True, return the sum of the weights, optional.
x : ndarray
Sample points
w : ndarray
Weights
mu : float
Sum of the weights

`gegenbauer`(n, alpha, monic=False)

rGegenbauer (ultraspherical) polynomial.

Defined to be the solution of

for ; is a polynomial of degree .

n : int
Degree of the polynomial.
monic : bool, optional
If True, scale the leading coefficient to be 1. Default is False.
C : orthopoly1d
Gegenbauer polynomial.

The polynomials are orthogonal over with weight function .

`roots_chebyt`(n, mu=False)

rGauss-Chebyshev (first kind) quadrature.

Computes the sample points and weights for Gauss-Chebyshev quadrature. The sample points are the roots of the n-th degree Chebyshev polynomial of the first kind, . These sample points and weights correctly integrate polynomials of degree or less over the interval with weight function .

n : int
mu : bool, optional
If True, return the sum of the weights, optional.
x : ndarray
Sample points
w : ndarray
Weights
mu : float
Sum of the weights

`chebyt`(n, monic=False)

rChebyshev polynomial of the first kind.

Defined to be the solution of

is a polynomial of degree .

n : int
Degree of the polynomial.
monic : bool, optional
If True, scale the leading coefficient to be 1. Default is False.
T : orthopoly1d
Chebyshev polynomial of the first kind.

The polynomials are orthogonal over with weight function .

chebyu : Chebyshev polynomial of the second kind.

`roots_chebyu`(n, mu=False)

rGauss-Chebyshev (second kind) quadrature.

Computes the sample points and weights for Gauss-Chebyshev quadrature. The sample points are the roots of the n-th degree Chebyshev polynomial of the second kind, . These sample points and weights correctly integrate polynomials of degree or less over the interval with weight function .

n : int
mu : bool, optional
If True, return the sum of the weights, optional.
x : ndarray
Sample points
w : ndarray
Weights
mu : float
Sum of the weights

`chebyu`(n, monic=False)

rChebyshev polynomial of the second kind.

Defined to be the solution of

is a polynomial of degree .

n : int
Degree of the polynomial.
monic : bool, optional
If True, scale the leading coefficient to be 1. Default is False.
U : orthopoly1d
Chebyshev polynomial of the second kind.

The polynomials are orthogonal over with weight function .

chebyt : Chebyshev polynomial of the first kind.

`roots_chebyc`(n, mu=False)

rGauss-Chebyshev (first kind) quadrature.

Computes the sample points and weights for Gauss-Chebyshev quadrature. The sample points are the roots of the n-th degree Chebyshev polynomial of the first kind, . These sample points and weights correctly integrate polynomials of degree or less over the interval with weight function .

n : int
mu : bool, optional
If True, return the sum of the weights, optional.
x : ndarray
Sample points
w : ndarray
Weights
mu : float
Sum of the weights

`chebyc`(n, monic=False)

rChebyshev polynomial of the first kind on .

Defined as , where is the nth Chebychev polynomial of the first kind.

n : int
Degree of the polynomial.
monic : bool, optional
If True, scale the leading coefficient to be 1. Default is False.
C : orthopoly1d
Chebyshev polynomial of the first kind on .

The polynomials are orthogonal over with weight function .

chebyt : Chebyshev polynomial of the first kind.

 [1] Abramowitz and Stegun, “Handbook of Mathematical Functions” Section 22. National Bureau of Standards, 1972.
`roots_chebys`(n, mu=False)

rGauss-Chebyshev (second kind) quadrature.

Computes the sample points and weights for Gauss-Chebyshev quadrature. The sample points are the roots of the n-th degree Chebyshev polynomial of the second kind, . These sample points and weights correctly integrate polynomials of degree or less over the interval with weight function .

n : int
mu : bool, optional
If True, return the sum of the weights, optional.
x : ndarray
Sample points
w : ndarray
Weights
mu : float
Sum of the weights

`chebys`(n, monic=False)

rChebyshev polynomial of the second kind on .

Defined as where is the nth Chebychev polynomial of the second kind.

n : int
Degree of the polynomial.
monic : bool, optional
If True, scale the leading coefficient to be 1. Default is False.
S : orthopoly1d
Chebyshev polynomial of the second kind on .

The polynomials are orthogonal over with weight function .

chebyu : Chebyshev polynomial of the second kind

 [1] Abramowitz and Stegun, “Handbook of Mathematical Functions” Section 22. National Bureau of Standards, 1972.
`roots_sh_chebyt`(n, mu=False)

rGauss-Chebyshev (first kind, shifted) quadrature.

Computes the sample points and weights for Gauss-Chebyshev quadrature. The sample points are the roots of the n-th degree shifted Chebyshev polynomial of the first kind, . These sample points and weights correctly integrate polynomials of degree or less over the interval with weight function .

n : int
mu : bool, optional
If True, return the sum of the weights, optional.
x : ndarray
Sample points
w : ndarray
Weights
mu : float
Sum of the weights

`sh_chebyt`(n, monic=False)

rShifted Chebyshev polynomial of the first kind.

Defined as for the nth Chebyshev polynomial of the first kind.

n : int
Degree of the polynomial.
monic : bool, optional
If True, scale the leading coefficient to be 1. Default is False.
T : orthopoly1d
Shifted Chebyshev polynomial of the first kind.

The polynomials are orthogonal over with weight function .

`roots_sh_chebyu`(n, mu=False)

rGauss-Chebyshev (second kind, shifted) quadrature.

Computes the sample points and weights for Gauss-Chebyshev quadrature. The sample points are the roots of the n-th degree shifted Chebyshev polynomial of the second kind, . These sample points and weights correctly integrate polynomials of degree or less over the interval with weight function .

n : int
mu : bool, optional
If True, return the sum of the weights, optional.
x : ndarray
Sample points
w : ndarray
Weights
mu : float
Sum of the weights

`sh_chebyu`(n, monic=False)

rShifted Chebyshev polynomial of the second kind.

Defined as for the nth Chebyshev polynomial of the second kind.

n : int
Degree of the polynomial.
monic : bool, optional
If True, scale the leading coefficient to be 1. Default is False.
U : orthopoly1d
Shifted Chebyshev polynomial of the second kind.

The polynomials are orthogonal over with weight function .

`roots_legendre`(n, mu=False)

Computes the sample points and weights for Gauss-Legendre quadrature. The sample points are the roots of the n-th degree Legendre polynomial . These sample points and weights correctly integrate polynomials of degree or less over the interval with weight function .

n : int
mu : bool, optional
If True, return the sum of the weights, optional.
x : ndarray
Sample points
w : ndarray
Weights
mu : float
Sum of the weights

`legendre`(n, monic=False)

rLegendre polynomial.

Defined to be the solution of

is a polynomial of degree .

n : int
Degree of the polynomial.
monic : bool, optional
If True, scale the leading coefficient to be 1. Default is False.
P : orthopoly1d
Legendre polynomial.

The polynomials are orthogonal over with weight function 1.

Generate the 3rd-order Legendre polynomial 1/2*(5x^3 + 0x^2 - 3x + 0):

```>>> from scipy.special import legendre
>>> legendre(3)
poly1d([ 2.5,  0. , -1.5,  0. ])
```
`roots_sh_legendre`(n, mu=False)

Computes the sample points and weights for Gauss-Legendre quadrature. The sample points are the roots of the n-th degree shifted Legendre polynomial . These sample points and weights correctly integrate polynomials of degree or less over the interval with weight function .

n : int
mu : bool, optional
If True, return the sum of the weights, optional.
x : ndarray
Sample points
w : ndarray
Weights
mu : float
Sum of the weights

`sh_legendre`(n, monic=False)

rShifted Legendre polynomial.

Defined as for the nth Legendre polynomial.

n : int
Degree of the polynomial.
monic : bool, optional
If True, scale the leading coefficient to be 1. Default is False.
P : orthopoly1d
Shifted Legendre polynomial.

The polynomials are orthogonal over with weight function 1.