numerical-stability

数值稳定性

目标

提供一份可重复使用的检查清单和由脚本驱动的检查，以确保时变模拟稳定且有据可依。

要求

• Python 3.10+
• NumPy（用于 matrix_condition.py 和 von_neumann_analyzer.py）
• 依赖项请参阅 scripts/requirements.txt

需收集的输入

输入	说明	示例
网格间距 dx	空间离散	0.01 m
时间步长 dt	时间离散	1e-4 s
速度 v	对流速度	1.0 m/s
扩散率 D	热扩散率／质量扩散率	1e-5 m²/s
反应速率 k	一阶速率常数	100 s⁻¹
维数	1D、2D 或 3D	2
格式类型	显式或隐式	explicit

决策指南

选择显式还是隐式

Is the problem stiff (fast + slow dynamics)?
├── YES → Use implicit or IMEX scheme
│         └── Check conditioning with matrix_condition.py
└── NO → Is CFL/Fourier satisfied with reasonable dt?
    ├── YES → Use explicit scheme (cheaper per step)
    └── NO → Consider implicit or reduce dx

稳定性限制速查表

物理过程	数值	显式限制（1D）	公式
对流	CFL	C ≤ 1	C = v·dt/dx
扩散	Fourier	Fo ≤ 0.5	Fo = D·dt/dx²
反应	反应数	R ≤ 1	R = k·dt

多维修正：对于 d 个维度，扩散限制为 Fo ≤ 1/(2d)。

脚本输出（JSON 字段）

脚本	主要输出
scripts/cfl_checker.py	metrics.cfl、metrics.fourier、recommended_dt、stable
scripts/von_neumann_analyzer.py	results.max_amplification、results.stable
scripts/matrix_condition.py	results.condition_number、results.is_symmetric
scripts/stiffness_detector.py	results.stiffness_ratio、results.real_part_stiffness_ratio、results.imag_dominated、results.stiff、results.recommendation、results.warning

工作流

1. 识别主导物理过程（对流与扩散或反应）
2. 使用 scripts/cfl_checker.py 运行 CFL 检查器
3. 与限制值比较，并在需要时调整 dt
4. 如果存在多个尺度，使用 scripts/stiffness_detector.py 检查刚性
5. 使用 scripts/von_neumann_analyzer.py 分析自定义格式
6. 对隐式求解使用 scripts/matrix_condition.py 检查条件数
7. 记录稳定性判断和建议的时间步长

对话工作流示例

用户：我的相场模拟在运行 100 步后发散了。我使用的是显式 Euler，dx=0.01、dt=1e-4，扩散率 D=1.0。

Agent 工作流：

1. 检查稳定性准则：

   python3 scripts/cfl_checker.py --dx 0.01 --dt 1e-4 --diffusivity 1.0 --dimensions 2 --json
   

2. 解读结果：
   • Fourier 数：Fo = 1.0 × 1e-4 / (0.01)² = 1.0
   • 2D 限制：Fo ≤ 0.25
   • 超出限制：Fo = 1.0 > 0.25，不稳定！
3. 建议修复方法：
   • 将 dt 减小到 2.5e-5（工具给出的 recommended_dt，此时 Fo = 0.25）
   • 或增大 dx，或改用隐式格式

模拟前稳定性检查清单

• 识别主导物理过程和无量纲群
• 使用 cfl_checker.py 计算 CFL 数／Fourier 数／反应数
• 如果使用显式格式且超出限制，则减小 dt 或更换格式
• 如果刚性比 > 1000，则选择隐式求解器／刚性积分器
• 对于自定义格式，验证放大因子 ≤ 1
• 记录稳定性判断依据及输入和输出

CLI 示例

# Check CFL/Fourier for 2D diffusion-advection
python3 scripts/cfl_checker.py --dx 0.1 --dt 0.01 --velocity 1.0 --diffusivity 0.1 --dimensions 2 --json

# Von Neumann analysis for custom 3-point stencil
python3 scripts/von_neumann_analyzer.py --coeffs 0.2,0.6,0.2 --dx 1.0 --nk 128 --json

# Detect stiffness from eigenvalue estimates
python3 scripts/stiffness_detector.py --eigs=-1,-1000 --json

# Check matrix conditioning for implicit system
python3 scripts/matrix_condition.py --matrix A.npy --norm 2 --json

错误处理

错误	原因	解决方法
dx and dt must be positive	值为零或负数	提供有效的正数
No stability criteria applied	缺少速度／扩散率	至少提供一个物理参数
Matrix not found: <path>	路径无效	检查矩阵文件是否存在
Could not compute eigenvalues	矩阵奇异或格式不正确	检查矩阵是否有效

解读指南

场景	含义	操作
stable: true	所有已检查的标准均满足	继续进行模拟
stable: false	至少有一项限制被违反	减小 dt 或更改方案
stable: null	无法应用任何标准	提供更多物理参数
刚度比 > 1000	问题是刚性的	使用隐式积分器
条件数 > 10⁸	条件较差（status: poorly-conditioned）	可能需要预条件
条件数 > 10¹⁰	病态（status: ill-conditioned）	使用缩放／预条件

条件阈值假定采用 IEEE 双精度：求解会损失大约 log10(κ) 位有效数字，且当 κ ≈ 1/eps ≈ 4.5e15 时，矩阵会变为数值奇异。> 10⁸／> 10¹⁰ 的阈值（与 matrix_condition.py 的 status 一致）留有充足余量；对于双精度而言条件良好的 FEM 矩阵（κ 最高约为 10⁶–10⁷）会报告 status: ok。

验证清单

在满足以下所有适用项，并从脚本中记录相应数值之前，不要将某个方案／时间步声明为“稳定”。

• 运行了 cfl_checker.py --json，并记录 metrics.cfl、metrics.fourier 和／或 metrics.reaction 与报告的 limits.* 的对比（显式默认值：CFL ≤ 1、Fo ≤ 1/(2d)、R ≤ 1），且 stable: true，预期采用的标准已出现在 criteria_applied 中。
• 确认 criteria_applied 非空且 stable 不为 null——也就是说，确实提供了至少一个物理参数（--velocity／--diffusivity／--reaction-rate），因此执行了实际检查，而非静默地什么也没做。
• --dx 使用所有方向中最小的网格间距（各向异性网格：最小的 dx／dy／dz），并在每次网格细化后重新运行 cfl_checker.py，因为 Fo ∝ dt/dx²，限值对 dx 高度敏感。
• 如果所选 dt 低于限值，记录工具给出的 recommended_dt（以及使用的 --safety 因子），以明确并确保稳定边界余量可复现。
• 对于任何自定义／非标准更新模板，运行了 von_neumann_analyzer.py --json，并确认 results.max_amplification ≤ 1（stable: true）；记录 k_at_max，并处理所有关于偶数长度模板的 warning。
• 对于多尺度系统，运行了 stiffness_detector.py --json，并记录 real_part_stiffness_ratio、imag_dominated 和 stiff；仅当实部刚度比显示 stiff: true 且没有 warning 时，才选择 BDF／Radau，而非仅依据幅值。
• 对于隐式求解，运行了 matrix_condition.py --json，并记录 condition_number 和 status；将 poorly-conditioned（>1e8）／ill-conditioned（>1e10）视为应先进行缩放／预条件，再信任求解结果的标志。

常见陷阱与错误理由

看似省事的做法	错在哪里／应该怎么做
“方案是隐式的，所以 dt 取什么值都行。”	无条件稳定不等于准确。对于 --scheme implicit，cfl_checker.py 会在放宽（无限制）的限值下报告 stable: true，甚至会添加“检查准确性”的说明——较大的 dt 仍会破坏时间误差。设置 dt 时应考虑准确性，而不仅仅是稳定性。
“运行了 100 步都没崩溃，所以设置是稳定的。”	舍入误差、守恒量丢失或临界 Fo 都可能导致后期爆炸。运行完成 ≠ 结果正确——记录 metrics.fourier/metrics.cfl 与 limits.* 的对比，并确认 stable: true 后再相信运行结果。
“一维傅里叶数限值是 0.5，所以 Fo ≤ 0.5 就安全。”	显式扩散的限值是 Fo ≤ 1/(2d)——二维为 0.25，三维为 0.167。传入实际的 --dimensions；cfl_checker.py 会自动收紧 diffusion_limit，而在二维／三维中使用 0.5 会导致不稳定。
“刚度比很大，所以用 BDF/Radau。”	对于振荡／平流／哈密顿谱，stiffness_ratio 的大小容易误导。检查 imag_dominated 和 real_part_stiffness_ratio：如果 imag_dominated: true，检测器会返回 stiff: false 并附带 warning——应优先选择辛格式／leapfrog，或步长按 CFL 限制的 A 稳定格式，而不是隐式刚性求解器。
“我缩小了 dt，所以细化网格后，旧网格的 dt 仍然适用。”	Fo ∝ dt/dx²：dx 减半会使 Fo 增至四倍。沿用网格细化前的 dt 会再次违反限值。每次更改网格后，都应使用 cfl_checker.py 重新计算。
“矩阵求解成功了，所以它的条件数没问题。”	求解器可能返回数值，却悄然损失约 log10(κ) 位有效数字。记录 matrix_condition.py 返回的 condition_number/status；如果结果是 poorly-conditioned/ill-conditioned，应先进行缩放或预处理，再相信结果。

安全性

输入验证

• 在任何计算开始前，都会将 dx、dt 和 safety 验证为有限正数；如果提供了 velocity、diffusivity 和 reaction_rate，则会将它们验证为有限数值
• --dimensions 仅限于 {1, 2, 3}（否则 cfl_checker.py 会报错并以状态码 2 退出）
• 逗号分隔的特征值列表（--eigs）最多允许 10,000 项，并会验证每项均为有限数值
• 模板系数列表（--coeffs）最多允许 10,000 项，并会验证每项均为有限浮点数

文件访问

• matrix_condition.py 只读取由 --matrix 指定的单个矩阵文件（.npy 或文本）；不会遍历给定路径之外的目录
• 如果矩阵／雅可比文件超过 500 MB，则会在解析前拒绝处理（matrix_condition.py、stiffness_detector.py）
• 调用 np.load() 时会设置 allow_pickle=False，以防止通过特制的 .npy 文件执行任意代码（matrix_condition.py、stiffness_detector.py）
• 脚本只向 stdout 写入内容（JSON 输出）；除非 Agent 明确使用 Write 工具，否则不会创建文件

工具限制

• Read：用于检查脚本源代码、参考资料和用户配置文件
• Bash：用于通过显式参数列表执行四个 Python 分析脚本（cfl_checker.py、von_neumann_analyzer.py、matrix_condition.py、stiffness_detector.py）
• Write：用于保存分析结果或生成的报告；写入范围仅限于用户的工作目录
• Grep/Glob：用于定位相关文件和搜索参考资料

安全措施

• 不使用 eval()、exec() 或动态代码生成
• 所有子进程调用都使用显式参数列表（不使用 shell=True）
• matrix_condition.py 中的矩阵维度上限（每边 100,000）可防止内存耗尽
• JSON 输出模式（--json）会生成结构化、可解析的结果，不包含可由 shell 解释的内容

限制

• CFL/Fourier 检查仅适用于显式格式（隐式格式无条件稳定）
• Von Neumann 分析假设线性、常系数、周期性边界条件
• 刚性检测需要用户提供特征值估计。stiff 判定依据是真正衰减的模态（Re(λ) < 0）之间的尺度分离；仍会报告幅值比 stiffness_ratio。由虚轴主导的谱（振荡／对流／Hamilton 系统）会通过 imag_dominated 和 warning 标记，且不会被归类为刚性问题——对于这类问题，优先选择辛格式／leapfrog，或依据 CFL 限制确定步长的 A 稳定格式，而非 BDF/Radau。

参考资料

• references/stability_criteria.md - 判定阈值与公式
• references/common_pitfalls.md - 常见失效模式与修正方法
• references/scheme_catalog.md - 常见格式的稳定性属性

版本历史

• v1.2.2（2026-06-24）：新增“验证清单”（7 项基于证据、与 4 个脚本的 JSON 输出对应的检查项）和“常见陷阱与自我合理化”表格（6 行），让稳定性判定可复现，并防止因走捷径而产生虚假的信心
• v1.2.0（2026-06-23）：修正相场示例（D=1.0 时 Fo=1.0 确实不稳定）；修正 eval case 1 的 Fourier 值；统一条件数阈值和错误消息文档与脚本之间的内容；基于实部的刚性分类，并检测虚轴主导情况；落实文档规定的输入、文件大小和维度防护措施
• v1.1.0（2024-12-24）：完善文档、判定指南和示例
• v1.0.0：首次发布，包含 4 个稳定性分析脚本