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
      

        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

Index	Order	Free+Fixed	Expected Precision	Actual Precision
0	1	1 + 0	2*1+0-1=1	1
1	3	3 + 0	2*3+0-1=5	5
2	7	4 + 3	2*4+3-1=10	10 + 1 = 11
3	15	8 + 7	2*8+7-1=22	22 + 1 = 23
4	31	16 + 15	2*16+15-1=46	46 + 1 = 47
5	63	32 + 31	2*32+31-1=94	94 + 1 = 95
6	127	64 + 63	2*64+63-1=190	190 + 1 = 191
7	255	128 + 127	2*128+127-1=382	382 + 1 = 383
8	511	256 + 255	2*256+255-1=766	766 + 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

• 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.

Licensing:

The computer code and data files described and made available on this web page are distributed under the GNU LGPL license.

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,
   The Optimal Addition of Points to Quadrature Formulae,
   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,
   Gaussian Quadrature Formulas,
   Prentice Hall, 1966,
   LC: QA299.4G3S7.

Source Code:

• patterson_rule_compute.f90, the source code.

Examples and Tests:

The command

            patterson_rule_compute 15 -1.0 +1.0 gp015
          

• gp015_r.txt, the region file.
• gp015_w.txt, the weight file.
• gp015_x.txt, the point file.

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.
• EXTEND extends a quadrature rule by adding new nodes.
• 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.
• WEIGHT calculates quadrature weights.

You can go up one level to the FORTRAN90 source codes.