SciPy 科研实战教程

环境：Windows + Miniconda，Python 3.14，SciPy 1.17.x，VS Code。
本文所有示例均在 SciPy 1.17.x 真实 API 下验证通过（本机实测 1.18.0，覆盖范围内 API 与 1.17 完全一致）。凡标注"以官方文档为准"处，请以 docs.scipy.org 为准。

1. 简介与定位：SciPy 与 NumPy 的关系

一句话定位：NumPy 提供"数组 + 基础数值运算"，SciPy 在 NumPy 之上提供"科学计算算法库"（优化、积分、插值、信号、统计、稀疏矩阵等）。二者共用同一个 ndarray 类型，SciPy 函数基本都能接受 NumPy 数组。

import numpy as np   # 数组、基础线性代数
import scipy         # 高级算法（构建在 numpy 之上）

SciPy 子模块总览（与你的科研场景对应的重点已加粗）：

子模块用途计算化学 / 催化典型场景
scipy.constants物理常数与单位玻尔兹曼常数、eV/Hartree 换算、kBT
scipy.optimize拟合、寻根、最小化Arrhenius 活化能拟合、Langmuir 等温线、过渡态几何优化辅助
scipy.integrate数值积分、解 ODE反应动力学、Langmuir-Hinshelwood 速率方程
scipy.interpolate插值谱数据插值、势能面 (PES) 插值
scipy.linalg稠密线性代数分子振动 / 哈密顿量对角化
scipy.sparse稀疏矩阵大体系哈密顿量、格点模型
scipy.stats分布、假设检验、误差误差棒、置信区间、Arrhenius 图线性回归
scipy.signal峰检测、滤波、平滑XRD / 色谱 / DOS 寻峰、S-G 平滑
scipy.fft快速傅里叶变换谱分析、信号频谱
scipy.special特殊函数球谐函数、误差函数
scipy.ndimage多维图像处理电子密度网格平滑
scipy.spatial空间算法 / 距离原子坐标最近邻、凸包

2. 安装与版本检查

import scipy, numpy
print("scipy", scipy.__version__)
print("numpy", numpy.__version__)

Miniconda 下安装 / 升级：

conda install scipy            # 推荐：conda 会处理好 BLAS/LAPACK 依赖
# 或
pip install --upgrade scipy

版本检查的常用套路（确认某个函数在不在当前版本）：

from scipy import constants, optimize, integrate
print(hasattr(constants, "physical_constants"))   # True
print(hasattr(optimize, "curve_fit"))             # True

提示：SciPy 1.x 已用 scipy.integrate.solve_ivp 取代旧的 odeint、用 scipy.optimize.root 取代旧的 fsolve，写新代码时优先用新 API。

3. constants 物理常数模块

scipy.constants 收录 CODATA 国际推荐值，全部以 SI 基本单位 存储（J、kg、m、s、K）。计算化学里高频的是这几个：

from scipy import constants as C

# 玻尔兹曼常数（J/K，CODATA 2018 精确值）
print(C.Boltzmann)         # 1.380649e-23
print(C.k)                 # 同义别名

# 阿伏伽德罗常数（1/mol）
print(C.Avogadro)          # 6.02214076e+23
print(C.N_A)               # 同义别名

# 基本物理量
print(C.eV)                # 电子伏特 = 1.602176634e-19 J
print(C.e)                 # 元电荷 = 1.602176634e-19 C
print(C.h)                 # 普朗克常数 = 6.62607015e-34 J·s
print(C.hbar)              # 约化普朗克常数
print(C.c)                 # 光速 m/s
print(C.R)                 # 摩尔气体常数 = 8.314462618 J/(mol·K)
print(C.atomic_mass)       # 原子质量单位 u = 1.66053906892e-27 kg

3.1 能量单位换算（计算化学高频）

constants 里没有 C.hartree 这个属性。Hartree 能量要从 CODATA 字典里取：

from scipy import constants as C

Eh_J   = C.physical_constants["Hartree energy"][0]        # 4.359744722206e-18 J
Eh_eV  = C.physical_constants["Hartree energy in eV"][0]  # 27.211386245981 eV

# 常用换算系数（直接算出来，不必死记）
eV_to_J        = C.eV                          # 1 eV = 1.602176634e-19 J
J_to_eV        = 1.0 / C.eV                    # 1 J  = 6.2415e18 eV
Hartree_to_eV  = Eh_eV                         # 1 Ha = 27.211386245981 eV
Hartree_to_kJ  = Eh_J * C.N_A / 1000           # 1 Ha = 2625.4996 kJ/mol
eV_to_kJ_mol   = C.eV * C.N_A / 1000           # 1 eV = 96.4853 kJ/mol
kcal_to_eV     = 4184.0 / C.N_A / C.eV         # 1 kcal/mol = 0.043364 eV

print(f"1 Hartree = {Hartree_to_eV:.4f} eV = {Hartree_to_kJ:.2f} kJ/mol")
print(f"1 eV      = {eV_to_kJ_mol:.4f} kJ/mol")

3.2 室温热运动能 kBT（催化 / 动力学常用）

T = 298.15
kBT_eV = C.k * T / C.eV          # 0.0256926 eV ≈ 25.7 meV
kBT_J  = C.k * T                 # 4.116e-21 J
print(f"298.15 K 时 kBT = {kBT_eV*1000:.1f} meV")

记住这个量级：室温 kBT ≈ 25.7 meV，约 2.48 kJ/mol。判断吸附能 / 活化能是否跨越热扰动，全靠它。

3.3 按名字搜索常数

不确定一个常数叫什么时，用 C.find() 模糊搜索：

from scipy import constants as C
print(C.find("hartree")[:5])     # 列出所有含 "hartree" 的常数名
print(C.physical_constants["Boltzmann constant"])  # (值, 单位, 不确定度)

4. optimize：拟合、寻根、最小化

这是催化 / 动力学数据处理的核心模块。

4.1 curve_fit —— 重点掌握

curve_fit(f, xdata, ydata, p0=..., sigma=..., absolute_sigma=..., bounds=...) 用最小二乘拟合任意非线性模型，返回 (popt, pcov)：popt 是最优参数，pcov 是协方差矩阵（对角开方即参数标准误）。

完整签名（1.17）：

scipy.optimize.curve_fit(f, xdata, ydata, p0=None, sigma=None,
                         absolute_sigma=False, check_finite=None,
                         bounds=(-inf, inf), method=None, jac=None, ...)

示例 A：Arrhenius 方程拟合活化能

k = A \exp\left(-\frac{E_a}{RT}\right)

import numpy as np
from scipy.optimize import curve_fit

R = 8.314462618   # J/(mol·K)

# 1) 造一组"实验"数据：Ea = 40 kJ/mol, 指前因子 A = 1e7
T = np.array([300, 320, 340, 360, 380, 400])
k_true = 1e7 * np.exp(-40000 / (R * T))
rng = np.random.default_rng(42)
k = k_true * (1 + 0.03 * rng.normal(size=T.shape))   # 加 3% 相对噪声

# 2) 定义模型（参数顺序：A, Ea）
def arrhenius(T, A, Ea):
    return A * np.exp(-Ea / (R * T))

# 3) 拟合：p0 给一个合理的初值
popt, pcov = curve_fit(arrhenius, T, k, p0=[1e7, 40000])

A_fit, Ea_fit = popt
perr = np.sqrt(np.diag(pcov))          # 参数标准误
print(f"A  = {A_fit:.3e} ± {perr[0]:.3e}")
print(f"Ea = {Ea_fit:.1f} ± {perr[1]:.1f} J/mol = {Ea_fit/1000:.2f} kJ/mol")

示例 B：Langmuir 吸附等温线拟合

q = q_{\max}\frac{K P}{1 + K P}

# 1) 合成数据：qmax = 5 mmol/g, K = 0.02 (1/kPa)
P = np.linspace(0, 200, 20)
q_true = 5 * 0.02 * P / (1 + 0.02 * P)
q = q_true + 0.05 * rng.normal(size=P.shape)

def langmuir(P, qmax, K):
    return qmax * K * P / (1 + K * P)

# 2) 拟合，给 bounds 防止 qmax、K 变负（物理上必须 > 0）
popt, pcov = curve_fit(langmuir, P, q, p0=[5, 0.02],
                       bounds=(0, np.inf))
print("qmax =", popt[0], " K =", popt[1])

关键点：p0（初值）给错会收敛到错误解或直接发散。Langmuir 常用线性化初值：画 1/q 对 1/P 得 qmax、K 的粗估计（见第 13 节"常见坑"）。

示例 C：带权重拟合 + 置信区间

当每个数据点误差不同时，用 sigma 传误差、absolute_sigma=True 让协方差反映真实误差：

from scipy import stats

y_err = 0.05 * q          # 假设每点相对误差 5%
popt, pcov = curve_fit(langmuir, P, q, p0=[5, 0.02],
                       sigma=y_err, absolute_sigma=True)

# 95% 置信区间（用 t 分布，自由度 = n - 参数个数）
n = len(P); p = len(popt)
dof = n - p
t_val = stats.t.ppf(0.975, dof)
perr = np.sqrt(np.diag(pcov))
for name, v, e in zip(["qmax", "K"], popt, perr):
    print(f"{name} = {v:.4f} ± {t_val*e:.4f}  (95% CI)")

4.2 minimize —— 通用最小化（几何 / 参数寻优）

minimize(fun, x0, method=..., bounds=..., constraints=...) 返回 OptimizeResult，结果在 .x、.fun（目标函数值）、.success。

from scipy.optimize import minimize

# 例：用最小二乘目标函数拟合（等价于 curve_fit，但更灵活）
def residual(params):
    qmax, K = params
    return np.sum((langmuir(P, qmax, K) - q)**2)

res = minimize(residual, x0=[5, 0.02], method="Nelder-Mead")
print("最优参数:", res.x, " 目标值:", res.fun, " 成功:", res.success)

4.3 root —— 解方程 / 求根

解非线性方程（例如求反应平衡转化率、求交叉点）。推荐用 root，旧 fsolve 是它的包装：

from scipy.optimize import root

# 例：一级反应达到 90% 转化所需时间 t，解 1 - exp(-k t) = 0.9
k = 0.5
def f(t):
    return 1 - np.exp(-k * t) - 0.9

sol = root(f, x0=5.0)      # x0 是初值猜测
print("t =", sol.x[0], "s,  验证:", 1 - np.exp(-k * sol.x[0]))

4.4 least_squares —— 带边界的非线性最小二乘

比 curve_fit 更底层、更强（支持残差向量、边界、鲁棒损失）。拟合时若需要硬边界 + 更精细控制，用这个：

from scipy.optimize import least_squares

def residual_vec(params):          # 返回残差向量（不是平方和）
    qmax, K = params
    return langmuir(P, qmax, K) - q

res = least_squares(residual_vec, x0=[5, 0.02], bounds=(0, np.inf))
print("qmax, K =", res.x)

5. integrate：积分与常微分方程

5.1 quad —— 定积分

quad(func, a, b) 返回 (积分值, 误差估计)。适合算配分函数、归一化、期望值。

from scipy.integrate import quad

# 例：玻尔兹曼分布归一化 —— 积分 exp(-E/kBT) dE 从 0 到无穷
# 这里直接验证 ∫_0^∞ exp(-x) dx = 1
val, err = quad(lambda x: np.exp(-x), 0, np.inf)
print(val, err)   # 1.0 1e-9 量级

# 例：Maxwell-Boltzmann 速率分布的平均动能 <(1/2)mv^2>，略复杂可自行尝试
# 例：径向分布函数 g(r) 在壳层内的积分（配位数）
val, err = quad(lambda r: 4 * np.pi * r**2 * np.exp(-r), 0, np.inf)
print(val)        # 8π ≈ 25.13

quad 还能处理带参数的被积函数（args=）和奇异点（points=），见官方文档。

5.2 solve_ivp —— 解常微分方程（反应动力学）

scipy.integrate.solve_ivp(fun, t_span, y0, method='RK45', t_eval=..., rtol=..., atol=...) 返回 OdeResult：.t 是时间点，.y 是解数组，形状 (变量数, 时间点数)（注意第一维是变量，不是时间）。

这是新版推荐 API，取代旧的 odeint（odeint 已弃用，新代码请用 solve_ivp）。

示例：一级反应 A → B

\frac{d[A]}{dt} = -k[A]

from scipy.integrate import solve_ivp

k = 0.5                     # 速率常数
A0 = 1.0
t_eval = np.linspace(0, 10, 100)

def dA_dt(t, A):
    return -k * A

sol = solve_ivp(dA_dt, [0, 10], [A0], t_eval=t_eval)
t, A = sol.t, sol.y[0]      # y 的第 0 行是 [A]

import matplotlib.pyplot as plt
plt.plot(t, A, label="数值解")
plt.plot(t, A0 * np.exp(-k * t), "--", label="解析解")
plt.xlabel("t / s"); plt.ylabel("[A]")
plt.legend(); plt.show()

示例：二级反应 2A → products

\frac{d[A]}{dt} = -2k[A]^2

k2 = 0.1
def dA_dt2(t, A):
    return -2 * k2 * A**2

sol2 = solve_ivp(dA_dt2, [0, 50], [A0], t_eval=np.linspace(0, 50, 100))
# 解析解 1/[A] = 1/[A]0 + 2 k t，可自行验证
print("末浓度:", sol2.y[0, -1])

示例：Langmuir-Hinshelwood 速率方程（多组分耦合）

表面反应速率 r = k \frac{K C}{1 + K C}，组分消耗 \frac{dC}{dt} = -r：

k, K = 0.8, 0.05

def dC_dt(t, C):
    r = k * K * C / (1 + K * C)     # Langmuir-Hinshelwood 形式
    return -r

sol3 = solve_ivp(dC_dt, [0, 20], [10.0], t_eval=np.linspace(0, 20, 200))
print("末浓度:", sol3.y[0, -1])

solve_ivp 常用选项

选项含义
method'RK45'(默认，非刚性) / 'Radau' / 'BDF'(刚性)
t_eval指定输出时间点（不设则自动选点）
rtol / atol相对 / 绝对误差容限（精度关键参数）
dense_output=True返回连续可插值的解，可求任意时刻
events事件检测（如浓度降到阈值时停止）

刚性体系（快慢时间尺度差异大，如含快速平衡的机理）要用 method="BDF" 或 "Radau"，否则 RK45 会极慢。

6. interpolate：插值（谱数据、势能面）

6.1 interp1d（经典一维插值）

from scipy.interpolate import interp1d

# 例：实验谱只有离散点，插到统一波长网格上
wavelength = np.array([200, 300, 400, 500, 600])
absorbance = np.array([0.1, 0.4, 0.9, 0.6, 0.2])

f = interp1d(wavelength, absorbance, kind="cubic", fill_value="extrapolate")
w_new = np.linspace(200, 600, 400)
a_new = f(w_new)

interp1d 是经典旧接口，kind="cubic" 用的是分段三次（非自然样条）。需要真正的样条时用下面的 CubicSpline。

6.2 CubicSpline（三次样条，更推荐）

CubicSpline(x, y, bc_type="not-a-knot") 返回可调用对象，能求导 .derivative()、积分 .antiderivative()：

from scipy.interpolate import CubicSpline

cs = CubicSpline(wavelength, absorbance)      # 自然边界样条
print(cs(350.0))                               # 插值
print(cs.derivative()(350.0))                  # 一阶导（谱图一阶导可判峰位）

6.3 griddata（散点 → 规则网格，势能面插值）

计算势能面 (PES) 时，往往只有一批散点 (r1, r2, E)，要插成规则网格画等高线图。griddata(points, values, xi, method=...)：

from scipy.interpolate import griddata

# 散点：20 个 (r1, r2) 坐标对应的能量
rng = np.random.default_rng(0)
pts = rng.random((40, 2)) * 5
E = np.sin(pts[:, 0]) * np.cos(pts[:, 1])     # 假想的 PES

# 目标网格
gx = np.linspace(0, 5, 50)
gy = np.linspace(0, 5, 50)
GX, GY = np.meshgrid(gx, gy)
xi = np.column_stack([GX.ravel(), GY.ravel()])

E_grid = griddata(pts, E, xi, method="cubic").reshape(GX.shape)

plt.contourf(GX, GY, E_grid, levels=30)
plt.scatter(pts[:, 0], pts[:, 1], c="k", s=8)   # 原始散点
plt.colorbar(label="E"); plt.show()

method 可选 "linear" / "nearest" / "cubic"。规则网格（各维独立取样）用 RegularGridInterpolator 更快。

7. linalg 与 sparse：线性代数与稀疏矩阵

7.1 scipy.linalg 与 np.linalg 的分工

  • np.linalg：基础分解，够用但功能少。
  • scipy.linalg：功能更全（广义特征值 eigh、带状矩阵、Schur 分解等），且默认调 LAPACK，精度和速度更好。
import numpy as np
from scipy import linalg

# 例：分子振动 —— 质量加权 Hessian 矩阵对角化，本征值对应振动频率
H = np.array([[2.0, -0.5],
              [-0.5, 1.5]])          # 简化的 Hessian

w, v = linalg.eigh(H)                # 厄米矩阵用 eigh（Hessian 是对称的）
print("本征值(力常数):", w)
print("本征矢量(正则模式):\n", v)

# 解线性方程组 Ax = b（比 np.linalg.solve 更稳，支持更多类型）
A = np.array([[4., 1.], [1., 3.]])
b = np.array([1., 2.])
print(linalg.solve(A, b))

7.2 sparse —— 稀疏矩阵（大体系 / 格点模型）

当矩阵绝大多数元素为 0（如紧束缚 / 格点哈密顿量、大分子的邻接矩阵），用稀疏存储能省内存、加速。

from scipy import sparse
from scipy.sparse.linalg import spsolve, eigsh

# 1) 构造稀疏矩阵（COO 三列格式，最适合逐项填充）
rows = [0, 0, 1, 1, 2]
cols = [0, 1, 0, 1, 2]
data = [4.0, 1.0, 1.0, 3.0, 2.0]
A_sp = sparse.coo_matrix((data, (rows, cols))).tocsr()   # 转 CSR 便于运算

# 2) 稀疏求解 Ax = b
x = spsolve(A_sp, np.array([1.0, 2.0, 3.0]))
print(x)

# 3) 稀疏矩阵求最小本征值（大体系哈密顿量基态能量）
w_min, v_min = eigsh(A_sp, k=1, which="SA")   # SA = 最小代数本征值
print("最小本征值:", w_min[0])

常用稀疏格式：coo_matrix（构造）、csr_matrix（矩阵乘法、行切片快）、csc_matrix（列切片快）、lil_matrix（逐项改值快）。

8. stats：分布、拟合、假设检验、误差

8.1 常用分布对象

scipy.stats.norm 是正态分布对象，提供 pdf / cdf / ppf（分位数）/ rvs（采样）/ fit（估计参数）：

from scipy import stats

x = np.linspace(-4, 4, 200)
plt.plot(x, stats.norm.pdf(x, loc=0, scale=1), label="N(0,1)")
plt.plot(x, stats.norm.cdf(x, loc=0, scale=1), label="CDF")
plt.legend(); plt.show()

# 95% 双侧临界值
print(stats.norm.ppf(0.975))     # 1.96

# 从数据估计正态参数
data = np.array([2.1, 2.3, 1.9, 2.0, 2.2, 1.8])
loc, scale = stats.norm.fit(data)
print("均值估计:", loc, " 标准差估计:", scale)

8.2 线性回归（Arrhenius 图 / van't Hoff 图）

stats.linregress(x, y) 返回含 slope、intercept、rvalue（相关系数）、pvalue、stderr、intercept_stderr 的结果对象：

# Arrhenius 线性化：ln k = ln A - (Ea/R) * (1/T)
R = 8.314462618
T = np.array([300, 320, 340, 360, 380, 400])
k = 1e7 * np.exp(-40000 / (R * T))

res = stats.linregress(1 / T, np.log(k))
Ea = -res.slope * R          # 斜率 = -Ea/R
lnA = res.intercept
print(f"Ea = {Ea:.1f} J/mol, lnA = {lnA:.2f}, R² = {res.rvalue**2:.4f}")
print(f"斜率标准误 = {res.stderr:.3e}")

8.3 误差与置信区间（误差棒 / 误差分析）

设做了 n 次平行实验：

data = np.array([1.92, 2.03, 1.98, 2.01, 1.95, 2.00, 1.97, 2.04])
n = len(data)
mean = np.mean(data)
sd   = np.std(data, ddof=1)            # 样本标准差（ddof=1 无偏）
sem  = stats.sem(data)                 # 标准误 = sd / sqrt(n)

# 95% 置信区间（小样本用 t 分布，不是 1.96）
t_crit = stats.t.ppf(0.975, df=n - 1)
ci = t_crit * sem
print(f"均值 = {mean:.3f}, 标准差 = {sd:.3f}, 标准误 = {sem:.3f}")
print(f"95% 置信区间: [{mean - ci:.3f}, {mean + ci:.3f}]")

画误差棒：

plt.errorbar([1], [mean], yerr=[sem], fmt="o", capsize=5)
plt.ylim(0, 3); plt.show()

8.4 正态性检验

判断一组数据是否服从正态分布（决定用 t 检验还是非参数检验）：

# Shapiro-Wilk：小样本首选
W, p = stats.shapiro(data)
print("Shapiro p =", p, "(p>0.05 则不能拒绝正态假设)")

# 综合检验（偏度+峰度）
k2, p2 = stats.normaltest(data)
print("normaltest p =", p2)

其他常用：stats.ttest_ind（两独立样本 t 检验）、stats.ttest_rel（配对）、stats.kstest（分布拟合优度）。

9. signal：峰检测、滤波、平滑（谱图 / XRD / 色谱）

9.1 find_peaks —— 寻峰

find_peaks(x, height=..., threshold=..., distance=..., prominence=...) 返回 (峰索引数组, 属性字典)：

from scipy.signal import find_peaks

# 例：模拟双峰谱（如 XRD 两个衍射峰 / DOS 两个特征峰）
x = np.linspace(0, 10, 1000)
spectrum = (np.exp(-((x - 3)**2) / 0.2)
            + 0.5 * np.exp(-((x - 7)**2) / 0.15))

# height=最小峰高, distance=两峰最小间隔(点数)
peaks, props = find_peaks(spectrum, height=0.3, distance=50)

print("峰位置 x =", x[peaks])
print("峰高 =", props["peak_heights"])

plt.plot(x, spectrum)
plt.plot(x[peaks], spectrum[peaks], "ro")
plt.show()

常用属性：prominence（显著度，判断"峰"是否明显高于基线）、width（半高宽）、left_bases / right_bases。prominence 是排除噪声小凸起的利器。

9.2 savgol_filter —— Savitzky-Golay 平滑（保峰形）

对谱数据平滑同时尽量保留峰的高度和宽度（比简单移动平均好，不易把峰抹平）：

from scipy.signal import savgol_filter

# 加噪声的谱
noisy = spectrum + 0.02 * rng.normal(size=spectrum.shape)

# window_length 必须为奇数，polyorder < window_length
smooth = savgol_filter(noisy, window_length=51, polyorder=3)

plt.plot(x, noisy,  alpha=0.4, label="含噪")
plt.plot(x, smooth, label="S-G 平滑")
plt.legend(); plt.show()

deriv=1 / deriv=2 可直接求平滑后的一阶 / 二阶导（一阶导零点找峰位，二阶导找拐点）。

9.3 butter + filtfilt —— 数字滤波（去基线漂移 / 高频噪声）

from scipy.signal import butter, filtfilt

# 设计低通 Butterworth 滤波器
fs = 100.0                 # 采样频率 Hz
cutoff = 10.0              # 截止频率 Hz
b, a = butter(N=4, Wn=cutoff / (fs / 2), btype="low")  # Wn 归一化到 Nyquist

filtered = filtfilt(b, a, noisy)   # 零相位滤波（不引入相位畸变）

色谱 / 电化学数据常见基线漂移，可用 scipy.signal.detrend 去线性 / 常数基线。

10. fft：快速傅里叶变换（谱分析）

新代码用 scipy.fft（scipy.fftpack 是旧接口）。

from scipy.fft import rfft, rfftfreq

# 例：时域信号 = 5 Hz + 20 Hz 两个正弦叠加
fs = 100.0
t = np.arange(0, 2, 1 / fs)
signal = np.sin(2 * np.pi * 5 * t) + 0.6 * np.sin(2 * np.pi * 20 * t)

# 实数信号用 rfft（只算正频率一半，效率高）
spec = rfft(signal)
freq = rfftfreq(len(signal), d=1 / fs)
power = np.abs(spec)**2

peaks, _ = find_peaks(power, height=power.max() * 0.1)
print("检测到的频率(Hz):", freq[peaks])

plt.plot(freq, power)
plt.xlabel("Frequency / Hz"); plt.ylabel("Power")
plt.show()

要点：rfftfreq(n, d=采样间隔) 返回对应频率轴；next_fast_len(n) 可把长度补到 2 的幂提升 FFT 速度。

11. special：特殊函数（简要）

scipy.special 提供物理 / 数学特殊函数，计算化学常见：

from scipy import special

# 误差函数（高斯积分 / 扩散问题）
print(special.erf(1.0), special.erfc(1.0))

# 伽马函数与对数伽马（配分函数 / 阶乘，gammaln 更稳）
print(special.gamma(5), special.gammaln(100))

# 组合数 / 阶乘
print(special.comb(10, 3), special.factorial(10))

# 球谐函数（注意：SciPy 1.17 中是 sph_harm_y，不是 sph_harm）
# sph_harm_y(n, m, theta, phi)：n=角量子数, m=磁量子数, theta=极角, phi=方位角
Y = special.sph_harm_y(n=2, m=0, theta=0.5, phi=0.3)   # 返回复数
print("Y_2^0(0.5, 0.3) =", Y)

注：sph_harm_y 的相位约定与部分教材（如 Arfken）略有差异，含 Condon-Shortley 相位。做量子化学 / 轨道可视化时，请以 SciPy 官方文档的约定为准。

12. 科研实战综合示例

12.1 Arrhenius 活化能拟合（完整流程）

import numpy as np
from scipy.optimize import curve_fit
from scipy import stats

R = 8.314462618

# --- 1) 原始数据：不同温度 T 下的速率常数 k ---
T = np.array([298, 310, 325, 340, 355, 370])
k = np.array([0.012, 0.034, 0.11, 0.31, 0.82, 1.9])   # 单位 s^-1

# --- 2) 线性化初值估计（ln k vs 1/T）---
res = stats.linregress(1 / T, np.log(k))
Ea_guess = -res.slope * R
A_guess  = np.exp(res.intercept)
print(f"线性化初值: Ea = {Ea_guess:.0f} J/mol, A = {A_guess:.3e}")

# --- 3) 非线性拟合（用初值保证收敛）---
def arr(T, A, Ea):
    return A * np.exp(-Ea / (R * T))

popt, pcov = curve_fit(arr, T, k, p0=[A_guess, Ea_guess])
A, Ea = popt
perr = np.sqrt(np.diag(pcov))

# --- 4) 报告结果 ---
t_crit = stats.t.ppf(0.975, len(T) - 2)
print(f"指前因子 A = {A:.3e} ± {perr[0]:.3e}")
print(f"活化能   Ea = {Ea/1000:.2f} ± {t_crit*perr[1]/1000:.2f} kJ/mol")

# --- 5) 可视化 ---
T_fine = np.linspace(T.min(), T.max(), 100)
plt.plot(T, k, "o", label="实验")
plt.plot(T_fine, arr(T_fine, A, Ea), "-", label="拟合")
plt.xlabel("T / K"); plt.ylabel("k / s$^{-1}$"); plt.legend(); plt.show()

12.2 Langmuir 等温线拟合（完整流程）

# --- 1) 吸附等温线数据：平衡压力 P(kPa) vs 吸附量 q(mmol/g) ---
P = np.array([1, 2, 5, 10, 20, 50, 100, 150])
q = np.array([0.10, 0.19, 0.45, 0.83, 1.43, 2.50, 3.33, 3.85])

def langmuir(P, qmax, K):
    return qmax * K * P / (1 + K * P)

# --- 2) 线性化初值（1/q vs 1/P）---
res = stats.linregress(1 / P, 1 / q)
# 1/q = 1/qmax + 1/(qmax*K) * 1/P
qmax_guess = 1 / res.intercept
K_guess    = res.intercept / res.slope

popt, pcov = curve_fit(langmuir, P, q, p0=[qmax_guess, K_guess],
                       bounds=(0, np.inf))
qmax, K = popt
print(f"饱和吸附量 qmax = {qmax:.2f} mmol/g")
print(f"吸附平衡常数 K = {K:.4f} (1/kPa)")

P_fine = np.linspace(0, 150, 200)
plt.plot(P, q, "o", label="实验")
plt.plot(P_fine, langmuir(P_fine, qmax, K), "-", label="Langmuir 拟合")
plt.xlabel("P / kPa"); plt.ylabel("q / mmol·g$^{-1}$"); plt.legend(); plt.show()

12.3 一级 / 二级反应动力学数值积分

from scipy.integrate import solve_ivp

# 一级：d[A]/dt = -k1[A]
k1 = 0.4
sol1 = solve_ivp(lambda t, A: -k1 * A, [0, 20], [2.0],
                 t_eval=np.linspace(0, 20, 200))

# 二级：d[A]/dt = -2 k2 [A]^2
k2 = 0.15
sol2 = solve_ivp(lambda t, A: -2 * k2 * A**2, [0, 20], [2.0],
                 t_eval=np.linspace(0, 20, 200))

plt.plot(sol1.t, sol1.y[0], label="一级 k=%.2f" % k1)
plt.plot(sol2.t, sol2.y[0], label="二级 k=%.2f" % k2)
plt.xlabel("t / s"); plt.ylabel("[A]"); plt.legend(); plt.show()

12.4 DOS / 谱数据平滑与寻峰

from scipy.signal import find_peaks, savgol_filter

# --- 模拟态密度 DOS：三个展宽峰 + 噪声 ---
E = np.linspace(-10, 10, 2000)
dos = (2.0 * np.exp(-((E + 4)**2) / 0.8)      # 宽峰
       + 3.0 * np.exp(-((E - 1)**2) / 0.3)    # 尖峰
       + 1.5 * np.exp(-((E - 6)**2) / 1.0))   # 宽峰
dos_noisy = dos + 0.1 * rng.normal(size=dos.shape)

# --- 平滑 + 寻峰 ---
dos_smooth = savgol_filter(dos_noisy, window_length=41, polyorder=3)
peaks, props = find_peaks(dos_smooth, height=0.5, prominence=0.3)

print("峰位能量 E =", E[peaks])
print("峰高 =", props["peak_heights"])

plt.plot(E, dos_noisy, alpha=0.35, label="含噪 DOS")
plt.plot(E, dos_smooth, label="S-G 平滑")
plt.plot(E[peaks], dos_smooth[peaks], "ro", label="峰")
plt.xlabel("E / eV"); plt.ylabel("DOS"); plt.legend(); plt.show()

13. 常见坑（务必阅读）

坑 1：拟合初值选择不当

curve_fit 是局部优化，初值离真值太远会收敛到局部极小或发散。

  • 解法：先用线性化得到粗初值（Arrhenius 用 ln k 对 1/T，Langmuir 用 1/q 对 1/P），再非线性拟合。
  • 加 bounds 约束物理范围（速率常数、吸附量必须 ≥ 0）。
  • 试多组初值，看是否收敛到同一组参数。

坑 2：单位不一致

这是计算化学最常见的错误来源。 例：Arrhenius 里 R = 8.314 对应 Ea 用 J/mol；若活化能写成 kJ/mol 而不除以 1000，结果错 3 个数量级。eV、Hartree、kJ/mol、kcal/mol 混用必出错。建议：

  • 一个脚本里统一用 SI 单位，只在输出时换算。
  • 用第 3 节的换算系数显式转换，别心算。

坑 3：过拟合

多项式阶数太高、或模型参数太多，会"完美拟合噪声"，外推能力反而变差。

  • 解法：优先用有物理意义的模型（Arrhenius、Langmuir 等），而不是高阶多项式。
  • 报告参数的标准误 / 置信区间，看参数是否显著。
  • 用交叉验证 / 留一法评估泛化能力；比较模型时可用 AIC / BIC。

坑 4：solve_ivp 的 y 形状

solve_ivp 返回的 .y 形状是 (变量数, 时间点数)，第一维是变量不是时间。取第 i 个变量用 sol.y[i]，别写成 sol.y[:, i]（那是某一时刻所有变量）。

坑 5：把数值积分误差当真实误差

quad 返回的第二个值是数值误差估计，不是实验误差。拟合误差要看 pcov，不要混淆。

坑 6：find_peaks 默认会捡到一堆噪声小峰

find_peaks 不加 height / prominence 会把每个噪声凸起都当峰。务必设置 prominence（最有效）或 height，必要时配合 distance。

坑 7：用旧 API

新代码避免 scipy.integrate.odeint（用 solve_ivp）、scipy.optimize.fsolve（用 root）、scipy.interpolate.interp1d 的 cubic（用 CubicSpline）、scipy.fftpack（用 scipy.fft）。这些旧接口要么已弃用，要么功能弱于新接口。

14. 延伸资源

  • SciPy 官方文档：https://docs.scipy.org/doc/scipy/（最权威，含每个函数的完整签名、示例、注意事项）
  • 官方教程（SciPy Tutorial）：https://docs.scipy.org/doc/scipy/tutorial/（每个子模块有循序渐进的官方教程）
  • API 参考：https://docs.scipy.org/doc/scipy/reference/（遇到不确定的参数 / 签名，直接查这里）
  • 在 VS Code / 交互式环境里 help(scipy.optimize.curve_fit) 或 ?scipy.optimize.curve_fit 快速查签名。
  • 计算化学经典参考：
    • Atkins, Physical Chemistry（反应动力学、吸附等温线的理论背景）
    • NIST CODATA 常数表：https://physics.nist.gov/cuu/Constants/（scipy.constants 的数据来源）
  • 数据可视化配合 Matplotlib 官方文档：https://matplotlib.org/stable/

本文档中所有代码示例均可直接复制到 .py 文件运行；涉及画图的部分需先 import matplotlib.pyplot as plt。如有 API 细节疑问，一律以 SciPy 1.17 官方文档为准。