To get a better intuition for how the Fast Multipole Method (FMM) approximates n-body interactions, I implemented the equations in Mathematica along with the exact solution. Below are some plots for different situations. The thing to keep in mind is that the FMM is valid only a certain distance from the center of the particle cluster, so it’s expected that the approximation will not be accurate in that region.

The first row shows contour plots of clusters of 3 and 40 points and a difference between them. With the difference plot, I was trying to get a sense of how the potential field varies under differing conditions.
The second row is again for 3 and 40 points with a difference. The cutoff spike in the difference surface is an artifact of numerical evaluation. Since the potential in 2D is logarithmic, near the particles it tends to blow up. The coloring of the surfaces follows the phase of the potential. With 3 particles, there potential wave has 3 clear periods while with 40 particles, it starts to get noisier (an artifact of sampling), but there are likewise 40 periods. The reason the number of particles matches the number of cycles is from the linear summing of the phase term.
The last row shows the electrical field as derived from the complex potential. The two images are for 3 particles with the one on the left being exact and on the right calculated from the FMM method. These plots only show the calculation of the multipole expansion. To actually implement the FMM as a solution of the n-body problem, someadditional techniques for manipulating the multipole expansion are required.
The FMM simplifies the n-body problem by relating clusters of particles to other clusters of particles. As the multipole expansion shows, the affect on the potential field a cluster of particles generates can be easily calculated. The trick now is to figure out how to incorporate the multipole expansion into a hierarchical data structure without having to recalculate it at every single level. Instead, what we’d like to do is calculate the multipole expansion at the lowest level and accumulate the expansions up the hierarchy through a simple summing of terms. The problem is that each low-level expansion is centered at a different point, so we need to be able to translate the expansions to the center of higher-level (parent) cells.

The translation lemma shows how to accumulate the multipole expansion into coarser clusters, which can dramatically speed up the calculation for faraway clusters. This is the bottom up step of the FMM technique. As has been emphasized a number of times, however, the multipole expansion is only valid for well-separated clusters that are a certain distance from each other. This means that only clusters that are faraway can be treated this way or, alternatively, we can reduce the size of the clusters in order to apply the technique to closer in clusters. In order to do this, each box in the hierarchy is given an interaction list.. The interaction list says which boxes interact with each other. The interaction list consists of all of the child nodes of the neighbors of a box’s parent node not bordering it.

At each level, the interaction list covers a smaller and more refined area of the space. Any area not covered by the interaction list, which is essentially those boxes at the finest level bordering the box under consideration, are calculated using the O(N^2) n-body method directly, but these calculations should be quite minimal.
In the diagram above, the pink box affects affects the blue boxes. In order to make use of the multipole expansion, the expansion within the blue boxes has to be recentered on the pink box since expansions can only be summed in a simple way if they are centered on the same location. This process is described as converting a multipole expansion into a local expansion where local is the pink box.

Once this calculation is done, we have a description for the effects of every particle outside the box, except for those nearest neighbors, and can evaluate their effects simply by plugging in the locations of the particles in the box into the sum of all of the multipole expansions. This sounds like a lot of work compared to the O(N^2) formulation, but for thousands or millions of particles, it’s vastly more scalable. The trickiest part is all of the bookkeeping that has to be done to shift the multipole expansions around the space and up and down the spatial hierarchy. The payoff though, once it’s encoded in a easy to use software package should be significant and we won’t have to think about the hairy details so much anymore and can get to figuring out what massively detailed particle simulations can be used for in terms of spatial composition.
Note: The lemmas posted above were taken from Greengard and Beatson, A Short Course in Fast Multipole Methods.
Download: FMM Mathematica Files