常微分方程（五）：级数解法与特殊函数

有些 ODE 的解，根本写不成熟悉的初等函数。 Bessel 方程描述圆柱里的热传导和鼓面的振动，Legendre 方程出现在球坐标的每一处分离变量，Airy 方程刻画量子隧穿；它们的解定义了全新的"特殊函数"。本章给出找到这些解的统一方法——幂级数与 Frobenius 法——并解释为什么同一小撮特殊函数会反复出现在物理与工程之中。

本章要点

常点处的幂级数解，及其收敛半径的几何意义
正则奇点与 Frobenius 方法
Bessel 函数：振动鼓膜与圆柱热传导
Legendre 多项式：球谐函数与多极展开
Hermite 多项式：量子谐振子
Airy 函数：量子力学中的转折点
把上述所有特殊函数串起来的 Sturm-Liouville 框架

前置知识

Taylor 级数与收敛半径
第一至三章：二阶线性 ODE
基本的复数概念（我们要用复平面读出收敛半径）

1. 为什么需要幂级数？

$$x^2 y'' + x y' + (x^2 - n^2) y = 0,$$

也就是 Bessel 方程。前几章学过的招数——常系数特征方程、积分因子、参数变易——一个都用不上：系数是$x$的函数而不是常数；更糟糕的是，$x = 0$处方程奇异，但圆柱的轴恰恰落在那里，物理上必须给出有意义的解。

对策。 假设解是一个级数，代回方程，让 ODE 自己来逐项决定系数。我们需要两种风味：

在常点处用普通 Taylor 级数$\sum a_k (x - x_0)^k$。
在正则奇点处用 Frobenius 级数$(x - x_0)^r \sum a_k (x - x_0)^k$，指数$r$也由方程决定。

同一套机器可以解出 Bessel、Legendre、Hermite、Airy、Chebyshev、Laguerre、超几何函数——古典特殊函数的整片家族。

2. 常点处的幂级数解

2.1 算法

$$y'' + P(x)\,y' + Q(x)\,y = 0.$$

如果$P, Q$都在$x_0$处解析，那么$x_0$是常点。算法如下：

验证$x_0$是常点。
假设$y = \sum_{k=0}^{\infty} a_k (x - x_0)^k$。
逐项求导后代入 ODE。
对每个幂次$(x - x_0)^k$把系数置零。
解出$a_k$关于$a_0, a_1$的递推式（这两个自由系数对应两个初值条件）。

Fuchs 定理保证：得到的级数在以$x_0$为中心、不包含$P, Q$任何奇点的最大开圆盘内收敛——而且收敛半径正好等于复平面上到最近奇点的距离。这就是为什么即使 ODE 是实系数，也得用复变的眼光看它。

在常点$x_0$附近，幂级数解收敛于一个最大圆盘，其边界由系数函数$P, Q$在复平面上的最近奇点决定。例如$P(x) = 1/(1 + x^2)$的奇点在$x = \pm i$，所以围绕$x_0 = 0$的收敛半径恰好$R = 1$。

2.2 热身：重新发现$e^x$用幂级数解$y' = y$。设$y = \sum a_k x^k$，则$y' = \sum (k+1) a_{k+1} x^k$，匹配系数得

$$(k+1)\,a_{k+1} = a_k, \qquad a_k = \frac{a_0}{k!},$$$$y = a_0 \sum_{k=0}^{\infty} \frac{x^k}{k!} = a_0\,e^x.$$

级数法"重新发现"了指数函数。意义在于：哪怕已经知道答案，递推也是机械可执行、不会失败的程序。

2.3 Airy 方程

$$y'' - x\,y = 0.$$

凡是波（光、量子概率幅、水波）遇到"局部波数从实变虚"的转折点处，都会出现这个方程；最经典的例子是量子力学中线性势垒附近的隧穿。

$$a_{k+2} = \frac{a_{k-1}}{(k+1)(k+2)} \qquad (k \geq 1), \qquad a_2 = 0.$$

系数被劈成两条独立链，分别由$a_0, a_1$起始；得到的两个线性独立解就是 Airy 函数$\text{Ai}(x)$与$\text{Bi}(x)$。Ai 在$x > 0$（经典禁区）指数衰减，在$x < 0$（经典允许区）振荡——量子隧穿的图像清清楚楚。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import airy

x = np.linspace(-15, 5, 1000)
Ai, Aip, Bi, Bip = airy(x)

plt.figure(figsize=(10, 5))
plt.plot(x, Ai, 'b-', linewidth=2, label='Ai(x)')
plt.plot(x, Bi, 'r-', linewidth=2, label='Bi(x)')
plt.axhline(0, color='k', linewidth=0.5)
plt.title("Airy 函数：y'' - xy = 0 的解")
plt.legend(); plt.ylim(-0.6, 1.2)
plt.grid(True, alpha=0.3); plt.tight_layout()
plt.show()

3. 正则奇点与 Frobenius 方法

3.1 普通 Taylor 级数会失败

$$4 x^2 y'' + y = 0.$$

直接验证可知$y = \sqrt{x}$是它的一个解；但$\sqrt{x}$在原点切线竖直、不解析，任何形如$\sum a_k x^k$的级数都不可能复现它。必须让"分数次幂"进入待定式。

左：过原点的多项式在$x = 0$处斜率有限，永远无法匹配$\sqrt{x}$的尖点。右：Frobenius 待定式$y = x^{1/2}\sum a_k x^k$精确还原它。$4x^2 y'' + y = 0$的指标方程恰有重根$r = 1/2$。

3.2 正则奇点的判别

$$(x - x_0)\,P(x), \qquad (x - x_0)^2\,Q(x)$$

都在$x_0$解析，那么$x_0$叫正则奇点。直观地说：$P$至多是单极，$Q$至多是二阶极。

3.3 Frobenius 待定式

$$y = (x - x_0)^r \sum_{k=0}^{\infty} a_k (x - x_0)^k, \qquad a_0 \neq 0,$$$$r(r - 1) + p_0\,r + q_0 = 0,$$

其中$p_0 = \lim_{x \to x_0}(x-x_0) P(x), q_0 = \lim_{x \to x_0}(x-x_0)^2 Q(x)$。两根$r_1 \geq r_2$决定如何拼出两个独立解：

情形	两个独立解
$r_1 - r_2 \notin \mathbb{Z}$	两个 Frobenius 级数，各对应一个根
$r_1 = r_2$（重根）	一个 Frobenius 级数，加上一个含$\ln(x - x_0)$的解
$r_1 - r_2 \in \mathbb{Z}_{>0}$	一个 Frobenius 级数；第二解可能含对数项

更高次幂的系数匹配则给出$a_k$的递推，逐项求解即可。

4. Bessel 函数

4.1 Bessel 方程

$$x^2 y'' + x y' + (x^2 - n^2) y = 0.$$

这是圆柱物理里最重要的方程，至少出现在三处：

圆形鼓面振动——极坐标下亥姆霍兹方程$\nabla^2 u + k^2 u = 0$的径向部分。
圆柱波导——金属管内的电磁模式。
圆柱体内的热扩散——金属棒的瞬态冷却。

$x = 0$是正则奇点，$p_0 = 1, q_0 = -n^2$，指标方程$r^2 - n^2 = 0$给出$r = \pm n$。

4.2 第一类 Bessel 函数

$$\boxed{\; J_n(x) = \sum_{m=0}^{\infty} \frac{(-1)^m}{m!\,\Gamma(m + n + 1)} \left(\frac{x}{2}\right)^{2m+n}.\;}$$

后面会反复用到的三条性质：

递推关系：$\displaystyle J_{n-1}(x) + J_{n+1}(x) = \frac{2n}{x}\,J_n(x)$，可在不同阶之间穿梭。
大$x$渐近：$\displaystyle J_n(x) \approx \sqrt{\frac{2}{\pi x}}\cos\!\left(x - \frac{n\pi}{2} - \frac{\pi}{4}\right)$，故远端就是衰减余弦。
零点：$J_n$的正零点$j_{n,1} < j_{n,2} < \dots$间距约为$\pi$，量子化了圆膜的共振频率。

前五个 Bessel 函数$J_0, \dots, J_4$。标出的零点$j_{0,1} \approx 2.405, j_{0,2} \approx 5.520, \ldots$正是单位圆膜对称振动模态的径向波数。

4.3 从零点到鼓面模态

$$u_{n,m}(r,\varphi,t) = J_n(j_{n,m}\,r)\,\cos(n\varphi)\,\cos(c\,j_{n,m}\,t),$$

其中$j_{n,m}$是$J_n$的第$m$个正零点。整数$n$是角向节线条数，$m$是径向节圆个数（含边界）。

单位圆鼓四个本征模态的瞬态图。黑色虚线是节线（不动的位置）。$(0,1)$是基础"呼吸"模态；$(0,2)$多了一道径向节；$(1,1)$多了一道角向节；$(2,1)$多了两道角向节。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import jv, jn_zeros

x = np.linspace(0.01, 20, 1000)
plt.figure(figsize=(8, 4))
for n in range(5):
    plt.plot(x, jv(n, x), linewidth=2, label=f'$J_{n}(x)$')
plt.axhline(0, color='k', linewidth=0.5)
plt.title('第一类 Bessel 函数')
plt.legend(fontsize=9); plt.grid(True, alpha=0.3)
plt.tight_layout(); plt.show()

4.4 第二类 Bessel 函数

当$n$是非负整数时取负根$r = -n$，递推会失效，必须引入对数项；得到的线性独立伙伴叫第二类 Bessel 函数$Y_n(x)$。$Y_n$在原点附近行为是$Y_n(x) \sim -(2/\pi)(n-1)!(x/2)^{-n}$，会发散。所以：定义域包含轴心的问题（实心鼓、实心圆柱）只能取$J_n$；定义域是圆环域的问题（同轴电缆）则需$J_n$与$Y_n$都用上。

5. Legendre 多项式

5.1 Legendre 方程

$$(1 - x^2)\,y'' - 2x\,y' + n(n+1)\,y = 0.$$

凡是球坐标下做分离变量都会遇到它：拉普拉斯算子的角向部分、氢原子波函数、引力与静电的多极展开、天线辐射方向图。这里$x$扮演的是$\cos\theta$的角色，所以自然定义域是$[-1, 1]$。

$x = \pm 1$（首项系数$1 - x^2$的零点）是正则奇点，$x = 0$是常点。

5.2 为什么是多项式？

$$a_{k+2} = \frac{k(k+1) - n(n+1)}{(k+1)(k+2)}\,a_k.$$

当$n$取非负整数时，分子在$k = n$处归零，级数截断为多项式。另一支线性独立解仍是无穷级数，并且在$x = \pm 1$发散；要求解在球的两极正则，就把参数$n$量子化为非负整数。

按$P_n(1) = 1$归一化：

$n$	$P_n(x)$
0	$1$
1	$x$
2	$\frac{1}{2}(3x^2 - 1)$
3	$\frac{1}{2}(5x^3 - 3x)$
4	$\frac{1}{8}(35x^4 - 30x^2 + 3)$

5.3 Rodrigues 公式与正交性

$$P_n(x) = \frac{1}{2^n n!}\,\frac{d^n}{dx^n}(x^2 - 1)^n.$$$$\int_{-1}^{1} P_m(x)\,P_n(x)\,dx = \frac{2}{2n+1}\,\delta_{mn}.$$

所以$[-1, 1]$上的合理函数都能展成$f(x) = \sum_n c_n P_n(x)$，系数$c_n = \frac{2n+1}{2}\int_{-1}^{1} f(x) P_n(x)\,dx$——这就是电磁学和引力中多极展开的数学骨架。

Legendre 多项式$P_0, \dots, P_5$。归一化保证它们都过$(1, 1)$，奇偶性保证它们都过$(-1, (-1)^n)$。每个$P_n$在$(-1, 1)$内恰有$n$个简单零点——Gauss-Legendre 数值积分用的就是这些根。

6. Hermite 多项式与量子谐振子

6.1 Hermite 方程

$$y'' - 2x\,y' + 2n\,y = 0.$$$$H_0 = 1, \quad H_1 = 2x, \quad H_2 = 4x^2 - 2, \quad H_3 = 8x^3 - 12x, \ldots$$

正交性：$\int_{-\infty}^{\infty} H_m H_n e^{-x^2} dx = \sqrt{\pi}\,2^n n!\,\delta_{mn}$，权函数恰是高斯——一会儿就成为波函数的包络。

6.2 量子谐振子

$$-\tfrac{1}{2}\psi'' + \tfrac{1}{2} x^2 \psi = E\,\psi.$$$$\psi_n(x) = \frac{1}{\sqrt{2^n n!}\,\pi^{1/4}}\,H_n(x)\,e^{-x^2/2}, \qquad E_n = n + \tfrac{1}{2}\quad(\text{以 }\hbar\omega\text{ 为单位}).$$

两条物理结论自动落地：

能量量子化。 只有$n = 0, 1, 2, \ldots$给出可归一化解，所以能级是分立的。
零点能。 基态$E_0 = \tfrac{1}{2}\hbar\omega \neq 0$——抛物势中的量子粒子不可能完全静止。

量子谐振子前六个本征态$\psi_n$，画在各自能量$E_n = n + \tfrac{1}{2}$处，共同叠在抛物势$V(x) = \tfrac{1}{2} x^2$之上。$\psi_n$恰有$n$个节点，就是底层 Hermite 多项式留下的指纹。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import hermite
from math import factorial

x = np.linspace(-4, 4, 500)
plt.figure(figsize=(8, 6))
for n in range(5):
    Hn = hermite(n)
    norm = 1 / np.sqrt(2**n * factorial(n)) * (1/np.pi)**0.25
    psi = norm * Hn(x) * np.exp(-x**2 / 2)
    plt.plot(x, psi + n, linewidth=2, label=f'n={n}')
    plt.fill_between(x, n, psi + n, alpha=0.2)
plt.title('量子谐振子波函数')
plt.grid(True, alpha=0.3); plt.tight_layout(); plt.show()

7. Sturm-Liouville 统一框架

$$\frac{d}{dx}\!\left[p(x)\,\frac{dy}{dx}\right] + \bigl[q(x) + \lambda\,w(x)\bigr]\,y = 0,$$

配合恰当的边界条件。换一组$p, q, w$与区间，就换出一族特殊函数：

函数	$p(x)$	$w(x)$	区间
Legendre$P_n$	$1 - x^2$	$1$	$[-1, 1]$
Bessel$J_n$	$x$	$x$	$[0, a]$
Hermite$H_n$	$e^{-x^2}$	$e^{-x^2}$	$(-\infty, \infty)$
Laguerre$L_n$	$x e^{-x}$	$e^{-x}$	$[0, \infty)$
Chebyshev$T_n$	$\sqrt{1 - x^2}$	$1/\sqrt{1 - x^2}$	$[-1, 1]$

每一行都自动继承三件事：

本征值是实的——算子在权$w$下对称。
正交性——$\int y_m y_n\,w\,dx = 0$（$m \neq n$）。
完备性——一般函数都能展成本征函数的广义 Fourier 级数。

这就是为什么"展成球谐函数"“展成 Bessel 模态"“展成 Hermite 态"感觉如出一辙：它们确实是同一条定理换权重。

8. Python：调用特殊函数

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
from scipy import special
import numpy as np

# Bessel
print(f"J_0(5) = {special.jv(0, 5):.6f}")
print(f"前 5 个 J_0 零点: {special.jn_zeros(0, 5)}")

# Legendre
print(f"P_3(0.5) = {special.eval_legendre(3, 0.5):.6f}")

# Hermite
print(f"H_4(1.0) = {special.eval_hermite(4, 1.0):.1f}")

# Airy
Ai, Aip, Bi, Bip = special.airy(2.0)
print(f"Ai(2) = {Ai:.6f}")

本文所有插图由脚本 scripts/figures/ode/05-power-series.py 生成；拉取代码后重跑该脚本，会把全部 PNG 同时写入 EN 与 ZH 资源目录。

总结

概念	核心思想
常点幂级数	$y = \sum a_k x^k$，求$a_k$的递推
收敛半径	等于复平面上到$P, Q$最近奇点的距离
正则奇点	Frobenius 法：$y = x^r \sum a_k x^k$
指标方程	关于$r$的二次方程，定首项指数
Bessel 函数	圆柱几何；零点决定鼓膜共振频率
Legendre 多项式	球几何；多极展开的基底
Hermite 多项式	量子谐振子波函数
Sturm-Liouville	一个框架，多族正交本征函数

练习题

基础题

用幂级数解$y'' + y = 0$，验证得到$\sin x$与$\cos x$。
证明$x = 0$是$2x^2 y'' + 3x y' - (1 + x) y = 0$的正则奇点，并求指标方程。
用 Rodrigues 公式算出$P_4(x), P_5(x)$，并验证递推$(n+1) P_{n+1}(x) = (2n+1) x P_n(x) - n P_{n-1}(x)$。

进阶题

对 Airy 方程$y'' = xy$写出递推，分别给出两个独立解的前 8 个非零系数。
由$J_n$的级数定义出发，证明 Bessel 递推$J_{n-1}(x) + J_{n+1}(x) = (2n/x) J_n(x)$。
把$f(r) = 1 - r^2$在$[0, 1]$上展成$J_0$的 Fourier-Bessel 级数，并数值验证前几个系数。

编程题

自己实现$J_n(x)$的级数求和，与 scipy.special.jv 在$x \in [0, 30], n = 0, 1, 2$上对比。
复现鼓面图：把模态$(n, m) \in \{(0,1),(0,2),(1,1),(2,1),(3,1)\}$画出来，并对其中一个做一个周期的时间动画。
画出量子谐振子前 10 个概率密度$|\psi_n|^2$叠在抛物势上，并标出经典转折点$x = \pm\sqrt{2n+1}$。

参考资料

Arfken, Weber & Harris, Mathematical Methods for Physicists, Academic Press (2012)
Bender & Orszag, Advanced Mathematical Methods, Springer (1999)
NIST DLMF: dlmf.nist.gov

系列导航


上一章	第四章：拉普拉斯变换
当前	第五章：级数解法与特殊函数
下一章	第六章：线性微分方程组