In simulating interstellar dust clouds astrophysicists need high accuracy integration formulae for functions on the sphere. There are several known approaches to construct such cubature formulae. Aims are to obtain nodes inside the region of integration, to get positive weights to assure convergence results and to get a high degree of exactness.
Probably the best known approach are product formulae where univariate quadratures are composed to get multivariate formulae. Positive weights can be assured by using univariate Gauss formulae for this construction. For the sphere the main disadvantage is that lattice points lie much denser near the poles than in the equatorial region. In applications this point distribution is often not suitable.
Another approach are cubature formulae constructed with the aid of invariant theory (see e. g. Cools [6]) where the number of moment equations occuring in the calculation of nodes and weights are drastically reduced with the aid of consistency conditions. Here, classes of points with equal weights are used. The formulae in general do not have only positive weights and it cannot be decided in advance wether a formula of a certain degree exists or not.
In contrast to this, ideal theory makes it possible to decide the question of existence or nonexistence of a formula. Lower bounds for the number of nodes needed for a cubature formula of specific degree can also be obtained (see e. g. [6] and the references therein). Nodes are common zeros of certain orthogonal polynomials. However, the calculation of these zeros is problematic.
None of the formulae mentioned above is a multivariate Gauss formula which in the univariate case are known to be the formulae with highest degree of exactness and with general convergence due to their positive weights.
Of course there has also been some research concerning multivariate Gauss quadrature formulae (see e. g. Bos [3], Reimer [17, 18, 19]). It turned out that the existence of multivariate Gauss formulae for the sphere is closely related to the existence of cubature formulae with equal weights. Multivariate Gauss formulae for the sphere can even be characterized by this property [17]. Bos [3] has proved that nodes yielding equal weights also solve the Fejér problem, that is the square-sum of the Lagrange polynomials is equal to one for these nodes. This property leads to results concerning the existence or nonexistence of multivariate Gauss quadrature formulae [3].
Nodal systems yielding cubature formulae with equal weights on the sphere are also known as spherical designs. Here, the number of nodes is not restricted to the dimension of the corresponding polynomial space. There are several results on the existence and characterization of spherical designs (Delsarte, Goethals and Seidel [8], Seidel [21], Bannai and Damerell [1, 2]), but no algorithm for the calculation of these designs is known.
An estimation of the integration error for equally weighted cubature formulae was recently given by Cui and Freeden [7]. However, convergence for this approach has up to now only be proved for spherical designs.
After having finished our numerical calculations we discovered that Zhou
[25] has partially worked on the same subject. For his thesis he
minimized numerically several objective functions using quasi-Newton methods.
Zhou puts an emphasis on geometric considerations but does not inspect
cubature formulae. He restricts his computations to facility numbers
.
Steinacker, Thamm and Maier [23] computed quadrature formulae for the sphere, yielding better results then the approach of Reimer and Sündermann [20], but the algorithm for the node distribution is not an optimal one.
In this paper the terms facilities, nodes and points are used for the same object without further notion.
We propose a two-stage algorithm using first a stochastical algorithm and an improved quasi-Newton method to compute nearly equidistant nodes on the sphere and then to calculate corresponding weights which differ only little from equal weights. Nodes as well as weights are computed with high precision and the resulting cubature formulae are tested for several functions. We computed point systems for up to 1600 nodes and calculated cubature weights for up to 900 points.
This article is organized as follows. Section 2 summarizes results about spherical designs and multivariate Gauss quadrature formulae. In Section 3 our technique for computing nodes and weights is described. Section 4 contains implementation details and Section 5 presents numerical results.
For 
the unit sphere in 
is denoted by 
.
For 
let 
denote the linear space
of polynomial restrictions on in r variables of degree 
and let 
denote the corresponding space of homogeneous polynomials.
A system 
of nodes 
is called a
fundamental system with respect to 
,
if the corresponding evaluation functionals form a basis in its dual space.
This property ensures that the
corresponding
Lagrange fundamental polynomials 
are well defined
and that 
holds
(
the Kronecker symbol).
Further let 
be the surface
integral of .
The reproducing kernel of is denoted by 
and has the
representation
where 
is the surface area of and 
denotes the Gegenbauer polynomial of degree m with index r/2.
The reproducing kernel of is denoted by 
and has the representation
