# Earth Cubed

## Recursive Numeric Integration, The FFT and Spherical Harmonics

I was briefly exploring various methods of numeric integration,  because if any spectral techniques are used (e.g. Spherical Harmonics) then, the computation of there coefficients involves an inner product between the basis function, and the properties of the earth who’s spectral expansion is being computed. The computation of the inner product involves integration.

In order to reduce computations, some kind of recursive integration may be desirable. Examples include: Clenshaw Curtis_Quadratures, Gaussian Quadrature, Romberg method.

These are of course  one dimensional techniques. An interesting concept for multiple dimensional integrals is as follows:

Recursive Integration Methodologies with Statistical Applications
A. J. Hayter, Industrial & Systems Engineering
Georgia Tech, Atlanta, USA
Summary

This paper shows how recursive integration methodologies can be used to evaluate high dimensional integral expressions. This has applications to many areas of statistical inference where probability calculations and critical point evaluations often require such high dimensional integral evaluations. Recursive integration can allow an integral expression of a given dimension to be evaluated by a series of calculations of a smaller dimension. This significantly reduces the computation time. The application of the recursive integration methodology is illustrated with several examples.

http://www2.isye.gatech.edu/statistics/papers/04-27.pdf

As for the single dimension methods of integration, Gaussian quadrature works well, when the  function being ingegrated is well approximated by a polynomial.  Romberg is quite interesting, it will work for uneven step sizes. It is derived by applying Richardson extrapolation (Richardson 1910) repeatedly on the trapezium rule.

Clenshaw Curtis_quadratures is particularly interesting. In this method a change of variables is done by letting x=cos(theata). This reduces the problem of integration to the problem of finding the Fourier series, which can be efficiently by using the Fast Fourier Transform. The Fast Fourier Transform has the advantage that the computational complexity grows with the logarithm of the number of the points.

In climatology spherical harmonics are used to represent the spectrum instead of Fourier Transforms. It is worth checking to see if there are efficient ways to compute these transforms. Some relevant links include:

Abstract
Spherical Harmonic Transforms (SHTs) which are essentially Fourier transforms on the sphere are critical in global geopotential and related applications. Among the best known strategies for discrete SHTs are Chebychev quadratures and least squares. The numerical evaluation of the Legendre functions are especially challenging for very high degrees and orders which are required for advanced geocomputations. The computational aspects of SHTs and their inverses using both quadrature and least-squares estimation methods are discussed with special emphasis on numerical preconditioning that guarantees reliable results for degrees and orders up to 3800 in REAL*8 or double precision arithmetic. These numerical results of spherical harmonic synthesis and analysis using simulated spectral coefficients are new and especially important for a number of geodetic, geophysical and related applications with ground resolutions approaching 5 km.

and:

SpharmonicKit is a collection of routines, written in C, which implement discrete Legendre and spherical harmonic transforms by a number of different algorithms. For certain algorithms, code for the inverse transform is also provided. Included as well are routines for spherical convolutions.The Legendre transforms in the Kit require that a function of bandwidth B be sampled on the 2B-many Chebyshev points and not the B-many Gaussian points. Therefore, the number of samples needed is twice the function’s bandwidth. This implies that for the spherical transforms in the Kit, a function of bandwidth B is sampled on the equiangular 2Bx2B grid on the sphere.

Some of the algorithms implemented in this package are based on the work of Driscoll and Healy. A preprint that offers a detailed development and description of the algorithms implemented in SpharmonicKit is here.

Caveat emptor – this is research code only and has not been hardened. All the code works quite well on DEC Alpha workstations using OSF1 Versions 3.2 and 4.0. Some code has also been tested and successfully run on Pentium-based GNU/Linux workstations, SGI workstations using IRIX 5.3 and IRIX 6.4, Sun workstations using SunOS 4.1.2 1, an HP Exemplar X-Class running SPPUX 5.2.1, and even a NeXTstation Color Turbo running NeXTStep 3.3!

We would also like to mention that there exists other software for performing spherical transforms. Details on other possibilities may be found here.

http://www.cs.dartmouth.edu/~geelong/sphere/

August 17, 2009 - Posted by | Math

1. It appears, that the spherical harmonic kit, gives you the option to install a modified fftpack which will use the fast Fourier transform to speed up computations in the toolkit.

http://www.cs.dartmouth.edu/~geelong/sphere/modified_fftpack.html

Comment by s243a | August 17, 2009 | Reply

2. I just noticed the following at the bottom of the wikipedia page on the FFT:

“Other generalizations

An O(N5/2 log N) generalization to spherical harmonics on the sphere S2 with N2 nodes was described by Mohlenkamp (1999), along with an algorithm conjectured (but not proven) to have O(N2 log2 N) complexity; Mohlenkamp also provides an implementation in the libftsh library. A spherical-harmonic algorithm with O(N2 log N) complexity is described by Rokhlin and Tygert (2006).

Various groups have also published “FFT” algorithms for non-equispaced data, as reviewed in Potts et al. (2001). Such algorithms do not strictly compute the DFT (which is only defined for equispaced data), but rather some approximation thereof (a Non-uniform discrete Fourier transform, or NDFT, which itself is often computed only approximately).”
http://en.wikipedia.org/wiki/Fast_Fourier_transform#Other_generalizations

Comment by s243a | August 17, 2009 | Reply

3. Searching for “FFT Spherical harmonics also seems to yield some results” (perhaps not as good):

“Abstract. Spherical Harmonic Transforms (SHTs) which are essentially Fourier transforms on the sphere are critical in global geopotential and related applications. Among the best known strategies for discrete SHTs are Chebychev quadratures and least squares. The numerical evaluation of the Legendre functions are especially challenging for very high degrees and orders which are required for advanced geocomputations. The computational aspects of SHTs and their
inverses using both quadrature and least-squares estimation methods are discussed with special emphasis on numerical preconditioning that guarantees reliable results for degrees and orders up to 3800 in REAL*8 or double precision
arithmetic. These numerical results of spherical harmonic synthesis and analysis using simulated spectral coefficients are new and especially important for a number of geodetic, geophysical and related applications with ground resolutions approaching 5 km.”
http://www.ucalgary.ca/~blais/Blais2008b.pdf

Comment by s243a | August 17, 2009 | Reply

4. The following paper derives Romberg method.

http://www.stanford.edu/~lambers/math105a/lecture13.pdf

Comment by s243a | August 18, 2009 | Reply

5. Here is an interesting PDF describing the connection between polynomials and spherical harmonics.

http://www.regonaudio.com/SphericalHarmonics.pdf

Basically, you can represent the spherical harmonics you can to a transformation, between spherical harmonics and polynomials which are linearly independent and orthogonal over the sphere.

Here is a document covering a lot about orthogonal polynomials.

http://www.emis.de/journals/SAT/papers/3/3.pdf

Comment by s243a | August 18, 2009 | Reply

6. I find the following aslo interesting to this topic:

“spherical harmonics and related Fourier transforms
Agata BEZUBIK † ∗, Agata DĄBROWSKA † ∗∗ and Aleksander STRASBURGER ‡

Abstract
This article summarizes a new, direct approach to the determination of the expansion into spherical harmonics of the exponential ei(x|y) with x, y ∈ Rd. It is elementary in
the sense that it is based on direct computations with the canonical decomposition of homogeneous polynomials into harmonic components and avoids using any integral identities. The proof makes also use of the standard representation theoretic properties of spherical harmonics and the explicit form of the reproducing kernels for these spaces by means of classical Gegenbauer polynomials. In the last section of the paper a new method of computing the Fourier transforms of SO(d)-finite functions on the unit
sphere is presented which enables us to reobtain both the classical Bochner identity and generalizations of it due to one of the present authors and F. J. Gonzalez Vieli.”