  We present a method for computing nearly singular integrals that occur when
single or double layer surface integrals, for harmonic potentials or Stokes
flow, are evaluated at points nearby. Such values could be needed in solving an
integral equation when one surface is close to another or to obtain values at
grid points. We replace the singular kernel with a regularized version having a
length parameter $\delta$ in order to control discretization error. Analysis
near the singularity leads to an expression for the error due to regularization
which has terms with unknown coefficients multiplying known quantities. By
computing the integral with three choices of $\delta$ we can solve for an
extrapolated value that has regularization error reduced to $O(\delta^5)$,
uniformly for target points on or near the surface. In examples with $\delta/h$
constant and moderate resolution we observe total error about $O(h^5)$ close to
the surface. For convergence as $h \to 0$ we can choose $\delta$ proportional
to $h^q$ with $q < 1$ to ensure the discretization error is dominated by the
regularization error. With $q = 4/5$ we find errors about $O(h^4)$. For
harmonic potentials we extend the approach to a version with $O(\delta^7)$
regularization; it typically has smaller errors but the order of accuracy is
less predictable.
