Certified ODE twins¶

Twindex's polynomial ODE backend checks two parameter vectors against the same bounded-noise observations. A valid pair proves a lower bound L on the target parameter diameter. No upper bound U is claimed. A failed search proves no uniqueness. Initial values are exact polynomials in the parameters; uncertain initial values must be included as parameters.

The numerical generator needs pip install 'twindex[general]'. Replay needs only Python's standard library:

twindex check examples/study/ode/sir_pair.json examples/study/ode/pk_pair.json
# Intentional negative example: exits 1
twindex check examples/study/ode/pk_incompatible.json
# Rebuild the synthetic observations, certificates and timing table
PYTHONPATH=. python studies/ode_gallery.py

The SIR model observes only I, with S(0)=0.99 and I(0)=0.01:

S' = -beta S I
I' = beta S I - gamma I

Eight observations through t=4, with absolute noise 0.001, accept both (beta,gamma)=(0.5,0.2) and (0.72,0.4178). The certified beta gap is 11/50. The hidden removed compartment is determined by conservation. These parameters share the initial growth rate; their trajectories still need exact replay to prove compatibility at every listed observation.

The linear two-compartment pharmacokinetic model observes the central compartment C only, with C(0)=1 and P(0)=0:

C' = -(k10+k12) C + k21 P
P' = k12 C - k21 P

Eight observations through t=1, with absolute noise 0.01, accept both (k10,k12,k21)=(0.2,0.5,0.3) and (0.3,0.4,0.3). The certified k10 gap is 1/10. With later observations through t=4, the proposed second vector is proved incompatible. That result excludes this witness; it does not prove uniqueness or bound the diameter of the remaining compatible set.

Observations are synthetic reference-polynomial samples rounded to ten decimal places. They are stored as exact decimal rationals; the reference is itself checked against those same observations. These are mathematical regression examples, without claims about clinical or epidemiological data.

Candidate generation and limits¶

candidate_ode integrates with DOP853, samples Chebyshev nodes, rounds to dyadic coefficients, and converts exactly to monomials in each segment's local time. fit_witness_ode reuses the dense numerical solution while trying exact tube radii, then increases segment count and polynomial degree up to explicit limits. scales='auto' additionally tries positive rational state-magnitude weights. These weights are heuristic; the exact checker verifies their effect.

find_pair_ode adapts the model to audit_model, tries target endpoints and numerical null directions, rationalizes proposals, then certifies them independently. Its numerical audit remains numerical evidence even when a pair passes replay. Integration errors raise; exhausted sufficient bounds return INCONCLUSIVE. The implementation makes no stiff-system performance claim.

examples/study/ode/benchmark.json records generation plus replay wall time, segment counts, degree and success for these two examples. Timings depend on the machine and are a small regression benchmark, not a cross-method comparison.

A local sample run (four segments per witness, degree 10):

Case	Horizon	Generation + replay	Result	Certified L
Partially observed SIR	4	0.296 s	VALID pair	11/50
Central-only PK	1	0.021 s	VALID pair	1/10
PK later-observation witness	4	not timed	INVALID	none

Both proposed pairs succeeded (2/2); the intentionally incompatible witness was rejected (1/1). This table measures only these fixed proposals, not pair-search success probability. Regeneration records fresh timings in the JSON table.

Scaling and nonlinear stress cases¶

The same four-segment, degree-10 candidate for a three-state scaled cascade is INCONCLUSIVE with unit weights and VALID with automatic weights (1/1000, 1, 1/1000):

x' = -k x
y' = 1000 k x - k y
z' = k y/1000 - k z
(x,y,z)(0) = (0.001, 1, 0.001), k = 1

This is a change of state units in a linear cascade, not a stiff kinetic model. At the fixed mesh the unit-weight log norm fails mu*h < 1; diagonal scaling closes the sufficient bound without changing the polynomial pieces or data. Both certificate files are saved, including the honest INCONCLUSIVE result.

A FitzHugh–Nagumo example uses v' = v-v^3/3-w+drive, w' = 0.08*(v+0.7-0.8*w) with initial (-1,1). Four synthetic v observations through t=1 and epsilon=0.05 accept drive=0.5 and 0.51, proving L=1/100. This short horizon and loose noise do not establish performance on long trajectories, oscillations or stiff systems.

A local run measured 0.010 s for the unit-scale failure, 0.015 s for automatic scaling, and 0.050 s for the FHN pair, including candidate fitting and fresh replay. Reference/data generation separately took 0.004 s per model. Automatic scaling includes the failed unit-weight trials. See benchmark.json for the exact timings and fixed budgets; replay tests forbid numerical imports for all six gallery certificates.