# Reproduce a two-soliton KdV collision

This finite-time experiment supports `/learn/lab-kdv-collision/`. It evolves
the exact two-soliton initial field through a collision and compares the entire
computed field with the exact answer. Separate refinements test time, space and
the artificial periodic interval. Endpoint peak measurements distinguish
numerical error from the finite-separation bias in an asymptotic phase shift.

## Files and execution

Use Python 3.9–3.12 and the pinned NumPy dependency. The recorded environment is
Python 3.9.6, NumPy 2.0.2, IEEE 754 binary64. From a site checkout:

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

With NumPy installed, `python3 public/computations/kdv-collision/experiment.py --check`
is sufficient. On Windows use the virtual environment's `Scripts/python.exe`.

The complete download at `/computations/kdv-collision.zip` preserves the required
layout, including the shared solver. Extract it, enter
`integrable-kdv-collision/kdv-collision`, install `requirements.txt` in a virtual
environment, and run `python3 experiment.py --check` with that environment's
Python. The archive's root README gives the exact commands.

For individual downloads, preserve this small directory layout:

```text
computations/
  kdv/
    experiment.py          # companion from /computations/kdv/experiment.py
  kdv-collision/
    experiment.py
    inputs.json
    results.json
    requirements.txt
    README.md
```

The companion is the existing one-soliton program, loaded as a library. This
experiment reuses its Fourier grid, nonlinear term, integrating-factor RK4 step,
quadratures, error norms and independent method checks. It does not execute that
program's main routine, read its inputs/results, or change its baseline. Direct
source loading avoids writing Python bytecode caches. No renamed solver copy is
needed. The ZIP and all six individual files are linked from the lab.

`--check` reruns the experiment, validates mathematical and numerical checks,
then compares results without changing any artifact. No option prints fresh
JSON. `--write` deliberately regenerates `results.json`; inspect changed inputs,
results, teaching tables and figure together. Ordinary site builds do not run
this solver. The saved figures and numerical tables have different provenance:
the profile figure is an exact formula evaluation, not output from this solver.

## Exact target and its independent checks

On the real line use `u_t + 6*u*u_x + u_xxx = 0`, with zero background. For
`k1>k2>0`, set

```text
A = ((k1-k2)/(k1+k2))^2
eta_j = 2*k_j*(x - 4*k_j^2*t - x_j)
tau = 1 + exp(eta1) + exp(eta2) + A*exp(eta1+eta2)
u = 2*partial_x^2 log(tau).
```

The inputs take `k=(1,0.5)`, `A=1/9`, and `x_j=log(A)/(4*k_j)`. These phase
parameters enforce `u(-x,-t)=u(x,t)` and are not both incoming intercepts.
The integration interval is `[-4,4]`. At time zero, independently expanding
`tau=(1+exp(x))^3` gives `u=(3/2)*sech(x/2)^2`.

Numerical evaluation uses four log weights, subtracting their maximum before
exponentiation. If `w_s` are the normalized positive weights and
`r=(0,2*k1,2*k2,2*(k1+k2))`, then `u=2*sum(w_s*(r_s-mean(r))^2)`.
The centered variance avoids cancellation in a difference of two large moments;
the exponential rescaling avoids overflow. Extremely small tails may underflow
to zero in binary64. They are not assigned a nonzero numerical resolution.

Analytic derivatives use centered moments: `u_x=2*mu3`,
`u_xxx=2*(mu5-10*mu3*mu2)`, and
`u_t=2*E[(r-E[r])^2*(s-E[s])]`, with
`s=(0,-8*k1^3,-8*k2^3,-8*(k1^3+k2^3))`.
The program checks their PDE residual at 1,001 positions and nine times. A
separate exact calculation uses Python `Fraction` to expand every coefficient
of the bilinear polynomial `(D_x D_t + D_x^4) tau . tau`: all nine exponent
classes vanish for the prescribed rational kappas. Setting `A=1` fails this test.

The moderate-exponent direct quotient `2*(tau_xx/tau-(tau_x/tau)^2)` is compared
with the variance. Independently, for `(k1,k2)=(1,0.5),(1.2,0.7)` and moderate
`x,t`, construct `G_jl=delta_jl+w_j*w_l/(k_j+k_l)`,
`w_j=sqrt(c_j(t))*exp(-k_j*x)`,
`c_j(0)=2*k_j*exp(2*k_j*x_j)/A`, `c_j(t)=c_j(0)*exp(8*k_j^3*t)`.
Solving `G*z=w` gives `S=w.z` and `u=4*(k*w).z-2*S^2`. This independent
Marchenko reconstruction is compared with the stable tau evaluator; it is not
used far to the left where this direct matrix evaluation loses conditioning.

640-point Gauss–Legendre quadrature on `[-48,48]` at times `-4,0,4` checks
the line integrals `M=4*sum(k)`, `P=16*sum(k^3)/3`, `E=32*sum(k^5)/5`.
Their exact values are `(6,6,6.6)`. Quadrature errors include finite tails and
floating-point evaluation; these are numerical checks, not exact equalities.

## Numerical boundary problem and norms

The solver is periodic on `[-length/2,length/2)`. It starts from Fourier-projected
samples of the full tau field at `t=-4`, not from two separately moving pulses.
The continuous line solution does not satisfy exact periodic boundary conditions.
At 161 uniformly spaced times the report records both endpoint magnitudes, their
value mismatch, and their derivative mismatch, as maxima over those samples.
These sampled diagnostics are not a rigorous bound between sample times.

Use `N` equispaced points and retain the Fourier modes with integer `|m|<N/3`.
This strict two-thirds cutoff prevents aliasing into the retained quadratic
nonlinearity. It also resolves the cubic quadrature used for the projected
field's energy: three retained wave numbers cannot sum to a nonzero multiple
of `N`. The Fourier equation is `v_t=i*k^3*v -3*i*k*FFT(u^2)`.
The companion applies classical RK4 to the local integrating-factor variable;
this is not ETDRK4. Independent transformed-stage, explicit-convolution and
linear-Airy tests are rerun before accepting results.

At each of the **nine specified checkpoint times**, errors compare the whole
spatial field, without shifting or aligning it:

```text
relative L2   = sqrt(dx*sum(abs(u_num-u_exact)^2)) / sqrt(dx*sum(abs(u_exact)^2))
relative Linf = max(abs(u_num-u_exact)) / max(abs(u_exact)).
```

The displayed maximum is over those nine times, not all intermediate integration
steps. Initial projection error is the first row. Invariant drift is checked at
every integration step and is `max_n abs(I_num(n)-I_num(0))/abs(I_line)` for each
of `M,P,E`. Initial error relative to exact line integrals is stored separately.

- Time refinement fixes length 80 and N=768 and halves the step three times.
- Space refinement fixes length 80 and 12,800 steps while changing N.
- Domain refinement fixes spacing 1/8 and 12,800 steps while changing length and N.
- The reference uses length 80, N=768 and 51,200 steps. Its smaller step is stated
  separately; the space/domain tables share a coarser time-error floor.

The observed time orders describe the saved finite sequence. They are not an
independent proof of the order-four method or a rigorous continuum error bound.
Once other errors dominate, further refinement need not improve the result.

## Phase measurements and an incorrect control

For each kappa the incoming/outgoing intercepts are
`(x1,x2-log(A)/(2*k2))` and `(x1-log(A)/(2*k1),x2)`. The asymptotic shifts are
`(+log(3),-2*log(3))`. The matrix-free analytic derivative locates an isolated
maximum near each predicted center using a sign-changing bracket of radius
`1/kappa` and 60 bisection steps. This is done only at separated endpoint times.
No two pulse identities or maxima are followed through the overlap.

For the numerical reference, the derivative of the Fourier interpolant locates
endpoint maxima in the same brackets. Compare each center first with the exact
finite-time maximum. A shift estimate is
`center_out-center_in-4*kappa^2*(t_out-t_in)`; its numerical error and its
finite-time bias from the analytic asymptotic shift are reported separately.
An independent exact endpoint-time sequence `T=2,3,4,6` shows the latter bias
shrinking. The slow pulse's broad tail perturbs the fast peak more strongly.

The negative control adds two independent one-soliton fields with the correct
incoming intercepts. Its PDE residual is `6*partial_x(u1*u2)`. It is initially
close when the pulses are separated but has the wrong collision and no outgoing
phase shifts. Saved full-field errors and residuals expose this failure.

## Saved comparison and figure

`--check` rejects any changed inputs before running, requires finite outputs and
all acceptance checks, and compares every scientific output field with relative
tolerance 2e-5 and absolute tolerance 3e-10. Integer counts and metadata are exact;
runtime environment version strings may differ. The absolute tolerance permits
roundoff in near-zero checks; it is not the experiment's accuracy. No figure is
rewritten by this command. Comparisons deliberately do not use whole-file hashes
as a scientific criterion. Small differences within the stated numerical
tolerances are accepted; missing fields or material numerical changes fail.

The original plot source is `figures-src/kdv/collision-profiles.py`, with output
`public/figures/kdv/collision-profiles.svg`. It shows exact profiles at `t=-4,0,4`
on shared linear scales; dashed endpoint lines are asymptotic predictions.
The corresponding lesson supplies a semantic table and a caption explaining
the finite-time distinction. To regenerate in a checkout with Matplotlib 3.9.4:

```sh
python3 figures-src/kdv/collision-profiles.py
python3 figures-src/kdv/collision-profiles.py --check
```

The second command compares an in-memory regenerated SVG without changing files.
The figure check verifies the independent time-zero profile and signed shifts;
it does not rerun numerical evolution. Plot serialization can differ with other
Matplotlib versions without changing the mathematics.

## Sources

- Tuncay Aktosun, “Inverse Scattering Transform and the Theory of Solitons,” in
  *Encyclopedia of Complexity and Systems Science*, Springer, 2009, pp. 4960–4971.
  DOI: https://doi.org/10.1007/978-0-387-30440-3_295 . Author version
  https://arxiv.org/html/0905.4746v1 , § IX, equations (9.1)–(9.3). The source potential
  is the negative of this positive KdV field. Our Library derivation reduces the
  two-by-two determinant to tau and fixes the growing-exponential phase mapping.
- Aly-Khan Kassam and Lloyd N. Trefethen, “Fourth-order time-stepping for stiff
  PDEs,” *SIAM Journal on Scientific Computing* 26(4), 2005, pp. 1214–1233.
  DOI: https://doi.org/10.1137/S1064827502410633 . Author PDF:
  https://people.maths.ox.ac.uk/trefethen/publication/PDF/2005_111.pdf .
  Printed pp. 1215–1216, equations (1.5)–(1.9) support the integrating-factor RK4
  construction; their later ETDRK4 construction is a different method.

Reported measurements are original outputs of this experiment, not numbers
quoted from either source. Exact algebra establishes this two-soliton family;
finite numerical agreement does not establish a general inverse-scattering
theorem, completeness, or a long-time error bound.