Документ взят из кэша поисковой машины. Адрес оригинального документа : http://shamolin2.imec.msu.ru/art-104-1.pdf Дата изменения: Sat Jul 21 09:55:19 2012 Дата индексирования: Mon Oct 1 20:35:59 2012 Кодировка:

ISSN 0027-1322, Moscow University Mathematics Bul letin, 2010, Vol. 65, No. 4, pp. 135­139 . c Al lerton Press, Inc., 2010. Original Russian Text c R.R. Aidagulov, M.V. Shamolin, 2010, published in Vestnik Moskovskogo Universiteta, Matematika. Mekhanika, 2010, Vol. 65, No. 4, pp. 3­7.

Integration Formulas of Tenth Order and Higher
R. R. Aidagulov and M. V. Shamolin
Moscow State University, Institute of Mechanics, Leninskie Gory, Moscow, 119991 Russia Received Decemb er 18, 2009, in final form, Septemb er 22, 2009 Abstract--Nowadays, due to the considerable growth of computer capacities, the development of more efficient quadrature formulas may seem unnecessary. However, if the calculation of each integrand value requires much computational time or we have to study the dep endence of the integral on a large numb er of parameters the integrand is determined through, then it is necessary to use more efficient formulas.

DOI: 10.3103/S0027132210040017
a

Quadrature formulas of 8th order for calculation of integrals of the form
b

f (x) dx are presented (and

used) in reference-bo oks and program packages. In fact, it is easy to calculate all the co efficients of quadrature formulas of Gauss type for the order 10 and higher. In particular, for "symmetric weights" we have
a

(x)f (x) dx, where (x) is an even weight function; the degree of the polynomial for the calculation of
-a

weights is decreased by two times, more precisely, P (x) = x f (x2 ), where P (x) is the polynomial of degree m = 2k + , k = [m/2], for determination of the weights, f (x) is a polynomial of degree k (here [... ] is the integral part). In all practically interesting cases of approximate calculation of integrals, one can cho ose an even weight function so that this expression is a Gauss integral formula with this weight. In the case a = we assume that all integrals for power functions are determined. Thus, taking into account the evenness of the weight function, we get that integrals of odd functions (odd powers are such functions) are equal to zero. Since the general theory does not depend on the evenness, at the beginning we present our arguments for an arbitrary weight function and search for the integration formula in the following form:
a m

I (f ) =
b

(x)f (x) dx = Sm (f )+ Rm (f ),

Sm (f ) =
j =1

a

mj

f (xmj ).

The polynomial Pm (x) =
j

(x - xmj ) is said to be the polynomial of no des. We perform the choice of

2m variables am1 ,... ,amm , xm1 ,...,xmm so that for the function f (x) = xj , j = 0, 1,... , 2m - 1, the error of the integral sum Rm (f ) vanishes. Such formulas are called Gauss quadrature formulas [1; 2, Section 2]. In this case the error order of a quadrature function is estimated by the value M2m m - aj x2m /(2m)!, j
j

where M

2m

= sup |f
x

(2m)

(x)|.

This implies that the polynomials of nodes are mutually orthogonal in the integral metrics:
a

(x)Pm (x)Pn (x) dx = 0,
b

m = n.

If (x) is a nonnegative function with a positive integral, then the orthogonal polynomials with the coefficient 1 at the highest degree are uniquely determined by the orthogonalization metho d (the necessity of the condition P0 (x) 1 follows from the consideration of integrals of f (x) = Pm (x)). The theory of orthogonal polynomials implies that these polynomials have exactly m roots in the studied interval and hence they are polynomials of no des. Taking the function f (x) from the condition f (x) = Pm (x)/(x - xmj ) =
i=j

(x - xmi ),

135


136 we determine the weights

AIDAGULOV, SHAMOLIN

a

a

mj

=

1 fm (xmj )

(x)
b i=j

(x - xmi ) dx.

(1)

In order to calculate the weights amj , one can obtain a polynomial formula of degree m - 1: amj = m-1 (xmj ). From the computational viewpoint, it is convenient to express the values of no des and weights through the moments up to the order 2m. The metho d of determination of Gauss formulas was perfectly well described for the case (x) 1 in [1; 2, Section 2]. Here we bring the corresponding calculations to particular formulas for an arbitrary symmetric weight function (x). Since in the even case all moments of odd order are equal to zero and we consider only them in the sequel, we determine the moments by the formulas
a

2 n =
-a

(x)x2n dx,

n = 0, 1,... .

Note that, due to the evenness of the and antisymmetric no des xj + xm+1-j shifting the old one so that the positive x-1 = -x1 ,...,x-k = -xk , additionally This gives a j = 2 0 ,
j

weight function (x), we obtain symmetric weights aj = am+1-j = 0. Therefore, it is better to introduce a new enumeration by no des were represented as x1 ,x2 ,... ,xk and the negative ones as assume x0 = 0 in the case of odd m. n = 1, 2,... ,m - 1,

yj = x2 , j
j k

n aj y j = n ,

k = [m/2].

(2)

Using the cases n = 1, 2,... ,k for o dd m and n = 0, 1,... ,k - 1 for even m, we come to the following expression for aj : (-1)i+j ij , (3) i-1+ m-2k aj =1+ = yj i where =
i
(yj - yi ) is the Vandermonde determinant, ij is the determinant obtained from by

cancellation of the ith row and j th column. Substituting expressions (3) into equations (2), we have
l l-1+

cj = l

j -1+

,

cj = l

j l ,

(4)

where the determinant j is obtained from the Vandermonde determinant by replacement of the lth row, l j j j l l l namely, the row y1-1 ,y2-1 ,...,yk-1 , is replaced by the j th row y1-1 ,y2-1 ,...,yk-1 . It is clear that cj = j / for l k , j k , and so equations (4) give nothing. For j = k + 1,... , 2k l l we obtain a system of k Diophantine equations for determination of the unknowns y1 ,... ,yk . It is also clear that cj is a symmetric polynomial of degree j - l of the variables y1 ,...,yk and the maximal degree l at the variable y1 is equal to j - k , and the co efficient at the same variable is (-1)k-l . Correspondingly, ck+1 = (-1)k-l k+1-l , cj = (-1)k-1 k cj -1 , where k are some numbers to be determined. For example, 1 k l from these equations we may obtain the symmetric functions 1 = y1 + y2 + ... + yk , 2 =
i
yi yj ,

... ,

k =
j

yj .

Using linear combinations, from (4) we get the following linear equations with respect to j : (-1)j
0j k k+j -1+ j

= 0,

m = k + j - 1+ = k +1,... , 2k.

MOSCOW UNIVERSITY MATHEMATICS BULLETIN

Vol. 65

No. 4

2010


INTEGRATION FORMULAS OF TENTH ORDER AND HIGHER

137

This gives the following equation for the determination of squares of no des: 1 y ... y k 1 2 ... k+1 = 0, ... ... ... ... k k+1 ... 2k 1 y ... y k 0 1 ... k ... ... ... ... k-1 k ... 2k-

= 1,

= 0,
1

= 0.

Calculating the roots of this equation, from (2)­(4) we get the values of the weights a1 ,... ,ak and a0 = 0 - a1 - ... - ak . Also note that setting the weights of all yi equal to one and the weights of i equal to i, we obtain everywhere uniform formulas, and the weights ai are equal to zero. There also exists a symmetry for the change of the weight function (x) by (x)x2 with the increment of weights of the variables i by one. Consider small odd values in the general case in more detail. The case k = 0 is reduced to the following formula: (x)f (x) dx 0 f (x). For k = 1 we come to formula (1) with the following values: y1 = 2 , 1 x1 = y1 , a1 = 1 2 = 1, y1 2 a0 = 0 - 2 1 =0 1 2
1 2

2 .

Further, for k = 2 we have (2 - 1 3 )y 2 - (3 2 - 1 4 )y +(2 - 2 4 ) = 0, 2 3 D = 2 2 - 61 2 3 4 - 32 2 +43 4 +41 3 , 14 23 2 3 a1 = y 2 1 - 2 , y1 a2 = 2 - y 1 1 , y2 0 1 a3 = 1 2 2 3
2 3 4

2 3

3 4

.

The case k = 3 results in the following values: 1 2 2 3 3 4 a1 =
3 4 5

1 2 y 3 - 2 3 3 4 a2 =

4 5 6

1 3 y 2 + 2 4 3 5

4 5 6

2 3 y - 3 4 4 5

4 5 6

= 0, 1 y1 y2 - 2 (y1 + y2 )+ 3 , y3 (y3 - y1 )(y3 - y2 )
4 5 6

1 y2 y3 - 2 (y2 + y3 )+ 3 , y2 (y2 - y1 )(y3 - y1 )

-1 y1 y3 + 2 (y1 + y3 ) - 3 , y2 (y2 - y1 )(y3 - y2 )
0 1 2 3

a3 =

1 2 - 2 1 + 3 a0 = 0 - = 3



1 2 3 4



2 3 4 5



3 4 5 6

2 3 3 4 4 5

.

Now consider particular weight functions and calculate the formulas of the 10th and 14th orders of accuracy. For the 14th order of accuracy the weights and no des are ro ots of cubic equations and can be expressed by Cardano formulas. Since those expressions are bulky, we present here only approximate values. 1) In some applications we have to calculate the integrals


f (x)exp(-x2 /2) dx.
-

In this case it is convenient to take (x) = exp(-x2 /2), a = . The corresponding orthogonal polynomials (normalized polynomials of no des) are called Sonin polynomials, and the polynomials forhe weight function t (x) = exp(-x2 ) connected to the previous ones by the contraction of the axis x by 2 times are called Hermitian or Chebyshev­Hermitian [3, Ch. 1; 4, Ch. 4]. The equality n = /2(2n - 1)!! is valid. For m = 5 (k = 2) we obtain the no des and weights for the 10th order of accuracy: 7+2 10 x1 = 5 - 10 = -x-1 , a1 = a-1 = , 30 2
MOSCOW UNIVERSITY MATHEMATICS BULLETIN Vol. 65 No. 4 2010


138 x2 = 5+ 10,

AIDAGULOV, SHAMOLIN

a2 = a

-1

7 - 2 10 = 30

, 2

a0 =

8 15

, 2

and for m = 7 (k = 3) the corresponding polynomial has the form y 3 - 21y 2 + 105y - 105 = 0, and in this case the nodes and weights are x1 = 1.154405394739968, x2 = 2.366759410734541, x3 = 3.750439717725742,

a0 = 0.5898718274019, a1 = 0.601899548885595, a2 = 0.05124934362178949, a3 = 0.01029341740621618. 2) The following integrals of hypergeometric form often o ccur in practice [4, Ch. 4]


Re
-

(x - a1 )(a2 - x) ... (an - x) e (b1 - x)(b2 - x) ... (bn - x)
-cx

bi -cx

dx =
ia i

fi (x) dx (x - ai )(bi - x)

,

fi (x) = (x - ai ) e

(bi - x) (x - ai )

j

aj - x , bj - x

a1 < b1 < a2 < b2 < ... < an < bn .
1 1-x

In order to calculate such integrals efficiently, one should use formulas for the weight n = (2n - 1)!! 0 = . 2 (2n)!! 2

2

, a = 1. Then

This case is the simplest from the computational viewpoint because the quadrature formula takes the simple form
1

-1

f (x) dx = 2 m 1-x

f (xj ),
j

xl = - cos

(2l - 1) . 2m

(5)

The consideration of this case is reduced to the following equations
1 1 2 2

I (f ) =
-1

f (x) dx = 2 1 - x2 S (f ) = m

sin2k
-1

x dx = 2 2
0

cos2k

x dx = 2k+1 2 2
j C2k

j C2k 0j k 0

cos

2(k - j )x dx, 2

cos2k
0lm

(2l - 1) = 2k 2m 2

0j k

1 m

cos
0lm

(k - j )(2l - 1) . m

Therefore, it remains to prove that
2

1 2
0

cos(n x) dx =

1 m

cos
0lm

n(2l - 1) , 2m

n < m.

(6)

The validity of equality (6) follows from the consideration of the following cases: 1) for n = 0 in the left- and right-hand sides of equality (6) we have 1; 2) for 0 < n < m in the left- and right-hand sides of equality (6) n we get 0. The latter takes place due to the auxiliary multiplication of both the sides by cos 2m = 0. Actually, it is not difficult to calculate the difference between the integral and the approximative sum for k m: I (x2k ) - S (x2k ) = 22k
k- C2k 1lk/m lm

.

In this case (the same) weights and no des are represented by simple formulas (5) and we have no need to consider separate cases. These normalized polynomials of nodes are called Chebyshev polynomials of the first kind. 3) (x) 1, a = 1. In this case n = 1/(2n + 1) and the polynomials of nodes (which are Legendre polynomials up to normalization) have the form n! dn (x2 - 1)n . (2n)! dxn
MOSCOW UNIVERSITY MATHEMATICS BULLETIN Vol. 65 No. 4 2010


INTEGRATION FORMULAS OF TENTH ORDER AND HIGHER

139

Correspondingly, for k = 2 we get x1 = 35 - 280 , 63 x2 = 35 + 280 , 63 322 + 13 70 a1 = , 900 322 - 13 70 a2 = , 900 a0 = 128 ; 225

for k = 3 the equation for no des has the form 35 - 315y + 639y 2 - 429y 3 = 0 and the no des and weights are x1 = 0.4058451513773965, x2 = 0.741531185599394, x3 = 0.949107912342758,

a0 = 0.0794161895996775, a2 = 0.133350981823048378, For calculation of the integral R0 0.72676,
1 -1

a1 = 0.3818300505051187, a3 = 0.405402778072155431. these formulas give the following errors: R2 7.021 · 10
-8

cos

2

x dx =

4

R1 8.842 · 10

-4

,

,

R3 < 10

-12

.

In principle. the latter estimates allow us to use the integration formulas of the 10th and 14th orders on a computer (and also, as before, on programmable calculators) for the latter two cases. ACKNOWLEDGMENTS The work was supported by the Russian Foundation for Basic Research (pro ject no. 08­01­00231-). REFERENCES
1. 2. 3. 4. S. A. G. A. M. Nikolskii, Quadrature Formulas (Fizmatgiz, Moscow, 1958) [in Russian]. N. Krylov, Lectures on Approximate Calculations (Moscow, 1950) [in Russian]. Szego, Orthogonal Polynomials (AMS, Providence, RI, 1975; Mir, Moscow 1962). F. Nikiforov and V. B. Uvarov, Special Functions of Mathematical Physics (Nauka, Moscow, 1976) [in Russian].

Translated by V. Valedinskii

MOSCOW UNIVERSITY MATHEMATICS BULLETIN

Vol. 65

No. 4

2010