gaussQuadRule¶

Computes a Gauss, Gauss-Radau, or Gauss-Lobatto quadrature rule with various classical weight functions.

Synopsis¶

gaussQuadRule (quadpts, weights, points)

Required Arguments¶

int quadpts (Input): Number of quadrature points.
float weights[] (Output): Array of length n containing the quadrature weights.
float points[] (Output): Array of length n containing quadrature points. The default action of this routine is to produce the Gauss Legendre points and weights.

Optional Arguments¶

chebyshevFirst

Compute the Gauss points and weights using the weight function

\[1 / \sqrt{1 - x^2}\]

on the interval (−1, 1).

chebyshevSecond

Compute the Gauss points and weights using the weight function

\[\sqrt{1-x^2}\]

on the interval (−1, 1).

hermite

Compute the Gauss points and weights using the weight function \(\exp(-x^2\)) on the interval (−∞, ∞).

cosh

Compute the Gauss points and weights using the weight function \(1 ∕ (\cosh(x))\) on the interval (−∞, ∞).

jacobi, float alpha, float beta (Input)

Compute the Gauss points and weights using the weight function \((1 − x)^a (1 + x)^b\) on the interval (−1, 1).

genLaguerre, float (Input)

Compute the Gauss points and weights using the weight function \(\exp(-x)x^a\) on the interval (0, ∞).

fixedPoint, float (Input)

Compute the Gauss-Radau points and weights using the specified weight function and the fixed point a. This formula will integrate polynomials of degree less than \(2n − 1\) exactly.

twoFixedPoints, float a, float b (Input)

Compute the Gauss-Lobatto points and weights using the specified weight function and the fixed points a and b. This formula will integrate polynomials of degree less than \(2n − 2\) exactly.

Description¶

The function gaussQuadRule produces the points and weights for the Gauss, Gauss-Radau, or Gauss-Lobatto quadrature formulas for some of the most popular weights. The default weight is the weight function identically equal to 1 on the interval (−1, 1). In fact, it is slightly more general than this suggests, because the extra one or two points that may be specified do not have to lie at the endpoints of the interval. This function is a modification of the subroutine GAUSSQUADRULE (Golub and Welsch 1969).

In the default case, the function returns points in x = points and weights in w = weights so that

\[\int_a^b f(x)w(x)dx = \sum_{i=1}^{N} f\left(x_i\right)w_i\]

for all functions f that are polynomials of degree less than 2n.

If the keyword fixedPoint is specified, then one of the above \(x_i\) is equal to a. Similarly, if the keyword twoFixedPoints is specified, then two of the components of x are equal to a and b. In general, the accuracy of the above quadrature formula degrades when n increases. The quadrature rule will integrate all functions f that are polynomials of degree less than 2n − F, where F is the number of fixed points.

Examples¶

Example 1¶

The three-point Gauss Legendre quadrature points and weights are computed and used to approximate the integrals

\[\int_{-1}^{1} x^i dx \phantom{...} i = 0, \ldots, 6\]

Notice that the integrals are exact for the first six monomials, but that the last approximation is in error. In general, the Gauss rules with k points integrate polynomials with degree less than 2k exactly.

from __future__ import print_function
from numpy import *
from pyimsl.math.gaussQuadRule import gaussQuadRule

QUADPTS = 3
POWERS = 7
weights = []
points = []
s = zeros(POWERS, dtype=double)
gaussQuadRule(QUADPTS, weights, points)
for i in range(0, POWERS, 1):
    s[i] = 0.0
    for j in range(0, QUADPTS, 1):
        s[i] = s[i] + weights[j] * pow(points[j], i)

print("The integral from -1 to 1 of pow(x, i) is")
print("Function            Quadrature    Exact\n")
z = 1.0
for i in range(0, POWERS, 1):
    z = (1 - i % 2) * 2. / (i + 1.)
    print("pow(x, %d)        %10.3f  %10.3f" % (i, s[i], z))

Output¶

The integral from -1 to 1 of pow(x, i) is
Function            Quadrature    Exact

pow(x, 0)             2.000       2.000
pow(x, 1)            -0.000       0.000
pow(x, 2)             0.667       0.667
pow(x, 3)            -0.000       0.000
pow(x, 4)             0.400       0.400
pow(x, 5)            -0.000       0.000
pow(x, 6)             0.240       0.286

Example 2¶

The three-point Gauss Laguerre quadrature points and weights are computed and used to approximate the integrals

\[\int_{0}^{\infty} x^i xe^{-x} dx = i! \phantom{...} i = 0, \ldots, 6\]

Notice that the integrals are exact for the first six monomials, but that the last approximation is in error. In general, the Gauss rules with k points integrate polynomials with degree less than 2k exactly.

from __future__ import print_function
from numpy import *
from pyimsl.math.gaussQuadRule import gaussQuadRule

QUADPTS = 3
POWERS = 7
weights = []
points = []
s = zeros(POWERS, dtype=double)
gaussQuadRule(QUADPTS, weights, points, genLaguerre=1.0)

for i in range(0, POWERS, 1):
    s[i] = 0.0
    for j in range(0, QUADPTS, 1):
        s[i] = s[i] + weights[j] * pow(points[j], i)

print("The integral from 0 to infinity of pow(x, i)*x*exp(x) is")
print("Function            Quadrature    Exact\n")
z = 1.0
for i in range(0, POWERS, 1):
    z = z * (i + 1)
    print("pow(x, %d)        %10.3f  %10.3f" % (i, s[i], z))

Output¶

The integral from 0 to infinity of pow(x, i)*x*exp(x) is
Function            Quadrature    Exact

pow(x, 0)             1.000       1.000
pow(x, 1)             2.000       2.000
pow(x, 2)             6.000       6.000
pow(x, 3)            24.000      24.000
pow(x, 4)           120.000     120.000
pow(x, 5)           720.000     720.000
pow(x, 6)          4896.000    5040.000