Numerical Stability

Goal

Provide a repeatable checklist and script-driven checks to keep time-dependent simulations stable and defensible.

Requirements

• Python 3.8+
• NumPy (for matrix_condition.py and von_neumann_analyzer.py)
• See

  scripts/requirements.txt

  for dependencies

Inputs to Gather

dx

0.01 m

dt

1e-4 s

v

1.0 m/s

D

1e-5 m²/s

k

100 s⁻¹

2

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

C = v·dt/dx

Fo = D·dt/dx²

R = k·dt

Multi-dimensional correction: For d dimensions, diffusion limit is

Fo ≤ 1/(2d)

Script Outputs (JSON Fields)

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

results.recommendation

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=1e-3.

Agent workflow:

1. Check stability criteria:

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

2. Interpret results:
   • Fourier number:

     Fo = 1e-3 × 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

     (to get 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

dx and dt must be positive

No stability criteria applied

Matrix file not found

Could not compute eigenvalues

Interpretation Guidance

stable: true

stable: false

stable: null

Security

Input Validation

• All numeric parameters (

  dx

  ,

  dt

  ,

  velocity

  ,

  diffusivity

  ,

  dimensions

  ) are validated as finite positive numbers before any computation

• --dimensions

  is restricted to

  {1, 2, 3}

• Comma-separated eigenvalue lists (

  --eigs

  ) are capped at 10,000 entries and validated as finite numbers
• Stencil coefficient lists (

  --coeffs

  ) are length-limited and validated as finite floats

File Access

• matrix_condition.py

  reads a single matrix file (

  .npy

  format) specified by

  --matrix

  ; no directory traversal beyond the given path
• Matrix files are rejected if they exceed 500 MB before parsing

• np.load()

  is called with

  allow_pickle=False

  to prevent arbitrary code execution via crafted

  .npy

  files
• 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 limits (100,000 per dimension) prevent 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

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