On numerical improvement of the first kind Gauss–Chebyshev quadrature rules

Abstract

One of the integration methods of the equality type is Gauss–Chebyshev quadrature rule, which is in the following form:

$$\underset{-1}{\overset{1}{\int }}\frac{f(x)}{\sqrt{1-x^2}}\mathrm{d}x=\frac{\pi }{n}\sum _{k=1}^{n}f\left(\cos \left(\frac{(2k-1)\pi }{2n}\right)\right)+\frac{2\pi }{2^{2n}(2n)!}{f}^{(2n)}(\eta )\text{,}-1<\eta <1\text{.}$$

According to Gauss quadrature rules, the precision degree of above formula is the highest, i.e. 2n − 1. Hence, it is not possible to increase the precision degree of Gauss–Chebyshev integration formulas anymore. In this way, we present a matrix proof for this matter. But, on the other hand, we claim that we can improve the above formula numerically. To do this, we consider the integral bounds as two unknown variables. This causes to numerically be extended the monomial space f(x) = x^j from j = 0, 1, … , 2n − 1 to j = 0, 1, … , 2n + 1. This means that we have two monomials more than Gauss–Chebyshev integration method. In other words, we give an approximate formula as:

$$\underset{a}{\overset{b}{\int }}\frac{f(x)}{\sqrt{1-x^2}}\mathrm{d}x\simeq \sum _{i=1}^n{w}_{i}f(x_i)\text{,}$$

in which a, b and w₁, w₂, … , w_n and x₁, x₂, … , x_n are all unknowns and the formula is almost exact for the monomial basis f(x) = x^j, j = 0, 1, … , 2n + 1. Some important examples are finally given to show the excellent superiority of the proposed nodes and weights with respect to the usual Gauss–Chebyshev nodes and weights. Let us add that in this part we have also some wonderful 2-point formulas that are comparable with 71-point formulas of Gauss–Chebyshev quadrature rules in average.

Introduction

It is known that the general form of Gauss quadrature rules are given by:

$$\underset{a}{\overset{b}{\int }}f(x)\mathrm{d}w(x)=\sum _{j=1}^n{w}_{j}f(x_j)+\sum _{k=1}^m{v}_{k}f(z_k)+R[f]\text{,}$$

where the weights $[w_j{]}_{j=1}^n\text{,}[v_k{]}_{k=1}^m$ and nodes $[x_j{]}_{j=1}^n$ are unknowns and the nodes $[z_k{]}_{k=1}^m$ are predetermined. w is also a positive measure on [a, b] (see [5], [8], [9], [10]).

The residue R[f] is determined (see for instance [15]) by:

$$R[f]=\frac{f^{(2n+m)}(\eta )}{(2n+m)!}\underset{a}{\overset{b}{\int }}\prod _{k=1}^m(x-z_k){\left[\prod _{j=1}^n(x-x_j)\right]}^2\mathrm{d}w(x)\text{,}a<\eta <b\text{.}$$

By selecting $\mathrm{d}w(x)=\frac{\mathrm{d}x}{\sqrt{1-x^2}}$, a = −1, b = 1 and m = 0 we reach the Gauss–Chebyshev formula. In other words we have:

$$\underset{-1}{\overset{1}{\int }}\frac{f(x)}{\sqrt{1-x^2}}\mathrm{d}x=\sum _{j=1}^n{w}_{j}f(x_j)+R_{\mathrm{T}}(f)\text{,}$$

where

$$R_{\mathrm{T}}(f)=\frac{f^{(2n)}(\eta )}{(2n)!}\underset{-1}{\overset{1}{\int }}{\left[\prod _{i=1}^M(x-x_i)\right]}^2(1-x^2{)}^{-\frac{1}{2}}\mathrm{d}x\text{.}$$

Calculating $[w_j{]}_{j=1}^n\text{,}[x_k{]}_{k=1}^n$ and R_T[f] in (3) yields (see [7], [12]):

$$\underset{-1}{\overset{1}{\int }}\frac{f(x)}{\sqrt{1-x^2}}\mathrm{d}x=\frac{\pi }{n}\sum _{k=1}^{n}f\left(\cos \left(\frac{(2k-1)\pi }{2n}\right)\right)+\frac{2\pi }{2^{2n}(2n)!}{f}^{(2n)}(\eta )\text{,}-1<\eta <1\text{.}$$

As regards, the integration formula (5) has precision degree 2n − 1. Note that 2n − 1 is the highest precision in all integration formulas. This subject can be proved by the following theorem:

Theorem 1

In the general n-point integration formula:

$$\underset{a}{\overset{b}{\int }}w(x)f(x)\mathrm{d}x=\sum _{j=1}^n{w}_{j}f(x_j)+E[f]\text{,}$$

the highest precision degree is 2n − 1, which w(x) is the weight function and E[f] is the error of integration formula.

Proof

It is clear that a ≠ b in relation (6) and x_i ≠ x_j if i ≠ j.

Now suppose that the integration formula (6) has the precision degree 2n (contrary hypothesis), then we have:

$$\underset{a}{\overset{b}{\int }}w(x)f(x)\mathrm{d}x=\sum _{j=1}^n{w}_{j}f(x_j)\text{,}f(x)=x^i\text{,}i=0\text{,}1\text{,}\dots \text{,}2n\text{.}$$

The expansion of above relation in form of a nonlinear system is as:

$$\left\{\begin{array}{c}{w}_1+w_2+\cdots +w_n=\underset{a}{\overset{b}{\int }}w(x)\mathrm{d}x\text{,}\hfill \\ w_1{x}_1+w_2{x}_2+\cdots +w_n{x}_n=\underset{a}{\overset{b}{\int }}\mathit{xw}(x)\mathrm{d}x\text{,}\hfill \\ w_1{x}_1^2+w_2{x}_2^2+\cdots +w_n{x}_n^2=\underset{a}{\overset{b}{\int }}{x}^{2}w(x)\mathrm{d}x\text{,}\hfill \\ ⋮\hfill \\ ⋮\hfill \\ w_1{x}_1^{2n-2}+w_2{x}_2^{2n-2}+\cdots +w_n{x}_n^{2n-2}=\underset{a}{\overset{b}{\int }}{x}^{2n-2}w(x)\mathrm{d}x\text{,}\hfill \\ w_1{x}_1^{2n-1}+w_2{x}_2^{2n-1}+\cdots +w_n{x}_n^{2n-1}=\underset{a}{\overset{b}{\int }}{x}^{2n-1}w(x)\mathrm{d}x\text{,}\hfill \\ w_1{x}_1^{2n}+w_2{x}_2^{2n}+\cdots +w_n{x}_n^{2n}=\underset{a}{\overset{b}{\int }}{x}^{2n}w(x)\mathrm{d}x\text{.}\hfill \end{array}\right.$$

But if we define the following matrices:

$$V_1≔\left[\begin{array}{cccc}\hfill 1\hfill & \hfill 1\hfill & \hfill \dots \hfill & \hfill 1\hfill \\ \hfill x_1^2\hfill & \hfill x_2^2\hfill & \hfill \dots \hfill & \hfill x_n^2\hfill \\ \hfill x_1^4\hfill & \hfill x_2^4\hfill & \hfill \dots \hfill & \hfill x_n^4\hfill \\ \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill \hfill & \hfill ⋮\hfill \\ \hfill x_1^{2n-2}\hfill & \hfill x_2^{2n-2}\hfill & \hfill \dots \hfill & \hfill x_n^{2n-2}\hfill \end{array}\right]\text{,}$$

$$V_2≔\left[\begin{array}{cccc}\hfill x_1\hfill & \hfill x_2\hfill & \hfill \dots \hfill & \hfill x_n\hfill \\ \hfill x_1^3\hfill & \hfill x_2^3\hfill & \hfill \dots \hfill & \hfill x_n^3\hfill \\ \hfill x_1^5\hfill & \hfill x_2^5\hfill & \hfill \dots \hfill & \hfill x_n^5\hfill \\ \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill \hfill & \hfill ⋮\hfill \\ \hfill x_1^{2n-1}\hfill & \hfill x_2^{2n-1}\hfill & \hfill \dots \hfill & \hfill x_n^{2n-1}\hfill \end{array}\right]\text{,}$$

$$V_3≔\left[\begin{array}{cccc}\hfill x_1^2\hfill & \hfill x_2^2\hfill & \hfill \dots \hfill & \hfill x_n^2\hfill \\ \hfill x_1^4\hfill & \hfill x_2^4\hfill & \hfill \dots \hfill & \hfill x_n^4\hfill \\ \hfill x_1^6\hfill & \hfill x_2^6\hfill & \hfill \dots \hfill & \hfill x_n^6\hfill \\ \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill \hfill & \hfill ⋮\hfill \\ \hfill x_1^{2n}\hfill & \hfill x_2^{2n}\hfill & \hfill \dots \hfill & \hfill x_n^{2n}\hfill \end{array}\right]\text{,}$$

$$D≔\left[\begin{array}{cccc}\hfill x_1\hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill \\ \hfill \hfill & \hfill x_2\hfill & \hfill \hfill & \hfill \hfill \\ \hfill \hfill & \hfill \hfill & \ddots \hfill & \hfill \hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill x_n\hfill \end{array}\right]\text{,}$$

$$b_1≔\left[\begin{array}{c}\hfill \underset{a}{\overset{b}{\int }}w(x)\mathrm{d}x\hfill \\ \hfill \underset{a}{\overset{b}{\int }}{x}^{2}w(x)\mathrm{d}x\hfill \\ \hfill ⋮\hfill \\ \hfill \underset{a}{\overset{b}{\int }}{x}^{2n-2}w(x)\mathrm{d}x\hfill \end{array}\right]\text{,}$$

$$b_2≔\left[\begin{array}{c}\hfill \underset{a}{\overset{b}{\int }}\mathit{xw}(x)\mathrm{d}x\hfill \\ \hfill \underset{a}{\overset{b}{\int }}{x}^{3}w(x)\mathrm{d}x\hfill \\ \hfill ⋮\hfill \\ \hfill \underset{a}{\overset{b}{\int }}{x}^{2n-1}w(x)\mathrm{d}x\hfill \end{array}\right]\text{,}$$

$$b_3≔\left[\begin{array}{c}\hfill \underset{a}{\overset{b}{\int }}{x}^{2}w(x)\mathrm{d}x\hfill \\ \hfill \underset{a}{\overset{b}{\int }}{x}^{4}w(x)\mathrm{d}x\hfill \\ \hfill ⋮\hfill \\ \hfill \underset{a}{\overset{b}{\int }}{x}^{2n}w(x)\mathrm{d}x\hfill \end{array}\right]\text{,}$$

$$W≔\left[\begin{array}{c}\hfill w_1\hfill \\ \hfill w_2\hfill \\ \hfill ⋮\hfill \\ \hfill w_n\hfill \end{array}\right]\text{.}$$

Then according to the properties of the weight function and this fact that a ≠ b, it is simple to show that b₁ ≠ b₃. Because, let p_n(x) be an orthogonal polynomials with respect to the weight function w(x) on [a, b] (see [16]) that is:

$$\underset{a}{\overset{b}{\int }}{p}_n(x)p_m(x)w(x)\mathrm{d}x=\left\{\begin{array}{cc}0\text{,}\hfill & \mathrm{if}n\ne m\text{,}\hfill \\ \underset{a}{\overset{b}{\int }}{p}_0^2(x)w(x)\mathrm{d}x>0\text{,}\hfill & \mathrm{if}n=m\text{.}\hfill \end{array}\right.$$

Now, it can be verified that:

$$x^n=\sum _{i=0}^n{c}_i{p}_i(x)\text{,}\mathrm{and}{x}^n-x^{n-2}=\sum _{i=0}^n{b}_i{p}_i(x)\text{,}\mathrm{such}\mathrm{that}{c}_n\ne 0\mathrm{and}{b}_n\ne 0\text{.}$$

Therefore, we should have:

$$\underset{a}{\overset{b}{\int }}{x}^n(x^n-x^{n-2})w(x)\mathrm{d}x=\underset{a}{\overset{b}{\int }}\left(\sum _{i=0}^n{c}_i{p}_i(x)\right)\left(\sum _{i=0}^n{b}_i{p}_i(x)\right)w(x)\mathrm{d}x=\sum _{j=0}^n{c}_j{b}_j\underset{a}{\overset{b}{\int }}{p}_j^2(x)w(x)\mathrm{d}x\ne 0\text{.}$$

Hence b₁ ≠ b₃, which proves our claim.

Moreover, one can replace the expanded nonlinear system (8) to the following system:

$$\left\{\begin{array}{c}{V}_{1}W=b_1\text{,}\hfill \\ V_{2}W=b_2\text{,}\hfill \\ V_{3}W=b_3\text{.}\hfill \end{array}\right.$$

On the other hand we have:

$$V_2=V_{1}D\text{,}$$

$$V_3=V_1{D}^2\text{,}$$

where V₁ is an inversable Vandermonde matrix, because if i ≠ j then x_i ≠ x_j(see [11]).

Now noting to relations (20), (21) and also this matter that V₁ is inversable yields:

$$V_1{\mathit{DV}}_1^{-1}{b}_2=b_1\text{.}$$

Also by combining (20), (21) we get:

$$V_1(I-D)W=b_1-b_2\text{.}$$

Hence the following result will be derived:

$$(I-D)W=V_1^{-1}(b_1-b_2)\text{.}$$

On the other hand, combining (20), (22) yields:

$$V_{1}D(I-D)W=b_2-b_3\text{.}$$

This leads to have the following result:

$$D(I-D)W=V_1^{-1}(b_2-b_3)\text{.}$$

If we put the value (I − D)W from (25) in the relation (27) then:

$${\mathit{DV}}_1^{-1}(b_1-b_2)=V_1^{-1}(b_2-b_3)\text{.}$$

Consequently

$$V_1{\mathit{DV}}_1^{-1}(b_1-b_2)=(b_2-b_3)\text{.}$$

But, V₂ = V₁D, therefore:

$$V_2{V}_1^{-1}{b}_1-V_2{V}_1^{-1}{b}_2=b_2-b_3\text{.}$$

Now, noting to (20) results:

$$V_2{V}_1^{-1}{b}_2=b_3\text{.}$$

And (21), (23), (31) give us finally:

$$b_3=b_1\text{,}$$

which is a contradiction. So the highest precision degree is at most 2n − 1 and the theorem is proved.  □

Thus, for Chebyshev integration formula we can write:

$$\underset{-1}{\overset{1}{\int }}\frac{f(x)}{\sqrt{1-x^2}}\mathrm{d}x=\sum _{j=1}^n{w}_{j}f(x_j)\text{,}f(x)=x^i\text{,}i=0\text{,}1\text{,}\dots \text{,}2n-1\text{.}$$

By applying above relation one can reach a nonlinear system with 2n equations and 2n unknowns, so that the unknowns are $[x_j{]}_{j=1}^n\mathrm{and}[w_j{]}_{j=1}^n$. Hence, solving this nonlinear system results the integration formula (5). Now if we suppose a and b (i.e. bounds of integration formula (33)) to be two additional unknown variables, then by extending the basis space {1, x, x², … , x²ⁿ⁻¹} to {1, x, x², … , x²ⁿ⁺¹} we have:

$$\underset{a}{\overset{b}{\int }}\frac{f(x)}{\sqrt{1-x^2}}\mathrm{d}x=\sum _{j=1}^n{w}_{j}f(x_j)\text{,}f(x)=x^i\text{,}i=0\text{,}1\text{,}\dots \text{,}2n+1\text{.}$$

Certainly, the nonlinear system (34) does not have any analytical solution, because as we said in the Theorem 1, the highest precision degree in general integration formulas is 2n − 1.

On the other hand, let us add that this nonlinear system can be solved numerically, which causes to increase the monomial space two degrees. This would be the main approach behind our work. In the next section, some examples are given to show the numerical superiority of the presented rules with respect to the numerical results of Gauss–Chebyshev quadrature rules.

But, one of the applications of Gauss–Chebyshev integration rule is the integration of trigonometric functions, because by considering the relation:

$$\underset{0}{\overset{\pi }{\int }}f(\cos (\theta ))\mathrm{d}\theta =\underset{-1}{\overset{1}{\int }}\frac{f(x)}{\sqrt{1-x^2}}\mathrm{d}x\text{,}$$

one can transform the trigonometric functions integration rules to Gauss–Chebyshev quadrature formulas (see [14]). This also will be shown in the numerical examples of the next sections.