显式 vs 隐式
显式方法简单但有稳定性限制；隐式方法可取更大时间步但需求解线性系统。
CS267 · 第 6 课 · James Demmel
偏微分方程 (PDEs) 是科学计算的核心。本课探讨 PDEs 的显式与隐式求解方法、并行化策略，以及不同算法的复杂度与适用场景。
上一课介绍了四类模拟：离散事件系统、粒子系统、ODEs，以及本课的 PDEs。
共性问题：
新需求：可重复性 (Reproducibility)
| 类型 | 依赖关系 | 示例 |
|---|---|---|
| 椭圆型 (Elliptic) | 稳态，全局空间依赖 | 静电/引力势：\(V(\text{position})\) |
| 双曲型 (Hyperbolic) | 时间依赖，局部空间依赖 | 声波：\(P(\text{position}, \text{time})\) |
| 抛物型 (Parabolic) | 时间依赖，全局空间依赖 | 热流：\(T(\text{position}, \text{time})\) |
考虑一根均匀材料的杆，两端绝缘。设 \(u(x,t)\) 为位置 \(x\) 时间 \(t\) 的温度：
\[\frac{\partial u(x,t)}{\partial t} = C \cdot \frac{\partial^2 u(x,t)}{\partial x^2}\]
离散化时间和空间：
// 用有限差分近似
(u(x,t+dt) - u(x,t))/dt = C * [u(x-h,t) - 2u(x,t) + u(x+h,t)] / h²
// 解出 u(x,t+dt)
u(x,t+dt) = u(x,t) + z * [u(x-h,t) - 2u(x,t) + u(x+h,t)]
// 其中 z = C*dt/h²简化（设 \(x = j \cdot h\), \(t = i \cdot dt\)）：
\[u[j, i+1] = z \cdot u[j-1, i] + (1-2z) \cdot u[j, i] + z \cdot u[j+1, i]\]
\[u[:, i+1] = T \cdot u[:, i]\]
其中 T 是三对角矩阵：
\[T = I - zL, \quad L = \begin{pmatrix} 2 & -1 \\ -1 & 2 & -1 \\ & \ddots & \ddots & \ddots \\ & & -1 & 2 \end{pmatrix}\]
L 称为 Laplacian（1D）。对于 2D 网格（5-point stencil），L 是五对角的。
数值不稳定性：对于热方程（抛物型），显式方法要求 \(z = C \cdot dt / h^2 \leq 0.5\)。
当 h 很小时，\(dt < 0.5 \cdot h^2 / C\)，需要极小的时间步！
// Backward Euler
(u(x,t+dt) - u(x,t))/dt = C * [u(x-h,t+dt) - 2u(x,t+dt) + u(x+h,t+dt)] / h²
// 整理得
(I + zL) * u[:, i+1] = u[:, i]
// 需要求解稀疏线性系统更好的数值性质：
\[(I + (z/2) \cdot L) \cdot u[:, i+1] = (I - (z/2) \cdot L) \cdot u[:, i]\]
椭圆型 PDE 的特例（稳态）：
\[\frac{\partial^2 V}{\partial x^2} + \frac{\partial^2 V}{\partial y^2} + \frac{\partial^2 V}{\partial z^2} = 0\]
与引力、静电势密切相关：
对于 2D 网格（5-point stencil）：
\[L = \begin{pmatrix} 4 & -1 & -1 \\ -1 & 4 & -1 & -1 \\ & \ddots & \ddots & \ddots & \ddots \\ & & -1 & -1 & 4 \end{pmatrix}\]
| 算法 | 串行复杂度 | 并行复杂度 |
|---|---|---|
| Dense LU | \(N^3\) | \(N\) |
| Band LU | \(N^2\) | \(N\) |
| Jacobi | \(N^2\) | \(N\) |
| Explicit Inverse | \(N\) | \(\log N\) |
| Conjugate Gradient | \(N^{3/2}\) | \(N^{1/2} \log N\) |
| Red/Black SOR | \(N^{3/2}\) | \(N^{1/2}\) |
| Sparse LU | \(N^{3/2}\) | \(N^{1/2}\) |
| FFT | \(N \log N\) | \(\log N\) |
| Multigrid | \(N\) | \(\log^2 N\) |
| Lower Bound | \(N\) | \(\log N\) |
Multigrid 是最快的！串行复杂度 \(O(N)\)，接近理论下界。
核心思想：在粗网格上消除低频误差，在细网格上消除高频误差。
1024-CPU NCUBE-2 上的对流扩散问题：
| Solver | CPU Time (s) | Mflop/s |
|---|---|---|
| Jacobi | 2124 | 1800 |
| Gauss-Seidel | 885 | 1365 |
| Multigrid | 7 | 318 |
结论：算法选择比峰值性能更重要！
| 类型 | 描述 |
|---|---|
| Regular | 1D/2D/3D 规则网格，基础构建块 |
| Composite | 多个"弯曲"规则网格在边缘连接 |
| Unstructured | 任意网格点和连接性 |
| Adaptive (AMR) | 求解过程中动态调整分辨率 |
在求解过程中动态调整网格分辨率：
不同目标需要不同的矩阵排序：
| 目标 | 排序策略 |
|---|---|
| 并行ism | 最小化非对角非零元（图分区） |
| 减少 fill-in | 最小度排序 (Minimum Degree) |
| 改善数值稳定性 | 列主元 (Column Pivoting) |
Q1：热方程的显式方法为什么需要极小的时间步？
Q2：Multigrid 算法的串行复杂度是多少？
Q3：隐式方法相比显式方法的主要优点是什么？
Q4：Poisson 方程属于哪种类型的 PDE？
Q5：非结构网格并行求解的主要挑战不包括什么？
显式方法简单但有稳定性限制；隐式方法可取更大时间步但需求解线性系统。
对于椭圆型 PDE，Multigrid 达到 \(O(N)\) 串行复杂度，接近理论下界。
实际中 Multigrid 比 Jacobi 快 300 倍，尽管峰值 Mflop/s 更低。
需要专门的工具处理网格生成、分区、求解器设计。