04- 快速多极展开

1. 相关概念

1.1 快速多极算法简介

快速多极展开算法（FMM）最初是用来处理大量粒子之间相互影响的问题，也可以用来计算基于 位势理论 的拉普拉斯方程的数值解，同时也可用来处理泊松方程。

同理，我们也可以用相关场的求解过程来辅助理解 FMM 求解泊松方程的这一过程，该处以多体系统的静电场为例。

1.2 静电场

1.2.1 公式的离散形式

对于空间 N 个带电粒子构成的系统，第 i 个粒子所在位置的静电势 $\Phi(x_i)$ ：

$\Phi(x_i) = \sum^N_{j=1, j\neq i} \cfrac{m_j}{r_{ij}}$

所受到的静电力 $E(x_i)$ ：

$E(x_i) = \sum^N_{j=1, j\neq i} m_j \cfrac{\bold{x}_i - \bold{x}_j}{r_{ij}^3} = \sum^N_{j=1, j\neq i} m_j \cfrac{\bold{x}_i - \bold{x}_j}{\Vert \bold{x}_i - \bold{x}_j \Vert^3}$

Note ：

• $\bold{x}_i$ ：表示第 i 个粒子所在位置；
• $m_j$ ：表示第 j 个粒子所占权重，与带电量 $q_j$ 成正比；
• $r_{ij}$ ：表示第 i 个粒子和第 j 个粒子之间的距离；
• $\Vert \cdot \Vert$ ：表示两点间的距离；
• 求解 $E(x_i)$ 相当于对 i 粒子与 i 粒子以外的所有粒子的相互力做一个累和，复杂度为 $O(N)$ ，同理，如果对全部 N 个粒子都使用该方式进行求解，复杂度高达 $O(N^2)$ 。

1.2.2 公式的积分形式

对于空间连续介质问题，假设粒子是连续分布在空间曲面 $S$ 上，则离散的多粒子系统可以过渡到连续介质系统。在点 $\bold{x}$ 处受到整个曲面的势为：

$\Phi(\bold{x}) = \int_S \cfrac{m(\bold{y})}{\vert \bold{x} - \bold{y} \vert} \mathrm{d}S(\bold{y}),\bold{y} \in S$

点 $\bold{x}$ 处受到整个曲面的力为：

$E(\bold{x}) = \int_S m(\bold{y}) \cfrac{\bold{x} - \bold{y}}{\vert \bold{x} - \bold{y} \vert^3} \mathrm{d}S(\bold{y}),\bold{y} \in S$

Note ：

• $\vert \cdot \vert$ ：表示两点间的距离。

1.2.3 公式的复势形式

课程中为了式子尽可能的直观，主要是以该表达方式为例。

对于单个粒子产生的场

位于复平面上的点 $\bold{z}_j$ 的平面点电荷 $q_j$ 在点 $\bold{z}$ 的复势：

$f_j(\bold{z}) = \cfrac{-q_j}{2\pi \epsilon_0} \ln(\bold{z} - \bold{z}_j)$

其中的实部为势函数：

$\phi_j(\bold{z}) = \mathrm{Re}f_j(\bold{z}) = \cfrac{q_j}{2\pi \epsilon_0} \mathrm{Re}(\ln(\bold{z} - \bold{z}_j))$

Note ：

• 虚部为通量函数，由于该处我们仅讨论势函数，故忽略。

对于 j 个粒子组成的系统产生的场

点 $\bold{c}$ 附近的 $m$ 个电荷对点 $\bold{z}$ 的势（直接累和）：

$\phi(\bold{z}) = \cfrac{1}{2\pi \epsilon_0} \mathrm{Re} \sum^m_{j = 1} q_j(\ln(\bold{z} - \bold{z}_j))$

该处简化为：

$\phi(\bold{z}) = \sum^m_{j = 1} q_j\ln(\bold{z} - \bold{z}_j) \tag{1}$

将 $\ln(\bold{z} - \bold{z}_j)$ 在点 $\bold{z}$ 处做泰勒级数展开，即：

$\ln(\bold{z} - \bold{z}_j) = \ln \bold{z} - \sum^{\infin}_{k=1} \cfrac{\bold{z}^k_j}{k \bold{z}^k} \tag{2}$

将公式（2）代入公式（1），得到：

$\phi(\bold{z}) = \sum^m_{j = 1} q_j \ln \bold{z} - \sum^m_{j = 1} \sum^{\infin}_{k=1} \cfrac{q_j \bold{z}^k_j}{k \bold{z}^k} \tag{3}$

不妨定义一个总电荷变量 $Q$ 、 $Q_k$ 如下：

$\begin{cases} Q = \sum^m_{j = 1} q_j \\[2ex] Q_k = - \sum^m_{j = 1} \cfrac{q_j \bold{z}^k_j}{k} \end{cases} \tag{4}$

将（4）式代入（3）式，得到：

$\phi(\bold{z}) = Q \ln \bold{z} - \sum^{\infin}_{k=1} \cfrac{Q_k}{\bold{z}^k} \tag{5}$

对于 $\ln(\bold{z} - \bold{z}_j)$ 的泰勒展开：

$\begin{aligned} \ln(\bold{z} - \bold{z}_j) &= \cfrac{\ln \bold{z}}{0!} + (-\bold{z}_j) \cfrac{\cfrac{1}{\bold{z}}}{1!} + (-\bold{z}_j)^2 \cfrac{- \cfrac{1}{\bold{z}^2}}{2!} + (-\bold{z}_j)^3 \cfrac{(-1)^2 2 \cfrac{1}{\bold{z}^3}}{3!} + \cdots \\[2ex] &+ (-\bold{z}_j)^k \cfrac{(-1)^{k-1} (k-1)! \cfrac{1}{\bold{z}^k}}{k!} + \cdots \\[2ex] &= \ln \bold{z} - \sum^{\infin}_{k=1} \cfrac{\bold{z}^k_j}{k \bold{z}^k} \end{aligned}$

1.3 相关方程

1.3.1 拉普拉斯方程（调和函数）

如果一个函数满足拉普拉斯方程：

$\nabla^2 \Phi = \cfrac{\partial^2 \Phi}{\partial x^2} + \cfrac{\partial^2 \Phi}{\partial y^2} + \cfrac{\partial^2 \Phi}{\partial z^2} = 0$

则称它为调和函数。

• 位势理论 ：满足拉普拉斯方程的函数的理论。

1.3.2 球面调和函数（球谐函数）

球坐标下的三维拉普拉斯方程：

$\cfrac{1}{r^2} \cfrac{\partial}{\partial r}(r^2 \cfrac{\partial \Phi}{\partial r}) + \cfrac{1}{r^2 \sin{\theta}} \cfrac{\partial}{\partial \theta}(\sin{\theta} \cfrac{\partial \Phi}{\partial \theta}) + \cfrac{1}{r^2 \sin^2{\theta}} \cfrac{\partial^2 \theta}{\partial \phi^2} = 0$

用分离变量法解出的该方程的标准解为：

$\Phi = \sum^{\infty}_{n=0} \sum^n_{m=-n}(L^m_n r^n + \cfrac{M^m_n}{r^{n+1}}) Y^m_n(\theta, \phi)$

Note ：

• n 次球面调和函数 ： $Y^m_n(\theta, \phi)r^n$ ；
• -(n+1) 次球面调和函数 ：$Y^m_n(\theta, \phi)/r^{n+1}$ ；
   
• 系数 $L^m_n$ 、 $M^m_n$ 称为 展开的矩 ，它们仅与电荷的本质属性有关（电荷的位置和带电量）；
  • 在多极展开中，系数 $L^m_n = 0$ ；
  • 在局部展开中，系数 $M^m_n = 0$ ；

 

• 勒让德多项式 ：

$Y^m_n(\theta, \phi) = \sqrt{\cfrac{(n - \vert m \vert)!}{(n + \vert m \vert)!}} \cdot P^{\vert m \vert}_n(\cos{\theta})e^{im\phi}$

• 其中 $P^{\vert m \vert}_n(\cos{\theta})$ 的值可通过以下递推式求得：

$\begin{cases} (2n+1) x P^m_n(x) = (n-m+1)P^m_{n+1}(x) + (n+m)P^m_{n-1} \\[2.5ex] P^{m+2}_n(x) + 2(m+1) \cfrac{x}{\sqrt{x^2 - 1}} P^{m+1}_n(x) = (n-m)(n+m+1)P^m_n(x) \end{cases}$

2. 快速多级算法

2.1 关键思路

• 有一个指定的可接受的计算精度 $\epsilon$ ；
• 将空间分层细分为 源（sources） 的 平面（panels） 或 簇（clusters） ；
• 一个 核函数（kernel） $k(\bold{x}, \bold{y})$ 的远场展开，其中源（source）和评价点（evaluation points）的影响是分开的；
• 远场展开到局部展开的转换。

2.2 网格分层构建

• 选取一个能够将所有粒子都包含进去的正方形平面作为计算的定义域，对其进行分成；
• 规定该定义域为第 0 层网格（即最大的单个方格）；
• 第 L+1 层网格是第 L 层网格进行 2×2 的划分，同时，从第 L 层中的某个网格 box 中划分出来的 4 个网格是网格 box 的子结点；
• 显然，对于第 L 层，所有网格的个数为 $4^L$ （如果是三维情况，八叉树情况下网格个数为 $8^L$ ）。

2.3 网格之间的关系

• 以下图片来自 《快速多极边界元法在动力学问题中的若干应用_姚振汉》

2.3.1 邻接网格（nearest box）

邻接网格（nearest box） ：对于第 $L$ 层的某个网格 $S_1$ ，其 “邻接网格” 是同第 $L$ 层且有共同边界点的网格，图中的 $S_2$ 即为 $S_1$ 的 “邻接网格” 。

2.3.2 分离网格对（well separated pairs）

分离网格对（well separated pairs） ：两个网格在同一层，且不是 “邻接网格” ，则称为两个网格是 分离关系（well separated） 。

2.3.3 交互列表（interaction list）

交互列表（interaction list） ：对于第 $L$ 层的某个网格 $S_1$ ， $S_1$ 的父结点的 “邻接网格” 的子结点的集合称之为 $S_1$ 的 “交互列表” ，同时集合中的网格与 $S_1$ 是 “分离关系” 。

• “交互列表” 是一组结点的集合，这组结点不是 $S_1$ 的 “邻接网格” ，是至少间隔一个网格的；
• 在某层，对于其中的每一个网格，均存在一组对应的 “交互列表” ；
• “交互列表” 中的网格，是当前层待多极展开计算的网格。

3. 算法流程概述

• 图片来自 《基于快速多极展开算法的直升机旋翼尾迹计算分析_丁国华》

3.1 Tree Code（NlogN）

3.1.1 递归起始条件

• 首先，从第 0 层到第 1 层，是一个 2×2 的网格，不存在 “分离网格对”；
• 在第 2 层，存在若干个 “分离网格对” ，对其使用多极展开计算交互；
• 计算分两个部分：
  • 计算网格 $\bold{X}$ 内部粒子之间的交互；
  • 计算白色网格（与网格 $\bold{X}$ 是 “分离关系” 的网格）的多极展开；
  • 网格 $\bold{X}$ 以外的灰色网格（网格 $\bold{X}$ 的“邻接网格”）不计算。
• 在计算过程中，需要考虑误差的界限，对于精度 $\epsilon$ ，展开项数 $p$ 为 $p = \log_2(1/ \epsilon)$ 。

3.1.2 递归过程

• 显然，在第 2 层中，每个网格的粒子簇和其邻接网格中粒子簇的交互还需要去计算。

接下来讨论第 3 层：

• 在第 3 层中，父网格的邻接网格以外的网格对网格 $\bold{X}$ 的交互已经计算完毕（图中亮灰色网格，对应的第 2 层的白色网格）；
• 在当前层中，网格 $\bold{X}$ 与其邻接网格的交互无法计算（图中深灰色网格）；
• 排除以上两个部分，剩下的网格即是上述定义的 “交互列表” （图中白色网格），此时是 “分离关系”，可以使用多极展开计算交互。

Note ：

• 第 3 层中的网格 $\bold{X}$ 可以看成是第 2 层中的网格 $\bold{X}$ 进一步的细分，目标粒子同时存在于这两个网格中；
• 第 3 层中的亮灰色网格（“已计算网格”）是在第 2 层中已经计算完毕的白色网格（“交互列表”，或者形象的称之为“待计算网格”）；
• 第 2 层中的灰色网格（“邻接网格”），有一部分在第 3 层中变成了白色（“交互列表”），此时，变成白色的部分可以在当前层进行多极展开计算；
• 通过这种方式进行递归，将未计算交互的部分（即上一层的 “邻接网格” 部分）一层层的定义成白色（“交互列表”）进行计算，对于目标粒子来说，来自 “远场” 的交互，像是一个从外向内，越来越薄的洋葱皮。

3.1.3 算法分析

整体角度 ：

• 递归树的高度为 $\log N$ ；
• 每一层的工作量约为 $O(N)$ ，对于每个多极展开的项数为 $p$ ，创建所有展开式的操作数为 $Np$ ；

单个粒子的角度 ：

• 最多有 27 个网格需要进行展开计算（27 是 “交互列表” 的最大网格数），因此操作数为 $27Np$ ：
  • 图中红色网格代表的无法计算的 “邻接网格” ；
  • 在第 $L+1$ 层中，白色网格部分是当前层需要进行计算交互的 “交互列表”，同时也是第 $L$ 层中的 “邻接网格” 的一部分；
  • 可见当网格 $\bold{X}$ 处于中间方位的时候，白色网格部分最大值为 27 。

递归出口 ：

• 在最下面一层，大概划分出了 $4^{\log_4 N} = N$ 个网格；
• 还剩下目标粒子所在网格的 “邻接网格” 未计算交互；
• 根据均匀性假设，每个网格的粒子数为 $O(1)$ ，最后一步是目标粒子所在的网格与周围 8 个 “邻接网格” 进行直接计算，操作数为 $8N$ 。

总的时间花费为：

$28 N p \log N + 8 N$

3.1.4 自适应算法

当粒子的分布是不均匀的情况下，我们可以使用一种自适应的数据结构，确保每个网格都含有粒子，且粒子数目适中。

在网格划分的过程中，检查每个网格是否包含有粒子：

• 若存在粒子，则继续划分该网格；
• 若不存在粒子，将该网格及其兄弟网格从树中剪枝，回溯到父网格，提前结束递归，不参与后续的划分操作。

3.2 快速多极展开

向上遍历 ：

• 多极展开计算 P2M ：计算第 $n$ 层中每个网格中的所有粒子在网格中心的电势的多极展开。

• 多极展开的平移 M2M（子结点到父结点） ：从第 $n-1$ 层到第 $0$ 层，计算每个网格中心点势能在其父网格中心点的电势的多极展开。

向下遍历 ：

• 多极展开向局部展开的转换 M2L ：

• 局部展开的平移 L2L （父结点到子结点）：

• 局部展开的计算 L2P ：
• 直接计算 P2P ：

3. 多极展开计算（P2M）

计算一个远场区域内的 $m$ 个电荷对点 $\bold{z}$ 的势：

$\phi(\bold{z}) = \sum^m_{j = 1} q_j \ln \bold{z} - \sum^m_{j = 1} \sum^{\infin}_{k=1} \cfrac{q_j \bold{z}^k_j}{k \bold{z}^k} \tag{3}$

令远场中心点的电荷为远场的总电荷 $Q$ ：

$\begin{cases} Q = \sum^m_{j = 1} q_j \\[2ex] Q_k = - \sum^m_{j = 1} \cfrac{q_j \bold{z}^k_j}{k} = - Q \cfrac{\sum^m_{j = 1} \bold{z}^k_j}{k} \end{cases} \tag{4}$

用式（4）替换式（3）中的对应系数：

$\phi(\bold{z}) = Q \ln \bold{z} + \sum^{\infin}_{k=1} \cfrac{Q_k}{\bold{z}^k} \tag{5}$

当然，在实际实现多极展开时，需要确定多极展开的项数 $p$ ， $p$ 的选取决定了计算结果的精度，这里直接表示在式中泰勒展开的级数：

$\phi(\bold{z}) = Q \ln \bold{z} + \sum^{p}_{k=1} \cfrac{Q_k}{\bold{z}^k} \tag{6}$

4. 多极展开的平移(M2M)

设在中心 $\bold{c}_1$ 的半径圆内，存在一组强度为 $q_1, q_2, \cdots , q_j$ 的 $j$ 个电荷引起的势的多极展开，即 $M_1$ 网格，其在点 $\bold{z}$ 产生的势为：

$\phi_{\bold{c}_1}(\bold{z}) = Q \ln (\bold{z} - \bold{c}_1) + \sum^{p}_{k=1} \cfrac{a_k}{(\bold{z} - \bold{c}_1)^k} \tag{7}$

其中的高阶项 $a_k$ ：

$a_k = - Q \cfrac{\sum^m_{j = 1} (\bold{z}_j - \bold{c}_1)^k}{k}$

这里我们需要寻找到一个 $M_1 \mapsto M_2$ 的变换，在网格 $M_2$ 处找到一个展开式 $\phi_{C_1}(\bold{z})$ 使得 $\phi_{C_1}(\bold{z}) = \phi_{C_2}(\bold{z})$ ，相当于用 $\phi_{C_1}(\bold{z})$ 在点 $\bold{c}_2$ 处进行泰勒展开：

$\phi_{\bold{c}_2}(\bold{z}) = Q \ln (\bold{z} - \bold{c}_2) + \sum^{p}_{k=1} \cfrac{b_k}{\bold{z}^k} \tag{8}$

其关键在于高阶项 $b_k$ ， $b_k$ 是对 $Q_k$ 的一个泛化：

$b_k = - \cfrac{Q(\bold{c}_1 - \bold{c}_2)^k}{k} + \sum^k_{i=1} Q_i(\bold{c}_1 - \bold{c}_2)^{k-i} \begin{pmatrix} k-1 \\ i-1 \end{pmatrix} \tag{9}$

• 对一个分布若干粒子的平面进行网格划分；
• 对每个格子的中心点进行计算，主要计算 $Q$ ， $Q_k$ ， $\phi_{\bold{c}}(\bold{z})$ ；
• 对于离格子较远处的某个点 $\bold{z}$ （换句话说，该格子对于这个点 $\bold{z}$ 来说是个 远场 ），该点受到的格子产生的势的影响 $\phi_{\bold{c}}(\bold{z})$ ，可被 $Q$ ， $Q_k$ 近似表示，见式（6）；
• 对于每一层，均执行一次向上抽象，即计算 $Q$ ， $Q_k$ ， $\phi_{\bold{c}}(\bold{z})$ 。

4.2 多极展开公式（multipole ）

5. 多极展开向局部展开变换(M2L)

6. 局部展开的平移(L2L)

7. 局部展开计算

球坐标系

球坐标系向直角坐标系的转换

球坐标 $(r, \theta, \varphi) \mapsto$ 直角坐标 $(x, y, z)$

$\begin{cases} x = r \sin{\theta} \cos{\varphi} \\ y = r \sin{\theta} \sin{\varphi} \\ z = r \cos{\theta} \end{cases}$

球坐标系下两直线的夹角

若已知球坐标系中两个方向分别为 $\vec{r_1}(1,\theta_1, \phi_1)$ 和 $\vec{r_2}(1,\theta_2, \phi_2)$ ，求它们之间的夹角 $\alpha$ 。

• 先计算两个单位向量的直角坐标：

$\begin{cases} \vec{r_1} = (\sin{\theta_1} \cos{\phi_1}, \sin{\theta_1} \sin{\phi_1}, \cos{\theta_1}) \\[2ex] \vec{r_2} = (\sin{\theta_2} \cos{\phi_2}, \sin{\theta_2} \sin{\phi_2}, \cos{\theta_2}) \end{cases}$

• 再计算两个单位向量之间的夹角 $\alpha$ ：

$\begin{aligned} \cos{\alpha} &= \vec{r_1} \cdot \vec{r_2} \\ &= \sin{\theta_1} \cos{\phi_1} \sin{\theta_2} \cos{\phi_2} + \sin{\theta_1} \sin{\phi_1} \sin{\theta_2} \sin{\phi_2} + \cos{\theta_1} \cos{\theta_2} \\ &= \sin{\theta_1} \sin{\theta_2} (\cos{\phi_1} \cos{\phi_2} + \sin{\phi_1} \sin{\phi_2}) + \cos{\theta_1} \cos{\theta_2} \\ &= \sin{\theta_1} \sin{\theta_2} \cos{(\phi_2 - \phi_1)} + \cos{\theta_1} \cos{\theta_2} \end{aligned}$

参考资料 ：

• 二维电荷分布所产生电势的多极展开的简易方法
• 基于快速多极展开算法的直升机旋翼尾迹计算分析_丁国华
• 数学物理方法（梁昆森编）
• 球坐标系两方向的夹角