Our pseudoparticle method and Anderson's method are quite similar. Both approximate the multipole expansion by a function on a circle (two dimensions) or a sphere (three dimensions). The difference is that the value of potential is used in Anderson's method and the mass distribution itself is used in our method.
Both the potential and mass distribution are formally defined by the inverse transform from the multipole expansion. The only difference lies in the multiplication factor applied in the inverse transform.
One practical advantage of our method is that the calculation of the M2L part, which is known to be the dominant part of the computation, is significantly simpler for our method. Therefore, the overall calculation speed for the same expansion order would be faster for our method.
In the case of the non-adaptive tree, one can pre-compute the translation matrix for M2L part. In this case, the calculation cost of our method and that of Anderson's method is essentially the same. In this case, theoretically, both scheme would be slightly slower than the original FMM, since the number of points needed to express expansions up to order p is somewhat larger than the number of terms, .
If we view Anderson's method as one way of expressing the multipole expansion, it seems clear that there is a room for improvement in the original formulation by Anderson. In the first transform from physical particles to the values of potential on the outer ring, he used the (or ) potential without truncation. This treatment naturally introduces fictitious high-frequency terms in the potential on the ring. This is the reason why he had to carefully chose the radius of the ring to suppress the high-order terms. These high-order terms contaminate the values of potential itself through aliasing, unless we make the radius of the ring sufficiently large. On the other hand, if we make the radius of the ring too large, the solution inside the ring tends to be quite inaccurate.
If we use the truncated form of the potential for the first conversion, the values of potential on the ring always represent exact values of the multipole expansion, for any choice of the radius of the ring. Thus, such choice should improve the accuracy of Anderson's method significantly, without increasing the calculation cost.
In our method, untruncated potential is used in M2L part. This also introduce the high-order terms which are not in the original multipole expansion. However, this does not cause much degradation in accuracy, since here it is guaranteed that the contribution of the high-order terms is small.
Recently, several methods which reduce the cost of the translation of the multipole expansion from to have been proposed [16,8]. White and Head-Gordon [16] described a very clever transformation, in which they rotate the coordinates so that the axis of symmetry becomes parallel to the line which connects the two centers of expansion. The cost of the rotation is , and the translation along the symmetry axis is also . Greengard and Rokhlin [8] described a more complex method with theoretically better scaling.
Timing results reported [16,8] indicate that these new methods, like the method described in [5], are advantageous only for pretty large values of p. For example, White and Head-Gordon [16] reported the maximum speedup over FMM of a factor of 2.5 for p=21. Typical speedup for p=21 was around a factor of two.
In practice, for many problems expansions of order 2 to 8 give sufficient accuracy. For these problems, our method will be competitive with these new methods. Of course, a careful and detailed timing comparison under realistic circumstance will be necessary to draw a firm conclusion.
Our group has developed a series of special-purpose computers for N-body problem [14,12]. The basic function of these machines is to evaluate and accumulate the gravitational interaction between particles. Though a modified version of the tree algorithm has been implemented on these machines [10], previously only the monopole (effectively dipole) approximation could be used, since the hardware could only calculate the interaction between point particles.
In our pseudoparticle method, the high order expansion is expressed by means of particles, which means we can use the special-purpose hardware to evaluate high order expansion. Our pseudoparticle method combined with special-purpose hardware will provide a very large speed advantage over FMM or tree algorithms on general-purpose computer.
As we discussed in section 3.1, the present implementation uses rather large number of particles to represent the multipoles. This is because we put a stringent restriction to the placement of particles: We fix the positions and only vary the mass of particles. This restriction makes it possible to obtain the mass of particles using linear convolution. However, this is certainly not optimal. For example, it is clear that a monopole can be exactly expressed by one particle, while our method require two. In the case of gravitational potential, a dipole term can also be expressed by one particle placed at the center of mass, while our method requires four. A quadrupole can be expressed by four particles, while out method requires 12. For these low-order expansions, it would be relatively easy to obtain the location of particles without actually solving the nonlinear equation.
For orders higher than 2, the calculation cost of solving the nonlinear equation would be too high.