Sperical splines


Project description



In general terms, a spline is a real-valued function f on some domain D which is defined by parts: f(p) = fi(p) iff p belongs to Di, where G = {D1, D2, ... Dn} is some fixed partition of D, and each fi is a function defined on Di.

Splines are often used for function approximation (modeling some function that cannot be easily computed, such as the solution of a differential equation); and for data interpolation (modeling a function which is known only at certain scattered points). In such applications, the partial functions fi (the pieces) of f) are easily computed formulas, such as polynomials of a fixed maximum degree g; and the subsets D i (the tiles) have a simple geometry such as intervals, rectangles, triangles, etc.. Also, the spline is often required to be continuous, and to have continuous derivatives up to some order r, even across the boundaries between the tiles Di.

Spherical splines

A spherical function is a real-valued function whose domain is the d-dimensional sphere Sd (the set of all unit-length vectors of Rd+1).

A spherical spline is a piecewise-polynomial function f defined on the sphere Sd, in terms of a fixed triangulation T of the latter. Specifically, within each triangle ti of T, the value f(x) is defined by a polynomial fi(x0,.. xd) whose arguments are the d+1 Cartesian coordinates of the point x.

Why splines?

As alternatives to other "global" approximation and interpolation models --- such as polynomials and trigonometric series --- splines have a number of advantages, such as flexibility, numerical stability, efficient evaluation, analytic tractability, and more. Splines are competitive alternatives to other "local" schemes such as radial bases and wavelets.

Spherical spline spaces

The spaces Prg and Hrg

In this project we study the space Prg of all spherical splines with a given triangulation T and a given maximum degree g within each triangle, which are also continuous and differentiable up to a given order r.

We also study the subspace Hrg of Prg, consisting of the splines whose pieces fi are homogeneous polynomials of degree exactly g. This subspace was exhaustively studied by P. Alfeld, M. Neamtu and L. Schumaker, who provided an algorithm to construct explicit basis with local support (the ANS basis) --- provided that g is at least 2r+1.

Bases for Prg

We have shown [??] that, for dimension d = 2 and higher, the space Prg is the direct sum of Hrg and Hrg-1. Therefore, one can obtain a basis for the former by simply concatenating the ANS bases for the latter two spaces. This construction works as long as g is at least 2r+2.

Least-squares approximation

Once a basis B for Prg has been constructed, we can use it to approximate a given spherical function f by a spline s of Prg, with least-squares error. For that purpose, we need the matrix M, whose generic entry Mij is the scalar product <Bi|Bj> of the basis elements Bi and Bj; and the vector y such that yi = <Bi|f>. The approximation s is then SUMj xj Bj, where x is the solution of the linear equation system M x = y.

The function s thus obtained minimizes the squared norm of the approximation error, ||s-f||2 = <s-f|s-f>. For the scalar product <f|h>, we use the integral of f(p)h(p) over the whole sphere (the L2 metric), optionally multiplied by a weight function w(p).

Since, in our case, each basis element Bi has finite support, confined to a few triangles, the matrix M is quite sparse --- it has only a constant number of nonzero entries per row. Thus the system M x = y can be solved fairly efficiently, even for moderately large triangulations. We have been using a direct solution (via Cholesky factorization of the matrix M) which seems adequate for bases with up to ?? elements. Beyond that limit the factorization is prone to fail; presumably iterative methods like Gauss-Jordan can exceed that limit.

PDE integration

We have also investigated the suitability of the spline spaces Prg and Hrg for integration of partial differential equations (PDEs) on the sphere, such as those that arise in global weather modeling.

The general problem

The generic PDE has the form K(f)(p) = h(p), where f is the unknown function to be determined, h is is a given function independent of f, and K is a known differential operator. h(p). A standard example for the sphere is the Helmholtz equation, (Lf)(p) - Cf(p) = h(p) where L denotes the spherical Laplacian operator and C is a given real constant.

Numerical PDE integration

When solving (integrating) the PDE numerically, we look for an optimal approximation s to the true solution f in some finite-dimensional space S of approximating functions.

In principle, when K is a linear operator, there is an explicit formula for approximate solution s. The computation becomes much easier if we use Galerkin's criterion of optimality -- namely, make the equation's residual K(s)(p) - h(p) orthogonal to the approximation space S.

With this criterion, if B={B1,.. Bn} is a basis for S, the PDE integration problem reduces to solving the linear equation system N x = y where Nij is <K(Bi)|Bj> yi = <Bi|h>, and the unknowns xj are the coefficients of Bj in the approximate solution s.

The finite-element approach

In our case, S is a space of spherical splines, either Prg or Hrg, over some fixed triangulation T of the sphere S2, and the basis B is obtained as discussed above. As in least-squares approximation, the matrix N is quite sparse.

When the basis elements Bi are splines, one must worry about K(Bi) becoming infinite at tile boundaries. Such infinities would cause problems in the computation of <K(Bi)|Bj>. Thus, if K includes first-order derivatives, the elements Bi must be at least continuous (C0).

However, when K is a second-order operator, we can use Green's theorem can to reduce <K(Bi)|Bj> to a combination of <Bi|Bj> and <G(Bi)|G(Bj)> where G is the spherical gradient operator. Thus we can avoid second derivatives, and use C0 instead of C1 splines.

Solving non-linear PDEs

When the operator K is not linear, there is no general direct method; in fact, Galerkin's criterion may not be enough to ensure a unique solution.

If we are lucky, the solution may be found by the following iterative method. We start with an initial guess s := s0 = SUMj xj Bj. At each iteration k > 0, we replace the operator K by a first-order approximation K ??



Last edited on 2003-06-11 01:26:37 by stolfi