Extension to 3D

Because the Laplacian of a function is an invariant, the method should work the same way in 3D. Computing the Laplacian of a function of two arguments in the nodes of an unstructured mesh is an interesting mathematical problem. The simplest approach, analoguous to the one used for Cartesian grids, is to interpolate a surface through local sets of points and to compute the Laplacian analytically from its coefficients:

 

z = a + bx + cy + dx² + exy + fy² + gx² y²,

 

 450#450 (183)

Since the surface has seven terms, we must know the function z=f(x,y) at seven points in order to find the Laplacian at one of the seven. This means solving a 451#451 linear system to find the coefficients a, b, c, d, e, f and g, given three pairs of (x, y, z) points. It seems frustrating that we have to compute all the seven coefficients while using only the values of d, f, and g. Fortunately, we do not have to solve innumerable amounts of 451#451 systems for each downward continuation step: the x and y values depend only on the geometry of the spatial mesh. For each point, the matrix inversion can be done only once: in the beginning, and after that for each point. We only need to multiply the vector of z values with three rows of the precomputed inverted matrix in order to find out the values of d, f, and g.

On a Cartesian grid, only five values are needed to compute the Laplacian. In that case, we also deal with two extra ``hidden'' relationships that state the particular geometrical relationships between the five points. In the general case, we do not have that information and, therefore, need more points.

The Laplacian may be found in other ways as well; perhaps interpolating with splines or other basis functions the entire wavefield - not just local neighbourhoods - at each depth step. The fastest and most elegant approach would nevertheless not involve finding a complete analytical expression of the wavefield function, but only its curvature information represented by the Laplacian.