A grid-based boundary integral method for elliptic problems in three dimensions

We develop a simple, efficient numerical method of boundary integral type for solving an elliptic partial differential equation in a three-dimensional region using the classical formulation of potential theory. Accurate values can be found near the boundary using special corrections to a standard quadrature. We treat the Dirichlet problem for a harmonic function with a prescribed boundary value in a bounded three-dimensional region with a smooth boundary. The solution is a double layer potential, whose strength is found by solving an integral equation of the second kind. The boundary surface is represented by rectangular grids in overlapping coordinate systems, with the boundary value known at the grid points. A discrete form of the integral equation is solved using a regularized form of the kernel. It is proved that the discrete solution converges to the exact solution with accuracy O(h(p)), p < 5, depending on the smoothing parameter. Once the dipole strength is found, the harmonic function can be computed from the double layer potential. For points close to the boundary, the integral is nearly singular, and accurate computation is not routine. We calculate the integral by summing over the boundary grid points and then adding corrections for the smoothing and discretization errors using formulas derived here; they are similar to those in the two-dimensional case given by [J. T. Beale and M.-C. Lai, SIAM J. Numer. Anal., 38 ( 2001), pp. 1902-1925]. The resulting values of the solution are uniformly of O(h(p)) accuracy, p < 3. With a total of N points, the calculation could be done in essentially O(N) operations if a rapid summation method is used.