A kernel-independent adaptive fast multipole algorithm in two and three dimensionsIntroductionSynopsis of the new methodRelated workOrganization of the paperReview of the fast multipole methodThe new algorithmDensity translationsEquivalent densities and check potentialsM2M translationM2L translationL2L translationSummaryDiscretizationRegularizationSurface geometry and discretization2D case3D caseSummaryThe complete algorithmDefinitionsImplementation issuesAcceleration techniquesSVD-based acceleration (2D)FFT-based acceleration (3D)HomogeneitySymmetryLazy computationError analysisFMM factorizationDiscretization errorSingle step errorNumerical resultsAccuracy on the equivalent density approximation2D case3D caseOverall approximation errorConclusions and future workKernelsCoefficients of fast multipole methodProofs of lemmasReferencesA kernel-independent adaptive fast multipole algorithmin two and three dimensionsqLexing Ying, George Biros*, Denis ZorinCourant Institute of Mathematical Sciences, New York University, New York 10012, USAReceived 17 July 2003; received in revised form 6 November 2003; accepted 11 November 2003AbstractWe present a new fast multipole method for particle simulations. The main feature of our algorithm is that it doesnot require the implementation of multipole expansions of the underlying kernel, and it is based only on kernelevaluations. Instead of using analytic expansions to represent the potential generated by sources inside a box of thehierarchical FMM tree, we use a continuous distribution of an equivalent density on a surface enclosing the box. Tofind this equivalent density, we match its potential to the potential of the original sources at a surface, in the far field, bysolving local Dirichlet-type boundary value problems. The far-field evaluations are sparsified with singular value de-composition in 2D or fast Fourier transforms in 3D. We have tested the new method on the single and double layeroperators for the Laplacian, the modified Laplacian, the Stokes, the modified Stokes, the Navier, and the modifiedNavier operators in two and three dimensions. Our numerical results indicate that our method compares very well withthe best known implementations of the analytic FMM method for both the Laplacian and modified Laplacian kernels.Its advantage is the (relative) simplicity of the implementation and its immediate extension to more general kernels. 2003 Elsevier Inc. All rights reserved.Keywords: Fast multipole methods; Fast solvers; Integral equations; Single-layer potential; Double-layer potential; Particle methods;N-body problems1. IntroductionMany methods in computational physics (e.g., vortex methods, molecular dynamics) are based on theevolution of particle systems with pairwise interactions corresponding to potentials related to the funda-mental solution of elliptic partial differential equations (PDEs). The most important among these kernels isthe single-layer Laplacian. Other kernels include the kernels of the Stokes and Navier ope rators, theirmodified versions, and their derivatives (double-layer and hypersingular kernels).www.elsevier.com/locate/jcpJournal of Computational Physics 196 (2004) 591–626qThis work is supported by the National Science FoundationÕs Knowledge and Distributed Intelligence (KDI) program throughGrant DMS-9980069.*Corresponding author. Tel.: +1-212-998-3358; fax: +1-212-995-4122.E-mail addresses: [email protected] (L. Ying), [email protected], [email protected] (G. Biros), [email protected] (D. Zorin).0021-9991/$ - see front matter  2003 Elsevier Inc. All rights reserved.doi:10.1016/j.jcp.2003.11.021Particle formulations result in dense linear algebraic systems because all pairwise interactions have to becomputed. This is a significant bottleneck since for N particles and results in a OðN2Þ computation. In orderto make large scale problems tractable it is essential to efficiently compute these interactions. A number ofalgorithms have been proposed for this purpose. The fast multipole method (FMM) has been one of themost successful, especially for nonuniform particle distributions.In this paper, we present a new kernel-independent FMM-like algori thm. Our algorithm has thestructure of the adaptive FMM algorithm [12] but requ ires only the kernel evaluations, and it does notsacrifice the efficiency of the original algorithm. The crucial element of our approach is to replace theanalytic expansions and translations with equivalent density representations. These repres entations arecomputed by solving local exterior and interior problems on circles (2D), spheres or cub es (3D) using theintegral equation formulations. We demonstrate the efficiency of our method in both 2D and 3D for manykernels: the single and double layer potentials of the Laplacian, the modified Laplacian, the Navier, theStokes, and their modified variants. Our method has OðNÞ asymptotic complexity, and, like analytic FMM,works well for nonuniform particle distributions.1.1. Synopsis of the new methodThe basic structure of our method follows [14], the original fast multipole method, which we brieflyreview in Section 2. FMM consists of the following steps:1. generation of a hierarchical tree partitioning of the computational domain;2. accumulation of the multipole expansions for the far field by a postorder traversal of the tree;3. translation of the multipole moments to the local expansions;4. construction of local expansions by a pre order traversal of the tree;5. evaluation of the far-field action on the particles using local expansions;6. evaluation of the near field interactions.The same steps are used in our algorithm. However in the postorder traversal of the tree, the multipoleexpansion construction is replaced by solving local exterior inverse problems. To represent the potentialgenerated by pa rticles inside a box, we use a continuous distribution of an equivalent density on a surfaceenclosing the box. To find this equivalent den sity on the surface, we match its potential to the potential ofthe original sources at another surface in the far field. The translations are done by direct evaluation on thefar field, sparsified with SVD or FFT. During the preorder traversal of the tree, we evaluate the far-fieldinteraction on a surface enclosing a target box, and solve an interior Dirichlet-type integral equation tocompute an equivalent density. Then we use this density to represent the potential inside a target