# Three-site open Toda experiment

This small experiment supports `/learn/project-three-site-toda/`. It compares
velocity Verlet with an independently evaluated QR solution and, for the symmetric
initial data, an elementary exact solution. It is a finite-time numerical check,
not a proof of Liouville integrability or a benchmark for very long times.

## Reproduce

Use Python 3.9–3.12 for the pinned NumPy dependency, as documented in the
[NumPy 2.0.2 release notes](https://github.com/numpy/numpy/releases/tag/v2.0.2).
The recorded run used Python 3.9.6,
NumPy 2.0.2, and IEEE 754 binary64 arithmetic. From the repository root:

```sh
python3 -m venv .tmp/toda-venv
.tmp/toda-venv/bin/python -m pip install -r public/computations/toda/requirements.txt
.tmp/toda-venv/bin/python public/computations/toda/experiment.py --check
```

If NumPy is already available, `python3 public/computations/toda/experiment.py --check`
is sufficient. Downloaded copies of `experiment.py`, `inputs.json`, and
`results.json` also work together in any directory; paths are relative to the script.
The command recomputes the experiment, checks the equations and convergence, and
compares saved results using relative tolerance 1e-6 and absolute tolerance 2e-9.
It writes no artifacts. Platform-dependent last digits and near-zero residuals are
allowed; whole-file hashes do not decide scientific agreement.

With no option, the script prints freshly computed results without comparing the
saved file. `--write` explicitly replaces `results.json`; review changed inputs,
results, article numbers, and the figure together. No solver runs during an ordinary
site build. Retain the original inputs if performing an extension in a copied folder.

## Problem and two independent evaluations

All quantities are dimensionless. There are three real positions and momenta, with
canonical Poisson brackets, unit masses, and **two** exponential bonds:

`H = (p1^2+p2^2+p3^2)/2 + exp(q1-q2) + exp(q2-q3)`.

The force is `(-e1, e1-e2, e2)` with `e1=exp(q1-q2)` and
`e2=exp(q2-q3)`. There is no periodic bond or endpoint spring. The Lax matrix has
diagonal `p` and off-diagonal `a_i=exp((q_i-q_(i+1))/2)`. Its skew matrix has
upper entries `-a_i/2` and lower entries `a_i/2`; the evolution is `dL/dt=[B,L]`.
The direct scalar invariants checked against matrix traces are
`P=sum(p)`, `H=tr(L^2)/2`, and
`I3=sum(p^3)/3 + (p1+p2)*e1 + (p2+p3)*e2`.

- Velocity Verlet applies half a force kick, a full position drift, and a second
  half kick. All grid points through `T=6` enter maximum-error measurements.
- Independently factor `exp(-t*L0/2)=Q*R` with positive diagonal of `R`, and set
  `L(t)=Q.T*L0*Q`. Diagonalizing the initial symmetric matrix evaluates the matrix
  exponential; no time integration is used. Reconstruct position differences as
  `q_i-q_(i+1)=2*log(a_i)`, and fix the common position shift by
  `mean(q(t))=mean(q0)+t*P/3`. The matrix exponential is scaled by a positive scalar
  to avoid overflow; this does not change `Q`. Severe conditioning at long times
  is not cured by this scaling and is outside the verified interval.
- For `q0=(0,0,0), p0=(1,0,-1)`, also compare to the closed form
  `q=(x,0,-x), p=(v,0,-v)`, where
  `s=sqrt(3)*t/2-atanh(1/sqrt(3))`,
  `x=log(3/2)-2*log(cosh(s))`, `v=-sqrt(3)*tanh(s)`.
  This reference uses neither the force implementation nor the Lax matrix.

The second initial condition is deliberately nonsymmetric:
`q0=(0.2,-0.3,0.1), p0=(0.7,-0.4,0.2)`. Its cubic invariant is nonzero.
For the first case, preserved reflection symmetry forces `P=I3=0` even in a
trajectory with substantial integration error; their apparent conservation alone
would be a weak test. In the second case both bonds and all three momenta evolve.
The `symmetric` case name is tied to the stated closed-form initial data: when
changing that input in an extension, adapt or remove its analytic oracle as well.
The separate `shift_boost_check` function tests the exact transformation
`q'=q-0.2+0.4*t`, `p'=p+0.4` against the transformed closed-form solution. Its
invariant and trajectory residuals are also retained in the result file.

## Evidence and limitations

`results.json` preserves input values, selected trajectories (every tenth step of
the finest run), maximum absolute invariant/eigenvalue drifts, trajectory errors,
observed refinement orders, and independent equation checks. Trajectory error is
the maximum of all six absolute component errors over **every** time step of each
run. Invariant drift is measured from the initial scalar value; spectral drift
uses ascending sorted eigenvalues. No relative drift is defined for zero invariants.

The algebra check uses nonuniform bonds and compares analytic `dL/dt` with `[B,L]`,
scalar invariants with traces, and force with a centered energy derivative of step
1e-5. The QR ODE check uses fourth-order centered derivatives at times
`0,0.5,1.25,3,6`, at derivative spacings `0.01` and `0.005`. Each derivative samples
up to two spacings beyond its central time; these small endpoint extensions are
included only in the differential residual test. The fine-spacing residual is
required below 1e-7. The exact symmetric reference checks QR to 1e-9. Step halving
must give observed order in `(1.95,2.05)`, and the finest trajectory error must be
below 1e-4 for both prescribed cases. These are regression tolerances, not rigorous
error bounds. The saved run resolves the primary trajectory error at about 1e-5,
while QR versus the closed form differs by less than 3e-12.

## Figure brief and regeneration

- Reader question: does smaller time step improve the full trajectory at the
  expected rate, including when symmetry is broken?
- Takeaway: halving the Verlet step reduces maximum state error by about four in
  both fixed finite-time experiments.
- Objects: two log–log series of maximum `q,p` error versus step size. Filled
  circles/solid line identify symmetric data; hollow squares/dashed line identify
  asymmetric data. Axes are dimensionless; data are numerical, not schematic.
- Caption/alt: both cases show second-order trajectory convergence on `0<=t<=6`.
  The article provides the numerical values and explains the norm next to the plot.
- Provenance: original computed data in `results.json`, no third-party figure.
- Independent check: every plotted number is an all-grid maximum against QR;
  QR is checked against Hamilton's equations and the symmetric closed form. The
  expected factor of four follows from the independent second-order method result.

The editable source is `figures-src/toda/plot_experiment.py`; its SVG output is
`public/figures/toda/three-site-convergence.svg`. Install the optional plotting
dependency in the same environment and run:

```sh
.tmp/toda-venv/bin/python -m pip install matplotlib==3.9.4
.tmp/toda-venv/bin/python figures-src/toda/plot_experiment.py
.tmp/toda-venv/bin/python figures-src/toda/plot_experiment.py --check
```

The figure uses a light canvas and redundant grayscale line/marker encodings.
The renderer is headless and omits timestamp metadata. Plot verification compares
the regenerated SVG with the saved figure; it is a source-to-output check, not a
new scientific certification. Different Matplotlib versions may change the SVG
serialization without changing the calculation.

## Sources and convention mapping

- Anthony M. Bloch and Steven N. Karp, *Symmetric Toda, gradient flows, and
  tridiagonalization*, arXiv:2304.10697v1 (2023), pp. 1–2.
  https://arxiv.org/abs/2304.10697v1 . Their general symmetric matrix equation
  `dot M=[M,pi_u(M)]` has solution obtained from `exp(t*M0)`. Substituting `M=-L/2`
  gives this experiment's `exp(-t*L0/2)` and `dot L=[B,L]`. This uses the general
  symmetric-matrix flow, so the negative off-diagonal entries of `M` are allowed.
  The project also derives the sign directly by differentiating QR.
- Ernst Hairer, *Geometric Numerical Integration, Lecture 2: Symplectic
  integrators*, TU München, January–February 2010, §1, Theorem 2, p. 2.
  https://www.unige.ch/~hairer/poly_geoint/week2.pdf . The separable-Hamiltonian
  specialization gives the velocity Verlet steps used here and their order two.

These sources support method/convention choices. Reported numerical values are
recomputed by the supplied original code; they are not quoted from either source.