Quadrature formulas for integration of multivariate trigonometric [PDF]

Quadrature formulas for integration of multivariate trigonometric polynomials on spherical triangles J. Beckmann∗, H. N. Mhaskar†, and J. Prestin‡

Abstract We describe an explicit construction of quadrature rules exact for integrating multivariate trigonometric polynomials of a given coordinatewise degree on a spherical triangle. The theory is presented in the more general setting of quadrature formulas on a compact subset of the unit hypersphere, Sq , embedded in the Euclidean space Rq+1 . The number of points at which the polynomials are sampled is commensurate with the dimension of the polynomial space.

1

Introduction

High accuracy quadrature formulas for integration on the unit sphere embedded in a Euclidean space are necessary in many such practical applications as geophysics and the numerical solutions of differential and integral equations [1, 2]. For example, many of these applications involve the solution of a pseudo– differential equation. Many software tools are available to do these tasks if the data is available in the form of spherical harmonic coefficients. In practice, however, the data is available only in the form of the values of the function in question. A straightforward use of Monte Carlo integration is not satisfactory from a theoretical point of view [9]. Therefore, a careful construction of quadrature formulas exact for computing integrals of spherical harmonics of high degree is desirable. The existence of quadrature formulas on the sphere, exact for polynomials of high degrees, is established in [10]. A novelty of this result is that the quadrature formulas are based on sufficiently dense data sets consisting of scattered sites, i.e., where no specific choice of the sites is required. The investigation was continued in [11, 12] where the existence of similar formulas was established for integration on spherical caps and other compact subsets of the Euclidean sphere. In contrast to the results in [10], the focus in these later papers is on reducing the diameter of the compact subset, holding the degree of the polynomials fixed, as in spline approximation. In order to construct such quadrature formulas numerically, one needs to know first the exact values of the integrals over the compact set in question for suitable polynomials forming a basis for the space of spherical polynomials of the desired degree. Once such integrals are known exactly for the basis polynomials, one can use the procedures outlined in [9] to obtain quadrature formulas based on scattered data. For the canonical orthonormal basis for the space of spherical polynomials, it is easy to compute analytically the integrals over the whole sphere and spherical caps. To the best of our knowledge, there are no such analytic formulas available for spherical triangles. Therefore, Larry Schumaker proposed that we should look for quadrature formulas known to be exact to evaluate the necessary integrals. The objective of this paper is to provide an explicit construction for integration on spherical triangles. This construction can be used to evaluate the exact integrals for basic polynomials of a fixed (low) degree on the triangles. We observe that it is not important for this purpose to use the minimal number of ∗ Institute of Mathematics, University of L¨ ubeck, Ratzeburger Allee 160, 23562, L¨ ubeck, Germany, email: [email protected]. † Department of Mathematics, California State University, Los Angeles, CA 90032 and California Institute of Technology, Pasadena, California, 91125, USA, email: [email protected]. The research of this author was supported, in part, by grant DMS-0908037 from the National Science Foundation and grant W911NF-09-1-0465 from the U.S. Army Research Office. ‡ Institute of Mathematics, University of L¨ ubeck, Ratzeburger Allee 160, 23562, L¨ ubeck, Germany, email: [email protected].

1

quadrature nodes, nor is it critically important to come up with numerically stable algorithms. Since the focus is on reducing the diameter of the domain of integration, we have developed quadrature formulas for integrating all multivariate trigonometric polynomials of a given degree on a triangle, based on nodes belonging to the spherical cap circumscribing the triangle. It is theoretically possible to obtain the quadrature formula based on a smaller number of nodes if one restricts oneself to those trigonometric polynomials which are functions on the sphere. In the case of S2 , a complete characterization for such polynomials as well as an explicit form for the reproducing kernel canRbe found in [4]. The idea behind the paper is the following. We wish to calculate K T (ξ)dµ∗q (ξ) where T is a multivariate trigonometric polynomial, K is a compact subset of the sphere Sq , and µ∗q is the area element of the sphere. First, we observe that any (multivariate) trigonometric polynomial can be expressed in terms of a reproducing kernel as an integral over the spherical cap circumscribing K (Theorem 3.1). This integral, being a tensor product integral, can be discretized using appropriate quadrature formulas, thereby expressing the trigonometric polynomials as a linear combination of reproducing kernels with finitely many “centers”. An analogous construction for algebraic polynomials is given by Hesse and Womersley in [8]. These centers are our quadrature nodes. The quadrature weights are obtained by integrating these kernels over K (Theorem 3.1). To evaluate these integrals, we no longer need to restrict ourselves to use centers on a spherical cap, but need to focus on the domain K of integration. In the case of a spherical triangle on S2 , we may express this integral as a double integral. The specific form of the reproducing kernel ensures that the inner integral will be a trigonometric polynomial multiplied by one of the four weight functions described in (4.7)–(4.10). We can then fix quadrature formulas to evaluate the outer integral, and with each of the nodes in these quadrature formulas, evaluate the quadrature formulas for the inner integral. This yields the desired integrals for the reproducing kernels, and hence, for the original problem. In Section 2 we describe the representation of a univariate trigonometric polynomial as a discrete integral over an arc of the circle, involving a suitable reproducing kernel. This enables us to express the quadrature formulas for integration over K in terms of the quadrature formulas for integration of reproducing kernels over K in Section 3. We describe this construction in the context of a spherical triangle on S2 in Section 4. This is summarized in the form of an algorithm in Section 5. Numerical results are presented in Section 6.

2

Integration on a circular arc

The goal of this section is to obtain an integral representation for all univariate trigonometric polynomials as a discrete integral over an arc, using a suitable reproducing kernel. For an integer n ≥ 0, let Hn (respectively, Πn ) denote the space of all trigonometric (respectively, algebraic) polynomials of degree at most n. We find it convenient to extend this notation also when n is not an integer; thus, for example, Πx is Π⌊x⌋ . If T ∈ Hn and γ ∈ [−π, π], then there exist unique algebraic polynomials Q ∈ Πn and R ∈ Πn−1 such that T (t + γ) = Q(cos t) + (sin t)R(cos t). (2.1) Therefore, we will obtain reproducing kernels for T in terms of those for algebraic polynomials with respect to different weights. R1 If w : [−1, 1] → [0, ∞) is a weight function, i.e., −1 |x|n w(x)dx < ∞ for n = 0, 1, . . ., it is well known that there exists a unique system of polynomials pk (w; ◦) with degree k and positive leading coefficients, such that for integers k, ℓ = 0, 1, . . ., Z

1

pk (w; x)pℓ (w; x)w(x)dx =

−1



1, if k = ℓ, 0, otherwise.

Writing Kn (w; t, u) =

n X

pk (w; cos t)pk (w; cos u),

k=0

2

(2.2)

it is obvious that for P ∈ Πn , Z π Z π w(cos u) P (cos u)Kn (w; t, u) P (cos u)Kn (w; t, u)w(cos u)(sin u)du = P (cos t). (2.3) | sin u|du = 2 −π 0 Proposition 2.1 Let −π ≤ α < β ≤ π,  (1 − t2 )−1/2 , if t ∈ (cos((β − α)/2), 1), w(α, β; t) := 0, otherwise,

(2.4)

and w(α, ˜ β; t) = (1 − t2 )w(α, β; t). Let γ = (α + β)/2 and Kn∗ (α, β; t, u) := Kn (w(α, β); t − γ, u − γ) + sin(t − γ) sin(u − γ)Kn−1 (w(α, ˜ β); t − γ, u − γ) . 2 Then T (t) =

β

Z

T (u)Kn∗ (α, β; t, u)du,

α

T ∈ Hn .

(2.5)

Proof. In this proof only, let ψ = (β − α)/2. We choose Q and R as in (2.1). In view of (2.3) and the fact that Kn (α, β; ±t, ±u) = Kn (α, β; t, u), we obtain that ψ

Z

Q(cos t) =

Q(cos u)Kn (w(α, β); t, u)du

0

Z

=

ψ

T (u + γ) + T (−u + γ) Kn (w(α, β); t, u)du 2

Z

ψ

0

1 2

=

T (u + γ)Kn (w(α, β); t, u)du.

(2.6)

−ψ

Similarly, (sin t)R(cos t) = = =

Z

Z

ψ

R(cos u)(sin t)Kn−1 (w(α, ˜ β); t, u)(sin2 u)du

0 ψ

T (u + γ) − T (−u + γ) (sin t)(sin u)Kn−1 (w(α, ˜ β); t, u)du 2

Z

ψ

0

1 2

T (u + γ)(sin t)(sin u)Kn−1 (w(α, ˜ β); t, u)du.

(2.7)

We deduce (2.5) from (2.1), (2.6), and (2.7), with an obvious change of variables.

2

−ψ

Next, we turn our attention to expressing the integral in (2.5) as a finite sum. It is well known that for each integer n ≥ 1, the polynomial pn (w) has exactly n distinct zeros in the interior of the smallest interval containing the support of w. We denote these zeros by xk,n (w) = cos(χk,n ),

χk,n = χk,n (w) ∈ [0, π],

k = 1, . . . , n.

(2.8)

We will denote by νn (w) the measure (on the “circle”) that associates with each t ∈ [−π, π] for which pn (w; cos t) = 0 the mass −1   n−1 X p2j (w; cos t) . (2.9) λk,n = λk,n (w) :=   j=0

It is well known [3, §I.3, §I.4] that Z Z P (cos t)dνn (w; t) =

π

P (cos t)w(cos t)| sin t|dt,

−π

3

P ∈ Π2n−1 .

By symmetry, it is clear that Z Z (sin t)P (cos t)dνn (w; t) =

π

(sin t)P (cos t)w(cos t)| sin t|dt = 0,

−π

P ∈ Π2n−1 .

Since any T ∈ H2n−1 can be written in the form P1 (cos t) + (sin t)P2 (cos t) for some P1 ∈ Π2n−1 , P2 ∈ Π2n−2 , we conclude that Z Z π T (t)dνn (w; t) = T (t)w(cos t)| sin t|dt, T ∈ H2n−1 . (2.10) −π

Theorem 2.1 Let n ≥ 1 be an integer, n1 be the integer part of (n + 2)/2, −π ≤ α < β ≤ π. For an ∗ integer m ≥ 1, let νm (α, β) be defined by Z Z ∗ f (t)dνm (α, β; t) = f (t + γ)dνm (w(α, β); t) (2.11) then for T ∈ Hn , we have

Z

β

T (t)dt =

α

and T (t) =

Z

Z

T (u)dνn∗1 (α, β; u)

∗ T (u)Kn∗ (α, β; t, u)dνn+1 (α, β; u),

(2.12)

T ∈ Hn .

(2.13)

Proof. In this proof only, let ψ = (β − α)/2. We observe that the measure νn∗1 (α, β) is supported on the arc (−ψ, ψ). Since n ≤ 2n1 − 1, we obtain from (2.10) that Z ψ Z β Z Z T (t + γ)dt = T (t)dt. T (t)dνn∗1 (α, β; t) = T (t + γ)dνn1 (w(α, β); t) = −ψ

α

This proves (2.12). Since Kn∗ (α, β; t, u) is a trigonometric polynomial of degree at most n in both t and u, the formula (2.13) follows from the reproducing formula (2.5) using the quadrature formula (2.12) with 2n in place of n. 2

3

Quadrature formula on compact sets

In this section, we will describe how the kernels described above can be used to construct quadrature rules on an arbitrary compact subset of the sphere. The main result (Theorem 3.1) is proved by expressing a multivariate trigonometric polynomial using the reproducing kernels, followed by a change of the order of integration. In the sequel, let Hqn denote the class of all trigonometric polynomials in q variables with coordinatewise degree at most n. By a spherical cap of radius ρ, center x0 , we mean the set Sqρ (x0 ) := {x ∈ Sq : x · x0 ≥ cos ρ} = {x ∈ Sq : kx − x0 k ≤ 2 sin(ρ/2)}. We assume the standard parameterization of Sq in terms of the angles θ1 , . . . , θq , where −π ≤ θ1 ≤ π and 0 ≤ θk ≤ π for k = 2, . . . , q. If x ∈ Sq , then the kth component of x is given by  Qq θj , k = 1,   j=1 sinQ q cos θk−1 j=k sin θj , 2 ≤ k ≤ q, xk = (3.1)   cos θq , k = q +1. The area measure µ∗q on Sq can be expressed in these coordinates as dµ∗q (x)

=

q Y

sink−1 (θk )dθk .

k=1

4

Note that dµ∗q = sinq−1 (θq )dθq dµ∗q−1 . If K is a subset of Sq , K ◦ ⊆ [−π, π] × [0, π]q−1 is the corresponding set of parameters, and f is an integrable function on K ◦ , we will abuse the notation to write Z

K

f dµ∗q

=

Z

f (θ1 , . . . , θq )

K◦

q Y

sink−1 (θk )dθk .

k=1

We will also identify the sets K and K ◦ . R In the next theorem, we express the integral K T dµ∗q , where K is a compact subset of Sq and T ∈ Hqn as a discrete integral based on points in a spherical cap containing K. In the sequel, we may assume that the center of this cap is x0 = n := (0, . . . , 0, 1), and write Sqρ := Sqρ (n). Theorem 3.1 Let q ≥ 2 be an integer, 0 < ρ ≤ π, and K ⊆ Sqρ be a compact subset. For ζ = (θ1 , . . . , θq ), ξ = (φ1 , . . . , φq ), and an integer n ≥ 1, we use the notation as in (2.11) to define the measure ∗ (ρ; ζ) dνn,q

:=

(

q Y

k−1

sin

)

∗ dνn+q (0, ρ; θq )

(θk )

k=1

q−1 Y

!

∗ ∗ (0, π; θk )dνn+1 (−π, π; θ1 ) dνn+k

k=2

,

and the kernel ∗ ∗ Kn,q (ρ; ζ, ξ) := Kn+q−1 (0, ρ; θq , φq )

q−1 Y

∗ Kn+k−1 (0, π; θk , φk )Kn∗ (−π, π; θ1 , φ1 ).

(3.2)

k=2

Then for T ∈ Hqn , we have Z

K

T dµ∗q =

Z

T (ζ)

(Z

K

∗ Kn,q (ρ; ζ, ξ)

q Y

dφk

k=1

)

∗ dνn,q (ρ; ζ).

(3.3)

Proof. Using the discrete reproduction formula (2.13) with each coordinate separately, we deduce that T (ξ)

q Y

k−1

sin

(φk ) =

k=1

Z

Sqρ

∗ ∗ T (ζ)Kn,q (ρ; ζ, ξ)dνn,q (ρ; ζ).

The formula (3.3) follows by an application of Fubini’s theorem.

2

We observe that formula (3.3) is a quadrature formula for integration on K, based on nodes located in Sqρ . In particular, if q = 2, then ∗ ∗ ∗ dνn,2 (ρ; ζ) = sin(θ2 )dνn+2 (0, ρ; θ2 )dνn+1 (−π, π; θ1 ). ∗ The quadrature weights involved in this formula are given for each ξ in the support of dνn,q by

Z

K

∗ Kn,q (ρ; ζ, ξ)

q Y

dφk .

(3.4)

k=1

Clearly, an explicit computation of these weights will depend upon the geometry of K. In the next section, we do this for a spherical triangle in S2 .

5

4

Integration on a spherical triangle

In this section, we focus on evaluating the quadrature weights (3.4) in the case when K is a spherical triangle ∆∗ in S2 (Fig. 1). It is no longer necessary to have any restriction on the placement of the points involved in these formulas, since we are only evaluating the integrals of specific polynomials. Since the measure µ∗q is rotation invariant, we may assume a standardization of the triangle ∆∗ , which we now describe. First, we observe that the triangle ∆∗ can be decomposed into three disjoint triangles, each with one vertex at the centroid x0 of ∆∗ . Let ρ ∈ (0, π) be chosen so that the cap circumscribing ∆∗ is S2ρ (x0 ). Since a quadrature formula on ∆∗ can then be computed as the sum of the corresponding formulas on each of these disjoint triangles, in our discussion we may replace the original ∆∗ by one of these component subtriangles. By appropriate rotations, we may then assume that the vertices of ∆∗ are given by x0 = n = (0, 0, 1), a = (sin ρ cos α, sin ρ sin α, cos ρ), b = (sin ρ cos α, − sin ρ sin α, cos ρ)

(4.1)

for some α ∈ (0, π) (cf. Figs. 2, 3). (In the context of the original triangle, different rotations are needed for the different components to accomplish this normalization for each component triangle.) By further subdivisions of ∆∗ , it is not a loss of generality to assume that ρ, α ∈ (0, π/2).

n 1.0

p1

0.5

0.0

z

p2

m

- 0.5

-1.0 1.0

p3

1.0 0.5

0.5 0.0

0.0

y

x - 0.5

- 0.5 -1.0

-1.0

Figure 1: An arbitrary spherical triangle and its circumcircle with center m.

In this section, it is convenient to write the parametric representation of S2 in the notation (x, y, z) = (sin θ cos φ, sin θ sin φ, cos θ) =: p(θ, φ),

θ ∈ [0, π], φ ∈ [−π, π],

i.e., instead of the notation introduced in (3.1) write φ in place of θ1 , θ in place of θ2 , and (x, y, z) in place of (x1 , x2 , x3 ). Since the equation of the plane passing through the origin of R3 and the points a, b is given by x cos ρ = z sin ρ cos α, i.e., tan θ cos φ = tan ρ cos α =: L, we deduce that for any integrable f : ∆∗ → C, Z

∆∗

f (θ, φ)dθdφ =

Z

α

−α

Z

arctan(L sec φ)

f (θ, φ)dθdφ.

0

In our new notation, the reproducing kernel in (3.2) takes the form ∗ ˜ φ), ˜ p(θ, φ)) = K ∗ (0, ρ; θ, ˜ θ)K ∗ (−π, π; φ, ˜ φ) Kn,2 (ρ; p(θ, n+1 n

6

1.0

1.0

1.0

1.0

0.5

0.5

0.0 z - 0.5 -1.0

0.0 z - 0.5 -1.0

a

0.5

0.5

c

a

c

y 0.0

y 0.0

n

n

b b

- 0.5

- 0.5

-1.0

-1.0 - 0.5

-1.0

0.0

0.5

1.0

- 0.5

-1.0

x

0.0

0.5

1.0

x

Figure 2: The north pole n is element of the spherical triangle ∆∗ . Decomposition of ∆∗ into three disjoint triangles.

Figure 3: The north pole n lies outside the spherical triangle ∆∗ . Decomposition of ∆∗ into three disjoint triangles.

and we need to evaluate (cf. (3.4)) Z

∆∗

∗ ˜ θ)K ∗ (−π, π; φ, ˜ φ)dθdφ = Kn+1 (0, ρ; θ, n

Z

α

−α

˜ φ) Kn∗ (−π, π; φ,

(Z

arctan(L sec φ)

0

∗ ˜ θ)dθ Kn+1 (0, ρ; θ,

)

dφ.

∗ ˜ ◦) is a trigonometric polynomial of degree at most n+ 1 with For each θ˜ ∈ [0, π], the function Kn+1 (0, ρ; θ, ˜ ˜ ◦). Therefore, coefficients dependent on θ, and a similar observation holds for the function Kn∗ (−π, π; φ, it is sufficient to develop a quadrature formula to compute an integral of the form

Z

α

−α

P (φ)

Z

arctan(L sec φ)

Q(θ)dθdφ

(4.2)

0

for arbitrary P ∈ Hn and Q ∈ Hn+1 . The following proposition summarizes the nature of the inner integral. Proposition 4.1 Let n ≥ 1 be an integer, m0 , m1 be the integer parts of (n+1)/2, (n+2)/2, respectively, and X iℓθ ˆ Q(ℓ)e , θ ∈ [−π, π]. Q(θ) = |ℓ|≤n+1

Then there exist trigonometric polynomials Q0 ∈ H2m0 and Q1 ∈ H2m1 −1 such that, for φ ∈ [0, π/2), Z

arctan(L sec φ)

ˆ arctan(L sec φ) Q(θ)dθ = Q(0)

0

Z arctan(L sec φ) ˆ Q(θ) − Q(θ + π) Q(θ) + Q(θ + π) − 2Q(0) dθ + dθ + 2 2 0 0 Z π 1 Q(θ) − Q(θ + π) ˆ Q(0) arctan(L sec φ) + dθ 2 0 2 Q1 (φ) Q0 (φ) + 2 . + 2 2 m 0 (L + cos φ) (L + cos2 φ)m1 −1/2 Z

=

arctan(L sec φ)

7

(4.3)

Remark. We note that if n is even, then Q0 ∈ Hn , Q1 ∈ Hn+1 , while if n is odd, then Q0 ∈ Hn+1 and Q1 ∈ Hn . Proof of Proposition 4.1. In this proof only, we write X iℓθ ˆ Q(θ) = Q(ℓ)e , ℓ∈Z

ˆ where Q(ℓ) = 0 if |ℓ| > n + 1. Then Q(θ)

ˆ = Q(0) +

2kiθ ˆ Q(2k)e +

X

ˆ Q(2k − 1)e(2k−1)iθ

X k∈Z

k∈Z\{0}

ˆ Q(θ) − Q(θ + π) Q(θ) + Q(θ + π) − 2Q(0) ˆ + . = Q(0) + 2 2 Therefore, Z arctan(L sec φ)

Q(θ)dθ

ˆ arctan(L sec φ) + Q(0)

=

0

ˆ Q(2k) (exp (2ki arctan(L sec φ)) − 1) 2ki

X

k∈Z\{0}

+

X Q(2k ˆ − 1) k∈Z

It is easy to verify that Z So, (4.4) implies that Z arctan(L sec φ)

π 0

Q(θ)dθ

i(2k − 1)

exp ((2k − 1)i arctan(L sec φ)) −

X Q(2k ˆ − 1) k∈Z

i(2k − 1)

. (4.4)

X Q(2k ˆ − 1) Q(θ) − Q(θ + π) dθ = −2 . 2 i(2k − 1) k∈Z

=

ˆ arctan(L sec φ) + Q(0)

Z

π

0

0

+

X

k∈Z\{0}

+ =

ˆ Q(2k) (exp (2ki arctan(L sec φ)) − 1) 2ki

X Q(2k ˆ − 1) k∈Z

Q(θ) − Q(θ + π) dθ 4

i(2k − 1)

exp ((2k − 1)i arctan(L sec φ))

ˆ arctan(L sec φ) + Q(0)

Z

0

π

Q(θ) − Q(θ + π) dθ + Se + So . 4

(4.5)

It is elementary to check that exp(i arctan(L/ cos φ)) =

(cos φ + iL) . (L2 + cos2 φ)1/2

ˆ Since Q(ℓ) = 0 if |ℓ| > n + 1, we deduce that Se

=

m0 m0 X X ˆ ˆ Q(2k) Q(2k) (cos φ + iL)2k − 2 2 k 2ik (L + cos φ) 2ik k=−m k=−m 0

0

k6=0

=

k6=0

m0 ˆ X Q(2k) (cos φ + iL)2k (L2 + cos2 φ)m0 −k

k=1

(L2 + cos2 φ)m0

2ik

−

m0 X ˆ (L2 + cos2 φ)m0 Q(2k) 2ik (L2 + cos2 φ)m0 k=−m

+

m0 ˆ X Q(−2k) (cos φ − iL)2k (L2 + cos2 φ)m0 −k k=1

(−2ik)

(L2 + cos2 φ)m0

0

k6=0

=

Q0 (φ) , (L2 + cos2 φ)m0

(4.6) 8

where Q0 ∈ H2m0 . Similarly, So

=

m1 X

k=−m1 +1

=

ˆ Q(2k − 1) (cos φ + iL)2k−1 (2k − 1)i (L2 + cos2 φ)k−1/2

m1 ˆ X Q(2k − 1) (cos φ + iL)2k−1 (L2 + cos2 φ)m1 −k (2k − 1)i (L2 + cos2 φ)m1 −1/2 k=1

+

m1 ˆ X Q(−2k + 1) (cos φ − iL)2k−1 (L2 + cos2 φ)m1 −k k=1

=

(L2

(L2 + cos2 φ)m1 −1/2

(−2k + 1)i

Q1 (φ) , + cos2 φ)m1 −1/2

where Q1 ∈ H2m1 −1 . Together with (4.5) and (4.6), this completes the proof.

2

In view of Proposition 4.1, the evaluation of the integral in (4.2) involves the evaluation of integrals with four weight functions. Let  (1 − t2 )−1/2 , if t ∈ [cos α, 1], W (α; t) = (4.7) 0, otherwise, W0 (ρ, α; t) :=

W1 (ρ, α; t) := Ws (ρ, α; t) :=

(L2 + t2 )−m0 W (α; t), 2

2 −m1 +1/2

(L + t ) W (α; t), arctan(L/t)W (α; t).

(4.8)

(4.9) (4.10)

The evaluation of the integral in (4.2) in terms of the values of the polynomials P and Q is summarized in the following theorem. Theorem 4.1 Let n ≥ 1 be an integer, mp be the integer part of (n + 1 + p)/2, p = 0, 1, 2, P ∈ Hn , and Q ∈ Hn+1 . We use the notation (cf. (2.10)) ∗ ∗ ∗ = νn+1 (W0 ), σn,1 = νn+1 (W1 ), = νm1 (Ws ), σn,0 σn∗ = νm1 (W ), σn,s

  n+1 1 X 2kπ C0 := , Q n+2 n+2 k=0

and (cf. (2.11)) C1 :=

Z

Q(θ) − Q(θ + π) ∗ dνm2 (0, π; θ). 4

Let C0 , C1 be the supports of σ0∗ , σ1∗ , respectively. Then Z α Z arctan(L sec φ) Z ∗ P (φ)Q(θ)dθdφ = C0 P (φ)dσn,s (φ) −α 0 Z +C1 P (φ)dσn∗ (φ) Z Z Q(θ) + Q(θ + π) − 2C0 ∗ ∗ (φ) dνm2 (0, arctan(L/ cos φ); θ)dσn,0 + P (φ)(L2 + cos2 φ)m0 2 C0 Z Z + P (φ)(L2 + cos2 φ)m1 −1/2 C1   Q(θ) − Q(θ + π) ∗ ∗ (0, arctan(L/ cos φ); θ)dσn,1 (φ). − C1 dνm × 2 2

9

(4.11)

Proof. We evaluate the inner integral first, using Proposition 4.1. Since Q ∈ Hn+1 , the equation ˆ Q(0) =

  n+1 1 X 2kπ = C0 Q n+2 n+2

(4.12)

k=0

follows from the orthogonality relations    n+1 1 X 2kjπi 1, = exp 0, n+2 n+2 k=0

The equation Z

0

π

if j = 0, if 0 < |j| ≤ n + 1.

Q(θ) − Q(θ + π) dθ = C1 4

(4.13)

follows from the quadrature formula (2.12). Hence, (4.3) shows that Z

α

P (φ)

−α

Z

arctan(L sec φ)

Q(θ)dθdφ 0

Z α ˆ = Q(0) P (φ) arctan(L sec φ)dφ −α Z π Z α Q(θ) − Q(θ + π) + dθ P (φ)dφ 4 0 −α Z α Z α P (φ)Q0 (φ) P (φ)Q1 (φ) dφ. + dφ + 2 2 m0 2 2 m1 −1/2 −α (L + cos φ) −α (L + cos φ)

It is easy to verify that the integer parts of (n + 2m0 + 2)/2 and (n + 2m1 + 1)/2 are both equal to n + 1. Since P Q0 ∈ Hn+2m0 and P Q1 ∈ Hn+2m1 −1 , we may use the quadrature formula (2.12), (4.12), and (4.13) to complete the proof of (4.11). 2

5

Algorithm

In this chapter, we summarize the construction of the quadrature formula described above on a spherical triangle ∆∗ with the edges in (4.1) in the R form of an algorithm. For any T ∈ H2n , we will compute ∆∗ T dµ∗2 using Theorem 3.1. Using the notation in Section 4, formula (3.3) takes the form (Z (Z ) ) Z Z α arctan(L sec φ) ∗ ∗ ∗ ˜ φ)) ˜ ˜ φ) ˜ φ)), ˜ ˜ θ)dθ dφ dν ∗ (p(θ, T (p(θ, T dµ = K (−π, π; φ, K (0, ρ; θ, ∆∗

2

C

−α

n

0

n+1

n,2

(5.1) ∗ ˜ φ) ˜ ∈ C, we will compute the inner where C ⊂ S2ρ is the support of the discrete measure νn,2 . For each p(θ, integral on the right hand side of (5.1) using Theorem 4.1. All the quadrature measures involved in the use of all these formulas are computed using the method of discretization of the weight functions which ∗ Gautschi describes in [5] to [7]. To evaluate the kernel function Kn+1 (0, ρ; ◦, ◦) in the corresponding nodes we use the three term recurrence relation of the orthogonal polynomials (see (2.2)) with respect to the weight function w(0, ρ; ◦) in (2.4) and w ˜ = (1 − ◦2 )w(0, ρ; ◦). The Gautschi algorithm yields approximations of the recurrence coefficients. The kernel Kn∗ (−π, π; ◦, ◦) is associated to the normalized Chebyshev polynomials of the first and second kind. Algorithm Assumption : The triangle is given by the vertices as in (4.1). Input : Polynomial degree n, radius ρ, angle α. Calculate the weights and the nodes from the Gaussian quadratures with the following weight functions, where N denotes the number of the weights and nodes (cf. (2.8) and (2.9)): 10

(1)

:=

(2)

:=

(3)

:=

(1)

:= χj,N (W0 (ρ, α)), λj

(2)

:= χj,N (W1 (ρ, α)), λj

(3)

:= χj,N (Ws (ρ, α)), λj

1. With weight function W0 (ρ, α) as in (4.8) and N = 2(n + 1), find χj λj,N (W0 (ρ, α)). 2. With weight function W1 (ρ, α) as in (4.9) and N = 2(n + 1), find χj λj,N (W1 (ρ, α)).

3. With weight function Ws (ρ, α) as in (4.10) and N = 2⌊ n+2 2 ⌋, find χj λj,N (Ws (ρ, α)).

(4)

4. With weight function W (α) as in (4.7) and N = 2⌊ n+2 2 ⌋, find χj λj,N (W (α)).

(4)

(5)

(6)

(5)

=

(6)

=

= χj,N (w(0, ρ)), λj

5. With weight function w(0, ρ) as in (2.4) and N = 2(n + 2), find χj λj,N (w(0, ρ)). 6. With weight function w(0, π) as in (2.4) and N = 2⌊ n+3 2 ⌋, find χj λj,N (w(0, π)).

:=

:= χj,N (W (α)), λj

= χj,N (w(0, π)), λj

(1)

(7)

(2)

(8)

7. With weight function w1 = w(0, arctan(L/ cos χj )) as in (2.4) and N = 2⌊ n+3 2 ⌋, find χj,k = (7)

(1)

χj,k,N (w1 ), λj,k = λj,k,N (w1 ) for every node χj .

8. With weight function w2 = w(0, arctan(L/ cos χj )) as in (2.4) and N = 2⌊ n+3 2 ⌋, find χj,k = (8)

(2)

χj,k,N (w2 ), λj,k = λj,k,N (w2 ) for every node χj . (9)

(9)

(9)

(9)

9. For s = 1, 2, . . . , n + 1, set χs = −χ−s = arccos (2s−1)π 2(n+1) , λs = λ−s =

π n+1 .

Evaluate C0 (χ(5) r ) =

  n+1 2πk 1 X ∗ , , Kn+1 0, ρ; χ(5) r n+2 n+2 k=0

⌊ n+3 2 ⌋

(5) (6) ∗ (6) Kn+1 (0, ρ; χr , χk )

X

C1 (χ(5) r ) =

λk

k=−⌊ n+3 2 ⌋,k6=0

(5)

(6)

∗ − Kn+1 (0, ρ; χr , χk + π) , 4

and the polynomials (1)

Q0 (χ(5) r , χj ) =

(1)

(L2 + cos2 χj )⌊(n+1)/2⌋ ⌊ n+3 2 ⌋

× (2)

Q1 (χ(5) r , χj ) =

(5) (7) ∗ (7) Kn+1 (0, ρ; χr , χj,k )

X

λj,k

(5)

(7)

(5)

∗ + Kn+1 (0, ρ; χr , χj,k + π) − 2C0 (χr )

2

k=−⌊ n+3 2 ⌋,k6=0

,

(2)

(L2 + cos2 χj )⌊(n+2)/2⌋−1/2 ⌊ n+3 2 ⌋

×

X

(5) (8) ∗ (8) Kn+1 (0, ρ; χr , χj,k ) λj,k

k=−⌊ n+3 2 ⌋,k6=0

(5)

(8)

∗ − Kn+1 (0, ρ; χr , χj,k + π)

2

and the summands in (4.11) ⌊ n+2 2 ⌋ (9) S1 (χ(5) r , χs ))

=

C0 (χ(5) r )

X

j=−⌊ n+2 2 ⌋,j6=0 ⌊ n+2 2 ⌋ (9) S2 (χ(5) r , χs )) =

C1 (χ(5) r )

X

j=−⌊ n+2 2 ⌋,j6=0

11

  (3) (3) , , χ λj Kn∗ −π, π; χ(9) s j   (4) (4) , λj Kn∗ −π, π; χ(9) s , χj

!

− C1 (χ(5) r ) ,

n+1 X

(9) S3 (χ(5) r , χs )) =

j=−(n+1),j6=0 (9) S4 (χ(5) r , χs )) =

n+1 X

j=−(n+1),j6=0

  (1) (1) (1) (9) ∗ , −π, π; χ , χ )K , χ λj Q0 (χ(5) s n r j j   (2) (2) (2) ∗ (9) . λj Q1 (χ(5) r , χj )Kn −π, π; χs , χj

Output: The quadrature weights are (9) (5) (9) (5) (9) (5) (9) (5) (9) wr,s = λ(5) r λs (S1 (χr , χs ) + S2 (χr , χs ) + S3 (χr , χs ) + S4 (χr , χs )),

(5.2)

the nodes are (9) (χ(5) r , χs ), r = −(n + 2), . . . , −1, 1, . . . , n + 2, s = −(n + 1) . . . , −1, 1, . . . , n + 1.

6

(5.3)

Numerical examples

We model the surface of the earth by the sphere, and construct a triangulation T using the algorithm given by Renka [13], i.e., by choosing the geographical centers of the countries as the vertices of the triangles involved. This yields 476 spherical triangles ∆k , k = 1, 2, . . . , 476, which have different geometrical properties. Figure 4 shows some of the different kinds of spherical triangles and their radii ρk thus generated.

log10 Ρ k 0.0

100

200

300

400

D -Index k

-0.5

-1.0

-1.5

-2.0

Figure 4: Triangulation of the surface of the earth modeled as the sphere. The vertices of the triangles ∆k , k = 1, 2, . . . , 476, are the geographical centers of the countries. The triangles have different geometrical properties, for example different radii ρk .

Decomposing into isosceles spherical triangles, calculating the weights in (5.2) and the nodes in (5.3) for each subtriangle numerically as well as summarizing the weights in a suitable way lead to the weights (w(n,k) )j := wn,k,j and the nodes xn,k,j , j = 1, 2, . . . , 4(n + 2)(n + 1), of a triangle ∆k , k = 1, 2, . . . , 476, given the polynomial degree n. Unfortunately, the results of the Gautschi algorithm and the following calculation of the sum in (5.2) are not stable. But we can handle this cancellation problem by applying the computer algebra system Mathematica which can deal with numbers of high precision. On the one 12

hand, this decreases the speed of calculation. On the other hand, the precision of the weights in (5.2) and the nodes in (5.3) can be arbitrarily high, depending on the input data. If we merge all the weights R and nodes of the√ triangles in T , we will obtain a quadrature rule on the sphere. We observe that S2 Yl (x)dµ∗2 (x) = 4πδl,0 . To show that the computed quadrature weights and nodes for a given n are really exact for polynomials of degree at most n, we compose the matrices (Y(r,k) )l,j = Yl (xn,k,j ), l = 1, 2, . . . , (r + 1)2 , j = 1, 2, . . . , 4(n + 2)(n + 1), k = 1, 2, . . . , 476, and the (r) vectors e1 := (1, 0, . . . , 0)T of length (r + 1)2 . We calculate the error values √ P476 (r) k k=1 Y(r,k) w(n,k) − 4πe1 k2 (n,r) , r = 0, 1, . . . , 33, E2 := √ (r) k 4πe1 k2 n = 5, 10, 15, 20, 25, 30. In Figure 5 the error values are plotted for the degrees from n = 5 up to n = 30 and prove the expected exactness of the quadrature. The number of digits lost depends on the polynomial degree n and the geometrical property of the spherical triangle for which the quadrature is calculated. Let p(n,k) be the precision of the weights vector w(n,k) , i.e., the smallest value of precision of the entries. The entries of 476 1 X (n,k) p is nearly linear in the degree n. Table 1 show that the average loss of precision 476 k=1

H n ,rL

log10 E 2 0

5

10

15

20

25

30

35

r

-5 n =30

-10

n =25 -15 n =20 -20

n =15

-25

n =10

-30

n =5

-35 (n,r)

Figure 5: Error values log10 E2 of the quadrature rule on the sphere with given polynomial degree n = 5, 10, 15, 20, 25, 30.

n Loss of Precision

5 40.15

10 71.34

Table 1: Average loss of precision 1/476

7

15 104.07

P476

k=1

20 132.10

25 159.92

30 185.20

p(n,k) depending on polynomial degree n.

Conclusions

Many applications in geophysics involve approximation of functions. In [9], we have developed a technique to accomplish this in a localized manner, yielding results far superior to the traditional least square fit in the case of several benchmark functions. In turn, our procedure involves the computation of integrals of functions over the whole sphere. In the spirit of localized approximation, it is desirable to be able to compute the integrals over subsets of the sphere, rather than the whole sphere. In this paper, we have described a method to compute such integrals exact for high degree polynomials on compact subsets of a 13

sphere embedded in a Euclidean space Rq+1 , q ≥ 2. Specializing to the case of integration over spherical triangles on S2 , this method is developed further into a specific algorithm. We have tested this algorithm on triangles of various shapes, based on data on the surface of the sphere. Further work in progress uses our algorithm to construct similar quadrature formulas for integration based on values of the integrand at scattered sites; i.e., in situations where one has no control on the sites where the integrand needs to be evaluated.

Acknowledgement We would like to thank Larry Schumaker for proposing the problem and many subsequent fruitful discussions. Moreover, we are grateful to the referees for their many helpful comments.

References [1] W. Freeden, O. Glockner, and M. Schreiner, Spherical Panel Clustering and Its Numerical Aspects, AGTM Report No. 183, University of Kaiserlautern, Geomathematics Group, 1997. [2] W. Freeden and U. Windheuser, Spherical Wavelet Transform and Its Discretization, AGTM Report No. 125, University of Kaiserlautern, Geomathematics Group, 1995. [3] G. Freud, “Orthogonal Polynomials”, Acad´emiai Kiado, Budapest, 1971. [4] M. Ganesh and H. N. Mhaskar, Matrix–free interpolation on the sphere, SIAM J. Numer. Analysis 44 (3) (2006), pp. 1314–1331. [5] W. Gautschi, On Generating Orthogonal Polynomials, SIAM J. Sci. Stat. Comput. 3 (3) (1982), 289–317. [6] W. Gautschi, Orthogonal Polynomials: Computation and Approximation, Oxford University Press, Oxford, 2004. [7] W. Gautschi, Orthogonal polynomials (in Matlab), J. Comput. Appl. Math. 178 (1–2) (2005), 215– 234. [8] K. Hesse and R. S. Womersley, Numerical integration with polynomial exactness over a spherical cap, Adv. in Comp. Math., available online Sept. 23, 2011. [9] Q. T. Le Gia and H. N. Mhaskar, Localized linear polynomial operators and quadrature formulas on the sphere, SIAM J. Numer. Anal. 47 (1) (2008), 440–466. [10] H. N. Mhaskar, F. J. Narcowich, and J. D. Ward, Spherical Marcinkiewicz-Zygmund Inequalities and positive quadrature, Math. Comp. 70 (2001), 1113-1130. [11] H. N. Mhaskar, Local quadrature formulas on the sphere, J. Complexity, 20 (2004), 753–772. [12] H. N. Mhaskar, Local quadrature formulas on the sphere, II, in “Advances in Constructive Approximation” (M. Neamtu and E. B. Saff eds), Nashboro Press, Nashville, 2004, pp. 333–344. [13] R. J. Renka, Algorithm 772: STRIPACK: Delaunay Triangulation and Voronoi Diagram on the Surface of a Sphere, ACM Trans. Math. Softw. 23 (3) (1997), 416-434.

14