# PATTERSON_RULE_COMPUTE Compute Gauss-Patterson Quadrature Rules

PATTERSON_RULE_COMPUTE is a FORTRAN90 program which computes the points and weights of a 1D Gauss-Patterson quadrature rule of order 1, 3, 7, 15, 31, 63, 127, 255 or 511, to approximate integrals in the interval [A,B], as specified by the user.

The rule is written to three files for easy use as input to other programs.

The program is based on Patterson's ACM TOMS algorithm 672, which is actually designed for the more general problem of creating a more accurate quadrature rule from any given one by adding intermediate points.

This program differs from the patterson_rule program in that it actually computes the points and weights directly, rather than simply looking them up in a table.

The Gauss-Patterson quadrature is a nested family which begins with the Gauss-Legendre rules of orders 1 and 3, and then succesively inserts one new abscissa in each subinterval. Thus, after the second rule, the Gauss-Patterson rules do not have the super-high precision of the Gauss-Legendre rules. They trade this precision in exchange for the advantages of nestedness. This means that Gauss-Patterson rules are only available for orders of 1, 3, 7, 15, 31, 63, 127, 255 or 511.

The standard Gauss-Patterson quadrature rule is used as follows:

Integral ( A <= x <= B ) f(x) dx

is to be approximated by
Sum ( 1 <= i <= order ) w(i) * f(x(i))

The polynomial precision of a Gauss-Patterson rule can be checked numerically by the INT_EXACTNESS_LEGENDRE program. We should expect
IndexOrderFree+FixedExpected PrecisionActual Precision
011 + 02*1+0-1=11
133 + 02*3+0-1=55
274 + 32*4+3-1=1010 + 1 = 11
3158 + 72*8+7-1=2222 + 1 = 23
43116 + 152*16+15-1=4646 + 1 = 47
56332 + 312*32+31-1=9494 + 1 = 95
612764 + 632*64+63-1=190190 + 1 = 191
7255128 + 1272*128+127-1=382382 + 1 = 383
8511256 + 2552*256+255-1=766766 + 1 = 767
where the extra 1 degree of precision comes about because the rules are symmetric, and can integrate any odd monomial exactly. Thus, after the first rule, the precision is 3*2^index - 1.

### Usage:

patterson_rule_compute order a b filename
where
• order is the number of points in the quadrature rule. Acceptable values are 1, 3, 7, 15, 31, 63, 127, 255 or 511.
• a is the left endpoint;
• b is the right endpoint;
• filename specifies the output filenames: filename_w.txt, filename_x.txt, and filename_r.txt, containing the weights, abscissas, and interval limits.

### Languages:

PATTERSON_RULE_COMPUTE is available in a FORTRAN90 version.

### Related Data and Programs:

CCN_RULE, a FORTRAN90 program which defines a nested Clenshaw Curtis quadrature rule.

CHEBYSHEV1_RULE, a FORTRAN90 program which can compute and print a Gauss-Chebyshev type 1 quadrature rule.

CHEBYSHEV2_RULE, a FORTRAN90 program which can compute and print a Gauss-Chebyshev type 2 quadrature rule.

CLENSHAW_CURTIS_RULE, a FORTRAN90 program which defines a Clenshaw Curtis quadrature rule.

GEGENBAUER_RULE, a FORTRAN90 program which can compute and print a Gauss-Gegenbauer quadrature rule.

GEN_HERMITE_RULE, a FORTRAN90 program which can compute and print a generalized Gauss-Hermite quadrature rule.

GEN_LAGUERRE_RULE, a FORTRAN90 program which can compute and print a generalized Gauss-Laguerre quadrature rule.

HERMITE_RULE, a FORTRAN90 program which can compute and print a Gauss-Hermite quadrature rule.

INT_EXACTNESS_LEGENDRE, a FORTRAN90 program which checks the polynomial exactness of a Gauss-Legendre quadrature rule.

JACOBI_RULE, a FORTRAN90 program which can compute and print a Gauss-Jacobi quadrature rule.

KRONROD, a FORTRAN90 library which can compute a Gauss and Gauss-Kronrod pair of quadrature rules of arbitrary order, by Robert Piessens, Maria Branders.

LAGUERRE_RULE, a FORTRAN90 program which can compute and print a Gauss-Laguerre quadrature rule.

LEGENDRE_RULE, a FORTRAN90 program which can compute and print a Gauss-Legendre quadrature rule.

LEGENDRE_RULE_FAST, a FORTRAN90 program which uses a fast (order N) algorithm to compute a Gauss-Legendre quadrature rule of given order.

LOGNORMAL_RULE, a FORTRAN90 program which can compute and print a quadrature rule for functions of a variable whose logarithm is normally distributed.

PATTERSON_RULE, a FORTRAN90 program which returns the points and weights of a 1D Gauss-Patterson quadrature rule of order 1, 3, 7, 15, 31, 63, 127, 255 or 511.

TOMS699, a FORTRAN77 library which implements a new representation of Patterson's quadrature formula;
this is ACM TOMS algorithm 699.

TRUNCATED_NORMAL_RULE, a FORTRAN90 program which computes a quadrature rule for a normal distribution that has been truncated to [A,+oo), (-oo,B] or [A,B].

### Reference:

1. Milton Abramowitz, Irene Stegun,
Handbook of Mathematical Functions,
National Bureau of Standards, 1964,
ISBN: 0-486-61272-4,
LC: QA47.A34.
2. Philip Davis, Philip Rabinowitz,
Methods of Numerical Integration,
Second Edition,
Dover, 2007,
ISBN: 0486453391,
LC: QA299.3.D28.
3. Gene Golub, Thomas Robertson,
A generalized Bairstow Algorithm,
Communications of the ACM,
Volume 10, Number 6, June 1967, pages 371-373.
4. Thomas Patterson,
Mathematics of Computation,
Volume 22, Number 104, October 1968, pages 847-856.
5. Thomas Patterson,
An algorithm for generating interpolatory quadrature rules of the highest degree of precision with preassigned nodes for general weight functions,
Transactions on Mathematical Software,
Volume 15, Number 2, June 1989, pages 123-136.
6. Thomas Patterson,
Algorithm 672: EXTEND: generation of interpolatory quadrature rules of the highest degree of precision with preassigned nodes for general weight functions,
Transactions on Mathematical Software,
Volume 15, Number 2, June 1989, pages 137-143.
7. Arthur Stroud, Don Secrest,
Prentice Hall, 1966,
LC: QA299.4G3S7.

### Examples and Tests:

The command

patterson_rule_compute 15 -1.0 +1.0 gp015

creates the following "region", "weight" and "point" files:

### List of Routines:

• MAIN is the main program for PATTERSON_RULE_COMPUTE.
• TEST01: extension of 3 point Gauss-Legendre rule.
• ASSIGN generates the polynomial whose roots are the preassigned nodes.
• BAIR seeks roots of a polynomial.
• CH_CAP capitalizes a single character.
• CHECK tests a computed quadrature rule.
• CH_EQI is a case insensitive comparison of two characters for equality.
• CH_TO_DIGIT returns the integer value of a base 10 digit.
• DGEFA factors a real general matrix.
• DGESL solves a real general linear system A * X = B.
• EPROD expands a product of two orthogonal polynomials.
• GENER calculates the polynomial defining the optimal new nodes.
• GET_UNIT returns a free FORTRAN unit number.
• IDAMAX indexes the array element of maximum absolute value.
• LFACT removes a linear factor from a polynomial expansion.
• NEWTON applies Newton's method for a single root of a polynomial.
• ORDER_CHECK checks the value of ORDER.
• QFACT divides a polynomial by a quadratic factor.
• R8MAT_WRITE writes an R8MAT file.
• RECURA is the recurrence used for Gauss, Lobatto and Radau rules.
• RESCALE rescales a Legendre quadrature rule from [-1,+1] to [A,B].
• ROOTS calculates roots of a quadratic factor.
• RSORT carries out a simple ripple sort.
• RULE_WRITE writes a quadrature rule to a file.
• SOLVE calculates roots of an orthogonal polynomial expansion.
• S_TO_I4 reads an I4 from a string.
• S_TO_R8 reads an R8 value from a string.
• TIMESTAMP prints the current YMDHMS date as a time stamp.
• TRANSF scales a polynomial expansion with respect to the moments.