Monday, November 6, 2017

Python: Orthogonal Polynomials and Generalized Gauss Quadrature

A new (to me) python library for easily computing families of orthogonal polynomials.

Getting standard (generalized) Gauss quadrature schemes is extremely simple. For example to get 13 nodes and weights for Gauss-Laguerre integration, correct up to 50 decimal places:

pts,wts = orthopy.schemes.laguerre(13, decimal_places=50)

The numpy Polynomial package provides similar functionality:

pts, wts = numpy.polynomial.laguerre.laggauss(13)

A nice feature (besides arbitrary precision) is that you can derive custom orthogonal polynomials and quadrature rules. All you need to provide is a weight function and domain of the polynomials. From the project webpage:
import orthopy
moments = orthopy.compute_moments(lambda x: x**2, -1, +1, 20)
alpha, beta = orthopy.chebyshev(moments)
points, weights = orthopy.schemes.custom(alpha, beta, decimal_places=30)
This generates a 10-point scheme for integrating functions over the interval [-1, 1], with weight function \(w(x) = x^2\).

No comments: