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

有些常微分方程（ODE）的解无法用初等函数表示。贝塞尔（Bessel）方程、勒让德（Legendre）方程、艾里（Airy）方程——它们都自然地出现在物理学中（如圆柱体内的热传导、行星引力场、量子隧穿）。这些方程的解本身定义了全新的函数。本章将教你如何借助幂级数求解，为何在奇点处必须使用弗罗贝尼乌斯（Frobenius）方法，以及为何同一组“特殊函数”会反复出现在物理与工程的各个角落。

本章学习目标#

• 常点处的幂级数解法及其收敛半径的几何意义
• 正则奇点与 Frobenius 方法
• 贝塞尔函数：鼓面振动与圆柱热传导
• 勒让德多项式：球谐函数与多极展开
• 埃尔米特（Hermite）多项式：量子谐振子
• 艾里函数：量子力学中的转折点
• 统一上述所有特殊函数的 Sturm-Liouville 框架

前置知识#

• 泰勒（Taylor）级数与收敛半径
• 第 1 至 3 章：二阶线性 ODE
• 复数基础（我们将借助复平面确定收敛半径）

为何需要幂级数？#

这就是著名的贝塞尔方程。我们之前学过的方法——常系数特征方程、积分因子、参数变易法——全都失效了，因为这里的系数是 $x$ 的函数，而非常数。更棘手的是，该方程在 $x = 0$ 处存在奇点，而这一点恰恰对应圆柱的轴线，物理上却必须有良好定义的解。

解决策略是：假设解具有级数形式，代入方程，让微分方程自己“生成”系数。为此，我们需要两种级数形式：

1. 在常点处使用普通泰勒级数 $\sum a_k (x - x_0)^k$ ；
2. 在正则奇点处使用 Frobenius 级数 $(x - x_0)^r \sum a_k (x - x_0)^k$ ，其中指数 $r$ 本身也由方程决定。

这套统一的机制能解锁贝塞尔、勒让德、埃尔米特、艾里、切比雪夫（Chebyshev）、拉盖尔（Laguerre）、超几何函数——整个经典特殊函数家族。

常点处的幂级数解#

方法#

若 $P(x)$ 和 $Q(x)$ 在点 $x_0$ 处解析，则称 $x_0$ 为常点。求解步骤如下：

1. 验证 $x_0$ 是否为常点；
2. 设 $y = \sum_{k=0}^{\infty} a_k (x - x_0)^k$ ；
3. 逐项求导后代入原方程；
4. 合并同类幂次 $(x - x_0)^k$ 的系数并令其为零；
5. 解出关于 $a_k$ 的递推关系，通常以 $a_0$ 和 $a_1$ 为自由参数（对应两个初始条件）。

福克斯（Fuchs）定理保证：所得级数在以 $x_0$ 为中心、且不包含 $P$ 或 $Q$ 任何奇点的最大开圆盘内收敛——收敛半径恰好延伸至复平面上最近的奇点。这正是为何即使处理实系数 ODE，复分析依然至关重要。

热身：重新“发现” $e^x$ #

幂级数方法“重新发现”了指数函数。关键在于：即便我们已知答案，这种递推过程仍是机械且可靠的。

艾里方程#

此方程出现在波（光波、量子概率幅、水波）遇到“局部波数从实数平滑过渡到虚数”的转折点时——典型例子是线性势垒附近的量子隧穿。

系数分裂为两条独立链，分别由 $a_0$ 和 $a_1$ 启动，最终得到两个线性无关的解，即艾里函数 $\text{Ai}(x)$ 与 $\text{Bi}(x)$ 。其中 $\text{Ai}(x)$ 在 $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()

正则奇点与 Frobenius 方法#

为何普通泰勒级数会失效#

容易验证，$y = \sqrt{x}$ 是它的一个解。但 $\sqrt{x}$ 在原点处切线垂直，并非解析函数，因此任何形式为 $\sum a_k x^k$ 的级数都无法再现它。我们必须在试探解中引入 $x$ 的“分数次幂”。

正则奇点的定义#

在 $x_0$ 处均解析，则称 $x_0$ 为正则奇点。直观地说，这意味着 $P$ 的奇异性至多是单极点，$Q$ 至多是二阶极点。

Frobenius 试探式#

其中 $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$ ，则根据根的差异，可构造两个线性无关解：

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

贝塞尔函数#

贝塞尔方程#

该方程在物理学中频繁出现，典型场景包括：

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

此处 $x = 0$ 是正则奇点，且 $p_0 = 1$ 、$q_0 = -n^2$ ，故指标方程 $r^2 - n^2 = 0$ 给出 $r = \pm 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$ 的正零点 $j_{n,1} < j_{n,2} < \dots$ 间距约为 $\pi$ ，直接决定了圆形膜的共振频率。

从零点到鼓面振动模态#

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

 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()

第二类贝塞尔函数#

当 $n$ 为非负整数时，若取负根 $r = -n$ ，递推关系会失效，此时必须引入对数项。由此得到的线性无关解称为第二类贝塞尔函数 $Y_n(x)$ 。它在原点附近的行为为 $Y_n(x) \sim -(2/\pi)(n-1)!\,(x/2)^{-n}$ ，因此在 $x = 0$ 处发散。对于定义域包含轴心的问题（如实心鼓、实心导线），物理上只能使用 $J_n$ ；而对于环形区域的问题（如同轴电缆），则需同时使用 $J_n$ 和 $Y_n$ 。

勒让德多项式#

勒让德方程#

每当在球坐标系中分离变量时，都会遇到此方程：它是拉普拉斯算子的角向部分，出现在氢原子波函数、引力与静电的多极展开、天线辐射方向图等问题中。此处 $x$ 实际代表 $\cos\theta$ ，因此自然定义域为 $[-1, 1]$ 。

点 $x = \pm 1$ 是正则奇点（因首项系数 $1 - x^2$ 在此处为零），而 $x = 0$ 是常点。

为何解是多项式？#

当 $n$ 为非负整数时，分子在 $k = n$ 处为零，级数截断为多项式。另一线性无关解仍为无穷级数，且在 $x = \pm 1$ 处发散。若要求解在球极处正则，则必须选择多项式解，从而自然导致参数 $n$ 的量子化。

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

罗德里格斯（Rodrigues）公式与正交性#

这正是电磁学与引力理论中多极展开的数学基础。

埃尔米特多项式与量子谐振子#

埃尔米特方程#

其中的高斯权重 $e^{-x^2}$ 正好对应后续波函数的包络。

量子谐振子#

由此自然导出两个重要物理结论：

• 能量量子化：仅当 $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$ ，这是其背后埃尔米特多项式的直接体现。

 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()

Sturm-Liouville 统一框架#

配合适当的边界条件。通过选择不同的 $p(x)$ 、$q(x)$ 、权函数 $w(x)$ 及区间，即可复现各类特殊函数：

每族函数都继承以下三大性质：

1. 实本征值——算子在权 $w$ 下对称；
2. 正交性——$\int y_m y_n\,w\,dx = 0$ （当 $m \neq n$ ）；
3. 完备性——任意“合理”函数均可展开为这些本征函数的广义傅里叶级数。

这正是为何“球谐展开”、“贝塞尔模态展开”、“埃尔米特态展开”感觉如此相似：它们本质上是同一个定理，只是权函数不同而已。

Python：调用特殊函数#

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

from scipy import special
import numpy as np

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

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

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

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

本文插图由脚本 scripts/figures/ode/05-power-series.py 生成。拉取最新代码后重新运行该脚本，即可将所有 PNG 图像同步生成到英文和中文资源目录中。

总结#

练习题#

基础题

1. 用幂级数解 $y'' + y = 0$ ，验证结果为 $\sin x$ 和 $\cos x$ 。
2. 证明 $x = 0$ 是方程 $2x^2 y'' + 3x y' - (1 + x) y = 0$ 的正则奇点，并求其指标方程。
3. 利用 Rodrigues 公式计算 $P_4(x)$ 和 $P_5(x)$ ，并验证递推关系 $(n+1) P_{n+1}(x) = (2n+1) x P_n(x) - n P_{n-1}(x)$ 。

进阶题

4. 对艾里方程 $y'' = x y$ 写出递推关系，并分别计算两个独立解的前 8 个非零系数。
5. 从 $J_n(x)$ 的级数定义出发，证明贝塞尔递推公式 $J_{n-1}(x) + J_{n+1}(x) = (2n/x) J_n(x)$ 。
6. 将函数 $f(r) = 1 - r^2$ 在区间 $[0, 1]$ 上展开为以 $J_0$ 为基的 Fourier-Bessel 级数，并数值验证前几项系数。

编程题

7. 从头实现 $J_n(x)$ 的级数计算，并与 scipy.special.jv 在 $x \in [0, 30]$ 、$n = 0, 1, 2$ 范围内进行对比。
8. 复现鼓面振动模态图，模态取 $(n, m) \in \{(0,1),(0,2),(1,1),(2,1),(3,1)\}$ ，并对其中一个模态制作一个完整周期的时间动画。
9. 绘制前 10 个量子谐振子的概率密度 $|\psi_n|^2$ ，叠加在抛物势上，并标出经典转折点 $x = \pm\sqrt{2n+1}$ 。

下一步#

下一章从单变量升级到多变量，研究 ODE 系统。一个核心观察是：任何高阶 ODE 都可以写成一阶系统。这让我们能用矩阵语言统一处理：特征值决定稳定性，特征向量决定模态。线性代数和微分方程在这里第一次真正握手。

参考文献#

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