2. Python科学计算基础

在计算物理中，我们面对的往往不是单个标量，而是成千上万个离散点上的场量：一维的波形、二维的温度分布、三维的势场，甚至随时间演化的状态序列。Python 原生列表虽然灵活，但逐元素循环在解释器层面执行，速度慢、代码冗长。NumPy 的核心价值就在于把“对每个元素做同一件事”交给编译好的底层循环完成，这就是向量化。本节以一个热身任务贯穿：计算一维谐振子在若干离散位置上的势能与受力，并统计随机初相位下的平均能量。任务虽小，却覆盖了数组创建、索引、广播、数学函数、随机数与性能对比六个知识点。

一、数组创建与索引：把物理量装进 ndarray

先定义计算域。设谐振子劲度系数 k = 2.0 N/m，在 x ∈ [−3, 3] m 上取 N = 7 个等距采样点。用 np.linspace 生成坐标，用 np.zeros、np.ones、np.full 预分配结果数组：

import numpy as np

k = 2.0
x = np.linspace(-3.0, 3.0, 7)      # [-3. -2. -1.  0.  1.  2.  3.]
U = 0.5 * k * x**2                  # 势能，逐元素平方
F = -k * x                          # 恢复力，方向指向原点
mask = np.abs(x) <= 1.0             # 布尔索引：|x| ≤ 1 的区域
print(x[mask], U[mask])

这里体现了 NumPy 索引的三种常用形态：整数索引 x[0] 取端点、切片 x[1:-1] 取内部点、布尔索引 x[mask] 按物理条件筛选。注意 x**2 不是对数组对象求幂，而是对每个元素求平方，返回同形状新数组。切片返回的是视图而非副本，若后续要修改子数组又不想影响原数组，应显式调用 .copy()。这一点在迭代更新场量时极易出错。

二、广播机制：不同形状数组的自动对齐

广播解决的是“形状不同但语义可对齐”的运算。规则可概括为：从尾部维度开始逐维比较，维度相等或其中之一为 1 即可广播，否则报错。热身任务中，若要计算多个劲度系数下同一组位置的势能，可构造形状 (3, 1) 的 k 列向量与形状 (7,) 的 x 行向量相乘，结果自动扩展为 (3, 7)：

k_vec = np.array([1.0, 2.0, 4.0])[:, None]   # 形状 (3, 1)
U_mat = 0.5 * k_vec * x[None, :]**2           # 形状 (3, 7)
print(U_mat.shape)                            # (3, 7)

广播不复制数据，只在计算时逻辑扩展，因此内存开销小。但要警惕隐式扩展带来的误算：若把形状 (7,) 与 (3,) 直接相加会因尾部维度 7≠3 而报错，此时应明确用 [:, None] 或 reshape 表达意图。物理上，广播最适合“参数扫描 × 空间网格”这类笛卡尔积结构。

三、向量化运算与常用数学函数

NumPy 的 np.sin、np.exp、np.sqrt、np.log 等通用函数（ufunc）都按元素作用，并支持 out= 参数复用内存。计算总能量 E = U + K 时，若速度场 v 与 x 同形，直接 E = U + 0.5*m*v**2 即可，无需循环。聚合函数 np.sum、np.mean、np.max 可配合 axis 参数按行或按列统计，例如 U_mat.mean(axis=1) 得到每个 k 下的平均势能。

四、随机数生成：可复现的抽样

统计平均能量需要随机初相位。推荐使用新式接口 np.random.default_rng(seed)，它比旧的 np.random.seed 更隔离、更可控。设相位 φ 在 [0, 2π) 均匀分布，抽样 10000 次，计算 x = A·cos(φ) 处的势能均值：

rng = np.random.default_rng(2024)
phi = rng.uniform(0.0, 2*np.pi, size=10000)
A = 1.5
x_s = A * np.cos(phi)
U_s = 0.5 * k * x_s**2
print(U_s.mean())   # 理论值 k*A²/4 = 1.125

固定种子保证结果可复现，便于调试与教学对照。示例中均值应接近理论值 kA²/4 = 1.125，偏差随样本量增大而减小，这正是蒙特卡罗方法收敛性的直观体现。

五、性能对比：循环与向量化

下面用同一物理计算对比两种写法。对 10⁶ 个点逐元素求势能，纯 Python 循环与 NumPy 向量化的耗时差异通常在两个数量级以上：

import time
n = 1_000_000
xs = np.linspace(-3, 3, n)

t0 = time.perf_counter()
U_loop = [0.5 * k * xi**2 for xi in xs]
t1 = time.perf_counter()

t2 = time.perf_counter()
U_vec = 0.5 * k * xs**2
t3 = time.perf_counter()

print(f"循环 {(t1-t0)*1e3:.1f} ms, 向量化 {(t3-t2)*1e3:.1f} ms")

具体毫秒数依机器而异，但量级关系稳定。向量化快的根本原因是循环发生在 C 层、避免了 Python 对象装箱与解释开销。需要提醒的是，向量化并非万能：当计算存在强前后依赖（如逐时步迭代）或需要复杂条件分支时，可考虑 np.where、掩码运算，或退回到 Numba、Cython 等方案。

六、把热身任务串起来

将上述片段组合，可得到一份可直接运行的脚本：生成坐标网格 → 计算势能与受力 → 用布尔索引标出 |x| ≤ 1 的束缚区 → 用广播扫描三组 k → 用随机相位估计平均势能 → 用计时器对比循环与向量化。运行后应观察到：U_mat 形状为 (3, 7)，束缚区包含 x = −1、0、1 三个点，随机均值接近 1.125，向量化耗时显著低于循环。若均值偏离较大，先检查样本量与种子；若广播报错，先打印 .shape 定位维度不匹配。这个闭环虽小，却把 NumPy 的六项基础能力落到了同一个物理对象上，后续的数值积分、微分方程求解与场可视化都将建立在这一套数组思维之上。