Fast numerical solution of the vector Molodensky problem

When standard boundary element methods (BEM) are used in order to solve the linearized vector Molodensky problem we are confronted with two problems: (i) the absence of O (|x|(-2)) terms in the decay condition is not taken into account, since the single layer ansatz, which is commonly used as representation of the disturbing potential, is of the order O (|x|(-1)) as x --> infinity. This implies that the standard theory of Galerkin BEM is applicable since the injectivity of the integral operator fails; (ii) the N x N stiffness matrix is dense, with N typically of the order 10(5). Without fast algorithms, which provide suitable approximations to the stiffness matrix by a sparse one with O (N . log(s) N), s greater than or equal to 0, non-zero elements, high-resolution global gravity field recovery is not feasible. We propose solutions to both problems. (i) A proper variational formulation taking the decay condition into account is based on some closed subspace of codimension 3 of L-2(T). Instead of imposing the constraints directly on the boundary element trial space, we incorporate them into a variational formulation by penalization with a Lagrange multiplier. The conforming discretization yields an augmented linear system of equations of dimension (N + 3) x (N + 3). The penalty term guarantees the well-posedness of the problem, and gives precise information about the incompatibility of the data. (ii) Since the upper left submatrix of dimension N x N of the augmented system is the stiffness matrix of the standard BEM, the approach allows to use all techniques to generate sparse approximations to the stiffness matrix such as wavelets, fast multipole methods, panel clustering etc. without any modification. We use a combination of panel clustering and fast multipole method in order to solve the augmented linear system of equations in O(N) operations. The method is based on an approximation of the kernel function of the integral operator by a degenerate kernel in the far field, which is provided by a multipole expansion of the kernel function. In order to demonstrate the potential of the method we solve a Robin problem on the sphere with a nullspace of dimension 3. For N = 65538 unknowns the matrix assembly takes about 600 s and the solution of the sparse linear system using GMRES without any preconditioning takes about 8 s. 30 iterations are sufficient to make error smaller than the discretization error.