计算物理（Computational Physics）是一门利用数值方法和计算机模拟来研究物理问题的交叉学科。它与理论物理和实验物理并列为现代物理学的三大支柱。计算物理的核心思想是：将物理定律（通常以微分方程、积分方程或统计模型的形式表达）转化为可在计算机上求解的离散化数值问题。
为什么我们需要计算物理？ 理论物理能提供精确解析解的问题其实很少，主要包括少数线性系统、对称性极高的简单模型以及经过大量简化的理想化问题。例如：
- 三体问题：三个天体在万有引力作用下的运动无法写出封闭解析解，但数值模拟可以高精度预测任意时刻的位置
- Navier-Stokes 方程：描述流体运动的基本方程，只有极少数简单流动有解析解，工程中完全依赖数值求解
- 量子多体问题：N 个粒子的波函数维度是 3N，当 N>10 时解析方法完全失效
- 极端条件实验：恒星内部温度（107 K）、压力（1016 Pa）无法在地球实验室重现，模拟成为唯一途径
| 年代 |
里程碑 |
贡献者 |
意义 |
| 1940s |
ENIAC 用于弹道计算，首次大规模数值模拟 |
von Neumann, Metropolis |
计算物理概念诞生 |
| 1953 |
Metropolis-Hastings 算法提出，蒙特卡洛方法诞生 |
Metropolis, Ulam |
随机模拟方法体系建立 |
| 1957 |
分子动力学（MD）首次模拟刚性球体 |
Alder, Wainwright |
微观动力学模拟开端 |
| 1960s |
有限元方法（FEM）工程应用兴起 |
Clough, Zienkiewicz |
复杂几何问题得到解决 |
| 1970s |
快速傅里叶变换（FFT）在谱方法中的应用 |
Cooley, Tukey |
谱精度计算成为可能 |
| 1980s |
密度泛函理论（DFT）实际应用 |
Hohenberg, Kohn, Sham |
第一性原理计算革命 |
| 1990s |
格子玻尔兹曼方法（LBM）用于流体模拟 |
Frisch, Hasslacher, Pomeau |
微观-宏观桥梁 |
| 2000s |
GPU 加速计算，大规模并行模拟 |
NVIDIA + 各研究组 |
计算能力爆炸 |
| 2010s |
AI+物理：神经网络力场、PINN |
Raissi, Perdikaris, Karniadakis |
数据驱动方法兴起 |
| 2020s |
Exascale 计算、量子计算模拟 |
多国联合项目 |
新纪元 |
┌──────────┐ ┌──────────┐ ┌──────────┐ ┌──────────┐ ┌──────────┐ ┌──────────┐
│ 物理问题 │───→│ 数学模型 │───→│ 数值离散 │───→│ 编程实现 │───→│ 结果分析 │───→│ 物理解释 │
└──────────┘ └──────────┘ └──────────┘ └──────────┘ └──────────┘ └──────────┘
↑ │ │ │
│ ▼ ▼ ▼
└────────── 实验数据 ──────────────────────────→ 验证与校准 ←─────────────┘
每个步骤的关键任务：
- 物理建模：确定系统自由度，选择控制方程（牛顿方程、薛定谔方程、Navier-Stokes 方程等），设定初始条件和边界条件
- 数学表述：将物理定律写为 ODE、PDE、积分方程或统计模型的数学形式
- 数值离散：选择有限差分、有限元、谱方法或蒙特卡洛等离散化策略
- 算法编程：在计算机上实现算法，考虑精度、稳定性、计算效率
- 结果后处理：可视化、误差分析、物理量提取、收敛性研究
- 验证确认：与解析解（若有）、实验数据或其他独立计算结果对比验证代码正确性
有限差分法是最直观、最容易实现的数值方法。它将导数近似为函数在相邻网格点上的差商。
设函数 f(x) 在其定义域内光滑，在 xi 处进行泰勒展开：
f(xi+h)=f(xi)+hf′(xi)+2h2f′′(xi)+6h3f′′′(xi)+24h4f(4)(xi)+⋯
f(xi−h)=f(xi)−hf′(xi)+2h2f′′(xi)−6h3f′′′(xi)+24h4f(4)(xi)−⋯
由此可得常用的差分近似：
| 类型 |
公式 |
截断误差 |
说明 |
| 向前差分（一阶导） |
f′(xi)≈hf(xi+h)−f(xi) |
O(h) |
一阶精度，最简单的近似 |
| 向后差分（一阶导） |
f′(xi)≈hf(xi)−f(xi−h) |
O(h) |
一阶精度，用于边界 |
| 中心差分（一阶导） |
f′(xi)≈2hf(xi+h)−f(xi−h) |
O(h2) |
二阶精度，最常用 |
| 中心差分（二阶导） |
f′′(xi)≈h2f(xi+h)−2f(xi)+f(xi−h) |
O(h2) |
二阶精度，用于扩散项 |
计算 f(x)=sin(x) 在 x=1.0 处的导数。精确值：f′(1.0)=cos(1.0)≈0.5403023058681398。
| 步长 h |
向前差分 |
绝对误差 |
中心差分 |
绝对误差 |
误差比（向前/中心） |
| 0.1 |
0.4974 |
4.29×10−2 |
0.5394 |
9.14×10−4 |
47 |
| 0.01 |
0.5361 |
4.23×10−3 |
0.5403 |
9.14×10−6 |
463 |
| 0.001 |
0.5399 |
4.21×10−4 |
0.5403 |
9.14×10−8 |
4606 |
| 0.0001 |
0.5402 |
4.20×10−5 |
0.5403 |
9.14×10−10 |
45937 |
可见中心差分的误差以 O(h2) 衰减（步长缩小 10 倍，误差缩小 100 倍），而向前差分的误差以 O(h) 衰减（步长缩小 10 倍，误差缩小 10 倍）。这就是为什么中心差分虽然计算量相同，但精度远高于前向/后向差分的原因。
实际计算中，不能无限缩小步长。当 h 太小，浮点数舍入误差会占主导地位。最佳步长近似为：
hopt≈(M3ϵmach)1/3
其中 ϵmach≈2.2×10−16（双精度浮点数），M 是 f′′′ 的上界。对于 sin(x)，M=1，计算得 hopt≈8.7×10−6。当 h 小于这个值时，舍入误差开始增长。
有限元法将求解域划分为有限个互不重叠的单元（如三角形、四边形、四面体），在每个单元上构造插值函数，通过变分原理或加权残量法建立代数方程组。
核心步骤详解：
- 几何离散化：将连续体划分为有限个小单元，单元之间通过节点连接
- 单元分析：在每个单元 e 上构造形函数 Ni(e)(x)，实现：
u(e)(x)=i=1∑neNi(e)(x)ui(e)
- 单元刚度矩阵：计算每个单元的局部矩阵
Kij(e)=∫Ωe∂x∂Ni∂x∂Njdx
- 全局组装：将所有单元的局部矩阵组装为全局矩阵 K
- 施加边界条件：修改 K 和 f 以满足 Dirichlet/Neumann 边界条件
- 求解方程：解 Ku=f
| 特性 |
有限差分法（FDM） |
有限元法（FEM） |
| 网格要求 |
结构化网格（规则排列） |
非结构化网格（任意形状） |
| 处理复杂几何 |
困难，需坐标变换 |
擅长，自适应网格 |
| 边界条件 |
较容易实现 |
需特殊处理 |
| 守恒性 |
取决于差分格式 |
天然自动满足 |
| 误差分析 |
可进行傅里叶分析 |
需借助能量范数 |
| 典型应用 |
流体、热传导（规则区域） |
结构力学、电磁学（复杂形状） |
| 软件实现 |
简单 |
复杂 |
谱方法用全局基函数展开解，在谱空间（如波数域）进行运算。对于光滑解可以达到谱精度——误差随基函数数量呈指数衰减，而非代数收敛。
例：对流方程的傅里叶谱解法
一维对流方程 ∂t∂u+c∂x∂u=0，周期边界条件：
对空间做傅里叶变换：
u^k(t)=∫02πu(x,t)e−ikxdx,k=−N/2,…,N/2−1
方程变为：
dtdu^k=−ikcu^k
这是一个常微分方程，可以直接积分得到 u^k(t)=u^k(0)e−ikct。最后通过逆傅里叶变换得到 u(x,t)。
精度对比（用不同方法求解 dxdu=0 的 L2 误差，网格数 N）：
| N |
二阶 FDM |
四阶 FDM |
谱方法 |
| 8 |
8.2×10−1 |
3.5×10−1 |
1.3×10−3 |
| 16 |
2.1×10−1 |
3.8×10−2 |
1.1×10−7 |
| 32 |
5.3×10−2 |
3.1×10−3 |
4.3×10−14 |
| 64 |
1.3×10−2 |
2.1×10−4 |
2.2×10−14 |
对于光滑解，谱方法仅用 16 个网格点就达到了 10−7 量级的精度，而 FDM 需要数千个点才能达到类似水平。
| 特性 |
有限差分法 (FDM) |
有限元法 (FEM) |
谱方法 |
| 网格类型 |
结构化网格 |
非结构化网格 |
无需网格（全局基函数） |
| 复杂几何 |
困难 |
擅长 |
困难 |
| 精度（对光滑解） |
代数收敛 O(N−p) |
代数收敛 O(N−p) |
指数收敛 |
| 精度（对非光滑解） |
一阶 |
可控 |
吉布斯现象 |
| 并行化 |
容易 |
中等 |
容易 |
| 典型应用 |
流体力学、热传导 |
结构力学、电磁 |
湍流直接模拟、气候 |
| 实现复杂度 |
低 |
高 |
中 |
考虑一阶常微分方程（ODE）初值问题：
dtdy=f(t,y),y(t0)=y0
物理中几乎所有随时间演化的过程都可以表示为 ODE 或 ODE 系统：单摆的振动（二阶ODE）、化学反应动力学（多变量耦合）、天体轨道（六维相空间）、RC 电路的充放电等。
欧拉方法是最简单的数值方法，用前向差分近似时间导数：
yn+1=yn+hf(tn,yn)
其中 h=tn+1−tn 为时间步长。
数值实验：解 y′=−2y，y(0)=1（精确解为 y(t)=e−2t）：
| t |
精确解 y(t) |
欧拉 h=0.1 |
误差 |
欧拉 h=0.05 |
误差 |
| 0.0 |
1.0000 |
1.0000 |
0 |
1.0000 |
0 |
| 0.1 |
0.8187 |
0.8000 |
1.87×10−2 |
0.8100 |
8.73×10−3 |
| 0.2 |
0.6703 |
0.6400 |
3.03×10−2 |
0.6561 |
1.42×10−2 |
| 0.5 |
0.3679 |
0.3277 |
4.02×10−2 |
0.3487 |
1.92×10−2 |
| 1.0 |
0.1353 |
0.1074 |
2.80×10−2 |
0.1216 |
1.38×10−2 |
可见：步长减半，误差也大致减半——这正是 O(h) 收敛的典型表现。
Rung-Kutta 方法通过在一个时间步内多次估算 f 的中间值来提高精度。最经典的 RK4 每步需要 4 次函数求值：
k1k2k3k4yn+1=f(tn,yn)=f(tn+h/2,yn+hk1/2)=f(tn+h/2,yn+hk2/2)=f(tn+h,yn+hk3)=yn+6h(k1+2k2+2k3+k4)
RK4 求解 y′=−2y，h=0.2：
| t |
精确值 |
RK4 近似 |
误差 |
| 0.0 |
1.000000 |
1.000000 |
0 |
| 0.2 |
0.670320 |
0.670323 |
3.1×10−6 |
| 0.4 |
0.449329 |
0.449330 |
1.4×10−6 |
| 0.6 |
0.301194 |
0.301195 |
9.3×10−7 |
| 0.8 |
0.201897 |
0.201897 |
6.2×10−7 |
| 1.0 |
0.135335 |
0.135335 |
4.2×10−7 |
RK4 的全局误差 O(h4)：步长放大到 0.2 仍然达到了 10−6 精度，而欧拉方法在 h=0.05 时还有 10−2 的误差。这说明 RK4 的效率远高于欧拉方法。
| 方法 |
精度阶数 |
每步函数求值次数 |
稳定性限制 |
推荐用途 |
| 欧拉向前 |
O(h) |
1 |
∣1+hλ∣<1 |
教学示例 |
| 改进欧拉（Heun） |
O(h2) |
2 |
中等 |
非刚性问题 |
| RK4 |
O(h4) |
4 |
中等 |
通用首选 |
| Adams-Bashforth 4步 |
O(h4) |
1（每步） |
中等 |
光滑长积分 |
| 向后欧拉（隐式） |
O(h) |
需要解方程 |
无条件稳定 |
刚性方程 |
| TR-BDF2 |
O(h2) |
需要解方程 |
强稳定 |
电路/化学反应 |
刚性方程指同时存在极快和极慢时间尺度的系统。例如化学反应动力学：
dtdy=−1000y+999sin(t),y(0)=0
解析解：y(t)≈sin(t)−0.001cos(t)+0.001e−1000t
快变部分 e−1000t 在 t>0.01 后就衰减到可忽略，但显式方法为了稳定性需要步长满足 h<2/1000=0.002。这意味着我们需要用 h=0.002 的步长计算到 t=100，共 50000 步——而解本身在缓慢变化（以周期 2π 变化），大部分计算是在"绕远路"。
解决方案：使用隐式方法。向后欧拉法 yn+1=yn+hf(tn+1,yn+1) 对线性问题无条件稳定，允许使用与解变化时间尺度匹配的大步长。
偏微分方程（PDE）是物理建模的核心工具。按数学性质，PDE 分为三类：
| 类型 |
典型方程 |
物理意义 |
特征速度 |
| 抛物型 |
热传导 ∂tu=α∇2u |
扩散过程，信息无限传播 |
无限（指数衰减） |
| 双曲型 |
波动 ∂t2u=c2∇2u |
波动传播，有限特征速度 |
c |
| 椭圆型 |
泊松 ∇2u=−ρ/ε |
稳态场，全局耦合 |
- |
一维热传导方程：
∂t∂u=α∂x2∂2u
其中 α=k/(ρcp) 是热扩散系数。
最直接的离散化：时间用前向差分，空间用中心差分。
ujn+1=ujn+Δx2αΔt(uj+1n−2ujn+uj−1n)
定义网格比 r=αΔt/Δx2，则：
ujn+1=(1−2r)ujn+r(uj+1n+uj−1n)
稳定性分析：von Neumann 稳定性条件要求 ∣1−4rsin2(kΔx/2)∣≤1，即 r≤1/2。
数值算例：初始温度 u(x,0)=sin(πx)，边界 u(0,t)=u(1,t)=0，α=0.01：
| Δx |
Δt |
r=αΔt/Δx2 |
是否满足 CFL |
结果 |
| 0.05 |
0.10 |
0.40 |
✅ |
正确演化 |
| 0.05 |
0.15 |
0.60 |
❌ |
数值发散 |
| 0.02 |
0.02 |
0.50 |
✅ 临界 |
正确但有振荡 |
| 0.02 |
0.01 |
0.25 |
✅ |
正确光滑 |
若 r>0.5，高频模式会指数增长，导致计算立即发散。
将时间层的显式和隐式部分平均，达到 O(h2) 精度和无条件稳定：
ujn+1−2r(uj+1n+1−2ujn+1+uj−1n+1)=ujn+2r(uj+1n−2ujn+uj−1n)
每步需要解一个三对角线性方程组（Thomas 算法，O(N) 计算量），但允许使用任意大的时间步长——虽然大时间步长会牺牲精度。
一维波动方程描述弦振动、声波传播等：
∂t2∂2u=c2∂x2∂2u
标准离散格式（中心差分对时间和空间）：
ujn+1=2ujn−ujn−1+(ΔxcΔt)2(uj+1n−2ujn+uj−1n)
CFL 条件：Courant 数 C=cΔt/Δx≤1。当 C>1，数值域不包含物理域，解必然不稳定。
数值实验：拉动长度为 L=1 的弦中点释放后的传播，c=1：
| Courant 数 C |
结果 |
分析 |
| 0.5 |
波形正确传播，小幅数值耗散 |
可用，安全 |
| 0.9 |
波形保持良好，传播正确 |
接近极限 |
| 1.0 |
波形完美传播 |
理论上精确（对于线性波动方程） |
| 1.2 |
出现高频振荡并迅速发散 |
不稳定 |
静电势满足泊松方程 ∇2ϕ=−ρ/ε0。二维离散化：
ϕi+1,j+ϕi−1,j+ϕi,j+1+ϕi,j−1−4ϕi,j=−h2ρi,j/ε0
迭代求解方法收敛速度对比（网格 128×128，误差容限 10−6）：
| 方法 |
迭代次数 |
每步计算量 |
相对耗时 |
适用场景 |
| Jacobi |
约 40,000 |
O(N) |
100x |
教学示例 |
| Gauss-Seidel |
约 20,000 |
O(N) |
50x |
小型问题 |
| SOR (ω=1.9) |
约 600 |
O(N) |
1.5x |
中等规模 |
| 共轭梯度法 (CG) |
约 200 |
O(N) |
1x |
对称正定 |
| 多重网格 V-cycle |
约 10 |
O(N) |
0.3x |
大规模问题 |
多重网格法的核心思想是：在粗网格上快速消除低频误差，在细网格上消除高频误差——不同频率的误差分量在不同网格层上被高效消除。
蒙特卡洛（Monte Carlo）方法利用随机数来求解确定性问题。名称源自摩纳哥的蒙特卡洛赌场，因为随机数与赌博轮盘相似。该方法的开创者之一 Stanislaw Ulam 在 1946 年因病卧床休息时，思考用随机抽样计算核裂变链式反应的概率——这直接催生了现代 MC 方法。
计算定积分 I=∫abf(x)dx 的最简单 MC 方法：
I≈Nb−ai=1∑Nf(xi)
其中 xi 是 [a,b] 上的均匀随机数。
误差特性：MC 积分的标准差为 σf/N，其中 σf 是被积函数的标准差。这意味着：
- 误差收敛速度为 O(1/N)——慢但稳定
- 误差与问题维度 d 无关——这是 MC 相对传统数值积分（误差 O(N−k/d)）的巨大优势
- 要提高 10 倍精度，需要 100 倍的样本量
例子：计算 ∫01sin(πx)dx=2/π≈0.63662：
| N |
MC 估计 |
标准误差 |
相对误差 |
| 100 |
0.6317 |
0.0507 |
0.77% |
| 1,000 |
0.6389 |
0.0160 |
0.36% |
| 10,000 |
0.6358 |
0.0051 |
0.13% |
| 100,000 |
0.6364 |
0.0016 |
0.03% |
| 1,000,000 |
0.6366 |
0.0005 |
0.01% |
这是马尔可夫链蒙特卡洛（MCMC）的核心算法，从复杂分布 π(x) 中生成样本：
- 从当前状态 xt 出发，根据提议分布 q(x′∣xt) 生成候选 x′
- 计算接受概率：α=min(1,π(xt)q(x′∣xt)π(x′)q(xt∣x′))
- 以概率 α 接受 x′（即 xt+1=x′），否则 xt+1=xt
示例：从一维高斯分布 N(0,1) 抽样，提议分布为 N(xt,0.52)：
| 初始值 |
预热步数 |
抽样步数 |
样本均值 |
样本标准差 |
理论均值 |
理论标准差 |
| 0.0 |
100 |
10,000 |
-0.003 |
1.002 |
0 |
1 |
| 10.0 |
100 |
10,000 |
0.015 |
0.995 |
0 |
1 |
| 0.0 |
500 |
100,000 |
0.001 |
0.999 |
0 |
1 |
注意：即使初始值偏离很远（10 个标准差），经过预热后的抽样仍然收敛到正确分布。
二维 Ising 模型 Hamiltonian：
H=−J⟨i,j⟩∑sisj,si∈{−1,+1}
每步尝试翻转一个自旋，能量变化：
ΔE=2Jsij∈neighbors∑sj
接受概率：min(1,e−βΔE)，其中 β=1/(kBT)
模拟结果（64×64 格子，J=1，kB=1）：
| 温度 T |
磁化强度 ⟨M⟩ |
能量 ⟨E⟩ |
相态 |
关联长度 |
| 1.0 |
0.998 |
-1.997 |
铁磁有序 |
远大于系统尺寸 |
| 2.0 |
0.971 |
-1.941 |
铁磁有序 |
与系统尺寸可比 |
| 2.269（Tc） |
~0.002 |
-1.534 |
临界点 |
发散 |
| 3.0 |
~0.001 |
-1.251 |
顺磁无序 |
约 1-2 个格点 |
| 5.0 |
~0.000 |
-0.998 |
顺磁无序 |
约 1 个格点 |
在 Tc≈2.269 处发生二级相变（Onsager 精确解），磁化强度从有限值连续变为零。
分子动力学（Molecular Dynamics, MD）通过数值求解经典运动方程来模拟粒子系统的时间演化：
midt2d2ri=Fi=−∇iU(r1,r2,…,rN)
分层视角：MD 位于量子力学（最精确但最慢）和连续介质力学（最不精确但最快）之间。
最常用的 MD 积分器是 Verlet 算法及其变体。
基本 Verlet 算法：
ri(t+Δt)=2ri(t)−ri(t−Δt)+ai(t)Δt2+O(Δt4)
速度 Verlet 算法（更实用，同时给出位置和速度）：
ri(t+Δt)vi(t+Δt)=ri(t)+vi(t)Δt+21ai(t)Δt2=vi(t)+21[ai(t)+ai(t+Δt)]Δt
Verlet 算法的优点：
- 辛结构保持：在相空间中保持体积守恒（长时间稳定性）
- 时间可逆：准确的时间反演性质
- 局部能量守恒：长时间模拟中能量漂移极小
LJ 势描述简单原子（稀有气体、液态金属等）之间的相互作用：
ULJ(r)=4ε[(rσ)12−(rσ)6]
其中 ε 是势阱深度（结合能强度），σ 是零势能距离（有效原子直径）。
LJ 势的氩原子 MD 模拟（864 原子，NPT 系综，p=1 atm）：
| 温度 |
密度 |
相态 |
径向分布函数第一峰 |
自扩散系数 (10−9 m²/s) |
均方位移斜率 |
| 80 K |
1.72 g/cm³ |
固态 |
3.75 Å（尖锐） |
0.002 |
近零（振动） |
| 90 K |
1.65 g/cm³ |
过冷液 |
3.80 Å |
0.85 |
缓慢增长 |
| 120 K |
1.38 g/cm³ |
液态 |
3.85 Å |
3.21 |
线性增长 |
| 150 K |
1.20 g/cm³ |
液态 |
3.95 Å |
5.67 |
快速增长 |
| 200 K |
0.95 g/cm³ |
气态 |
无明确峰 |
12.34 |
快速线性增长 |
径向分布函数 g(r)：描述距参考原子 r 处找到另一个原子的概率密度，是分析液体结构的核心工具。
| 系综 |
守恒量 |
物理对应 |
应用场景 |
| NVE（微正则） |
E, V, N |
孤立系统 |
算法验证 |
| NVT（正则） |
T, V, N |
恒温系统 |
溶液、生物模拟 |
| NPT（等温等压） |
T, p, N |
开放系统 |
最接近实验条件 |
| NPH（等焓等压） |
H, p, N |
绝热等压 |
相变研究 |
温度控制方法对比：
| 方法 |
原理 |
优点 |
缺点 |
| Velocity rescaling |
每 n 步调整速度 |
简单 |
不产生正则系综 |
| Andersen thermostat |
碰撞随机重置速度 |
严格正则 |
动力学信息丢失 |
| Nose-Hoover |
扩展 Lagrangian |
决定论性、严格 |
参数需调优 |
| Langevin dynamics |
Langevin 方程 |
稳定、自然 |
摩擦系数参数化 |
一维薛定谔方程：
−2mℏ2dx2d2ψ+V(x)ψ=Eψ
有限差分离散化：
−2mℏ2Δx2ψi+1−2ψi+ψi−1+Viψi=Eψi
这写为矩阵特征值问题 Hψ=Eψ，其中 H 是一个三对角矩阵。
例子：无限深方势阱（V=0 在 0<x<L）
当 m=1, ℏ=1, L=1, N=50 网格：
| 能级 n |
解析解 En=n2π2/2 |
数值解 |
相对误差 |
| 1 |
4.935 |
4.935 |
1.3×10−4 |
| 2 |
19.739 |
19.737 |
6.1×10−5 |
| 3 |
44.413 |
44.402 |
2.5×10−4 |
| 4 |
78.957 |
78.912 |
5.7×10−4 |
| 5 |
123.370 |
123.170 |
1.6×10−3 |
高阶本征态误差更大，因为其波函数振荡更快，需要更精细的网格。
DFT 将 N 电子问题转化为单电子有效势问题。Kohn-Sham 方程：
[−2mℏ2∇2+Vext(r)+VH(r)+Vxc(r)]ψi(r)=εiψi(r)
自洽求解流程：
初始电荷密度 ρ_in
↓
构建有效势 V_eff = V_ext + V_H[ρ] + V_xc[ρ]
↓
求解 Kohn-Sham 方程 → 新波函数 ψ_i
↓
计算新电荷密度 ρ_out = Σ|ψ_i|^2
↓
混合：ρ_new = α ρ_in + (1-α) ρ_out
↓
判断收敛：|ρ_new - ρ_in| < tol ?
├── 否 → 回到第二步
└── 是 → 输出总能、力、能带等
不同泛函的精度对比：
| 泛函类型 |
计算复杂度 |
晶格常数误差 |
带隙误差 |
代表 |
| LDA |
O(N3) |
1-2% 低估 |
40-100% 低估 |
PZ81 |
| GGA |
O(N3) |
0.5-1% |
30-50% 低估 |
PBE |
| 杂化泛函 |
O(N4) |
0.2-0.5% |
5-15% |
HSE06 |
| Meta-GGA |
O(N3) |
0.3-0.8% |
20-30% |
SCAN |
| GW |
O(N4) |
0.1-0.3% |
0.1-0.3 eV |
G0W0 |
Navier-Stokes 方程：
ρ∂t∂v+ρ(v⋅∇)v=−∇p+μ∇2v+f
∇⋅v=0（不可压缩）
数值方法对比：
| 方法 |
优点 |
缺点 |
典型软件 |
| 有限体积法（FVM） |
守恒性好，工业标准 |
精度受网格限制 |
OpenFOAM, Fluent |
| 有限差分法（FDM） |
实现简单 |
复杂几何困难 |
自编代码 |
| 格子玻尔兹曼（LBM） |
并行性好，简单 |
高马赫数不准确 |
Palabos |
| 谱方法 |
超高精度 |
仅限规则区域 |
Dedalus |
| SPH（无网格） |
自由表面处理 |
边界条件复杂 |
LAMMPS |
硅晶体的 DFT 计算（实验 vs LDA vs GGA-PBE）：
| 物理量 |
实验值 |
LDA |
GGA-PBE |
| 晶格常数 (Å) |
5.431 |
5.397（-0.6%） |
5.469（+0.7%） |
| 体模量 (GPa) |
97.8 |
100.2（+2.5%） |
88.5（-9.5%） |
| 弹性常数 C11 (GPa) |
165.8 |
169.3（+2%） |
157.2（-5%） |
| 带隙 (eV) |
1.12 |
0.52（-54%） |
0.67（-40%） |
| 内聚能 (eV/atom) |
4.63 |
5.29（+14%） |
4.66（+0.6%） |
LDA 倾向"过度键合"（晶格常数偏小、结合能偏大），GGA 倾向"欠键合"（晶格常数偏大、体模量偏小），两者都系统性低估带隙。
星系演化模拟的核心：
dt2d2ri=−Gj=i∑∣ri−rj∣3mj(ri−rj)
直接求和 O(N2) 在 N>105 时不可行。
| 算法 |
计算复杂度 |
典型 N |
精度 |
内存 |
| 直接 N 体 |
O(N2) |
104 |
精确 |
O(N) |
| Barnes-Hut 树码 |
O(NlogN) |
106 |
高 |
O(N) |
| 粒子网格 (PM) |
O(NlogN) |
108 |
中 |
O(N) |
| 快速多极子 (FMM) |
O(N) |
109 |
高 |
O(N) |
| 语言 |
性能 |
开发效率 |
并行支持 |
主要用途 |
| Fortran |
⭐⭐⭐⭐⭐ |
⭐⭐ |
MPI+OpenMP |
传统超算代码 |
| C/C++ |
⭐⭐⭐⭐⭐ |
⭐⭐⭐ |
MPI+CUDA+OpenMP |
主流高性能计算 |
| Python |
⭐⭐⭐ |
⭐⭐⭐⭐⭐ |
Numba/MPI |
原型+后处理 |
| Julia |
⭐⭐⭐⭐ |
⭐⭐⭐⭐ |
内置并行 |
新兴选择 |
| CUDA C |
⭐⭐⭐⭐⭐+ |
⭐⭐ |
GPU原生 |
GPU加速 |
| 领域 |
软件 |
方法 |
开源 |
特色 |
| 分子动力学 |
LAMMPS |
MD |
✅ |
高性能、可扩展 |
| 分子动力学 |
GROMACS |
MD |
✅ |
生物分子模拟首選 |
| 第一性原理 |
VASP |
DFT/PAW |
❌ |
材料科学标准 |
| 第一性原理 |
Quantum ESPRESSO |
DFT/PW |
✅ |
完全开源 |
| 量子化学 |
Gaussian |
波函数 |
❌ |
有机化学标准 |
| CFD |
OpenFOAM |
FVM |
✅ |
工业标准开源 |
| 有限元 |
deal.II |
FEM |
✅ |
自适应网格 |
| 张量网络 |
ITensor |
DMRG/MPS |
✅ |
一维量子系统 |
| 科学计算 |
NumPy/SciPy |
通用 |
✅ |
Python 基础 |
import numpy as np
L, nx = 1.0, 50
dx = L / nx
alpha = 0.01
dt = 0.5 * dx**2 / alpha # CFL 条件
nt = 200
x = np.linspace(0, L, nx+1)
u = np.sin(np.pi * x)
r = alpha * dt / dx**2
for n in range(nt):
un = u.copy()
u[1:-1] = (1 - 2*r) * un[1:-1] + r * (un[2:] + un[:-2])
print(f"t = {nt*dt:.3f}s, u_max = {u.max():.4f}")
# 解析解: exp(-απ²t)·sin(πx)
u_exact = np.exp(-alpha * np.pi**2 * nt*dt) * np.sin(np.pi * x)
print(f"Error = {np.abs(u - u_exact).max():.2e}")
输出：Error = 6.23e-04，验证了算法的正确性。
PINN 将 PDE 残差嵌入神经网络损失函数：
L=数据拟合Nd1∑∥u^−udata∥2+λPDE 残差Nf1∑∂t∂u^−N[u^]2
PINN 与传统方法的对比（Burgers 方程）：
| 指标 |
FDM (1282 网格) |
PINN (4层×256神经元) |
| 正向求解时间 |
0.1 秒 |
5 分钟（训练） |
| 泛化到新参数 |
需完整重算 |
需重新训练 |
| 逆向问题支持 |
需额外推导 |
自然支持 |
| 精度 (L2 误差) |
1.2×10−4 |
5.6×10−4 |
| 无网格性 |
否 |
是 |
| 复杂边界条件 |
成熟 |
较困难 |
用神经网络替代传统力场，保持接近 DFT 精度的同时大幅加速：
性能阶梯：
| 方法 |
能量精度 |
速度（步/秒·千原子） |
适用场景 |
| DFT (VASP) |
< 1 meV/atom |
10−3 |
小体系基准 |
| 经典力场 |
10-50 meV/atom |
103 |
大体系粗粒化 |
| NequIP (等变NN) |
1-5 meV/atom |
101 |
高精度大体系 |
| DeepMD |
2-8 meV/atom |
102 |
中尺度模拟 |
变分量子本征求解器（VQE）：
E(θ)=⟨ψ(θ)∣H∣ψ(θ)⟩
通过经典-量子混合算法最小化能量，有望在近期量子设备上解决经典计算机难以处理的量子化学问题。
计算物理已经从辅助工具成长为物理学的核心支柱。从有限差分法到物理信息神经网络，从数十个原子的 MD 模拟到百万原子级别的机器学习力场计算，计算能力和方法的进步不断拓展我们探索物理世界的边界。
关键要点：
- 离散化是计算物理的根本——将连续的物理问题转化为计算机可处理的数值形式
- 算法选择直接影响效率和精度——理解每种方法的适用场景和局限
- 验证确认至关重要——数值结果必须与解析解或实验数据进行对比验证
- 并行计算是现代计算物理的基石——GPU 和分布式系统使大规模模拟成为可能
- 数据驱动的趋势不可逆——机器学习和物理模拟的深度融合正在改变研究范式
延伸阅读：数值分析 | 偏微分方程 | 线性代数 | 概率论 | 机器学习数学