  Integral equation methods for the solution of partial differential equations,
when coupled with suitable fast algorithms, yield geometrically flexible,
asymptotically optimal and well-conditioned schemes in either interior or
exterior domains. The practical application of these methods, however, requires
the accurate evaluation of boundary integrals with singular, weakly singular or
nearly singular kernels. Historically, these issues have been handled either by
low-order product integration rules (computed semi-analytically), by
singularity subtraction/cancellation, by kernel regularization and asymptotic
analysis, or by the construction of special purpose "generalized Gaussian
quadrature" rules. In this paper, we present a systematic, high-order approach
that works for any singularity (including hypersingular kernels), based only on
the assumption that the field induced by the integral operator is locally
smooth when restricted to either the interior or the exterior. Discontinuities
in the field across the boundary are permitted. The scheme, denoted QBX
(quadrature by expansion), is easy to implement and compatible with fast
hierarchical algorithms such as the fast multipole method. We include accuracy
tests for a variety of integral operators in two dimensions on smooth and
corner domains.
