High degree quadrature rules with pseudorandom rational nodes

TL;DR

This paper introduces a method to construct high-degree quadrature rules with pseudorandom rational nodes by combining companion rules and mean rules, optimizing approximation accuracy.

Contribution

It presents a novel framework for creating high-degree quadrature rules using linear combinations of lower-degree rules with rational pseudorandom nodes.

Findings

Construction of higher-degree rules from lower-degree ones

Mean rule as the best least-squares approximation

Explicit examples of quadrature rule construction

Abstract

After introducing the definitions of positive, negative and companion rules, from a given pair of companion rules we construct a new rule with higher degree of precision The scheme is generalized giving rise to a transformation which we call the mean rule. We show that the mean rule is the best approximation, in the sense of least-squares, obtained from a linear combination of two rules of the same degree of precision. Finally, we show that a rule of degree 2k+1 can be constructed as linear combination of k+1 rules of degree one and rational pseudorandom nodes. Several worked examples are presented.

Tables

Table 1. Table 1: Rule ( 13 ) for t 1 = 0 subscript 𝑡 1 0 t_{1}=0 and t 2 = − t 0 subscript 𝑡 2 subscript 𝑡 0 t_{2}=-t_{0} .

Table 2. Table 2: Errors of A ​ ( g ) 𝐴 𝑔 A(g) , B ​ ( g ) 𝐵 𝑔 B(g) and the mean rule of degree 7.

Equations

$$\sigma (t)=a+\frac{b-a}{2}(t+1),-1\le t\le 1$$

$${\mu }_{m+1}={\int }_{-1}^1{\varphi }_{m+1}(t)dt.$$

$$\left\{\begin{array}[]{ll}Q(\phi_{0})&=\mu_{0}\\ Q(\phi_{1})&=\mu_{1}\\ &\vdots\\ Q(\phi_{m})&=\mu_{m},\end{array}\right.$$

$$Q({\varphi }_{m+1})\ne {\mu }_{m+1}.$$

$$\mu ={\int }_{-1}^1{\varphi }_{m+1}(t)dt=\frac{2}{m+2}.$$

$${\gamma }_Q=\mu -Q({\varphi }_{m+1}).$$

$${\mu }_A=A({\varphi }_{m+1})\ne {\mu }_B=B({\varphi }_{m+1}),$$

$$Y(g)=\frac{\mu -{\mu }_B}{{\mu }_A-{\mu }_B}A(g)+\frac{{\mu }_A-\mu }{{\mu }_A-{\mu }_B}B(g),$$

$${\mu }_B<\mu <{\mu }_A\text{mbox}or{\mu }_A<\mu <{\mu }_B,$$

$$A(g)\le Y(g)\le B(g)\text{mbox}orB(g)\le Y(g)\le A(g).$$

$$\begin{array}[]{ll}Y(\phi_{j})&=\displaystyle\frac{\mu-\mu_{B}}{\mu_{A}-\mu_{B}}\,A(\phi_{j})+\displaystyle\frac{\mu_{A}-\mu}{\mu_{A}-\mu_{B}}\,B(\phi_{j}),\\ \\ &=\displaystyle\frac{\mu-\mu_{B}+\mu_{A}-\mu}{\mu_{A}-\mu_{B}}\,\,I(\phi_{j})=I(\phi_{j})\ .\end{array}$$

$$\begin{array}[]{ll}Y(\phi_{m+1})&=\displaystyle\frac{\mu-\mu_{B}}{\mu_{A}-\mu_{B}}\,\mu_{A}+\displaystyle\frac{\mu_{A}-\mu}{\mu_{A}-\mu_{B}}\,\mu_{B},\\ \\ &=\displaystyle\frac{\mu\mu_{A}-\mu\,\mu_{B}}{\mu_{A}-\mu_{B}}=\mu=I(\phi_{m+1})\ .\end{array}$$

$${\mu }_A<\mu <{\mu }_B.$$

$$Y(g)=\frac{\alpha A(g)+\beta B(g)}{\alpha +\beta },$$

$$c_1=\frac{\alpha }{\alpha +\beta }>0\text{mbox}and{c}_2=\frac{\beta }{\alpha +\beta }>0,$$

$$M(g)=2g(0),T(g)=g(-1)+g(1).$$

$$\mu ={\int }_{-1}^1{\varphi }_2(t)dt=2/3.$$

$${\gamma }_M=\mu -M({\varphi }_2)=2/3>0\text{mbox}and{\gamma }_T=\mu -T({\varphi }_2)=-1/3<0.$$

$$\begin{array}[]{ll}Y(g)&=\displaystyle{\frac{\mu-\mu_{T}}{\mu_{M}-\mu_{T}}M(g)+\frac{\mu_{M}-\mu}{\mu_{M}-\mu_{T}}T(g)}=\displaystyle\frac{2}{3}\,M(g)+\displaystyle\frac{1}{3}\,T(g).\end{array}$$

$$Y(g)=\frac{1}{3}\left[g(-1)+4g(0)+g(1)\right],$$

$$A(g)=g\left(-\frac{\sqrt{3}}{3}\right)+g\left(\frac{\sqrt{3}}{3}\right),$$

$$S(g)=\frac{1}{3}\left(g(-1)+4g(0)+g(1)\right),$$

$$\begin{array}[]{l}\mu_{A}=A(t^{4})=2/9\quad\Longrightarrow\quad\gamma_{A}=\mu-A(\phi_{4})=8/45>0\\ \mu_{S}=S(t^{4})=2/3\quad\Longrightarrow\quad\gamma_{S}=\mu-S(\phi_{4})=-4/15<0,\\ \end{array}$$

$$Y(g)=\frac{{\gamma }_S}{{\mu }_A-{\mu }_S}A(g)-\frac{{\gamma }_A}{{\mu }_A-{\mu }_S}S(g)=\frac{3A(g)+2S(g)}{5}.$$

$$Y(g)=\frac{1}{15}\left[2g(-1)+9g(-\sqrt{3}/3)+8g(0)+9g(\sqrt{3}/3)+2g(1)\right].$$

$$\begin{array}[]{l}Y(\phi_{4})=2/5\quad\mbox{and}\quad I(\phi_{4})=2/5,\\ Y(\phi_{5})=0\quad\mbox{and}\quad I(\phi_{5})=0,\\ Y(\phi_{6})=14/45\quad\neq\quad I(\phi_{6})=2/7\ .\\ \end{array}$$

$$I(g)={\int }_{-1}^1\frac{2}{1+t^2}dt=\pi .$$

$$E_A(g)=I(g)-A(g)\simeq 0.14,E_S(g)=I(g)-S(g)\simeq -0.19,$$

$$Q(g)=A_{0}g(t_0)+A_{1}g(t_1),$$

$${\mathcal{L}}_2:1/3+t_0{t}_1=0.$$

Full text

Key words: Positive rule; Negative rule; Companion rules; Combined rule; Midpoint; Trapezoidal; Simpson rule; Mean rule; Pseudorandom node.

2010 Mathematics Subject Classification: 65-05, 65D30, 65D32.

1 Introduction

Given two quadrature rules $A(g)$ and $B(g)$ with the same degree of precision $m\geq 0$, we begin by proposing a scheme to construct a new rule $Y(g)$ of higher degree. One desirable assumption is that the rules $A(g)$ and $B(g)$ are companion, in the sense that one can assign opposite signals to the respective error. In particular, we show that the basic quadrature rules known as midpoint, trapezoidal, and Simpson rules can all be obtained as linear combinations of companion rules of lower degree of precision. The theoretical background applied in this work relies on the method of undetermined coefficients ([2], p. 565), ([4] p. 170).

We generalize the referred scheme by presenting a rule transformation $W(g)$ (defined in the set ${\cal Q}$ of rules of degree $m$), which to a pair of rules $\left(A(g),B(g)\right)$ of ${\cal Q}$, assigns a new rule of greater degree. This leads to an algorithm to obtain quadrature rules of arbitrary order of precision, as suggested in the worked examples. The rule $W(g)$ is called the mean rule since it is a weighted mean of $A(g)$ and $B(g)$. It can be seen as a least-squares best approximation as discussed in paragraph 4.1.

The main results of this paper are discussed in Section 5. We first show that if one takes a set of $k+1$ of open rules of degree one, $Q_{0}(g),\ldots,Q_{k}(g)$, whose first member is the midpoint rule $Q_{0}(g)=2\,g(0)$ and the other members have two symmetrical rational nodes, there exists a unique linear combination of the rules such that the combined rule $W_{k}(g)$ has degree $2\,k+1$ (see Proposition 3). As an illustration we apply the composite version of the rule $W_{5}(g)$ to obtain approximations of $\pi$ with 60 significant digits (Example 8) using the model function $g(t)=2/(1+t^{2}),\,\-1\leq t\leq 1$.

Finally we show that one can expand the scheme considering combined rules where the nodes of the $1$-degree starting rules are pseudorandom rational numbers (paragraph 5.1). In particular, we use a pseudorandom $151$-degree combined rule giving an approximation of $\pi$ with more than $500$ significant digits as detailed in Example 9.

2 Notation and definitions

A quadrature rule $Q(f)$ is an approximation of the integral $\displaystyle\int_{a}^{b}f(x)dx$ obtained using values of $f$ (and/or its derivatives) on a discrete set of points in $[a,b]$ (see for instance [3]).

Without loss of generality, we consider $[-1,1]$ to be the interval of integration and we will denote by $Q(g)$ a general quadrature rule to approximate the integral $I(g)=\displaystyle\int_{-1}^{1}g(t)\,dt$.

Note that a rule $Q(g)$, defined in $[-1,1]$, can be rewritten for the interval $[a,b]$ as $Q(f)$, through a change of variable defined by the bijection

[TABLE]

where $g(t)=f(\sigma(t))\ .$

The monomials $1,t,t^{2},\ldots$, are denoted by $\phi_{j}(t)=t^{j}$, for $j=0,1,\ldots$.

Definition 2.1**.**

(Degree $m$ of rule)

Let $m\geq 0$ be an integer, $\mu_{0}=2$, $\mu_{j}=0$ for odd $j$ and $\mu_{j}=\displaystyle\int_{-1}^{1}\phi_{j}(t)\,dt=2/(j+1)$, for even $j$, and let

[TABLE]

A rule $Q(g)$ is of degree $m\geq 0$ if it is exact for $\phi_{j}(t)$, with $0\leq j\leq m$, but not exact for $\phi_{m+1}(t)$, that is,

[TABLE]

but

[TABLE]

In what follows, $\mu_{m+1}$ (a.k.a. the principal moment) plays a key role and we shorten its notation to

[TABLE]

Definition 2.2**.**

(Sign of a rule)

Let $Q_{m}(f)$ be a quadrature rule of degree $m\,(m\geq 0)$, such that

[TABLE]

The rule is positive (resp. negative) if $\gamma_{Q}>0$ (resp. $\gamma_{Q}<0$).

A pair of rules of the same degree whose respective value $\gamma$ have opposite signs are rules of particular interest. Such rules will be called companion rules.

Definition 2.3**.**

(Companion rules)

Two rules $A(g)$ and $B(g)$, of the same degree $m$, such that $\gamma_{A}$ and $\gamma_{B}$ have opposite signs are called companion rules.

3 Linear combination of two rules

We now address the problem of combining a pair of companion rules and show that the resulting rule is not only a weighted mean of the two rules considered but also it has higher degree of precision. An obvious advantage of this approach is that at a minor computational cost, the value of the new rule can be closer to the integral than the values of the rules of the starting pair (see Example 2).

Proposition 1**.**

Let $A(g)$ and $B(g)$ be two rules both of degree $m\geq 0$, such that

[TABLE]

and consider the linear combination of the rules

[TABLE]

where $\mu$ is the principal moment in (1). Then, the degree of the rule $Y(g)$ is at least $m+1$.

Corollary 1**.**

If $A(g)$ and $B(g)$ are companion rules of degree $m\geq 0$, such that

[TABLE]

then the combined rule (3) has degree at least $m+1$ and

[TABLE]

Proof of Proposition 1.

Both rules are exact for $\phi_{j}(t)=t^{j}$, with $j=0,1,\ldots,m$, that is, $A(\phi_{j})=B(\phi_{j})=I(\phi_{j})$. Therefore,

[TABLE]

So, the rule $Y(g)$ has degree at least $m$. It remains to show that this rule is exact for $\phi_{m+1}(t)$ which implies that its degree is at least $m+1$. From (2), we have

[TABLE]

∎

Proof of Corollary 1.

We assume that $A(g)$ is positive and $B(g)$ is negative being the proof in the other case completely analogous. By definition of sign of a rule, we have $\mu>\mu_{A}$ and $\mu<\mu_{B}$, i.e.

[TABLE]

Let $\alpha=\mu-\mu_{B}<0$ and $\beta=\mu_{A}-\mu<0$ which implies $\alpha+\beta=\mu_{A}-\mu_{B}<0$. The rule $Y(g)$ can be written as

[TABLE]

where

[TABLE]

So, the rule $Y(g)$ is a linear combination of the rules $A(g)$ and $B(g)$ with positive coefficients $c_{1}$ and $c_{2}$. Consequently, the inequalities in (4) hold, and the rule $Y(g)$ is a weighted mean of the companion rules $A(g)$ and $B(g)$. ∎

Example 1**.**

($S=2/3\,M+1/3\,T$)

The well-konown midpoint and trapezoidal rules are, respectively,

[TABLE]

For $\phi_{2}(t)=t^{2}$, the principal moment is

[TABLE]

Since $M(\phi_{2})=0$ and $T(\phi_{2})=1$, we have

[TABLE]

Thus, by Definition 2.2, $M(g)$ is positive, while $T(g)$ is negative. Moreover, as $\gamma_{M}\ .\gamma_{T}<0$, the rules $M$ and $T$ are companion rules (see Definition 2.3).

Let us show that the linear combination (3) of these two rules of degree 1 coincides with the Simpson rule, which is a rule of degree $m=3$.

As $\mu_{M}=M(\phi_{2})=0$ and $\mu_{T}=T(\phi_{2})=2$, from (3) it follows

[TABLE]

Taking into account the expressions (5), we get

[TABLE]

which coincides with Simpson rule $S(g)$ (see next example).

Example 2**.**

(A combined rule of degree 5 using the Simpson rule )

One easily verifies that in $[-1,1]$ the (open) rule

[TABLE]

and the Simpson rule (closed)

[TABLE]

are both of degree $3$. Let us show that they are companion rules, and verify that the respective combined rule (3) has degree $m=5$.

The principal moment is $\mu=\displaystyle\int_{-1}^{1}t^{4}\,dt=2/5$, and

[TABLE]

and so $A(f)$ is positive while $S(f)$ is negative. We have $\mu_{A}-\mu_{S}=2/9-2/3=-4/9$ and the rule $Y(g)$ has the form

[TABLE]

The explicit expression of the weighted mean (6) is

[TABLE]

The rule $Y(g)$ has degree $m=5$ since

[TABLE]

The result in (7) shows the dependence of $Y(g)$ on 5 nodes. In paragraph 3.2 we will construct another rule of degree 5 using only 3 nodes.

In computational terms, given two rules $A(g)$ and $B(g)$ one does not need to use the expression of the rule $Y(g)$ in terms of the nodes like in (7) but just the linear combination (6).

In order to observe the numerical improvement one can get passing from a pair of companion rules to the respective combined rule $Y(g)$, let us approximate the following integral which will be used as a test model in the subsequent examples:

[TABLE]

For the first rule we obtain $A(g)=3$, for the Simpson rule $S(g)=10/3\simeq 3.333\ldots$ and for the combined rule $Y(g)=\left(3\,A(g)+2\,S(g)\right)/5=47/15\simeq 3.1333\ldots$. Thus, the respective errors are

[TABLE]

while $E_{Y}(g)=I(g)-Y(g)\simeq 0.008.$ Once computed $A(g)$ and $S(g)$ the cost for computing $Y(g)$ is only an addition, two multiplications and a division. In this example an $\simeq 5\verb+$ relative error in $A(g)$ and $S(g)$ yields to a relative error of $\simeq 0.3\verb+$ of the rule $Y(g)$, meaning that the combined rule is approximately 15 times more accurate than the two referred companion rules.

3.1 Families of companion two-point rules

Let $t_{0},t_{1}$ be two distinct nodes belonging to $[-1,1]$, and the two-point quadrature rule

[TABLE]

where the parameters $A_{0}$, $A_{1}$ will be determined in order that the rule has a degree of precision at least 1.

Facts:

1. (i)

There are infinite choices of points $(t_{0},t_{1})$ in the square region $D=[-1,1]\times[-1,1]$, for which the rule $Q(g)$ (considered as a function of $(t_{0},t_{1})$) is either positive or negative and simultaneously has degree 1; 2. (ii)

There are more positive rules than negative ones; 3. (iii)

The points $(t_{0},t_{1})\in D$ for which $Q(g)$ has exact degree $m=2$, lie on the hyperbole ${\cal L}_{2}$,

[TABLE] 4. (iv)

There is exactly one (two-point) rule $Q(g)$ with degree $m=3$, namely for $t_{0}=-\sqrt{3}/3$, $t_{1}=\sqrt{3}/3$, that is given by

[TABLE]

In Figure 1 the points in the square region $D$ for which the rule $Q(g)$ is either positive or negative are displayed in gray and white respectively. The boundary of this region contains points of the hyperbole ${\cal L}_{2}$ in (9). For $(t_{0},t_{1})$ belonging to the hyperbole ${\cal L}_{2}$ the rule $Q(g)$ has exactly degree one, meaning that there is an infinite number of rules of $1$-degree each one with nodes $t_{0},t_{1}$ for which the point $(t_{0},t_{1})$ belongs to one of the two branches of the hyperbole ${\cal L}_{2}$ in Figure 1. Moreover, one immediately sees that the two-point closed Newton-Cotes rule ($t_{0}=-1$ and $t_{1}=1$) is a negative rule of order one, as it was assumed in Example 1. It is also worth to recall that any closed Newton-Cotes rule of any degree $m>1$ is negative as well, whereas the open Newton-Cotes rules are all positive.

The mathematical aspects behind the facts $(i)$-$(iv)$ above and the geometry of figures 1-2 can be explained as follows.

Instead of the canonical polynomial basis of monomials of degree $\leq 3$, we consider the basis

[TABLE]

Imposing the condition that the rule $Q(f)$ is of degree at least 1, the parameters $A_{0},A_{1}$ in (8) are computed. That is, considering $Q(\Psi_{0})=I(\Psi_{0})$ and $Q(\Psi_{1})=I(\Psi_{1})$, we obtain the following triangular system in the unknowns $A_{0},A_{1}$:

[TABLE]

The solution of this system is $A_{1}=-2\,t_{0}/(t_{1}-t_{0})$ and $A_{0}=2\,t_{1}/(t_{1}-t_{0})$. Consequently, the two-point rule

[TABLE]

has degree at least one.

The rule applied to the nodal polynomial $\Psi_{2}(t)$ gives $Q(\Psi_{2})=0$, whereas the moment $I(\Psi_{2})=2/3+2\,t_{0}\,t_{1}$. So,

[TABLE]

Thus, the rule is of degree one and positive for the points $(t_{0},t_{1})$ displayed in gray color in Figure 1 and negative for the white points.

Remark 1**.**

One advantage of the polynomial basis $\Psi_{0}(t),\Psi_{1}(t),\ldots$ is that from the nodal polynomial onwards the rule is null and so the value of the parameter $\gamma_{Q}$ coincides with $\mu=I(\Psi_{3})$. This is the reason why we call $\mu$ the principal moment of the rule which also gives the sign of a rule in the sense of Definition 2.2.

The rule $Q(g)$ has degree $\geq 3$ if and only if (9) holds and it has null moment $\mu=\displaystyle\int_{-1}^{1}\Psi_{3}(t)\,dt$. The rule has degree 2 and is positive when $\mu>0$ and negative otherwise. The condition $\mu=0$ is given by the equation

[TABLE]

In Figure 2, the points for which $\mu>0$ (positive rule of degree 2) are displayed in dark gray and in light gray the points with $\mu<0$ (negative rule). The boundary of the respective domains is shown in dark bold. This boundary is the algebraic curve ${\cal L}_{3}$ defined by a cubic polynomial as in (11).

The intersection points of the curves ${\cal L}_{2}$ and ${\cal L}_{3}$ are symmetrically distributed with respect to the origin and are given by $t_{0}=\pm\sqrt{3}/3$, with $t_{1}=-t_{0}$. In figures 1-2 these intersection points are the center of the circles in blue.

This fully justify the facts $(i)$-$(iv)$. Analogous procedures enable us to obtain the geometry of a family of rules with three distinct nodes. Indeed, following the same lines as above, we will obtain a famous (positive) rule of degree 5, known as Gauss-Legendre rule (see for instance [1], Ch. 6) with three nodes, which might be be combined with the 5 degree rule $Y(g)$ given in Example 2.

3.2 Families of companion three-point rules

In $D=\left\{(t_{0},t_{2},t_{3}):\,-1\leq t_{0}\leq 1,\,-1\leq t_{1}\leq 1,\,-1\leq t_{2}\leq 1\right\}$, we consider the set of 3-point quadrature rules

[TABLE]

Firstly, we obtain the weights $A_{0},A_{1},A_{2}$ in order that all the rules in the set have degree 2 and are either positive or negative. This will enable us to find candidates to pairs of companion rules of degree 2.

The weights $A_{i}$ can be written as functions $A_{i}=\Omega_{i}:D\subset\mathbb{R}^{3}\mapsto\mathbb{R}$, for $i=0,1,2$. In fact, taking the polynomials

[TABLE]

the rule has degree $\geq 2$ if and only if it is exact for $\Psi_{0},\Psi_{1}$ and $\Psi_{2}$. That is, the weights are solution to the triangular system

[TABLE]

Since we are assuming $t_{0}\neq t_{1}\neq t_{2}$, this system has a unique solution,

[TABLE]

Thus,

[TABLE]

where $\Omega_{i}$ are given by (12). After some simplifications the expressions in (12) become

[TABLE]

Now, let us consider the polynomials

[TABLE]

Remark 2**.**

*Note that the rule (13) is null when applied to $\Psi_{j}$, for $j\geq 3$, and so the degree of the rule is the same of the first index $j\geq 3$ for which $I(\Psi_{j})\neq 0$. *

By construction, the rule $Q(g)$ has order $\geq 2$. Therefore, for any polynomial of the vector space ${\cal P}_{2}$ (polynomials of degree $\leq 2$) the rule is exact. In particular, it is exact for the elements of the canonical basis of ${\cal P}_{2}$, that is for $\phi_{j}(t)=t^{j}$, with $0\leq j\leq 2$.

Since $I(\phi_{3})=\displaystyle\int_{-1}^{1}t^{3}\,dt=0$, we have $\gamma_{3}=-Q(\phi_{3})$, where $\gamma_{j}=I(\phi_{j})-Q(\phi_{j}).$ Thus, $Q(g)$ has order exactly 2 if $Q(\phi_{3})\neq 0$ and order $\geq 3$ if $Q(\phi_{3})=0$.

Let $P$ be the following polynomial

[TABLE]

where $\Omega_{i}$ are given by (14).

The rule (13) is:

1. (i)

positive if $P(t_{0},t_{1},t_{2})<0$; 2. (ii)

negative if $P(t_{0},t_{1},t_{2})>0$; 3. (iii)

of degree $\geq 3$ if $(t_{0},t_{1},t_{2})$ is a root of the polynomial equation $P(t_{0},t_{1},t_{2})=0$.

Specific instances of the rule (13) with weights (14) are given in the following example.

Example 3**.**

Let $t_{0}=-15/16,t_{1}=-7/8$, $t_{2}=-3/4$ (see Figure 3 where the point $A$ has coordinates111After defining the polynomial expression (15), the coordinates of the point A have been found using the predicate $\{P(t_{0},t_{1},t_{2})<0,-1\leq t_{0}\leq 1,-1\leq t_{1}\leq 1,-1\leq t_{2}\leq 1,t_{0}\neq t_{1}\neq t_{2}\}$ as argument to the Mathematica [6] command FindInstance. $(t_{0},t_{1},t_{2})$) and $P$ as in (15). Computing the $\Omega_{i}$’s, given by (14) we get $P(t_{0},t_{1},t_{2})<0$. Therefore the rule

[TABLE]

is positive, with degree 2.

Considering now $t_{0}=-3/4,t_{1}=-7/8$ and $t_{2}=-3/4$ (point $B$ in Figure 3) we have $P(t_{0},t_{1},t_{2})>0$. Thus, the rule

[TABLE]

is negative, with degree 2. Taking $\phi_{3}(t)=t^{3}$, the coefficients (3) of the combined rule are:

[TABLE]

As $A(g)$ and $B(g)$ have errors of opposite sign they are companion rules (see Definition 2.3). The combined rule is

[TABLE]

which coincides with the (open) 4-point rule

[TABLE]

It can be verified that $Y(\phi_{j})=Y(\phi_{j})$, for $j=0,1,2,3$, and

[TABLE]

Thus the combined rule $Y(g)$ is negative and of degree $m=3$.

As an exercise, the interested reader can verify that for any of the $6$ points referred in Figure 3 the corresponding rule $Q(g)$ in (13) has degree 3. Thus, one may conclude that there exists an infinite set of rules of the referred type which are negative and of degree 3. Consequently, as all rules belonging to the set of open Newton-Cotes rules, ${\cal N}{\cal C}$, are positive, so they are good candidates to use in pairs of companion rules in order to obtain combined rules of arbitrary order.

Example 4**.**

($3$-point Gauss-Legendre rule of degree 5)

A simple choice of nodes $t_{0},t_{1}$ and $t_{2}$ for the weights (14) of the rule (13) is $t_{1}=0$, and $t_{2}=-t_{0}$, giving

[TABLE]

For the polynomials $\phi_{j}(t)=t^{j}$, with $j\geq 3$, Table 1 displays the errors $\gamma_{j}=I(\phi_{j})-Q(\phi_{j})$. The rule has degree $\geq 4$ if and only if $\gamma_{4}=0$, that is, for $t_{0}=\sqrt{3/5}$. In this case the rule in (16) becomes

[TABLE]

Taking into account the values of $\gamma_{5}=0$ and $\gamma_{6}>0$, the rule (17) has degree 5 and is positive. This is the Gauss-Legendre rule with 3 nodes. In Figure 3, the points $L_{1}$ and $L_{2}$ correspond to the nodes $(t_{0},t_{1},t_{2})=\pm(-\sqrt{3/5},0,\sqrt{3/5})$ for which the rule has degree 5.

4 The mean transformation rule $W(g)$

In what follows, ${\cal Q}$ denotes the set of rules with degree $m\geq 0$. In Proposition 1, a new rule has been assigned to a pair of companion rules belonging to ${\cal Q}$. In order to generalize this scheme, let us define a transformation

[TABLE]

where $W(g)$ will enjoy analogous properties of the linear combination $Y(g)$ in (3), in the sense that the rule $W(g)$ has degree greater than those of the arguments $A(g)$ and $B(g)$. The rule $W(g)$ will be called mean rule.

As before, for a given rule $Q(g)\in{\cal Q}$, we compute the quantities $\mu$ and $\mu_{Q}$, defined by

[TABLE]

where $\phi_{j}(t)=t^{j},\,\,j=0,1,\ldots$.

Definition 4.1**.**

(Mean rule)

Let ${\cal Q}$ be the set of rules of degree $m\geq 0$ and $A(g),B(g)$ belonging to ${\cal Q}$. The mean rule of $A(g)$ and $B(g)$ is

[TABLE]

where $\mu_{A}$ and $\mu_{B}$ are as in (18).

Proposition 2**.**

The mean rule (19) has degree at least $m+1$.

Proof.

As the rules $A(g)$ and $B(g)$ have degree $m\geq 0$, we know that $A(\phi_{j})=B(\phi_{j})=0$, for odd $j$, and

[TABLE]

In the case $\mu_{A}=\mu_{B}$, from (20), it follows

[TABLE]

So, $W(g)$ has degree $\geq m+1$.

When $\mu_{B}-\mu_{A}\neq 0$, let us show that there exists a unique pair $(\alpha,\beta)$ such that

[TABLE]

has degree $\geq m+1$.

For $j=0,1,\ldots,m,m+1$, substituting $g$ in $W(g)$ by each $\phi_{j}(t)$ and applying (20), we obtain the linear system

[TABLE]

Since $\mu_{A}\neq\mu_{B}$ this system has the unique solution

[TABLE]

Thus, $W(g)$ has degree at least $m+1$. ∎

Example 5**.**

We now construct the mean rule of (10) and the Simpson’s rule (both rules of degree 3) and show that the mean rule has degree 5. The starting rules will be denoted by $A(g)$ and $S(g)$, respectively:

[TABLE]

We have,

[TABLE]

Thus,

[TABLE]

As $W(\phi_{5})=I(\phi_{5})=0$ the rule has degree $m=5$.

We note that $W(g)$ is a companion rule of the positive Gauss-Legendre rule (17) since $\gamma_{6}=I(\phi_{6})-W(\phi_{6})=-8/315<0$. We may also obtain the mean rule of (17) and (21). That is,

[TABLE]

This rule has degree $m=7$ and is negative since $\gamma_{8}=I(\phi_{8})-\widetilde{W}(\phi_{8})=-16/1575$.

Let $g(t)=2/(1+t^{2})$ and $I(g)=\int_{-1}^{1}g(t)\,dt=\pi$. We have

[TABLE]

Recall that the error with Simpson’s rule is approximately $-0.19$ (see Example 2) and so the mean rule $\widetilde{W}$ leads to a remarkable gain in accuracy.

4.1 The mean rule $W(g)$ as a least-squares approximation

We now show that given two rules $A(g)$ and $B(g)$, its mean rule $W(g)$ is the least-squares approximation to the vector of moments

[TABLE]

where, as before, the moments are: $\mu_{j}=\displaystyle\int_{-1}^{1}\phi_{j}(t)\,dt$ (for $j=0,\ldots,m$) and $\mu=\displaystyle\int_{-1}^{1}\phi_{m+1}(t)\,dt$.

We know that $A(\phi_{j})=B(\phi_{j})=\mu_{i}$, for $i=0,\ldots m$ and $A(\phi_{m+1})=\mu_{A}$, $B(\phi_{m+1})=\mu_{B}$ are two distinct numbers (the case $\mu_{A}=\mu_{B}$ is trivial since the arithmetic mean is a least-squares approximation of $h$ by rules of the type (23) below).

Consider the linear independent ($\mu_{A}\neq\mu_{B}$) vectors of $\mathbb{R}^{m+2}$:

[TABLE]

The least-squares approximation of (22) by quadrature rules of the form

[TABLE]

is equivalent to the least-squares approximation of $\mathbf{h}$ by vectors of the form

[TABLE]

That is, the minimizer of the function

[TABLE]

Denoting by $s$ the number $s=\sum_{i=0}^{m}\mu_{i}^{2}$, the minimum of $F$ is the solution of the system of normal equations

[TABLE]

whose solution is

[TABLE]

Thus, the rule that is the best approximation of (22), in the sense of least-squares, coincides with the mean rule (19). For other connections of quadrature with least-squares approximations see [5].

Example 6**.**

(A mean rule of degree 7)

Consider the Gauss-Legendre rule (17)

[TABLE]

and the open Newton-Cotes rule with 5 nodes:

[TABLE]

The rules $A(g)$ and $B(g)$ have degree $m=5$ and are both positive with

[TABLE]

The respective mean rule of $A(g)$ and $B(g)$ is

[TABLE]

This rule is positive of degree $m=7$:

[TABLE]

Note that $\gamma_{W}<\gamma_{B}<\gamma_{A}$ and $\gamma_{W}/(\gamma_{A}+\gamma_{B})/2)\simeq 0.36$, suggesting that the absolute error of $W(g)$ given by (24) is approximately $1/3$ of the arithmetic mean of errors of the rules $A(g)$ and $B(g)$.

For $I(g)=\displaystyle\int_{-1}^{1}2/(1+t^{2})\,dt=\pi$, we compare in Table 3 the errors of $A(g)$, $B(g)$ with the error of the mean rule $W(g)$.

Dividing the interval $[-1,1]$ into $n\geq 2$ equal parts, and considering the composite rules $A_{n}(g)$, $B_{n}(g)$ and $W_{n}(g)$, the gain of accuracy of the mean rule $W_{n}(g)$ relatively to the two rules of degree $5$ is numerically illustrated by the Table 3, where the respective errors are displayed. For subintervals of length $h=2/1024\simeq 0.0020$ the composite rule $W_{1024}(g)$ produces an approximation of $\pi$ with 33 significant digits.

5 Open rules of arbitrary degree with pseudorandom rational nodes

The usual approach to approximate $I(g)=\displaystyle\int_{-1}^{1}g/t)\,dt$ is to construct a quadrature rule by mean of the interpolating polynomial of a given set of nodes in $[-1,1]$. It is well-known that the resulting rules can be highly unstable when one increases the number of nodes. For instance, this is the case of the closed Newton-Cotes rules with a number of equally spaced nodes greater than $10$. In order to overcoming such instability we are going to purpose the construction of rules of high degree by combinations of starting rules of degree one.

For any integer $k\geq 1$, suppose it is given $k+1$ open rules

[TABLE]

where the symmetrical nodes $t_{i}$, for $i=1,\ldots,k$, belong to $(-1,1)$, are nonzero distinct rational numbers. Consider the combined rule

[TABLE]

In general the rules (25) are numerically stable and the assumption of the rationality of the nodes $t_{i}$, has the advantage of obtaining values free of rounding errors provided exact computation is performed. This is possible with any symbolic language system whose arithmetic is exact for rational numbers such as the Mathematica. Thus, a rational linear combination of the rules (25) will be numerically more interesting than the consideration of rules computed directly from an interpolatory rule of the $2\,k+1$ nodes.

Due to the fact we are considering nodes symmetrically distributed in $(-1,1)$, for any odd $j\geq 1$ we have

[TABLE]

Moreover, all the rules in (25) have degree one since

[TABLE]

The rules (25) cannot be of degree $2$, unless $t_{i}^{2}=1/3$, which is never the case since we are assuming the rationality of the nodes. If the combined rule (26) has a certain even degree $d$, then $W_{k}(g)$ has degree at least $d+1$ due to (27). Thus, the degree of the combined rule is always odd.

In the following proposition we prove that there exist unique rational coefficients $a_{0},a_{1},\ldots,a_{k}$ such that the linear combination (26) has odd degree $m=2\,k+1$. This result suggests that an efficient algorithm can be designed to obtain combined rules of arbitrary degree, starting from rules of degree one.

Proposition 3**.**

Let be given $k+1$ rules as in (25) with $t_{i}\in\mathbb{Q}$, and consider the linear combination (26) with coefficients $a_{0},\ldots,a_{k}$. Then,

1. (i)

the weights $a_{0},a_{1},\ldots,a_{k}$ exist and are unique; 2. (ii)

$a_{0}+a_{1}+\ldots+a_{k}=1$; 3. (iii)

the weights are rational numbers; 4. (iii)

the degree of the combined rule $W_{k}(g)$ in (26) is $m=2\,k+1$.

Proof.

For the sake of simplicity we just prove the statements for $k=1$ and $k=2$ but the result for any other $k\geq 1$ will follow by induction on $k$.

The combined rule for $k=1$ is

[TABLE]

Due to (27), one has $W_{1}(\phi_{1})=0$ and $I(\phi_{1})=0$. So, $W_{1}(g)$ has degree $d\geq 2$ if and only if the following two conditions hold:

[TABLE]

that is,

[TABLE]

The first equation is just $(ii)$. Moreover, the above system has the unique solution $a_{1}=1/(3\,t_{1}^{2}),\,\,a_{0}=1-a_{1}$. Since $t_{1}\in\mathbb{Q}$, then $a_{0},a_{1}$ are also rational and $(i)$ and $(iii)$ hold. Finally, by (27) we have

[TABLE]

which implies that the rule has degree $m\geq 3$. However, as

[TABLE]

the equation $2/3\,t_{1}^{2}-2/5=0$ does not have solution in $\mathbb{Q}$ and so the rule $W_{1}(g)$ cannot have degree $4$. So, $m=3=2\,k+1$.

Let $k=2$. Consider the combined rule

[TABLE]

Due to (27), trivially $W_{1}(\phi_{1})=0$ and $I(\phi_{3})=0$. So $W_{2}(g)$ has degree $d\geq 4$ if and only if the following three conditions hold:

[TABLE]

that is,

[TABLE]

Equivalently,

[TABLE]

The first equation of the above triangular system is just $(ii)$. Also, as the solution of the this system consists of sums, products and quotients of nonzero rationals, then $a_{0},a_{1},a_{2}\in\mathbb{Q}$ and so $(i)$ is true. Finally, by (27) we have

[TABLE]

and so that the rule has degree $m\geq 5$. However,

[TABLE]

which implies that $W_{2}(g)$ cannot have degree $6$. So, its degree is $m=5=2\,k+1$ and $(iii)$ holds. ∎

Remark 3**.**

*The system (28) is almost singular if one takes the nodes $t_{1}\neq t_{2}$ very close to the central node [math]. Therefore for non exact arithmetic one expects that the combined rule $W_{2}(g)$ will be numerically unstable for $t_{1},t_{2}\simeq 0$. Thus, one might prefer a combined rule with distinct nodes $t_{1},t_{2}$ closer to 1 rather than 0. However, for exact computations on the rationals the solution $a_{0},a_{1},a_{2}$ is exact and so, by construction, the rule $W_{2}(g)$ is stable, assuming that the function $g$ is sufficiently smooth in $[-1,1]$. *

Example 7**.**

(A combined open rule of degree 7)

Consider the rules

[TABLE]

and the respective combined rule

[TABLE]

As predicted by Proposition 3, this rule has degree $m=2\,k+1=7$. Indeed, the weights $a_{i}$ satisfy the conditions $W_{3}(\phi_{i})=I(\phi_{i})$, for $i=0,2,4,6=2\,k$, that is, they are solutions of the system

[TABLE]

Solving this system, we get

[TABLE]

One may also confirm that the sum of the weights in (29) is 1. As

[TABLE]

the combined rule as degree $m=7$ (and is a positive rule).

In Example 6 we have constructed another rule of degree 7, the mean rule $W(g)$ in (24). For the model function $g$, comparing the parameter $\gamma_{{}_{W_{3}}}\simeq 0.1202$ with the corresponding parameter $\gamma_{W}\simeq 0.014$, it is expected that the rule (24) will perform better than the rule (30) (particularly when we considerer the composite versions as in Table 3). However, the rule $W_{3}(g)$ uses only rational nodes whereas the rule (24) does not. So, due to Remark 3, we prefer the rule $W_{3}(g)$.

Composite rules of a certain degree $m$ for which the respective value of the parameter $\gamma$ is very close to zero are particularly useful. In fact, in this case it means that the nodes $t_{1},t_{2},\ldots$ have been chosen such that the combined rule is almost optimal in the sense that it acts like a perturbed rule of the optimal rule of maximum degree. In the following example, from a set of starting rules of degree $1$ of type (25)-(26), we construct a $11$-degree combined rule whose parameter $\gamma$ is small.

Example 8**.**

(A combined open rule of degree 11)

Consider the Legendre polynomial of degree 10

[TABLE]

whose roots belong to $(-1,1)$. It is well known that an interpolatory rule taking as nodes the 10 zeros of the polynomial $P_{10}(t)$ (know as a 10-point Gauss-Legendre rule) is a rule of maximal degree $m=17$.

A combined rule of the type (26), based on degree 1 rules (25) whose nodes are (a certain number of) rational approximations of the zeros of $P_{10}$ is expected to be almost as accurate as the maximal degree rule of Gauss-Legendre type. For instance, consider for the positive nodes the following rational approximations of the 5 positive roots of $P_{10}(t)$ (computed with an error222Using the Mathematica command $Rationalize[argument,10^{-16}]$ . $\leq 10^{-16}$),

[TABLE]

These nodes define the rules $Q_{0}(g),Q_{1}(g),\ldots,Q_{5}(g)$ in (25). Let us now compute the weights of the corresponding combined rule

[TABLE]

The conditions $W_{5}(\phi_{i})=I(\phi_{i}),\,\,i=0,2,4,6,8,10$ give rise to a system in the unknowns $a_{0},\ldots,a_{5}$ whose solution is given in Figure 4

The rule has degree $m=11$ (see Proposition 3) and its parameter $\gamma_{{}_{W_{5}}}$ is

[TABLE]

Comparing with the parameter $\gamma_{W}\simeq 0.014$ of the rule (24), in Example 6, (or (30) for a another rule of degree 7), one can predict that the composite rule $W_{5}(g)_{n}$, for $n$ subintervals of $[-1,1]$, will produce better numerical results than the rule $W$.

For $I(g)=\displaystyle\int_{-1}^{1}2/(1+t^{2})\,dt=\pi$, the Table 4 shows the errors of the composite rule $W_{5}(g)_{n}$, for $n=2,4,8,\ldots 1024$ subintervals. Comparing the values in this table with the ones in Table 3 we observe a real improvement of the accuracy relatively to the approximations obtained with the rule (24). Note that the last computed value of $W_{5}(g)_{1024}$ gives an approximation of $\pi$ with $60$ significant digits.

In spite of the starting rules $Q_{0}(g),\ldots Q_{5}(g)$ being rules of degree 1 only, we invite the reader to verify that a rule such as the composite Simpson rule (degree 3) is totally unable to give an approximation of $\pi$ with the referred high precision of the value $W_{5}(g)_{1024}$.

Thus, we see that the rules $Q_{0}(g),\ldots,Q_{5}(g)$ are like a basis for the combined rules $W_{5}(g)$ of degree $m=11$. It is interesting to observe what happens when we replace the element $Q_{0}(g)$ (midpoint rule) by a new one, say $\tilde{Q}_{0}(g)=g(-1)+g(1)$ (trapezoidal rule) and consider the new combined rule $\tilde{W}_{5}(g)$. The same code that produced (31) and the error values in the Table 4 gives for $\tilde{W}_{5}(g)$:

[TABLE]

and the respective errors are displayed in Table 5.

Thus $W_{5}(g)$ and $\tilde{W}_{5}(g)$ are companion rules. Consequently, for $g(t)=2/(1+t^{2})$, we have

[TABLE]

In fact, rounding $\tilde{W}_{5}(g)_{1024}$ to 61 decimal places we get

[TABLE]

whose error is

[TABLE]

5.1 High degree pseudorandom combined rules

Considering $k$ pseudorandom rational numbers in $(0,1)$, the 1-degree starting rules $Q_{i}(g)$ and the respective combined rule $W_{k}(g)$

[TABLE]

the Proposition 3 is obviously valid and so combined rules with pseudorandom nodes preserve the properties discussed before.

Example 9**.**

(A pseudorandom rule of degree $m=151$)

In order to test the stability property of the combined rule when $k$ is large, we take $k=75$. Using the Mathematica command $\verb+SeedRandom[2020]+$ we generate 76 pseudorandom rational numbers333Recurring to the command $Rationalize[RandomReal[{0,1}],10^{-4}]$.. Then, we compute the weights of the combined rule $W_{75}(g)$ in (32). Finally, we use the composite version $W_{75}(g)_{n}$, for $n$ subintervals of $[-1,1]$, applied to the function $g(t)=2/(1+t^{2})$.

The combined rule $W_{75}(g)$ has degree $m=151$ and its error parameter is

[TABLE]

We show in Table 6 the errors for the composite rule with $n=2,4,8,\ldots,1024$ subintervals, where $I(g)=\pi$. Thus, $W_{75}(g)_{1024}$ produces and approximation of $\pi$ with $507$ significant digits.