Equidistribution on the Sphere 1 Introduction - CiteSeerX
Equidistribution on the Sphere
Jianjun Cui and Willi Freeden Laboratory of Technomathematics, Geomathematics Group University of Kaiserslautern, 67653 Kaiserslautern Email: [email protected] Key Words: sphere, pointsets, pseudodi erential operators, generalized discrepancy,
equidistribution, low discrepancy (quasi-Monte-Carlo) method, approximate integration. AMS subject classi cation: 33C55, 40A10, 47G30, 62E25, 86-08
Abstract A concept of generalized discrepancy, which involves pseudodi erential operators to give a criterion of equidistributed pointsets, is developed on the sphere. A simply structured formula in terms of elementary functions is established for the computation of the generalized discrepancy. With the help of this formula ve kinds of point systems on the sphere, namely lattices in polar coordinates, transformed 2-dimensional sequences, rotations on the sphere, triangulation, and "sum of three squares sequence", are investigated. Quantitative tests are done, and the results are compared with each other. Our calculations exhibit di erent orders of convergence of the generalized discrepancy for di erent types of point systems.
1 Introduction Of practical importance is the problem of generating equidistributed pointsets on the sphere. For that reason a concept of generalized discrepancy, which involves pseudodifferential operators to give a quantifying criterion of equidistributed pointsets, is of great interest. In this paper an explicit formula in terms of elementary functions is developed for the generalized discrepancy. Essential tools are Sobolev space structures and pseudodi erential operator techniques. It is mentioned that an optimal pointset may be obtained by minimizing the generalized discrepancy. But, in spite of the elementary representation of the generalized discrepancy, this is a non-linear optimization problem which will not be discussed here. Our investigations, however, show that there are many promising ways to generate point systems on the sphere such that the discrepancy becomes small. To be speci c, we distinguish ve kinds of point systems on the sphere: lattices in polar coordinates, transformed 2-dimensional sequences, rotations on the sphere, triangulations, and "sum of three squares sequence". By use of our developed formulas, the ve classes of point systems are described, and their discrepancies are 1
explicitly calculated for increasing numbers of points. The results show di erent orders of convergence indicated by the generalized discrepancy. Furthermore, our computations enable us to give a quantitative comparison between the di erent point systems. It is somehow surprising that certain types of transformed sequences yield the best results. Nevertheless, there are special other pointsets which provide us with better results for comparable numbers of points. For instance, the soccer ball (C60) leads us to the best result in all our considered pointsets of about 60 points. The problem of generating a large number of "equidistributed points" on the sphere has many applications in various elds of computation, particularly in geoscience, medicine. The advantage of equidistributed point systems lies in the fact that relatively few sampling of the data is needed, and approximate integration can be simply performed by computation of a mean value, i.e., the arithmetical mean. The outline of the paper is as follows: Section 2 begins with a discussion of Sobolev spaces and (invariant) pseudodi erential operators on the sphere. The generalized discrepancy is developed in Section 3. The announced point systems are studied in Section 4. The low discrepancy method, also called quasi-Monte-Carlo method, for numerical integration on the sphere is discussed in Section 5 for three types of test examples.
2 Sobolev Spaces and Pseudodi erential Operators on the Sphere Let IR3 denote the three-dimensional Euclidean space. We use x; y; z : : : to represent the elements of IR3 and the Greek alphabet ; : : : to represent the vectors of the unit sphere
in IR3. denotes the Beltrami operator on the unit sphere. A function F : 7! IR possessing k continuous derivatives on is said to be of class C k ( ). C ( ) = C 0( ) is the class of real continuous scalar-valued functions on . By L2( ) we denote the space of Lebesgue square-integrable scalar functions on . The spherical harmonics of degree n are de ned as the everywhere on in nitely di erentiable eigenfunctions of the Beltrami-operator corresponding to the eigenvalues ()n = ?n(n +1); n = 0; 1; : : : ;. As it is well-known (cf. e.g. [14]), the linear space Harmn of all spherical harmonics of degree n is of dimension 2n+1. We denote fYn;j (); n = 0; 1; : : : ; j = 1; 2; : : : ; 2n +1)g to be an orthonormalized basis of L2( ), where n is called degree and j order of the spherical harmonics. The space Harm0;:::;m of all spherical harmonics of degree m possesses the dimension (m + 1)2. The Legendre polynomials Pn are the in nitely di erentiable eigenfunctions of the Legendre operator, satisfying the normalization condition Pn (1) = 1. It should be noted that jPn (t)j 1 and jPn0 (t)j n(n2+1) for all t 2 [?1; +1] and all 2
n 2 IN0. The well-known addition theorem of spherical harmonics states (cf. [5], [14]) 2X n+1 1 P ( ) ; (; ) 2 : Yn;j ()Yn;j () = 2n4+ (1) n j =1 In particular, we have s 2X n+1 2 n + 1 jYn;j ()j2 = 4 ; jYn;j ()j 2n4+ 1 ; 2 : j =1
(2)
For s 2 IR consider the space 1 2X n+1 X 2 n^ 2s < 1g; Fn;j
E s ( ) = fF 2 C 1( )j
n=0 j =1
where Fn;j are the spherical harmonic coecients of F and by de nition, 8 < n^ = : 1 if n = 0 n elsewhere: On E s ( ) we are able to impose an inner product as follows: n+1 1 2X X Fn;j Gn;j n^2s: < F; G >H s = n=0 j =1
(3)
(4)
(5)
The Sobolev space H s( ) is the completion of E s( ) with respect to the topology kkH s = (; )1H=s2. The following relationships between Sobolev spaces are valid: H s1 ( ) H s2 ( ), whenever s1 > s2. Sobolev spaces and classical function spaces C k ( ) are related via imbedding theorems.
Lemma 2.1 (Imbedding Theorem). H s ( ) is a subspace of C k ( ), if s > k + 1. Proof. It is sucient to prove Lemma 2.1 for the case k = 0. For s > 1, applying the Cauchy-Schwarz inequality we obtain n+1 1 2X n+1 1 2X 1 2n+1 X X 2 n^ 2s ) 21 ( X X (Yn;j ( ))2 n^ ?2s ) 21 : Fn;j jFn;j Yn;j ()j ( n=0 j =1
n=0 j =1
n=0 j =1
From the addition theorem of the spherical harmonics it follows that 1 1 2X n+1 X X 1 < 1; (Yn;j ())2 n^?2s = 24nn+ C= 2 s ^ n=0 n=0 j =1 provided that s > 1. Therefore the series is absolutely convergent. But this means that the F corresponds to a continuous function on . 3
From Lemma 2.1 (Imbedding Theorem) we are able to conclude that the spherical harmonic expansion of any function in H s ( ) converges uniformly provided that s > 1. This is of signi cance, since there are functions in C ( ) which do not allow a uniformly convergent spherical harmonic expansion (cf. [7]). By aid of the Sobolev space structure introduced above we make some comments on (invariant) pseudodi erential operators on the sphere, since these operators do not belong to the traditional equipment. De nition 2.1 (Pseudodi erential Operator on the Sphere). Let fAng, n=0,1,: : : , be a sequence of real numbers An satisfying jAnj nlim !1 (n + 12 )t = const 6= 0 for some t 2 IR. Then the operator A : H s( ) 7! H s?t ( ) de ned by 1 2X n+1 X AnFn;j Yn;j ; F 2 H s ( ) (6) AF = n=0 j =1
is called pseudodi erential operator of order t. fAn g is called spherical symbol of A. Moreover, if jAnj = 0 lim n!1 (n + 1 )t 2 s for all t 2 IR, then the operator A:H ( ) 7! C 1 ( ) is called pseudodi erential operator of order -1 (note that the equation (6) is understood in H s?t ( )-topology).
The spherical symbol has many appealing properties. For example, it is readily seen that (A0 + A00)n = A0n + A00n ; (A0A00)n = A0nA00n for every choice of pseudodi erential operators A0; A00. Obviously we have
AYn;j = AnYn;j ; n = 0; 1; : : : ; j = 1; : : : ; n + 1 (i.e. all pseudodi erential operators considered here are invariant operators). Using the kernel 1 n+1 1 2X X X 1 A P ( ); AnYn;j ()Yn;j () = 2n4+ KA ( ) = n n n=0
n=0 j =1
AF is obtainable by convolution as follows Z AF = KA F = F ()KA()d!() : The kernel satis es the property KA 2 H ?(t+) ( ) for all > 0 provided that t > ?1. In the case t = ?1 we have KA 2 C 1( ). In the sequel we list some particularly important pseudodi erential operators: 4
(1) The Beltrami operator is not invertible, since ()n = 0 for n = 0, but the operator ? + 41 has the symbol f(n + 21 )2g; n = 0; 1; : : : ; and, hence, has an inverse (? + 41 )?1 which is a pseudodi erential of order ?2. More generally, (? + 14 )s is a pseudodi erential operator of order 2s and has the spherical symbol f(n + 21 )2sg; n = 0; 1; : : :. Using the kernel K(+ 14 )?2s ; s > 1, we obtain Z F () = ( + 14 )s K(+ 14 )?2s ( )( + 41 )sF ()d!()
for all F 2 H s ( ) and 2 . This is easily seen by the Parseval identity for the system of L2( )-orthonormal spherical harmonics. (2) The integral operator of the single layer potential on given by Z AU () = 21 jU?() j d!() ; 2 ; (7) is a pseudodi erential operator of order -1 possessing the symbol fAng with An = (n + 21 )?1. This can be seen by inserting the expansion of j ? j?1 in terms of Legendre polynomials. A consequence is that A is invertible. Other types of geodetically relevant integral operators can be found in [17]. (3) The operator A = (?2)(? + 14 ) 21 has the symbol fAng, where An is given by An = (2n + 1)n(n + 1); n = 0; 1; : : :. A is not invertible, but the related operator B with symbol fBng given by 8 < Bn = : 1; n = 0 (2n + 1)(n(n + 1)); n = 1; 2; : : : is invertible. Since (cf. e.g. [12])
0 s 1 1 P (t) = 1 ? 2 ln @1 + 1 ? t A ; t 2 [?1; +1] n 2 n=1 n(n + 1) 1 X
we nd for all ; 2 .
1 0 s 1 1 ? 1 A KB?1 ( ) = 2 ? 2 ln @1 + 2
(4) An = rn; 0 r < 1 (or An = e 21 n(n+1)h; h > 0) de ne pseudodi erential operators of order ?1. 5
To show the relationship between the pseudodi erential operators and the Sobolev spaces we give an equivalent reformulation of the Sobolev space H s ( ). Let A be a pseudodi erential operator of order s. Then, the Sobolev space H s ( ) can be equivalently expressed in the following form:
H s( ) = fF jF : 7! IR; AF 2 L2( )g:
(8)
3 Equidistribution on the Sphere In this section we discuss the problem of equidistribution. This problem can be stated as follows: nd a pointset f1; : : : ; N g such that the remainder term Z N X j G()d!() ? 4N G(i)j ; (9) i=1 for any function G 2 H s ( ) converges to 0 as N ! 1. Let us begin the discussion with an a priori estimate of the remainder. The result is an analogue to the Koksma-Hlawka inequality in Euclidean spaces (see also, [8,9,10] and many other references therein). Theorem 3.1 . Let A be a pseudodi erential operator of order s, s > 1, with symbol fAng satisfying An 6= 0 n 1. Then, for any function AG 2 L2( ), we have the estimate "X # 12 Z 1 2n + 1 N X N X N X 1 1 1 j 4 G()d!() ? N G(i)j N 2 Pn (i t ) kAGkL2 : (10) t=1 i=1 n=1 4An i=1 Proof. From (8) we know that AG 2 L2( ) is equivalent to G 2 H s ( ). Because of the assumption s > 1 it follows from Lemma 2.1 (Imbedding Theorem) that the spherical harmonic expansion 1 2X n+1 X Gn;j Yn;j (); 2 ; G() = n=0 j =1
converges absolutely and uniformly. We rewrite the spherical harmonic expansion of G as follows: Z n+1 Z AG( )Y ( ) 1 2X X 1 n;j G() = 4 G()d!() + d!()Yn;j (): An
n=1 j =1
By setting = i and taking the sum over all indices i; 1 i N , we obtain " X # Z N N n+1 Z AG( )Y ( ) 1 2X X 1 1 1X n;j d!() N Yn;j (i) : N i=1 G(i) = 4 G()d!() + n=1 j=1
An i=1 6
Hence, the integral error can then be estimated as follows Z N X j 41 G()d!() ? N1 G(i )j 1 2n+1 N i=1 Z 1 X X X Yn;j (i) A G ( ) Y N n;j ( )d! ( ) An
nZ=1 j=1 i=11 2n+1 N X X X Yn;j (i)Yn;j () = N1 AG() d! ( ) : An
n=1 j =1 i=1 By using the Cauchy-Schwarz inequality and the addition theorem of the spherical harmonics we get Z N X j 41 G()d!() ? N1 G(i )j i=1 2 0 3 21 12 1 Z Z N 1 2 n +1 X X X Yn;j (i)Yn;j () A N1 (AG())2d!() 2 64 @ d!()75 A n n=1 j =1 i=1 2 1 2n+1 PN !23 12 X X Y ( ) 1 n;j i i=1 5 = N kAGkL2 4 A n n=1 j =1 3 21 2 N N 1 2n+1 X X X X Y ( ) Y ( ) 1 n;j i n;j t 5 = N kAGkL2 4 2 A n i=1 t=1 n=1 j =1 "X # 21 1 2n + 1 N X N X 1 = N kAGkL2 4A2 Pn(i t) : i=1 t=1 n=1
n
This completes the proof. Theorem 3.1 extends a result given for the Beltrami operator (cf. [5]) to the general context of pseudodi erential operators. It is obvious that the right-hand side of the estimate (10) consists of two parts. The rst part depends only on the pointset, while the second part depends merely on the function G under consideration. This gives rise to the de nition of generalized discrepancy.
De nition 3.1 (Generalized discrepancy). Let A be a pseudodi erential operator of order s, s > 1, with symbol fAng, An = 6 0 for n 1. Then the generalized discrepancy associated to a pseudodi erential operator A is de ned by
"X # 21 N X 1 2n + 1 N X 1 (11) D(f1 ; : : :; N g; A) = N 2 Pn (i t) : i=1 t=1 n=1 4An The generalized discrepancy characterizes "how well a pointset is equidistributed". This can be explained by a heuristic geometrical interpretation as follows: for any subset 7
B of and any pointset f1; : : :; N g on we de ne the counting function, N X #B = B (i) ; i=1
(12)
where B is the characteristic function of B . Thus #B indicates the number n, 1 n N , of points i contained in B . Let B 2 H s, s > 1 be a "smooth" approximation of the characteristic function B. Then, for any non-empty Lebesgue-measurable subset B of
, we have an integral inequality of the following form Z Z N X j #NB ? 41 Bd!()j = j N1 B (i) ? 41 Bd!()j (13) i=1 Z N X j N1 B (i) ? 41 Bd!()j D(f1; : : : ; N g; A)kA B kL2 : i=1 Therefore, the generalized discrepancy gives a quantitative description of how well a pointset is equidistributed. Remark. It is also worth to translate the Wozniakowski idea of averaging the integration error (cf. [13]) into the spherical context. But this is a challenge for future work. For di erent pseudodi erential operators we have di erent generalized discrepancies. It is of interest to understand the relationship between di erent generalized discrepancies associated to di erent pseudodi erential operators. Lemma 3.1 . Let A, B be two pseudodi erential operators of order s1, s2, (s1 > 1, s2 > 1), and with symbols fAng, fBng satisfying An 6= 0 ; Bn 6= 0 for n 1, respectively. Then, the following relations are valid for the generalized discrepancies associated to the pseudodi erential operators A and B: 1) If the two pseudodi erential operators A and B are of the same order, i.e. s1 = s2, then there exist two positive constants C1 and C2 such that
C1D(f1; : : : ; N g; A) D(f1; : : : ; N g; B) C2D(f1; : : : ; N g; A) :
(14)
2) If s1 > s2, then there exists a constant C3 such that
D(f1; : : : ; N g; A) < C3D(f1 ; : : :; N g; B) :
Proof. Applying the addition theorem of spherical harmonics we have
n+1 Y ( )Y ( ) 1 2X 1 2n + 1 N X N X N X N X X X n;j i n;j t P ( ) = n i t 2 A2n i=1 t=1 n=1 j =1 i=1 t=1 n=1 4An n+1 PN Y ( ) !2 1 2X X i=1 n;j i : = An n=1 j =1
8
(15)
Thus the generalized discrepancy can be rewritten as follows: 2 1 2n+1 PN !23 12 X X Y ( ) i=1 n;j i 5 : D(f1; : : : ; N g; A) = 4 An n=1 j =1 If s1 = s2, then it follows from our de nition of pseudodi erential operators that there exist two constants C1 and C2 such that 1 2 1 2 1 2 C1 A B C2 A ; n 1: n n n Therefore it is clear that C1D(f1 ; : : :; N g; A) D(f1 ; : : :; N g; B) C2D(f1 ; : : :; N g; A) : The part (2) can be proved in the same manner. From Lemma 3.1 we deduce that the generalized discrepancy is dependent only on the order of the pseudodi erential operators. Using the generalized discrepancy we can study the property of point systems in quantitative way. De nition 3.2 (Equidistribution in H s ( )). A point system f1N ; : : :; NN g; N = 1; 2; : : : is called A-equidistributed in H s ( ), s > 1, if the generalized discrepancy associated to a pseudodi erential operator A of order s, s > 1, satis es lim D(f1N ; : : :; NN g; A) = 0 : (16) It should be noted that if a point system is equidistributed in H s0 ( ), s0 > 1, then the point system is equidistributed in H s( ), s > s0. This shows us that we should use s0, s0 > 1, as small as possible for our quantitative investigations. However, in order to compute the generalized discrepancy as illustrated above, one has to evaluate a series expansion in terms of Legendre polynomials. This summation is very time consuming, in particular, for relatively small s. Therefore it is not advisable to use this series expansion for the computation of the generalized discrepancy. Fortunately, for certain pseudodi erential operators, the series is representable in terms of elementary functions. One example should be discussed in more detail. Consider theqpseudodi erential operator D = (?2) 21 (? + 14 ) 41 . Its symbol is fDn g with Dn = (2n + 1)n(n + 1), its order is 32 . Thus,pthe condition DG 2 L2( ) is equivalent to G 2 H 32 ( ), more precisely, kDGkL2 6kGkH 23 . It is worth noting that the functions in the Sobolev space H 32 ( ) are "slightly smoother" than the continuous functions, and its elements generally are not Lipschitz-continuous. Observing our particular choice the generalized discrepancy is then given by "X # 21 1 N X N X 1 1 Pn (i t) : D(f1; : : : ; N g; D) = N i=1 t=1 n=1 4n(n + 1) N !1
9
As mentioned above, the series on the right hand side is expressible in terms of elementary functions 2N N 3 12 s X X 1 1 ? i t 5 (1 ? 2 ln(1 + D(f1; : : : ; N g; D) = 2pN 4 (17) 2 )) : i=1 t=1
To keep an appropriate balance between computation and the requirement that the order should be as small as possible, we simply speak in the following of equidistribution on the sphere, when f1N ; : : : ; NN g, N=1,2,: : : , is D-equidistributed, i.e. lim D(f1N ; : : : ; NN g; D) = 0:
N !1
In connection with Theorem 3.1 we then obtain Z N p X j 41 G()d!() ? N1 G(i)j 6D(f1 ; : : :; N g; D)kGkH 23 ; i=1
(18)
where D(f1; : : : ; N g; D) is given by (17) and G 2 H 23 ( ). To get an optimal point system on the sphere we have to minimize the expression D(f1 ; : : :; N g; D) within a set of N points 1; : : : ; N on . This is unfortunately a non-linear optimization problem. It is a challenge for further investigation. Another possibility for construction of pointsets is to require that sum of some low degree spherical harmonics vanish, that is to say, moments up to a particular order should vanish.
De nition P3.3 (m-design). A pointset f1; : : :; N g is called an m-design if the so-called Weyl sums
N Y ( ) i=1 n;j i N X i=1
vanish for n = 1; : : : ; m. To be speci c,
Yn;j (i) = 0 for n = 1; 2; : : : m; j = 1; : : : 2n + 1 :
Equivalently a pointset f1; : : :; N g is an m-design if and only if N 1 Z Y ()d!() = 1 X 4 n N i=1 Yn (i) ;
(19)
(20)
holds for all spherical harmonics Yn of degree n m. It can be shown that a pair of antipodal points, the vertices of a regular tetrahedron, the regular octahedron, and the regular icosahedron give 1-,2-,3-,and 5-designs, respectively. In [15] the existence of spherical m-designs is proved, but only for suciently large N. A construction of mdesigns with number of points N = m4:5 can be found in [2]. In this paper we only deal with the possibility of an estimate of the generalized discrepancy in case of an m-design f1; : : :; N g. 10
Lemma 3.2 . Let A be a pseudodi erential operator of order s, s > 1, with symbol fAng, An 6= 0 for n 1. If a pointset f1; : : :; N g is an m-design, then the generalized discrepancy satis es
D(f1; : : : ; N g; A) mCs?1 ; where the constant C depends only on s.
Proof. From the de nitions of m-design and generalized discrepancy it follows that D(f1; : : :; N g; A) "X # 21 1 2n + 1 N X N X 1 = N 2 Pn (i t ) i=1 t=1 n=1 4An 2 1 N N 3 12 X X X 2 n + 1 1 5 = N4 2 Pn (i t ) n=m+1 i=1 t=1 4An 3 12 2 1 N N X X X 2 n + 1 1 5 N4 4A2 : n=m+1 i=1 t=1
n
By using the well-known inequality
1 1 X Cz;s ; (21) s ms?1 n=m+1 n where the Zeta constant is given by C = P1 1 , we have the desired estimate. k=1 ks
z;s
For large values s the m-design gives a good rate of convergence. Suppose that the number of points is N = mt (for t 4:5 we know the existence of an m-design, but the minimal value of t is unknown), then the rate of convergence amounts to N (s?1)=t. However, for s = 32 and t = 1, we obtain a rate of convergence of order N 12 . Therefore this type of an m-design does not seem to be well-suited.
4 Point Systems on the Sphere and Their Discrepancy Of practical importance is the basic problem of generating a large number of points on the sphere having a small generalized discrepancy. The geometrical interpretation of a pointset with small generalized discrepancy means that the points in the set are in a well-separated and distributed manner. A point system is called hierarchical, if +1g for all N. In our approach we distinguish ve kinds f1N ; : : : ; NN g f1N +1; : : :; NN+1 of point systems on the sphere: lattices in polar coordinates, transformed 2-dimensional 11
sequences, rotations on the sphere, triangulations, and the "sum of three squares sequence". In this section brie y describe how these point systems are generated. Then we numerically compute the generalized discrepancies and compare their values. To investigate the order of convergence in quantative way a parameter determined by 1 ; : : :; N g; D) ; D(f1; : : : ; N g; D) = N ? ; i.e., = ? log D(flog (22) N is introduced, where D(f1 ; : : :; N g; D) is given by (17). A large number of means that the point system is well-structured. The results of computation exhibit the di erent orders of convergence for di erent point systems.
4.4.1 Conventional Polar Coordinate Lattices
The lattices are generated by taking an equidistant division of the latitude and the longitude. Two examples are considered: (a) The "simple lattice" f(m; n); m = 1; : : : ; P ; n = 1; : : : ; Qg :
m = m=P ; m = 1; 2; : : : P ? 1 ; n = 2n=Q ; n = 1; 2; : : : Q ? 1 : This is the most well-known point system in geosciences. Unfortunately, for large numbers P and Q, there is strong concentration of points near the poles. Table 1 gives the computed values of the generalized discrepancy and the factor for this "simple lattice". As expected the point system does not give good results.
Tab.1: The generalized discrepancy of "simple lattices" No. of Pts. D(f1 ; : : :; N g; D) No. of Pts. D(f1 ; : : :; N g; D) 4 16 36 64
0.0972 0.0445 0.0354 0.0328
1.6817 1.1222 0.9320 0.8218
121 169 225 289
0.0315 0.0312 0.0310 0.0309
0.7211 0.6762 0.6415 0.6138
(b) To avoid the concentration of points near the poles we keep equidistant laterals
m = m=P ; m = 1; 2; : : : ; P ? 1 ; and divide each lateral circle into
n(m) = integer(2P sin(m=P )) parts, i.e.
n = 2n=n(m) n = 1; 2; : : : ; n(m) : 12
Hence, this "improved lattice" does not show strong concentration of points near the poles, since the division of lateral circles depends on the variable . The total number of points in the set f(m; n); m = 1; 2; : : : ; P ? 1 ; n = 1; 2; : : : ; n(m)g amounts to PX ?1 (23) NP = n(m) : m=1
Table 2 gives the computed values of the generalized discrepancy and the factor for the "improved lattices". Actually the rate of improvement is slow.
Tab.2 Generalized discrepancy of "improved lattice" No. of Pts. D(f1 ; : : :; N g; D) No. of Pts. D(f1 ; : : :; N g; D) 6 20 66 120 162 204
0.0812 0.0280 0.0179 0.0162 0.0162 0.0158
1.4016 1.194 0.9600 0.8610 0.8099 0.7797
234 276 296 350 394 442
0.0159 0.0158 0.0158 0.0158 0.0158 0.0158
0.7593 0.7382 0.7294 0.7079 0.6943 0.6811
4.4.2 Transformed Sequences
We again follow the strategy: (i) consider equidistributed pointsets on the rectangle [?1; 1] [0; 2) (ii) transform this system via the mapping
p
p
(t; ) 7! = ( 1 ? t2 cos ; 1 ? t2 sin ; t)T : to the unit sphere . As a matter of fact, there are many ways to generate equidistributions in a rectangle (for example, see [8,9,10,13] and many other references). (a) Hammersley system: the p-adic Van der Corput sequence is de ned by r with n = a0 + : : : + ar pr ; (24) x(np) = ap0 + : : : + par+1 where p is an integer greater than or equal 2. 2n?1 Let tn = x(2) n be the Van der Corput sequence with basis p=2 and n = 2N . Then, the Hammersley system on the sphere is de ned by (tn; n) n = 1; 2 : : : N . The Hammersley system is not hierarchical. A hierarchical point system can be generated by two p-adic Van der Corport sequences (tn; n) = (x(np1); x(np2)); n = 1; 2; : : :, with di erent prime numbers p1 and p2. Table 3 gives the computed values of the generalized discrepancy and the factor for the Hammersley system. The values of show a signi cant improvement compared with the former pointsets.
Tab.3: Generalized discrepancy of Hammersley systems 13
No. of Pts. D(f1 ; : : :; N g; D) 5 0.0813 23 0.0262 45 0.0157 85 0.0097 125 0.0072 165 0.0060
No. of Pts. D(f1 ; : : :; N g; D) 1.5590 218 0.0048 1.1612 250 0.0043 1.0908 296 0.0038 1.0427 350 0.0034 1.0218 394 0.0031 1.0035 440 0.0028
0.9919 0.9874 0.9779 0.9712 0.9671 0.9637
(b) Uniformly pseudo-random generation: we use the one-dimensional pseudo-random number generator to dertermine two sequences tj ; l j; l = 1; 2 : : :. For example, we choose a classical congruential generator, xn+1 = axn mod N with a = 31415821, N = 108 and x0 = 10000. Table 4 gives the computed values of the generalized discrepancy and the factor for the pseudo-random sequence on the sphere.
Tab.4: Generalized discrepancy of pseudo-random sequences No. of Pts. D(f1 ; : : :; N g; D) No. of Pts. D(f1 ; : : :; N g; D) 22 62 102 142 202
0.0542 0.0341 0.0296 0.0194 0.0189
0.9430 0.8185 0.7614 0.7957 0.7480
222 262 302 342 402
0.0214 0.0174 0.0202 0.0152 0.0142
0.7113 0.7272 0.6833 0.7173 0.7092
4.4.3 Rotations on the Sphere
This so-called operator sequence is described in [11]. Let A, B, and C be rotations (orthogonal transformations), respectively, about the X, Y and Z-axes, each through an angle. Let Rs be the set of non-trivial reduced words in A, B, C, A?1, B ?1, C ?1 of length less than or equal to s (by reduced we mean all the obvious cancellations such as AA?1 have been carried out), i.e. Rs = fA; B; C; A?1; B ?1; C ?1; AA; AB; AC; AB?1; AC ?1; : : :g. Rs contains 2N= 23 (5s ? 1) elements of rotations. A stereographic projection can be used to relate points in the complex plane C to points on the sphere. To every point = (1; 2; 3)T on the unit sphere, except the North pole (0; 0; 1)T , we associate a complex number i2 : (25) z = 11 + ? 3
Under this stereographic projection the fractional linear transformations in the complex plane correspond to rotations on the sphere. Taking this into account it is convenient to actually compute with points in C and fractional linear transformations instead of 14
rotations. In particular, we consider the following fractional linear transformations 0 0 1 1 0 1 1 + 2 i 0 1 2 1 2 i 1 1 1 A ;B = p @ A ;C = p @ A: A= p @ 1 ? 2i 5 0 5 ?2 1 5 2i 1 Table 5 gives the computed values of the generalized discrepancy and the factor with initial complex number i as starting point for rotations.
Tab.5: Generalized discrepancy of rotations No. of Pts. D(f1 ; : : :; N g; D) No. of Pts. D(f1 ; : : :; N g; D) 12 52 92 132 192
0.0577 0.0443 0.0419 0.0450 0.0416
1.1478 0.7887 0.7015 0.6353 0.6048
212 252 292 332 412
0.0396 0.0373 0.0357 0.0371 0.0390
4.4.4 Triangulations
0.6026 0.5948 0.5872 0.5676 0.5387
One can generate point systems on the sphere by triangulation: the triangulation are initiated with well-known polyhedral corner points of number 4, 8 or 20, which are, respectively, tetrahedron, octahedron, or icosahedron. To re ne the triangulation we proceed as follows: connect the midpoints of sides of each triangle, then project the midpoints on the sphere in radial direction, this will form four new triangular elements. The points are taken by the centre of each new triangular elements. Table 6,7,8 give the computed values of the generalized discrepancy and the factor for triangulations. To our surprise the results are not very convincing for larger pointsets.
Tab.6: Generalized discrepancy of tetrahedral tringulation No. of Pts. D(f1 ; : : :; N g; D) 4 16 64 256 1024
0.1115 0.1038 0.0723 0.0712 0.0693
1.5827 0.8172 0.6315 0.4766 0.3851
Tab.7: Generalized discrepancy of octahedral tringulation
15
No. of Pts. 8 32 128 512 2048
D(f1 ; : : :; N g; D) 0.0621 0.1743 0.1659 0.1705 0.1753
1.3367 0.5041 0.3703 0.2836 0.2284
Tab.8: Generalized discrepancy of icosahedral triangulation No. of Pts. D(f1 ; : : :; N g; D) 20 80 320 1280
0.0462 0.0883 0.0759 0.0792
1.0262 0.5539 0.4469 0.3545
4.4.5 Sum of Three Squares Sequence
This sequence arises in Number Theory, and a description can be found in the paper [1]. Taking the projection of the integer solutions of the equations x21+x22+x23 = n; n 2 IN, on the unit sphere, we are led to sets of type (26) X^N = f p1n (x1; x2; x3)jx21 + x22 + x23 = n; x1; x2; x3 2 ZZ; n = 1; 2; : : : ; mg : Removing the redundancy in the set X^N we obtain the so-called "sum of three squares sequence", denoted by XN . This point system is hierarchical, since XN XN +1 for all N. Table 9 gives the computed values of the generalized discrepancy and the factor for the sum of three squares sequence.
Tab.9: Generalized discrepancy of sum of three squares sequence No. of Pts. D(f1 ; : : :; N g; D) No. of Pts. D(f1 ; : : :; N g; D) 12 52 92 132 172
0.1067 0.0520 0.0408 0.0384 0.0364
4.4.6 Conclusion
0.9006 0.7484 0.7077 0.6678 0.6325
212 252 292 332 392
0.0378 0.0354 0.0372 0.0345 0.0350
0.6116 0.6041 0.5797 0.5800 0.5613
The results of the generalized discrepancy for di erent types of point systems on the sphere will be plotted below. It is surprising that certain types of transformed sequences yield the best results in our computation. It is also worth mentioning that the triangu16
lations do not take an exceptional role in equidistribution on the sphere. Furthermore it should be noted that there exist some well-distributed pointsets with special number of points on the unit sphere. For example, the soccer ball (the carbon-60 molecule) has the smallest generalized discrepancy (=0.0012...) compared with all our pointsets with about 60 elements. 0.05 0.045
Rotations
generalized discrepancy
0.04 0.035
Sum of Three Squares Simple Lattice
0.03 0.025 0.02 0.015
Improved Lattice Random
0.01
Hammersley
0.005 0 150
200
250
300 Number of Points
350
400
450
5 Low Discrepancy Method The Monte-Carlo method is very simple, namely the integral is replaced by its arithmetic mean value over a nite set of points. The quasi-Monto-Carlo method, sometimes called low discrepancy method, is understood to be a Monte-Carlo method, but based on pointsets forming an equidistribution. The proper justi cation of the low discrepancy method must not be based on the randomness of the procedure, which is spurious, but on equidistribution of the sets of points at which the values of the integrand are computed. Besides the simplicity, the other advantage of the low discrepancy method is its exibility of adaptive increasing number of points in order to guarantee a given accuracy. The transition from N points to N+1 points for a hierachical system needs only one further operation. The drawback of the low discrepancy method is that quite a large number of sampling may be needed, since the convergence tends to be quite slow. 17
From Section 4, we know empirically the factor . In the following test, the Hammersley point system is used, which turned out in our study to have a good value of . Example 1: rstly we discuss the integration of a discontinuous function, namely the characteristic function given by 8 < (27) h( 0) = : 1 if h 0 1 0 elsewhere. for xed 0 2 . The exact value of the integral is easily calculable: Z h( 0)d!() = 2(1 ? h) :
(28)
For our tests we choose the parameter h = cos(0:1) and 0 = (0; 0; 1)T , the exact value of the integral is equal to 0.30752... . Example 2: the second problem deals with integration of the function (cf. [5]) Z cos(r0 )d!() = 4 sin(r r) : (29)
This function is smooth for small r. For the numerical tests, r = 1 and 0 = (0; 0; 1)T , R are taken, therefore we obtain cos(r0 )d!() = 0:84147:::. Example 3: the third problem is concerned with the computation of a highly oscillating integrand. We choose the square of a spherical harmonic with degree n = 100 and order j = 50. Of course, the integral value of this square of the spherical harmonic is equal to 1.
Absolute Error of Low Discrepancy Method No. of Points 100 300 500 700 900 1100 1300 1500
Example 1 5.6194e-002 1.4306e-002 5.9281e-003 2.3377e-003 3.4303e-004 9.2630e-004 1.8051e-003 2.4495e-003
Example 2 4.3443e-004 4.9351e-005 3.2776e-005 1.4992e-005 1.1899e-005 5.7437e-006 4.1731e-006 3.8074e-006
Example 3 2.4112e-002 3.0263e-002 3.5989e-002 2.5217e-003 3.2383e-003 1.5361e-002 1.9815e-003 1.9823e-003
18
References [1] Arenstorf, R. F. and Johnson, D. (1979) Uniform Distribution of Integral Points on 3-Dimensional Sphere via Modular Forms, Journal of Number Theory 11, 218-238. [2] Bajnok, B. (1991) Construction of Design on 2-Sphere, Europ. J. Combinatorics, (12), 377-382. [3] Cui, J. (1995) Finite Pointset Methods on the Sphere and Their Application in Physical Geodesy, Ph. D. Thesis, University of Kaiserslautern, Geomathematics Group, Germany. [4] Delsarte, P., Goethals, J. M. and Seidel, J. J. (1977) Spherical Codes and Design , Geom. Dedicata 6, 363-388. [5] Freeden, W. (1980) On Integral Formulas of the Unit Sphere and Their Application to Numerical Computation of Integrals, Computing 25, 131-146. [6] Grabner, P. J. and Tichy, R. F., (1993) Spherical Designs, Discrepancy and Numerical Integration, Mathematics of Computation, Vol. 60, No. 201, 327-336. [7] Gronwall, T. (1914) On the Degree of Convergence of Laplace's Series. Trans. Amer. Math. Soc., 15 pp. 1-30. [8] Hlawka, E. (1982) Gleichverteilung auf Produkten von Spharen. J. Reine Angew. Math. 30, 1-30. [9] Hua, L., Wang, Y. (1981) Application of Number Theory to Numerical Analysis, Springer Verlag and Science Press, Beijing. [10] Kuipers, L. and Niederreiter, H. (1974) Uniform Distribution of Sequences, Wiley, New York. [11] Lubotzky, A., Phillips, R., and Sarnak, P. (1986) Hecke Operators and Distributing Points on the Sphere I, Comm. on Pure and Appl. Math., Vol. XXXIX, 149-186. [12] Magnus, W., Oberhettinger, F., Soni, R. (1966) Formulas and Theorems for the Special Functions of Mathematical Physics, Grundlehren der Mathematischen Wissenschaften in Einzeldarstellungen, Bd. 52 Springer Verlag. [13] Moroko , W. and Ca isch, R. E, (1994) Quasi-Random Sequences and Their Discrepancies, SIAM J. Sci. Comput., Vol 15, No. 6, 1251-1279. 19
[14] Muller, C.(1966) Spherical Harmonics, Lecture Notes in Mathematics 17, Springer Verlag. [15] Seymour, P. D. and Zaslavsky, T. (1984) Averaging Set: a Generalization of Mean Values and Spherical Designs, Adv. Math. (52), 213-240. [16] Sloane, N. (1983) Encrypting by Random Rotations, Cryptograph, Lecture Notes in Computer Science, 149, T. Beth editor, Berlin. [17] Svensson, S. (1983) Pseudodi erential Operators - A New Approach to the Boundary Problems of Physical Geodesy, Manuscripta Geodeaetica, Vol.8, 1-40. [18] Tichy, R. F., (1990) Random Points in the Cube and on the Sphere with Applications to Numerical Analysis, Journal of Comp. and Appl. Math., 31, 191-197. [19] Wahba, G. (1981) Spline Interpolation and Smoothing on the Sphere, SIAM J. Sci. Comput., 2 1-15.
20
Jianjun Cui and Willi Freeden Laboratory of Technomathematics, Geomathematics Group University of Kaiserslautern, 67653 Kaiserslautern Email: [email protected] Key Words: sphere, pointsets, pseudodi erential operators, generalized discrepancy,
equidistribution, low discrepancy (quasi-Monte-Carlo) method, approximate integration. AMS subject classi cation: 33C55, 40A10, 47G30, 62E25, 86-08
Abstract A concept of generalized discrepancy, which involves pseudodi erential operators to give a criterion of equidistributed pointsets, is developed on the sphere. A simply structured formula in terms of elementary functions is established for the computation of the generalized discrepancy. With the help of this formula ve kinds of point systems on the sphere, namely lattices in polar coordinates, transformed 2-dimensional sequences, rotations on the sphere, triangulation, and "sum of three squares sequence", are investigated. Quantitative tests are done, and the results are compared with each other. Our calculations exhibit di erent orders of convergence of the generalized discrepancy for di erent types of point systems.
1 Introduction Of practical importance is the problem of generating equidistributed pointsets on the sphere. For that reason a concept of generalized discrepancy, which involves pseudodifferential operators to give a quantifying criterion of equidistributed pointsets, is of great interest. In this paper an explicit formula in terms of elementary functions is developed for the generalized discrepancy. Essential tools are Sobolev space structures and pseudodi erential operator techniques. It is mentioned that an optimal pointset may be obtained by minimizing the generalized discrepancy. But, in spite of the elementary representation of the generalized discrepancy, this is a non-linear optimization problem which will not be discussed here. Our investigations, however, show that there are many promising ways to generate point systems on the sphere such that the discrepancy becomes small. To be speci c, we distinguish ve kinds of point systems on the sphere: lattices in polar coordinates, transformed 2-dimensional sequences, rotations on the sphere, triangulations, and "sum of three squares sequence". By use of our developed formulas, the ve classes of point systems are described, and their discrepancies are 1
explicitly calculated for increasing numbers of points. The results show di erent orders of convergence indicated by the generalized discrepancy. Furthermore, our computations enable us to give a quantitative comparison between the di erent point systems. It is somehow surprising that certain types of transformed sequences yield the best results. Nevertheless, there are special other pointsets which provide us with better results for comparable numbers of points. For instance, the soccer ball (C60) leads us to the best result in all our considered pointsets of about 60 points. The problem of generating a large number of "equidistributed points" on the sphere has many applications in various elds of computation, particularly in geoscience, medicine. The advantage of equidistributed point systems lies in the fact that relatively few sampling of the data is needed, and approximate integration can be simply performed by computation of a mean value, i.e., the arithmetical mean. The outline of the paper is as follows: Section 2 begins with a discussion of Sobolev spaces and (invariant) pseudodi erential operators on the sphere. The generalized discrepancy is developed in Section 3. The announced point systems are studied in Section 4. The low discrepancy method, also called quasi-Monte-Carlo method, for numerical integration on the sphere is discussed in Section 5 for three types of test examples.
2 Sobolev Spaces and Pseudodi erential Operators on the Sphere Let IR3 denote the three-dimensional Euclidean space. We use x; y; z : : : to represent the elements of IR3 and the Greek alphabet ; : : : to represent the vectors of the unit sphere
in IR3. denotes the Beltrami operator on the unit sphere. A function F : 7! IR possessing k continuous derivatives on is said to be of class C k ( ). C ( ) = C 0( ) is the class of real continuous scalar-valued functions on . By L2( ) we denote the space of Lebesgue square-integrable scalar functions on . The spherical harmonics of degree n are de ned as the everywhere on in nitely di erentiable eigenfunctions of the Beltrami-operator corresponding to the eigenvalues ()n = ?n(n +1); n = 0; 1; : : : ;. As it is well-known (cf. e.g. [14]), the linear space Harmn of all spherical harmonics of degree n is of dimension 2n+1. We denote fYn;j (); n = 0; 1; : : : ; j = 1; 2; : : : ; 2n +1)g to be an orthonormalized basis of L2( ), where n is called degree and j order of the spherical harmonics. The space Harm0;:::;m of all spherical harmonics of degree m possesses the dimension (m + 1)2. The Legendre polynomials Pn are the in nitely di erentiable eigenfunctions of the Legendre operator, satisfying the normalization condition Pn (1) = 1. It should be noted that jPn (t)j 1 and jPn0 (t)j n(n2+1) for all t 2 [?1; +1] and all 2
n 2 IN0. The well-known addition theorem of spherical harmonics states (cf. [5], [14]) 2X n+1 1 P ( ) ; (; ) 2 : Yn;j ()Yn;j () = 2n4+ (1) n j =1 In particular, we have s 2X n+1 2 n + 1 jYn;j ()j2 = 4 ; jYn;j ()j 2n4+ 1 ; 2 : j =1
(2)
For s 2 IR consider the space 1 2X n+1 X 2 n^ 2s < 1g; Fn;j
E s ( ) = fF 2 C 1( )j
n=0 j =1
where Fn;j are the spherical harmonic coecients of F and by de nition, 8 < n^ = : 1 if n = 0 n elsewhere: On E s ( ) we are able to impose an inner product as follows: n+1 1 2X X Fn;j Gn;j n^2s: < F; G >H s = n=0 j =1
(3)
(4)
(5)
The Sobolev space H s( ) is the completion of E s( ) with respect to the topology kkH s = (; )1H=s2. The following relationships between Sobolev spaces are valid: H s1 ( ) H s2 ( ), whenever s1 > s2. Sobolev spaces and classical function spaces C k ( ) are related via imbedding theorems.
Lemma 2.1 (Imbedding Theorem). H s ( ) is a subspace of C k ( ), if s > k + 1. Proof. It is sucient to prove Lemma 2.1 for the case k = 0. For s > 1, applying the Cauchy-Schwarz inequality we obtain n+1 1 2X n+1 1 2X 1 2n+1 X X 2 n^ 2s ) 21 ( X X (Yn;j ( ))2 n^ ?2s ) 21 : Fn;j jFn;j Yn;j ()j ( n=0 j =1
n=0 j =1
n=0 j =1
From the addition theorem of the spherical harmonics it follows that 1 1 2X n+1 X X 1 < 1; (Yn;j ())2 n^?2s = 24nn+ C= 2 s ^ n=0 n=0 j =1 provided that s > 1. Therefore the series is absolutely convergent. But this means that the F corresponds to a continuous function on . 3
From Lemma 2.1 (Imbedding Theorem) we are able to conclude that the spherical harmonic expansion of any function in H s ( ) converges uniformly provided that s > 1. This is of signi cance, since there are functions in C ( ) which do not allow a uniformly convergent spherical harmonic expansion (cf. [7]). By aid of the Sobolev space structure introduced above we make some comments on (invariant) pseudodi erential operators on the sphere, since these operators do not belong to the traditional equipment. De nition 2.1 (Pseudodi erential Operator on the Sphere). Let fAng, n=0,1,: : : , be a sequence of real numbers An satisfying jAnj nlim !1 (n + 12 )t = const 6= 0 for some t 2 IR. Then the operator A : H s( ) 7! H s?t ( ) de ned by 1 2X n+1 X AnFn;j Yn;j ; F 2 H s ( ) (6) AF = n=0 j =1
is called pseudodi erential operator of order t. fAn g is called spherical symbol of A. Moreover, if jAnj = 0 lim n!1 (n + 1 )t 2 s for all t 2 IR, then the operator A:H ( ) 7! C 1 ( ) is called pseudodi erential operator of order -1 (note that the equation (6) is understood in H s?t ( )-topology).
The spherical symbol has many appealing properties. For example, it is readily seen that (A0 + A00)n = A0n + A00n ; (A0A00)n = A0nA00n for every choice of pseudodi erential operators A0; A00. Obviously we have
AYn;j = AnYn;j ; n = 0; 1; : : : ; j = 1; : : : ; n + 1 (i.e. all pseudodi erential operators considered here are invariant operators). Using the kernel 1 n+1 1 2X X X 1 A P ( ); AnYn;j ()Yn;j () = 2n4+ KA ( ) = n n n=0
n=0 j =1
AF is obtainable by convolution as follows Z AF = KA F = F ()KA()d!() : The kernel satis es the property KA 2 H ?(t+) ( ) for all > 0 provided that t > ?1. In the case t = ?1 we have KA 2 C 1( ). In the sequel we list some particularly important pseudodi erential operators: 4
(1) The Beltrami operator is not invertible, since ()n = 0 for n = 0, but the operator ? + 41 has the symbol f(n + 21 )2g; n = 0; 1; : : : ; and, hence, has an inverse (? + 41 )?1 which is a pseudodi erential of order ?2. More generally, (? + 14 )s is a pseudodi erential operator of order 2s and has the spherical symbol f(n + 21 )2sg; n = 0; 1; : : :. Using the kernel K(+ 14 )?2s ; s > 1, we obtain Z F () = ( + 14 )s K(+ 14 )?2s ( )( + 41 )sF ()d!()
for all F 2 H s ( ) and 2 . This is easily seen by the Parseval identity for the system of L2( )-orthonormal spherical harmonics. (2) The integral operator of the single layer potential on given by Z AU () = 21 jU?() j d!() ; 2 ; (7) is a pseudodi erential operator of order -1 possessing the symbol fAng with An = (n + 21 )?1. This can be seen by inserting the expansion of j ? j?1 in terms of Legendre polynomials. A consequence is that A is invertible. Other types of geodetically relevant integral operators can be found in [17]. (3) The operator A = (?2)(? + 14 ) 21 has the symbol fAng, where An is given by An = (2n + 1)n(n + 1); n = 0; 1; : : :. A is not invertible, but the related operator B with symbol fBng given by 8 < Bn = : 1; n = 0 (2n + 1)(n(n + 1)); n = 1; 2; : : : is invertible. Since (cf. e.g. [12])
0 s 1 1 P (t) = 1 ? 2 ln @1 + 1 ? t A ; t 2 [?1; +1] n 2 n=1 n(n + 1) 1 X
we nd for all ; 2 .
1 0 s 1 1 ? 1 A KB?1 ( ) = 2 ? 2 ln @1 + 2
(4) An = rn; 0 r < 1 (or An = e 21 n(n+1)h; h > 0) de ne pseudodi erential operators of order ?1. 5
To show the relationship between the pseudodi erential operators and the Sobolev spaces we give an equivalent reformulation of the Sobolev space H s ( ). Let A be a pseudodi erential operator of order s. Then, the Sobolev space H s ( ) can be equivalently expressed in the following form:
H s( ) = fF jF : 7! IR; AF 2 L2( )g:
(8)
3 Equidistribution on the Sphere In this section we discuss the problem of equidistribution. This problem can be stated as follows: nd a pointset f1; : : : ; N g such that the remainder term Z N X j G()d!() ? 4N G(i)j ; (9) i=1 for any function G 2 H s ( ) converges to 0 as N ! 1. Let us begin the discussion with an a priori estimate of the remainder. The result is an analogue to the Koksma-Hlawka inequality in Euclidean spaces (see also, [8,9,10] and many other references therein). Theorem 3.1 . Let A be a pseudodi erential operator of order s, s > 1, with symbol fAng satisfying An 6= 0 n 1. Then, for any function AG 2 L2( ), we have the estimate "X # 12 Z 1 2n + 1 N X N X N X 1 1 1 j 4 G()d!() ? N G(i)j N 2 Pn (i t ) kAGkL2 : (10) t=1 i=1 n=1 4An i=1 Proof. From (8) we know that AG 2 L2( ) is equivalent to G 2 H s ( ). Because of the assumption s > 1 it follows from Lemma 2.1 (Imbedding Theorem) that the spherical harmonic expansion 1 2X n+1 X Gn;j Yn;j (); 2 ; G() = n=0 j =1
converges absolutely and uniformly. We rewrite the spherical harmonic expansion of G as follows: Z n+1 Z AG( )Y ( ) 1 2X X 1 n;j G() = 4 G()d!() + d!()Yn;j (): An
n=1 j =1
By setting = i and taking the sum over all indices i; 1 i N , we obtain " X # Z N N n+1 Z AG( )Y ( ) 1 2X X 1 1 1X n;j d!() N Yn;j (i) : N i=1 G(i) = 4 G()d!() + n=1 j=1
An i=1 6
Hence, the integral error can then be estimated as follows Z N X j 41 G()d!() ? N1 G(i )j 1 2n+1 N i=1 Z 1 X X X Yn;j (i) A G ( ) Y N n;j ( )d! ( ) An
nZ=1 j=1 i=11 2n+1 N X X X Yn;j (i)Yn;j () = N1 AG() d! ( ) : An
n=1 j =1 i=1 By using the Cauchy-Schwarz inequality and the addition theorem of the spherical harmonics we get Z N X j 41 G()d!() ? N1 G(i )j i=1 2 0 3 21 12 1 Z Z N 1 2 n +1 X X X Yn;j (i)Yn;j () A N1 (AG())2d!() 2 64 @ d!()75 A n n=1 j =1 i=1 2 1 2n+1 PN !23 12 X X Y ( ) 1 n;j i i=1 5 = N kAGkL2 4 A n n=1 j =1 3 21 2 N N 1 2n+1 X X X X Y ( ) Y ( ) 1 n;j i n;j t 5 = N kAGkL2 4 2 A n i=1 t=1 n=1 j =1 "X # 21 1 2n + 1 N X N X 1 = N kAGkL2 4A2 Pn(i t) : i=1 t=1 n=1
n
This completes the proof. Theorem 3.1 extends a result given for the Beltrami operator (cf. [5]) to the general context of pseudodi erential operators. It is obvious that the right-hand side of the estimate (10) consists of two parts. The rst part depends only on the pointset, while the second part depends merely on the function G under consideration. This gives rise to the de nition of generalized discrepancy.
De nition 3.1 (Generalized discrepancy). Let A be a pseudodi erential operator of order s, s > 1, with symbol fAng, An = 6 0 for n 1. Then the generalized discrepancy associated to a pseudodi erential operator A is de ned by
"X # 21 N X 1 2n + 1 N X 1 (11) D(f1 ; : : :; N g; A) = N 2 Pn (i t) : i=1 t=1 n=1 4An The generalized discrepancy characterizes "how well a pointset is equidistributed". This can be explained by a heuristic geometrical interpretation as follows: for any subset 7
B of and any pointset f1; : : :; N g on we de ne the counting function, N X #B = B (i) ; i=1
(12)
where B is the characteristic function of B . Thus #B indicates the number n, 1 n N , of points i contained in B . Let B 2 H s, s > 1 be a "smooth" approximation of the characteristic function B. Then, for any non-empty Lebesgue-measurable subset B of
, we have an integral inequality of the following form Z Z N X j #NB ? 41 Bd!()j = j N1 B (i) ? 41 Bd!()j (13) i=1 Z N X j N1 B (i) ? 41 Bd!()j D(f1; : : : ; N g; A)kA B kL2 : i=1 Therefore, the generalized discrepancy gives a quantitative description of how well a pointset is equidistributed. Remark. It is also worth to translate the Wozniakowski idea of averaging the integration error (cf. [13]) into the spherical context. But this is a challenge for future work. For di erent pseudodi erential operators we have di erent generalized discrepancies. It is of interest to understand the relationship between di erent generalized discrepancies associated to di erent pseudodi erential operators. Lemma 3.1 . Let A, B be two pseudodi erential operators of order s1, s2, (s1 > 1, s2 > 1), and with symbols fAng, fBng satisfying An 6= 0 ; Bn 6= 0 for n 1, respectively. Then, the following relations are valid for the generalized discrepancies associated to the pseudodi erential operators A and B: 1) If the two pseudodi erential operators A and B are of the same order, i.e. s1 = s2, then there exist two positive constants C1 and C2 such that
C1D(f1; : : : ; N g; A) D(f1; : : : ; N g; B) C2D(f1; : : : ; N g; A) :
(14)
2) If s1 > s2, then there exists a constant C3 such that
D(f1; : : : ; N g; A) < C3D(f1 ; : : :; N g; B) :
Proof. Applying the addition theorem of spherical harmonics we have
n+1 Y ( )Y ( ) 1 2X 1 2n + 1 N X N X N X N X X X n;j i n;j t P ( ) = n i t 2 A2n i=1 t=1 n=1 j =1 i=1 t=1 n=1 4An n+1 PN Y ( ) !2 1 2X X i=1 n;j i : = An n=1 j =1
8
(15)
Thus the generalized discrepancy can be rewritten as follows: 2 1 2n+1 PN !23 12 X X Y ( ) i=1 n;j i 5 : D(f1; : : : ; N g; A) = 4 An n=1 j =1 If s1 = s2, then it follows from our de nition of pseudodi erential operators that there exist two constants C1 and C2 such that 1 2 1 2 1 2 C1 A B C2 A ; n 1: n n n Therefore it is clear that C1D(f1 ; : : :; N g; A) D(f1 ; : : :; N g; B) C2D(f1 ; : : :; N g; A) : The part (2) can be proved in the same manner. From Lemma 3.1 we deduce that the generalized discrepancy is dependent only on the order of the pseudodi erential operators. Using the generalized discrepancy we can study the property of point systems in quantitative way. De nition 3.2 (Equidistribution in H s ( )). A point system f1N ; : : :; NN g; N = 1; 2; : : : is called A-equidistributed in H s ( ), s > 1, if the generalized discrepancy associated to a pseudodi erential operator A of order s, s > 1, satis es lim D(f1N ; : : :; NN g; A) = 0 : (16) It should be noted that if a point system is equidistributed in H s0 ( ), s0 > 1, then the point system is equidistributed in H s( ), s > s0. This shows us that we should use s0, s0 > 1, as small as possible for our quantitative investigations. However, in order to compute the generalized discrepancy as illustrated above, one has to evaluate a series expansion in terms of Legendre polynomials. This summation is very time consuming, in particular, for relatively small s. Therefore it is not advisable to use this series expansion for the computation of the generalized discrepancy. Fortunately, for certain pseudodi erential operators, the series is representable in terms of elementary functions. One example should be discussed in more detail. Consider theqpseudodi erential operator D = (?2) 21 (? + 14 ) 41 . Its symbol is fDn g with Dn = (2n + 1)n(n + 1), its order is 32 . Thus,pthe condition DG 2 L2( ) is equivalent to G 2 H 32 ( ), more precisely, kDGkL2 6kGkH 23 . It is worth noting that the functions in the Sobolev space H 32 ( ) are "slightly smoother" than the continuous functions, and its elements generally are not Lipschitz-continuous. Observing our particular choice the generalized discrepancy is then given by "X # 21 1 N X N X 1 1 Pn (i t) : D(f1; : : : ; N g; D) = N i=1 t=1 n=1 4n(n + 1) N !1
9
As mentioned above, the series on the right hand side is expressible in terms of elementary functions 2N N 3 12 s X X 1 1 ? i t 5 (1 ? 2 ln(1 + D(f1; : : : ; N g; D) = 2pN 4 (17) 2 )) : i=1 t=1
To keep an appropriate balance between computation and the requirement that the order should be as small as possible, we simply speak in the following of equidistribution on the sphere, when f1N ; : : : ; NN g, N=1,2,: : : , is D-equidistributed, i.e. lim D(f1N ; : : : ; NN g; D) = 0:
N !1
In connection with Theorem 3.1 we then obtain Z N p X j 41 G()d!() ? N1 G(i)j 6D(f1 ; : : :; N g; D)kGkH 23 ; i=1
(18)
where D(f1; : : : ; N g; D) is given by (17) and G 2 H 23 ( ). To get an optimal point system on the sphere we have to minimize the expression D(f1 ; : : :; N g; D) within a set of N points 1; : : : ; N on . This is unfortunately a non-linear optimization problem. It is a challenge for further investigation. Another possibility for construction of pointsets is to require that sum of some low degree spherical harmonics vanish, that is to say, moments up to a particular order should vanish.
De nition P3.3 (m-design). A pointset f1; : : :; N g is called an m-design if the so-called Weyl sums
N Y ( ) i=1 n;j i N X i=1
vanish for n = 1; : : : ; m. To be speci c,
Yn;j (i) = 0 for n = 1; 2; : : : m; j = 1; : : : 2n + 1 :
Equivalently a pointset f1; : : :; N g is an m-design if and only if N 1 Z Y ()d!() = 1 X 4 n N i=1 Yn (i) ;
(19)
(20)
holds for all spherical harmonics Yn of degree n m. It can be shown that a pair of antipodal points, the vertices of a regular tetrahedron, the regular octahedron, and the regular icosahedron give 1-,2-,3-,and 5-designs, respectively. In [15] the existence of spherical m-designs is proved, but only for suciently large N. A construction of mdesigns with number of points N = m4:5 can be found in [2]. In this paper we only deal with the possibility of an estimate of the generalized discrepancy in case of an m-design f1; : : :; N g. 10
Lemma 3.2 . Let A be a pseudodi erential operator of order s, s > 1, with symbol fAng, An 6= 0 for n 1. If a pointset f1; : : :; N g is an m-design, then the generalized discrepancy satis es
D(f1; : : : ; N g; A) mCs?1 ; where the constant C depends only on s.
Proof. From the de nitions of m-design and generalized discrepancy it follows that D(f1; : : :; N g; A) "X # 21 1 2n + 1 N X N X 1 = N 2 Pn (i t ) i=1 t=1 n=1 4An 2 1 N N 3 12 X X X 2 n + 1 1 5 = N4 2 Pn (i t ) n=m+1 i=1 t=1 4An 3 12 2 1 N N X X X 2 n + 1 1 5 N4 4A2 : n=m+1 i=1 t=1
n
By using the well-known inequality
1 1 X Cz;s ; (21) s ms?1 n=m+1 n where the Zeta constant is given by C = P1 1 , we have the desired estimate. k=1 ks
z;s
For large values s the m-design gives a good rate of convergence. Suppose that the number of points is N = mt (for t 4:5 we know the existence of an m-design, but the minimal value of t is unknown), then the rate of convergence amounts to N (s?1)=t. However, for s = 32 and t = 1, we obtain a rate of convergence of order N 12 . Therefore this type of an m-design does not seem to be well-suited.
4 Point Systems on the Sphere and Their Discrepancy Of practical importance is the basic problem of generating a large number of points on the sphere having a small generalized discrepancy. The geometrical interpretation of a pointset with small generalized discrepancy means that the points in the set are in a well-separated and distributed manner. A point system is called hierarchical, if +1g for all N. In our approach we distinguish ve kinds f1N ; : : : ; NN g f1N +1; : : :; NN+1 of point systems on the sphere: lattices in polar coordinates, transformed 2-dimensional 11
sequences, rotations on the sphere, triangulations, and the "sum of three squares sequence". In this section brie y describe how these point systems are generated. Then we numerically compute the generalized discrepancies and compare their values. To investigate the order of convergence in quantative way a parameter determined by 1 ; : : :; N g; D) ; D(f1; : : : ; N g; D) = N ? ; i.e., = ? log D(flog (22) N is introduced, where D(f1 ; : : :; N g; D) is given by (17). A large number of means that the point system is well-structured. The results of computation exhibit the di erent orders of convergence for di erent point systems.
4.4.1 Conventional Polar Coordinate Lattices
The lattices are generated by taking an equidistant division of the latitude and the longitude. Two examples are considered: (a) The "simple lattice" f(m; n); m = 1; : : : ; P ; n = 1; : : : ; Qg :
m = m=P ; m = 1; 2; : : : P ? 1 ; n = 2n=Q ; n = 1; 2; : : : Q ? 1 : This is the most well-known point system in geosciences. Unfortunately, for large numbers P and Q, there is strong concentration of points near the poles. Table 1 gives the computed values of the generalized discrepancy and the factor for this "simple lattice". As expected the point system does not give good results.
Tab.1: The generalized discrepancy of "simple lattices" No. of Pts. D(f1 ; : : :; N g; D) No. of Pts. D(f1 ; : : :; N g; D) 4 16 36 64
0.0972 0.0445 0.0354 0.0328
1.6817 1.1222 0.9320 0.8218
121 169 225 289
0.0315 0.0312 0.0310 0.0309
0.7211 0.6762 0.6415 0.6138
(b) To avoid the concentration of points near the poles we keep equidistant laterals
m = m=P ; m = 1; 2; : : : ; P ? 1 ; and divide each lateral circle into
n(m) = integer(2P sin(m=P )) parts, i.e.
n = 2n=n(m) n = 1; 2; : : : ; n(m) : 12
Hence, this "improved lattice" does not show strong concentration of points near the poles, since the division of lateral circles depends on the variable . The total number of points in the set f(m; n); m = 1; 2; : : : ; P ? 1 ; n = 1; 2; : : : ; n(m)g amounts to PX ?1 (23) NP = n(m) : m=1
Table 2 gives the computed values of the generalized discrepancy and the factor for the "improved lattices". Actually the rate of improvement is slow.
Tab.2 Generalized discrepancy of "improved lattice" No. of Pts. D(f1 ; : : :; N g; D) No. of Pts. D(f1 ; : : :; N g; D) 6 20 66 120 162 204
0.0812 0.0280 0.0179 0.0162 0.0162 0.0158
1.4016 1.194 0.9600 0.8610 0.8099 0.7797
234 276 296 350 394 442
0.0159 0.0158 0.0158 0.0158 0.0158 0.0158
0.7593 0.7382 0.7294 0.7079 0.6943 0.6811
4.4.2 Transformed Sequences
We again follow the strategy: (i) consider equidistributed pointsets on the rectangle [?1; 1] [0; 2) (ii) transform this system via the mapping
p
p
(t; ) 7! = ( 1 ? t2 cos ; 1 ? t2 sin ; t)T : to the unit sphere . As a matter of fact, there are many ways to generate equidistributions in a rectangle (for example, see [8,9,10,13] and many other references). (a) Hammersley system: the p-adic Van der Corput sequence is de ned by r with n = a0 + : : : + ar pr ; (24) x(np) = ap0 + : : : + par+1 where p is an integer greater than or equal 2. 2n?1 Let tn = x(2) n be the Van der Corput sequence with basis p=2 and n = 2N . Then, the Hammersley system on the sphere is de ned by (tn; n) n = 1; 2 : : : N . The Hammersley system is not hierarchical. A hierarchical point system can be generated by two p-adic Van der Corport sequences (tn; n) = (x(np1); x(np2)); n = 1; 2; : : :, with di erent prime numbers p1 and p2. Table 3 gives the computed values of the generalized discrepancy and the factor for the Hammersley system. The values of show a signi cant improvement compared with the former pointsets.
Tab.3: Generalized discrepancy of Hammersley systems 13
No. of Pts. D(f1 ; : : :; N g; D) 5 0.0813 23 0.0262 45 0.0157 85 0.0097 125 0.0072 165 0.0060
No. of Pts. D(f1 ; : : :; N g; D) 1.5590 218 0.0048 1.1612 250 0.0043 1.0908 296 0.0038 1.0427 350 0.0034 1.0218 394 0.0031 1.0035 440 0.0028
0.9919 0.9874 0.9779 0.9712 0.9671 0.9637
(b) Uniformly pseudo-random generation: we use the one-dimensional pseudo-random number generator to dertermine two sequences tj ; l j; l = 1; 2 : : :. For example, we choose a classical congruential generator, xn+1 = axn mod N with a = 31415821, N = 108 and x0 = 10000. Table 4 gives the computed values of the generalized discrepancy and the factor for the pseudo-random sequence on the sphere.
Tab.4: Generalized discrepancy of pseudo-random sequences No. of Pts. D(f1 ; : : :; N g; D) No. of Pts. D(f1 ; : : :; N g; D) 22 62 102 142 202
0.0542 0.0341 0.0296 0.0194 0.0189
0.9430 0.8185 0.7614 0.7957 0.7480
222 262 302 342 402
0.0214 0.0174 0.0202 0.0152 0.0142
0.7113 0.7272 0.6833 0.7173 0.7092
4.4.3 Rotations on the Sphere
This so-called operator sequence is described in [11]. Let A, B, and C be rotations (orthogonal transformations), respectively, about the X, Y and Z-axes, each through an angle. Let Rs be the set of non-trivial reduced words in A, B, C, A?1, B ?1, C ?1 of length less than or equal to s (by reduced we mean all the obvious cancellations such as AA?1 have been carried out), i.e. Rs = fA; B; C; A?1; B ?1; C ?1; AA; AB; AC; AB?1; AC ?1; : : :g. Rs contains 2N= 23 (5s ? 1) elements of rotations. A stereographic projection can be used to relate points in the complex plane C to points on the sphere. To every point = (1; 2; 3)T on the unit sphere, except the North pole (0; 0; 1)T , we associate a complex number i2 : (25) z = 11 + ? 3
Under this stereographic projection the fractional linear transformations in the complex plane correspond to rotations on the sphere. Taking this into account it is convenient to actually compute with points in C and fractional linear transformations instead of 14
rotations. In particular, we consider the following fractional linear transformations 0 0 1 1 0 1 1 + 2 i 0 1 2 1 2 i 1 1 1 A ;B = p @ A ;C = p @ A: A= p @ 1 ? 2i 5 0 5 ?2 1 5 2i 1 Table 5 gives the computed values of the generalized discrepancy and the factor with initial complex number i as starting point for rotations.
Tab.5: Generalized discrepancy of rotations No. of Pts. D(f1 ; : : :; N g; D) No. of Pts. D(f1 ; : : :; N g; D) 12 52 92 132 192
0.0577 0.0443 0.0419 0.0450 0.0416
1.1478 0.7887 0.7015 0.6353 0.6048
212 252 292 332 412
0.0396 0.0373 0.0357 0.0371 0.0390
4.4.4 Triangulations
0.6026 0.5948 0.5872 0.5676 0.5387
One can generate point systems on the sphere by triangulation: the triangulation are initiated with well-known polyhedral corner points of number 4, 8 or 20, which are, respectively, tetrahedron, octahedron, or icosahedron. To re ne the triangulation we proceed as follows: connect the midpoints of sides of each triangle, then project the midpoints on the sphere in radial direction, this will form four new triangular elements. The points are taken by the centre of each new triangular elements. Table 6,7,8 give the computed values of the generalized discrepancy and the factor for triangulations. To our surprise the results are not very convincing for larger pointsets.
Tab.6: Generalized discrepancy of tetrahedral tringulation No. of Pts. D(f1 ; : : :; N g; D) 4 16 64 256 1024
0.1115 0.1038 0.0723 0.0712 0.0693
1.5827 0.8172 0.6315 0.4766 0.3851
Tab.7: Generalized discrepancy of octahedral tringulation
15
No. of Pts. 8 32 128 512 2048
D(f1 ; : : :; N g; D) 0.0621 0.1743 0.1659 0.1705 0.1753
1.3367 0.5041 0.3703 0.2836 0.2284
Tab.8: Generalized discrepancy of icosahedral triangulation No. of Pts. D(f1 ; : : :; N g; D) 20 80 320 1280
0.0462 0.0883 0.0759 0.0792
1.0262 0.5539 0.4469 0.3545
4.4.5 Sum of Three Squares Sequence
This sequence arises in Number Theory, and a description can be found in the paper [1]. Taking the projection of the integer solutions of the equations x21+x22+x23 = n; n 2 IN, on the unit sphere, we are led to sets of type (26) X^N = f p1n (x1; x2; x3)jx21 + x22 + x23 = n; x1; x2; x3 2 ZZ; n = 1; 2; : : : ; mg : Removing the redundancy in the set X^N we obtain the so-called "sum of three squares sequence", denoted by XN . This point system is hierarchical, since XN XN +1 for all N. Table 9 gives the computed values of the generalized discrepancy and the factor for the sum of three squares sequence.
Tab.9: Generalized discrepancy of sum of three squares sequence No. of Pts. D(f1 ; : : :; N g; D) No. of Pts. D(f1 ; : : :; N g; D) 12 52 92 132 172
0.1067 0.0520 0.0408 0.0384 0.0364
4.4.6 Conclusion
0.9006 0.7484 0.7077 0.6678 0.6325
212 252 292 332 392
0.0378 0.0354 0.0372 0.0345 0.0350
0.6116 0.6041 0.5797 0.5800 0.5613
The results of the generalized discrepancy for di erent types of point systems on the sphere will be plotted below. It is surprising that certain types of transformed sequences yield the best results in our computation. It is also worth mentioning that the triangu16
lations do not take an exceptional role in equidistribution on the sphere. Furthermore it should be noted that there exist some well-distributed pointsets with special number of points on the unit sphere. For example, the soccer ball (the carbon-60 molecule) has the smallest generalized discrepancy (=0.0012...) compared with all our pointsets with about 60 elements. 0.05 0.045
Rotations
generalized discrepancy
0.04 0.035
Sum of Three Squares Simple Lattice
0.03 0.025 0.02 0.015
Improved Lattice Random
0.01
Hammersley
0.005 0 150
200
250
300 Number of Points
350
400
450
5 Low Discrepancy Method The Monte-Carlo method is very simple, namely the integral is replaced by its arithmetic mean value over a nite set of points. The quasi-Monto-Carlo method, sometimes called low discrepancy method, is understood to be a Monte-Carlo method, but based on pointsets forming an equidistribution. The proper justi cation of the low discrepancy method must not be based on the randomness of the procedure, which is spurious, but on equidistribution of the sets of points at which the values of the integrand are computed. Besides the simplicity, the other advantage of the low discrepancy method is its exibility of adaptive increasing number of points in order to guarantee a given accuracy. The transition from N points to N+1 points for a hierachical system needs only one further operation. The drawback of the low discrepancy method is that quite a large number of sampling may be needed, since the convergence tends to be quite slow. 17
From Section 4, we know empirically the factor . In the following test, the Hammersley point system is used, which turned out in our study to have a good value of . Example 1: rstly we discuss the integration of a discontinuous function, namely the characteristic function given by 8 < (27) h( 0) = : 1 if h 0 1 0 elsewhere. for xed 0 2 . The exact value of the integral is easily calculable: Z h( 0)d!() = 2(1 ? h) :
(28)
For our tests we choose the parameter h = cos(0:1) and 0 = (0; 0; 1)T , the exact value of the integral is equal to 0.30752... . Example 2: the second problem deals with integration of the function (cf. [5]) Z cos(r0 )d!() = 4 sin(r r) : (29)
This function is smooth for small r. For the numerical tests, r = 1 and 0 = (0; 0; 1)T , R are taken, therefore we obtain cos(r0 )d!() = 0:84147:::. Example 3: the third problem is concerned with the computation of a highly oscillating integrand. We choose the square of a spherical harmonic with degree n = 100 and order j = 50. Of course, the integral value of this square of the spherical harmonic is equal to 1.
Absolute Error of Low Discrepancy Method No. of Points 100 300 500 700 900 1100 1300 1500
Example 1 5.6194e-002 1.4306e-002 5.9281e-003 2.3377e-003 3.4303e-004 9.2630e-004 1.8051e-003 2.4495e-003
Example 2 4.3443e-004 4.9351e-005 3.2776e-005 1.4992e-005 1.1899e-005 5.7437e-006 4.1731e-006 3.8074e-006
Example 3 2.4112e-002 3.0263e-002 3.5989e-002 2.5217e-003 3.2383e-003 1.5361e-002 1.9815e-003 1.9823e-003
18
References [1] Arenstorf, R. F. and Johnson, D. (1979) Uniform Distribution of Integral Points on 3-Dimensional Sphere via Modular Forms, Journal of Number Theory 11, 218-238. [2] Bajnok, B. (1991) Construction of Design on 2-Sphere, Europ. J. Combinatorics, (12), 377-382. [3] Cui, J. (1995) Finite Pointset Methods on the Sphere and Their Application in Physical Geodesy, Ph. D. Thesis, University of Kaiserslautern, Geomathematics Group, Germany. [4] Delsarte, P., Goethals, J. M. and Seidel, J. J. (1977) Spherical Codes and Design , Geom. Dedicata 6, 363-388. [5] Freeden, W. (1980) On Integral Formulas of the Unit Sphere and Their Application to Numerical Computation of Integrals, Computing 25, 131-146. [6] Grabner, P. J. and Tichy, R. F., (1993) Spherical Designs, Discrepancy and Numerical Integration, Mathematics of Computation, Vol. 60, No. 201, 327-336. [7] Gronwall, T. (1914) On the Degree of Convergence of Laplace's Series. Trans. Amer. Math. Soc., 15 pp. 1-30. [8] Hlawka, E. (1982) Gleichverteilung auf Produkten von Spharen. J. Reine Angew. Math. 30, 1-30. [9] Hua, L., Wang, Y. (1981) Application of Number Theory to Numerical Analysis, Springer Verlag and Science Press, Beijing. [10] Kuipers, L. and Niederreiter, H. (1974) Uniform Distribution of Sequences, Wiley, New York. [11] Lubotzky, A., Phillips, R., and Sarnak, P. (1986) Hecke Operators and Distributing Points on the Sphere I, Comm. on Pure and Appl. Math., Vol. XXXIX, 149-186. [12] Magnus, W., Oberhettinger, F., Soni, R. (1966) Formulas and Theorems for the Special Functions of Mathematical Physics, Grundlehren der Mathematischen Wissenschaften in Einzeldarstellungen, Bd. 52 Springer Verlag. [13] Moroko , W. and Ca isch, R. E, (1994) Quasi-Random Sequences and Their Discrepancies, SIAM J. Sci. Comput., Vol 15, No. 6, 1251-1279. 19
[14] Muller, C.(1966) Spherical Harmonics, Lecture Notes in Mathematics 17, Springer Verlag. [15] Seymour, P. D. and Zaslavsky, T. (1984) Averaging Set: a Generalization of Mean Values and Spherical Designs, Adv. Math. (52), 213-240. [16] Sloane, N. (1983) Encrypting by Random Rotations, Cryptograph, Lecture Notes in Computer Science, 149, T. Beth editor, Berlin. [17] Svensson, S. (1983) Pseudodi erential Operators - A New Approach to the Boundary Problems of Physical Geodesy, Manuscripta Geodeaetica, Vol.8, 1-40. [18] Tichy, R. F., (1990) Random Points in the Cube and on the Sphere with Applications to Numerical Analysis, Journal of Comp. and Appl. Math., 31, 191-197. [19] Wahba, G. (1981) Spline Interpolation and Smoothing on the Sphere, SIAM J. Sci. Comput., 2 1-15.
20
