Computing with Functions in Spherical and Polar Geometries I. The Sphere

Abstract

Abstract. A collection of algorithms is described for numerically computing with smooth functions defined on the unit sphere. Functions are approximated to essentially machine precision by using a structure-preserving iterative variant of Gaussian elimination together with the double Fourier sphere method. We show that this procedure allows for stable differentiation, reduces the oversampling of functions near the poles, and converges for certain analytic functions. Operations such as function evaluation, differentiation, and integration are particularly efficient and can be computed by essentially one-dimensional algorithms. A highlight is an optimal complexity direct solver for Poisson's equation on the sphere using a spectral method. Without parallelization, we solve Poisson's equation with 100 million degrees of freedom in 1 minute on a standard laptop. Numerical results are presented throughout. In a companion paper (part II) we extend the ideas presented here to computing with functions on the disk.Key words. low rank approximation, Gaussian elimination, functions, approximation theory AMS subject classification. 65D05 DOI. 10.1137/15M10458551. Introduction. Spherical geometries are universal in computational science and engineering, arising in weather forecasting and climate modeling [11,13,17,24,28,31,35], geophysics [16,48], and astrophysics [3,8,39]. At various levels these applications all require the approximation of functions defined on the surface of the unit sphere. For such computational tasks, a standard approach is to use longitude-latitude coordinates (λ, θ) ∈ [−π, π] × [0, π], where λ and θ denote the azimuthal and polar angles, respectively. Thus, computations with functions on the sphere can be conveniently related to analogous tasks involving functions defined on a rectangular domain. This is a useful observation that, unfortunately, also has many severe disadvantages due to artificial pole singularities introduced by the coordinate transform.In this paper, we synthesize a classic technique known as the double Fourier sphere (DFS) method [6,18,28,31,50] together with new algorithmic techniques in low rank function approximation [5,42]. This alleviates many of the drawbacks inherent with standard coordinate transforms. Our approximants have several attractive properties: (1) no artificial pole singularities, (2) a representation that allows for fast algorithms, (3) a structure so that differentiation is stable, and (4) an underlying interpolation grid that rarely oversamples functions near the poles.

Rhodonea Curves as Sampling Trajectories for Spectral Interpolation on the Unit Disk

“…The reflection symmetry corresponds to the invariance of g under the reflection operator i † on J (m) , while the glide reflection symmetry is described by the operator i * on the subset K (m) . In [33,35], this glide reflection symmetry is referred to as block-mirror centrosymmetric (BMC) structure. determined by a discrete Fourier transform.…”

Tensor calculus in spherical coordinates using Jacobi polynomials. Part-I: Mathematical analysis and derivations

“…Several authors have employed Chebyshev polynomials as a basis for Reg(ℓ) in spherical problems; see [27,3,17,29,50]. It is possible to use Chebyshev polynomials and achieve reasonably good results.…”

“…The unit ball has three separate coordinate singularities: disk-like singularities at the north and south pole, and a third at the centre of the domain. Spherical geometry is important for a large number of two-and three-dimensional applications; see e.g., [1,46,50,43,42,48,29,17,22]. Astrophysical and planetary applications provide many obvious examples with stars and planets.…”

A spectral interpolation scheme on the unit sphere based on the nodes of spherical Lissajous curves

“…. Now, based on the formulas (18), (19) and (20) as well as Proposition 5, the statements about the orthogonality and the norms of the basis functions can be derived similarly as in Theorem 6. As a template for the entire procedure, we calculate the norm χ…”

“…where we used the product formulas (18), (19) and (20) to manipulate the function terms in the integral. Next, we check in which cases the condition (14) given in Proposition 5 is satisfied and determine in this way the value of the norm.…”