# A periodic oscillator with no D-unstable child selection

Example for **A reaction-deletion-minimal mass-action oscillator without D-unstable child selections**, manuscript entry 37, SHA-256 `9bdeae52d9666dd64bb1f20cacb59053d81494d3a2102b340139b8f6d9a3bb95`.

The full reaction system can oscillate even though every eligible square child-selection matrix remains nonunstable under every positive diagonal scaling. This package reconstructs the literal mass-action model, verifies exact structural and crossing identities, freshly evaluates the rational interval certificate for the first Lyapunov coefficient, and computes a numerical periodic illustration. The exact sign and crossing support a conventional local Hopf theorem; the displayed decimal orbit is not an interval-validated witness.

## Run and editable study inputs

Python 3.11 or newer:

```sh
python -m venv .venv
# Windows: .venv\Scripts\activate
# macOS/Linux: source .venv/bin/activate
python -m pip install -r requirements.txt
python -m unittest -v test_example.py
python example.py --output outputs
```

The top of `example.py` exposes orbit amplitude, midpoint refinement counts, reference parameter for exact rate export, spectral range and sampling/tolerance settings. The default run takes about ten seconds. It writes a continuous numerical orbit, a spectral sweep, exact deletion covectors, rates, midpoint diagnostics, the regenerated rational Lyapunov enclosure and two figure pairs. The package is self-contained; no research workspace or Lean installation is required.

The local shooting amplitude is guarded to (0,0.05], and each nonlinear/linear solve must meet its stated criterion. An unsuccessful solve or positive-domain exit raises an error. This implementation guard is not a mathematically certified existence interval.

## Literal reactions and numerical representation

Species order is X1,X2,X3,X4, with the five irreversible reactions

```text
r1: X1 + 400 X4 -> 4 X3 + 398 X4
r2: 4 X2 + 2 X3 -> empty
r3: 375 X3 -> 371 X3 + 2 X4
r4: 2 X4 -> X1 + 6 X2 + X3
r5: empty -> 2 X3 + 2 X4
```

Classical mass action means v_j=k_j*product(x_i^Y_ij), using the full reactant coefficients, including species regenerated on the product side. The maximum reactant molecularity is 401. These are mathematical witness reactions, not realistic elementary chemistry, measured kinetics or a thermodynamically closed realization. Reaction r5 supplies inflow and r2 removes material.

`MassActionSource` owns reactant/product/net matrices, the positive circuit and exact falling-factorial derivative tensors. `RelativeReactor` owns a stationary state and flux, reconstructs log rates, evaluates the full nonlinear field/Jacobian, and integrates. `LogRates` accepts arbitrary positive paper-network rates represented by their logarithms and returns the unique equilibrium. `ChildCertificates`, `HopfCrossing`, the interval module and `PeriodicShooter` provide separate structural, bifurcation, rigorous arithmetic and numerical tasks.

The family is xbar(t)=(2,600,1500,8/t), vbar=(2,3,2,2,2), k_j=vbar_j/xbar^Y_j. In relative concentrations z=x/xbar, time is unchanged and

    dz/ds = diag(1/xbar) S diag(vbar) z^Y.

Because the high-order rate constants underflow in ordinary floating point, the model stores log(k_j). It never replaces an underflowed rate by zero. Dynamics evaluate the balanced difference using log1p and expm1. For z=1+r*u, the simulator integrates the smooth amplitude-scaled field F(1+r*u)/r; at r=0 it returns the linear field A*u. All four species remain present.

At arbitrary positive rates the equilibrium can be computed without a root search, in logarithms:

    log x4 = (log k5 - log k4)/2
    log x3 = (log k5 - log k3)/375
    log x2 = (log(3/2) + log k5 - log k2 - 2 log x3)/4
    log x1 = log k5 - log k1 - 400 log x4.

This follows from the positive flux circuit and the constant fifth reaction, not a numerical uniqueness assumption. The determinant is exactly 216000*k5^4/(x1*x2*x3*x4)>0, so no positive equilibrium has a zero eigenvalue. Uniqueness and nonsingularity do not imply attraction.

Reuse example:

```python
import numpy as np
from example import MassActionSource, RelativeReactor, LogRates

model = RelativeReactor.family(MassActionSource.paper(), t=1.0)
logs = np.array(model.log_rates.values)
logs[3] += np.log(1.01)  # independently increase k4 by one percent
changed = RelativeReactor.from_log_rates(LogRates(tuple(logs)))
print(changed.state, np.linalg.eigvals(changed.jacobian))
times, u = changed.integrate(
    initial=np.array([1., 0., 0., 0.]), period=10., amplitude=1e-4)
physical_concentrations = changed.state * (1 + 1e-4*u)
```

That new rate assignment need not lie in the attracting-orbit region. The exact certificate is bound to the paper's source and one-parameter family; changing kinetic orders must not silently reuse it. The interval routine checks the literal characteristic identities before accepting its eigenvector argument.

`reference_rates.csv` contains exact rational rates at the separately editable reference parameter (default t=1), and log rates at the numerical orbit in a clearly separate column. Those columns do not describe the same parameter value.

## Structural certificates

Child selections injectively assign each chosen species a reaction in which it is a reactant. They see net stoichiometry and eligibility, but not the sizes of kinetic orders. Removing 398 copies of X4 from both sides of r1 and 371 copies of X3 from both sides of r3 gives an unpadded skeleton with identical net stoichiometry and child selections. The plotted skeleton is numerically stable along this family; its universal mass-action stability is the companion paper's theorem.

All 24 nonempty children are covered by six matrices. Five have exact diagonal-weight energy certificates -(Pi*M+M^T*Pi)>=0, checked through every principal minor. The last has scaled characteristic polynomial lambda*(lambda+4*d1)*(lambda+4*d2+2*d3). Restrictions inherit noninstability. Singular children retain zero eigenvalues, so the claim is D-nonunstability, not strict D-stability.

The net matrix has rank four and kernel spanned by the positive vector (2,3,2,2,2). Every proper column subset is independent. For each deleted reaction j, the code constructs a rational covector w satisfying w*S_l=1 for every retained reaction. Along a positive trajectory of the deleted system,

    d(w*x)/ds = sum of retained reaction fluxes.

This is strictly positive if at least one retained rate is positive, ruling out positive equilibria, returns and recurrence under arbitrary nonnegative retuning. If all rates vanish, every state is stationary and every solution is constant. The mixed-sign covectors are not positive mass bounds and do not establish global existence. Deletion removes balanced positive operation itself; minimality is relative to these five reactions, not all possible oscillators.

## Crossing and the rational Lyapunov sign

The characteristic quartic has positive coefficients and third Hurwitz determinant 3*h(t)/7812500, where

    h(t)=2825760+242612196*t+3570069381*t^2-4001047375*t^3.

Exact signs, Descartes' rule and the derivative identity isolate its unique positive simple root tH in (1/2,1). At tH the quartic factors into a simple imaginary pair and a stable real quadratic. The exact Hurwitz derivative identity gives positive crossing speed. Numerically tH≈0.956453654736546, omegaH≈0.100998757759762, and alpha'(tH)≈0.004291967. The original and relative Jacobians are similar; normalization does not change stability.

`interval_certificate.py` adapts the arithmetic approach of the paper workspace's `amplitude_enclosure.py` (source SHA-256 recorded in that module). It bisects the crossing with exact rational signs, encloses frequency by squared endpoint inequalities, and evaluates the literal first/second/third derivative tensors from falling factorials. Real operations round outward to a 2^-100 dyadic grid. Complex intervals and all 14 Gaussian-elimination pivots have exact denominator-separation checks; a nonzero left/right pairing normalizes the eigenvectors. The exact quartic factorization and invertible leading block justify the omitted eigenvector equations through the Schur complement.

The first Lyapunov coefficient uses the invariant expression with derivative tensors B,C (no factorials absorbed). It is initially computed with q4=1, then divided by ||q||^2 to obtain the unit-norm convention. The rational endpoints satisfy

    -23/1000 < l1 < -22/1000,

with numerical center about -0.02239969743507. All acceptance decisions use rationals, including source binding, root signs, pivots and the final enclosure. Together with the exact crossing and a conventional nondegenerate Hopf theorem, the negative sign gives an attracting hyperbolic positive periodic orbit for t>tH sufficiently close and persistence on some open set of positive rates. The theorem supplies neither an explicit endpoint of that set nor a certified orbit at a printed decimal parameter. This rational computation and Hopf application are not Lean-compiled.

## Numerical orbit and mean identities

`PeriodicShooter` fixes an initial center amplitude in a numerical real eigenbasis, then solves four return equations for two stable coordinates, parameter and period. Its implicit midpoint refinements reproduce the paper's periods 62.2226143528, 62.2132421634 and 62.2108996453 at 128, 256 and 512 steps. The period error decreases at the expected second-order scale.

A separate continuous Radau shooting solve gives t≈0.956506669872689 and period≈62.2101188529, with return residual about 2.55e-13 in the shooting basis. BDF independently follows the returned initial state and differs in relative concentration by about 6.83e-11. The minimum relative concentration is about 0.980019. These are numerical diagnostics, not a validated continuous orbit or a numerical proof of stability.

For every true positive periodic orbit, integrating S*v over a period forces its mean flux to equal the equilibrium flux: (k5,3*k5/2,k5,k5,k5). Consequently the four normalized reaction monomials have mean one, including mean(z4^2)=1. A nonconstant periodic orbit must have varying z4, so mean(z4)^2+Var(z4)=1 implies mean(z4)<1. Products of means are not substituted for means of products.

The continuous numerical quadrature gives mean(z4)-1≈-6.41367e-10 and variance≈1.28273e-9, with mean-square identity defect around -3.84e-18. The tiny mean shift is computed from deviations, avoiding subtraction of two rounded near-unit averages. Midpoint mean-flux identities are also a consequence of the conservative algebra of a returned discrete scheme; they are useful implementation checks, not independent existence evidence. The orbit figure scales the fourth coordinate by 10^6 and the third by 100 only for visibility.

## Validation and scope

Seven test groups check literal padding and arbitrary log-rate equilibrium reconstruction, every proper stoichiometric subset and deletion identity, all child certificates, exact crossing identities and inertia, all ordered jets through degree three against symbolic derivatives, outward arithmetic and the Lyapunov enclosure including a rejected altered-source control, and numerical refinement/return/mean consistency. Clean-package replay compares scientific JSON/CSV. Lean is not rerun. No global attraction, basin, minimal kinetic order, low-molecularity realization or experimental calibration is asserted.

Licensing: MIT is proposed for this example, pending owner confirmation; this statement does not grant a license. The manuscript retains its own terms.