Example code
Every reaction in this four-species model consumes at most two molecules. A child selection chooses one distinct consuming reaction for each species in a subset. None of the 25 child matrices develops a growing mode when its columns are multiplied by positive factors. Nevertheless, small concentration disturbances can grow in the full mass-action system at its positive equilibrium.
The example makes that separation explicit: a reusable reaction model reconstructs rates from a stationary flux and composition, a structural layer certifies every child, and a kinetic layer evaluates the full Jacobian (the matrix governing small concentration disturbances). Changing operating conditions can switch stability while leaving all the child certificates unchanged.
At the manuscript's default rates , the equilibrium is . An exact shifted-polynomial calculation places an unstable eigenvalue's real part between 2 and 4. Numerical evaluation gives approximately . At the stable control , the same source has a strictly stable equilibrium. Four child matrices are singular, so their certified property is noninstability rather than strict stability.
The operating map varies circulation relative to return flux and the inverse concentration of the fast species . The code reproduces all allowable positive stationary fluxes and the exact unstable ray , . A Schur reduction explains the large- threshold near : the full Jacobian (the matrix governing small concentration disturbances) combines reaction sensitivities that no individual child selection records. The shaded grid is a numerical illustration, not a certified bifurcation diagram.
The package also supplies a feed-and-dilution reactor, keeping the chosen equilibrium while shifting the Jacobian by . At the default witness, every remains unstable. The code checks all 121 child selections of this twelve-reaction completion using the paper's triangular-block transport argument. A separate interface reconstructs the moving equilibrium after independent changes to the original six rate constants.
The short trajectory illustrates local departure, checked against the linear mode and a second stiff solver. It does not establish a periodic orbit. These are dimensionless abstract-species kinetics: the active source has no positive mass vector, and one reaction produces nine molecules. The example therefore demonstrates reactant-bimolecular instability, not an atom-balanced mechanism or a two-sided bimolecular result. Seven scientific test groups and saved exact certificates accompany the reusable model; Lean is not rerun.
Python source
"""Low-order mass-action instability invisible to every child selection.
Exact structural and spectral checks are separate from numerical trajectories.
"""
from __future__ import annotations
# EDITABLE INPUTS: dimensionless manuscript operating conditions --------------
CIRCULATION_RATIO = '100' # T>1
INVERSE_D_CONCENTRATION = '100' # L>0; equilibrium (1,1,1,1/L)
THROUGHPUT = '1' # s>0; rescales all original rates and time
DILUTION_RATE = '1' # completion keeps the selected equilibrium
PERTURBATION_SIZE = 1e-7
TRAJECTORY_END = 0.8
TRAJECTORY_POINTS = 401
RATE_PERTURBATION_FRACTION = 0.001
PARAMETER_T_RANGE = (1.05, 120.0)
PARAMETER_L_RANGE = (0.1, 200.0)
MANUSCRIPT_SHA256 = '12479870603ad4aafcc2a7cef8d0c32f22a78f07e210ba83c4f904c249063dbf'
# ---------------------------------------------------------------------------
import argparse
from dataclasses import dataclass
import hashlib
import itertools
import json
from pathlib import Path
import platform
import csv
import numpy as np
import sympy as sp
from scipy.integrate import solve_ivp
def Q(x):
return sp.Rational(str(x))
def vec(xs):
return sp.Matrix([sp.sympify(x) if isinstance(x, sp.Basic) else Q(x) for x in xs])
def encoded(xs):
return [str(sp.simplify(x)) for x in xs]
@dataclass(frozen=True)
class Reaction:
name: str
reactant: tuple[int, ...]
product: tuple[int, ...]
class Source:
"""Literal integer complexes; deterministic homodimer event rate is k*x**2."""
def __init__(self, species, reactions):
self.species = tuple(species); self.reactions = tuple(reactions)
self.n = len(self.species); self.m = len(self.reactions)
if not self.n or not self.m or len(set(self.species)) != self.n:
raise ValueError('Distinct species and at least one reaction required.')
for r in self.reactions:
if len(r.reactant) != self.n or len(r.product) != self.n or any(
not isinstance(v, int) or v < 0 for v in (*r.reactant, *r.product)):
raise ValueError('Complexes must be matching nonnegative integer vectors.')
self.Y = sp.ImmutableMatrix.hstack(*[sp.ImmutableMatrix(r.reactant) for r in self.reactions])
self.P = sp.ImmutableMatrix.hstack(*[sp.ImmutableMatrix(r.product) for r in self.reactions])
self.S = self.P-self.Y
def monomials(self, x):
return sp.Matrix([sp.prod(x[i]**self.Y[i, j] for i in range(self.n)) for j in range(self.m)])
def reconstruct(self, equilibrium, flux):
x = vec(equilibrium); v = vec(flux)
if len(x) != self.n or len(v) != self.m or any(a <= 0 for a in (*x, *v)):
raise ValueError('Positive matching equilibrium and flux vectors required.')
if (self.S*v).applyfunc(sp.simplify) != sp.zeros(self.n, 1):
raise ValueError('Flux must be stationary, S*v=0.')
mono = self.monomials(x)
return MassAction(self, [sp.simplify(v[j]/mono[j]) for j in range(self.m)], x)
def molecularity(self):
return {'reactants': [int(sum(self.Y[:, j])) for j in range(self.m)],
'products': [int(sum(self.P[:, j])) for j in range(self.m)],
'catalyst_free': all(self.Y[i, j]*self.P[i, j] == 0 for i in range(self.n) for j in range(self.m))}
class MassAction:
def __init__(self, source, rates, equilibrium=None):
self.source = source; self.k = vec(rates)
if len(self.k) != source.m or any(k <= 0 for k in self.k):
raise ValueError('One positive rate constant per reaction required.')
self.equilibrium = None if equilibrium is None else vec(equilibrium)
self.x = sp.Matrix(sp.symbols(f'x0:{source.n}'))
self.symbolic_flux = sp.matrix_multiply_elementwise(self.k, source.monomials(self.x))
self.symbolic_field = source.S*self.symbolic_flux
self.symbolic_jacobian = self.symbolic_field.jacobian(self.x)
self._field = sp.lambdify([tuple(self.x)], self.symbolic_field, 'numpy')
self._jac = sp.lambdify([tuple(self.x)], self.symbolic_jacobian, 'numpy')
if self.equilibrium is not None:
if len(self.equilibrium) != source.n or any(v <= 0 for v in self.equilibrium):
raise ValueError('Equilibrium must be positive and match species.')
if self.exact_field(self.equilibrium).applyfunc(sp.simplify) != sp.zeros(source.n, 1):
raise ValueError('Supplied state is not an exact equilibrium.')
def exact_field(self, x):
return self.symbolic_field.subs(dict(zip(self.x, x)))
def jacobian(self, x=None):
state = self.equilibrium if x is None else vec(x)
if state is None:
raise ValueError('Supply a state for the Jacobian.')
return self.symbolic_jacobian.subs(dict(zip(self.x, state))).applyfunc(sp.simplify)
def field(self, x):
return np.asarray(self._field(np.asarray(x, float)), float).reshape(-1)
def integrate(self, initial, times, method='Radau'):
initial = np.asarray(initial, float); times = np.asarray(times, float)
if initial.shape != (self.source.n,) or np.any(initial <= 0) or not np.all(np.isfinite(initial)):
raise ValueError('Positive finite initial concentrations required.')
if len(times) < 2 or times[0] != 0 or not np.all(np.diff(times) > 0):
raise ValueError('Strictly increasing times starting at zero required.')
solution = solve_ivp(lambda t, x: self.field(x), (0, times[-1]), initial, t_eval=times,
jac=lambda t, x: np.asarray(self._jac(x), float), method=method,
rtol=1e-11, atol=1e-14)
if not solution.success or np.min(solution.y) <= 0:
raise ArithmeticError('Numerical integration failed or left the positive region.')
return solution.y.T
class OperatingFamily:
def __init__(self):
self.source = Source(('A', 'B', 'C', 'D'), (
Reaction('A+D -> 2C', (1,0,0,1), (0,0,2,0)),
Reaction('B+C -> empty', (0,1,1,0), (0,0,0,0)),
Reaction('2C -> D', (0,0,2,0), (0,0,0,1)),
Reaction('2D -> A+4B+4C', (0,0,0,2), (1,4,4,0)),
Reaction('empty -> A', (0,0,0,0), (1,0,0,0)),
Reaction('empty -> D', (0,0,0,0), (0,0,0,1))))
def flux(self, T, s=1):
T, s = Q(T), Q(s)
if T <= 1 or s <= 0:
raise ValueError('Positive stationary cone requires T>1 and s>0.')
return s*sp.Matrix([T, 4, T, 1, T-1, 2])
def model(self, T=100, L=100, s=1):
L = Q(L)
if L <= 0: raise ValueError('L must be positive.')
return self.source.reconstruct([1, 1, 1, 1/L], self.flux(T, s))
def from_rates(self, rates):
"""Closed-form positive equilibrium for any six positive original rates.
This follows from the complete stationary flux cone, rather than a
nonlinear root solver. Stability still needs a separate calculation.
"""
k = vec(rates)
if len(k) != 6 or any(v <= 0 for v in k): raise ValueError('Six positive rates required.')
s = k[5]/2; T = 1+k[4]/s
D = sp.sqrt(s/k[3]); C = sp.sqrt(s*T/k[2])
state = sp.Matrix([s*T/(k[0]*D), 4*s/(k[1]*C), C, D])
return MassAction(self.source, k, state)
@staticmethod
def coefficients(T, L):
return ((T+4)*L+5*T+8, (14*T+32)*L+4*T*(T+6), 112*T*L+16*T*T, 64*T*T*L)
def exact_family_certificate(self):
T, L, u = sp.symbols('T L u', positive=True)
a = self.coefficients(T, L); H = sp.expand(Quartic.hurwitz(a))
# Derive the Jacobian from flux and the mass-action derivative identity.
v = sp.Matrix([T, 4, T, 1, T-1, 2])
J = self.source.S*sp.diag(*v)*self.source.Y.T*sp.diag(1,1,1,L)
if any(sp.expand(x-y) != 0 for x, y in zip(J.charpoly().all_coeffs()[1:], a)):
raise AssertionError('Family characteristic polynomial mismatch.')
ray = sp.Poly(H.subs({T:100, L:100+u}), u)
if not all(v < 0 for v in ray.all_coeffs()): raise AssertionError('Ray sign certificate failed.')
reduced = (J[:3,:3]-J[:3,3]*(1/J[3,3])*J[3,:3]).applyfunc(sp.factor)
if any(L in x.free_symbols for x in reduced): raise AssertionError('Fast reduction retained L.')
cubic = reduced.charpoly().all_coeffs()[1:]
threshold = (41+sp.sqrt(2577))/4
return {'family_coefficients': encoded(a), 'H': str(H),
'negative_ray_coefficients_descending': encoded(ray.all_coeffs()),
'large_L_leading_coefficient': str(sp.factor(sp.Poly(H,L).LC())),
'fast_D_matrix': [[str(x) for x in reduced.row(i)] for i in range(3)],
'fast_D_cubic_coefficients': encoded(cubic),
'fast_D_hurwitz': str(sp.factor(cubic[0]*cubic[1]-cubic[2])),
'large_L_threshold': str(threshold), 'large_L_threshold_numeric': float(threshold)}
class Quartic:
@staticmethod
def hurwitz(coefficients):
a,b,c,d = coefficients
return a*b*c-c*c-a*a*d
def __init__(self, matrix):
if matrix.shape != (4,4): raise ValueError('Quartic check needs a four-species matrix.')
self.z = sp.Symbol('z')
self.p = matrix.charpoly(self.z).as_expr()
self.a = sp.Poly(self.p, self.z).all_coeffs()[1:]
def classification(self):
if not all(v > 0 for v in self.a):
return {'status': 'unknown', 'reason': 'Positive-coefficient quartic criterion not applicable.'}
H = sp.simplify(self.hurwitz(self.a))
status = 'strictly_stable' if H > 0 else ('unstable_complex_pair' if H < 0 else 'imaginary_pair_boundary')
return {'status': status, 'coefficients': encoded(self.a), 'H': str(H),
'second_hurwitz': str(sp.simplify(self.a[0]*self.a[1]-self.a[2]))}
def shifted_root_certificate(self, lo=2, hi=4):
lo, hi = Q(lo), Q(hi)
if not 0 < lo < hi: raise ValueError('Positive ordered shift bracket required.')
t = sp.Symbol('t')
a,b,c,d = sp.Poly(sp.expand(self.p.subs(self.z, self.z+t)), self.z).all_coeffs()[1:]
H = sp.expand(self.hurwitz((a,b,c,d)))
# Positive polynomial coefficients imply A(t), C(t)>0 throughout t>=lo.
positive = all(all(v >= 0 for v in sp.Poly(expr.subs(t,t+lo), t).all_coeffs()) and expr.subs(t,lo)>0 for expr in (a,c))
endpoints = [H.subs(t,lo), H.subs(t,hi)]
certified = positive and endpoints[0]*endpoints[1] < 0
return {'certified': bool(certified), 'real_part_bracket': encoded([lo,hi]),
'H_at_endpoints': encoded(endpoints), 'shifted_coefficients': encoded([a,b,c,d]),
'H_shift': str(H), 'root_construction': 'At an intervening zero of H(t), omega=sqrt(C(t)/A(t)); p(t+i*omega)=0.'}
class ChildLattice:
"""Enumerate injective consuming assignments, including the empty selection."""
def __init__(self, source):
self.source = source
def patterns(self, budget=100000):
choices = [[None]+[j for j in range(self.source.m) if self.source.Y[i,j]>0] for i in range(self.source.n)]
if np.prod([len(c) for c in choices], dtype=object) > budget:
raise RuntimeError('Enumeration budget exceeded; no absence-of-core conclusion.')
return [p for p in itertools.product(*choices) if len([x for x in p if x is not None]) == len({x for x in p if x is not None})]
def matrix(self, pattern):
selected = [i for i,r in enumerate(pattern) if r is not None]
if not selected: return sp.zeros(0,0)
return self.source.S.extract(selected, [pattern[i] for i in selected])
def padded(self, pattern):
if self.source.n != 4: raise ValueError('This quartic certificate route requires four species.')
M = -sp.eye(4)
for j,r in enumerate(pattern):
if r is not None:
for i,q in enumerate(pattern):
M[i,j] = self.source.S[i,r] if q is not None else 0
return M
def certificates(self):
scales = sp.symbols('a b c d', positive=True)
records = []
def nonnegative(expr, strict=False):
coefficients = sp.Poly(sp.expand(expr), *scales).coeffs()
return all(v >= 0 for v in coefficients) and (not strict or any(v > 0 for v in coefficients))
for pattern in self.patterns():
padded = self.padded(pattern)
a,b,c,d = (padded*sp.diag(*scales)).charpoly().all_coeffs()[1:]
H = sp.expand(Quartic.hurwitz((a,b,c,d)))
passed = nonnegative(a,True) and nonnegative(b) and nonnegative(c,True) and nonnegative(d) and nonnegative(H)
records.append({'assignment': list(pattern), 'species': [i for i,r in enumerate(pattern) if r is not None],
'padded_coefficients': encoded([a,b,c,d]), 'hurwitz_polynomial': str(H),
'hurwitz_monomials': [{'exponents': list(m), 'coefficient': str(v)} for m,v in sp.Poly(H,*scales).terms()],
'certified_D_nonunstable': bool(passed), 'singular': d == 0,
'scope': 'All positive diagonal column scalings, by polynomial nonnegativity and the quartic exclusion lemma.' if passed else 'Criterion inconclusive; no instability conclusion.'})
return records
def maximal_patterns(self):
patterns = self.patterns()
return [p for p in patterns if not any(p!=q and all(r is None or r==q[i] for i,r in enumerate(p)) for q in patterns)]
class DilutionCompletion:
"""Four active reactions, four adjusted feeds, and four first-order losses."""
def __init__(self, original, delta):
self.original = original; self.delta = Q(delta)
if self.delta <= 0 or original.equilibrium is None:
raise ValueError('Positive dilution and an exactly stationary original model required.')
s = original.source
if s.n != 4 or s.m != 6 or s.Y[:,4:] != sp.zeros(4,2):
raise ValueError('Completion expects this six-reaction source with its two terminal inflows.')
empty = (0,0,0,0); basis = [tuple(int(i==j) for i in range(4)) for j in range(4)]
reactions = list(s.reactions[:4])+[Reaction('feed '+s.species[i], empty, basis[i]) for i in range(4)]
reactions += [Reaction('loss '+s.species[i], basis[i], empty) for i in range(4)]
feeds = s.S[:,4:]*original.k[4:,:]+self.delta*original.equilibrium
self.feed_concentrations = feeds/self.delta
self.model = MassAction(Source(s.species,reactions), list(original.k[:4,:])+list(feeds)+[self.delta]*4, original.equilibrium)
def transport_certificates(self, original_records):
old = {tuple(r['assignment']): r['certified_D_nonunstable'] for r in original_records}
lattice = ChildLattice(self.model.source); rows = []
for p in lattice.patterns():
losses = [i for i,r in enumerate(p) if r is not None and r>=8]
oldpattern = tuple(None if i in losses else r for i,r in enumerate(p))
if any(p[i] != 8+i for i in losses): raise AssertionError('Non-diagonal loss assignment.')
retained = [i for i,r in enumerate(oldpattern) if r is not None]
# Verify the exact triangular block: loss columns vanish on every other species.
if any(self.model.source.S[j,p[i]] != -int(i==j) for i in losses for j in range(4)):
raise AssertionError('Loss transport block failed.')
if retained and self.model.source.S.extract(retained,[p[i] for i in retained]) != self.original.source.S.extract(retained,[p[i] for i in retained]):
raise AssertionError('Retained old-child block changed.')
rows.append({'assignment':list(p),'loss_species':losses,'old_assignment':list(oldpattern),
'certified_D_nonunstable':old.get(oldpattern,False)})
return rows
def local_departure(model, size=PERTURBATION_SIZE, end=TRAJECTORY_END):
if model.equilibrium is None or size <= 0 or end <= 0: raise ValueError('Positive perturbation, duration and reference equilibrium required.')
J = np.array(model.jacobian(),float)
eigenvalues, vectors = np.linalg.eig(J)
index = max(range(len(eigenvalues)),key=lambda i:(eigenvalues[i].real,eigenvalues[i].imag))
lam = eigenvalues[index]; w = vectors[:,index]
pivot = int(np.argmax(np.abs(w))); w = w/w[pivot] # deterministic phase and maximum modulus 1
equilibrium = np.array(model.equilibrium,float).reshape(-1)
times = np.linspace(0,end,TRAJECTORY_POINTS)
nonlinear = model.integrate(equilibrium+size*w.real,times)
linear = size*np.real(np.exp(lam*times[:,None])*w[None,:])
deviations = nonlinear-equilibrium
independent = model.integrate(equilibrium+size*w.real,times,method='BDF')
scale = float(np.max(np.linalg.norm(linear,axis=1)))
return times,deviations,linear,{'eigenvalue_real':float(lam.real),'eigenvalue_imag':float(lam.imag),
'relative_linear_error':float(np.max(np.linalg.norm(deviations-linear,axis=1))/scale),
'independent_solver_error_relative_to_mode':float(np.max(np.linalg.norm(nonlinear-independent,axis=1))/scale),
'scope':'Short-time numerical departure only; no periodic-orbit or long-time claim.'}
def write_csv(out,name,headers,rows):
with (out/name).open('w',newline='') as stream:
writer=csv.writer(stream);writer.writerow(headers);writer.writerows(rows)
def main():
parser=argparse.ArgumentParser(description=__doc__);parser.add_argument('--output',type=Path,default=Path('outputs'))
out=parser.parse_args().output;out.mkdir(parents=True,exist_ok=True)
family=OperatingFamily();model=family.model(CIRCULATION_RATIO,INVERSE_D_CONCENTRATION,THROUGHPUT)
children=ChildLattice(family.source);records=children.certificates()
if not all(r['certified_D_nonunstable'] for r in records): raise AssertionError('Original child exclusion not established.')
completion=DilutionCompletion(model,DILUTION_RATE);transport=completion.transport_certificates(records)
spectral=Quartic(model.jacobian());source_witness=family.model()
source_root=Quartic(source_witness.jacobian()).shifted_root_certificate()
times,deviations,linear,dynamics=local_departure(model)
cert=family.exact_family_certificate()
result={'inputs':{'T':CIRCULATION_RATIO,'L':INVERSE_D_CONCENTRATION,'s':THROUGHPUT,'dilution':DILUTION_RATE},
'molecularity':family.source.molecularity(),'rates':encoded(model.k),'equilibrium':encoded(model.equilibrium),
'flux':encoded(family.flux(CIRCULATION_RATIO,THROUGHPUT)),
'jacobian':[[str(v) for v in model.jacobian().row(i)] for i in range(4)],
'operating_classification':spectral.classification(),
'source_100_100_root_certificate':source_root,'family_certificate':cert,
'stable_control':Quartic(family.model(2,1).jacobian()).classification(),
'child_count':len(records),'singular_children':sum(r['singular'] for r in records),
'maximal_assignments':[list(p) for p in children.maximal_patterns()],
'completion':{'reaction_count':12,'feed_concentrations':encoded(completion.feed_concentrations),
'child_count':len(transport),'all_children_transported':all(r['certified_D_nonunstable'] for r in transport),
'spectral_classification':Quartic(completion.model.jacobian()).classification(),
'source_margin_applies':bool(Q(CIRCULATION_RATIO)==100 and Q(INVERSE_D_CONCENTRATION)==100 and Q(DILUTION_RATE)<2*Q(THROUGHPUT))},
'local_dynamics':dynamics,
'mass_obstruction':'2C->D and 2D->A+4B+4C imply mass_A+4*mass_B=0, impossible with positive masses.',
'scope':'Kinetic abstract-species counterexample; products can have nine molecules. No atom-balanced mechanism, Hopf theorem or Lean rerun.'}
(out/'results.json').write_text(json.dumps(result,indent=2)+'\n')
(out/'child_certificates.json').write_text(json.dumps(records,indent=2)+'\n')
(out/'completion_children.json').write_text(json.dumps(transport,indent=2)+'\n')
write_csv(out,'local_trajectory.csv',['time',*[f'nonlinear_d{x}' for x in family.source.species],*[f'linear_d{x}' for x in family.source.species]],
[[t,*d,*lin] for t,d,lin in zip(times,deviations,linear)])
perturbations=[]
for reaction in range(6):
for direction in (-1,1):
rates=model.k.copy();rates[reaction]*=1+direction*Q(RATE_PERTURBATION_FRACTION)
perturbed=family.from_rates(rates)
eig=np.linalg.eigvals(np.array(perturbed.jacobian(),float))
perturbations.append([reaction,direction*RATE_PERTURBATION_FRACTION,*[float(x) for x in perturbed.equilibrium],max(eig.real)])
write_csv(out,'rate_perturbations.csv',['reaction','relative_change','A_star','B_star','C_star','D_star','numeric_spectral_abscissa'],perturbations)
Ts=np.linspace(*PARAMETER_T_RANGE,141);Ls=np.geomspace(*PARAMETER_L_RANGE,111)
TT,LL=np.meshgrid(Ts,Ls);HH=Quartic.hurwitz(family.coefficients(TT,LL))
write_csv(out,'operating_grid.csv',['T','L','numeric_H'],zip(TT.flat,LL.flat,HH.flat))
lines=[f'Operating rates: {encoded(model.k)}; equilibrium: {encoded(model.equilibrium)}.',
f'Operating quartic: {spectral.classification()}.',
f'{len(records)} original children certified for all positive column scalings; {result["singular_children"]} singular.',
f'Source root bracket [2,4]: H endpoints {source_root["H_at_endpoints"]}.',
f'Stable control (T,L)=(2,1): H={result["stable_control"]["H"]}.',
f'Dilution completion: {len(transport)} transported children; feed={encoded(completion.feed_concentrations)}.',
f'Numerical local mode: {dynamics["eigenvalue_real"]:.9g} + {dynamics["eigenvalue_imag"]:.9g} i.',
f'Local nonlinear/linear relative error: {dynamics["relative_linear_error"]:.6g}.',
'Unstable linearization is not a sustained oscillation or atom-resolved mechanism.']
(out/'console.txt').write_text('\n'.join(lines)+'\n');print('\n'.join(lines))
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
fig,ax=plt.subplots(figsize=(8,4.8),layout='constrained')
ax.contourf(TT,LL,(HH<0).astype(float),levels=[-.5,.5,1.5],colors=['#edf1f3','#bdd8e7'])
ax.contour(TT,LL,HH,levels=[0],colors='#216b91',linewidths=1)
ax.axvline(cert['large_L_threshold_numeric'],color='#777',ls='--',label='Large-L threshold (not finite-L boundary)')
ax.plot([100,100],[100,200],color='#b95024',lw=3,label='Certified ray: T=100, L>=100')
ax.scatter([2],[1],color='black',s=30,label='Certified stable control')
ax.scatter([float(Q(CIRCULATION_RATIO))],[float(Q(INVERSE_D_CONCENTRATION))],marker='x',color='#b95024',label='Selected operating point')
ax.text(62,35,'H < 0: unstable\n(numerical grid)',ha='center',color='#184a66')
ax.text(62,.65,'H > 0: strictly stable\n(numerical grid)',ha='center')
ax.set(yscale='log',xlabel='Circulation ratio T',ylabel='Inverse D concentration L',title='Equilibrium stability versus circulation\nand concentration',xlim=PARAMETER_T_RANGE,ylim=PARAMETER_L_RANGE)
ax.legend(fontsize=8,loc='lower right')
fig.savefig(out/'operating_region.png',dpi=180);fig.savefig(out/'operating_region.svg');plt.close(fig)
fig,axes=plt.subplots(1,2,figsize=(10,4),layout='constrained')
axes[0].plot(times,deviations[:,0]/PERTURBATION_SIZE,label='Nonlinear mass action',color='#216b91')
axes[0].plot(times,linear[:,0]/PERTURBATION_SIZE,'--',label='Linear mode',color='#b95024')
axes[0].set(xlabel='Dimensionless time',ylabel='(A - A*) / perturbation size',title='Short-time local departure');axes[0].legend(fontsize=8)
axes[1].plot(deviations[:,0]/PERTURBATION_SIZE,deviations[:,2]/PERTURBATION_SIZE,color='#216b91')
axes[1].scatter([deviations[0,0]/PERTURBATION_SIZE],[deviations[0,2]/PERTURBATION_SIZE],color='#b95024',label='Initial state')
axes[1].set(xlabel='(A - A*) / perturbation size',ylabel='(C - C*) / perturbation size',title='Short-time concentration deviations in the\nA–C plane');axes[1].legend(fontsize=8)
for ax in axes:ax.grid(alpha=.2)
fig.savefig(out/'local_dynamics.png',dpi=180);fig.savefig(out/'local_dynamics.svg');plt.close(fig)
digest=lambda p:hashlib.sha256(p.read_bytes()).hexdigest()
(out/'run_metadata.json').write_text(json.dumps({'manuscript_sha256':MANUSCRIPT_SHA256,
'source_sha256':digest(Path(__file__)),'python':platform.python_version(),
'output_sha256':{p.name:digest(p) for p in sorted(out.iterdir()) if p.is_file() and p.name!='run_metadata.json'}},indent=2)+'\n')
if __name__=='__main__':main()
Run output
Operating rates: ['10000', '4', '100', '10000', '99', '2']; equilibrium: ['1', '1', '1', '1/100'].
Operating quartic: {'status': 'unstable_complex_pair', 'coefficients': ['10908', '185600', '1280000', '64000000'], 'H': '-5025252352000000', 'second_hurwitz': '2023244800'}.
25 original children certified for all positive column scalings; 4 singular.
Source root bracket [2,4]: H endpoints ['-2133097221521408', '2676722101092352'].
Stable control (T,L)=(2,1): H=626688.
Dilution completion: 121 transported children; feed=['100', '1', '1', '201/100'].
Numerical local mode: 2.99974676 + 15.6895257 i.
Local nonlinear/linear relative error: 4.92017e-06.
Unstable linearization is not a sustained oscillation or atom-resolved mechanism.