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
拡散	フーリエ数	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ステップ後に発散します。dx=0.01、dt=1e-4、拡散係数 D=1.0 の陽的オイラー法を使用しています。

エージェントのワークフロー：

1. 安定性の基準を確認する：

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

2. 結果を解釈する：
   • フーリエ数：Fo = 1.0 × 1e-4 / (0.01)² = 1.0
   • 2D の限界：Fo ≤ 0.25
   • 違反：Fo = 1.0 > 0.25 のため、不安定です！
3. 修正を提案する：
   • dt を 2.5e-5（Fo = 0.25 となる、ツールの recommended_dt）に減らす
   • または dx を大きくするか、陰解法に切り替える

シミュレーション前の安定性チェックリスト

• 支配的な物理現象と無次元数を特定する
• cfl_checker.py で CFL 数／フーリエ数／反応数を計算する
• 陽解法で限界に違反した場合は、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	速度または拡散係数が未指定	物理パラメーターを少なくとも1つ指定してください
Matrix not found: <path>	パスが無効	行列ファイルが存在することを確認してください
Could not compute eigenvalues	行列が特異、または形式が不正	行列が有効であることを確認してください

解釈の指針

シナリオ	意味	対応
stable: true	すべての確認対象の基準を満たしています	シミュレーションを続行します
stable: false	1 つ以上の制限に違反しています	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 でないことを確認します。つまり、実際に検査が行われるよう、少なくとも 1 つの物理パラメーター（--velocity／--diffusivity／--reaction-rate）が指定されており、何も実行されない状態ではないことを確認します。
• --dx には全方向の格子間隔のうち最小の値を使用します（異方性格子では最小の dx／dy／dz）。Fo ∝ dt/dx² により上限が dx に非常に敏感なため、メッシュを細分化した後は cfl_checker.py を再実行します。
• 選択した 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 の場合にのみ BDF／Radau を選択し、warning がないことを確認します。
• 陰解法で解く場合は、matrix_condition.py --json を実行し、condition_number と status を記録します。poorly-conditioned（>1e8）／ill-conditioned（>1e10）は、解法の結果を信頼する前にスケーリング／前処理を行うべき警告として扱います。

よくある落とし穴と誤った正当化

つい選びがちな近道	それが誤りである理由／すべきこと
「暗黙的なスキームだから、どんな dt でも問題ありません。」	無条件の安定性は精度を意味しません。cfl_checker.py は --scheme implicit に対して緩和された（無限の）上限で stable: true を報告し、「精度を確認してください」という注記も追加します。大きな dt では時間誤差が依然として大きくなります。安定性だけでなく、精度を基準に dt を決めてください。
「クラッシュせずに100ステップ実行できたので、設定は安定しています。」	丸め誤差、保存則の損失、または限界的な Fo による後半での発散はよくあります。完了 ≠ 正しさです。実行結果を信頼する前に、metrics.fourier／metrics.cfl と limits.* を記録し、stable: true を確認してください。
「1D のフーリエ数の上限は 0.5 なので、Fo ≤ 0.5 なら安全です。」	陽的拡散の上限は Fo ≤ 1/(2d) で、2D では 0.25、3D では 0.167 です。実際の --dimensions を指定してください。cfl_checker.py は diffusion_limit を自動的に厳しくします。2D／3D で 0.5 を使うと不安定になります。
「剛性比が非常に大きいので、BDF／Radau を使いましょう。」	stiffness_ratio の値は、振動的／移流的／ハミルトン系のスペクトルでは誤解を招きます。imag_dominated と real_part_stiffness_ratio を確認してください。imag_dominated: true の場合、検出器は warning とともに stiff: false を返します。陰的な剛性ソルバーではなく、シンプレクティック法／リープフロッグ法、または CFL に合わせた A 安定スキームを優先してください。
「dt を小さくしたので、メッシュを細分化した後も以前のメッシュの dt を使えます。」	Fo ∝ dt/dx² です。dx を半分にすると Fo は4倍になります。細分化前の 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）
• 細工された .npy ファイルを介した任意コード実行を防ぐため、np.load() は allow_pickle=False を指定して呼び出されます（matrix_condition.py、stiffness_detector.py）
• スクリプトは標準出力（JSON 出力）にのみ書き込みます。エージェントが明示的に Write ツールを使用しない限り、ファイルは作成されません

ツールの制限

• Read：スクリプトのソース、参考資料、ユーザー設定ファイルの確認に使用します
• Bash：4つの 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）では、シェルで解釈可能な内容を含まない、構造化された解析可能な結果を生成します

制限事項

• CFL/Fourier チェックは陽解法のみです（陰解法は無条件安定です）
• フォン・ノイマン解析では、線形、定係数、周期境界条件（BCs）を仮定します
• スティフネスの検出には、ユーザーが指定する固有値の推定値が必要です。stiff 判定は、実際に減衰するモード（Re(λ) < 0）間のスケール分離に基づきます。大きさの比は引き続き stiffness_ratio として報告されます。虚軸が支配的なスペクトル（振動／移流／ハミルトン系）は imag_dominated と warning で示され、スティフとは分類されません。このような場合は、BDF/Radau ではなく、シンプレクティック法／リープフロッグ法、または CFL 制限に合わせて刻み幅を設定した A 安定なスキームを優先してください。

参考資料

• references/stability_criteria.md - 判定しきい値と公式
• references/common_pitfalls.md - よくある失敗モードと対処法
• references/scheme_catalog.md - 一般的なスキームの安定性特性

バージョン履歴

• v1.2.2（2026-06-24）：安定性判定を再現可能にし、誤った確信に基づく近道を防ぐため、「検証チェックリスト」（4 つのスクリプトの JSON 出力に関連付けられた、エビデンスに基づく 7 項目）と「よくある落とし穴と合理化」テーブル（6 行）を追加
• v1.2.0（2026-06-23）：相場場の計算例を修正（D=1.0 で実際に不安定な Fo=1.0）；eval ケース 1 の Fourier 値を修正；条件数のしきい値とエラーメッセージのドキュメントをスクリプトに合わせて調整；実部に基づくスティフネス分類と虚軸支配の検出を追加；文書化された入力／ファイルサイズ／次元の安全対策を適用
• v1.1.0（2024-12-24）：ドキュメント、判定ガイダンス、例を拡充
• v1.0.0：4 つの安定性解析スクリプトを含む初回リリース