Bulletin T.CXXXVII de l’Académie serbe des sciences et des arts − 2008 Classe des Sciences mathématiques et naturelles Sciences mathématiques, N\underline{o} 33

QUADRATURE PROCESSES – DEVELOPMENT AND NEW DIRECTIONS

G. V. MILOVANOVIĆ

(Presented at the 7th Meeting, held on October 26, 2007)

A b s t r a c t. We present a survey on quadrature processes, beginning with Newton’s idea of approximate integration and Gauss’ discovery of his famous quadrature method, as well as significant contributions of Jacobi and Christoffel. Beside the stable construction of Gauss-Christoffel quadratures for classical and non-classical weights we give some recent applications in the summation of slowly convergent series and moment-preserving spline approximation. Also, we consider quadratures of the maximal degree of precision with multiple nodes, as well as a more general concept of orthogonality with respect to a given linear moment functional and corresponding quadratures of Gaussian type. A short account of non-standard quadratures of Gaussian type is also included. Finally, we mention the Gaussian integration which is exact on the space of Müntz polynomials.

AMS Mathematics Subject Classification (2000): 41A55, 33C45, 33C47, 42C05, 65D30, 65D32

Key Words: Quadrature process; Newton-Cotes formula, Gauss-Christo-ffel quadrature formula; Orthogonal polynomials; Moments; Moment functional; Three-term recurrence relation; Weight; Node; Multiple nodes; Summation of series; Moment-preserving spline approximation; Non-standard quadratures; Müntz polynomials

1. Introduction

In this paper we give an account of the classical Newton-Cotes and Gauss-Christoffel quadratures, as well as several their generalizations and extensions. We consider only one dimensional cases.

The paper is organized as follows. In Sections 2 and 3, beside the basic ideas of Isaac Newton (1643–1727) and Carl Friedrich Gauss (1777–1855), including significant contributions of Carl Gustav Jacob Jacobi (1804–1851) and Elwin Bruno Christoffel (1829–1900), we give several recent results on the Gauss-Christoffel quadrature rules for classical and non-classical weight functions, where the theory of orthogonal polynomials is a basic tool. We also mention some recent applications of such quadratures in the summation of slowly convergent series and the moment-preserving spline approximation.

In Section 4 we consider the quadratures with multiple nodes introduced fifty years ago by Chakalov, Popoviciu and Turán, but intensively studied using the s- and σ-orthogonality up to date.

Section 5 is devoted to the concept of orthogonality with respect to a given linear moment functional and to the corresponding quadrature rules of Gaussian type. As illustrations we mention only orthogonality on the semicircle introduced in our joint papers with Gautschi and Landau and recent results on orthogonal polynomials for oscillatory weights obtained with Cvetković. In the first case, the corresponding quadratures of Gaussian type can be applied for calculating the Cauchy principal value integrals and for numerical differentiation of analytic functions. In the second one, they can be applied for numerical calculation of integrals involving highly oscillatory integrands.

In Section 6 we give some considerations on the so-called non-standard quadratures of Gaussian type, which have been introduced and investigated recently. Such quadratures are based on operator values for a general family of linear operators, acting of the space of algebraic polynomials, such that the degrees of polynomials are preserved.

Finally, in Section 7 we consider the Gaussian integration which is exact on the space of Müntz polynomials. The extensions to trigonometric systems and the other non-polynomial systems are not included.

2. Theory of Newton and Gauss

The major advance in integration came in the 17th century with the independent discovery of the fundamental theorem of calculus by Newton and Leibniz. This theorem gives a connection between integration and differentiation. Formalizing integrals (B. F. Riemann) and formulating different definitions of integrals (H. Lebesgue, T. J. Stieltjes, S. Bochner, …), founded in measure theory, integration becomes an important modern theory.

Numerical integration also begins by Newton’s idea from 1676 for finding the weight coefficients A_{1}, A_{2}, …, A_{n} in the so-called n-point quadrature formula

I(f) = ∫_{a}^{b} f(t) dt ≈ Q_{n}(f) = A_{1}f(τ_{1}) + A_{2}f(τ_{2}) + ⋯ + A_{n}f(τ_{n}), (1)

for given (usually equidistantly) n points (nodes) τ_{1}, τ_{2}, …, τ_{n}, such that (1) is exact for all algebraic polynomials of degree at most n − 1, i.e., for each f ∈ 𝒫_{n−1}. In modern terminology, given n distinct points τ_{k} and corresponding values f(τ_{k}), Newton constructs the unique polynomial P ∈ 𝒫_{n−1}, which at the points τ_{k} assumes the same values as f, i.e., P (τ_{k}) = f(τ_{k}), k = 1, 2, …, n, expressing this interpolation polynomial in terms of divided differences. Evidently, for the interpolation error r_{n}(f; t) := f(t) − P (t) we have r_{n}(f; t) = 0 for all f ∈ 𝒫_{n−1}. Subsequently, integrating this polynomial over [a, b], Newton obtains (1). Here, we give it in a more convenient (Lagrange) form. Namely, the interpolation polynomial P can be expressed in terms of the so-called fundamental polynomials

ℓ_{k}(t) = \frac{ω_{n}(t)}{ω′_{n}(τ_{k})(t − τ_{k})}, k = 1, …, n,  (2)

where ω_{n}(t) = ∏_{k=1}^{n}(t − τ_{k}), and therefore A_{k} = I(ℓ_{k}), k = 1, 2, …, n, and R_{n}(f) = I(f) − Q_{n}(f) = I(r_{n}(f; ·)). Obviously, the remainder term R_{n}(f) = 0 for each f ∈ 𝒫_{n−1}. The quadrature formula obtained in this way is known as interpolatory and it has degree of exactness at least n − 1. We write d = d(Q_{n}) ≥ n − 1.

Starting from the work of Newton and Cotes and combining it with his earlier work on the hypergeometric series, Gauss in 1814 (see [23]) develops his famous method of numerical integration which dramatically improves the earlier method of Newton and Cotes. Namely, Gauss’ question was what is the maximum degree of exactness that can be achieved if the nodes τ_{k} are free. He started with the conjecture that max_{τ_{k},A_{k}} d(Q_{n}) = 2n − 1, since there are 2n unknowns in (1): τ_{k}, A_{k}, k = 1, …, n, and 2n conditions R_{n}(t^{ν}) = 0, ν = 0, 1, …, 2n − 1. To simplify consideration, we put (a, b) = (−1, 1).

Since

R_{n} (\frac{1}{z − ·}) = ∑_{ν=0}^{+∞} \frac{R_{n}(t^{ν})}{z^{ν+1}}, |z| > 1 ≥ |t|,

the problem can be reduced to

R_{n} (\frac{1}{z − ·}) = ∑_{ν=2n}^{+∞} \frac{R_{n}(t^{ν})}{z^{ν+1}} = O (\frac{1}{z^{2n+1}} ) , z → ∞.

On other hand

I (\frac{1}{z − ·}) = ∫_{−1}^{1} \frac{dt}{z − t} = log \frac{1 + z^{−1}}{1 − z^{−1}}

can be expressed as a continued function

I (\frac{1}{z − ·}) = \frac{2}{z − \frac{1/3}{z − \frac{2^{2}/(3 · 5)}{z − \frac{3^{2}/(5 · 7)}{z − \frac{4^{2}/(7 · 9)}{z − _{...}}}}}} ,

which was well known to Gauss (see [23]). Also, Gauss knew that the n-th rational convergent of the previous continued function, R_{n−1,n}(z) = p_{n−1}(z)/q_{n}(z) (dg p_{n−1} = n − 1, dg q_{n} = n), satisfies

I (\frac{1}{z − ·}) = R_{n−1,n}(z) + O (\frac{1}{z^{2n+1}} ) , z → ∞.

Expanding R_{n−1,n}(z) in partial fractions and taking poles (zeros of q_{n}(z)) as the nodes of the quadrature formula (1), Gauss obtains

R_{n−1,n}(z) = ∑_{k=1}^{n} \frac{A_{k}}{z − τ_{k}} , with A_{k} = Res_{z=τ_{k}} R_{n−1,n}(z), k = 1, …, n.

By such a quadrature formula Q_{n}, defined by Q_{n}(1/(z − ·)) := R_{n−1,n}(z), it follows that

R_{n} (\frac{1}{z − ·}) = I (\frac{1}{z − ·}) − Q_{n} (\frac{1}{z − ·}) = I (\frac{1}{z − ·}) − R_{n−1,n}(z),

i.e.,

R_{n} (\frac{1}{z − ·}) = O (\frac{1}{z^{2n+1}} ) , z → ∞.

Thus, d(Q_{n}) = 2n − 1. Gauss expressed the denominator and numerator polynomials (q_{n} and p_{n−1}) in terms of his hypergeometric series. Otherwise, these polynomial are known now as Legendre polynomials of the first and second kind.

No doubt that Gauss’ discovery is the most significant event of the 19th century in the field of numerical integration and perhaps in all of numerical analysis. An elegant alternative derivation of this famous method was provided by Jacobi, and a significant generalization to arbitrary measures was given by Christoffel. The error term and convergence were proved by Markov and Stieltjes, respectively. Today these formulae with maximal degree of precision are known as the Gauss-Christoffel quadrature formulae. A nice survey of Gauss-Christoffel quadrature formulae was written by Gautschi [24].

In modern terminology, the formulation of this classical theory can be given in the following form: Let dμ be a finite positive Borel measure on the real line ℝ such that its support is an infinite set, and all its moments μ_{k} = ∫_{ℝ} t^{k} dμ(t), k = 0, 1, …, exist and are finite. Then, for each n ∈ ℕ, there exists the n-point Gauss-Christoffel quadrature formula

∫_{ℝ} f(t) dμ(t) = ∑_{k=1}^{n} A_{k}f(τ_{k}) + R_{n}(f),  (3)

which is exact for all algebraic polynomials of degree ≤ 2n−1, i.e., R_{n}(f) = 0 for each f ∈ P_{2n−1}.

The Gauss-Christoffel quadrature formula can be characterized as an interpolatory formula for which its node polynomial ω_{n}(t) = ∏_{k=1}^{n}(t − τ_{k}) is orthogonal to 𝒫_{n−1} with respect to the inner product defined by

(f, g) = ∫_{ℝ} f(t)g(t) dμ(t) (f, g ∈ L^{2}(dμ)).  (4)

Thus, the nodes in (3) are zeros of the monic orthogonal polynomial ω_{n}(t) = π_{n}(dμ; t). The corresponding weights A_{k} (Christoffel numbers) can be expressed by the so-called Christoffel function λ_{n}(dμ; t) (cf. [53, Chapters 2 & 5]) in the form A_{k} = λ_{n}(dμ; τ_{k})(> 0), k = 1, …, n. In the special case dμ(t) = dt on [−1, 1] considered by Gauss, the nodes are zeros of the Legendre polynomial P_{n}.

If μ is an absolutely continuous function, then we say that μ′(t) = w(t) is a weight function and the measure dμ can be expressed as dμ(t) = w(t) dt, where the weight function w is non-negative and measurable in Lebesgue’s sense for which all moments exist and μ_{0} > 0. If the support of the measure (weight) is supp (w) = [a, b], where −∞ < a < b < +∞, we have Gaussian quadratures (and the corresponding orthogonal polynomials) on the finite interval [a, b].

In the general case, the function μ can be written in the form μ = μ_{ac} + μ_{s} + μ_{j}, where μ_{ac} is absolutely continuous, μ_{s} is singular, and μ_{j} is a jump function.

One of the most important uses of orthogonal polynomials is in the construction of quadrature formulas of the maximal, or nearly maximal, algebraic degree of exactness for integrals involving a positive measure dμ. The following theorem is due to Jacobi [44] (cf. [53, p. 297]).

Theorem 1. Given a positive integer m (≤ n), the quadrature formula (3) has degree of exactness d = n − 1 + m if and only if the following conditions are satisfied:

1^{∘} Formula (3) is interpolatory; 2^{∘} The node polynomial ω_{n} satisfies

(∀p ∈ 𝒫_{m−1}) (p, ω_{n}) = ∫_{ℝ} p(t)q_{n}(t) dμ(t) = 0.

According to Theorem 1, an n-point quadrature formula (3) with respect to the positive measure dμ(t) has the maximal algebraic degree of exactness 2n − 1, i.e., m = n is optimal. The higher m (> n) is impossible. Indeed, according to 2^{∘}, the case m = n + 1 requires the orthogonality (p, ω_{n}) = 0 for all p ∈ 𝒫_{n}, which is impossible when p = ω_{n}.

The cases m = n − 1 and m = n − 2 lead to the Gauss-Radau (one of the endpoints a or b is included in the set of nodes) and Gauss-Lobatto formulas (τ_{1} = a and τ_{n} = b), respectively.

3. Gauss-Christoffel quadrature rules for classical and non-classical weights

We start this section with monic polynomials

π_{k}(dμ; t) = π_{k}(t) = t^{k} + lower degree terms,

k = 0, 1, …, which are orthogonal with respect to the inner product (4). Because of the property (tf, g) = (f, tg), these polynomials satisfy the threeterm recurrence equation

π_{k+1}(t) = (t − α_{k})π_{k}(t) − β_{k}π_{k−1}(t),  k = 0, 1, 2, …,  (5)

π_{0}(t) = 0, π_{−1}(t) = 0,

where (α_{k}) = (α_{k}(dμ)) and (β_{k}) = (β_{k}(dμ)) are sequences of recursion coefficients. The coefficient β_{0} which is multiplied by π_{−1}(x) = 0 in (5) may be arbitrary. Sometimes, it is convenient to define it by β_{0} = μ_{0} = ∫_{ℝ} dμ(x). Then the norm of π_{n} can be expressed in the form ∥π_{n}∥ = \sqrt{(π_{n}, π_{n})} = \sqrt{β_{0}β_{1} ⋯ β_{n}}.

For generating Gauss-Christoffel quadrature rules there are numerical methods, which are computationally much better than a computation of nodes by using Newton’s method and then a direct application of Christoffel’s expressions for the weights (see e.g. Davis and Rabinowitz [18]). The characterization of the Gauss-Christoffel quadratures via an eigenvalue problem for the Jacobi matrix has become the basis of current methods for generating these quadratures. The most popular of them is one due to Golub and Welsch [43]. Their method is based on determining the eigenvalues and the first components of the eigenvectors of a symmetric tridiagonal Jacobi matrix.

Theorem 2. The nodes τ_{k} in the Gauss-Christoffel quadrature rule (3), with respect to a positive measure dμ, are the eigenvalues of the n-th order Jacobi matrix

J_{n}(dμ) = [ \begin{matrix}α_{0} & \sqrt{β_{1}} & & & O \\ \sqrt{β_{1}} & α_{1} & \sqrt{β_{2}} \\ & \sqrt{β_{2}} & α_{2} & ... \\ & & ... & ... & \sqrt{β_{n−1}} \\ O & & & \sqrt{β_{n−1}} & α_{n−1} \end{matrix} ] ,  (6)

where α_{ν} and β_{ν}, ν = 0, 1, …, n − 1, are the coefficients in the three-term recurrence relation (5) for the monic orthogonal polynomials π_{ν}(dμ; ·). The weights A_{k} are given by

A_{k} = β_{0}v_{k,1}^{2}, k = 1, …, n,

where β_{0} = μ_{0} = ∫_{ℝ} dμ(t) and v_{k,1} is the first component of the normalized eigenvector v_{k} corresponding to the eigenvalue x_{k},

J_{n}(dμ)v_{k} = x_{k}v_{k},  v_{k}^{T}v_{k} = 1,  k = 1, …, n.

Simplifying QR algorithm so that only the first components of the eigenvectors are computed, Golub and Welsch [43] gave an efficient procedure for constructing the Gaussian quadrature rules. This procedure was implemented in several programming packages including the most known ORTHPOL given by Gautschi [31] (see also our package OrthogonalPolynomials realized in Mathematica [13]).

According to Theorem 2, we need the recursion coefficients α_{k} and β_{k}, k ≤ m − 1, for the monic polynomials π_{ν}(dμ; ·), in order to construct the n-point Gauss-Christofell quadrature formula (3), with respect to a positive measure dμ, for each n ≤ m. In the case of the classical orthogonal polynomials, i.e., Jacobi, generalized Laguerre and Hermite polynomials (cf. [53, Section 2.3]), these coefficients are known explicitly^{1} and the construction problem of Gaussian quadratures is completely solved by Theorem 2. However, in the case of strong non-classical polynomials (see [53, Subsection 2.4.7]), we need an additional numerical construction of recursion coefficients (see [53, Subsection 2.4.8]).

In the sequel we mention some non-classical measures dμ(x) = w(x) dx for which the recursion coefficients α_{k}(dμ), β_{k}(dμ), k = 0, 1, …, n − 1, have been provided in the literature and used in the construction of Gaussian quadratures.

1^{∘} One-sided Hermite weight w(x) = exp(−x^{2}) on [0, c], 0 < c ≤ +∞. This distribution w(x) dx is known as the Maxwell (velocity) distribution. The cases c = 1, n = 10 and c = +∞, n = 15 were considered by Steen, Byrne and Gelbard [96] (see also Gautschi [30]).

2^{∘} Logarithmic weight w(x) = x^{α} log(1/x), μ > −1 on (0, 1). Piessens and Branders [79] considered cases when α = 0, ±1/2, ±1/3, −1/4, −1/5 (see also Gautschi [29]).

3^{∘} Airy weight w(x) = exp(−x^{3}/3) on (0, +∞). The inhomogeneous Airy functions Hi(x) and Gi(x), arise in theoretical chemistry (e.g. in harmonic oscillator models for large quantum numbers) and their integral representations (see Lee [50]) are given by

Hi(t) = \frac{1}{π} ∫_{0}^{+∞} w(x)e^{xt} dx,
Gi(t) = −\frac{1}{π} ∫_{0}^{+∞} w(x)e^{−xt/2} cos(\frac{\sqrt{3}}{2} xt + \frac{2π}{3} ) dx.

These functions can effectively be evaluated by the Gaussian quadrature relative to the Airy weight w(x). It needs orthogonal polynomials with

1Also, there are a few non-classical cases for which the recursion coefficients are known explicitly (see [53, Section 2.4]).

respect to this weight. Gautschi [27] computed the recursion coefficients for n = 15 with 16 decimal digits after the decimal point (D).

4^{∘} Reciprocal gamma function w(x) = 1/Γ(x) on (0, +∞). Gautschi [26] determined the recursion coefficients for n = 40 with 20 significant decimal digits (S). This function could be useful as a probability density function in reliability theory (see Fransén [21]).

5^{∘} Einstein’s and Fermi’s weight functions on (0, +∞),

w_{1}(x) = ε(x) = \frac{x}{e^{x} − 1}  and  w_{2}(x) = φ(x) = \frac{1}{e^{x} + 1} .

These functions arise in solid state physics. Integrals with respect to the measure dμ(x) = ε(x)^{r} dx, r = 1 and r = 2, are widely used in phonon statistics and lattice specific heats and occur also in the study of radiative recombination processes. Similarly, integrals with φ(x) are encountered in the dynamics of electrons in metals. For w_{1}(x), w_{2}(x), w_{3}(x) = ε(x)^{2} and w_{4}(x) = φ(x)^{2}, Gautschi and Milovanović [36] gave the first systematic investigation on the derivation of quadrature rules with high-precision, determined the recursion coefficients α_{k} and β_{k}, k < n = 40, and gave an application of the corresponding Gauss-Christoffel quadratures to the summation of slowly convergent series, whose general term is expressible in terms of a Laplace transform or its derivative. We call such a summation the method of Laplace transform.

In the numerical construction for the measure dμ(x) = [ε(x)]^{r} dx on (0, +∞), r ≥ 1, we used the discretized Stieltjes procedure based on the Gauss-Laguerre quadratures, so that

∫_{0}^{+∞} P (x) dμ(x) = \frac{1}{r} ∫_{0}^{+∞} P (x/r) (\frac{x/r}{1 − e^{−x/r}} )^{r} e^{−x} dx

≈ ∑_{k=1}^{N} \frac{A_{k}^{L}}{r} (\frac{x_{k}^{L}/r}{1 − e^{−x_{k}^{L}/r}} )^{r} P (x_{k}^{L}/r),

where P ∈ P and N ≫ n. The Gauss-Laguerre nodes x_{k}^{L} (zeros of the standard Laguerre polynomial L_{N}(x)) and the weights A_{k}^{L} can be easily computed for an arbitrary N by the Golub-Welsch algorithm.

6^{∘} The hyperbolic weights on (0, +∞),

w_{1}(x) = \frac{1}{cosh^{2} x}  and  w_{2}(x) = \frac{sinh x}{cosh^{2} x} .

The recursion coefficients α_{k}, β_{k}, for k < 40 were obtained by Milovanović [56]. The main application of these quadratures is the summation of slowly convergent series, with the general term a_{k} = f(k). Such a method known as the method of contour integration was given in [56] (for further applications see [57]).

If F is an integral of f and

Φ(x, y) = −\frac{1}{2}[F (x + iy) + F (x − iy)],

then, under certain conditions,

T_{m} = ∑_{k=m}^{+∞} a_{k} = ∫_{0}^{+∞} Φ(m − 1/2, t/π)w_{1}(t) dt

= ∑_{ν=1}^{n} A_{ν}Φ(m − 1/2, τ_{ν}/π) + R_{n}(Φ).

As an illustration, we consider only the series

T_{1}(a) = ∑_{k=1}^{+∞} \frac{1}{\sqrt{k}(k + a)} ,  (7)

which appeared (for a = 1) in a study of spirals (see Davis [17]) and defines the “Theodorus constant.” The first 1 000 000 terms of the series T_{1}(1) give the result 1.8580 …, i.e., T_{1}(1) ≈ 1.86 (only 3-digit accuracy). Using the method of Laplace transform, as a increases, the convergence of the Gauss quadrature formula (with Einstein’s weight) slows down considerably. For example, when a = 8, we get results with relative errors 1.4×10^{−1}, 2.3×10^{−2}, 1.5 × 10^{−3}, 1.9 × 10^{−4}, 2.5 × 10^{−5}, 2.1 × 10^{−6}, 2.5 × 10^{−7}, and 2.6 × 10^{−8}, for n = 5(5)40, respectively.

Now, we directly apply the method of contour integration to (7) with

F (z) = \frac{2}{\sqrt{a}} (arctan \sqrt{\frac{z}{a}} − \frac{π}{2} ) ,

where the integration constant is taken so that F (∞) = 0.

We can represent (7) in the form

T_{1}(a) = ∑_{k=1}^{m−1} \frac{1}{\sqrt{k}(k + a)} + T_{m}(a),  T_{m}(a) = ∑_{k=m}^{+∞} \frac{1}{\sqrt{k}(k + a)}, (8)

and then use the Gaussian quadrature formula (with the hyperbolic weight w_{1}) to calculate T_{m}(a). Relative errors in approximations for T_{1}(a), when a ≥ 1/2 and m = 4 are less than 10^{−18} with only n = 10 nodes. The method is very efficient. Moreover, its convergence is slightly faster if the parameter a is larger. Taking n = 25 nodes in our quadratures, we obtain the exact sum T_{1}(1) = 1.86002507922119030718069591572 … (to 30 significant digits).

An approach to summation formulas due to Plana, Lindelöf and Abel was recently given by Dahlquist [14, 15, 16].

Here we also mention an interesting application of Gaussian type formulas in the so-called moment-preserving spline approximation of a given function f on [0, +∞) (or on a finite interval, e.g. [0, 1]). Such kind of problems appeared in the physics literature, for example in the approximation of the Maxwell velocity distribution by a linear combination of Dirac δ-functions (see Laframboise and Stauffer [49]) or in the corresponding approximation by a linear combination of Heaviside step functions (see Calder and Laframboise [8]). The authors used some classical methods, which are very sensitive to rounding errors, and therefore their computations require a high-precision arithmetic. In order to get a stable method for this kind of approximation, Gautschi and Milovanović [38] found new applications of Gaussian type of quadratures.

Let f be a given function defined on the positive real line ℝ_{+} = [0, +∞) and s_{n,m} be a spline of the form

s_{n,m}(t) = ∑_{ν=1}^{n} a_{ν}(t_{ν} − t)_{+}^{m}, 0 ≤ t < +∞,  (9)

where the plus sign on the right is the cutoff symbol, u_{+} = u if u > 0 and u_{+} = 0 if u ≤ 0, 0 < t_{1} < ⋯ < t_{n}, a_{ν} ∈ ℝ. We considered the moment-preserving spline approximation f(t) ≈ s_{n,m}(t) such that

∫_{0}^{+∞} s_{n,m}(t)t^{j} dV = ∫_{0}^{+∞} f(t)t^{j} dV, j = 0, 1, …, 2n − 1,

where dV is the volume element depending on the geometry of the problem. In some concrete applications in physics, up to unimportant numerical factors, dV = t^{d−1} dt, where d = 1, 2, and 3 for rectilinear, cylindric, and spherical geometry, respectively.

For fixed n, m ∈ ℕ, d ∈ {1, 2, 3} and certain conditions on f we proved in [38] that the spline function s_{n,m} exists uniquely if and only if the measure

dλ_{m}(t) = \frac{(−1)^{m+1}}{m!} t^{m+d}f^{(m+1)}(t) dt on [0, +∞)

admits an n-point Gauss-Christoffel quadrature formula

∫_{0}^{+∞} g(x) dλ_{m}(x) = ∑_{ν=1}^{n} λ_{ν}^{(n)}g(τ_{ν}^{(n)}) + R_{n}(g; dλ_{m}),  (10)

with distinct positive nodes τ_{ν}^{(n)}, where R_{n}(g; dλ_{m}) = 0 for all g ∈ 𝒫_{2n−1}. In that event, the knots t_{ν} and weights a_{ν} in (9) are given by

t_{ν} = τ_{ν}^{(n)}, a_{ν} = t_{ν}^{−(m+d)}λ_{ν}^{(n)}, ν = 1, …, n.

The approximation error e_{n,m}(t) := f(t) − s_{n,m}(t) can be expressed by the remainder term of (10), i.e., e_{n,m}(t) = R_{n}(σ_{t}; dλ_{m}), where σ_{t}(x) = x^{−(m+d)}(x − t)_{+}^{m}, x, t > 0.

Note that for m = 0, (9) reduces to linear combination of Heaviside step functions. The case with Dirac δ-functions can be formally obtained putting m = −1.

Approximation on a compact interval was considered by Frontini, Gautschi and Milovanović [22].

On the end of this section we mention that there are several other extensions and generalizations, e.g. the so-called optimal quadratures, Radau and Lobatto quadratures, Kronrod quadratures, etc. An account of the role played by moments and modified moments in the construction of interpolatory quadrature rules, especially weighted Newton-Cotes and Gaussian rules, is given by Gautschi [33].

4. Quadratures with multiple node

More than hundred years after the famous Gauss method of approximate integration, there appeared the idea of numerical integration involving multiple nodes Chakalov [9, 10, 11], Turán [99], Popoviciu [81], Ghizzetti and Ossicini [40, 41], etc.

Let η_{1}, …, η_{m} (η_{1} < ⋯ < η_{m}) be given fixed (or prescribed) nodes, with multiplicities m_{1}, …, m_{m}, respectively, and τ_{1}, …, τ_{n} (τ_{1} < ⋯ < τ_{n}) be free nodes, with given multiplicities n_{1}, …, n_{n}, respectively. Interpolation quadrature formulae of a general form

I(f) = ∫_{ℝ} f(t) dμ(t) ≅ Q(f),

where

Q(f) = ∑_{ν=1}^{n} ∑_{i=0}^{n_{ν}−1} A_{i,ν}f^{(i)}(τ_{ν}) + ∑_{ν=1}^{m} ∑_{i=0}^{m_{ν}−1} B_{i,ν}f^{(i)}(η_{ν}),  (11)

with an algebraic degree of exactness at least M + N − 1, were investigated by D. D. Stancu [88, 91, 98].

Using fixed and free nodes we introduce two polynomials

q_{M}(t): = ∏_{ν=1}^{m}(t − η_{ν})^{m_{ν}} and Q_{N}(t): = ∏_{ν=1}^{n}(t − τ_{ν})^{n_{ν}},  (12)

where M = ∑_{ν=1}^{m} m_{ν} and N = ∑_{ν=1}^{n} n_{ν}. Choosing the free nodes to increase the degree of exactness leads to so-called Gaussian type of quadratures. If the free (or Gaussian) nodes τ_{1}, …, τ_{n} are such that I(f) = Q(f) for each f ∈ 𝒫_{M+N+n−1}, the corresponding quadrature Q we call the Gauss-Stancu formula. D.D. Stancu [93] proved that τ_{1}, …, τ_{n} are the Gaussian nodes if and only if

∫_{ℝ} t^{k}Q_{N}(t)q_{M}(t) dμ(t) = 0  (13)

for k = 0, 1, …, n − 1. Under some restrictions of node polynomials q_{M}(t) and Q_{N}(t) on the support interval of the measure dμ(t) we can give sufficient conditions for the existence of Gaussian nodes (cf. Stancu [93] and [42]). For example, if the multiplicities of the Gaussian nodes are odd, e.g., n_{ν} = 2s_{ν} + 1, ν = 1, …, n, and if the polynomial with fixed nodes q_{M}(t) does not change its sign in the support interval of the measure dμ(t), then, in this interval, there exist real distinct nodes τ_{ν}, ν = 1, …, n.

The last condition for the polynomial q_{M}(t) means that the multiplicities of the internal fixed nodes must be even. Defining a new (nonnegative) measure dμ̂(t) by dμ̂(t) = γq_{M}(t) dμ(t) (γ = sgn (q_{M}(t)), the “orthogonality conditions” (13) can be expressed in the simpler form

∫_{ℝ} t^{k}Q_{N}(t) dμ̂(t) = 0, k = 0, 1, …, n − 1.

This means that the general quadrature problem (11), under these conditions, can be reduced to a problem with only Gaussian nodes, but with respect to another modified measure. Computational methods for this purpose are based on Christoffel’s theorem and described in details in [32] (see also [42] and [47]).

Let π_{n}(t): = ∏_{ν=1}^{n}(t − τ_{ν}). Since Q_{N}(t)/π_{n}(t) = ∏_{ν=1}^{n}(t − τ_{ν})^{2s_{ν}} ≥ 0 over the support interval, we can make an additional reinterpretation of the “orthogonality conditions” (13) in the form

∫_{ℝ} t^{k}π_{n}(t) dμ(t) = 0, k = 0, 1, …, n − 1,  (14)

where

dμ(t) = ( ∏_{ν=1}^{n}(t − τ_{ν})^{2s_{ν}} )dλ̂(t).

This means that π_{n}(t) is a polynomial orthogonal with respect to the new nonnegative measure dμ(t) and, therefore, all zeros τ_{1}, …, τ_{n} are simple, real, and belong to the support interval. As we see the measure dμ(t) involves the nodes τ_{1}, …, τ_{n}, i.e., the unknown polynomial π_{n}(t), which is implicitly defined (see Engels [20, pp. 214–226]). This polynomial π_{n}(t) belongs to the class of so-called σ-orthogonal polynomials {π_{n,σ}(t)}_{n∈ℕ_{0}}, which correspond to the sequence σ = (s_{1}, s_{2}, …) connected with multiplicities of Gaussian nodes. Namely, π_{n}(t) = π_{n,σ}(t). If σ = (s, s, …), these polynomials reduce to the s-orthogonal polynomials. (For details see Milovanović [59].)

Quadratures with only Gaussian nodes (m = 0),

∫_{ℝ} f(t) dλ(t) = ∑_{ν=1}^{n} ∑_{i=0}^{2s_{ν}} A_{i,ν}f^{(i)}(τ_{ν}) + R(f),  (15)

which are exact for all algebraic polynomials of degree at most d_{max} = 2 ∑_{ν=1}^{n} s_{ν} + 2n − 1, are known as Chakalov-Popoviciu quadrature formulas (see [9, 10, 11], [81]). A deep theoretical progress in this subject was made by Stancu (see [98] and [86]–[95]). In the special case of the Legendre measure on [−1, 1], when all multiplicities are mutually equal, these formulas reduce to the well-known Turán quadrature [99]. The case with a weight function dμ(t) = w(t)dt on [a, b] has been investigated by Italian mathematicians Ossicini, Ghizzetti, Guerra, Rosati, and also by Chakalov, Stroud, Stancu, Ionescu, Pavel, etc. (see Milovanović [59] for references).

At the Third Conference on Numerical Methods and Approximation Theory (Niš, August 18–21, 1987) (see [54]) we presented a stable method with quadratic convergence for numerically constructing s-orthogonal polynomials, which zeros are nodes of Turán quadratures. The basic idea for our method to numerically construct s-orthogonal polynomials with respect to the measure dμ(t) on the real line ℝ is a reinterpretation of the s-orthogonality in terms of implicitly defined standard orthogonality. Further progress in this direction was made by Gautschi and Milovanović [39]. After our survey [59] (published in series “Numerical Analysis in the 20th Century”), where we gave a connection between quadratures, s and σorthogonality and moment-preserving approximation with defective splines, the interest for this subject rapidly increases (Gautschi, Yang, Shi, Yang & Wang, Gout & Guessab, Shi & Xu, etc.). A very efficient method for constructing quadratures with multiple nodes was given recently by Milovanović, Spalević and Cvetković [74].

The error analysis and estimates of the remainder term in quadratures on [−1, 1] with multiple nodes have been treated recently in the class of analytic functions in the interior of a simple closed curve Γ in the complex plane surrounding the interval [−1, 1]. Two choices of the contour Γ have been widely used: 1^{∘} a circle with center at origin and radius ϱ (> 1), and 2^{∘} an ellipse with foci at ±1. In several papers Milovanović and Spalević (see [68, 69, 70, 71, 72, 73]) developed a few methods for error estimates and obtained very exact L^{1}- and L^{∞}-estimates of the remainder term for the generalized Chebyshev weights.

Another type of quadratures with multiple nodes are the so-called Birkhoff quadratures. Roughly speaking their quadrature sums do not include all derivatives. We mention here only a very special case of Birkhoff quadratures – the generalized (0, m) quadrature problem

∫_{ℝ} f(x) dμ(x) = ∑_{k=1}^{n}[v_{k}f(x_{k}) + w_{k}f^{(m)}(x_{k})] + R(f)  (16)

of highest degree of precision, which first was stated in 1974 by Turán for dμ(x) = dx on [−1, 1], m = 2, and nodes as zeros of the polynomial Π_{n}(x) := (1−x^{2})P′_{n−1}(x), where P_{k} is the Legendre polynomial of degree k (cf. [100]). For some particular results on (16) see [101], [19], [75], [51].

We mention also two nice recent books by Shi [83, 84].

Further extensions dealing with quadratures with multiple nodes for ET (Extended Tschebycheff) systems were given by [45], [2], [4], [3], etc.

5. Orthogonality with respect to a moment functional and corresponding quadratures

In the previous sections the inner product was always positive definite provided the existence of the corresponding orthogonal polynomials with real zeros in the support of the measure. Such zeros appeared as the nodes of the Gaussian formulas. However, there are a more general concept of orthogonality with respect to a given linear moment functional L on the linear space 𝒫 of all algebraic polynomials. Because of linearity, the value of the linear functional L at every polynomial is known if the values of L are known at the set of all monomials, i.e., if we know L(x^{k}) = μ_{k}, for each k ∈ ℕ_{0}. In that case we can introduce a system of orthogonal polynomials {π_{k}}_{k=0∈ℕ_{0}} with respect to the functional L if for all nonnegative integers k and n (cf. Chihara [12, p. 7]),

• π_{k}(x) is polynomial of degree k,

• L(π_{k}(x)π_{n}(x)) = 0, if k ≠ n,

• L(π_{n}^{2}(x)) ≠ 0.

If the sequence of orthogonal polynomials exists for a given linear functional L, then L is called quasi-definite or regular linear functional. Under the condition L(π_{n}^{2}(x)) > 0, n ∈ ℕ_{0}, the functional L is called positive definite. The necessary and sufficient conditions for the existence of a sequence of orthogonal polynomials with respect to the linear functional L are that for each n ∈ ℕ the Hankel determinants

Δ_{n} = | μ_{0} μ_{1} … μ_{n−1} ; μ_{1} μ_{2}  μ_{n} ; ⋮ ; μ_{n−1} μ_{n}  μ_{2n−2} | ≠ 0.  (17)

This concept of orthogonality enables us to consider also the quadratures of Gaussian type with respect to the functional L, even in the cases when the weight is a complex function. In the next subsections we mention only the two such cases.

5.1. Orthogonality on the semicircle and quadratures. Let w be a weight function which is positive and integrable on the open interval (−1, 1), though possibly singular at the endpoints, and which can be extended to a function w(z) holomorphic in the half disc D_{+} = {z ∈ ℂ : |z| < 1, Im z > 0}. Consider the following “inner product”

(f, g) = ∫_{Γ} f(z)g(z)w(z)(iz)^{−1} dz = ∫_{0}^{π} f(e^{iθ})g(e^{iθ})w(e^{iθ}) dθ,  (18)

where Γ is the circular part of ∂D_{+} and all integrals are assumed to exist (possibly) as appropriately defined improper integrals. The “inner product” is not Hermitian; we deliberately did not conjugate the second factor g and did not integrate with respect to the measure |w(e^{iθ})| dθ. The existence of the corresponding orthogonal polynomials {π_{n}}_{n∈ℕ_{0}}, therefore, is not guaranteed.

The case w = 1 was considered by Gautschi and Milovanović [37]. The existence and uniqueness of polynomials orthogonal on the unit semicircle were proved via moment determinants.

A more general case of the complex weight was considered by Gautschi, Landau and Milovanović [35]. Under the condition

Re (1, 1) = Re ∫_{0}^{π} w(e^{iθ}) dθ ≠ 0,

they proved that the orthogonal polynomials {π_{n}}_{n∈ℕ_{0}} exist uniquely and they can be represented in the form π_{n}(z) = p_{n}(z) − iθ_{n−1}p_{n−1}(z), where

θ_{n−1} = θ_{n−1}(w) = \frac{μ_{0}p_{n}(0) + iq_{n}(0)}{iμ_{0}p_{n−1}(0) − q_{n−1}(0)}.

Here, p_{k} are standard (real) polynomials orthogonal with respect to the inner product

[f, g] = ∫_{−1}^{1} f(x)\overline{g(x)}w(x) dx,

and q_{k} are the corresponding associated polynomials of the second kind,

q_{k}(z) = ∫_{−1}^{1} \frac{p_{k}(z) − p_{k}(x)}{z − x}w(x) dx.

In the previous mentioned papers, we obtained several results for the polynomials π_{n}, including the three-term recurrence relation, zero distribution, and a linear differential equation. Also, we gave a few applications in numerical differentiation and numerical integration.

In [55] the Gegenbauer case w(z) = (1 − z^{2})^{λ−1/2}, λ > −1/2, on [−1, 1], was investigated, for which μ_{0} = π. The explicit expressions for θ_{n} in terms of gamma function enabled us to propose a stable construction of Gauss-Christoffel quadratures on the semicircle

∫_{0}^{π} f(e^{iθ})w(e^{iθ}) dθ = ∑_{ν=1}^{n} σ_{ν}f(ζ_{ν}) + R_{n}(f).

The nodes ζ_{ν} = ζ_{ν}^{(n)} are precisely the zeros of π_{n}(z). In order to determine the weights σ_{ν} = σ_{ν}^{(n)}, we consider an eigenvalue problem for the Jacobi matrix

J_{n} = [ \begin{matrix}iα_{0} & θ_{0} & & & O \\ θ_{0} & iα_{1} & θ_{1} \\ & θ_{1} & iα_{2} & ... \\ & & ... & ... & θ_{n−2} \\ O & & & θ_{n−2} & iα_{n−1} \end{matrix} ] .

Using some matrix transformations (see [55]) we reduce this complex matrix to a real nonsymmetric tridiagonal matrix given by

A_{n} = −iD_{n}^{−1}J_{n}D_{n} = [ \begin{matrix}α_{0} & θ_{0} & & & O \\ −θ_{0} & α_{1} & θ_{1} \\ & −θ_{1} & α_{2} & ... \\ & & ... & ... & θ_{n−2} \\ O & & & −θ_{n−2} & α_{n−1} \end{matrix} ] ,

whose eigenvalues are η_{ν} = −iζ_{ν}. Denote by V_{n} the matrix of the eigenvectors of this matrix A_{n}, each normalized so that the first component is equal to 1. Putting v = [σ_{1} σ_{2} … σ_{n}]^{T}, we can get a system of linear equations for finding the weight coefficients

V_{n}v = πe_{1},

where e_{1} = [1 0 … 0]^{T} (the first coordinate vector).

The zeros ζ_{ν} of π_{n} are located symmetrically with respect to the imaginary axis and for n ≥ 2 all zeros are contained in D_{+}. For the general Jacobi weight w(z) = (1 − z)^{α}(1 + z)^{β}, α, β > −1, we conjectured that all zeros are also contained in the upper unit half disc D_{+} (see [35]). This was verified numerically for many values of parameters α and β.

5.2. Orthogonal polynomials for oscillatory weights. Let w be a given weight function on [−1, 1] and dμ(x) = xw(x)e^{iζx} dx, where ζ ∈ ℝ. In this subsection we consider the existence of the orthogonal polynomials π_{n} with respect to the functional

L(p) = ∫_{−1}^{1} p(x) dμ(x), μ_{k} = L(x^{k}), k ∈ ℕ_{0}.  (19)

Two cases are intensively studied.

1^{∘} Case w(x) = 1, ζ = mπ, m ∈ ℤ \ {0}.

The existence and uniqueness for polynomials π_{n} were proved via moment determinants (see [61]). Namely, for each integer m ≠ 0, the orthogonal polynomials with respect to the measure dμ(x) = xe^{imπx} dx on −1, 1] exist uniquely and satisfy the three-term recurrence relation

π_{n+1}(x) = (x − iα_{n})π_{n}(x) − β_{n}π_{n−1}(x), n ≥ 0,

with π_{0}(x) = 1 and π_{−1}(x) = 0.

From the Hankel determinants it is clear that the recursion coefficients are rational functions in ζ = mπ. Using our software package [13] we can generate coefficients even in symbolic form for some reasonable values of n (e.g. n ≤ 20) and state the following conjecture: Let a_{n}(z) and c_{n}(z) be algebraic polynomials with integer coefficients of degree r_{n} and s_{n}, respectively, i.e., a_{n}(z) = A_{n}z^{r_{n}} +⋯ and c_{n}(z) = z^{s_{n}} +· · ·. If ζ = mπ and n ≥ 2, then

α_{n} = \frac{a_{n}(ζ^{2})}{ζ c_{n−1}(ζ^{2})c_{n}(ζ^{2})} , β_{n} = B_{n} \frac{c_{n−2}(ζ^{2})c_{n}(ζ^{2})}{ζ^{2} c_{n−1}(ζ^{2})^{2}} ,

where

A_{n} = \{ \begin{matrix}−\frac{n^{2} − 1}{4} & (n odd), \\ \frac{n^{2} + 10n + 8}{4} & (n even), \end{matrix} B_{n} = \{ \begin{matrix}1 & (n odd), \\ −n^{2} & (n even), \end{matrix}

and

r_{n} = \frac{n(n + 1)}{2} , s_{n} = \{ \begin{matrix}\frac{(n + 1)^{2}}{4} & (n odd), \\ \frac{n(n + 2)}{4} & (n even). \end{matrix}

The complexity of expressions for α_{n} and β_{n} increases exponentially with n. On the other side, there is an efficient algorithm for their numerical construction [61].

The corresponding Gaussian quadratures can also be constructed. A possible application of these quadratures is in numerical calculation of integrals involving highly oscillatory integrands, in particular for calculation of the Fourier coefficients,

C_{m}(f) + iS_{m}(f) = ∫_{−1}^{1} f(x)e^{imπx} dx = ∫_{−1}^{1} \frac{f(x) − f(0)}{x} dμ(x)

≈ ∑_{ν=1}^{n} \frac{A_{ν}^{n}}{x_{ν}^{n}} (f(x_{ν}^{n}) − f(0)),

where x_{ν}^{n} are zeros of π_{n}, and A_{k}^{n} are the corresponding weight coefficients. Taking f(x) = \frac{x}{x^{2} + 1/4} we have

S_{m}(f) = ∫_{−1}^{1} \frac{x}{x^{2} + 1/4} sin(mπx) dx ≈ G_{n}(f) = Im \{∑_{ν=1}^{n} \frac{A_{ν}^{n}}{(x_{ν}^{n})^{2} + 1/4}\} .

For example, for m = 100 and n = 10 we get

S_{100}(f) ≈ G_{10}(f) = −0.0050929580138121841037438653708.

If we increase m, the convergence of G_{n}(f) is rather faster, e.g., for m = 10^{6}, the relative error in G_{10}(f) is smaller than 10^{−60}. The integrand for m = 200 is presented in Figure .

Figure 1. Integrand for m = 200

2^{∘} Case w(x) = 1/\sqrt{1 − x^{2}}, ζ ∈ ℝ. This case was considered recently by Milovanović and Cvetković [62]. We showed that the corresponding moments can be expressed in terms of Bessel functions J_{0} and J_{1} as

μ_{k} = \frac{iπ}{(iζ)^{k}} (P_{k}(ζ^{2})J_{1}(ζ) + ζQ_{k}(ζ^{2})J_{0}(ζ)), k ∈ ℕ_{0},

where P_{k} and Q_{k} are algebraic polynomials with integer coefficients of degrees 2[k/2] and 2[(k − 1)/2], respectively, and satisfy the following recurrence relation

y_{k+2} = −(k + 2)y_{k+1} − ζ^{2}y_{k} − (k + 1)ζ^{2}y_{k−1},

with initial conditions

P_{0} = 1, P_{1} = −1, P_{2} = 2 − ζ^{2},
Q_{0} = 0, Q_{1} = 1, Q_{2} = −1.

When ζ is a positive zero of J_{0}, we proved that polynomials π_{n} orthogonal with respect to the functional (19) exist uniquely and satisfy the three-term recurrence relation, for which coefficients we found the asymptotic formulae. Moreover, we proved that the essential spectrum of the associated Jacobi operator created using three-term recurrence coefficients α_{k} and β_{k}, k ∈ ℕ_{0}, is [−1, 1].

6. Nonstandard quadratures of Gaussian type

All previous quadrature rules use the information on the integrand only at some selected points x_{k}, k = 1, …, n (the values of the function f and its derivatives in the cases of rules with multiple nodes). Such quadratures will be called the standard quadrature formulae. However, in many cases in physics and technics it is not possible to measure the exact value of the function f at points x_{k}, so that a standard quadrature cannot be applied. On the other side, some other information on f can be available, as

1^{∘} the averages

\frac{1}{2h_{k}} ∫_{I_{k}} f(x) dx

of this function over some non-overlapping subintervals I_{k}, with length of I_{k} equals 2h_{k}, and their union which is a proper subset of supp (dμ);

2^{∘} a fixed linear combination of the function values, e.g.

af(x − h) + bf(x) + cf(x + h)

at some points x_{k}, where a, b, c are constants and h is sufficiently small positive number, etc.

The last example is gained from communications. Namely, supposing that we are receiving a signal with an interference, then the measurements provide linear combinations of the function values rather than single function values. In another words, for signals (functions) which depend on time in a discrete sense, averaging can be given with respect to a discrete measure, with some jumps at x_{k} − h, x_{k}, x_{k} + h, in our simple example. The same problem appears when a signal f passes through a digital filter.

Suppose now that we want to find the value of the integral of the signal f with respect to some positive measure μ, i.e., to calculate ∫_{ℝ} f(x) dμ(x), and that we know only values of the signal after passing through some digital filter, i.e., we know the values g_{k} = af(x_{k} − h) + bf(x_{k}) + cf(x_{k} + h), k = 1, …, n.

In this very simple example, using the initial conditions, in general, we can state a system of linear equations and calculate the values f(x_{k}) from the values g_{k}, and then we apply the quadrature rule. Unfortunately, very often the system of linear equations is ill-conditioned and such a procedure significantly disturbs the final result of computations. Therefore, it is much better to use directly the values g_{k} instead of f(x_{k}) if it is possible. This idea leads to the so-called nonstandard quadratures (see [67]).

Thus, if the information data {f(x_{k})}_{k=1}^{n} in the standard quadrature is replaced by {(A^{h_{k}}f)(x_{k})}_{k=1}^{n}, where A^{h} is an extension of some linear operator A^{h} : P → P, h ≥ 0, we get a non-standard quadrature formula

∫_{ℝ} f(x) dμ(x) = ∑_{k=1}^{n} w_{k}(A^{h_{k}}f)(x_{k}) + R_{n}(f).  (20)

This kind of quadratures is based on operator values for a general family of linear operators A^{h}, acting of the space of algebraic polynomials, such that the degrees of polynomials are preserved.

As a typical example for such operators is the average (Steklov) operator mentioned before in 1^{∘}, i.e.,

(A^{h}p)(x) = \frac{1}{2h} ∫_{x−h}^{x+h} p(x) dx,  h > 0, p ∈ P.  (21)

The first idea on so-called interval quadratures, which are an example of non-standard quadrature rules, appeared a few decades ago. In 1976 Omladič, Pahor and Suhadolc [78] considered quadratures with the average operator (21) (see also Pitnauer and Reimer [80]). Some further investigations were given by Kuz’mina [48], Sharipov [82], Babenko [1], and Motornyi [76].

Let h_{1}, …, h_{n} be nonnegative numbers such that

a < x_{1} − h_{1} ≤ x_{1} + h_{1} < x_{2} − h_{2} ≤ x_{2} + h_{2} < ⋯ < x_{n} − h_{n} ≤ x_{n} + h_{n} < b,

(22)

and let w(x) be a given weight function on [a, b]. Using the previous inequalities it is obvious that we have 2(h_{1} + ⋯ + h_{n}) < b − a.

Recently, Bojanov and Petrov [5] proved that the Gaussian interval quadrature rule of the maximal algebraic degree of exactness 2n − 1 exists, i.e.,

∫_{a}^{b} f(x)w(x) dx = ∑_{k=1}^{n} \frac{w_{k}}{2h_{k}} ∫_{x_{k}−h_{k}}^{x_{k}+h_{k}} f(x)w(x) dx + R_{n}(f),  (23)

where R_{n}(f) = 0 for each f ∈ P_{2n−1}. Under conditions h_{k} = h, 1 ≤ k ≤ n, they also proved the uniqueness of (23). Moreover, in [6] Bojanov and Petrov proved the uniqueness of (23) for the Legendre weight (w(x) = 1) for any set of lengths h_{k} ≥ 0, k = 1, …, n, satisfying the condition (22). The question of the existence for bounded [a, b] is proved in [5] in much broader context for a given Chebyshev system of functions.

Recently in [60], using properties of the topological degree of non-linear mappings, it was proved that Gaussian interval quadrature formula is unique for the Jacobi weight function w(x) = (1−x)^{α}(1+x)^{β}, α, β > −1, on [−1, 1] and an algorithm for numerical construction was proposed. For the special case of the Chebyshev weight of the first kind and the special set of lengths an analytic solution can be given ([60]). Interval quadrature rules of Gauss-Radau and Gauss-Lobatto type with respect to the Jacobi weight functions are considered in [?].

Recently, Bojanov and Petrov [7] proved the existence and uniqueness of the weighted Gaussian interval quadrature formula for a given system of continuously differentiable functions, which constitute an ET system of order two on [a, b].

The cases with interval quadratures on unbounded intervals with the classical generalized Laguerre and Hermite weights have been recently investigated by Milovanović and Cvetković in [63] and [66].

In meantime we considered the nonstandard quadratures (20) with some special operators of the form (see [67])

(A^{h}p)(x) = \frac{1}{2h} ∫_{x−h}^{x+h} p(t) dt,

(A^{h}p)(x) = ∑_{k=−m}^{m} a_{k}p(x + kh) or (A^{h}p)(x) = ∑_{k=−m}^{m−1} a_{k}p (x + (k + \frac{1}{2})h) ,

and

(A^{h}p)(x) = ∑_{k=0}^{m} \frac{b_{k}h^{k}}{k!} 𝒟^{k}p(x),

where m is a fixed natural number and 𝒟^{k} = d^{k}/dx^{k}, k ∈ ℕ_{0}.

7. Gaussian quadrature for Müntz systems

Gaussian integration can be extended in a natural way to non-polynomial functions, taking a system of linearly independent functions

{P_{0}(x), P_{1}(x), P_{2}(x), …}  (x ∈ [a, b]),  (24)

usually chosen to be complete in some suitable space of functions. If dσ(x) is a given nonnegative measure on [a, b] and the quadrature rule

∫_{a}^{b} f(x)dσ(x) = ∑_{k=1}^{n} A_{k}f(x_{k}) + R_{n}(f)  (25)

is such that it integrates exactly the first 2n functions in (24), we call the rule (25) Gaussian with respect to the system (24). The existence and uniqueness of a Gaussian quadrature rule (25) with respect to the system (24), or shorter a generalized Gaussian formula, is always guaranteed if the first 2n functions of this system constitute a Chebyshev system on [a, b]. Then, all the weights A_{1}, …, A_{n} in (25) are positive. In the terms of moment spaces, the Gaussian rule corresponds to the unique lower principal representation of the measure dσ(x) (see Karlin and Studden [46]).

The generalized Gaussian quadratures for Müntz systems goes back to Stieltjes [97] in 1884. Taking P_{k}(x) = x^{λ_{k}} on [a, b] = [0, 1], where 0 ≤ λ_{0} < λ_{1} < ⋯, Stieltjes showed the existence of Gaussian formulae. In his short note he considered also Gauss-Radau formulae.

A numerical algorithm for constructing the generalized Gaussian quadratures was recently investigated by Ma, Rokhlin and Wandzura [52]. They take a Chebyshev system of functions {P_{0}, P_{1}, …, P_{2n−1}} on [a, b] with certain additional properties (Extended Hermite (EH) system). Their procedure requires the construction of the functions

ξ_{i}(x) = ∑_{j=1}^{2n} α_{ij}P_{j−1}(x),  η_{i}(x) = ∑_{j=1}^{2n} β_{ij}P_{j−1}(x)  (i = 1, …, n), (26)

such that

\{ \begin{matrix}ξ_{i}(x_{k}) = 0, \\ ξ′_{i}(x_{k}) = δ_{ik}, \end{matrix}  \{ \begin{matrix}η_{i}(x_{k}) = δ_{ik}, \\ η′_{i}(x_{k}) = 0, \end{matrix}  (27)

for all i = 1, …, n and all k = 1, …, n. The algorithm is ill conditioned (see [52, Remark 6.2]). In order to obtain the double precision results (REAL*8), the authors have performed the computations in extended precision (Qarithmetic – REAL*16) for generating Gaussian quadratures up to order 20, and in Mathematica (120 digits operations) for generating Gaussian quadratures of higher orders (n ≤ 40). In particular, they considered the following important cases of EH systems:

{1, log x, x, x log x, …, x^{n−1}, x^{n−1} log x}  (28)

and

{1, x^{s}, x, x^{s+1}, …, x^{n−1}, x^{s+n−1}}  (29)

for s = 1/3, s = −1/3, s = −2/3.

In [64], Milovanović and Cvetković presented an alternatively numerical method for constructing the generalized Gaussian quadrature (25) for Müntz polynomials, which is exact for each

f ∈ M_{2n−1}(Λ) = span{x^{λ_{0}}, x^{λ_{1}}, …, x^{λ_{2n−1}}}.

Beside the properties of orthogonal Müntz polynomials on (0, 1) and their connection with orthogonal rational functions, we gave also a method for numerical evaluation of such generalized polynomials (see [58]). Our method is rather stable and simpler than the previous one, since it is based on orthogonal Müntz systems. It performs calculations in double precision arithmetics in order to get double precision results.

As it is well-known the Gaussian quadrature rule is unique provided the measure σ has nonnegative absolutely continuous part and has finitely many atoms on [0, 1]. Such quadratures possess several properties of the classical Gaussian formulae (for polynomial systems), such as positivity of the weights, rapid convergence, etc. They can be applied to the wide class of functions, including smooth functions, as well as functions with endpoint singularities, such as those appeared in the boundary-contact value problems, integral equations, complex analysis, potential theory, and several other fields (see [52]).

REFERENCES

[1] V. F. B a b e n k o, On a certain problem of optimal integration, In: Studies on Contemporary Problems of Integration and Approximation of Functions and Their Applications, Collection of Research Papers, pp. 3–13, Dnepropetrovsk State University, Dnepropetrovsk, 1984. (Russian)

[2] D. L. B a r r o w, On multiple node Gaussian quadrature formulae, Math. Comp. 32 (1978), 431–439.

[3] B. B o j a n o v, Gaussian quadrature formulae for Tchebycheff systems, East J. Approx. 3 (1997), 71–88.

[4] B. D. B o j a n o v, D. B r a e s s, N. D y n, Generalized Gaussian quadrature formulas, J. Approx. Theory 48 (1986), 335–353.

[5] B. B o j a n o v, P. P e t r o v, Gaussian interval quadrature formula, Numer. Math. 87 (2001), 625–643.

[6] B. B o j a n o v, P. P e t r o v, Uniqueness of the Gaussian interval quadrature formula, Numer. Math. 95 (2003), 53–62.

[7] B. B o j a n o v, P. P e t r o v, Gaussian interval quadrature formula for Tchebycheff systems, SIAM J. Numer. Anal. 43 (2005), 787–795.

[8] A. C. C a l d e r, J. G. L a f r a m b o i s e, Multiple-waterbag simulation of inhomogeneous plasma motion near an electrode, J. Comput. Phys. 65 (1986), 18– 45.

[9] L. C h a k a l o v, Über eine allgemeine Quadraturformel, C. R. Acad. Bulgar. Sci. 1 (1948), 9–12.

[10] L. C h a k a l o v, General quadrature formulae of Gaussian type, Bulgar. Akad. Nauk Izv. Mat. Inst. 1 (1954), 67–84 (Bulgarian); English transl.: East J. Approx. 1 (1995), 261–276.

[11] L. C h a k a l o v, Formules générales de quadrature mécanique du type de Gauss, Colloq. Math. 5 (1957), 69–73.

[12] T. S. C h i h a r a, An Introduction to Orthogonal Polynomials, Gordon and Breach, New York, 1978.

[13] A. S. C v e t k o v i ć, G. V. M i l o v a n o v i ć, The Mathematica Package “OrthogonalPolynomials”, Facta Univ. Ser. Math. Inform. 19 (2004), 17-36.

[14] G. D a h l q u i s t, On summation formulas due to Plana, Lindelöf and Abel, and related GaussChristoffel rules, I, BIT 37 (1997), 256–295.

[15] G. D a h l q u i s t, On summation formulas due to Plana, Lindelöf and Abel, and related GaussChristoffel rules, II, BIT 37 (1997), 804-832.

[16] G. D a h l q u i s t, On summation formulas due to Plana, Lindelöf and Abel, and related GaussChristoffel Rules, III, BIT 39 (1999), 51-78.

[17] P. J. D a v i s, Spirals: from Theodorus to Chaos (With contributions by Walter Gautschi and Arieh Iserles), A.K. Peters, Wellesley, MA, 1993.

[18] P. J. D a v i s, P. R a b i n o w i t z, Methods of Numerical Integration (2nd edn.), Computer Science and Applied Mathematics, Academic Press Inc., Orlando, FL, 1984.

[19] D. K. D i m i t r o v, On a problem of Turń: (0, 2) quadrature formula with a high algebraic degree of precision, Aequationes Math. 41 (1991), 168–171.

[20] H. E n g e l s, Numerical Quadrature and Cubature, Academic Press, London, 1980.

[21] A. F r a n s é n, Accurate determination of the inverse gamma integral, BIT 19 (1979), 137–138.

[22] M. F r o n t i n i, W. G a u t s c h i, G. V. M i l o v a n o v i ć, Moment-preserving spline approximation on finite intervals, Numer. Math. 50 (1987), 503–518.

[23] C. F. G a u s s, Methodus nova integralium valores per approximationem inveniendi, Commentationes Societatis Regiae Scientarium Göttingensis Recentiores 3, 1814 [Werke III, pp. 163–196].

[24] W. G a u t s c h i, A survey of Gauss-Christoffel quadrature formulae, In: E.B. Christoffel – The Influence of his Work on Mathematics and the Physical Sciences (P. L. Butzer, F. Fehér, eds.), Birkhäuser, Basel, 1981, pp. 72–147.

[25] W. G a u t s c h i, On generating orthogonal polynomials, SIAM J. Sci. Statist. Comput. 3 (1982), 289-317.

[26] W. G a u t s c h i, Polynomials orthogonal with respect to the reciprocal gamma function, BIT 22 (1982), 387–389.

[27] W. G a u t s c h i,, How and how not to check Gaussian quadrature formulae, BIT 23 (1983), 209–216.

[28] W. G a u t s c h i, On some orthogonal polynomials of interest in theoretical chemistry, BIT 24 (1984), 473-483.

[29] W. G a u t s c h i, Computational aspects of orthogonal polynomials, In: Orthogonal Polynomials (Columbus, OH, 1989), NATO Adv. Sci. Inst. Ser. C: Math. Phys. Sci., 294 (P. Nevai, ed.), Kluwer, Dordrecht, 1990, 181–216.

[30] W. G a u t s c h i, Computational problems and applications of orthogonal polynomials, In: IMACS Annals on Computing and Applied Mathematics, Vol. 9: Orthogonal Polynomials and Their Applications (C. Brezinski, L. Gori and A. Ronveaux, eds.), IMACS, Baltzer, Basel, 1991, pp. 61–71.

[31] W. G a u t s c h i, Algorithm 726: ORTHPOL – a package of routines for generating orthogonal polynomials and Gauss-type quadrature rules, ACM Trans. Math. Software 20 (1994), 21–62.

[32] W. G a u t s c h i,, Orthogonal polynomials: applications and computation, Acta Numerica (1996), 45–119.

[33] W. G a u t s c h i, Moments in quadrature problems, Comput. Math. Appl. 33 (1997), 105–118.

[34] W. G a u t s c h i, Orthogonal Polynomials: Computation and Approximation, Clarendon Press, Oxford 2004.

[35] W. G a u t s c h i, H. L a n d a u, G. V. M i l o v a n o v i ć, Polynomials orthogonal on the semicircle, II, Constr. Approx. 3 (1987), 389–404.

[36] W. G a u t s c h i, G. V. M i l o v a n o v i ć, Gaussian quadrature involving Einstein and Fermi functions with an application to summation of series, Math. Comp. 44 (1985), 177–190.

[37] W. G a u t s c h i, G. V. M i l o v a n o v i ć, Polynomials orthogonal on the semicircle, J. Approx. Theory 46 (1986), 230-250.

[38] W. G a u t s c h i, G. V. M i l o v a n o v i ć, Spline approximations to spherically symmetric distributions, Numer. Math. 49 (1986), 111–121.

[39] W. G a u t s c h i, G. V. M i l o v a n o v i ć, s-orthogonality and construction of Gauss–Turán-type quadrature formulae, J. Comput. Appl. Math. 86 (1997), 205–218.

[40] A. G h i z z e t t i, A. O s s i c i n i, Quadrature Formulae, Akademie Verlag, Berlin, 1970.

[41] A. G h i z z e t t i, A. O s s i c i n i, Sull’ esistenza e unicità delle formule di quadratura gaussiane, Rend. Mat. (6) 8 (1975), 1–15.

[42] G. H. G o l u b, J. K a u t s k y, Calculation of Gauss quadratures with multiple free and fixed knots, Numer. Math. 41 (1983), 147–163.

[43] G. H. G o l u b, J. H. W e l s c h, Calculation of Gauss quadrature rules, Math. Comp. 23 (1969), 221–230.

[44] C. G. J. J a c o b i, Über Gaußs neue Methode, die Werte der Integrale näherungsweise zu finden, J. Reine Angew. Math. 30 (1826), 301–308.

[45] S. K a r l i n, A. P i n k u s, Gaussian quadrature formulae with multiple nodes, In: S. Karlin, C.A. Micchelli, A. Pinkus, I.J. Schoenberg (Eds.), Studies in Spline Functions and Approximation Theory, Academic Press, New York, 1976, pp. 113–141.

[46] S. K a r l i n, W. J. S t u d d e n, Tchebycheff Systems with Applications in Analysis and Statistics, Pure and Applied Mathematics, Vol. XV, John Wiley Interscience, New York, 1966.

[47] J. K a u t s k y, G. H. G o l u b, On the calculation of Jacobi matrices, Linear Algebra Appl. 52/53 (1983), 439–455.

[48] A. L. K u z m i n a, Interval quadrature formulae with multiple node intervals, Izv. Vuzov 7 (218) (1980), 39–44. (Russian)

[49] J. G. L a f r a m b o i s e, A. D. S t a u f f e r, Optimum discrete approximation of the Maxwell distribution, AIAA J. 7 (1969), 520–523.

[50] S.-Y. L e e, The inhomogeneous Airy functions Hi(z) and Gi(z), J. Chem. Phys. 72 (1980), 332–336.

[51] M. L é n á r d, Birkhoff quadrature formulae based on the zeros of Jacobi polynomials, Math. Comput. Modelling, 38 (2003), 917–927.

[52] J. M a, V. R o k h l i n, S. W a n d z u r a, Generalized Gaussian quadrature rules for systems of arbitrary functions, SIAM J. Numer. Anal., 33 (1996), 971–996.

[53] G. M a s t r o i a n n i, G. V. M i l o v a n o v i ć, Interpolation Processes – Basic Theory and Applications, Springer Verlag, Berlin – Heidelberg – New York, 2008.

[54] G. V. M i l o v a n o v i ć, Construction of s-orthogonal polynomials and Turán quadrature formulae, in: G.V. Milovanović (Ed.), Numerical Methods and Approximation Theory III (Niš, 1987), Univ. Niš, Niš, 1988, pp. 311–328.

[55] G. V. M i l o v a n o v i ć, Complex orthogonality on the semicircle with respect to Gegenbauer weight: theory and applications, In: Topics in Mathematical Analysis, A Volume Dedicated to the Memory of A. L. Cauchy (Th. M. Rassias, ed.), World Scientific, Singapore, 1989, pp. 695–722.

[56] G. V. M i l o v a n o v i ć, Summation of series and Gaussian quadratures, In: Approximation and Computation (R.V.M. Zahar, ed.), ISNM Vol. 119, Birkhäuser Verlag, Basel–Boston–Berlin, 1994, pp. 459–475.

[57] G. V. M i l o v a n o v i ć, Summation of series and Gaussian quadratures, II, Numer. Algorithms 10 (1995), 127-136.

[58] G. V. M i l o v a n o v i ć, Müntz orthogonal polynomials and their numerical evaluation, In: Applications and computation of orthogonal polynomials (W. Gautschi, G. H. Golub, G. Opfer, eds.), ISNM, Vol. 131, Birkhäuser, Basel, 1999, pp. 179–202.

[59] G. V. M i l o v a n o v i ć, Quadrature with multiple nodes, power orthogonality, and moment-preserving spline approximation, In: Numerical Analysis 2000, Vol. V, Quadrature and Orthogonal Polynomials (W. Gautschi, F. Marcellán, and L. Reichel, eds.), J. Comput. Appl. Math. 127 (2001), 267–286.

[60] G. V. M i l o v a n o v i ć, A. S. C v e t k o v i ć, Uniqueness and computation of Gaussian interval quadrature formula for Jacobi weight function, Numer. Math. 99 (2004), 141-162.

[61] G. V. M i l o v a n o v i ć, A. S. C v e t k o v i ć, Orthogonal polynomials and Gaussian quadrature rules related to oscillatory weight functions, J. Comput. Appl. Math. 179 (2005), 263-287.

[62] G. V. M i l o v a n o v i ć, A. S. C v e t k o v i ć, Orthogonal polynomials related to the oscillatory-Chebyshev weight function, Bull. Cl. Sci. Math. Nat. Sci. Math. 30 (2005), 47-60.

[63] G. V. M i l o v a n o v i ć, A. S. C v e t k o v i ć, Gauss-Laguerre interval quadrature rule, J. Comput. Appl. Math. 182 (2005), 433-446.

[64] G. V. M i l o v a n o v i ć, A. S. C v e t k o v i ć, Gaussian type quadrature rules for Müntz systems, SIAM J. Sci. Comput. 27 (2005), 893-913.

[65] G. V. M i l o v a n o v i ć, A. S. C v e t k o v i ć, Gauss-Radau and Gauss-Lobatto interval quadrature rules for Jacobi weight function, Numer. Math. 102 (2006), 523-542.

[66] G. V. M i l o v a n o v i ć, A. S. C v e t k o v i ć, Gauss-Hermite interval quadrature rule, Comput. Math. Appl. 54 (2007), 544-555.

[67] G. V. M i l o v a n o v i ć, A. S. C v e t k o v i ć, Nonstandard Gaussian quadrature formulae based on operator value, (to appear).

[68] G. V. M i l o v a n o v i ć, M. M. S p a l e v i ć, Error bounds for Gauss-Turn quadrature formulae of analytic functions, Math. Comp. 72 (2003), 1855–1872.

[69] G. V. M i l o v a n o v i ć, M. M. S p a l e v i ć, Error analysis in some Gauss-Turán-Radau and Gauss-Turán-Lobatto quadratures for analytic functions, J. Comput. Appl. Math. 165 (2004), 569–586.

[70] G. V. M i l o v a n o v i ć, M. M. S p a l e v i ć, An error expansion for some Gauss-Turán quadratures and L^{1}-estimates of the remainder term, BIT 45 (2005), 117–136.

[71] G. V. M i l o v a n o v i ć, M. M. S p a l e v i ć, Bounds of the error of Gauss-Turán-type quadratures, J. Comput. Appl. Math. 178 (2005), 333–346.

[72] G. V. M i l o v a n o v i ć, M. M. S p a l e v i ć, Gauss-Turán quadratures of Kronrod type for generalized Chebyshev weight functions, Calcolo 43 (2006), 171–195.

[73] G. V. M i l o v a n o v i ć, M. M. S p a l e v i ć, A note on the bounds of the error of Gauss-Turán-type quadratures, J. Comput. Appl. Math. 200 (2007), 276–282.

[74] G. V. M i l o v a n o v i ć, M. M. S p a l e v i ć, A. S. C v e t k o v i ć, Calculation of Gaussian type quadratures with multiple nodes, Math. Comput. Modelling 39 (2004), 325–347.

[75] G. V. M i l o v a n o v i ć, A. K. V a r m a, On Birkhoff (0, 3) and (0, 4) quadrature formulae, Numer. Funct. Anal. Optim. 18 (1997), 427–433.

[76] V. P. M o t o r n y i, On the best quadrature formulae in the class of functions with bounded r-th derivative, East J. Approx. 4 (1998), 459–478.

[77] M. O m l a d i č, Average quadrature formulas of Gauss type, IMA J. Numer. Anal. 12 (1992), 189–199.

[78] M. O m l a d i č, S. P a h o r, A. S u h a d o l c, On a new type quadrature formulas, Numer. Math. 25 (1976), 421–426.

[79] R. P i e s s e n s, M. B r a n d e r s, Tables of Gaussian quadrature formulas, Appl. Math. Progr. Div. University of Leuven, Leuven, 1975.

[80] Fr. P i t n a u e r, M. R e i m e r, Interpolation mit Intervallfunktionalen, Math. Z. 146 (1976), 7–15.

[81] T. P o p o v i c i u, Sur une généralisation de la formule d’intégration numérique de Gauss, Acad. R. P. Romîne Fil. Iaşi Stud. Cerc. Şti. 6 (1955), 29–57. (Romanian)

[82] R. N. S h a r i p o v, Best interval quadrature formulas for Lipschitz classes, In: Constructive theory of functions and functional analysis, No. IV, pp. 124–132, Kazan. Gos. Univ., Kazan’, 1983.

[83] Y. G. S h i, Theory of Birkhoff Interpolation, Nova Science Publishers, Inc., Hauppauge, NY, 2003.

[84] Y. G. S h i, Power Orthogonal Polynomials, Nova Science Publishers, Inc., New York, 2006.

[85] D. D. S t a n c u, On a class of orthogonal polynomials and on some general quadrature formulas with minimum number of terms, Bull. Math. Soc. Sci. Math. Phys. R. P. Romîne (N.S) 1 (49) (1957), 479–498.

[86] D. D. S t a n c u, On the interpolation formula of Hermite and some applications of it, Acad. R. P. Romîne Fil. Cluj Stud. Cerc. Mat. 8 (1957) 339–355 (Romanian).

[87] D. D. S t a n c u, Generalization of the quadrature formula of Gauss-Christoffel, Acad. R.P. Romîne Fil. Iaşi Stud. Cerc. Şti. Mat. 8 (1957), 1–18 (Romanian).

[88] D. D. S t a n c u, On a class of orthogonal polynomials and on some general quadrature formulas with minimum number of terms, Bull. Math. Soc. Sci. Math. Phys. R. P. Romîne (N.S) 1 (49) (1957), 479–498.

[89] D. D. S t a n c u, A method for constructing quadrature formulas of higher degree of exactness, Com. Acad. R. P. Romîne 8 (1958), 349–358 (Romanian).

[90] D. D. S t a n c u, On certain general numerical integration formulas, Acad. R. P. Romîne. Stud. Cerc. Mat. 9(1958), 209–216 (Romanian).

[91] D. D. S t a n c u, On certain general numerical integration formulas, Acad. R. P. Romîne. Stud. Cerc. Mat. 9 (1958), 209–216 (Romanian).

[92] D. D. S t a n c u, Sur quelques formules générales de quadrature du type Gauss-Christoffel, Mathematica (Cluj) 1 (24) (1959), 167–182.

[93] D. D. S t a n c u, Sur quelques formules générales de quadrature du type Gauss-Christoffel, Mathematica (Cluj) 1 (24) (1959), 167–182.

[94] D. D. S t a n c u, An extremal problem in the theory of numerical quadratures with multiple nodes, In: Proceedings of the Third Colloquium on Operations Research (Cluj-Napoca, 1978), Univ. “Babeş-Bolyai”, Cluj-Napoca, 1979, pp. 257–262.

[95] D. D. S t a n c u, A. H. S t r o u d, Quadrature formulas with simple Gaussian nodes and multiple fixed nodes, Math. Comp. 17 (1963), 384–394.

[96] N. M. S t e e n, G. D. B y r n e, E. M. G e l b a r d, Gaussian quadratures for the integrals ∫_{0}^{∞} exp(−x^{2})f(x) dx and ∫_{0}^{b} exp(−x^{2})f(x) dx, Math. Comp. 23 (1969), 661–671.

[97] T. J. S t i e l t j e s, Sur une généralisation de la théorie des quadratures mécaniques, C.R. Acad. Sci. Paris, 99 (1884), 850–851.

[98] A. H. S t r o u d, D. D. S t a n c u, Quadrature formulas with multiple Gaussian nodes, J. SIAM Numer. Anal. Ser. B 2 (1965), 129–143.

[99] P. T u r á n, On the theory of the mechanical quadrature, Acta Sci. Math. Szeged 12 (1950), 30–37.

[100] P. T u r á n, On some open problems of approximation theory, J. Approx. Theory 29 (1980), 23–85.

[101] A. K. V a r m a, On Birkhoff quadrature formulas, Proc. Amer. Math. Soc. 97 (1986), 38–40.

Department Mathematics Faculty of Electronic Engineering University of Niš P. O. Box 73 18000 Niš Serbia