3. Cahn-Hilliard方程半隐式傅里叶谱方法
半隐式谱方法的离散思路
Cahn-Hilliard方程 ∂c/∂t = M∇²(∂f/∂c − κ∇²c) 中，双阱势导数 f′(c)=c³−c 的线性部分 −c 与四阶项 −κ∇⁴c 在傅里叶空间分别对应增长率 −Mk² 和 −Mκk⁴。当 k 较大时，−Mκk⁴ 的量级迅速增大，显式格式的稳定步长被压缩到 Δt ≤ 2/(Mκk⁴_max)，网格越密越难推进。半隐式谱方法把线性刚性项 −Mκk⁴ 隐式处理，非线性项 c³ 显式处理，从而在保持谱精度离散的同时放宽时间步限制。
在周期边界下，序参量 c(x,t) 展开为傅里叶级数，谱空间离散后每个波数 k 独立演化。对线性部分采用隐式、非线性部分采用显式，得到逐波数更新的半隐式格式：
ĉ_k^{n+1} = (ĉ_k^n − Δt·M·k²·FFT[(c^n)³]_k) / (1 + Δt·M·κ·k⁴)其中 FFT 与 IFFT 用于在实空间计算非线性项 c³，再变换回谱空间。分母 1+Δt·M·κ·k⁴ 对高波数提供强阻尼，使格式对刚性项无条件稳定（非线性项仍受显式稳定性约束，但限制远小于纯显式）。
关键实现与参数选择
实现时需注意以下几点：
- 波数构造：对 N 个格点、周期长度 L，波数 k = 2π·[0,1,…,N/2−1,−N/2,…,−1]/L，保证 FFT 频率顺序正确。
- 非线性项处理：先在实空间计算 c³，再 FFT 回谱空间，避免谱空间卷积的高开销。
- 时间步长：半隐式格式中 Δt 不再受 κk⁴ 限制，但仍受非线性项显式稳定性约束，通常取 Δt ~ 10⁻³ 量级即可稳定推进。
- 守恒性：k=0 模对应总浓度，更新公式中 k=0 时 k²=0，ĉ_0 保持不变，自动满足质量守恒。
显式与半隐式效率对比
| 对比项 | 显式有限差分 | 半隐式傅里叶谱 |
|---|---|---|
| 空间离散 | 二阶中心差分 | 谱方法（指数收敛） |
| 时间步长限制 | Δt ≤ 2/(Mκk⁴_max) | 仅受非线性项约束 |
| N=128 时典型 Δt | ~10⁻⁷ | ~10⁻³ |
| 推进 t=0.5 步数 | ~4×10⁶ | ~500 |
| 每步计算量 | O(N) | O(N log N) |
| 总计算量（估算） | ~4×10⁶·N | ~500·N log N |
半隐式谱方法的核心优势在于：用谱方法获得空间高精度，用隐式处理线性刚性项解除时间步限制，两者结合使长时间旋节线分解模拟变得可行。实际使用时建议先用小 N 验证格式正确性，再逐步增大 N 和 Δt。
半隐式格式并非无条件稳定：非线性项 c³ 仍为显式处理，当 Δt 过大时可能因非线性项导致数值振荡。建议通过能量单调性或总浓度守恒性监测稳定性，必要时减小 Δt 或对非线性项做部分隐式处理。