【cg】【球谐光照】预备知识之拉普拉斯方程
前言
本文将从拉普拉斯方程的定义、表示和解来简单介绍下该方程的基本定义，为后面研究球谐函数做一下准备。
拉普拉斯方程
拉普拉斯方程(Laplace's equation)[1]是一种关于多元函数的二阶偏微分方程。
设三元函数\(u=f(x, y, z)\)满足拉普拉斯方程，则有
\[ \Delta u = \nabla \cdot \nabla u = \mathrm{div} \cdot \mathrm{grad} \ u = \frac{\partial^2 u}{\partial{x^2}} + \frac{\partial^2 u}{\partial{y^2}} + \frac{\partial^2 u}{\partial{z^2}} = 0 \tag{2.1} \]
前置知识
在公式\((2.1)\)中，\(\Delta\)便称为拉普拉斯算子。 \(\nabla\)为向量微分算子，即nabla算子。 \(grad\)表示数量场u的梯度，它的结果是一个向量场。 \(div\)表示向量场的散度，它的结果是一个标量场。
可以看出，理解该偏微分方程需要掌握一些前置知识。
微分方程
对于一元函数来说，微分方程是研究目标函数与其自变量，及其各阶导函数之间关系的方程。这类微分方程通常称为常微分方程。
\[ F(x, y, y', y'', y^{(n)}) = 0 \]
而对于多元函数来说，微分方程便是研究目标函数与其自变量，及对于每个自变量的各阶偏导数(包括混合偏导数)之间关系的方程。这类微分方程通常移为偏微分方程。
\[ F(x, y, z, u, u'_x, u'_y, u'_z, u''_{x^2}, u''{xy}, \ldots) = 0 \]
很显然拉普拉斯方程就是一个偏微分方程。
一般很难确定偏微分方程的解析解，通过使用分离变量等方法将求偏微分方程转换为求常微分方程的问题的处理。
多元函数微分法
多元函数微分法是研究关于多元函数微分相关的方法，包括偏导数、方向导数、梯度、多元复合函数求导法则等内容。
其中偏导数是指函数关于某一个自变量的导数，它表示该函数在该自变量表示折坐标轴上的导数。求偏导数只需要将其它自变量看作常量，使用一元函数微分法的规则求导即可。 同理，多元复合函数的求对某个自变量的偏导数时，也可以将其它自变量看作常量，使用一元复合函数的求导法则求出。
而方向导数是指函数沿某个方向\(\vec{l} = (\cos\alpha, \cos\beta, \cos\gamma)\)的导数，其值可由如下公式得到，可以理解为这个方向在各个坐标轴上投影的偏导数之和。
\[ \frac{\partial u}{\partial l} = \frac{\partial u}{\partial x} \cos\alpha + \frac{\partial u}{\partial y} \cos\beta + \frac{\partial u}{\partial z} \cos\gamma \]
梯度
梯度是一个向量，对于函数\(u = f(x, y, z)\)的某一个自变量\(P(x, y, z)\)，梯度是该自变量各个分量的偏导数组成的向量，即
\[ \nabla u = grad \ f(x, y, z) = \frac{\partial u}{\partial x} \boldsymbol{i} + \frac{\partial u}{\partial y} \boldsymbol{j} + \frac{\partial u}{\partial z} \boldsymbol{k} \Leftrightarrow ( \frac{\partial u}{\partial x}, \frac{\partial u}{\partial y}, \frac{\partial u}{\partial z} ) \]
所以它是一个从数量场到向量场的映射。
此时方向向量可表示为梯度向量与方向向量的内积，即
\[ \frac{\partial u}{\partial l} = \nabla u \cdot (\cos\alpha, \cos\beta, \cos\gamma) \]
由此可以看出梯度的方向即为该点方向导数取得最大值的方向，梯度向量的模长即为该点方向向量的最大值。
散度
对于某向量场\(A(x, y, z) = P(x, y, z) \boldsymbol{i} + Q(x, y, z) \boldsymbol{j} + R(x, y, z) \boldsymbol{k} \Leftrightarrow (P, Q, R)\)，则A的散度为
\[ div A = \frac{\partial P}{\partial{x}} + \frac{\partial Q}{\partial{y}} + \frac{\partial R}{\partial{z}} \]
散度可以看作是速度场在某位置的通量密度，即不可压缩流体在某点的源头强度。当div > 0时，表示流体从该点向外发散，此时该点称为源点。当div < 0时，表示流体向该点汇聚，此时该点称为汇点。而div = 0时表示流体在该点无源无汇，流进和流出的速率相同。若div处处为0，则称该向量场为无源场。
显然散度是一个从向量场到数量场的映射。
显然拉普拉斯方程就是一个数量场(原函数)的梯度(向量场)的散度(数量场)。
以上基本概念可参考同济版高数七、九、十一章，限于篇幅，不做赘述。
球面坐标系表示
公式\((2.1)\)是拉普拉斯方程的直角坐标系表示，但是我们更关心它的球面坐标系表示。
设\(f(x, y, z)\)的球坐标系表示为\(u(r, \theta, \varphi)\)，其中r为球面半径，\(\theta\)为极角(与正z轴的夹角)，\(\varphi\)为方位角(在xy平面上的投影与正x轴的夹角)。
那么将有球坐标系下的拉普拉斯方程如下。
\[ \Delta u = \frac{1}{r^2} \frac{\partial}{\partial r}(r^2 \frac{\partial u}{\partial r}) +\frac{1}{r^2 \sin\theta} \frac{\partial}{\partial \theta}(\sin\theta \frac{\partial u}{\partial \theta}) +\frac{1}{r^2 \sin^2\theta} \frac{\partial^2 u}{\partial \varphi^2} = 0 \tag{2.2} \]
我们可以通过多种方式得到这个结果。
使用多元复合函数的求导法则
对于\(f(x, y, z)\)在球坐标坐标系下有
\[ \left\{ \begin{array}{l} x = r\sin\theta \cos\varphi \\ y = r\sin\theta \sin\varphi \\ z = r\cos\theta \end{array} \right. \tag{2.3} \]
所以这是一个多元函数与多元函数的复合函数，利用多元复合函数的求导法则(同济高数9.4节)，我们可以轻松求出其关于\(r, \theta, \varphi\)的一阶和二阶偏导数。
对于x的一阶偏导数为
\[ \frac{\partial u}{\partial x(r, \theta, \varphi)} = \frac{\partial u}{\partial r} \frac{\partial r}{\partial x} + \frac{\partial u}{\partial \theta} \frac{\partial \theta}{\partial x} + \frac{\partial u}{\partial \varphi} \frac{\partial \varphi}{\partial x} \]
对于x的二阶偏导数为
\[ \frac{\partial^2 u}{\partial x^2} = \frac{\partial }{\partial x}(\frac{\partial u}{\partial x}) = \frac{\partial}{\partial x}( \frac{\partial u}{\partial r} \frac{\partial r}{\partial x} + \frac{\partial u}{\partial \theta} \frac{\partial \theta}{\partial x} + \frac{\partial u}{\partial \varphi} \frac{\partial \varphi}{\partial x} ) \]
这里只展开第一项，其他项类似
\[ \begin{array}{l} \frac{\partial}{\partial x}( \frac{\partial u}{\partial r} \frac{\partial r}{\partial x} ) &=& \frac{\partial^2 r}{\partial x^2} \frac{\partial u}{\partial r} + \frac{\partial r}{\partial x} \frac{\partial}{\partial x}(\frac{\partial u}{\partial r}) \\ &=& \frac{\partial^2 r}{\partial x^2} \frac{\partial u}{\partial r} + \frac{\partial r}{\partial x} ( \frac{\partial^2 u}{\partial r^2} \frac{\partial r}{\partial x} + \frac{\partial}{\partial \theta} (\frac{\partial u}{\partial r}) \frac{\partial \theta}{\partial x} + \frac{\partial}{\partial \varphi} (\frac{\partial u}{\partial r}) \frac{\partial \varphi}{\partial x} ) \\ &=& \frac{\partial^2 r}{\partial x^2} \frac{\partial u}{\partial r} + \frac{\partial r}{\partial x} \frac{\partial r}{\partial x} \frac{\partial^2 u}{\partial r^2} + \frac{\partial r}{\partial x} \frac{\partial \theta}{\partial x} \frac{\partial^2 u}{\partial \theta \partial r} + \frac{\partial r}{\partial x} \frac{\partial \varphi}{\partial x} \frac{\partial^2 u}{\partial \varphi \partial r} \end{array} \]
又由
\[ \left\{ \begin{array}{l} r &=& \sqrt{x^2+y^2+z^2} \\ \theta &=& \arctan(\sqrt{\frac{x^2+y^2}{z^2}}) \\ \varphi &=& \arctan(\frac{y}{x}) \end{array} \right. \]
可消去所有的\(r、\theta、\varphi\)关于\(x、y、z\)的一二阶偏导数，使得最终的结果只剩下关于\(r、\theta、\varphi\)的二阶偏导数的结果，如公式\((2.2)\)所示。
这里是通过建立直接坐标系与球坐标系互相转换的函数来做多元复合函数求导得出的结果，也可以通过先建立直接坐标系与柱坐标系的转换函数，再建立柱坐标系与球坐标系的转换函数的方式求出。
拉梅系数
拉梅系数是利用弧微分来计算与坐标系无关各数学物理量的通项，且可以很容易的转换为其他坐标系。
对于给定的曲线\(r = (q_1, q_2, q_3)\)，其弧微分为
\[ \mathrm{d} r = \frac{\partial r}{\partial q_1} \mathrm{d} q_1 + \frac{\partial r}{\partial q_2} \mathrm{d} q_2 + \frac{\partial r}{\partial q_3} \mathrm{d} q_3 \]
则拉梅系数的定义为
\[ H_i = |\frac{\partial r}{\partial q_i}| = \sqrt{ (\frac{\partial x}{\partial q_i})^2 +(\frac{\partial y}{\partial q_i})^2 +(\frac{\partial z}{\partial q_i})^2 } \tag{2.4} \]
那么由散度的定义可得拉普拉斯方程的拉梅系数表达式为
\[ \Delta u = \frac{1}{H_1H_2H_3}[ \frac{\partial}{\partial q_1} (\frac{H_2H_3}{H_1} \frac{u}{\partial q_1}) + \frac{\partial}{\partial q_2} (\frac{H_1H_3}{H_2} \frac{u}{\partial q_2}) + \frac{\partial}{\partial q_3} (\frac{H_1H_2}{H_3} \frac{u}{\partial q_3}) ] \tag{2.5} \]
将\((2.3)\)中的\(r、\theta、\varphi\)作为作为\(q_1、q_2、q_3\)代入\((2.4)\)求得其拉梅系数为
\[ \left\{ \begin{array}{l} H_1 = 1 \\ H_2 = r \\ H_3 = r \sin\theta \end{array} \right. \]
代入\((2.5)\)便可得拉普拉斯方程的球坐标表达式\((2.2)\)。
张量分析
个人理解，张量(tensor)是一种数学符号描述，它使得在不同的参考系(坐标系)下都可以按照某个相同的规则或者描述进行变换。
在张量分析我们使用爱因斯坦求和约定(Einstein summation convention)[2]来表示，这可以省掉很多求和符号，使公式更加简洁活泼。
对于三维空间中的坐标系\(\vec{r} = (x^1, x^2, x^3)\)，其向量的弧微分为
\[ \mathrm{d} \vec{r} = \frac{\partial \vec{r}}{\partial x^i} \mathrm{d} x^i = \vec{g_i} \mathrm{d} x^i \]
其中 \(\vec{g_i} = \frac{\partial \vec{r}}{\partial x^i}\)定义为协变基矢量
定义\(\vec{g}\)为所有\(g_i\)组成的矢量，即
\[ \vec{g} = (\vec{g_1}, \vec{g_2}, \vec{g_3}) \]
同时有一组与\(\vec{g_i}\)对偶的矢量\(\vec{g}^i\)，称为逆变基矢量，且有
\[ \vec{g}^j \cdot \vec{g_i} = \delta_i^j，(i,j = 1, 2, 3) \]
协变基矢量和逆变基矢量可以通过\([g_{ij}]\)矩阵互相转换。
\[ g_{ij} = \vec{g_i} \cdot \vec{g_i} = \begin{bmatrix} g_{11} & g_{12} & g_{13} \\ g_{21} & g_{22} & g_{23} \\ g_{31} & g_{32} & g_{33} \end{bmatrix} \]
其逆矩阵为
\[ g^{ij} = [g_{ij}]^{-1} \]
可通过如下公式互相转换。
\[ \left\{ \begin{array}{l} \vec{g_i} = g_{ij} \vec{g}^j \\ \vec{g}^i = g^{ij} \vec{g_j} \\ \end{array} \right. \]
当我们有了弧微分\(d\vec{r}\)，就可以写出弧长的表达式。
\[ \mathrm{d} s^2 = \mathrm{d} \vec{r} \cdot \mathrm{d} \vec{r} = g_{ij} \mathrm{d} x^i \mathrm{d} x^j \]
便可进一步求出\(g_{ij}\)、\(g^{ij}\)等，进一步求出目标坐标系的表达式。
对于拉普拉斯方程，其张量表达式为
\[ \Delta = \nabla^2 = \frac{1}{\sqrt{|g|}} \partial_{x^i} \cdot g^{ij} \sqrt{|g|} \partial_{x^j} \tag{2.6} \]
其中\(|g|\)为矩阵\([g_ij]\)的行列式。其中\(\partial_{x^i} = \frac{\partial}{\partial x^i}\)。
那么由\((2.3)\)可求得球坐标系下有
\[ \left\{ \begin{array}{l} \mathrm{d} \vec{r} &=& \mathrm{d} r \vec{e_r} + r \mathrm{d} \theta \vec{e_{\theta}} + r \sin\theta \mathrm{d} \varphi \vec{e_{\varphi}} \\ \mathrm{d} s^2 &=& \mathrm{d} r^2 + r^2 \mathrm{d} \theta^2 + r^2 \sin^2\theta \mathrm{d} \varphi^2 \end{array} \right. \]
所以有
\[ \left\{ \begin{array}{l} g_{ij} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & r^2 & 0 \\ 0 & 0 & r^2 \sin^2\theta \\ \end{bmatrix} \\ g^{ij} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & r^{-2} & 0 \\ 0 & 0 & r^{-2} \sin^{-2}\theta \\ \end{bmatrix} \\ |g| = r^4\sin^2\theta \end{array} \right. \]
将其代入\((2.6)\)，便可得公式\((2.2)\)。
变分原理
变分原理(variational principle)一开始提出是为了解决最速降曲线问题的，后来逐渐发展成求解极值问题的一种重要方法。
设\(I\)是关于\(f(y, y')\)满足如下关系的泛函。
\[ I = \int_a^b f(y, y') \mathrm{d} x \]
若I有极值，则\(f(y, y')\)满足如下方程，证明略。
\[ \frac{\partial f}{\partial y} - \frac{\mathrm{d}}{\mathrm{d} x}(\frac{\partial f}{\partial y'}) = 0 \]
该方程便是欧拉-拉格朗日方程(Euler-Lagrange equation)。
根据散度的定义来定义泛函I，并代入欧拉-拉格朗日方程，同样可以得到公式\((2.2)\)，此处限于篇幅不做展开。
拉普拉斯方程的解
拉普拉斯方程\((2.2)\)是一个多元的偏微分议程，一般使用分离变量的方法将该偏微分方程分解为几个常微分方程，然后求解其本征值问题。
分离半径和角度
设\(u(r, \theta, \varphi) = R(r)Y(\theta, \varphi)\)，代入\((2.2)\)，得
\[ \frac{Y}{r^2} \frac{\mathrm{d}}{\mathrm{d} r}(r^2 \frac{\mathrm{d} R}{\mathrm{d} r}) +\frac{R}{r^2 \sin\theta} \frac{\partial}{\partial \theta}(\sin\theta \frac{\partial Y}{\partial \theta}) +\frac{R}{r^2 \sin^2\theta} \frac{\partial^2 Y}{\partial \varphi^2} = 0 \tag{3.1} \]
方程两边同乘\(\frac{r^2}{YR}\)，得
\[ \frac{1}{R} \frac{\mathrm{d}}{\mathrm{d} r}(r^2 \frac{\mathrm{d} R}{\mathrm{d} r}) +\frac{1}{Y \sin\theta} \frac{\partial}{\partial \theta}(\sin\theta \frac{\partial Y}{\partial \theta}) +\frac{1}{Y \sin^2\theta} \frac{\partial^2 Y}{\partial \varphi^2} = 0 \tag{3.2} \]
即
\[ \frac{1}{R} \frac{\mathrm{d}}{\mathrm{d} r}(r^2 \frac{\mathrm{d} R}{\mathrm{d} r}) = -\frac{1}{Y \sin\theta} \frac{\partial}{\partial \theta}(\sin\theta \frac{\partial Y}{\partial \theta}) -\frac{1}{Y \sin^2\theta} \frac{\partial^2 Y}{\partial \varphi^2} = l(l+1) \tag{3.3} \]
其中\(l \in N\)(自然数0,1,2, ...)。此处\(l\)为该方程的一个本征值。
由此得到如下两个方程。
\[ \left\{ \begin{array}{lr} \frac{\mathrm{d}}{\mathrm{d} r}(r^2 \frac{\mathrm{d} R}{\mathrm{d} r}) - l(l+1)R = 0 & (3.4)\\ \frac{1}{\sin\theta} \frac{\partial}{\partial \theta}(\sin\theta \frac{\partial Y}{\partial \theta}) +\frac{1}{\sin^2\theta} \frac{\partial^2 Y}{\partial \varphi^2} + l(l+1)Y = 0 & (3.5) \end{array} \right. \]
其中\((3.4)\)是一个常微分方程，且是一个欧拉方程，可参考同济高数7.9节的方法，将自变量r换元为指数\(e^t\)求得结果为
\[ R(r) = Cr^l + \frac{D}{r^{l+1}}, l \in N \]
但其实我们并不怎么关心这个结果，我们更关心角度部分\(Y(\theta, \varphi)\)。
继续分离极角和方位角
\((3.5)\)显然还是一个偏微分方程，我们需要继续分离变量。
设\(Y(\theta, \varphi) = \Theta(\theta) \varPhi(\varphi)\)，代入\((3.5)\)，得
\[ \frac{\varPhi}{\sin\theta} \frac{\mathrm{d}}{\mathrm{d} \theta}(\sin\theta \frac{\mathrm{d} \Theta}{\mathrm{d} \theta}) +\frac{\Theta}{\sin^2\theta} \frac{\mathrm{d}^2 \varPhi}{\mathrm{d} \varphi^2} + l(l+1)\Theta\varPhi = 0 \tag{3.6} \]
方程两边同乘\(\frac{\sin^2\theta}{\Theta\varPhi}\)，得
\[ \frac{\sin\theta}{\Theta} \frac{\mathrm{d}}{\mathrm{d} \theta}(\sin\theta \frac{\mathrm{d} \Theta}{\mathrm{d} \theta}) +\frac{1}{\varPhi} \frac{\mathrm{d}^2 \varPhi}{\mathrm{d} \varphi^2} + l(l+1)\sin^2\theta = 0 \tag{3.7} \]
即
\[ \frac{\sin\theta}{\Theta} \frac{\mathrm{d}}{\mathrm{d} \theta}(\sin\theta \frac{\mathrm{d} \Theta}{\mathrm{d} \theta}) + l(l+1)\sin^2\theta = -\frac{1}{\varPhi} \frac{\mathrm{d}^2 \varPhi}{\mathrm{d} \varphi^2} = m^2 \tag{3.8} \]
其中\(m \in Z 且 |m| \leq l\)。
由此得到如下两个方程，这两个方程便都是常微分方程了。
\[ \left\{ \begin{array}{lr} \frac{\mathrm{d}^2 \varPhi}{\mathrm{d} \varphi^2} + m^2 \varPhi = 0 & (3.9) \\ \sin\theta \frac{\mathrm{d}}{\mathrm{d} \theta}(\sin\theta \frac{\mathrm{d} \Theta}{\mathrm{d} \theta}) + [l(l+1)\sin^2\theta - m^2] \Theta = 0 & (3.10) \end{array} \right. \]
求方位角项
对于自变量只有\(\varphi\)的常微分方程\((3.9)\)，显然其解为
\[ \varPhi(\varphi) = e^{im\varphi}, m \in Z \tag{3.11} \]
因为\(\varphi\)是以\(2\pi\)为周期的，所以m必须是整数。
求极角项
对于自变量只有\(\theta\)的常微分方程\((3.10)\)，进行如下换元。
令\(x = \cos\theta\)，并设\(y(x) = \Theta(\theta)\)，则有
\[ \left\{ \begin{array}{l} \frac{\mathrm{d} \Theta}{\mathrm{d} \theta} = \frac{\mathrm{d} y}{\mathrm{d} x} \frac{\mathrm{d} x}{\mathrm{d} \theta} = -\sin\theta \frac{\mathrm{d} y}{\mathrm{d} x} \\ \sin\theta \frac{\mathrm{d}}{\mathrm{d} \theta} (\sin\theta \frac{\mathrm{d} \Theta}{\mathrm{d} \theta}) = \sin^2\theta \frac{\mathrm{d}}{\mathrm{d} x}(\sin^2\theta \frac{\mathrm{d} y}{\mathrm{d} x}) = (1-x^2) \frac{\mathrm{d}}{\mathrm{d} x} [ (1-x^2) \frac{\mathrm{d} y}{\mathrm{d} x} ] \end{array} \right. \]
代入\((3.10)\)，得
\[ \frac{\mathrm{d}}{\mathrm{d} x} [ (1-x^2) \frac{\mathrm{d} y}{\mathrm{d} x} ] + [ l(l+1) - \frac{m^2}{1-x^2} ] y = 0 \tag{3.12} \]
方程\((3.12)\)便是伴随勒让德方程(Associated Legendre equation)[3]。
该方程仅当\(l \in N 且 m \in Z 且 |m| \leq l\)时，才在\(x \in [-1, 1]\)上有非奇异解，这个解称之为伴随勒让德多项式(ALP, Associated Legendre Polynomials)。对每一个满足条件的\(l、m\)都能确定一个伴随勒让德多项式记为\(P_l^m(x)\)。
所以有极角项的解为
\[ \Theta(\theta) = P_l^m(\cos\theta)，l \in N, m \in Z, |m| \leq l \tag{3.13} \]
至此，只要我们能求出\(P_l^m\)的解析式，便求出了拉普拉斯方程的解。
后记
关于伴随勒让德多项式，我们将在下一篇求解球谐函数时展开，此处只要认为是一系列(由l、m不同导致有多个)关于\(\cos\theta\)确定的多项式即可。