I'm trying to implement the Fast Multipole Method following the short course http://math.nyu.edu/faculty/greengar/shortcourse_fmm.pdf. Assume now that the task at hand is to compute the Coulomb interactions in a large particle system. If a point charge of strength $q$ is located at $P_0=\left(x_0, y_0, z_0\right)$, then the potential and electrostatic field due to this charge at a distinct point $P=(x, y, z)$ are given by $$\Phi=\frac{1}{\tilde{r}}$$ where $\tilde{r}$ denotes the distance between points $P_0$ and $P$. Theorem 5.2 (Multipole expansion) Suppose that $k$ charges of strengths $\left{q_i, i=1, \ldots, k\right}$ are located at the points $\left{Q_i=\left(\rho_i, \alpha_i, \beta_i\right), i=\right.$ $1, \ldots, k}$, with $\left|\rho_i\right|<a$. Then for any $P=(r, \theta, \phi) \in \mathbf{R}^3$ with $r>a$, the potential $\phi(P)$ is given by Equation (5-15): $$\phi(P)=\sum_{n=0}^{\infty} \sum_{m=-n}^n \frac{M_n^m}{r^{n+1}} \cdot Y_n^m(\theta, \phi)$$ where $$M_n^m=\sum_{i=1}^k q_i \cdot \rho_i^n \cdot Y_n^{-m}\left(\alpha_i, \beta_i\right)$$ Theorem 5.3 (Translation of a multipole expansion) Suppose that $l$ charges of strengths $q_1, q_2, \cdots, q_l$ are located inside the sphere $D$ of radius a with center at $Q=(\rho, \alpha, \beta)$, and that for points $P=(r, \theta, \phi)$ outside $D$, the potential due to these charges is given by the multipole expansion $$\Phi(P)=\sum_{n=0}^{\infty} \sum_{m=-n}^n \frac{O_n^m}{r^{\prime n+1}} \cdot Y_n^m\left(\theta^{\prime}, \phi^{\prime}\right)$$ where $P-Q=\left(r^{\prime}, \theta^{\prime}, \phi^{\prime}\right)$. Then for any point $P=(r, \theta, \phi)$ outside the sphere $D_1$ of radius $(a+\rho)$, $$\Phi(P)=\sum_{j=0}^{\infty} \sum_{k=-j}^j \frac{M_j^k}{r^{j+1}} \cdot Y_j^k(\theta, \phi)$$ and $$M_j^k=\sum_{n=0}^j \sum_{m=-n}^n \frac{O_{j-n}^{k-m} \cdot i^{|k|-|m|-|k-m|} \cdot A_n^m \cdot A_{j-n}^{k-m} \cdot \rho^n \cdot Y_n^{-m}(\alpha, \beta)}{A_j^k}$$ with $A_n^m$ defined by $$A_n^m=\frac{(-1)^n}{\sqrt{(n-m)!(n+m)!}}$$ and $$O_n^m=\sum_{i=1}^N q_i \cdot r_i^n \cdot Y_n^{-m}\left(\theta_i, \phi_i\right)$$ $O_n^m$ = Old Multipole (multipole expansion already computed). Its center is on the child center. $M_n^m$ = New Multipole. Multipole computed using the parent center but using the already computed $O_n^m$ considering the child center. I'm using the Sympy package to compute spherical harmonics http://balbuceosastropy.blogspot.com/2015/06/spherical-harmonics-in-python.html Spherical harmonics definition from sympy $$Y_n^m(\theta, \varphi):=\sqrt{\frac{(2 n+1)(n-m)!}{4 \pi(n+m)!}} \exp (i m \varphi) \mathrm{P}_n^m(\cos (\theta))$$ Spherical harmonics definition from short course $$Y_n^m(\theta, \phi)=\sqrt{\frac{(n-|m|)!}{(n+|m|)!}} \cdot P_n^{|m|}(\cos \theta) e^{i m \phi}$$ Note that in this case there is no $\frac{(2 n+1)}{4 \pi}$ in the root. Because this difference, this factor should appear in the implementaion of equation (5-15). It is in the function eval_potential in the line: term = (4 * np.pi / (2 * n + 1)) * (1 / (r_target**(n + 1))) * harmonic_val * multipole_moments[n, m + n] My implementation is OK for the Multipole expansion given by Equation (5-15). However, the Multipole translation, given by Equation (5-20) is not working properly. I checked all the code but I can't find the error. Can someone help me? Here is my code: from scipy.special import sph_harm, factorial import numpy as np def A_nm(n, m): """Compute A_n^m coefficient as defined in Greengard's paper.""" return (-1)**n / np.sqrt(factorial(n - abs(m)) * factorial(n + abs(m))) def M2M_theorem53(child_multipole, child_center, parent_center, p): """ Perform M2M translation based on Theorem 5.3 (Greengard, 1997). Translates multipole expansion from child center to parent center (origin). Parameters: child_multipole: np.array of shape (p+1, 2p+1), with child moments O_n^m at [n, m+n] child_center: Point object parent_center: Point object p: order of expansion Returns: parent_multipole: np.array of shape (p+1, 2p+1), with translated M_j^k """ parent_multipole = np.zeros((p + 1, 2 * p + 1), dtype=complex) # Vector from parent to child center dx = child_center.x - parent_center.x dy = child_center.y - parent_center.y dz = child_center.z - parent_center.z rho = np.sqrt(dx**2 + dy**2 + dz**2) if rho < 1e-14: return np.copy(child_multipole) # No translation needed alpha = np.arctan2(dy, dx) beta = np.arccos(dz / rho) for j in range(p + 1): for k in range(-j, j + 1): acc = 0.0j A_jk = A_nm(j, k) for n in range(j + 1): for m in range(-n, n + 1): j_minus_n = j - n k_minus_m = k - m if abs(k_minus_m) > j_minus_n: continue # Spherical harmonic is 0 outside its domain A_nm_val = A_nm(n, m) A_jn_km = A_nm(j_minus_n, k_minus_m) # Complex exponential factor: i^{|k| - |m| - |k - m|} phase = 1j ** (abs(k) - abs(m) - abs(k - m)) # Y_n^{-m}(alpha, beta) Y = sph_harm(-m, n, alpha, beta) term = (child_multipole[j_minus_n, k_minus_m + j_minus_n] * phase * A_nm_val * A_jn_km * (rho ** n) * Y / A_jk) acc += term parent_multipole[j, k + j] = acc return parent_multipole import numpy as np from scipy.special import sph_harm # from matplotlib import pyplot, rcParams # Not strictly needed for the core logic # %matplotlib inline # from mpl_toolkits.mplot3d import Axes3D # Not strictly needed for the core logic # --- Re-using your Point and Particle classes --- class Point(): """The class for a point.""" def __init__(self, coords=[], domain=1.0): if coords: assert len(coords) == 3, "the size of coords should be 3." self.x = coords[0] self.y = coords[1] self.z = coords[2] else: self.x = domain * np.random.random() self.y = domain * np.random.random() self.z = domain * np.random.random() def distance(self, other): return np.sqrt((self.x-other.x)**2 + (self.y-other.y)**2 + (self.z-other.z)**2) class Particle(Point): """The derived class for a particle.""" def __init__(self, coords=[], domain=1.0, m=1.0): Point.__init__(self, coords, domain) self.m = m # Mass of the particle (charge for potential) self.phi = 0. # Gravitational potential, initialized for direct sum # --- Your P2M (Particle-to-Multipole) function (with a small enhancement for mass and r=0 handling) --- def P2M(sources, center, p): """ Given sources and cell's information, return multipole moments. This computes M_n^m = sum q_i r_i^n * conj(Y_n^m(theta_i, phi_i)) """ xc, yc, zc = center.x, center.y, center.z # Coordinates of sources relative to the center x_rel = np.array([source.x for source in sources]) - xc y_rel = np.array([source.y for source in sources]) - yc z_rel = np.array([source.z for source in sources]) - zc masses = np.array([source.m for source in sources]) r_rel = np.sqrt(x_rel**2 + y_rel**2 + z_rel**2) phi_rel = np.arctan2(y_rel, x_rel) # Handle sources exactly at the center (r_rel = 0) # np.arccos(0/0) results in NaN, set theta to 0 for points at center theta_rel = np.arccos(z_rel / np.where(r_rel == 0, 1e-10, r_rel)) # Avoid division by zero theta_rel[r_rel == 0] = 0.0 # Standard convention for (0,0,0) in spherical coords multipole_moments = np.zeros((p + 1, 2 * p + 1), complex) for i in range(len(sources)): # Loop over each source particle for n in range(p + 1): for m in range(-n, n + 1): # q_i * r_i^n * conj(Y_n^m(theta_i, phi_i)) multipole_moments[n, m + n] += masses * (r_rel**n) * np.conj(sph_harm(m, n, phi_rel, theta_rel)) return multipole_moments # --- Your eval_potential (Multipole-to-Particle, M2P) function (corrected for multiple targets) --- def eval_potential(targets, multipole_moments, center, p): """ Given targets list, multipole moments and expansion center, return the array of target's potentials. This evaluates Phi(r') = sum (4pi/(2n+1)) * M_n^m / (r')^(n+1) * Y_n^m(theta', phi') """ potential_at_targets = np.zeros(len(targets), dtype=complex) for j, target in enumerate(targets): # Coordinates of target relative to the expansion center dx = target.x - center.x dy = target.y - center.y dz = target.z - center.z r_target = np.sqrt(dx**2 + dy**2 + dz**2) # This expansion is valid only if r_target > r_max_sources_in_cell. # For a target exactly at the center (r_target=0), the expansion is singular. # In a full FMM, P2P (direct) interaction handles targets close to sources/expansion center. if r_target < 1e-10: # If target is too close to center, expansion is invalid potential_at_targets[j] = np.nan # Indicate invalid result for this case continue phi_target = np.arctan2(dy, dx) theta_target = np.arccos(dz / r_target) current_potential_sum = 0.0j # Initialize sum for this specific target for n in range(p + 1): for m in range(-n, n + 1): # Y_n^m(theta', phi') harmonic_val = sph_harm(m, n, phi_target, theta_target) # Term for current (n, m) term = (4 * np.pi / (2 * n + 1)) * (1 / (r_target**(n + 1))) * harmonic_val * multipole_moments[n, m + n] current_potential_sum += term potential_at_targets[j] = current_potential_sum return potential_at_targets.real # Potentials are real # --- Your direct_sum for verification --- def direct_sum(sources, targets): """ Calculate the gravitational potential (target.phi) at each target particle using direct summation method. """ phi_values = np.zeros(len(targets)) for j, target in enumerate(targets): current_phi = 0.0 for source in sources: r = target.distance(source) if r < 1e-10: # Handle near-singularities / self-interaction current_phi += np.inf # For gravitational potential, it goes to infinity else: current_phi += source.m / r phi_values[j] = current_phi return phi_values # --- Test Setup for M2M --- if __name__ == "__main__": n_sources = 10 # number of particles in a child cell p = 8 # order of expansion # 1. Define child cell sources and center child_center_coords = [0.25, 0.25, 0.25] child_center = Point(child_center_coords) np.random.seed(42) # for reproducibility # Sources distributed within a small radius around the child center # Ensure they are within a sphere that defines the child cell radius_child = 0.1 child_source_coords = np.random.uniform(-radius_child, radius_child,(n_sources, 3)) + np.array(child_center_coords) child_sources = [Particle(coord.tolist(), m=1.0) for coord in child_source_coords] # 2. Compute child's multipole moments (P2M) child_multipole = P2M(child_sources, child_center, p) print(f"Child multipole moments computed for {n_sources} sources.") # 3. Define parent cell center # Parent center should be "outside" the child cell, making d non-zero parent_center_coords = [0.5, 0.5, 0.5] # E.g., a larger cell's center parent_center = Point(parent_center_coords) print(f"Parent center: {parent_center_coords}") print(f"Child center: {child_center_coords}") # 4. Perform M2M translation parent_multipole_direct = P2M(child_sources, parent_center, p) parent_multipole_M2M = M2M_theorem53(child_multipole, child_center, parent_center, p) print(f"M2M translation performed.") # 5. Define a distant target to test potentials # The target must be sufficiently far from *both* child and parent centers # for the multipole expansions to be accurate. target_coords = [[2.0, 2.0, 2.0]] # Very distant target targets = [Particle(coord) for coord in target_coords] print(f"Target for evaluation: {target_coords[0]}") # 6. Evaluate potential using child_multipole (centered at child_center) # The target must be outside the smallest sphere enclosing the sources centered at child_center phi_from_child_m = eval_potential(targets, child_multipole, child_center, p) print(f"Potential from child multipole: {phi_from_child_m[0]:.12e}") # 7. Evaluate potential using parent_multipole (centered at parent_center) # The target must be outside the smallest sphere enclosing the sources centered at parent_center phi_from_parent_m = eval_potential(targets, parent_multipole_M2M, parent_center, p) print(f"Potential from parent multipole: {phi_from_parent_m[0]:.12e}") # 8. Direct sum for comparison (ground truth) # This involves directly summing contributions from all original child sources to the target phi_direct = direct_sum(child_sources, targets) print(f"Potential from direct sum: {phi_direct[0]:.12e}") # 9. Check agreement error_child_vs_direct = np.abs(phi_from_child_m[0] - phi_direct[0]) error_parent_vs_direct = np.abs(phi_from_parent_m[0] - phi_direct[0]) error_m2m_consistency = np.abs(phi_from_parent_m[0] - phi_from_child_m[0]) print(f"\nErrors (absolute difference):") print(f" Child multipole vs Direct: {error_child_vs_direct:.2e}") print(f" Parent multipole vs Direct: {error_parent_vs_direct:.2e}") print(f" M2M consistency (Parent vs Child potential): {error_m2m_consistency:.2e}") if error_m2m_consistency < 1e-9: # A small tolerance for floating point comparisons print("\nM2M translation appears to be working correctly: " "Potential from parent's moments matches child's moments.") else: print("\nWarning: M2M translation might have an issue. " "Check the consistency error.") The result I obtained is: Child multipole moments computed for 10 sources. Parent center: [0.5, 0.5, 0.5] Child center: [0.25, 0.25, 0.25] M2M translation performed. Target for evaluation: [2.0, 2.0, 2.0] Potential from child multipole: 3.275435467795e+00 Potential from parent multipole: 1.123673185858e+00 Potential from direct sum: 3.275435467795e+00 Errors (absolute difference): Child multipole vs Direct: 8.26e-14 Parent multipole vs Direct: 2.15e+00 M2M consistency (Parent vs Child potential): 2.15e+00 Warning: M2M translation might have an issue. Check the consistency error. Continue reading...
