1. 一维Allen-Cahn方程显式有限差分解

问题设定与离散格式

Allen-Cahn 方程描述非保守序参量 η(x,t) 的弛豫过程：

∂η/∂t = -M (δF/δη) = M (ε² ∇²η - f'(η))
f(η) = (η² - 1)² / 4  →  f'(η) = η³ - η

其中 M 为迁移率，ε 为界面宽度参数。双阱势 f(η) 在 η = ±1 处取极小值，因此平衡态是 η = ±1 的均匀相，两相之间的过渡区由 ε 控制宽度。本节取 M = 1、ε = 0.03，计算域 x ∈ [0, 1]，周期边界条件。

显式有限差分离散：空间用二阶中心差分，时间用一阶显式 Euler：

η_i^{n+1} = η_i^n + Δt · [ (η_{i+1}^n - 2η_i^n + η_{i-1}^n) / Δx² · ε² - ( (η_i^n)³ - η_i^n ) ]

显式格式的稳定性由扩散项决定，需满足 Δt ≤ Δx² / (2ε²)。取 N = 200 个网格点，Δx = 1/200 = 0.005，则 Δt ≤ 0.005² / (2×0.03²) ≈ 0.0139。为留出安全裕度，取 Δt = 0.005，总步数 20000 步，对应物理时间 t = 100。

初始条件与演化机制

初值取围绕 η = 0 的小幅随机扰动，使体系处于双阱势的势垒顶部附近。由于 f'(η) 在 η = 0 附近为负（η³ - η < 0），扰动会被放大，η 逐渐向 ±1 靠拢；同时扩散项 ε²∇²η 会平滑界面，使相邻区域趋于一致。最终形成若干 η = +1 与 η = -1 的畴区，畴壁在曲率驱动下收缩，小畴消失，大畴长大。

完整可运行代码

以下 Python 脚本实现上述离散格式，输出三个时刻的 η(x) 曲线，并保存为 PNG 图像用于观察演化过程。

import numpy as np
import matplotlib.pyplot as plt

# ---------- 参数 ----------
N = 200                # 网格点数
L = 1.0                # 计算域长度
dx = L / N             # 空间步长
M = 1.0                # 迁移率
eps = 0.03             # 界面宽度参数
dt = 0.005             # 时间步长
nsteps = 20000         # 总步数
save_every = 2000      # 每隔多少步记录一次

# ---------- 初始化 ----------
x = np.linspace(0, L, N, endpoint=False)
np.random.seed(42)
eta = 0.05 * (np.random.rand(N) - 0.5)   # 围绕 0 的小幅随机扰动

# 周期边界：用 np.roll 实现
def laplacian(eta):
    return (np.roll(eta, -1) - 2*eta + np.roll(eta, 1)) / dx**2

# ---------- 时间推进 ----------
history = [(0, eta.copy())]
for n in range(1, nsteps + 1):
    lap = laplacian(eta)
    detadt = M * (eps**2 * lap - (eta**3 - eta))
    eta = eta + dt * detadt
    if n % save_every == 0:
        history.append((n * dt, eta.copy()))

# ---------- 绘图 ----------
fig, axes = plt.subplots(1, 3, figsize=(12, 3.5), sharey=True)
for ax, (t, e) in zip(axes, [history[0], history[len(history)//2], history[-1]]):
    ax.plot(x, e, 'b-', linewidth=1.5)
    ax.axhline(1, color='g', linestyle='--', linewidth=0.8, alpha=0.6)
    ax.axhline(-1, color='g', linestyle='--', linewidth=0.8, alpha=0.6)
    ax.set_title(f't = {t:.1f}')
    ax.set_xlabel('x')
    ax.set_ylim(-1.3, 1.3)
    ax.grid(True, alpha=0.3)
axes[0].set_ylabel('η')
plt.tight_layout()
plt.savefig('allen_cahn_1d.png', dpi=150)
plt.show()

结果判读与验证要点

• 初值阶段（t ≈ 0）：η 在 0 附近小幅波动，尚未形成明显畴区。
• 中间阶段（t ≈ 20）：η 快速向 ±1 靠拢，出现多个正负畴区，界面处 η 在 -1 与 +1 之间过渡，过渡宽度约为 ε 量级。
• 后期阶段（t ≈ 100）：小畴被大畴吞并，畴壁数量减少，最终趋于少数大畴或单一均匀态。这是非保守序参量曲率驱动收缩的典型特征。
• 稳定性验证：若将 dt 增大到 0.02（超过 0.0139），数值解会出现振荡发散，可据此验证稳定性条件。
• 能量单调下降：可额外计算总自由能 F = Σ [ε²/2 (∇η)² + f(η)] Δx，应随演化单调递减，作为正确性的辅助判据。