Journey 03 — Stiff chemistry (Robertson)

Integrate stiff Robertson kinetics with analytic Jacobian, non-negativity, and mass-conservation wrappers — the production stiff-ODE pattern.

Prerequisites: Journey 02 — you can read a result dict. Scenario: Robertson kinetics | rates span 1e9 | ndim=3

Key APIs

solvsrk.run(..., jac_fn=...)
cfg.handoff_nonneg
cfg.handoff_conservation_sum
cfg.cvode_enforce_nonneg

Concepts to internalize

Stiffness — rates span ~1e9; explicit RK45 stalls or fails (see journey 06).
Analytic Jacobian — stiff integrator uses it directly; finite differences are too expensive.
Physical wrappers — prevent negative concentrations and enforce sum(y)=1 at handoff and in CVODE.
Sanity checks — always verify conservation and non-negativity on success.

Important code

The guided tour script wraps this pattern. Copy the production shape into your pipeline:

def robertson_rhs(t, y):
    y1, y2, y3 = y
    return np.array([
        -0.04*y1 + 1e4*y2*y3,
        0.04*y1 - 1e4*y2*y3 - 3e7*y2**2,
        3e7*y2**2,
    ])
 
def robertson_jac(t, y):
    y1, y2, y3 = y
    return np.array([...])  # row-major 3x3
 
cfg = solvsrk.robertson_config()
result = solvsrk.run(cfg, 3, 0, 40, (1,0,0), robertson_rhs, jac_fn=robertson_jac)
assert abs(sum(result['y_final']) - 1.0) < 1e-6  # mass conserved

Act-by-act walkthrough

Act 1 — The system

Three chemical species. y1 leaks slowly to y2; y2 is consumed almost instantly by a fast quadratic reaction into y3. The mix of a very fast and a very slow mode is what 'stiff' means — and what defeats a plain explicit stepper. We integrate with SolvScout-tuned knobs (journey 04 shows untuned vs tuned); install SolvScout for live tokens.

Act 2 — Physical constraints as config

Concentrations cannot go negative and the three species must sum to 1 (mass is conserved). SolvSRK can enforce both at the internal handoff and inside the stiff integrator, so numerical noise never produces an unphysical state. These wrappers sit on top of the tuned knob set.

Field	Value	Notes
`handoff_nonneg`	1
`handoff_conservation_sum`	1.0
`cvode_enforce_nonneg`	1

Act 3 — Solve with an analytic Jacobian

For stiff problems an analytic Jacobian pays for itself — the stiff integrator uses it directly instead of finite-differencing. Pass jac_fn and SolvSRK marshals it row-major to the C core.

Field	Value	Notes
`knob_source`	SolvScout tuned (solvsrk, high)
`survived`	True
`failure_mode`	ok
`t_reached`	40
`y_final`	[0.715827, 9.18553e-06, 0.284164]
`total_rhs_evals`	795
`front_phase`	steps=0 rhs=0
`stiff_phase`	steps=613 rhs=795
`t_handoff`	0
`wall_s`	0.002629

Act 4 — Sanity check

When the integration succeeds, the three concentrations should still sum to ~1 — the conservation wrapper keeps mass balanced across the stiff tail. Check sum(y) and that every component stayed non-negative.

Field	Value	Notes
`y1`	0.715827
`y2`	9.18553e-06
`y3`	0.284164
`sum`	1.0000000000
`conservation_error`	8.882e-16

Guided tour output

Verified console capture from 03_stiff_chemistry.py:

Robertson kinetics | rates span 1e9 | ndim=3

-- Act 1 - The system --
  -> Three chemical species. y1 leaks slowly to y2; y2 is consumed almost instantly by a fast quadratic reaction into y3. The mix of a very fast and a very slow mode is what 'stiff' means — and what defeats a plain explicit stepper. We integrate with SolvScout-tuned knobs (journey 04 shows untuned vs tuned); install SolvScout for live tokens.

-- Act 2 - Physical constraints as config --
  -> Concentrations cannot go negative and the three species must sum to 1 (mass is conserved). SolvSRK can enforce both at the internal handoff and inside the stiff integrator, so numerical noise never produces an unphysical state. These wrappers sit on top of the tuned knob set.
  handoff_nonneg: 1
  handoff_conservation_sum: 1.0
  cvode_enforce_nonneg: 1

-- Act 3 - Solve with an analytic Jacobian --
  -> For stiff problems an analytic Jacobian pays for itself — the stiff integrator uses it directly instead of finite-differencing. Pass jac_fn and SolvSRK marshals it row-major to the C core.
  knob_source: SolvScout tuned (solvsrk, high)
  survived: True
  failure_mode: ok
  t_reached: 40
  y_final: [0.715827, 9.18553e-06, 0.284164]
  total_rhs_evals: 795
  front_phase: steps=0 rhs=0
  stiff_phase: steps=613 rhs=795
  t_handoff: 0
  wall_s: 0.002629

-- Act 4 - Sanity check --
  -> When the integration succeeds, the three concentrations should still sum to ~1 — the conservation wrapper keeps mass balanced across the stiff tail. Check sum(y) and that every component stayed non-negative.
  y1: 0.715827
  y2: 9.18553e-06
  y3: 0.284164
  sum: 1.0000000000
  conservation_error: 8.882e-16

-- Run complete --
  next_journey: 04_configuration.py

Run it

cd python
# Guided tour — narrated scenario walkthrough (recommended first run):
python examples/journeys/03_stiff_chemistry.py
 
# Production one-liner — same API call you ship:
python examples/journeys/03_stiff_chemistry.py

What comes next

Journey 04 compares untuned vs SolvScout-tuned knobs and demonstrates run_profile().

Navigation

← Journey 02 — Core path (Rössler) · Journey 04 — Configuration and SolvScout hand-off →

02 Core path (Rössler)04 Configuration and SolvScout hand-off