- Email: [email protected]
On the approximate calculation of multiple integrals
Journal of Complexity (
)
–
Contents lists available at ScienceDirect
Journal of Complexity journal homepage: www.elsevier.com/locate/jco
On the approximate calculation of multiple integrals Nikolai Sergeevich Bakhvalov Department of Computational Mathematics, Moscow State University, Russian Federation
article
info
Article history: Available online xxxx Keywords: Optimal quadrature formulas Monte Carlo methods
abstract When approximately calculating integrals of high dimension with the Monte Carlo method, one uses fewer values of the integrand than when calculating with the help of classical deterministic quadrature formulas. However, the error estimation for the Monte Carlo method does not depend on the smoothness of the integrand. This suggests that it is possible to obtain methods that give a better order of convergence in case of smooth functions. © 2015 Published by Elsevier Inc.
1. Lower estimation of the possible error We show a lower bound for the maximum error of evaluating integrals for functions of a given smoothness. Consider the following problem: suppose we have the class of all functions fj (n) (n = 1, 2, . . . , 2N) that take the values ±1, and one wants to evaluate the arithmetic mean of a function from this class 2N
Sj =
n =1
fj (n)
.
2N We have the following lemma. Lemma 1. If a deterministic method of approximately evaluating the arithmetic mean uses information about the function values at no more than N points, then there is a function from the class for which the difference between the computed arithmetic mean and the true value is at least 1/2. Indeed, suppose that the arithmetic mean was evaluated using information about values of f at points i1 , i2 , . . . , im (m ≤ N), and that f (ik ) = αk (k = 1, . . . , m). Consider functions ft (n) and fτ (n) such that ft (ik ) = fτ (ik ) = αk (k = 1, . . . , m) and ft (n) = −fτ (n) = 1 at the other points. http://dx.doi.org/10.1016/j.jco.2014.12.003 0885-064X/© 2015 Published by Elsevier Inc.
2
N.S. Bakhvalov / Journal of Complexity (
)
–
Since |St − Sτ | ≥ 1 and the calculated approximate values S t and S τ of the arithmetic means are the same, then max(|St − S t |, |Sτ − S τ |) ≥
1 2
.
Definition 1. A method of evaluating an integral (arithmetic mean) is called non-deterministic if the coordinates of each successive knot are random functions that depend on the integration domain and all values obtained previously in the computational process (coordinates of the knots, function values and its derivatives, the number of knots, etc.), and the approximate value of the integral (arithmetic mean) is a random function depending on all the obtained values. This definition applies to the Monte Carlo method, the method proposed in [4], methods described in Sections 2 and 4 of the present paper, and many others. Lemma 2. For any nondeterministic method of evaluating the arithmetic mean that uses information about function values at no more than N points, there is a function in the given class for which the average of the absolute error is at least C0 N −1/2 where C0 is an absolute constant. We will say that an event Kγ ,
α1 , . . . , αmγ γ = β1 , . . . , βmγ
,
mγ ≤ N, αk = ±1, took place in the process of evaluating the arithmetic mean, if information about values of the function at points β1 , . . . , βmγ was used, and f (βk ) = αk
(k = 1, . . . , mγ ).
Since the approximate value of the arithmetic mean (integral) can be written using a finite number of binary digits, it can take no more than countably many values Is (s = 1, 2, . . .). Denote by πγj the probability of the event Kγ in case the integral is evaluated for the function fj . γ
By σs we denote the probability that the approximate evaluation of the integral equals Is under the condition that the event Kγ took place. We associate with each function fj from the considered class the probability pj = 2−2N ; it is clear that
pj = 1.
j
Set
ρj =
γ
ρ=
πγ j
γ
σs |Is − Sj | ,
(1)
s
ρj pj .
(2)
j
Obviously, ρj is the average of the absolute value of the error when calculating the integral of the function fj . Plugging (1) into (2) we have
ρ=
γ
πγ pj j
s
j
where
qγ ,s =
j
πγj pj |Is − Sj | j
πγj pj
.
γ
σs qγ ,s ,
(3)
N.S. Bakhvalov / Journal of Complexity (
)
–
3
We will say that γ ∈ j if πγj ̸= 0. Consider those γ for which πγj ̸= 0 for at least one value of j (all πγj can be zero only when we do not use information in any of the points βk ). Obviously, the quantities πγj depend only on function values at points β1 , . . . , βmγ , and therefore, j
j
if γ ∈ j1 , γ ∈ j2 then πγ1 = πγ2 = πγ ̸= 0. Since πγj ̸= 0 only for 22N −mγ functions fj satisfying γ ∈ j, i.e., f (βk ) = αk (k = 1, . . . , mγ ), and all pj are equal to each other, then we have γ ∈j pj = 2−mγ . Hence qγ ,s = 2−(2N −mγ )
|Is − Sj |.
γ ∈j
All the functions fj satisfying the condition γ ∈ j can be divided into pairs fµk and fνk such that at points different from β1 , . . . , βγ is fµk (n) = −fνk (n). Obviously,
|Is − Sµk | + |Is − Sνk | ≥
1 γ γ |S µk | + |S νk | , 2N
where γ
Sj =
fj (n).
n̸=βk k=1,...,mγ
This implies qγ ,s ≥
t2N −mγ N
=
1
γ
N
−(2N −mγ )
2
|S j |
2
γ ∈j
.
The quantity t2N −mγ coincides with the average deviation of the binomial distribution with parameters (2N − mγ , 1/2, 1/2). By the central limit theorem, for large 2N − mγ we have
t2N −mγ =
2
π
2N − mγ + o
2N − mγ
.
Since mγ ≤ N, we conclude that for all N is t2N −mγ ≥ C0 N 1/2 . Therefore qγ ,s ≥ C0 N −1/2 . Using (3)
we obtain ρ ≥ C0 N −1/2 . As a consequence, there exists such j0 that ρj0 ≥ C0 N −1/2 , and this is what we were supposed to show. We will pass to estimating from below the maximum error of evaluating the integral I (f ) =
1
1
f (x1 , . . . , xm ) dx1 . . . dxm ,
··· 0
0
where f ∈ H (p, A, λ) and H (p, A, λ) is the space of all f ∈ C p ([0, 1]m ) such that all partial derivatives up to the order p ∈ N are bounded by A and the derivatives of order p are Hölder continuous with exponent 0 < λ ≤ 1 and constant A > 0, see [3]. The following considerations essentially coincide with the considerations that led to obtaining a lower bound for ε -entropy of the class of functions H (p, A, λ) [6]. Let N be an integer, n0 =
1 (2N ) m + 1,
p+λ) ψn0 (x1 , . . . , xm ) = Pp+λ,m n−( 0
m
[n0 xk (1 − n0 xk )]p+λ ,
k=1
where Pp+λ,m is a constant, such that ψn0 (x1 , . . . , xm ) ∈ H (p, 1, λ) for all n0 and 0 ≤ x1 , . . . , xm ≤ 1/n0 .
4
N.S. Bakhvalov / Journal of Complexity (
)
–
Let Q1 , . . . , Q2N be any 2N different cubes, nk /n0 ≤ xk ≤ (nk + 1)/n0 (k = 1, . . . , m; 0 ≤ nk < n0 ; nk an integer). Consider all the possible functions ϕj such that on each of the cubes
n1 nm ϕj = ±Aψn0 x1 − , . . . , xm − , n0
n0
and for the other points ϕj = 0. Obviously, for 0 ≤ x1 , . . . , xm ≤ 1, all the functions ϕj ∈ H (p, A, λ). There is a one-to-one correspondence between the functions ϕj and fj (n). Therefore, using Lemmas 1 and 2 we obtain the following theorem. Theorem 1. If the integral I (f ) is computed by a deterministic algorithm using information about the function values and its derivatives at no more than N points then there is a function f ∈ H (p, A, λ) for which the true value of the integral differs from the one computed with a given method by at least ′ −(p+λ)/m Cm , and if a nondeterministic method is used then there is a function f ∈ H (p, A, λ) for ,p+λ AN ′′
which the average of the absolute value of the error is larger than Cm,p+λ AN
1 − (p+λ)/m+ 2
.
Here Cm,p+λ , Cm,p+λ > 0 are some constants that depend on C0 and on the quantity ′
′′
1
1
ψ1 (x1 , . . . , xm ) dx1 . . . dxm .
··· 0
0
2. A method for the approximate evaluation of the integrals The classical Monte Carlo method relies on taking a sequence of independent random points from the domain of integration and computing the arithmetic mean of the function values at those points, see [1]. In case of smooth integrands, it seems better to take coordinates of the points depending on each other. An example of such a method for computation of integrals is the method proposed in [4]. Below we describe a method that is based on the same idea and allows, using function values at N points, to obtain the value of the integral with error of order O N
p+λ 1 + m 2
−
absolute value of the error of order O N
−
p+λ m
and with the average of the
.
Let t be the minimum integer such that 2t ≥ p + λ, 0 ≤ ξ ≤ 1, and let f ∈ H (p, A, λ) if −t /n0 ≤ x ≤ t /n0 . ξ Denote by Ω2t −1 (x) the Lagrange interpolation polynomial of degree 2t − 1 that agrees with the ξ ξ +t −1 function f at knots ± n , . . . , ± n . 0 0 It is known that ξ
Ω2t −1 (x) =
t −1
Wβ (x, ξ ) f
β +ξ n0
β=0
= f (x) + O
A p+λ
n0
+ W−β (x, ξ )f
β +ξ − n0
.
(4)
Integrating (4) in the intervals from −t /n0 to t /n0 we obtain t n0
−
t n0
f (x) dx =
t −1
β=0
πβ (ξ ) f
β +ξ n0
+f
−
β +ξ n0
+O
A 1+p+λ
n0
.
(5)
N.S. Bakhvalov / Journal of Complexity (
)
–
5
k k Let us perform a series of mnm 0 trials with outcomes ξn1 ...nm , where ξn1 ...nm are independent random quantities distributed uniformly on the interval [0, 1]. We will assume that f ∈ H (p, A, λ) on the whole of Rm and that f is periodic with respect to each variable with period 1. We set t −1
σn1 ...nm =
m
πβk ξnk1 ...nm
β1 ,...,βm =0 k=1
×
f
n1 + (−1)α1 (ξn11 ...nm + β1 ) n0
α1 ,...,αm =0,1
,...,
m nm + (−1)αm (ξn1 ...nm + βm )
.
n0
From (5) it follows that
σn1 ...nm =
n 1 +t n0
n m +t n0
···
n 1 −t n0
f dx1 . . . dxm + O
n m −t n0
A m+p+λ
n0
.
(6)
Taking the sum of (6) with respect to all nk (0 ≤ n1 , . . . , nm < n0 ) and dividing by (2t )m we obtain n 0 −1
σ =
n1 ...nm =0
σn1 ...nm
= I (f ) + O
(2t )m
A
.
p+λ
n0
(7)
From (6) we have
D σn1 ...nm = O
A2
.
2(m+p+λ)
n0
Due to independence of ξnk1 ...nm the values σn1 ...nm are also independent and therefore D(σ ) =
n 0 −1
1
(2t )2m
D σn1 ...nm = O
n1 ...nm =0
A2 m+2(p+λ)
n0
.
(8)
Now we show that M σ = I (f ). For short, we consider only m = 1. We have: M σn1 =
n1 +β+1 n0
t −1
n1 +β n0
β=0
−
t −1
β=0
− −
πβ (n0 x − n1 − β)f (x) dx
n1 −β−1 n0
n1 −β n0
πβ (|n0 x − n1 − β|) f (x) dx.
(9)
ξ
Since for f (x) = 1 is Ω2t −1 (x) = 1, then t −1
πβ (ξ ) =
β=0
t n0
.
(10)
Taking the sum of (9) with respect to n1 in the intervals from 0 to n0 − 1 and using (10) we obtain: Mσ =
n 0 −1 1
2t
n 1 =0
M σn1 = I (f ).
6
N.S. Bakhvalov / Journal of Complexity (
)
–
In this method of evaluating the integral one uses function values at N ≍ nm 0 points. Since
M (σ − I )2 ≥ M |σ − I |, from (7) and (8) we have
p+λ σ − I = O AN − m ,
M |σ − I | = O AN −
p+λ 1 m +2
.
By comparing these estimates with the results of Theorem 1 one can see that there is no method of calculating the integral whose error estimate has better order. 3. A theorem on existence of good quadrature formulas Let us consider the problem of evaluating the integral in case f ∈ Eγ (A), where Eγ (A) (γ > 1) is the class of periodic functions with period 1 with respect to each variable, such that the Fourier coefficients satisfy the condition
|⃗cn | ≤
A
|||⃗n|||γ
.
⃗ = (n1 , . . . , nm ), p = max(1, |p|), |||⃗n||| = k=1 nk , (⃗n, a⃗) = k=1 nk ak , Here and in the sequel n and σm , χm,γ are some constants. N.M. Korobov [7,8] dealt with derivation of formulas that are useful for computation of integrals for such functions. Using some other arguments he showed classes of quadrature formulas with error estimates that differ from the estimates obtained in this section by multipliers of the form lnα N. We will consider quadrature formulas of the form m
N −1
ΣNP
(f ) =
f (jP )
j =0
N
,
where P =
a
1
N
,...,
am N
m
, ak -integer.
We have [8]
RNP (f ) = ΣNP (f ) − I (f ) =
′ cn⃗ , (⃗n,⃗a)≡0 (mod N )
(11)
′
⃗ = 0⃗ is excluded. where means that in the summation n ⃗ we mean a vector with integer coordinates and we write In what follows, by n ⃗ (⃗n, P ) ≡ 0 (mod 1) = n⃗ ̸= 0: ⃗ (⃗n, a⃗) ≡ 0 (mod N ) . ⃗ ̸= 0: LP := n
We will prove the following theorem.
⃗ ∈ LP , then Theorem 2. Let M > 1. If |||⃗ n||| ≥ M for all points n m−1 P M R (f ) ≤ χm,γ · A ln N γ
M
for all f ∈ Eγ (A). For the proof of this theorem we need the following two lemmas. Lemma 3. Let m PM
:=
m
[ak , bk ] ⊂ R : ak < bk , k = 1, . . . , m, and m
k=1
m
(bk − ak ) < M
k=1
be the set of all axes-parallel hyperrectangles of volume less than M. Then, under the assumption of Theorem 2, we have
LP ∩ Ω ≤ 1 for all Ω ∈ P m . M
N.S. Bakhvalov / Journal of Complexity (
)
–
7
⃗, n⃗′ ∈ LP , then, obviously, we also have n⃗′′ = Proof. Assume that Ω contains 2 different points n ′ ′ ⃗ − n⃗ ∈ LP . By assumption, this means |||⃗n − n⃗||| ≥ M. But since Ω ∈ PMm and n⃗, n⃗′ ∈ Ω we know that n ⃗||| ≥ M: a contradiction. the volume of Ω is strictly smaller than M and larger than |||⃗ n′ − n Define the hyperbolic cross ΓMm by
ΓMm := n⃗: |||⃗n||| < M . We will show the following lemma. Lemma 4. Let M > 1 and ν ∈ N. The set ΓνmM can be covered using σm ν lnm−1 M ν hyperrectangles from m PM . Proof. We perform the proof using induction on m. It is clear that for m = 1 the assertion of the lemma holds with σm = 3. Suppose it is proven for some m. We split the set ΓνmM+1 into parts by the hyperplanes nm+1 = 0, nm+1 = ±2i ν , i = 0, 1, . . . , ⌊log2 M ⌋. m ⃗ ∈ ΓνmM+1 . Therefore the intersection of the set If 2i ν < |nm+1 | ≤ 2i+1 ν , then k=1 nk < M for n 2i
ΓνmM+1 with the hyperplane nm+1 = ν 2i + 1 can be covered using at most qi := σm lnm−1 (i)
(i)
M 2i
hyper-
m rectangles from PM , say Q1 , . . . , Qqi . We define the (m + 1)-dimensional hyperrectangles from /2i (i) m+1 i Pν M by Qℓ × [ν 2 , ν 2i+1 ], ℓ = 1, . . . , qi . It is clear that each such hyperrectangle can be covered m+1 using ν hyperrectangles of the type PM . As a consequence, all the points of the set ΓνmM+1 , such that m+1 2i ν < nm+1 ≤ 2i+1 ν , can be covered using ν qi hyperrectangles of the type PM .
If |nm+1 | ≤ ν , then, by inductive assumption, all points of the set ΓνmM+1 such that nm+1 = j ∈ N can be covered using q∗,j := σm
ν j
+ 1 lnm−1 (M ⌈ νj ⌉) (∗,j)
m hyperrectangles from PM , say Q1
ν
(∗,j)
Qℓ
j) , . . . , Qq(∗, ∗,j . Clearly, the collection of (m + 1)-dim. hyperrectangles
× [j − 1, j]: ℓ = 1, . . . , q∗,j
j =1 m+1 contains every point of ΓνmM+1 with 0 ≤ nm+1 ≤ ν , and each element is from PM .
From this it follows that the total number of hyperrectangles sufficient to cover the set ΓνmM+1 does not exceed 2
⌊log M ⌋ 2
νσ
m−1 M m ln 2i
+
ν j =1
i =0
σm
ν j
m−1
+ 1 ln
≤ 2 σm ⌊log2 M ⌋ν ln ′
m−1
M + 2νσm ln
m−1
Mν
ν
(M ⌈ j ⌉)
ν 1 j=1
j
≤ σm+1 ν lnm M ν. (The additional 2 comes from the absolute value in the restriction on nm+1 .) Thus the lemma is proven. We proceed with the proof of the theorem. We denote by Tν the number of solutions of the equation
(⃗n, a⃗) ≡ 0 (mod N ) that belong to ΓνmM . Combining Lemmas 3 and 4, we obtain Tν ≤ σm ν lnm−1 M ν.
8
N.S. Bakhvalov / Journal of Complexity (
)
–
From (11) it follows that ∞ P R (f ) ≤ | c | ≤ ⃗ n N ν=1
⃗∈LP n ∞
=
ATν+1
A
(M ν)γ
1
−
(M ν)γ
ν=1
(Tν+1 − Tν ) 1
(M (ν + 1))γ
−
AT1 Mγ
∞ ∞ ATν+1 (ν + 1)γ − ν γ A γ σm (ν + 1) lnm−1 (M (ν + 1)) ≤ Mγ ν γ (ν + 1)γ M γ ν=1 (ν + 1)ν γ ν=1
≤
≤ χm,γ A
lnm−1 M Mγ
.
The theorem is proven. We will show that if N = p is a prime then there exists sufficiently many points P corresponding to large values of M. ⃗ belonging to the set ΓMm . For 0 < θ < 1 let Mpθ be the largest Denote by SM the number of points n integer such that SM θ ≤ θ p. p
⃗ with at least one component being coprime with p, the number of solutions of For any fixed n ⃗) ≡ 0 (mod p) that satisfy the condition 0 ≤ a1 , . . . , am ≤ p − 1, equals pm−1 . the equation (⃗ n, a ⃗ |||⃗n||| < Mpθ , is not ⃗) ̸≡ 0 (mod p) for all n⃗ with n⃗ ̸= 0, Therefore the number of points P such that (⃗ n, a m less than (1 − θ )p . Consequently, such P can be obtained by randomly choosing arbitrary a1 , . . . , am ⃗) ̸≡ 0 (mod p) has solutions n⃗ ̸= 0⃗ such that |||⃗n||| < Mpθ . Obviously, and checking if the equation (⃗ n, a θ one has to check this condition 1−θ mp-times in expectation.
0
Let P =
a01 p
,...,
a0m p
be any point found that way. It is known that
SM ≤ χm M lnm−1 M for M > 1, therefore p
Mpθ ≍
(12)
lnm−1 p
and by Theorem 2
P RN0
(f ) = O A
ln(γ +1)(m−1) p
.
pγ
(13)
The estimation (13) can be made more accurate for most points P = Denote by
θ
⃗ a
a1 p
, . . . , apm .
⃗ (0 ≤ a1 , . . . , am ≤ p−1) such that (⃗n, a⃗) ̸≡ 0 (mod p) for |||⃗n||| < Mpθ , the sum of all a
⃗ We have ⃗= n ̸ 0. ′
sup RNP (f ) ≤
f ∈Eγ (A)
A
(⃗n,⃗a)≡0 (mod p)
|||⃗n|||γ
.
Hence
θ ⃗ a
max RNP (f ) ≤ A
f ∈Eγ (A)
||n⃗||>Mpθ
τn⃗ . |||⃗n|||γ
⃗) ≡ 0 (mod p) that satisfy the condition Here τn⃗ is the number of solutions of the equation (⃗ n, a 0 ≤ a1 , . . . , am ≤ p − 1.
N.S. Bakhvalov / Journal of Complexity (
)
–
9
It is clear that τn⃗ = pm−1 if at least one of the nk ’s cannot be divided by p, and τn⃗ = pm if all nk = tk p where tk are integers. Therefore
θ
max RNP (f ) ≤ pm σ ′ + pm−1 σ ′′ ,
(14)
f ∈Eγ (A)
⃗ a
where
σ′ = A
1
′
pγ |||⃗t |||
⃗t
, γ
1
σ ′′ = A
|||⃗n|||γ
||n⃗||>M θ
.
p
Obviously,
σ =O ′
A pγ
.
(15)
We have
|σ ′′ | ≤
∞ A s=Mpθ
sγ
(Ss+1 − Ss ).
Since Ss ≤ χm s lnm−1 s, we obtain in a similar way as above that ′′
|σ | ≤ A
∞
Ss+1
1
−
sγ
s=Mpθ
1
(s + 1)γ
≤ σm,γ A
lnm−1 Mpθ θ
(M p )γ −1
.
(16)
From (12), (14), (15), and (16) we get
θ ⃗ a
max RNP (f ) ≤ χm′ ,γ A
lnm−1 p
f ∈Eγ (A)
γ pm .
p
This implies that for no more than α pm points P is
χm′ ,γ A max RNP (f ) ≥ f ∈Eγ (A) α
lnm−1 p
γ
p
.
Since the summation is with respect to at least (1 − θ )pm points, then for (1 − θ − α)pm points P is
max RNP (f ) = O
A
f ∈Eγ (A)
lnm−1 p
γ
p
.
(17)
In case m = 2, using the apparatus of the theory of continued fractions one can immediately find points P such that
RNP
(f ) = O
A
Let a1 = 1, a2 = L, pµ qµ
ln N
Nγ L N
.
= [0, b1 , . . . , bn ],
= [0, b1 , . . . , bµ ],
p∗µ q∗µ
= [bµ , . . . , bn ].
10
N.S. Bakhvalov / Journal of Complexity (
)
–
From the theory of continued fractions it is known that |pµ−2 qn − pn qµ−2 | = p∗µ . Since a convergent is the best approximation of the second kind, then for integers ℓ and n2 , not simultaneously equal to 0, and |n2 | ≤ qµ−2 we have n2 NL − ℓ ≥ qµ−2 NL − pµ−2 . Hence for |n2 | ≤ qµ−2 we have |n2 L − ℓN | ≥ p∗µ , since pn = L and qn = N. Using that
(⃗n, a⃗) = n1 + n2 L ≡ 0 (mod N ) ⇐⇒ n1 = n2 L − ℓN for some ℓ ∈ Z, ⃗) ≡ 0 (mod N ) implies |n1 | ≥ p∗µ . Consequently, (⃗n, a⃗) ̸≡ we obtain that |n2 | ≤ qµ−2 and (⃗ n, a ⃗ ̸= 0⃗ with |||⃗n||| < minµ qµ−2 p∗µ = ML . 0 (mod N ) for all n We now show that ML ≥
N 2+b
,
where b = max bj . j
p∗ µ,s q∗ µ,s
:= [bµ , . . . , bs ], ts := (2 + bµ−1 )qµ−2 p∗µ,s and p∗µ,µ−1 := 1. By the well known formula qs = bs qs−1 + qs−2 we have Let
tµ−1 = (2 + bµ−1 )qµ−2 ≥ bµ−1 qµ−2 + qµ−3 = qµ−1 ; tµ = (2 + bµ−1 )qµ−2 bµ ≥ bµ qµ−1 + qµ−2 = qµ . Note that ts = bs ts−1 +ts−2 and hence ts ≥ qs for s ≥ µ+1. Therefore tn = (2+bµ−1 )qµ−2 p∗µ ≥ qn = N, which proves N ≤ ML (2 + b). For short, in all of the estimations above we omitted obvious cases µ = 2, n1 , n2 = 1. By Theorem 2,
P ln Nb R (f ) ≤ χ2,γ A bγ . N γ
(18)
N
The best quadrature formulas are apparently the formulas corresponding to L = ξk−2 , N = ξk and L = 2ξk−2 ,
N = 2ξk ,
where ξk is the successive element of the Fibonacci sequence 0, 1, 1, 2, 3, 5, 8, 13, 21, . . . whose (k + 1)st element is expressed by the previous elements according to the formula ξk+1 = ξk + ξk−1 . ⃗ = 0⃗ for For such L and N the equation n1 + n2 L ≡ 0 (mod N ) does not have solutions different from n
|||⃗n||| < ML ∼
2N 3+
√ . 5
Since E p+λ (A) ⊃ H (p, A′ , λ), by Theorem 1 the estimations (13), (17), (18) cannot be essentially m
improved, i.e., there are no quadrature formulas that in the considered case have the error of order o(N −γ ). 4. Computation of multiple integrals using the Monte Carlo method The now proposed method is a generalization of the method proposed by I. I. Shapiro-Piatetski [2] and allows computation of the integral of smooth functions with the higher accuracy. We denote by E γ (A) the class of functions satisfying the following conditions: |f | ≤ A and f periodic with respect to each variable with period 1, the absolute value of the derivative fxγ ...xγ and 1
all subordinate derivatives do not exceed A. If n1 . . . nk ̸= 0 then 1
cn1 ,...,nk ,0,...,0 =
1
f (x1 , . . . , xm ) · e
··· 0
2π ı
0
k
j=1
n j xj
1,...,1,0,...,0
dx1 . . . dxm =
dn1 ,...,nk ,0,...,0
(2π )kγ |||⃗n|||γ
,
m
N.S. Bakhvalov / Journal of Complexity (
)
–
11
1,...,1,0,...,0
where dn1 ,...,nk ,0,...,0 is the Fourier coefficient of the function fxγ ...xγ ; hence the Fourier coefficients of 1
the function f can be expressed as cn⃗ =
dn⃗
|||⃗n|||γ
,
where
k
d2n⃗ ≤ (σm′ )2 A2 .
⃗ n
Let
t
N
z
j
Nt
=
j=−N
µjN ,t z j .
j=−Nt
We put t = γ + 1 and Nt
SNP ,t (f ) =
j=−Nt
µjN ,t f (jP ) .
(2N + 1)t
(19)
Substituting
f =
2π ı
nk xk
k
cn⃗ e
in (19)
⃗ n
we obtain
SNP ,t (f ) =
⃗ n
cn⃗
sin (2N + 1)π
t
n x
k k k . (2N + 1) sin π nk x k k
Since I (f ) = c0⃗ we have
RNP ,t = SNP ,t (f ) − I (f )
t ′ dn⃗ sin (2N + 1)π k nk xk = . γ |||⃗n||| ⃗ n (2N + 1) sin π nk x k
(20)
k
Let γ , α > 0 and γ − α > 1/2. By Hölder and Bunyakovsky’s inequalities and the monotonicity of the ℓp -(quasi) norms for p > 0 we have:
P R (f ) ≤ N ,t
·
⃗ n
2 1/2 k 2 dn⃗ (2N + 1) sin π nk xk
sin (2N + 1)π
nk xk
k
⃗ n
1 γ
|||⃗n||| γ −α
γ −α γ γ −α sin (2N + 1)π nk x k k . (2N + 1) sin π n x k k k
12
N.S. Bakhvalov / Journal of Complexity (
)
–
Each of the expressions in the curly brackets has the form
1+ε sin (2N + 1)π nk x k k , S (x1 , . . . , xm ) = tn⃗ (2N + 1) sin π ⃗ n nk x k k where ε > 0 and the series n⃗ |tn⃗ | is convergent. As was shown by I. I. Shapiro-Piatetski, 1
... 0
cε
1
|tn⃗ |
⃗ n
S (x1 , . . . , xm ) dx1 . . . dxm ≤
,
N
0
and, consequently, the measure of the set of points P = (x1 , . . . , xm ) for which cε
|tn⃗ |
⃗ n
S (x1 , . . . , xm ) ≥
Nβ
is at most β . Since this applies to both expressions in the curly brackets, the measure of the set of points for which
P R (f ) > N ,t
1/2 2 c 1 dn⃗ c ⃗ n
βN
t /(t −α)
βN
⃗ n
γ −α
1
|||⃗n|||γ /(γ −α)
does not exceed 2β . From this it follows that if P is a randomly chosen point from the unit cube then with probability larger than 1 − δ is
RNP ,t (f ) = O
χ α,δ · A
N γ +1/2−α
.
Here α, δ > 0 are arbitrary small numbers. 5. On the set of points corresponding to good quadrature formulas Let f ∈ Eγ (A), γ > 1, and t = ⌈γ ⌉ be the smallest integer such that γ ≤ t and p = max(1, p). The following theorem holds. Theorem 3. Let N > 1 and α > 0. For almost all points P of the unit cube we have max RNP ,t (f ) ≤ A
f ∈Eγ (A)
χγ (P , α) Nγ
ln(m+α)(γ +1) N .
We first show a generalization of the theorem proven by A. Ya. Khinchin in [5] to the multidimensional case. Lemma 5. Let ϕ : N → (0, ∞). If the series
∞
1 s=0 sϕ(s)
is convergent, then for almost all points P =
(x1 , . . . , xm ) of the unit cube, and all n⃗ ̸= 0 we have cP nk x k ≥ m > 0. k [nk ϕ(nk )] k=1
Here ⟨⟨y⟩⟩ denotes the distance from y to the closest integer.
N.S. Bakhvalov / Journal of Complexity (
)
–
13
Proof. We prove that the number cP appearing in the lemma is almost everywhere positive. The ⃗ ̸= 0⃗ and ε < 1/2. measure of the set of points P satisfying the condition k nk xk ≤ ε equals 2ε , if n Consequently, the measure of the set of those points for which
≤
nk x k
k
η [nk ϕ(nk )]
k
⃗ ̸= 0⃗ does not exceed for any n 2η
[nk ϕ(nk )]
⃗ n
1
′
∞
1
s=−∞
sϕ(s)
= 2η
k
m
−1 .
Hence we get that the measure of the set of points P with the property that there is no η > 0 such ⃗ ̸= 0⃗ we have that for all n
>
nk x k
k
η m
[nk ϕ(nk )]
k=1
equals 0.
From now on we consider those points P for which
nk x k
cP
≥
k
m
nk log2 nk
1+ α > 0
(21)
m
k=1
⃗ From (20) we have ⃗ ̸= 0. for all n t sin (2N + 1)π nk x k ′ 1 P k . R (f ) ≤ A N ,t γ |||⃗n||| ⃗ n nk x k (2N + 1) sin π k
Let γ = 2γ /(γ − 1). We split the sum on the right hand side into two parts: ∗ ⃗ for which |||⃗n||| ≥ N γ , and ′′′ with respect to the remaining n⃗. Since n ∗
′′
N ı2ky e sin((2N + 1)y) k=−N (2N + 1) sin y = 2N + 1 ≤ 1, we have
′′ ≤A
||n⃗||≥N γ
1 ∗
|||⃗n|||γ
.
The last expression coincides with the expression contained in (16), and therefore
′′
= O A
2γ
ln N γ −1 N 2γ
m−1 −γ = O(AN ).
with respect to
14
N.S. Bakhvalov / Journal of Complexity (
)
–
For the second part note that, for 0 ≤ y ≤ 1 we have | sin π y| ≥ 2⟨⟨y⟩⟩, and hence,
′′′ ≤
A
(2N +
1
1)γ
||n⃗||
1
γ . |||⃗n|||γ nk xk 2
∗
(22)
k
Denote by which
′′′ s
cP 2s+2 sm+α
≤
⃗ for the part of the sum in the right hand side of (22) that corresponds to those n
cP
<
nk x k
2s+1 sm+α
k
In the following estimation of the sum Section 3. From (21) it follows that for |||⃗ n||| < 2s is
≥
nk xk
k
cP 2s sm+α
. ′′′ s
we essentially repeat the proof of Theorem 2 of
,
⃗ such and therefore in any hyperrectangle of the type P2ms (see Section 3) there is at most one point n that
<
nk xk
k
cP 2s+1 sm+α
.
This can be proven in the same way as Lemma 3. ⃗ such that Denoting by Tνs the number of points n
<
nk xk
k
cP
|||⃗n||| < ν 2s ,
and
2s+1 sm+α
we obtain from Lemma 4 that Tνs ≤ σm ν lnm−1 ν 2s . It is easily seen that
′′′ s≤
≤
A
(2N +
2s+1 sm+α
1)γ
cP
A
γ
cP
1
|||⃗n|||γ
cP nk xk < s+1 m+α 2 s k
2s+1 sm+α
(2N + 1)γ
γ ∞ s s Tν+ 1 − Tν . (ν 2s )γ ν=1
Using the same sequence of inequalities as in the proof of Theorem 2 we obtain
′′′ s ≤
χm′′ ,γ A
sm+α
(2N + 1)γ
γ
cP
sm−1 .
⃗ ̸= 0⃗ with |||⃗n||| < N γ we have Since for n ∗
k
nk x k
>
cP N γ (γ ∗ log2 N )m+α ∗
then
′′′ s
= 0 for s >
2γ
γ −1
log2 N .
,
N.S. Bakhvalov / Journal of Complexity (
)
–
15
Therefore
′′′ ≤ ≤
χm′′ ,γ A (2N +
1
1)γ
ξγ′ (P , α) A Nγ
γ
cP
2γ
γ −1
log2 N
(sm+α )γ sm−1
s=1
(ln N )(m+α)(γ +1) .
= O(AN −γ ) then P R (f ) ≤ χγ (P , α) A (ln N )(m+α)(γ +1) . N ,t γ
Since
′′
N
The methods described in Sections 3–5 can be directly applied only for periodic functions f (x1 , . . . , xm ). However, computation of integrals for functions satisfying all the conditions of being in the class E γ (A) with 1 < γ ≤ 2 except the periodicity condition, can be reduced to the computation of integrals
of functions from the class Eγ . As was shown by N. N. Chentsov, we have I (f ) = I (f˜ ), where
δ1 δm f (−1) x1 + δ1 , . . . , (−1) xm + δm ∈ Eγ .
f˜ = 2−m
δ1 ,...,δm =0,1
In practice, more convenient is, proposed by N. N. Chentsov and the author, the change I (f ) = I (fˆ ) where fˆ (x1 , . . . , xm ) = f
1 1 1 − 2 − x1 , . . . , 1 − 2 − xm ∈ Eγ . 2 2
Acknowledgment The author expresses his gratitude to N. N. Chentsov for his many valuable advises whilst accomplishing and writing down this work. References [1] G. Dahlquist, The Monte Carlo method, Nord. Mat. Tidskr. 2 (1954) 27–43 (in Swedish). [2] I.M. Gelfand, S.M. Feinberg, A.C. Frolov, N.N. Chentsov, On application of the method of random events (Monte Carlo method) to solving kinetic equations, in: Proc. of the 2nd Intern. Conf. OON on Using Atomic Energy for Piece Purposes (no year). [3] N.M. Günter, Potential Theory and its Applications to Basic Problems of Mathematical Physics, Gosudarstv. Izdat. Tehn.-Teor. Lit., Moscow, 1953. [4] J.M. Hammersley, J.G. Mauldon, General principles of antithetic variates, Proc. Cambridge Philos. Soc. 52 (1956) 476–481. [5] A.Ya. Khinchin, Metric problems in the theory of irrational numbers, Uspekhi Mat. Nauk (1) (1936) 7–32. [6] A.N. Kolmogorov, On certain asymptotic characteristics of completely bounded metric spaces, Dokl. Akad. Nauk SSSR (N.S.) 108 (1956) 385–388 (in Russian). [7] N.M. Korobov, Approximate calculation of repeated integrals by number-theoretical methods, Dokl. Akad. Nauk SSSR (N.S.) 115 (1957) 1062–1065 (in Russian). [8] N.M. Korobov, Computation of multiple integrals by the method of optimal coefficients, Vestnik Moskov. Univ. Ser. Mat. Meh. Astr. Fiz. Him. 1959 (1959) 19–25 (in Russian).
Further reading [1] A.Ya. Khinchin, Continued Fractions, Gosudarstv. Izdat. Tehn.-Teor. Lit., Moscow, Leningrad, 1949 (in Russian). [2] K.F. Roth, Rational approximations to algebraic numbers, Mathematika 2 (1955) 1–20. [3] P.A. Samet, The product of non-homogeneous linear forms, I & II, Proc. Cambridge Philos. Soc. 50 (1954) 372–379 & 380–390.
)
–
Contents lists available at ScienceDirect
Journal of Complexity journal homepage: www.elsevier.com/locate/jco
On the approximate calculation of multiple integrals Nikolai Sergeevich Bakhvalov Department of Computational Mathematics, Moscow State University, Russian Federation
article
info
Article history: Available online xxxx Keywords: Optimal quadrature formulas Monte Carlo methods
abstract When approximately calculating integrals of high dimension with the Monte Carlo method, one uses fewer values of the integrand than when calculating with the help of classical deterministic quadrature formulas. However, the error estimation for the Monte Carlo method does not depend on the smoothness of the integrand. This suggests that it is possible to obtain methods that give a better order of convergence in case of smooth functions. © 2015 Published by Elsevier Inc.
1. Lower estimation of the possible error We show a lower bound for the maximum error of evaluating integrals for functions of a given smoothness. Consider the following problem: suppose we have the class of all functions fj (n) (n = 1, 2, . . . , 2N) that take the values ±1, and one wants to evaluate the arithmetic mean of a function from this class 2N
Sj =
n =1
fj (n)
.
2N We have the following lemma. Lemma 1. If a deterministic method of approximately evaluating the arithmetic mean uses information about the function values at no more than N points, then there is a function from the class for which the difference between the computed arithmetic mean and the true value is at least 1/2. Indeed, suppose that the arithmetic mean was evaluated using information about values of f at points i1 , i2 , . . . , im (m ≤ N), and that f (ik ) = αk (k = 1, . . . , m). Consider functions ft (n) and fτ (n) such that ft (ik ) = fτ (ik ) = αk (k = 1, . . . , m) and ft (n) = −fτ (n) = 1 at the other points. http://dx.doi.org/10.1016/j.jco.2014.12.003 0885-064X/© 2015 Published by Elsevier Inc.
2
N.S. Bakhvalov / Journal of Complexity (
)
–
Since |St − Sτ | ≥ 1 and the calculated approximate values S t and S τ of the arithmetic means are the same, then max(|St − S t |, |Sτ − S τ |) ≥
1 2
.
Definition 1. A method of evaluating an integral (arithmetic mean) is called non-deterministic if the coordinates of each successive knot are random functions that depend on the integration domain and all values obtained previously in the computational process (coordinates of the knots, function values and its derivatives, the number of knots, etc.), and the approximate value of the integral (arithmetic mean) is a random function depending on all the obtained values. This definition applies to the Monte Carlo method, the method proposed in [4], methods described in Sections 2 and 4 of the present paper, and many others. Lemma 2. For any nondeterministic method of evaluating the arithmetic mean that uses information about function values at no more than N points, there is a function in the given class for which the average of the absolute error is at least C0 N −1/2 where C0 is an absolute constant. We will say that an event Kγ ,
α1 , . . . , αmγ γ = β1 , . . . , βmγ
,
mγ ≤ N, αk = ±1, took place in the process of evaluating the arithmetic mean, if information about values of the function at points β1 , . . . , βmγ was used, and f (βk ) = αk
(k = 1, . . . , mγ ).
Since the approximate value of the arithmetic mean (integral) can be written using a finite number of binary digits, it can take no more than countably many values Is (s = 1, 2, . . .). Denote by πγj the probability of the event Kγ in case the integral is evaluated for the function fj . γ
By σs we denote the probability that the approximate evaluation of the integral equals Is under the condition that the event Kγ took place. We associate with each function fj from the considered class the probability pj = 2−2N ; it is clear that
pj = 1.
j
Set
ρj =
γ
ρ=
πγ j
γ
σs |Is − Sj | ,
(1)
s
ρj pj .
(2)
j
Obviously, ρj is the average of the absolute value of the error when calculating the integral of the function fj . Plugging (1) into (2) we have
ρ=
γ
πγ pj j
s
j
where
qγ ,s =
j
πγj pj |Is − Sj | j
πγj pj
.
γ
σs qγ ,s ,
(3)
N.S. Bakhvalov / Journal of Complexity (
)
–
3
We will say that γ ∈ j if πγj ̸= 0. Consider those γ for which πγj ̸= 0 for at least one value of j (all πγj can be zero only when we do not use information in any of the points βk ). Obviously, the quantities πγj depend only on function values at points β1 , . . . , βmγ , and therefore, j
j
if γ ∈ j1 , γ ∈ j2 then πγ1 = πγ2 = πγ ̸= 0. Since πγj ̸= 0 only for 22N −mγ functions fj satisfying γ ∈ j, i.e., f (βk ) = αk (k = 1, . . . , mγ ), and all pj are equal to each other, then we have γ ∈j pj = 2−mγ . Hence qγ ,s = 2−(2N −mγ )
|Is − Sj |.
γ ∈j
All the functions fj satisfying the condition γ ∈ j can be divided into pairs fµk and fνk such that at points different from β1 , . . . , βγ is fµk (n) = −fνk (n). Obviously,
|Is − Sµk | + |Is − Sνk | ≥
1 γ γ |S µk | + |S νk | , 2N
where γ
Sj =
fj (n).
n̸=βk k=1,...,mγ
This implies qγ ,s ≥
t2N −mγ N
=
1
γ
N
−(2N −mγ )
2
|S j |
2
γ ∈j
.
The quantity t2N −mγ coincides with the average deviation of the binomial distribution with parameters (2N − mγ , 1/2, 1/2). By the central limit theorem, for large 2N − mγ we have
t2N −mγ =
2
π
2N − mγ + o
2N − mγ
.
Since mγ ≤ N, we conclude that for all N is t2N −mγ ≥ C0 N 1/2 . Therefore qγ ,s ≥ C0 N −1/2 . Using (3)
we obtain ρ ≥ C0 N −1/2 . As a consequence, there exists such j0 that ρj0 ≥ C0 N −1/2 , and this is what we were supposed to show. We will pass to estimating from below the maximum error of evaluating the integral I (f ) =
1
1
f (x1 , . . . , xm ) dx1 . . . dxm ,
··· 0
0
where f ∈ H (p, A, λ) and H (p, A, λ) is the space of all f ∈ C p ([0, 1]m ) such that all partial derivatives up to the order p ∈ N are bounded by A and the derivatives of order p are Hölder continuous with exponent 0 < λ ≤ 1 and constant A > 0, see [3]. The following considerations essentially coincide with the considerations that led to obtaining a lower bound for ε -entropy of the class of functions H (p, A, λ) [6]. Let N be an integer, n0 =
1 (2N ) m + 1,
p+λ) ψn0 (x1 , . . . , xm ) = Pp+λ,m n−( 0
m
[n0 xk (1 − n0 xk )]p+λ ,
k=1
where Pp+λ,m is a constant, such that ψn0 (x1 , . . . , xm ) ∈ H (p, 1, λ) for all n0 and 0 ≤ x1 , . . . , xm ≤ 1/n0 .
4
N.S. Bakhvalov / Journal of Complexity (
)
–
Let Q1 , . . . , Q2N be any 2N different cubes, nk /n0 ≤ xk ≤ (nk + 1)/n0 (k = 1, . . . , m; 0 ≤ nk < n0 ; nk an integer). Consider all the possible functions ϕj such that on each of the cubes
n1 nm ϕj = ±Aψn0 x1 − , . . . , xm − , n0
n0
and for the other points ϕj = 0. Obviously, for 0 ≤ x1 , . . . , xm ≤ 1, all the functions ϕj ∈ H (p, A, λ). There is a one-to-one correspondence between the functions ϕj and fj (n). Therefore, using Lemmas 1 and 2 we obtain the following theorem. Theorem 1. If the integral I (f ) is computed by a deterministic algorithm using information about the function values and its derivatives at no more than N points then there is a function f ∈ H (p, A, λ) for which the true value of the integral differs from the one computed with a given method by at least ′ −(p+λ)/m Cm , and if a nondeterministic method is used then there is a function f ∈ H (p, A, λ) for ,p+λ AN ′′
which the average of the absolute value of the error is larger than Cm,p+λ AN
1 − (p+λ)/m+ 2
.
Here Cm,p+λ , Cm,p+λ > 0 are some constants that depend on C0 and on the quantity ′
′′
1
1
ψ1 (x1 , . . . , xm ) dx1 . . . dxm .
··· 0
0
2. A method for the approximate evaluation of the integrals The classical Monte Carlo method relies on taking a sequence of independent random points from the domain of integration and computing the arithmetic mean of the function values at those points, see [1]. In case of smooth integrands, it seems better to take coordinates of the points depending on each other. An example of such a method for computation of integrals is the method proposed in [4]. Below we describe a method that is based on the same idea and allows, using function values at N points, to obtain the value of the integral with error of order O N
p+λ 1 + m 2
−
absolute value of the error of order O N
−
p+λ m
and with the average of the
.
Let t be the minimum integer such that 2t ≥ p + λ, 0 ≤ ξ ≤ 1, and let f ∈ H (p, A, λ) if −t /n0 ≤ x ≤ t /n0 . ξ Denote by Ω2t −1 (x) the Lagrange interpolation polynomial of degree 2t − 1 that agrees with the ξ ξ +t −1 function f at knots ± n , . . . , ± n . 0 0 It is known that ξ
Ω2t −1 (x) =
t −1
Wβ (x, ξ ) f
β +ξ n0
β=0
= f (x) + O
A p+λ
n0
+ W−β (x, ξ )f
β +ξ − n0
.
(4)
Integrating (4) in the intervals from −t /n0 to t /n0 we obtain t n0
−
t n0
f (x) dx =
t −1
β=0
πβ (ξ ) f
β +ξ n0
+f
−
β +ξ n0
+O
A 1+p+λ
n0
.
(5)
N.S. Bakhvalov / Journal of Complexity (
)
–
5
k k Let us perform a series of mnm 0 trials with outcomes ξn1 ...nm , where ξn1 ...nm are independent random quantities distributed uniformly on the interval [0, 1]. We will assume that f ∈ H (p, A, λ) on the whole of Rm and that f is periodic with respect to each variable with period 1. We set t −1
σn1 ...nm =
m
πβk ξnk1 ...nm
β1 ,...,βm =0 k=1
×
f
n1 + (−1)α1 (ξn11 ...nm + β1 ) n0
α1 ,...,αm =0,1
,...,
m nm + (−1)αm (ξn1 ...nm + βm )
.
n0
From (5) it follows that
σn1 ...nm =
n 1 +t n0
n m +t n0
···
n 1 −t n0
f dx1 . . . dxm + O
n m −t n0
A m+p+λ
n0
.
(6)
Taking the sum of (6) with respect to all nk (0 ≤ n1 , . . . , nm < n0 ) and dividing by (2t )m we obtain n 0 −1
σ =
n1 ...nm =0
σn1 ...nm
= I (f ) + O
(2t )m
A
.
p+λ
n0
(7)
From (6) we have
D σn1 ...nm = O
A2
.
2(m+p+λ)
n0
Due to independence of ξnk1 ...nm the values σn1 ...nm are also independent and therefore D(σ ) =
n 0 −1
1
(2t )2m
D σn1 ...nm = O
n1 ...nm =0
A2 m+2(p+λ)
n0
.
(8)
Now we show that M σ = I (f ). For short, we consider only m = 1. We have: M σn1 =
n1 +β+1 n0
t −1
n1 +β n0
β=0
−
t −1
β=0
− −
πβ (n0 x − n1 − β)f (x) dx
n1 −β−1 n0
n1 −β n0
πβ (|n0 x − n1 − β|) f (x) dx.
(9)
ξ
Since for f (x) = 1 is Ω2t −1 (x) = 1, then t −1
πβ (ξ ) =
β=0
t n0
.
(10)
Taking the sum of (9) with respect to n1 in the intervals from 0 to n0 − 1 and using (10) we obtain: Mσ =
n 0 −1 1
2t
n 1 =0
M σn1 = I (f ).
6
N.S. Bakhvalov / Journal of Complexity (
)
–
In this method of evaluating the integral one uses function values at N ≍ nm 0 points. Since
M (σ − I )2 ≥ M |σ − I |, from (7) and (8) we have
p+λ σ − I = O AN − m ,
M |σ − I | = O AN −
p+λ 1 m +2
.
By comparing these estimates with the results of Theorem 1 one can see that there is no method of calculating the integral whose error estimate has better order. 3. A theorem on existence of good quadrature formulas Let us consider the problem of evaluating the integral in case f ∈ Eγ (A), where Eγ (A) (γ > 1) is the class of periodic functions with period 1 with respect to each variable, such that the Fourier coefficients satisfy the condition
|⃗cn | ≤
A
|||⃗n|||γ
.
⃗ = (n1 , . . . , nm ), p = max(1, |p|), |||⃗n||| = k=1 nk , (⃗n, a⃗) = k=1 nk ak , Here and in the sequel n and σm , χm,γ are some constants. N.M. Korobov [7,8] dealt with derivation of formulas that are useful for computation of integrals for such functions. Using some other arguments he showed classes of quadrature formulas with error estimates that differ from the estimates obtained in this section by multipliers of the form lnα N. We will consider quadrature formulas of the form m
N −1
ΣNP
(f ) =
f (jP )
j =0
N
,
where P =
a
1
N
,...,
am N
m
, ak -integer.
We have [8]
RNP (f ) = ΣNP (f ) − I (f ) =
′ cn⃗ , (⃗n,⃗a)≡0 (mod N )
(11)
′
⃗ = 0⃗ is excluded. where means that in the summation n ⃗ we mean a vector with integer coordinates and we write In what follows, by n ⃗ (⃗n, P ) ≡ 0 (mod 1) = n⃗ ̸= 0: ⃗ (⃗n, a⃗) ≡ 0 (mod N ) . ⃗ ̸= 0: LP := n
We will prove the following theorem.
⃗ ∈ LP , then Theorem 2. Let M > 1. If |||⃗ n||| ≥ M for all points n m−1 P M R (f ) ≤ χm,γ · A ln N γ
M
for all f ∈ Eγ (A). For the proof of this theorem we need the following two lemmas. Lemma 3. Let m PM
:=
m
[ak , bk ] ⊂ R : ak < bk , k = 1, . . . , m, and m
k=1
m
(bk − ak ) < M
k=1
be the set of all axes-parallel hyperrectangles of volume less than M. Then, under the assumption of Theorem 2, we have
LP ∩ Ω ≤ 1 for all Ω ∈ P m . M
N.S. Bakhvalov / Journal of Complexity (
)
–
7
⃗, n⃗′ ∈ LP , then, obviously, we also have n⃗′′ = Proof. Assume that Ω contains 2 different points n ′ ′ ⃗ − n⃗ ∈ LP . By assumption, this means |||⃗n − n⃗||| ≥ M. But since Ω ∈ PMm and n⃗, n⃗′ ∈ Ω we know that n ⃗||| ≥ M: a contradiction. the volume of Ω is strictly smaller than M and larger than |||⃗ n′ − n Define the hyperbolic cross ΓMm by
ΓMm := n⃗: |||⃗n||| < M . We will show the following lemma. Lemma 4. Let M > 1 and ν ∈ N. The set ΓνmM can be covered using σm ν lnm−1 M ν hyperrectangles from m PM . Proof. We perform the proof using induction on m. It is clear that for m = 1 the assertion of the lemma holds with σm = 3. Suppose it is proven for some m. We split the set ΓνmM+1 into parts by the hyperplanes nm+1 = 0, nm+1 = ±2i ν , i = 0, 1, . . . , ⌊log2 M ⌋. m ⃗ ∈ ΓνmM+1 . Therefore the intersection of the set If 2i ν < |nm+1 | ≤ 2i+1 ν , then k=1 nk < M for n 2i
ΓνmM+1 with the hyperplane nm+1 = ν 2i + 1 can be covered using at most qi := σm lnm−1 (i)
(i)
M 2i
hyper-
m rectangles from PM , say Q1 , . . . , Qqi . We define the (m + 1)-dimensional hyperrectangles from /2i (i) m+1 i Pν M by Qℓ × [ν 2 , ν 2i+1 ], ℓ = 1, . . . , qi . It is clear that each such hyperrectangle can be covered m+1 using ν hyperrectangles of the type PM . As a consequence, all the points of the set ΓνmM+1 , such that m+1 2i ν < nm+1 ≤ 2i+1 ν , can be covered using ν qi hyperrectangles of the type PM .
If |nm+1 | ≤ ν , then, by inductive assumption, all points of the set ΓνmM+1 such that nm+1 = j ∈ N can be covered using q∗,j := σm
ν j
+ 1 lnm−1 (M ⌈ νj ⌉) (∗,j)
m hyperrectangles from PM , say Q1
ν
(∗,j)
Qℓ
j) , . . . , Qq(∗, ∗,j . Clearly, the collection of (m + 1)-dim. hyperrectangles
× [j − 1, j]: ℓ = 1, . . . , q∗,j
j =1 m+1 contains every point of ΓνmM+1 with 0 ≤ nm+1 ≤ ν , and each element is from PM .
From this it follows that the total number of hyperrectangles sufficient to cover the set ΓνmM+1 does not exceed 2
⌊log M ⌋ 2
νσ
m−1 M m ln 2i
+
ν j =1
i =0
σm
ν j
m−1
+ 1 ln
≤ 2 σm ⌊log2 M ⌋ν ln ′
m−1
M + 2νσm ln
m−1
Mν
ν
(M ⌈ j ⌉)
ν 1 j=1
j
≤ σm+1 ν lnm M ν. (The additional 2 comes from the absolute value in the restriction on nm+1 .) Thus the lemma is proven. We proceed with the proof of the theorem. We denote by Tν the number of solutions of the equation
(⃗n, a⃗) ≡ 0 (mod N ) that belong to ΓνmM . Combining Lemmas 3 and 4, we obtain Tν ≤ σm ν lnm−1 M ν.
8
N.S. Bakhvalov / Journal of Complexity (
)
–
From (11) it follows that ∞ P R (f ) ≤ | c | ≤ ⃗ n N ν=1
⃗∈LP n ∞
=
ATν+1
A
(M ν)γ
1
−
(M ν)γ
ν=1
(Tν+1 − Tν ) 1
(M (ν + 1))γ
−
AT1 Mγ
∞ ∞ ATν+1 (ν + 1)γ − ν γ A γ σm (ν + 1) lnm−1 (M (ν + 1)) ≤ Mγ ν γ (ν + 1)γ M γ ν=1 (ν + 1)ν γ ν=1
≤
≤ χm,γ A
lnm−1 M Mγ
.
The theorem is proven. We will show that if N = p is a prime then there exists sufficiently many points P corresponding to large values of M. ⃗ belonging to the set ΓMm . For 0 < θ < 1 let Mpθ be the largest Denote by SM the number of points n integer such that SM θ ≤ θ p. p
⃗ with at least one component being coprime with p, the number of solutions of For any fixed n ⃗) ≡ 0 (mod p) that satisfy the condition 0 ≤ a1 , . . . , am ≤ p − 1, equals pm−1 . the equation (⃗ n, a ⃗ |||⃗n||| < Mpθ , is not ⃗) ̸≡ 0 (mod p) for all n⃗ with n⃗ ̸= 0, Therefore the number of points P such that (⃗ n, a m less than (1 − θ )p . Consequently, such P can be obtained by randomly choosing arbitrary a1 , . . . , am ⃗) ̸≡ 0 (mod p) has solutions n⃗ ̸= 0⃗ such that |||⃗n||| < Mpθ . Obviously, and checking if the equation (⃗ n, a θ one has to check this condition 1−θ mp-times in expectation.
0
Let P =
a01 p
,...,
a0m p
be any point found that way. It is known that
SM ≤ χm M lnm−1 M for M > 1, therefore p
Mpθ ≍
(12)
lnm−1 p
and by Theorem 2
P RN0
(f ) = O A
ln(γ +1)(m−1) p
.
pγ
(13)
The estimation (13) can be made more accurate for most points P = Denote by
θ
⃗ a
a1 p
, . . . , apm .
⃗ (0 ≤ a1 , . . . , am ≤ p−1) such that (⃗n, a⃗) ̸≡ 0 (mod p) for |||⃗n||| < Mpθ , the sum of all a
⃗ We have ⃗= n ̸ 0. ′
sup RNP (f ) ≤
f ∈Eγ (A)
A
(⃗n,⃗a)≡0 (mod p)
|||⃗n|||γ
.
Hence
θ ⃗ a
max RNP (f ) ≤ A
f ∈Eγ (A)
||n⃗||>Mpθ
τn⃗ . |||⃗n|||γ
⃗) ≡ 0 (mod p) that satisfy the condition Here τn⃗ is the number of solutions of the equation (⃗ n, a 0 ≤ a1 , . . . , am ≤ p − 1.
N.S. Bakhvalov / Journal of Complexity (
)
–
9
It is clear that τn⃗ = pm−1 if at least one of the nk ’s cannot be divided by p, and τn⃗ = pm if all nk = tk p where tk are integers. Therefore
θ
max RNP (f ) ≤ pm σ ′ + pm−1 σ ′′ ,
(14)
f ∈Eγ (A)
⃗ a
where
σ′ = A
1
′
pγ |||⃗t |||
⃗t
, γ
1
σ ′′ = A
|||⃗n|||γ
||n⃗||>M θ
.
p
Obviously,
σ =O ′
A pγ
.
(15)
We have
|σ ′′ | ≤
∞ A s=Mpθ
sγ
(Ss+1 − Ss ).
Since Ss ≤ χm s lnm−1 s, we obtain in a similar way as above that ′′
|σ | ≤ A
∞
Ss+1
1
−
sγ
s=Mpθ
1
(s + 1)γ
≤ σm,γ A
lnm−1 Mpθ θ
(M p )γ −1
.
(16)
From (12), (14), (15), and (16) we get
θ ⃗ a
max RNP (f ) ≤ χm′ ,γ A
lnm−1 p
f ∈Eγ (A)
γ pm .
p
This implies that for no more than α pm points P is
χm′ ,γ A max RNP (f ) ≥ f ∈Eγ (A) α
lnm−1 p
γ
p
.
Since the summation is with respect to at least (1 − θ )pm points, then for (1 − θ − α)pm points P is
max RNP (f ) = O
A
f ∈Eγ (A)
lnm−1 p
γ
p
.
(17)
In case m = 2, using the apparatus of the theory of continued fractions one can immediately find points P such that
RNP
(f ) = O
A
Let a1 = 1, a2 = L, pµ qµ
ln N
Nγ L N
.
= [0, b1 , . . . , bn ],
= [0, b1 , . . . , bµ ],
p∗µ q∗µ
= [bµ , . . . , bn ].
10
N.S. Bakhvalov / Journal of Complexity (
)
–
From the theory of continued fractions it is known that |pµ−2 qn − pn qµ−2 | = p∗µ . Since a convergent is the best approximation of the second kind, then for integers ℓ and n2 , not simultaneously equal to 0, and |n2 | ≤ qµ−2 we have n2 NL − ℓ ≥ qµ−2 NL − pµ−2 . Hence for |n2 | ≤ qµ−2 we have |n2 L − ℓN | ≥ p∗µ , since pn = L and qn = N. Using that
(⃗n, a⃗) = n1 + n2 L ≡ 0 (mod N ) ⇐⇒ n1 = n2 L − ℓN for some ℓ ∈ Z, ⃗) ≡ 0 (mod N ) implies |n1 | ≥ p∗µ . Consequently, (⃗n, a⃗) ̸≡ we obtain that |n2 | ≤ qµ−2 and (⃗ n, a ⃗ ̸= 0⃗ with |||⃗n||| < minµ qµ−2 p∗µ = ML . 0 (mod N ) for all n We now show that ML ≥
N 2+b
,
where b = max bj . j
p∗ µ,s q∗ µ,s
:= [bµ , . . . , bs ], ts := (2 + bµ−1 )qµ−2 p∗µ,s and p∗µ,µ−1 := 1. By the well known formula qs = bs qs−1 + qs−2 we have Let
tµ−1 = (2 + bµ−1 )qµ−2 ≥ bµ−1 qµ−2 + qµ−3 = qµ−1 ; tµ = (2 + bµ−1 )qµ−2 bµ ≥ bµ qµ−1 + qµ−2 = qµ . Note that ts = bs ts−1 +ts−2 and hence ts ≥ qs for s ≥ µ+1. Therefore tn = (2+bµ−1 )qµ−2 p∗µ ≥ qn = N, which proves N ≤ ML (2 + b). For short, in all of the estimations above we omitted obvious cases µ = 2, n1 , n2 = 1. By Theorem 2,
P ln Nb R (f ) ≤ χ2,γ A bγ . N γ
(18)
N
The best quadrature formulas are apparently the formulas corresponding to L = ξk−2 , N = ξk and L = 2ξk−2 ,
N = 2ξk ,
where ξk is the successive element of the Fibonacci sequence 0, 1, 1, 2, 3, 5, 8, 13, 21, . . . whose (k + 1)st element is expressed by the previous elements according to the formula ξk+1 = ξk + ξk−1 . ⃗ = 0⃗ for For such L and N the equation n1 + n2 L ≡ 0 (mod N ) does not have solutions different from n
|||⃗n||| < ML ∼
2N 3+
√ . 5
Since E p+λ (A) ⊃ H (p, A′ , λ), by Theorem 1 the estimations (13), (17), (18) cannot be essentially m
improved, i.e., there are no quadrature formulas that in the considered case have the error of order o(N −γ ). 4. Computation of multiple integrals using the Monte Carlo method The now proposed method is a generalization of the method proposed by I. I. Shapiro-Piatetski [2] and allows computation of the integral of smooth functions with the higher accuracy. We denote by E γ (A) the class of functions satisfying the following conditions: |f | ≤ A and f periodic with respect to each variable with period 1, the absolute value of the derivative fxγ ...xγ and 1
all subordinate derivatives do not exceed A. If n1 . . . nk ̸= 0 then 1
cn1 ,...,nk ,0,...,0 =
1
f (x1 , . . . , xm ) · e
··· 0
2π ı
0
k
j=1
n j xj
1,...,1,0,...,0
dx1 . . . dxm =
dn1 ,...,nk ,0,...,0
(2π )kγ |||⃗n|||γ
,
m
N.S. Bakhvalov / Journal of Complexity (
)
–
11
1,...,1,0,...,0
where dn1 ,...,nk ,0,...,0 is the Fourier coefficient of the function fxγ ...xγ ; hence the Fourier coefficients of 1
the function f can be expressed as cn⃗ =
dn⃗
|||⃗n|||γ
,
where
k
d2n⃗ ≤ (σm′ )2 A2 .
⃗ n
Let
t
N
z
j
Nt
=
j=−N
µjN ,t z j .
j=−Nt
We put t = γ + 1 and Nt
SNP ,t (f ) =
j=−Nt
µjN ,t f (jP ) .
(2N + 1)t
(19)
Substituting
f =
2π ı
nk xk
k
cn⃗ e
in (19)
⃗ n
we obtain
SNP ,t (f ) =
⃗ n
cn⃗
sin (2N + 1)π
t
n x
k k k . (2N + 1) sin π nk x k k
Since I (f ) = c0⃗ we have
RNP ,t = SNP ,t (f ) − I (f )
t ′ dn⃗ sin (2N + 1)π k nk xk = . γ |||⃗n||| ⃗ n (2N + 1) sin π nk x k
(20)
k
Let γ , α > 0 and γ − α > 1/2. By Hölder and Bunyakovsky’s inequalities and the monotonicity of the ℓp -(quasi) norms for p > 0 we have:
P R (f ) ≤ N ,t
·
⃗ n
2 1/2 k 2 dn⃗ (2N + 1) sin π nk xk
sin (2N + 1)π
nk xk
k
⃗ n
1 γ
|||⃗n||| γ −α
γ −α γ γ −α sin (2N + 1)π nk x k k . (2N + 1) sin π n x k k k
12
N.S. Bakhvalov / Journal of Complexity (
)
–
Each of the expressions in the curly brackets has the form
1+ε sin (2N + 1)π nk x k k , S (x1 , . . . , xm ) = tn⃗ (2N + 1) sin π ⃗ n nk x k k where ε > 0 and the series n⃗ |tn⃗ | is convergent. As was shown by I. I. Shapiro-Piatetski, 1
... 0
cε
1
|tn⃗ |
⃗ n
S (x1 , . . . , xm ) dx1 . . . dxm ≤
,
N
0
and, consequently, the measure of the set of points P = (x1 , . . . , xm ) for which cε
|tn⃗ |
⃗ n
S (x1 , . . . , xm ) ≥
Nβ
is at most β . Since this applies to both expressions in the curly brackets, the measure of the set of points for which
P R (f ) > N ,t
1/2 2 c 1 dn⃗ c ⃗ n
βN
t /(t −α)
βN
⃗ n
γ −α
1
|||⃗n|||γ /(γ −α)
does not exceed 2β . From this it follows that if P is a randomly chosen point from the unit cube then with probability larger than 1 − δ is
RNP ,t (f ) = O
χ α,δ · A
N γ +1/2−α
.
Here α, δ > 0 are arbitrary small numbers. 5. On the set of points corresponding to good quadrature formulas Let f ∈ Eγ (A), γ > 1, and t = ⌈γ ⌉ be the smallest integer such that γ ≤ t and p = max(1, p). The following theorem holds. Theorem 3. Let N > 1 and α > 0. For almost all points P of the unit cube we have max RNP ,t (f ) ≤ A
f ∈Eγ (A)
χγ (P , α) Nγ
ln(m+α)(γ +1) N .
We first show a generalization of the theorem proven by A. Ya. Khinchin in [5] to the multidimensional case. Lemma 5. Let ϕ : N → (0, ∞). If the series
∞
1 s=0 sϕ(s)
is convergent, then for almost all points P =
(x1 , . . . , xm ) of the unit cube, and all n⃗ ̸= 0 we have cP nk x k ≥ m > 0. k [nk ϕ(nk )] k=1
Here ⟨⟨y⟩⟩ denotes the distance from y to the closest integer.
N.S. Bakhvalov / Journal of Complexity (
)
–
13
Proof. We prove that the number cP appearing in the lemma is almost everywhere positive. The ⃗ ̸= 0⃗ and ε < 1/2. measure of the set of points P satisfying the condition k nk xk ≤ ε equals 2ε , if n Consequently, the measure of the set of those points for which
≤
nk x k
k
η [nk ϕ(nk )]
k
⃗ ̸= 0⃗ does not exceed for any n 2η
[nk ϕ(nk )]
⃗ n
1
′
∞
1
s=−∞
sϕ(s)
= 2η
k
m
−1 .
Hence we get that the measure of the set of points P with the property that there is no η > 0 such ⃗ ̸= 0⃗ we have that for all n
>
nk x k
k
η m
[nk ϕ(nk )]
k=1
equals 0.
From now on we consider those points P for which
nk x k
cP
≥
k
m
nk log2 nk
1+ α > 0
(21)
m
k=1
⃗ From (20) we have ⃗ ̸= 0. for all n t sin (2N + 1)π nk x k ′ 1 P k . R (f ) ≤ A N ,t γ |||⃗n||| ⃗ n nk x k (2N + 1) sin π k
Let γ = 2γ /(γ − 1). We split the sum on the right hand side into two parts: ∗ ⃗ for which |||⃗n||| ≥ N γ , and ′′′ with respect to the remaining n⃗. Since n ∗
′′
N ı2ky e sin((2N + 1)y) k=−N (2N + 1) sin y = 2N + 1 ≤ 1, we have
′′ ≤A
||n⃗||≥N γ
1 ∗
|||⃗n|||γ
.
The last expression coincides with the expression contained in (16), and therefore
′′
= O A
2γ
ln N γ −1 N 2γ
m−1 −γ = O(AN ).
with respect to
14
N.S. Bakhvalov / Journal of Complexity (
)
–
For the second part note that, for 0 ≤ y ≤ 1 we have | sin π y| ≥ 2⟨⟨y⟩⟩, and hence,
′′′ ≤
A
(2N +
1
1)γ
||n⃗||
1
γ . |||⃗n|||γ nk xk 2
∗
(22)
k
Denote by which
′′′ s
cP 2s+2 sm+α
≤
⃗ for the part of the sum in the right hand side of (22) that corresponds to those n
cP
<
nk x k
2s+1 sm+α
k
In the following estimation of the sum Section 3. From (21) it follows that for |||⃗ n||| < 2s is
≥
nk xk
k
cP 2s sm+α
. ′′′ s
we essentially repeat the proof of Theorem 2 of
,
⃗ such and therefore in any hyperrectangle of the type P2ms (see Section 3) there is at most one point n that
<
nk xk
k
cP 2s+1 sm+α
.
This can be proven in the same way as Lemma 3. ⃗ such that Denoting by Tνs the number of points n
<
nk xk
k
cP
|||⃗n||| < ν 2s ,
and
2s+1 sm+α
we obtain from Lemma 4 that Tνs ≤ σm ν lnm−1 ν 2s . It is easily seen that
′′′ s≤
≤
A
(2N +
2s+1 sm+α
1)γ
cP
A
γ
cP
1
|||⃗n|||γ
cP nk xk < s+1 m+α 2 s k
2s+1 sm+α
(2N + 1)γ
γ ∞ s s Tν+ 1 − Tν . (ν 2s )γ ν=1
Using the same sequence of inequalities as in the proof of Theorem 2 we obtain
′′′ s ≤
χm′′ ,γ A
sm+α
(2N + 1)γ
γ
cP
sm−1 .
⃗ ̸= 0⃗ with |||⃗n||| < N γ we have Since for n ∗
k
nk x k
>
cP N γ (γ ∗ log2 N )m+α ∗
then
′′′ s
= 0 for s >
2γ
γ −1
log2 N .
,
N.S. Bakhvalov / Journal of Complexity (
)
–
15
Therefore
′′′ ≤ ≤
χm′′ ,γ A (2N +
1
1)γ
ξγ′ (P , α) A Nγ
γ
cP
2γ
γ −1
log2 N
(sm+α )γ sm−1
s=1
(ln N )(m+α)(γ +1) .
= O(AN −γ ) then P R (f ) ≤ χγ (P , α) A (ln N )(m+α)(γ +1) . N ,t γ
Since
′′
N
The methods described in Sections 3–5 can be directly applied only for periodic functions f (x1 , . . . , xm ). However, computation of integrals for functions satisfying all the conditions of being in the class E γ (A) with 1 < γ ≤ 2 except the periodicity condition, can be reduced to the computation of integrals
of functions from the class Eγ . As was shown by N. N. Chentsov, we have I (f ) = I (f˜ ), where
δ1 δm f (−1) x1 + δ1 , . . . , (−1) xm + δm ∈ Eγ .
f˜ = 2−m
δ1 ,...,δm =0,1
In practice, more convenient is, proposed by N. N. Chentsov and the author, the change I (f ) = I (fˆ ) where fˆ (x1 , . . . , xm ) = f
1 1 1 − 2 − x1 , . . . , 1 − 2 − xm ∈ Eγ . 2 2
Acknowledgment The author expresses his gratitude to N. N. Chentsov for his many valuable advises whilst accomplishing and writing down this work. References [1] G. Dahlquist, The Monte Carlo method, Nord. Mat. Tidskr. 2 (1954) 27–43 (in Swedish). [2] I.M. Gelfand, S.M. Feinberg, A.C. Frolov, N.N. Chentsov, On application of the method of random events (Monte Carlo method) to solving kinetic equations, in: Proc. of the 2nd Intern. Conf. OON on Using Atomic Energy for Piece Purposes (no year). [3] N.M. Günter, Potential Theory and its Applications to Basic Problems of Mathematical Physics, Gosudarstv. Izdat. Tehn.-Teor. Lit., Moscow, 1953. [4] J.M. Hammersley, J.G. Mauldon, General principles of antithetic variates, Proc. Cambridge Philos. Soc. 52 (1956) 476–481. [5] A.Ya. Khinchin, Metric problems in the theory of irrational numbers, Uspekhi Mat. Nauk (1) (1936) 7–32. [6] A.N. Kolmogorov, On certain asymptotic characteristics of completely bounded metric spaces, Dokl. Akad. Nauk SSSR (N.S.) 108 (1956) 385–388 (in Russian). [7] N.M. Korobov, Approximate calculation of repeated integrals by number-theoretical methods, Dokl. Akad. Nauk SSSR (N.S.) 115 (1957) 1062–1065 (in Russian). [8] N.M. Korobov, Computation of multiple integrals by the method of optimal coefficients, Vestnik Moskov. Univ. Ser. Mat. Meh. Astr. Fiz. Him. 1959 (1959) 19–25 (in Russian).
Further reading [1] A.Ya. Khinchin, Continued Fractions, Gosudarstv. Izdat. Tehn.-Teor. Lit., Moscow, Leningrad, 1949 (in Russian). [2] K.F. Roth, Rational approximations to algebraic numbers, Mathematika 2 (1955) 1–20. [3] P.A. Samet, The product of non-homogeneous linear forms, I & II, Proc. Cambridge Philos. Soc. 50 (1954) 372–379 & 380–390.