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.

http://www.springerlink.com/content/tn78433n5h53505h/

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 SPP–UX 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/