numerical-stability — Claude Code skill: PDE stability analysis

---
name: numerical-stability
description: >
  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:
  last_deterministic_eval: "2026-10-20"
  agent_evaluation_status: "not_revalidated"
  author: HeshamFS
  version: "1.2.3"
  security_tier: high
  security_reviewed: true
  tested_with:
    - claude-code
  last_evaluated: "2026-06-24"
  eval_cases: 7
  last_reviewed: "2026-10-19"
  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 dt (grid refinement requires a coupled timestep reduction)
```

### 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

1. **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 needed
4. **Check stiffness** with `scripts/stiffness_detector.py` if multiple scales
5. **Analyze custom schemes** with `scripts/von_neumann_analyzer.py`
6. **Check conditioning** with `scripts/matrix_condition.py` for implicit solves
7. **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-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.3** (2026-10-19): Use a combined explicit stability budget, reject nonfinite inputs, and serialize undefined diagnostics explicitly.

- **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


### Scope of the stability result

The explicit checker assumes forward Euler, upwind advection, uniform Cartesian
centered diffusion, and linear decay. For an unsplit update it tests the sum
of normalized component loads, not separate conditions alone. The scalar
velocity must represent the total absolute directional speed for multidimensional
advection. The default reaction bound ensures positivity and is conservative
relative to absolute stability. Custom limits require a separately justified
component discretization. Negative diffusion/growth and unspecified implicit
integrators return no certificate. Refining dx requires reducing dt.

For final explanations, read `component_checks` and `scientific_interpretation`.
Do not confuse a combined failure with an individual-limit failure: C=0.4 and
Fo=0.32 each satisfy their separate limits (1 and 0.5), while C+2Fo=1.04 fails
the combined monotonicity bound. 0.32 is less than 0.5.
For a final report that does not introduce model-generated arithmetic claims,
use the [deterministic explanation renderer](../../GROUNDED_EXPLANATIONS.md)
with trusted, input-verified tool JSON. Unsupported or missing evidence must
remain unresolved.