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."
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
