Report 1999-16

Fast numerical solution of the linearized Molodensky problem

R. Klees, M. van Gelderen, C. Lage and C. Schwab

Abstract: When standard boundary element methods (BEM) are used 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 perturbation potential, is of the order $O(|x|^{-1})$ as $x \to \infty$. This implies that the standard theory of Galerkin BEM is not applicable since the injectivity of the integral operator fails; (ii) the Nx N BEM 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 N)^s),\; s\ge 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 co-dimension 3 of $\LG$. 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. 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 on a workstation. 30 GMRES iterations are sufficient to make the error smaller than the discretization error.

Paper: Available as PDF (3,3 MB) or as hardcopy to order reports@sam.math.ethz.ch.