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

微分方程的世界里，有些方程无法用初等函数表示其解，但这并不意味着我们束手无策。幂级数解法打开了一扇新的大门——通过将解展开成无穷级数，我们不仅能求解更广泛的方程，还催生了贝塞尔函数、勒让德多项式等在物理学中无处不在的特殊函数。本章将带你领略这片精妙的数学世界。

从一个简单问题开始

想象你正在研究热传导问题，需要求解一个圆柱体内的温度分布。分离变量后，你会遇到这样一个方程：

这就是著名的贝塞尔方程。尝试用之前学过的方法——分离变量、积分因子、常数变易法——都行不通。这类方程的系数不是常数，而是关于的函数，传统方法失效了。

幸运的是，幂级数解法能够拯救我们。

幂级数解法的核心思想

基本原理

核心假设：假设解可以表示为幂级数

求解步骤： 1. 将的幂级数代入微分方程 2. 将所有项写成同一幂次的形式 3. 令各幂次的系数为零 4. 得到系数之间的递推关系 5. 用初始条件确定等

回顾：熟悉的例子

让我们用幂级数方法重新求解一个简单方程：

假设，则。

代入方程：

调整求和指标，令左边：

比较的系数：

由知，递推得：

所以：

幂级数方法"重新发现"了指数函数！

常点与奇点

定义

考虑二阶线性方程：

将其化为标准形式：

其中，。

定义： - 如果和在处解析（可展开为收敛幂级数），则称为常点（ ordinary point） - 否则，称为奇点（ singular point）

常点处的幂级数解

定理：如果是常点，则方程在附近有两个线性无关的幂级数解：

且收敛半径至少等于从到最近奇点的距离。

例子：艾里方程

艾里方程（ Airy equation）：

这是量子力学中 WKB 近似、光学中彩虹理论的核心方程。

分析：，，都是整函数，所以是常点。

求解：设，则。

代入方程：

调整指标使的幂次相同。对第一个和，令；对第二个和，令：

分离项：

令各项系数为零： -: -: 这个递推关系意味着系数分成三组： - 以为首的组： - 以为首的组： - 以为首的组：计算几项：

通解：

其中和是艾里函数，这就是两个特殊函数的诞生！

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)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

# 艾里函数
ax1.plot(x, Ai, 'b-', linewidth=2, label='Ai(x)')
ax1.plot(x, Bi, 'r-', linewidth=2, label='Bi(x)')
ax1.axhline(y=0, color='k', linestyle='-', linewidth=0.5)
ax1.axvline(x=0, color='k', linestyle='-', linewidth=0.5)
ax1.set_xlabel('x', fontsize=12)
ax1.set_ylabel('y', fontsize=12)
ax1.set_title('Airy Functions: Solutions of y\'\' - xy = 0', fontsize=14, fontweight='bold')
ax1.legend(fontsize=11)
ax1.grid(True, alpha=0.3)
ax1.set_ylim(-0.6, 1.2)

# 艾里函数在量子隧穿中的应用
# 势能阶梯 V(x) = V0 * x (线性势)
x_qm = np.linspace(-5, 3, 500)
V = x_qm  # 线性势（归一化）
psi, _, _, _ = airy(x_qm)  # 波函数正比于 Ai(x)
psi_prob = np.abs(psi)**2  # 概率密度

ax2.fill_between(x_qm, 0, V, alpha=0.3, color='gray', label='Potential V(x)')
ax2.plot(x_qm, psi_prob * 5, 'b-', linewidth=2, label='|ψ(x)|² (scaled)')
ax2.axvline(x=0, color='r', linestyle='--', linewidth=1.5, label='Classical turning point')
ax2.set_xlabel('Position x', fontsize=12)
ax2.set_ylabel('Energy / Probability', fontsize=12)
ax2.set_title('Quantum Tunneling Near Linear Potential Barrier', fontsize=14, fontweight='bold')
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('airy_functions.png', dpi=150, bbox_inches='tight')
plt.show()

print("Ai(0) =", airy(0)[0])
print("Bi(0) =", airy(0)[2])

正则奇点与 Frobenius 方法

什么是正则奇点？

并非所有奇点都一样"坏"。有些奇点可以通过改进的幂级数方法处理。

定义：如果是奇点，但 $和$

在处解析，则称为正则奇点（ regular singular point）。

否则称为非正则奇点（ irregular singular point）。

Frobenius 方法

思想：在正则奇点处，尝试广义幂级数解：

其中可以是任意实数（甚至复数）。

步骤： 1. 将广义幂级数代入方程 2. 令最低次幂的系数为零，得到指标方程（ indicial equation） 3. 解指标方程求出的可能值 4. 根据的值建立递推关系

指标方程

对于标准形式的方程：

设，。

指标方程：

设两个根为（约定）。

三种情况

根据的性质，有三种情况：

情况一：不是整数两个线性无关解：

情况二：

情况三：是正整数

其中可能为零。

贝塞尔方程与贝塞尔函数

贝塞尔方程

这是物理学中最重要的方程之一，出现在： - 圆柱坐标系中的拉普拉斯方程 - 圆形鼓膜的振动 - 圆柱形波导的电磁波传播 - 热传导问题

分析是正则奇点。将方程化为标准形式：

所以，。

在处：，。

指标方程：

解得。

第一类贝塞尔函数

取（），设。

代入方程，经过繁琐但直接的计算，得到递推关系：

通常取，。

结果：

第二类贝塞尔函数（诺伊曼函数）

当不是整数时，也是解，且与线性无关。

当是整数时，，需要构造第二个解。

诺伊曼函数的定义：

当为整数时，用极限定义。

贝塞尔函数的性质

递推关系：$

导数关系：$

渐近行为：当：

当：

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import jv, yv

# 贝塞尔函数可视化
x = np.linspace(0.01, 20, 1000)

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 第一类贝塞尔函数
ax1 = axes[0, 0]
for n in range(5):
    ax1.plot(x, jv(n, x), linewidth=2, label=f'$J_{n}(x)$')
ax1.axhline(y=0, color='k', linestyle='-', linewidth=0.5)
ax1.set_xlabel('x', fontsize=12)
ax1.set_ylabel('$J_n(x)$', fontsize=12)
ax1.set_title('Bessel Functions of the First Kind', fontsize=14, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)
ax1.set_ylim(-0.5, 1.1)

# 第二类贝塞尔函数
ax2 = axes[0, 1]
for n in range(4):
    ax2.plot(x, yv(n, x), linewidth=2, label=f'$Y_{n}(x)$')
ax2.axhline(y=0, color='k', linestyle='-', linewidth=0.5)
ax2.set_xlabel('x', fontsize=12)
ax2.set_ylabel('$Y_n(x)$', fontsize=12)
ax2.set_title('Bessel Functions of the Second Kind (Neumann)', fontsize=14, fontweight='bold')
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)
ax2.set_ylim(-1.5, 0.6)

# 贝塞尔函数的零点（鼓膜振动模式的关键）
ax3 = axes[1, 0]
from scipy.special import jn_zeros
for n in range(4):
    zeros = jn_zeros(n, 5)
    ax3.scatter(zeros, [n]*5, s=100, label=f'$J_{n}$zeros')
    for z in zeros:
        ax3.axvline(x=z, color=f'C{n}', alpha=0.3, linestyle='--')
ax3.set_xlabel('x', fontsize=12)
ax3.set_ylabel('n (order)', fontsize=12)
ax3.set_title('Zeros of Bessel Functions (Drum Vibration Modes)', fontsize=14, fontweight='bold')
ax3.legend(fontsize=10)
ax3.grid(True, alpha=0.3)

# 圆形鼓膜的振动模式
ax4 = axes[1, 1]
theta = np.linspace(0, 2*np.pi, 100)
r = np.linspace(0, 1, 50)
R, Theta = np.meshgrid(r, theta)
X = R * np.cos(Theta)
Y = R * np.sin(Theta)

# 某个振动模式 (n=1, m=1)
n, m = 1, 1
zero_nm = jn_zeros(n, m)[-1]
Z = jv(n, zero_nm * R) * np.cos(n * Theta)

contour = ax4.contourf(X, Y, Z, levels=20, cmap='RdBu')
ax4.set_xlabel('x', fontsize=12)
ax4.set_ylabel('y', fontsize=12)
ax4.set_title(f'Circular Drum Vibration Mode (n={n}, m={m})', fontsize=14, fontweight='bold')
ax4.set_aspect('equal')
plt.colorbar(contour, ax=ax4)

plt.tight_layout()
plt.savefig('bessel_functions.png', dpi=150, bbox_inches='tight')
plt.show()

勒让德方程与勒让德多项式

勒让德方程

或者标准形式：

这在球坐标系中无处不在： - 氢原子波函数 - 地球引力场展开 - 天线辐射模式

分析是常点，是正则奇点。

在附近展开：设，代入方程得递推关系：

勒让德多项式

当为非负整数时，级数截断为多项式！

因为当时，，后续所有偶（或奇）次系数都为零。

勒让德多项式是这些截断的解，经过归一化使得。

前几个勒让德多项式：$

罗德里格公式

这个公式非常优雅，直接给出任意阶的勒让德多项式。

正交性

勒让德多项式在上关于权重函数正交：

物理意义：任何在上的"合理"函数都可以展开为勒让德级数

这在多极展开中至关重要！

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import legendre

# 勒让德多项式可视化
x = np.linspace(-1, 1, 500)

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 前几个勒让德多项式
ax1 = axes[0, 0]
for n in range(6):
    Pn = legendre(n)
    ax1.plot(x, Pn(x), linewidth=2, label=f'$P_{n}(x)$')
ax1.axhline(y=0, color='k', linestyle='-', linewidth=0.5)
ax1.axvline(x=0, color='k', linestyle='-', linewidth=0.5)
ax1.set_xlabel('x', fontsize=12)
ax1.set_ylabel('$P_n(x)$', fontsize=12)
ax1.set_title('Legendre Polynomials', fontsize=14, fontweight='bold')
ax1.legend(fontsize=10, loc='lower right')
ax1.grid(True, alpha=0.3)

# 正交性验证
ax2 = axes[0, 1]
n_max = 6
ortho_matrix = np.zeros((n_max, n_max))
for i in range(n_max):
    for j in range(n_max):
        Pi = legendre(i)
        Pj = legendre(j)
        integral, _ = np.polynomial.legendre.leggauss(50)
        # 数值积分
        from scipy.integrate import quad
        result, _ = quad(lambda t: Pi(t) * Pj(t), -1, 1)
        ortho_matrix[i, j] = result

im = ax2.imshow(ortho_matrix, cmap='coolwarm', vmin=-0.5, vmax=2.5)
ax2.set_xlabel('n', fontsize=12)
ax2.set_ylabel('m', fontsize=12)
ax2.set_title('Orthogonality Matrix$\\int_{-1}^{1} P_m P_n dx$', fontsize=14, fontweight='bold')
plt.colorbar(im, ax=ax2)
ax2.set_xticks(range(n_max))
ax2.set_yticks(range(n_max))

# 函数展开示例：用勒让德多项式逼近阶跃函数
ax3 = axes[1, 0]
def step_function(x):
    return np.where(x >= 0, 1, -1)

# 计算展开系数
def legendre_coefficients(f, n_terms):
    coeffs = []
    for n in range(n_terms):
        Pn = legendre(n)
        integral, _ = quad(lambda t: f(t) * Pn(t), -1, 1)
        c_n = (2*n + 1) / 2 * integral
        coeffs.append(c_n)
    return coeffs

# 逼近
n_terms_list = [3, 7, 15, 31]
ax3.plot(x, step_function(x), 'k-', linewidth=3, label='Step function')

for n_terms in n_terms_list:
    coeffs = legendre_coefficients(step_function, n_terms)
    y_approx = np.zeros_like(x)
    for n, c in enumerate(coeffs):
        Pn = legendre(n)
        y_approx += c * Pn(x)
    ax3.plot(x, y_approx, linewidth=1.5, label=f'{n_terms} terms', alpha=0.7)

ax3.set_xlabel('x', fontsize=12)
ax3.set_ylabel('y', fontsize=12)
ax3.set_title('Approximating Step Function with Legendre Series', fontsize=14, fontweight='bold')
ax3.legend(fontsize=10)
ax3.grid(True, alpha=0.3)
ax3.set_ylim(-1.5, 1.5)

# 球谐函数可视化（勒让德多项式的应用）
ax4 = axes[1, 1]
theta = np.linspace(0, np.pi, 100)
phi = np.linspace(0, 2*np.pi, 200)
Theta, Phi = np.meshgrid(theta, phi)

# 球谐函数 Y_2^0(θ, φ) ∝ P_2(cos θ)
from scipy.special import sph_harm
l, m = 2, 0
Y = sph_harm(m, l, Phi, Theta)
Y_real = np.real(Y)

# 转换到笛卡尔坐标
r = np.abs(Y_real)
X = r * np.sin(Theta) * np.cos(Phi)
Y_coord = r * np.sin(Theta) * np.sin(Phi)
Z = r * np.cos(Theta)

ax4.contourf(Phi[:, 0] * 180/np.pi, Theta[0, :] * 180/np.pi, Y_real.T, 
             levels=20, cmap='RdBu')
ax4.set_xlabel('φ (degrees)', fontsize=12)
ax4.set_ylabel('θ (degrees)', fontsize=12)
ax4.set_title(f'Spherical Harmonic$Y_{l}^{m}$', fontsize=14, fontweight='bold')

plt.tight_layout()
plt.savefig('legendre_polynomials.png', dpi=150, bbox_inches='tight')
plt.show()

厄米方程与厄米多项式

量子谐振子

量子力学中，一维谐振子的薛定谔方程化简后得到厄米方程：

当为非负整数时，有多项式解——厄米多项式。

厄米多项式的定义

罗德里格公式：

前几个厄米多项式：$

正交性

厄米多项式关于权重函数在上正交：

量子谐振子波函数

谐振子的能量本征态：

能级：这就是著名的能量量子化！

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import hermite
from math import factorial

# 厄米多项式和量子谐振子
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

x = np.linspace(-4, 4, 500)

# 厄米多项式
ax1 = axes[0, 0]
for n in range(5):
    Hn = hermite(n)
    ax1.plot(x, Hn(x), linewidth=2, label=f'$H_{n}(x)$')
ax1.axhline(y=0, color='k', linestyle='-', linewidth=0.5)
ax1.axvline(x=0, color='k', linestyle='-', linewidth=0.5)
ax1.set_xlabel('x', fontsize=12)
ax1.set_ylabel('$H_n(x)$', fontsize=12)
ax1.set_title('Hermite Polynomials', fontsize=14, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)
ax1.set_ylim(-50, 50)

# 量子谐振子波函数
ax2 = axes[0, 1]
# 设 ℏ = m = ω = 1（自然单位）
def psi_n(n, x):
    Hn = hermite(n)
    norm = 1 / np.sqrt(2**n * factorial(n)) * (1/np.pi)**0.25
    return norm * Hn(x) * np.exp(-x**2 / 2)

for n in range(5):
    psi = psi_n(n, x)
    # 沿 y 轴偏移以便观看
    ax2.plot(x, psi + n, linewidth=2, label=f'n={n},$E_n$={n+0.5}')
    ax2.fill_between(x, n, psi + n, alpha=0.3)
    ax2.axhline(y=n, color='gray', linestyle='--', linewidth=0.5)

ax2.set_xlabel('Position x', fontsize=12)
ax2.set_ylabel('$\\psi_n(x)$+ Energy offset', fontsize=12)
ax2.set_title('Quantum Harmonic Oscillator Wavefunctions', fontsize=14, fontweight='bold')
ax2.grid(True, alpha=0.3)

# 概率密度
ax3 = axes[1, 0]
for n in range(4):
    psi = psi_n(n, x)
    prob = np.abs(psi)**2
    ax3.plot(x, prob, linewidth=2, label=f'n={n}')

# 添加经典概率分布对比（ n=10 的近似）
n_classical = 10
A = np.sqrt(2*n_classical + 1)  # 经典振幅
x_classical = np.linspace(-A, A, 500)
prob_classical = 1 / (np.pi * np.sqrt(A**2 - x_classical**2 + 0.01))
ax3.plot(x_classical, prob_classical / 10, 'k--', linewidth=2, 
         label='Classical (scaled)', alpha=0.7)

ax3.set_xlabel('Position x', fontsize=12)
ax3.set_ylabel('$|\\psi_n(x)|^2$', fontsize=12)
ax3.set_title('Probability Densities', fontsize=14, fontweight='bold')
ax3.legend(fontsize=10)
ax3.grid(True, alpha=0.3)

# 能级图
ax4 = axes[1, 1]
n_levels = 8
for n in range(n_levels):
    E = n + 0.5
    ax4.hlines(E, 0.2, 0.8, linewidth=3, color=f'C{n % 10}')
    ax4.text(0.85, E, f'n={n},$E_{n}$={E}', fontsize=10, va='center')

# 画势能曲线
x_pot = np.linspace(-3, 3, 100)
V = 0.5 * x_pot**2
ax4.plot(x_pot + 0.5, V, 'k-', linewidth=2, label='$V(x) = \\frac{1}{2}x^2$')

ax4.set_xlim(-0.5, 1.5)
ax4.set_ylim(-0.5, 8)
ax4.set_xlabel('Position / Level', fontsize=12)
ax4.set_ylabel('Energy', fontsize=12)
ax4.set_title('Energy Levels of Quantum Harmonic Oscillator', fontsize=14, fontweight='bold')
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('hermite_quantum.png', dpi=150, bbox_inches='tight')
plt.show()

实际应用：热传导与傅里叶-贝塞尔级数

问题设定

考虑一个半径为的圆柱形金属棒，初始温度分布为，表面保持零度。求温度随时间的变化。

数学模型

热传导方程（圆柱坐标，假设温度只依赖于和）：

边界条件：初始条件：

分离变量

设，代入得：

时间部分：空间部分：令，这正是零阶贝塞尔方程！

解

通解：由于在发散，取。

边界条件要求。

设的零点为，则。

最终解：

其中由初始条件和贝塞尔函数的正交性确定。

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import jv, jn_zeros
from scipy.integrate import quad

# 圆柱热传导问题求解
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 参数
a = 1.0  # 圆柱半径
k = 0.1  # 热扩散系数
n_terms = 10  # 级数项数

# J_0 的零点
zeros = jn_zeros(0, n_terms)

# 初始温度分布
def f(r):
    return 1 - r**2  # 抛物线分布

# 计算傅里叶-贝塞尔系数
def compute_coefficients(f, a, zeros, n_terms):
    coeffs = []
    for n in range(n_terms):
        alpha_n = zeros[n]
        
        # 分子：∫_0^a r f(r) J_0(α_n r/a) dr
        def integrand_num(r):
            return r * f(r) * jv(0, alpha_n * r / a)
        numerator, _ = quad(integrand_num, 0, a)
        
        # 分母：∫_0^a r [J_0(α_n r/a)]^2 dr = (a^2/2) [J_1(α_n)]^2
        denominator = (a**2 / 2) * jv(1, alpha_n)**2
        
        coeffs.append(2 * numerator / (a**2 * jv(1, alpha_n)**2))
    
    return coeffs

coeffs = compute_coefficients(f, a, zeros, n_terms)

# 温度分布
def u(r, t, coeffs, zeros, a, k):
    result = np.zeros_like(r)
    for n, (A_n, alpha_n) in enumerate(zip(coeffs, zeros)):
        result += A_n * jv(0, alpha_n * r / a) * np.exp(-k * (alpha_n/a)**2 * t)
    return result

r = np.linspace(0, a, 100)

# 不同时刻的温度分布
ax1 = axes[0, 0]
times = [0, 1, 5, 10, 20, 50]
for t in times:
    temp = u(r, t, coeffs, zeros, a, k)
    ax1.plot(r, temp, linewidth=2, label=f't = {t}')

ax1.set_xlabel('Radius r', fontsize=12)
ax1.set_ylabel('Temperature u(r, t)', fontsize=12)
ax1.set_title('Temperature Distribution in Cylinder Over Time', fontsize=14, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)

# 中心温度随时间变化
ax2 = axes[0, 1]
t_array = np.linspace(0, 50, 500)
center_temp = [u(np.array([0.0]), t, coeffs, zeros, a, k)[0] for t in t_array]

ax2.plot(t_array, center_temp, 'b-', linewidth=2)
ax2.set_xlabel('Time t', fontsize=12)
ax2.set_ylabel('Center Temperature u(0, t)', fontsize=12)
ax2.set_title('Temperature at Center vs Time', fontsize=14, fontweight='bold')
ax2.grid(True, alpha=0.3)

# 对数坐标查看衰减
ax3 = axes[1, 0]
ax3.semilogy(t_array, center_temp, 'b-', linewidth=2)
# 理论上主要由第一项主导
theoretical_decay = coeffs[0] * np.exp(-k * (zeros[0]/a)**2 * t_array)
ax3.semilogy(t_array, theoretical_decay, 'r--', linewidth=2, label='First term only')
ax3.set_xlabel('Time t', fontsize=12)
ax3.set_ylabel('Center Temperature (log scale)', fontsize=12)
ax3.set_title('Exponential Decay of Temperature', fontsize=14, fontweight='bold')
ax3.legend(fontsize=10)
ax3.grid(True, alpha=0.3)

# 3D 热图
ax4 = axes[1, 1]
R, T = np.meshgrid(r, t_array[:100])
U = np.zeros_like(R)
for i, t in enumerate(t_array[:100]):
    U[i, :] = u(r, t, coeffs, zeros, a, k)

contour = ax4.contourf(R, T, U, levels=20, cmap='hot')
ax4.set_xlabel('Radius r', fontsize=12)
ax4.set_ylabel('Time t', fontsize=12)
ax4.set_title('Heat Diffusion in Cylinder (Contour)', fontsize=14, fontweight='bold')
plt.colorbar(contour, ax=ax4, label='Temperature')

plt.tight_layout()
plt.savefig('heat_conduction.png', dpi=150, bbox_inches='tight')
plt.show()

print("傅里叶-贝塞尔系数（前 5 项）:")
for n, c in enumerate(coeffs[:5]):
    print(f"  A_{n+1} = {c:.6f}")

生成函数与渐近展开

贝塞尔函数的生成函数

这个公式连接了贝塞尔函数与复分析，在调频信号分析中非常有用。

勒让德多项式的生成函数

物理意义：这正是两个点电荷之间势能的展开！

如果有一个点电荷在原点，另一个在，它们之间的距离为。

当时：

这就是多极展开的基础！

渐近展开

当时，贝塞尔函数的行为：

物理意义：远离源点时，柱面波表现为衰减的正弦波。

特殊函数之间的联系

统一的框架： Sturm-Liouville 理论

所有这些特殊函数都是Sturm-Liouville 问题的本征函数： | 特殊函数 | 方程 | 权重| 区间 | |---------|------|------------|------| ||| 1 || ||| 1 || ||||| ||||| |||||

共同特征

正交性：不同本征值对应的本征函数正交
完备性：可以展开任意"合理"函数
递推关系：相邻阶的函数有简单关系
生成函数：紧凑的生成公式
渐近行为：在边界处有可预测的行为

数值方法与实现

幂级数的数值计算

直接计算幂级数可能有数值问题： 1. 项数多时效率低 2. 交替级数可能有消去误差

更好的方法： - 使用递推关系 - 利用渐近公式（大参数） - 查表与插值

Python 中的特殊函数

from scipy import special

# 贝塞尔函数
x = 5.0
print(f"J_0({x}) = {special.jv(0, x)}")
print(f"J_1({x}) = {special.jv(1, x)}")
print(f"Y_0({x}) = {special.yv(0, x)}")

# 贝塞尔函数的零点
zeros = special.jn_zeros(0, 5)
print(f"J_0 的前 5 个零点: {zeros}")

# 勒让德多项式
print(f"P_3(0.5) = {special.eval_legendre(3, 0.5)}")

# 球谐函数
l, m = 2, 1
theta, phi = np.pi/4, np.pi/3
Y = special.sph_harm(m, l, phi, theta)
print(f"Y_{l}^{m}(θ, φ) = {Y}")

# 厄米多项式
print(f"H_4(1.0) = {special.eval_hermite(4, 1.0)}")

# 艾里函数
Ai, Aip, Bi, Bip = special.airy(2.0)
print(f"Ai(2) = {Ai}, Ai'(2) = {Aip}")

自己实现贝塞尔函数

def bessel_j_series(nu, x, n_terms=50):
    """
    用级数计算第一类贝塞尔函数 J_ν(x)
    
    J_ν(x) = Σ_{m=0}^∞ (-1)^m / (m! Γ(m+ν+1)) * (x/2)^(2m+ν)
    """
    from scipy.special import gamma
    
    result = 0.0
    for m in range(n_terms):
        term = ((-1)**m / (np.math.factorial(m) * gamma(m + nu + 1)) 
                * (x/2)**(2*m + nu))
        result += term
        if abs(term) < 1e-15:  # 收敛检验
            break
    return result

# 测试
x = 3.0
for nu in [0, 1, 2, 0.5]:
    my_result = bessel_j_series(nu, x)
    scipy_result = special.jv(nu, x)
    print(f"J_{nu}({x}): 我的实现 = {my_result:.10f}, "
          f"scipy = {scipy_result:.10f}, "
          f"误差 = {abs(my_result - scipy_result):.2e}")

超几何函数：特殊函数的统一

高斯超几何函数

其中是 Pochhammer 符号。

惊人事实：许多特殊函数都是超几何函数的特例！

合流超几何函数

与厄米多项式的关系：

总结：级数解法的威力

本章要点

幂级数解法：将解展开为，通过递推关系求系数
常点与奇点：
- 常点：系数解析，幂级数解存在
- 正则奇点： Frobenius 方法，广义幂级数解
指标方程：决定幂次的代数方程
重要特殊函数：
- 贝塞尔函数：圆柱问题
- 勒让德多项式：球面问题
- 厄米多项式：量子谐振子
- 艾里函数：转折点问题
正交性与完备性：用特殊函数展开任意函数
统一框架： Sturm-Liouville 理论和超几何函数

下一章预告

在《线性微分方程组》中，我们将： - 学习矩阵指数和基解矩阵 - 分析耦合系统的行为 - 研究相空间中的轨迹 - 应用到生态学、电路和控制理论

练习题

基础题

用幂级数方法求解，，，并验证答案是。
验证是方程的正则奇点，并求出指标方程。
证明满足递推关系。
用罗德里格公式计算和。
验证厄米多项式满足。

进阶题

对于艾里方程：
- 写出的递推关系
- 计算前 8 个非零系数
- 估计级数的收敛半径
证明勒让德多项式满足递推关系：$$

(n+1)P_{n+1}(x) = (2n+1)xP_n(x) - nP_{n-1}(x)$$
证明满足微分方程的步骤：
- 将级数代入贝塞尔方程
- 验证递推关系给出正确的系数
用傅里叶-贝塞尔级数展开在上，满足。计算前 3 项系数。

编程题

编写函数计算第一类贝塞尔函数的级数展开，并与scipy.special.jv比较精度。
可视化量子谐振子的前 10 个能级波函数，标注节点位置。
模拟圆形鼓膜的振动：
- 使用和贝塞尔函数
- 制作动画展示不同振动模式
- 找出前 5 个共振频率
实现多极展开：
- 给定电荷分布，计算其电势
- 用勒让德多项式展开到 - 比较精确解和级数近似

思考题

为什么物理学中会出现如此多的特殊函数？它们与坐标系的选择有什么关系？
如果贝塞尔方程的参数是复数，解会有什么变化？
Sturm-Liouville 理论如何统一这些看似不同的特殊函数？

参考资料

Arfken, G. B., Weber, H. J., & Harris, F. E. (2012). Mathematical Methods for Physicists. Academic Press.
Bender, C. M., & Orszag, S. A. (1999). Advanced Mathematical Methods for Scientists and Engineers. Springer.
Abramowitz, M., & Stegun, I. A. (1964). Handbook of Mathematical Functions. Dover.
NIST Digital Library of Mathematical Functions: https://dlmf.nist.gov/
Watson, G. N. (1995). A Treatise on the Theory of Bessel Functions. Cambridge University Press.

本文是《常微分方程的世界》系列的第 5 章。