Hybrid Gauss-Trapezoidal Quadrature Rules.
A new class of quadrature rules for the integration of both regular and singular functions is constructed and analyzed. For each rule the quadrature weights are positive and the class includes rules of arbitrarily high-order convergence. The quadratures result from alterations to the trapezoidal rule, in which a small number of nodes and weights at the ends of the integration interval are replaced. The new nodes and weights are determined so that the asymptotic expansion of the resulting rule, provided by a generalization of the Euler--Maclaurin summation formula, has a prescribed number of vanishing terms. The superior performance of the rules is demonstrated with numerical examples and application to several problems is discussed.
r
2.2. Euler{Maclaurin summation formula. For a function f 2 Cp(R), p
1, the Euler{Maclaurin summation formula (see, for example, [15, (23.1.30)]) can be
derived by repeated integration by parts. We rst consider the interval
apply (5) to obtain
Z c+h Z 1
where we have derive the Euler{Maclaurin formula for the interval [a; b], we let h =(ba)=n and
rearrange terms to obtain

f(a) f(b)
a
The expression on the left-hand side of (8) is the well-known trapezoidal rule. Evaluation
of the expression on the right-hand side of (8) is simplied by the fact that
2.3. Generalized Riemann zeta-function. The generalized Riemann zeta-
function is dened by the formula
(v
This function has a continuation that is analytic in the entire complex s-plane, with
the exception of where it has a simple pole. In what follows, we shall be
concerned primarily with real s and v, with s<1 and v>0. We will use the
HYBRID GAUSS-TRAPEZOIDAL QUADRATURE RULES 1555
following representation derived from Plana's summation formula (see, for example,
[16, section 1.10 (7)]):
v1s vs Z 1 sin(s arctan t) dt
Equation can be used to derive the asymptotic expansion of as v !1.We
treat the integral as a sum of Laplace integrals, each with an asymptotic expansion
given by Watson's lemma (see, for example, [17, p. 263]), and obtain
as v !1, with an arbitrary positive integer. Equation (10) isa slight generalization of [16, section 1.18 (9)]. There is a direct connection between
the Bernoulli polynomials and ,
and generalizations of the dierence and dierentiation formulae hold:
(v
@(s; v)
2.4. Orthogonal polynomials and Gaussian quadrature. Suppose that !
is a positive continuous function on the interval (a; b) and ! is integrable on [a; b].
We dene the inner product with respect to ! of real-valued functions f and g by the
integral
a
There exist polynomials p0;p1;:::, of degree 0; 1;:::, respectively, such that (pn;pm)=
are unique up to the choice of leading coecients.With leading coecients one, they can be obtained recursively by the formulae (see,
for example, [18, p. 143])
where p1(x)=0andn;n are dened by the formulae
The zeros xn1 ;::: ;xn of pn are distinct and lie in the interval (a; b). There exist
positive numbers !1n;::: ;!n such that
a
whenever f is a polynomial of degree less than 2n. These Christoel numbers are
given by the formula (see, for example, [19, p. 48])
Moreover, if !(x)=(b  x) (x) with integrable on [a; b], then, with the denition
exist positive numbers 1n;::: ;n+1 such that
a
whenever f is a polynomial of degree less than or equal to 2n. These modied Christoffel
numbers are given by the formula> !in
a pn(b)
where !in is given by (16).
The summation in (15) is the n-node Gaussian quadrature with respect to !,
while that in (17) is an (n + 1)-node Gauss{Radau quadrature with respect to .
3. Hybrid Gauss-trapezoidal quadrature rules. In this section we introduce
new quadrature rules for regular integrands, singular integrands with a power
or logarithmic singularity, and improper integrals and determine their rate of convergence
as the number of quadrature nodes increases.
For notational convenience we generally consider quadratures on canonical inter-
vals, primarily [0; 1]. It is understood that these are readily transformed to quadratures
on any nite interval [a; b] by the appropriate linear transformation of the nodes
and weights.
3.1. Regular integrands. For j;n positive integers and a 2
0g, we dene a linear operator Tnja on C([0; 1]), depending on nodes x1;::: ;xj and
weights w1;::: ;wj, by the formula
where h =(n +2a  1)1 is chosen so that ah +(n  1)h =1 ah.
Theorem 3.1. Suppose f 2 Cp([0; 1]). The asymptotic expansion of Tnja(f) as
!1is given by the formula
HYBRID GAUSS-TRAPEZOIDAL QUADRATURE RULES 1557
Proof. We apply the Euler{Maclaurin formula (8) on the interval [ah; 1  ah]to
obtain
We now combine (19) and (21), the equality
Taylor expansion of all quantities about the Bernoulli polynomial expansion
formula (6), and dierence formula (4) to obtain (20).
Corollary 3.2. Suppose the nodes x1;::: ;xj and weights w1;::: ;wj satisfy
the equations
r +1
Then Tnja is a quadrature rule with convergence of order 2j +1 for f 2 Cp([0; 1]) with
2j. Moreover,
as n !1; provided and the remaining nodes x1;::: ;xj1 and
weights w1;::: ;wj satisfy the equations
2:
r +1
Then Tnja is a quadrature rule with convergence of order 2j for f 2 Cp([0; 1]) with
as n !1; provided f(2j1)(0)  f(2j1)(1) =0.We shall see below that (22) has a solution with the nodes and weights all positive
if a is suciently large and that numerical solution of (22) is equivalent to computing
the roots of a particular polynomial. This statement holds for (24) as well.
3.2. Singular integrands. For j;k;n positive integers and a; b 2 R+, we dene
a linear operator Snjkab on C((0; 1]), depending on nodes v1;::: ;vj;x1;::: ;xk and
weights u1;::: ;uj;w1;::: ;wk, by the formula
chosen so that ah +(n  1)h =1 bh. The following
theorem, which follows from a generalization of the Euler{Maclaurin formula due to
Navot [8] or a further generalization due to Lyness [20], presents a somewhat dierent
proof than the earlier ones.
Theorem 3.4. Suppose g(x)=xf(x); where >1 and f 2 Cp([0; 1]). The
asymptotic expansion of Snjkab(g) as n !1is given by the formula
Proof.Forc 2 R+, we dene polynomials pc0;pc1;:::, in analogy with the Bernoulli
polynomials, by the formula

r
Dierentiating, we verify that
d
dx
and, combining the dierence formula (11) with (28), we obtain
c;n=0;
Additionally, we dene functions q0c;q1c;::: by the
and observe that
d
qnc
dx
HYBRID GAUSS-TRAPEZOIDAL QUADRATURE RULES 1559
With these denitions, the proof follows that of the Euler{Maclaurin formula:
Z 1bh nX2 Z 1
Taylor expansion of f(r)(ah) about the denitions (28) and (29) for pcn and qnc ,
and the binomial theorem combine to yield
Likewise, Taylor expansion of f(r)(1  bh) about the denitions for pcn and qnc ,
the asymptotic expansion (10) for , the Bernoulli polynomial expansion formula (6),
and the binomial theorem combine to yield
We now combine (26) and (30){(32), the equality
expansion of the latter two integrals about and the dierence formula (4) for
the Bernoulli polynomials and (11) for to obtain (27).
Corollary 3.5. Suppose the nodes u1;::: ;uj and weights v1;::: ;vj satisfy the
equations
and the nodes x1;::: ;xj and weights w1;::: ;wj satisfy the equations
r +1
Then Snjjab is a quadrature rule with convergence of order 2j +1+minf0;g for
where g(x)=xf(x); with f 2 Cp([0; 1]) and p 2j. Moreover,
as n !1:
Corollary 3.6. Suppose the nodes u1;::: ;uj and weights v1;::: ;vj satisfy the
equations
r +1
and the nodes x1;::: ;xk and weights w1;::: ;wk satisfy the equations
Xk r Br+1(b)
r +1
Then Snjkab is a quadrature rule with convergence of order minfj +1;+j +1; 2k+1g
for minfj;2kg.
In Corollaries 3.5 and 3.6, an even number of constraints on the nodes and weights
at both ends of the interval are considered. Clearly, there are analogous quadrature
rules arising from an odd number of constraints at one or both ends; these are similar,
and explicit presentation of them is omitted. We now consider a dierent type of
singularity.
Theorem 3.7. Suppose g(x)=f(x) log x; where f 2 Cp([0; 1]). The asymptotic
expansion of Snjkab(g) as n !1is given by the formula
HYBRID GAUSS-TRAPEZOIDAL QUADRATURE RULES 1561
where 0 denotes the derivative of with respect to its rst argument.
Proof. This asymptotic expansion is derived from that of Theorem 3.4 by dierentiating
(27) with respect to  and evaluating the result at  =0.
Corollary 3.8. Suppose the nodes u1;::: ;uj and weights v1;::: ;vj satisfy the
equations
ui vir log
r +1
and the nodes x1;::: ;xk and weights w1;::: ;wk satisfy the equations
Xk r Br+1(b)
r +1
Then Snjkab is a quadrature rule with error of order O(minfhj+1 log h; h2k+1g) for
where g(x)=(x) log x minfj;2kg.
3.3. Improper integrals. For j;n positive integers, we dene a linear operator
Rnj on C([n; 1)), depending on nodes x1;::: ;xj and weights w1;::: ;wj, by the
Theorem 3.9. Suppose g(x)=eixf(x); where  2 R; =0; and f 2 Cp([1; 1));and that there exist positive constants ; 0;::: ;p; such that
x+r
The asymptotic expansion of Rnj (g) as n !1is given by the formula
Proof. We integrate by parts repeatedly to obtain
and in (41) we compute the Taylor expansion of f about
lies between 0 and xk for k =1;::: ;j. Now combining (42), (44), and (45),
we obtain (43).
Example. The function f dened by the formula
ar
x+r
with jarj < 1, satises the assumptions of Theorem 3.9 for every positive integer
p. We remark that Theorem 3.9 can, in some instances, be generalized to
the corresponding asymptotic expansion depends on a more detailed knowledge of f.
For f given by (46), for example, the quadrature nodes and weights for
on .
Corollary 3.10. Suppose f 2 Cp([1; 1));fsatises (42) for x 2 R; and f is
analytic in the half-plane Re(x) >afor some a 2 R. Suppose further that v1;::: ;vj
are the roots of the Laguerre polynomial Lj of degree j, that coecients u1;::: ;uj
satisfy the equations
and that the operator Rnj is dened with nodes xk =(i=)vk and weights wk =(i=)uk
for k =1;::: ;j. Suppose nally that T^mja is dened to be the quadrature rule Tmja with
nodes and weights satisfying (22) but translated and scaled to the interval [1;n]. Then
for p 2j, the expression T^(jna1)n(g)+Rnj (g) is an approximation for the integral
with error of order O(n2j) as n !1.Proof. This result is just a combination of the quadrature rule of Corollary 3.2,
for the interval [1;n], with the asymptotic expansion of Theorem 3.9, for the interval
[n; 1), provided
But (48) follows from (47), the equations
(see, for example, [15, (6.1.1) and (22.2.13)]), and the fact that Gaussian quadratures,
which are exact for polynomials of degree less than twice the number of nodes, have
nodes that coincide with the roots of the corresponding orthogonal polynomials (see
section 2.4).
We have completed the denition of the new quadratures, along with the demonstration
of their asymptotic performance. We shall see that the existence of these
rules, which depends on the solvability of the nonlinear systems of equations that
dene the nodes and weights, is assured by the theory of Chebyshev systems. The
uniqueness of the rules is similarly assured. These issues of existence and uniqueness
are treated next.
HYBRID GAUSS-TRAPEZOIDAL QUADRATURE RULES 1563
4. Existence and uniqueness.
4.1. Chebyshev systems. Material of this subsection is taken, with minor
alterations, from Karlin and Studden [21]. Suppose I is an interval of R, possibly
innite. A collection of n real-valued continuous functions f1;::: ;fn dened on I is
a Chebyshev system if any linear combination
with ai not all zero, has at most n  1 zeros on I. This condition is equivalent to the
statement that for distinct x1;::: ;xn in I,
det B@ . CA =0:f1(xn) fn(xn)
The Chebyshev property is a characteristic of the space, rather than the basis: if
f1;::: ;fn is a Chebyshev system, then so is any other basis of spanff1;::: ;fng.If
u is a continuous, positive function on I, then scaling by u preserves a Chebyshev
system. Finally, if u is strictly increasing and continuous on interval J with range I,
then f1 u;::: ;fn u is a Chebyshev system on J if and only if f1;::: ;fn is on I.
(Here fi u denotes the composition u followed by fi.)
The best-known example of a Chebyshev system is the set of polynomials
on any interval I R. We shall be concerned also with the Chebyshev systems
log x; x; xlog x;::: ;x(n1)=2;x(n1)=2 log x
on I =(0;a], where  2 RnZ and a>0. These systems are special cases of the system
of Muntz functions (see, for example, [22, p. 133])
on I =(0; 1);Pwhere 1;::: ;j are distinct real numbers and n1;::: ;nj are positive
integers with To see this is a Chebyshev system, suppose f 2 span M and
use induction in n on (d=d log x)[f(x) xj ], in combination with Rolle's theorem.
Another Chebyshev system that will arise is the system
on I =[0; 1), where P1;::: ;j are distinct positive real numbers and n1;::: ;nj are
positive integers with n. This is indeed a Chebyshev system, for if f 2 spanL,
Qj
then the function f(x) i=1(x + i)ni is a polynomial in x of degree n  1.
Suppose f1;::: ;fn is a Chebyshev system on the interval I. The moment space
Mn with respect to f1;::: ;fn is the set
where the measure ranges over the set of nondecreasing right-continuous functions
of bounded variation on I. It can be shown that Mn is the convex cone associated
with points in the curve C, where
In other words, Mn can be represented as

The index I(c)ofapointc of Mn is the minimum number of points of C that can be
used in the representation of c, under the convention that a point (f1(x);::: ;fn(x))
is counted as a half point if x is from the boundary of I and receives a full count
otherwise. The index of a quadrature involving x1;::: ;xp is determined by counting
likewise.
Proofs of the next three theorems are somewhat elaborate and are omitted here;
they can be found in Karlin and Studden [21].
Theorem 4.1. (See [21, p. 42].) Suppose I =[a; b] is a closed interval. A point
is a boundary point of Mn if and only if I(c) <n=2: Moreover, ifis a measure corresponding to a boundary point c 2Mn, then there is a unique
quadrature
I
Theorem 4.2. (See [21, p. 47].) Suppose I =[a; b] is a closed interval. Any
point c in the interior of Mn satises I(c)=n=2. Moreover, if is a measure
corresponding to c; then there are exactly two quadratures
I
of index n=2; where wi > 0;i=1;::: ;p. In particular, if n =2m; then
Theorem 4.3. (See [21, p. 65].) Let Pn denote the nonnegative linear combinations
of functions f1;::: ;fn,

HYBRID GAUSS-TRAPEZOIDAL QUADRATURE RULES 1565
The point c =(c1;::: ;cn) is an element of Mn if and only if
implies ai ci 0:
Moreover, c is in the interior of Mn if and only if
imply ai ci > 0:
Theorem 4.4. (See [21, p. 106].) Suppose fi(x)=xi1 for i =1;::: ;n; and
I =[a; b].Ifn =2m, then c =(c1;::: ;cn) is an element of Mn if and only if the
two quadratic forms
Xm Xm
are nonnegative denite. If n =2m+1; then c 2Mn if and only if the two quadratic
forms
are nonnegative denite. Moreover, for either parity of n, c is in the interior of Mn
if and only if the corresponding quadratic forms are both positive denite.
Proof. A theorem of Lukacs (see, for example, [19, p. 4]) states that a polynomial
f of degree n  1 that is nonnegative on [a; b] can be represented in the form
where p and q are polynomials such that the degree of each term in (64) does not
exceed n  1. The combination of (64) and Theorem 4.3 proves the theorem.
4.2. Muntz system quadratures. The systems of (22), (24), (36), and (39)
that dene the quadrature rules of section 3 are special cases of the system of equations
dXn=2e
for distinct real numbers 1;::: ;j and positive integers n1;::: ;nj with
Here (k) denotes the kth derivative of with respect to its rst argument. The existence
and uniqueness of the solution of (65) follow from the existence and uniqueness
of quadratures for Chebyshev systems, once it is established that there is a measure
a with
Z a
in other words, that the moment space Mn of the Chebyshev system of Muntz function

on (0;a] contains the point
We will show that this condition is satised provided that a is suciently large. It
would be convenient to have tight bounds for a, in particular for systems (22), (24),
(36), and (39), but it appears that such bounds are dicult to obtain. Even for the
regular cases (22) and (24), where by Theorem 4.4 the existence of a is equivalent
to the positive deniteness of two matrices, precise bounds for arbitrary j appear
dicult. (Numerical examples below provide evidence that a=j may be chosen as
small as 5/6.)
Theorem 4.5. Suppose 1;::: ;j are diPstinct real numbers, each greater than
1, and n1;::: ;nj are positive integers with n. Then for suciently large
a, there exists a measure a such that the system of (66) is satised and c dened by
(68) is in the interior of the moment space Mn.
Proof. We construct a continuous weight function a0 satisfying (66) and show
that for suciently large a, a0 (x) is positive for x 2 [0;a].
We linearly combine the equations of (66) to obtain the equivalent system
where we have used the binomial theorem to expand logk(x=a) = (log xlog a)k.We
dene the weight a0 by the formula
and combine (69), (70), and the equalities
x log xdx= ;>1;k=0; 1;::: ;
to obtain the equations in 0;a;::: ;n1;a,
This n-dimensional linear system is nonsingular, since the set of functions
forms a Chebyshev system on [0; 1), as established at (51). Thus (71) possesses a
unique solution 0;a;::: ;n1;a.
HYBRID GAUSS-TRAPEZOIDAL QUADRATURE RULES 1567
We now determine
The asymptotic expansion of (r)(i;a)asa !1can be derived by dierentiating
(10); the rst several terms are given by
Combining (71) and (73), changing the order of summation, and twice applying the
product dierentiation rule
dr  Xr r
dr s
we obtain
which immediately reduces to
The combination of (70), (72), and (74) gives
lim a0 (ax)=1;x2 [0; 1];
which implies that for a suciently large, a0 (x) > 0 for x 2 [0;a]. The point c
dened by (68) is in the interior of Mn, since small perturbations of c will preserve
the positivity of a0 .
Theorem 4.2 ensures the existence of Gaussian quadratures for a Chebyshev system
f1;::: ;fn dened on an interval I, under the assumption that I is closed, whereas
the system M of (67) is Chebyshev on I =(0;a]. As a consequence, we require the
following result.
Theorem 4.6. Suppose the collection of functions f1;::: ;fn forms a Chebyshev
system on I =(a; b] and each is integrable on [a; b] with respect to a measure corresponding
to a point c in the interior of Mn. Then there exists a unique quadrature
a
of index n=2; where wi > 0 and xi 2 I; i =1;::: ;p. In particular, if n =2m; then
Proof. The Chebyshev property implies that there exists with a<<bsuch
that fi is nonzero on (a; ], i =1;::: ;n. We dene the function u on I by the formula
and observe that u is continuous and positive on I and integrable on [a; b] with respect
to . Now we dene functions g1;::: ;gn on [a; b] by the formula< fi(x)
;x2 (a; b];
x!a
The system g1;:R::x;gn is a Chebyshev system on [a; b] and is integrable with respect
to the measure u(t) d(t). Theorem 4.2 therefore is applicable and ensures the
a
existence of exactly two quadratures of index n=2 for the interval [a; b], one of which
includes the point Our assumption xi 2 (a; b] excludes this case and we are
left with the single quadrature presented in (76) and (77).
The next theorem, which is the principal analytical result of this section, follows
directly from Theorems 4.5 and 4.6. The existence and uniqueness of the quadratures
dened in section 3 follow from it. It also hints at the existence of somewhat more
general quadratures, for singularities of the form x logk x, but we do not evaluate
these here.
Theorem 4.7. Suppose 1;::: ;j are diPstinct real numbers, each greater than
and n1;::: ;nj are positive integers with n. For suciently large a, the
system of equations
dXn=2e
has a unique solution w1;::: ;wdn=2e;x1;::: ;xdn=2e satisfying wi > 0 for
dn=2e and 0 <x1 < <xdn=2e a; with xdn=2e = a if n is odd.
5. Computation of the nodes and weights. The nodes and weights of the
quadratures dened in section 3 are computed by numerically solving the nonlinear
systems (22), (24), (36), and (39). Conventional techniques for this problem either are
overly cumbersome or converge too slowly to be practical. Recently, Ma, Rokhlin, and
Wandzura [14] addressed this need by developing a practical numerical algorithm that
is eective in a fairly general setting. They construct a simplied Newton's method
and combine it with a continuation (homotopy) method. We present their method in
an abbreviated form in section 5.2; the reader is referred to [14] for more detail.
The systems for regular integrands, however, can be solved even more simply, as
we see next.
5.1. Regular integrands. The classical theory of Gaussian quadratures for
polynomials, summarized in section 2.4, can be exploited to solve (22) and (24). In
particular, suppose that p0;::: ;pj are the orthogonal polynomials on [0;a], given by
HYBRID GAUSS-TRAPEZOIDAL QUADRATURE RULES 1569
the recurrence (13){(14), under the assumption that
Z a
r Br+1(a)
x
Then the roots x1;::: ;xj of pj and the corresponding Christoel numbers w1;::: ;wj
satisfy (22). The polynomials p1;::: ;pj can be calculated in symbolic form; their
coecients are rational if a is rational. The roots of pj can be computed by Newton
iteration and the Christoel numbers can be obtained using (16). Similar treatment
can be applied to the system (24) containing an odd number of equations, using the
interval [0;a 1], under the assumption
Z
r Br+1(a)
x 2:
The Gauss{Radau quadrature is computed using the formula (18) for the modied
Christoel numbers. For these tasks it is convenient to use a software system that
can manipulate polynomials with full-precision rational coecients. The author implemented
code for these computations in Pari/GP [23].
It should be noted that the proposed procedure is suitable for relatively small
values of j (less than, say, 20). It is neither very ecient nor very stable, but it
was quite adequate for our purposes. (Unlike the situation for standard Gaussian
quadratures, where the number of nodes depends on the size of the problem, here
depends only on the desired order of convergence.) If it is required to compute
the nodes and weights of (22) or (24) for large j, the reader may consider numerical
schemes for Gaussian quadrature proposed by other authors, for example, that of
Gautschi [24] or Golub and Welsch [25].
5.2. Singular integrands. The systems (36) and (39) for singular integrands
cannot be solved using methods for standard Gaussian quadrature, since the nodes
to be computed do not coincide with the roots of any closely related orthogonal
polynomials. We employ instead the algorithm for such systems developed by Ma,
Rokhlin, and Wandzura [14], which we now describe.
A collection of 2n real-valued continuously dierentiable functions f1;::: ;f2n
dened on an interval I =[a; b]isaHermite system if
for any choice of distinct x1;::: ;xn on I. A Hermite system that is also a Chebyshev
system is an extended Hermite system. The following theorem is a direct consequence
of the denition; the proofs of the subsequent two theorems are contained in [14].
Theorem 5.1. Suppose that the functions f1;::: ;f2n constitute a Hermite system
on the interval [a; b] and x1;::: ;xn are n distinct points on [a; b]. Then there
exist unique coecients ij;ij;i=1;::: ;n; j =1;::: ;2n; such that
for i =1;k=1;::: ;n; where the functions i;i are dened by the formulae
for i =1;::: ;n.
Theorem 5.2. (See [14, p. 979].) Suppose that the functions f1;::: ;f2n constitute
a Hermite system on [a; b]. Suppose further that S [a; b]n is the set of points
with distinct coordinates x1;::: ;xn. Suppose nally that the mapping F
dened by the formula
a a
with the functions 1;::: ;n dened by (80){(82). Then x1;::: ;xn are the Gaussian
nodes for the system of functions f1;::: ;f2n with respect to the weight ! if and only
if F(x1;::: ;xn)=0.
Theorem 5.3. (See [14, p. 983].) Suppose that the functions f1;::: ;f2n are
an extended Hermite system on [a; b];S [a; b]n is the set of points with distinct
coordinates, and the mapping G dened by the formula
a 1(x) !(x) dx a n(x) !(x) dx
a 1(x) !(x) dx a n(x) !(x) dx
with the functions 1;::: ;n and 1;::: ;n dened by (80){(83). Suppose further
that fi 2 C3((a; b)) for i =1;::: ;2n and the function F is dened by (84). Suppose
nally that x is the unique zero of F; that x0 is an arbitrary point of S; and that the
sequence x1; x2;::: is dened by the formula
Then there exists >0 and >0 such that the sequence x1; x2;::: generated by (86)
converges to x and
for any initial point x0 such that kx0  xk <.
The key feature of this theorem is the quadratic convergence indicated by (87).
The solution x is obtained by an iterative procedure; each step consists of computing
the coecients that determine 1;::: ;n and 1;::: ;n by inverting the matrix
of (79) then computing the integrals that dene G by taking linear combinations,
using these coecients, of integrals of f1;::: ;f2n. Theorem 5.3 ensures that with
appropriate choice of starting value x0, convergence is rapid and certain.
For quadrature nodes x =(x1;::: ;xn), the quadrature weights are given by the
integrals of i, namely,
a a
HYBRID GAUSS-TRAPEZOIDAL QUADRATURE RULES 1571
We note that Theorems 5.1{5.3 concern Gaussian quadratures with n nodes and
weights to integrate 2n functions on the interval [a; b] exactly. For Gauss{Radau
quadratures, in which node xn = b is xed and 2n1 functions are integrated exactly,
only a slight change is required. In particular, functions 1;::: ;n1 (without n)
and 1;::: ;n are dened as before by (80){(83), except that the summations in (82)
and (83) exclude f2n. Their coecients ij, i =1;::: ;n 1, j =1;::: ;2n  1, and
are obtained by inverting the matrix which results
from removing the last row and column from the matrix of (79). The revised mapping
components dened as the rst n  1 components in (85).
Finally, as before, the quadrature weights are given by (88).
In order to obtain a suciently good starting estimate for the solution of F(x)=0,
a continuation procedure can be used, as outlined in the following theorem.
Theorem 5.4. (See, for example, [14, p. 975].) Suppose that F :[0; 1]
Rn is a function with a unique solution xt to the equation F(t; x)=0
for all t 2 [0; 1]; suppose that xt is a continuous function of t; and suppose that
Finally, suppose x0 is given and for some >0 there is a procedure
P to compute xt for t 2 [0; 1]; given an estimate x~t with jx~t  xtj <. Then there
exists a positive integer m such that the following procedure can be used to compute
the solution of F(x)=0:
For use P to compute xi=m; given the estimate x(i1)=m.
The required solution of F(x)=0 is x1.
More typically, of course, and any bound on jxt+  xtj= depend on t and in a
practical implementation the step size is chosen adaptively.
To compute the solutions of (36) and (39), it is eective to use a continuation
procedure with respect to both j and a. Solutions for the rst few values of j are
readily obtained without requiring good initial estimates. Given a solution of (36) for
the interval [0;a] with nodes u1;::: ;uj and weights v1;::: ;vj, we choose an initial
estimate u~1;::: ;u~j+1,v~1;::: ;v~j+1 for and the interval [0;a+ 1] dened by the
ui;i=1;::: ;j; vi;i=1;::: ;j;
This choice exactly satises the equations
r +1
as follows immediately from the dierence formula (4) for Bn and (11) for , but the
corresponding equations for are not satised. Those equations are approximately
satised, however, and we can start with the actual values of the sums for
as the required values. These are then varied continuously, obtaining the corresponding
solutions, until they coincide with the intended values (  j;a +1)
and Bj+1(a +1)=(j + 1). This procedure can be used without alteration for (39).
Once the solution for and the interval [0;a + 1] is obtained, a can be
continuously varied, in a continuation procedure, to obtain solutions for dierent
intervals.
1572 BRADLEY K. ALPERT

Table

The minimum value of a, as a function of j, such that the point B1(a)=1;::: ;Bj(a)=j is in
the moment space M2j of the polynomials 1;x;::: ;x2j1 on the interval [0;a]. The moment space
is dened at (52).
j min a j min a j min a j min a
6 4.77448
9 7.21081
14 11.29815
6. Numerical examples. The procedures described in section 5 were implemented
in Pari/GP [23] for both the regular cases and the singular cases. The matrix
in (79), which must be inverted, is very poorly conditioned for many choices of n,
x1;::: ;xn, and f1;::: ;f2n. This diculty was met by using the extended precision
capability of Pari/GP.
6.1. Nodes and weights. The nodes and weights of (22), (24), (36), and (39)
that determine the quadratures of section 3 were computed for a range of values of the
parameter j. For each choice of j, a was chosen, by experiment, to be the smallest integer
leading to positive nodes and weights (see Theorem 4.7). For the regular case (22),
the characterization expressed in Theorem 4.4 was used to determine the minimum
value of a 2 R, for j =1;::: ;16, that satises
In particular, we obtained the minimum value of a such that the quadratic forms in
are nonnegative denite. This determination was made by calculating the determinant
of each corresponding matrix symbolically and solving for the largest root
of the resulting polynomial in a. These values are given in Table 1. This evidence
suggests that limj!1 j1 min a exists and is roughly 5=6, meaning that the number
of trapezoidal nodes displaced is less than the number of Gaussian nodes replacing
them, for quadrature rules of all orders. This relationship also appears to hold, to an
even greater extent, for the singular cases.
The values of selected nodes and weights, for the regular cases and for singularities
x1=2 and log x, are presented in an appendix. Of particular simplicity are the rst
two rules for regular integrands,


48 48

48 48
These rules are of third- and fourth-order
convergence, respectively. The rst is noteworthy for having the same weights as, but
higher order than, the trapezoidal rule; the second has asymptotic error 1=4 that of
Simpson's rule with the same number of nodes.
The lowest-order rule presented for logarithmic singularities,

HYBRID GAUSS-TRAPEZOIDAL QUADRATURE RULES 1573

Table
Relative errors in the computation of the integral in (94), for the regular case s(x)=0. Quadrature
rules with convergence of order 2; 4; 8; 16, and 32 were used with various numbers
of nodes. Here is the oversampling factor.
90 1.41
100 1.57
200 3.14
260 4.08

Table
Relative errors for the singular case s(x)=x1=2, for various numbers of nodes
and orders of convergence.
90 1.41
100 1.57
200 3.14
260 4.08
approximates g(x) dx with error of order O(h2 log h) for g(x)=(x) log x are regular functions on [0; 1]. The corresponding rule for the
singularity x1=2 is not quite as simple, for there the quantity 2 in (93) is replaced
with 4(1=2)2.
6.2. Quadrature performance. To demonstrate the performance of the quadrature
rules, they were used in a Fortran routine (with real*8 arithmetic) to numerically
compute the integrals
dx;for the functions s(x)=0,s(x)=x1=2, and log x. These integrals were also
obtained analytically and the relative error of the quadratures was computed. The
numerical integrations were computed for various orders of quadrature and various
numbers of nodes. Minimum sampling was taken to be two points per period of the
cosine (i.e., 200= 63:7 quadrature nodes). The accuracies were then compared for
various degrees of oversampling. The quadrature errors are listed in Tables 2{4 and
plotted, as a function of oversampling factor, in Figure 1. We note that the graphs are

Table
Relative errors for the singular case log x. Here the error is of order O(hl log h), where
l is shown.
90 1.41
100 1.57
200 3.14
260 4.08
-2-4-6Relative
Oversampling (f) Oversampling (f) Oversampling (f)
Fig. 1. The relative errors of the quadratures, shown in Tables 2{4, are plotted using logarithmic
scaling on both axes.
nearly straight lines (until the limit of machine precision is reached), as predicted from
the theoretical convergence rates. We remark also that excellent accuracy is attained
for even quite modest oversampling when quadratures with high-order convergence
are employed. For problems where the number of quadrature nodes is the major cost
therefore, one may benet by using the high-order quadratures even for modest
accuracy requirements.
HYBRID GAUSS-TRAPEZOIDAL QUADRATURE RULES 1575

Table
Relative errors in the computation of the integral in (95), for =1. Quadrature rules dened
in Corollary 3:10 with j =1; 2; 4; 8, and 16 were used with various numbers m of nodes.
We test the quadratures for improper integrals by numerically computing for
1 the integral
where H is the Heaviside step function. The integrand is oscillatory and decays
like x1 as x !1. The quadratures dened in Corollary 3.10 are employed, for
which the integral is split into a regular integral on a nite interval, chosen here
to be [5 m=4; 5 m=4], where m is the total number of quadrature nodes, and
two improper integrals in the imaginary direction, using Laguerre quadratures. The
quadrature errors are shown in Table 5.
7. Applications and summary. The chief motivation for the hybrid Gauss-
trapezoidal quadrature rules is the accurate computation of integral operators. We
dene an integral operator A by the formula
Z
where  is a regular, simple closed curve in the complex plane, the function f is
regular on , and the kernel regular function of its arguments,
except where they coincide; we assume
with and ^ regular on   and s regular on (0; 1), with an integrable singularity
at 0. A large variety of problems of classical physics can be formulated as integral
equations that involve such operators. When the operator occurs in an integral equa-
tion
some choice of discretization must be used to reduce the problem to a nite-dimensional
one for numerical solution. In the Nystrom method the integrals are replaced by
quadratures to yield the nite system of equations
Xm
1576 BRADLEY K. ALPERT
This linear system can be solved for f(x1);::: ;f(xm) by a variety of techniques. The
particular choice of xi and wij for ;mdetermines the order of convergence
(and therefore eciency) of the method.
For a curve parametrization :[0; , such as scaled arc length, the operator
A becomes
is convenient to use a uniform discretization 1=m; 2=m;::: ;1int and
How then is wij determined? We assume for the moment
that f is available at locations other than x1;::: ;xm. Continuing periodically with
period 1, and using the Gauss-trapezoidal quadratures, we obtain
Z 1+i=m
for dened by the formula
are determined for the singularity s of
K. Provided that the periodic continuation of is suciently regular, the quadrature
will converge to the integral with order greater than j as m !1, for i =1;::: ;m.
We relax the restriction that f be available outside x1;::: ;xm by using local Lagrange
interpolation of order equispaced nodes,
(j  1)=2c and
r  s
s=0;s=rNow wij for determined by combining (97){(101). The computation
of all m2 coecients requires m(m +2j  2a evaluations of the kernel K and
therefore order O(m2) operations. This cost can often be substantially reduced using
techniques that exploit kernel smoothness (see, for example, [6], [5]).
A slightly dierent application of the quadratures is the computation of Fourier
transforms of functions that fail to satisfy the assumptions usually made when using
the discrete Fourier transform. In particular, if a function decays slowly for large
argument or is compactly supported and singular at the ends of the support interval,
these quadratures can be used to compute its Fourier transform. One example of such
a function is that in (95). Since most of the nodes in these quadratures are equispaced,
with function values given equal weight, the fast Fourier transform can be used to
do the bulk of the computations; the overall complexity is O(n log n), where n is the
number of Fourier coecients to be computed. Other applications may include the
HYBRID GAUSS-TRAPEZOIDAL QUADRATURE RULES 1577
representation of functions for solving ordinary or partial dierential equations, when
high-order methods are required. In addition, an extension of these quadratures to
integrals on surfaces is under study.
In summary, the characteristics of the hybrid Gauss-trapezoidal quadrature rules
include
arbitrary order convergence for regular functions or functions with known
singularities of power or logarithmic type,
positive weights,
most nodes equispaced and most weights constant, and
invariant nodes and weights (aside from scaling) as the problem size increases.
The primary disadvantage of the quadrature rules, shared with other Gaussian quadratures
but exacerbated here by poor conditioning, is that the computation of the nodes
and weights is not trivial. Nevertheless, tabulation of nodes and weights for a given
order of convergence allows this issue to be avoided in the construction of high-order,
general-purpose quadrature routines.

Appendix

. Tables of quadrature nodes and weights.
Quadrature nodes and weights may also be obtained electronically from the author


Table

ThePnodes and weights for the quadrature rule
several choices of j and corresponding
minimum integer a.Forf a regular function, converges to f(x) dx as n !1withconvergence of order O.
Oa xi wi
HYBRID GAUSS-TRAPEZOIDAL QUADRATURE RULES 1579

Table
(Continued)
Oa xi wi
3.14968 50162 29433D01
1.39668 57813 42510D+00
2.17519 59032 06602D+00
3.06232 05758 80355D+00

Table
PThe nodes v1;:P:: ;vj and weights Pu1;::: ;uj for the quadrature rule Sn (g)=
g(x)=x1=2(x)+^(x), with and ^ regular functions. The nodes x1;::: ;xk and weights
w1;::: ;wk are found in Table 6.
Oa vi ui
1.00000 00000 00000D00 1.08019 20374 73384D+00
2.5 2 6.02387 37964 08450D02 2.85843 99904 20468D01
2.24632 55125 21893D01 4.87348 40566 46474D01
1.00000 00000 00000D+00 9.73575 20666 00344D01
2.69428 63467 92474D01 5.07743 45780 43636D01
2.57822 04347 38662D01 4.95598 17403 06228D01
2.00000 00000 00000D+00 1.00241 65465 50407D+00
8.28301 97052 96352D02 1.75524 44045 44475D01
4.13609 49257 26231D01 5.03935 05038 58001D01
1.08874 43736 88402D+00 8.26624 13396 80867D01
2.00648 21018 52379D+00 9.77306 58489 81277D01
3.00000 00000 00000D+00 9.99791 98099 47032D01
3.22395 27000 27058D02
1.79093 53836 49920D01
5.43766 38052 44631D01
1.17611 66283 96759D+00
2.03184 82107 16014D+00
4.00000 00000 00000D+00
6.19984 48842 97793D03
7.10628 67917 20044D02
2.40893 01044 10471D01
7.59244 65404 41226D01
9.32244 63996 14420D01
9.92817 14381 60095D01
9.99944 91256 89846D01
10.
HYBRID GAUSS-TRAPEZOIDAL QUADRATURE RULES 1581

Table
(Continued)
Oa vi ui
1.54042 43511 15548D02
8.83424 84071 96555D02
2.82446 20545 09770D01
6.57486 98923 05580D01
1.24654 10609 77993D+00
2.03921 84951 30811D+00
2.97933 34870 49800D+00
4.99724 08043 11428D+00
5.99986 87939 51190D+00
7.00000 00000 00000D+00
2.92101 89269 12141D03
3.43113 06112 56885D02
1.22466 94956 38615D01
2.76110 82420 22520D01
6.96655 56772 71379D01
8.79007 79419 72658D01
9.86862 24492 94327D01
1.00000 23977 96838D+00
14.0 9 3.41982 14602 49725D04
9.29659 34301 87960D03
4.21848 66056 53738D01
2.16099 75052 38153D+00
4.99973 27079 05968D+00
5.99987 51919 71098D+00
8.00000 00000 00000D+00
1.75095 72432 02047D03
2.08072 65842 87380D02
7.58683 06164 33430D02
3.20662 43620 72232D01
8.24495 90253 66557D01
9.84576 84431 63154D01
9.99285 27691 54770D01
1.00000 00814 05180D+00
5.89843 27437 09196D03
1.14558 64950 70213D01
2.79034 42188 56415D01
5.60011 37986 53321D01
9.81409 12428 83119D01
2.27017 91140 36658D+00
3.10823 46017 15371D+00
4.03293 08939 96553D+00

Table
PThe nodes v1;:P:: ;vj and weights Pu1;::: ;uj for the quadrature rule Sn (g)=
functions. The error is of order O(hl log h). The
nodes x1;::: ;xk and weights w1;::: ;wk are found in Table 6.
la vi ui
HYBRID GAUSS-TRAPEZOIDAL QUADRATURE RULES 1583

Table
(Continued)
la vi ui
2.44111 00950 09738D02
1.15385 12974 29517D01
7.32974 05318 07683D01
2.11435 87523 25948D+00
3.02608 45496 55318D+00
5.00014 11700 55870D+00
6.00000 10024 41859D+00
6.36419 07807 20557D03
4.72396 41432 87529D02
1584 BRADLEY K. ALPERT

Acknowledgments

. The author thanks Gregory Beylkin, Vladimir Rokhlin,
and Ronald Wittmann for helpful discussions and encouragement.


--R


On the Numerical Solution of One-Dimensional Integral and Dierential Equa- tions

A Nystro

Generalized Gaussian quadrature rules for systems of arbitrary functions
Handbook of Mathematical Functions (10th print- ing)
Higher Transcendental Functions
Advanced Mathematical Methods for Scientists and Engi- neers
Introduction to Numerical Analysis


Tchebyche Systems with Applications in Analysis and Statistics
Polynomials and Polynomial Inequalities
User's Guide to PARI-GP
On the construction of Gaussian quadrature rules from modied moments
Calculation of Gaussian quadrature rules
--TR

--CTR
Roberto Camassa , Jingfang Huang , Long Lee, Integral and integrable algorithms for a nonlinear shallow-water wave equation, Journal of Computational Physics, v.216 August 2006
Thomas Hagstrom , George Hagstrom, Grid stabilization of high-order one-sided differencing I: First-order hyperbolic systems, Journal of Computational Physics, v.223 n.1, p.316-340, April, 2007
Donna Calhoun, A Cartesian grid method for solving the two-dimensional streamfunction-vorticity equations in irregular regions, Journal of Computational Physics, v.176 n.2, p.231-275, March 1, 2002
Sheon-Young Kang , Israel Koltracht , George Rawitscher, Nystrm-Clenshaw-Curtis quadrature for integral equations with discontinuous kernels, Mathematics of Computation, v.72 n.242, p.729-756, 1 April
Bradley Alpert , Leslie Greengard , Thomas Hagstrom, Nonreflecting boundary conditions for the time-dependent wave equation, Journal of Computational Physics, v.180 n.1, p.270-296, July 20, 2002
high-order convergence;numerical integration;positive weights;singularity;gaussian quadrature;euler-maclaurin formula