球谐函数原理

参考：球谐函数极其应用

球谐函数是什么

球谐函数是傅里叶级数的高维类比，由一组表示球体表面的基函数构成。
同时，球谐函数是Laplace算子角动量在三个维度上的特征函数。

球坐标下的Laplace方程

Laplace方程(Laplace’s equation)

拉普拉斯算符(Laplacian) ${▽}^2$被定义为梯度的散度( $▽\cdot ▽$或$△$ )，有点类似于单变量函数的2阶导数。

对于函数$u(x_1,x_2,...)$， Laplace方程为：
$$\sum _{i=1}^n\frac{{\partial }^{2}u}{\partial x_i^2}=0$$

它表示一个函数在它各个变量的二阶偏导下都为0。
在一维情况下，在如下方程中可以理解为，当函数为Laplace方程时，加速度为0.
$$x(t)=vt+\frac{1}{2}a{t}^2$$

其他适用的场景包括：

• 室内均衡的气温的Laplacian等于0。
• 平衡水膜下受力的Laplacian等于0。
• 静电场中的电压场的Laplacian等于0。

综上，我们得出Laplacian等于0的函数可看作一种平衡状态，因此也被称为谐波函数（Harmonic）

三维直角坐标系下的Laplace方程

$${▽}^{2}u=\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}+\frac{{\partial }^{2}u}{\partial z^2}=0(1)$$

我们要做的就是求出三维Laplace方程的通解。

将其转化为球坐标系

将如下球面坐标表示代入$(1)$式
{ x = r sin ⁡ θ cos ⁡ φ y = r sin ⁡ θ sin ⁡ φ z = r cos ⁡ θ ( 2 ) \left\{

\right. \quad (2) ⎩ ⎨ ⎧​x=rsinθcosφy=rsinθsinφz=rcosθ​(2)
得到
$$\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2\frac{\partial f}{\partial r}\right)+\frac{1}{r^2\sin \theta }\frac{\partial }{\partial \theta }\left(\sin \theta \frac{\partial f}{\partial \theta }\right)+\frac{1}{r^2{\sin }^2\theta }\frac{{\partial }^{2}f}{\partial {\phi }^2}=0(3)$$
推导原理可参照：Laplacian in Spherical Coordinates

解 Laplacian 方程

将球坐标函数分离

$$u(r,\theta ,\psi )=R(r)Y(\theta ,\phi )$$

• $R(r)$表示距离部分
• $Y(\theta ,\phi )$就是球谐函数

分离后得到：

$$\frac{1}{R}\frac{d}{dr}\left(r^2\frac{dR}{dr}\right)=-\frac{1}{\sin (\theta )Y}\frac{\partial }{\partial \theta }\left(\sin (\theta )\frac{\partial Y}{\partial \theta }\right)-\frac{1}{Y}\frac{1}{{\sin }^2(\theta )}\frac{{\partial }^{2}Y}{\partial {\phi }^2}(4)$$
方程的左边是关于R的函数，右边是关于$\phi ，\theta$的函数，两边等式成立，需要两边都等于一个常数。

设这个常数为$$l(l+1)$$得到两个新方程
d d r ( r 2 ∂ Y ∂ θ ) − l ( l + 1 ) R = 0 1 sin ⁡ ( θ ) ∂ ∂ θ ( sin ⁡ ( θ ) ∂ Y ∂ θ ) + 1 sin ⁡ 2 ( θ ) ∂ 2 Y ∂ φ 2 + l ( l + 1 ) Y = 0

drd​(r2∂θ∂Y​)−l(l+1)R=0sin(θ)1​∂θ∂​(sin(θ)∂θ∂Y​)+sin2(θ)1​∂φ2∂2Y​+l(l+1)Y=0​

• 第一个方程是距离方程
• 第二个方程是角度方程

我们将角度方程继续分离
$$Y(\theta ,\psi )=\Theta (\theta )\Phi (\phi )$$
得到
$$\frac{1}{\sin \theta }\frac{\partial }{\partial \theta }\left(\sin \theta \frac{\partial [\Theta (\theta )\Phi (\phi )]}{\partial \theta }\right)+\frac{1}{{\sin }^2\theta }\frac{{\partial }^2[\Theta (\theta )\Phi (\phi )]}{\partial {\phi }^2}+l(l+1)[\Theta (\theta )\Phi (\phi )]=0$$

将上式最终化简为

设：等式两边都等于参数$\lambda$

化简后为：

求解$\Phi$

求解$\Theta$ ，设 $cos\theta =x$，有

$$\left(1-x^2\right)\frac{d^2\Theta }{d{x}^2}-2x\frac{d\Theta }{dx}+\left[l(l+1)-\frac{m^2}{1-x^2}\right]\Theta =0$$

这个方程是一个特殊的方程，叫做：l阶连带勒让德方程

参考链接：球谐函数的理解

勒让德多项式

勒让德方程：
$$\frac{\mathrm{d}}{\mathrm{d}x}\left[\left(1-x^2\right)\frac{\mathrm{d}}{\mathrm{d}x}{P}_l(x)\right]+l(l+1)P_l(x)=0$$

我们仅在区间 $x\in [-1,1]$ 上考虑 $l$ 为非负整数的情况，方程的解 $P_l(x)$ 是关于 $x$ 的 $l$ 阶多项式。
$$P_l(x)=\frac{1}{2^l}\sum _{s=0}^{[l/2]}\frac{(-1{)}^s(2l-2s)!}{s!(l-s)!(l-2s)!}{x}^{l-2s}$$
也可写作

这里列出了勒让德多项式的前6项结果，其他项可以根据公式自推。
也可参照勒让德多项式的基本介绍和MATLAB绘图
$$\begin{gathered} P_0(x)=1 \\ {P}_1(x)=x \\ {P}_2(x)=\frac{1}{2}(3x^2-1) \\ {P}_3(x)=\frac{1}{2}(5x^3-3x) \\ {P}_4(x)=\frac{1}{8}(35x^4-30x^2+3) \\ {P}_5(x)=\frac{1}{8}(63x^5-70x^3+15x) \end{gathered}$$

勒让德多项式也可以用罗德里格斯公式表示
$$P_l(x)=\frac{1}{2^{l}l!}\frac{d^l}{d{x}^l}(x^2-1{)}^l$$

连带勒让德多项式

连带勒让德方程来源于在球坐标系中使用分离变量法解拉普拉斯方程（球坐标系中的拉普拉斯方程）中
连带勒让德方程为
$$\frac{\mathrm{d}}{\mathrm{d}x}\left[\left(1-x^2\right)\frac{\mathrm{d}}{\mathrm{d}x}{P}_l^m(x)\right]+\left[l(l+1)-\frac{m^2}{1-x^2}\right]P_l^m(x)=0.$$

当$m=0$时，该方程变为勒让德方程。
$$\frac{\mathrm{d}}{\mathrm{d}x}\left[\left(1-x^2\right)\frac{\mathrm{d}}{\mathrm{d}x}{P}_l(x)\right]+l(l+1)P_l(x)=0$$

连带勒让德方程的解可用勒让德多项式$P_l(x)$生成
$$P_l^m=(-1{)}^m(1-x^2{)}^{m/2}\frac{d^m}{d{x}^m}{P}_l(x)$$

式中$(-1{)}^m$叫做 $Shortley-Condon$ 相位，通常只有在量子力学中会使用这个相位，不包含这个相位同样满足微分方程

将$P_l(x)$对于$x$的$m$次导数，记为$P_l^m(x)$。
$$P_l^m(x)=\Theta =(1-x^2{)}^{\frac{m}{2}}{P}_l^{[m]}(x)$$

• $P_l(x)$是连带勒让德函数，定义域为[-1,1]
• $l,m$分别是两带勒让德函数的阶(degree)和次(order)，$m\in (0,l)$
• $P_l^{[m]}(x)$是对$P_l^m(x)$求导。
• 广义勒让德多项式是复数，连带勒让德多项式是指实数部分。

连带勒让德多项式的性质

• 端点性、奇偶性、反对称性，求导；参照链接

• 特殊值：
  P m m ( x ) = { ( − 1 ) m ( 2 m − 1 ) ! ! ( 1 − x 2 ) m / 2 ( m > 0 ) 1 ( m = 0 ) P_{m}^{m}(x)=\left\{

  (−1)m(2m−1)!!(1−x2)m/2(m>0)1(m=0)" role="presentation" style="position: relative;">(−1)m(2m−1)!!(1−x2)m/21(m>0)(m=0)$$\begin{array}{ll}(-1{)}^m(2m-1)!!{\left(1-x^2\right)}^{m/2}& (m>0)\\ 1& (m=0)\end{array}$$
  \right. Pmm​(x)={(−1)m(2m−1)!!(1−x2)m/21​(m>0)(m=0)​

• 递归关系：
  

连带勒让德多项式的求解

根据式（12）~（14）可得：
$$p_m^m(x)=(-1{)}^m(2m-1)!!(1-x^2{)}^{\frac{m}{2}}，(m>0)$$

double pmm = 1.0;
if (m > 0) {
	double sign = (m % 2 == 0 ? 1 : -1);
	pmm = sign * DoubleFactorial(2 * m - 1) * pow(1 - x * x, m / 2.0);
}
1
2
3
4
5

$$p_{m+1}^m(x)=x(2m+1)p_m^m(x)$$

double pmm1 = x * (2 * m + 1) * pmm;
1

$$p_l^m=\frac{(2l-1)x{P}_{l-1}^m-(l+m-1)P_{l-2}^m}{l-m}$$

for (int n = m + 2; n <= l; n++) {
	double pmn = (x * (2 * n - 1) * pmm1 - (n + m - 1) * pmm) / (n - m);
	pmm = pmm1;
	pmm1 = pmn;
}
return pmm1;
1
2
3
4
5
6

根据上面公式即可得到任意关于 $P_l^m(x)$ 的函数值

将勒让德多项式代入球谐函数

将勒让德多项式代入如下球谐函数中，即可得到球坐标系下的解$Y(\theta ,\phi )$
$$Y(\theta ,\psi )=\Theta (\theta )\Phi (\phi )$$

这个解就是球谐函数
Y l m ( θ , φ ) = P l m ( c o s θ ) { s i n ( m φ ) , c o s ( m φ ) } = P l m ( c o s θ ⋅ e i m φ ) Y_l^m(\theta,\varphi) = P_l^m(cos\theta)\{sin(m\varphi),cos(m\varphi)\} = P_l^m(cos\theta \cdot e^{im\varphi}) Ylm​(θ,φ)=Plm​(cosθ){sin(mφ),cos(mφ)}=Plm​(cosθ⋅eimφ)

利用上面的连带勒让德函数的结论，可推导出前2阶球面谐波函数的表达式。

球谐函数的性质

正交性：基函数两两正交
完备性：线性叠加可以逼近任意一个球体函数
旋转不变性：只与$\theta ,\phi$有关，与半径无关。故球谐函数是旋转对称的。

球谐函数的使用

SH，球谐函数，归根到底只是一组基函数
有了基函数，就可以把任意一个函数，描述成几个基函数的加权和了。
例如
$$y\approx 0.1y_0+0.3y_1+0.8y_2+0.001y_3+...$$

球谐基函数构成了一个函数空间（Y0,Y1,Y2,…），对于给定球面函数$f(\theta ,\phi )$，我们可以将它投影到这个空间里去，求出其在此空间中的坐标（c0,c1,c2,…）,即系数。

$f(\theta ,\phi )$ 在球面S上的展开式为：
$$f(\theta ,\phi )=\sum _{l=0}^{\infty }\sum _{m=-l}^l{C}_l^m{Y}_l^m(\theta ,\phi )$$

其中 $Y_l^m(\theta ,\phi )$ 已知，即为球谐基函数。$f(\theta ,\phi )$为我们需要分解的函数。
需要我们计算的便是系数$C_l^m$
根据上式变换，得到系数计算公式：
$$C_l^m={\int }_{\Omega }f(s)Y_l^m(s)ds，\Omega 为球面积分$$
也就是频域上的系数$C_l^m$要求我们：求原函数 $f$与球谐函数 $Y_l^m$在球面上的乘积的积分。一般我们把这个求系数的过程叫做投影(Projection)。

在实际操作中，我们写程序时不可能会有对无穷级数进行储存和卷积的操作，一般展开项只能是有限项，也就是：
$$f(s)=\sum _{l=0}^{n-1}\sum _{m=-l}^l{C}_l^m{Y}_l^m(s)=\sum _{i=0}^{n^2}{c}_i{y}_i(s)$$

球谐函数的计算获取

double EvalSHSlow(int l, int m, double phi, double theta) {
		CHECK(l >= 0, "l must be at least 0.");
		CHECK(-l <= m && m <= l, "m must be between -l and l.");

		double kml = sqrt((2.0 * l + 1) * Factorial(l - abs(m)) /
			(4.0 * M_PI * Factorial(l + abs(m))));
		if (m > 0) {
			return sqrt(2.0) * kml * cos(m * phi) *
				EvalLegendrePolynomial(l, m, cos(theta));
		}
		else if (m < 0) {
			return sqrt(2.0) * kml * sin(-m * phi) *
				EvalLegendrePolynomial(l, -m, cos(theta));
		}
		else {
			return kml * EvalLegendrePolynomial(l, 0, cos(theta));
		}
	}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18