#!/usr/bin/env python """ 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: a1n f_n+1 (x) = (a2n + a3n x ) f_n (x) - a4n f_n-1 (x) The recursion relation of interest is P_n+1 (x) = (x - A_n) P_n (x) - B_n P_n-1 (x) where P has a different normalization than f. The coefficients can be found as: A_n = -a2n / a3n B_n = ( a4n / a3n sqrt(h_n-1 / h_n))**2 where h_n = int_a^b w(x) f_n(x)^2 See Numerical Recipies in C, page 156 and Abramowitz and Stegun p. 774, 782 """ from MLab import eig from MLab import * import cephes global p, q def gen_roots_and_weights(n,an_func,sqrt_bn_func,mu): """[x,w] = gen_roots_and_weights(n,an_func,sqrt_bn_func,mu) Returns the roots (x) of an nth order orthoganl polynomail, 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 """ global p,q nn = arange(1.0,n) sqrt_bn = sqrt_bn_func(nn) an = an_func(concatenate(([0],nn))) [x,v] = eig((diag(an)+diag(sqrt_bn,1)+diag(sqrt_bn,-1))) answer = [] sortind = argsort(x) answer.append(take(x,sortind)) answer.append(take(mu*v[:,0]**2,sortind)) return answer # Jacobi Polynomials 1 P^(alpha,beta)_n(x) def J_roots(n,alpha,beta): """[x,w] = J_roots(n,alpha,beta) Returns the roots (x) of the nth order Jacobi polynomial, P^(alpha,beta)_n(x) and weights (w) to use in Gaussian Quadrature over [-1,1] with weighting function (1-x)**alpha (1+x)**beta with alpha,beta > -1. """ if (alpha <= -1) or (beta <= -1): raise ValueError, "alpha and beta must be greater than -1." global p,q (p,q) = (alpha,beta) # from recurrence relation sbn_J = lambda k: 2.0/(2*k+p+q)*sqrt((k+p)*(k+q)*k*(k+p+q)/((2*k+p+q+1)*(2*k+p+q-1))) if (p == q): an_J = lambda k: 0*k else: an_J = lambda k: (q*q - p*p)/((2.0*k+p+q)*(2.0*k+p+q+2)) g = cephes.gamma mu0 = 2**(p+q+1)*g(p+1)*g(q+1)/(g(p+q+2)) return gen_roots_and_weights(n,an_J,sbn_J,mu0) # Jacobi Polynomials shifted G_n(p,q,x) def Js_roots(n,p1,q1): """[x,w] = Js_roots(n,p,q) Returns the roots (x) of the nth order shifted Jacobi polynomial, G_n(p,q,x), and weights (w) to use in Gaussian Quadrature over [0,1] with weighting function (1-x)**(p-q) x**(q-1) with p-q > -1 and q > 0. """ # from recurrence relation if not ( ( (p1 - q1) > -1 ) and ( q1 > 0 ) ): raise ValueError, "(p - q) > -1 and q > 0 please." global p,q (p,q) = (p1,q1) if (p == 0): sbn_Js = lambda k: sqrt((k+q-1)*(k-q))/((2*k-1)*2) else: sbn_Js = lambda k: sqrt((k*(k+q-1)*(k+p-1)*(k+p-q))/((2*k+p)*(2*k+p-2)))/(2*k+p-1) if (p == 1): an_Js = lambda k: 0*k + 0.5 else: an_Js = lambda k: (2*k*(k+p)+q*(p-1))/((2*k+p+1)*(2*k+p-1)) # integral of weight over interval g = cephes.gamma mu0 = g(q)*g(p-q+1)/g(p+1) return gen_roots_and_weights(n,an_Js,sbn_Js,mu0) # Generalized Laguerre L^(alpha)_n(x) def La_roots(n,alpha): """[x,w] = La_roots(n,alpha) Returns the roots (x) of the nth order generalized Laguerre polynomial, L^(alpha)(x,alpha), and weights (w) to use in Gaussian Quadrature over [0,inf] with weighting function exp(-x) x**alpha with alpha > -1. """ if not (alpha > -1): raise ValueError, "alpha > -1 please." global p,q (p,q) = (alpha,0.0) sbn_La = lambda k: -sqrt(k*(k + p)) # from recurrence relation an_La = lambda k: 2*k + p + 1 mu0 = cephes.gamma(alpha+1) # integral of weight over interval return gen_roots_and_weights(n,an_La,sbn_La,mu0) # Hermite 1 H_n(x) def H_roots(n): """[x,w] = H_roots(n) Returns the roots (x) of the nth order Hermite polynomial, H_n(x), and weights (w) to use in Gaussian Quadrature over [-inf,inf] with weighting function exp(-x**2). """ sbn_H = lambda k: sqrt(k/2) # from recurrence relation an_H = lambda k: 0*k mu0 = sqrt(pi) # integral of weight over interval return gen_roots_and_weights(n,an_H,sbn_H,mu0) # Hermite 2 He_n(x) def He_roots(n): """[x,w] = He_roots(n) Returns the roots (x) of the nth order Hermite polynomial, He_n(x), and weights (w) to use in Gaussian Quadrature over [-inf,inf] with weighting function exp(-(x/2)**2). """ sbn_He = lambda k: sqrt(k) # from recurrence relation an_He = lambda k: 0*k mu0 = sqrt(2*pi) # integral of weight over interval return gen_roots_and_weights(n,an_He,sbn_He,mu0) ## The remainder of the polynomials can be derived from the ones above. # Ultraspherical (Gegenbauer) C^(alpha)_n(x) def Cg_roots(n,alpha): """[x,w] = Cg_roots(n,alpha) Returns the roots (x) of the nth order Ultraspherical (Gegenbauer) polynomial, C^(alpha)_n(x), and weights (w) to use in Gaussian Quadrature over [-1,1] with weighting function (1-x**2)**(alpha-1/2) with alpha>-1/2. """ return J_roots(n,alpha-0.5,alpha-0.5) # Chebyshev of the first kind T_n(x) def T_roots(n): """[x,w] = T_roots(n) Returns the roots (x) of the nth order Chebyshev (of the first kind) polynomial, T_n(x), and weights (w) to use in Gaussian Quadrature over [-1,1] with weighting function (1-x**2)**(-1/2). """ return J_roots(n,-0.5,-0.5) # Chebyshev of the second kind U_n(x) def U_roots(n): """[x,w] = U_roots(n) Returns the roots (x) of the nth order Chebyshev (of the second kind) polynomial, U_n(x), and weights (w) to use in Gaussian Quadrature over [-1,1] with weighting function (1-x**2)**1/2. """ return J_roots(n,0.5,0.5) # Chebyshev of the first kind C_n(x) def C_roots(n): """[x,w] = C_roots(n) Returns the roots (x) of the nth order Chebyshev (of the first kind) polynomial, C_n(x), and weights (w) to use in Gaussian Quadrature over [-2,2] with weighting function (1-(x/2)**2)**1/2. """ [x,w] = J_roots(n,-0.5,-0.5) return [x*2,w] # Chebyshev of the second kind S_n(x) def S_roots(n): """[x,w] = S_roots(n) Returns the roots (x) of the nth order Chebyshev (of the second kind) polynomial, S_n(x), and weights (w) to use in Gaussian Quadrature over [-2,2] with weighting function (1-(x/2)**2)**1/2. """ [x,w] = J_roots(n,0.5,0.5) return [x*2,w] # Shifted Chebyshev of the first kind T^*_n(x) def Ts_roots(n): """[x,w] = Ts_roots(n) Returns the roots (x) of the nth order shifted Chebyshev (of the first kind) polynomial, T^*_n(x), and weights (w) to use in Gaussian Quadrature over [0,1] with weighting function (x-x**2)**(-1/2). """ return Js_roots(n,0.0,0.5) # Shifted Chebyshev of the second kind U^*_n(x) def Us_roots(n): """[x,w] = Us_roots(n) Returns the roots (x) of the nth order shifted Chebyshev (of the second kind) polynomial, U^*_n(x), and weights (w) to use in Gaussian Quadrature over [0,1] with weighting function (x-x**2)**1/2. """ return Js_roots(n,2.0,1.5) # Legendre def P_roots(n): """[x,w] = P_roots(n) Returns the roots (x) of the nth order Legendre polynomial, P_n(x), and weights (w) to use in Gaussian Quadrature over [-1,1] with weighting function 1. """ return J_roots(n,0.0,0.0) # Shifted Legendre P^*_n(x) def Ps_roots(n): """[x,w] = Ps_roots(n) Returns the roots (x) of the nth order shifted Legendre polynomial, P^*_n(x), and weights (w) to use in Gaussian Quadrature over [0,1] with weighting function 1. """ return Js_roots(n,1.0,1.0) # Laguerre L_n(x) def L_roots(n): """[x,w] = L_roots(n) Returns the roots (x) of the nth order Laguerre polynomial, L_n(x), and weights (w) to use in Gaussian Quadrature over [0,inf] with weighting function exp(-x). """ return La_roots(n,0.0)