Numerical Stability

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

HeshamFS Updated Sep 10, 2026

Install

Press C to copy

$ npx skillmds@latest add heshamfs/numerical-stability

Installs into ~/.claude/skills (or .claude/skills inside a project). Restart Claude Code to pick it up.

＋ Add to plugin

README badge

[![SkillMD](https://skillmd.com/badge/heshamfs/numerical-stability.svg)](https://skillmd.com/skills/heshamfs/numerical-stability)

Trust & safety

Audits pending

Not yet audited · commit 7cb44a9

Permissions this skill uses

Executes scripts
The instructions run shell commands or bundled script files.

Works with: Claude Code Claude.ai OpenAI Codex +

File contents

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

Script Outputs (JSON Fields)

Script	Key Outputs
`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`

Workflow

Identify dominant physics (advection vs diffusion vs reaction)
Run CFL checker with scripts/cfl_checker.py
Compare to limits and adjust dt if needed
Check stiffness with scripts/stiffness_detector.py if multiple scales
Analyze custom schemes with scripts/von_neumann_analyzer.py
Check conditioning with scripts/matrix_condition.py for implicit solves
Document the stability verdict and recommended time step

Conversational Workflow Example

User: My phase-field simulation is blowing up after 100 steps. I'm using explicit Euler with dx=0.01, dt=1e-4, and diffusivity D=1.0.

Agent workflow:

Check stability criteria:

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

Interpret results:
- Fourier number: Fo = 1.0 × 1e-4 / (0.01)² = 1.0
- 2D limit: Fo ≤ 0.25
- Violation: Fo = 1.0 > 0.25, unstable!
Recommend fix:
- 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.py status) 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
JSON output mode (--json) produces structured, parseable results without shell-interpretable content

Limitations

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
v1.1.0 (2024-12-24): Enhanced documentation, decision guidance, examples
v1.0.0: Initial release with 4 stability analysis scripts

---name: numerical-stabilitydescription: 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."---# Numerical Stability## GoalProvide 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)`.## Script Outputs (JSON Fields)| Script | Key Outputs ||--------|-------------|| `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` |## Workflow1. **Identify dominant physics** (advection vs diffusion vs reaction)2. **Run CFL checker** with `scripts/cfl_checker.py`3. **Compare to limits** and adjust `dt` if needed4. **Check stiffness** with `scripts/stiffness_detector.py` if multiple scales5. **Analyze custom schemes** with `scripts/von_neumann_analyzer.py`6. **Check conditioning** with `scripts/matrix_condition.py` for implicit solves7. **Document** the stability verdict and recommended time step## Conversational Workflow Example**User**: My phase-field simulation is blowing up after 100 steps. I'm using explicit Euler with dx=0.01, dt=1e-4, and diffusivity D=1.0.**Agent workflow**:1. Check stability criteria: ```bash python3 scripts/cfl_checker.py --dx 0.01 --dt 1e-4 --diffusivity 1.0 --dimensions 2 --json ```2. Interpret results: - Fourier number: `Fo = 1.0 × 1e-4 / (0.01)² = 1.0` - 2D limit: `Fo ≤ 0.25` - **Violation**: Fo = 1.0 > 0.25, unstable!3. Recommend fix: - 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```bash# Check CFL/Fourier for 2D diffusion-advectionpython3 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 stencilpython3 scripts/von_neumann_analyzer.py --coeffs 0.2,0.6,0.2 --dx 1.0 --nk 128 --json# Detect stiffness from eigenvalue estimatespython3 scripts/stiffness_detector.py --eigs=-1,-1000 --json# Check matrix conditioning for implicit systempython3 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.py` `status`) leave ample margin; well-conditioned-for-double FEM matrices (κ up to ~10⁶–10⁷) report `status: ok`.## Verification checklistDo 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- JSON output mode (`--json`) produces structured, parseable results without shell-interpretable content## Limitations- **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- **v1.1.0** (2024-12-24): Enhanced documentation, decision guidance, examples- **v1.0.0**: Initial release with 4 stability analysis scripts

HeshamFS/materials-simulation-skills/tree/main/skills/core-numerical/numerical-stability commit 7cb44a9e8e

Frequently asked questions

npx skillmds@latest add heshamfs/numerical-stability