Analyze numerical stability for time-dependent PDE simulations — check CFL and Fourier criteria, perform von Neumann stability analysis, detect stiffness, evaluate matrix conditioning, and recommend explicit vs implicit time-stepping schemes. Use when selecting time steps, diagnosing numerical blow-up or solver divergence, checking convergence criteria, or evaluating scheme stability for advection, diffusion, or reaction problems, even if the user doesn't explicitly mention "stability" or "CFL."
التثبيت
التثبيت باستخدام Codex أو Claude انسخ هذا Prompt والصقه في Codex أو Claude أو مساعد آخر ليراجع صفحة Skill ويثبّتها لك.
Analyze numerical stability for time-dependent PDE simulations — check CFL and Fourier criteria, perform von Neumann stability analysis, detect stiffness, evaluate matrix conditioning, and recommend explicit vs implicit time-stepping schemes. Use when selecting time steps, diagnosing numerical blow-up or solver divergence, checking convergence criteria, or evaluating scheme stability for advection, diffusion, or reaction problems, even if the user doesn't explicitly mention "stability" or "CFL."
allowed-tools
Read, Bash, Write, Grep, Glob
metadata
{"author":"HeshamFS","version":"1.2.2","security_tier":"high","security_reviewed":true,"tested_with":["claude-code"],"last_evaluated":"2026-06-24","eval_cases":4,"last_reviewed":"2026-06-23","standards":["Courant-Friedrichs-Lewy (CFL) condition (Courant, Friedrichs, Lewy 1928)","von Neumann (Fourier) stability analysis; amplification factor and Fourier limit Fo ≤ 1/(2d)","Dahlquist A-stability / L-stability theory for time integrators","BDF (Gear) and Radau IIA / Rosenbrock methods for stiff systems","IEEE 754 double precision (conditioning thresholds, κ ≈ 1/eps)"]}
Numerical Stability
Goal
Provide a repeatable checklist and script-driven checks to keep time-dependent simulations stable and defensible.
Requirements
Python 3.10+
NumPy (for matrix_condition.py and von_neumann_analyzer.py)
See scripts/requirements.txt for dependencies
Inputs to Gather
Input
Description
Example
Grid spacing dx
Spatial discretization
0.01 m
Time step dt
Temporal discretization
1e-4 s
Velocity v
Advection speed
1.0 m/s
Diffusivity D
Thermal/mass diffusivity
1e-5 m²/s
Reaction rate k
First-order rate constant
100 s⁻¹
Dimensions
1D, 2D, or 3D
2
Scheme type
Explicit or implicit
explicit
Decision Guidance
Choosing Explicit vs Implicit
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
Stability Limit Quick Reference
Physics
Number
Explicit Limit (1D)
Formula
Advection
CFL
C ≤ 1
C = v·dt/dx
Diffusion
Fourier
Fo ≤ 0.5
Fo = D·dt/dx²
Reaction
Reaction
R ≤ 1
R = k·dt
Multi-dimensional correction: For d dimensions, diffusion limit is Fo ≤ 1/(2d).
Reduce dt to 2.5e-5 (the tool's recommended_dt, giving Fo = 0.25)
Or increase dx, or switch to implicit
Pre-Simulation Stability Checklist
Identify dominant physics and nondimensional groups
Compute CFL/Fourier/Reaction numbers with cfl_checker.py
If explicit and limit violated, reduce dt or change scheme
If stiffness ratio > 1000, select implicit/stiff integrator
For custom schemes, verify amplification factor ≤ 1
Document stability reasoning with inputs and outputs
CLI Examples
# 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
Error Handling
Error
Cause
Resolution
dx and dt must be positive
Zero or negative values
Provide valid positive numbers
No stability criteria applied
Missing velocity/diffusivity
Provide at least one physics parameter
Matrix not found: <path>
Invalid path
Check matrix file exists
Could not compute eigenvalues
Singular or ill-formed matrix
Check matrix validity
Interpretation Guidance
Scenario
Meaning
Action
stable: true
All checked criteria satisfied
Proceed with simulation
stable: false
At least one limit violated
Reduce dt or change scheme
stable: null
No criteria could be applied
Provide more physics inputs
Stiffness ratio > 1000
Problem is stiff
Use implicit integrator
Condition number > 10⁸
Poorly-conditioned (status: poorly-conditioned)
Preconditioning likely needed
Condition number > 10¹⁰
Ill-conditioned (status: ill-conditioned)
Use scaling/preconditioning
Conditioning thresholds assume IEEE double precision: solving loses roughly log10(κ) significant digits, and a matrix becomes numerically singular near κ ≈ 1/eps ≈ 4.5e15. The > 10⁸ / > 10¹⁰ cutoffs (matching matrix_condition.pystatus) leave ample margin; well-conditioned-for-double FEM matrices (κ up to ~10⁶–10⁷) report status: ok.
Verification checklist
Do not declare a scheme/time step "stable" until each applicable item below is satisfied with a recorded value from the scripts, not a mental estimate.
Ran cfl_checker.py --json and recorded metrics.cfl, metrics.fourier, and/or metrics.reaction against the reported limits.* (explicit defaults: CFL ≤ 1, Fo ≤ 1/(2d), R ≤ 1), with stable: true and the intended criteria present in criteria_applied.
Confirmed criteria_applied is non-empty and stable is not null — i.e. at least one physics parameter (--velocity/--diffusivity/--reaction-rate) was actually supplied so a real check ran, not a silent no-op.
Used the smallest grid spacing across all directions for --dx (anisotropic grids: smallest dx/dy/dz) and re-ran cfl_checker.py after any mesh refinement, since Fo ∝ dt/dx² makes the limit highly sensitive to dx.
If the chosen dt is below the limit, recorded the tool's recommended_dt (and the --safety factor used) so the margin to the stability boundary is explicit and reproducible.
For any custom/non-standard update stencil, ran von_neumann_analyzer.py --json and confirmed results.max_amplification ≤ 1 (stable: true); noted k_at_max and resolved any even-length-stencil warning.
For multi-scale systems, ran stiffness_detector.py --json and recorded real_part_stiffness_ratio, imag_dominated, and stiff; only chose BDF/Radau when stiff: true on the real-part ratio (not on magnitude alone) and there is no warning.
For implicit solves, ran matrix_condition.py --json and recorded condition_number and status; treated poorly-conditioned (>1e8) / ill-conditioned (>1e10) as a flag to scale/precondition before trusting the solve.
Common pitfalls & rationalizations
Tempting shortcut
Why it's wrong / what to do
"Implicit scheme, so any dt is fine."
Unconditional stability is not accuracy. cfl_checker.py reports stable: true with relaxed (infinite) limits for --scheme implicit and even adds a note to "check accuracy" — a large dt still ruins temporal error. Size dt for accuracy, not just stability.
"It ran 100 steps without crashing, so the setup is stable."
Late-time blow-up from round-off, conservation loss, or marginal Fo is common. Completion ≠ correctness — record metrics.fourier/metrics.cfl vs limits.* and confirm stable: true before trusting the run.
"The 1D Fourier limit is 0.5, so Fo ≤ 0.5 is safe."
The explicit diffusion limit is Fo ≤ 1/(2d) — 0.25 in 2D, 0.167 in 3D. Pass the real --dimensions; cfl_checker.py tightens diffusion_limit automatically, and using 0.5 in 2D/3D is an instability.
"Stiffness ratio is huge, so use BDF/Radau."
The magnitude stiffness_ratio is misleading for oscillatory/advective/Hamiltonian spectra. Check imag_dominated and real_part_stiffness_ratio: if imag_dominated: true the detector returns stiff: false with a warning — prefer symplectic/leapfrog or a CFL-sized A-stable scheme, not implicit stiff solvers.
"I tightened dt, so the old mesh's dt still works after refining dx."
Fo ∝ dt/dx²: halving dx quadruples Fo. Reusing a pre-refinement dt reintroduces a violation. Recompute with cfl_checker.py after every mesh change.
"The matrix solved, so its conditioning is fine."
A solve can return numbers while silently losing ~log10(κ) digits. Record condition_number/status from matrix_condition.py; poorly-conditioned/ill-conditioned means scale or precondition before trusting the result.
Security
Input Validation
dx, dt, and safety are validated as finite positive numbers before any computation; velocity, diffusivity, and reaction_rate, when supplied, are validated as finite
--dimensions is restricted to {1, 2, 3} (cfl_checker.py raises and exits 2 otherwise)
Comma-separated eigenvalue lists (--eigs) are capped at 10,000 entries and validated as finite numbers
Stencil coefficient lists (--coeffs) are capped at 10,000 entries and validated as finite floats
File Access
matrix_condition.py reads a single matrix file (.npy or text) specified by --matrix; no directory traversal beyond the given path
Matrix/Jacobian files are rejected if they exceed 500 MB before parsing (matrix_condition.py, stiffness_detector.py)
np.load() is called with allow_pickle=False to prevent arbitrary code execution via crafted .npy files (matrix_condition.py, stiffness_detector.py)
Scripts write only to stdout (JSON output); no files are created unless the agent explicitly uses the Write tool
Tool Restrictions
Read: Used to inspect script source, references, and user configuration files
Bash: Used to execute the four Python analysis scripts (cfl_checker.py, von_neumann_analyzer.py, matrix_condition.py, stiffness_detector.py) with explicit argument lists
Write: Used to save analysis results or generated reports; writes are scoped to the user's working directory
Grep/Glob: Used to locate relevant files and search references
Safety Measures
No eval(), exec(), or dynamic code generation
All subprocess calls use explicit argument lists (no shell=True)
Matrix dimension limit (100,000 per side) in matrix_condition.py prevents memory exhaustion
Explicit schemes only for CFL/Fourier checks (implicit is unconditionally stable)
Von Neumann analysis assumes linear, constant-coefficient, periodic BCs
Stiffness detection requires eigenvalue estimates from user. The stiff verdict is based on scale separation among genuinely decaying modes (Re(λ) < 0); the magnitude ratio is still reported as stiffness_ratio. Imaginary-axis-dominated spectra (oscillatory/advection/Hamiltonian) are flagged via imag_dominated and a warning, and are NOT classified as stiff — for those, prefer symplectic/leapfrog or an A-stable scheme sized by the CFL limit, not BDF/Radau.
References
references/stability_criteria.md - Decision thresholds and formulas
references/common_pitfalls.md - Frequent failure modes and fixes
references/scheme_catalog.md - Stability properties of common schemes
Version History
v1.2.2 (2026-06-24): Added a "Verification checklist" (7 evidence-based items tied to the four scripts' JSON outputs) and a "Common pitfalls & rationalizations" table (6 rows) to make stability verdicts reproducible and guard against false-confidence shortcuts
v1.2.0 (2026-06-23): Corrected the phase-field worked example (genuinely unstable Fo=1.0 with D=1.0); fixed eval case 1 Fourier value; reconciled conditioning thresholds and error-message docs with the scripts; real-part-based stiffness classification with imaginary-dominated detection; enforced documented input/file-size/dimension safeguards
