# KdV pulse: independent error checks

This experiment solves a finite Fourier approximation to periodic
`u_t + 6 u u_x + u_xxx = 0` and compares it with the exact decaying-line pulse
`u = 2 kappa^2 sech^2(kappa (x - 4 kappa^2 t - x0))`.
The finite torus and real line are different boundary-value problems. The
separate domain study measures this distinction over the stated short interval.

## Run

Place `experiment.py`, `inputs.json`, and `results.json` together. Python 3.9–3.12
with NumPy 2.0.2 reproduces the recorded environment (Python 3.9.6, NumPy 2.0.2).
From this directory, after installing the supplied requirements:

```sh
python3 experiment.py --check
```

This is nonmutating. No flags prints fresh JSON; `--write` explicitly replaces
`results.json` after all scientific checks pass. Environment version strings are
reported but not compared. Numeric comparisons use the tolerances in inputs;
they allow platform-dependent roundoff rather than demanding identical bytes.
Input values and the shape of the saved result are also checked.

## Spatial discretization

Use an even number N of points on `[-L/2,L/2)`, with dx=L/N and Fourier wave
numbers k=2*pi*m/L. NumPy's FFT coefficients are unnormalized; its inverse FFT
includes 1/N. Retain only integer modes `abs(m) < N/3`, both initially and in
each nonlinear evaluation. Thus the grid has N collocation points but fewer
independent retained modes (341 for N=512). The Nyquist mode is excluded.

For a retained state v, the Fourier ODE is
`v_t = i*k^3*v - 3*i*k*FFT((IFFT(v).real)^2)` followed by the same projection.
This strict 2/3 rule prevents quadratic wraparound terms from entering retained
modes. The program compares this product with an explicit, non-circular Fourier
convolution on an independent random test. It is a Galerkin truncation; it does
not restore physical scales absent from the retained modes.

## Time stepping

Write v'=D v+G(v), D=i*k^3. On each step transform
`w(s)=exp(-s*D)*v(t+s)` and apply classical RK4 to
`w'=exp(-s*D)*G(exp(s*D)*w)` for 0<=s<=h. Reset the local clock each step.
The function `ifrk4_step` is the expanded update used by the solver. A separate
test evaluates ordinary RK4 in the transformed variable for a complex nonlinear
test system and compares the updates. With G=0, an independent Airy-wave check
verifies the sign of the exactly integrated linear phase.

This is integrating-factor RK4, not ETDRK4. No singular phi-function coefficients
or contour approximation is used. Treating dispersion exactly does not make
arbitrary time steps accurate or remove nonlinear stability limits.

Method source: Aly-Khan Kassam and Lloyd N. Trefethen, “Fourth-order time-stepping
for stiff PDEs,” *SIAM Journal on Scientific Computing* 26(4), 1214–1233 (2005),
[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), explain the IF transform and RK4.
The implementation, refinement choices and independent checks here are original.

## Error and conservation definitions

At T=1 compare the computed samples with the exact line pulse on the same grid:

- Relative L2 = sqrt(dx*sum(abs(error)^2))/sqrt(dx*sum(abs(exact)^2)).
- Relative Linf = max(abs(error))/max(abs(exact)).
- Discrete integrals use periodic trapezoidal quadrature; the derivative in the
  third integral is evaluated by Fourier differentiation.

The recorded integrals are M=int(u), P=int(u^2), E=int(u^3-u_x^2/2).
Their line-pulse values are 4*kappa, 16*kappa^3/3, 32*kappa^5/5. The model page
uses translation generator P/2 and Hamiltonian H=-E; these are normalization choices.
Each reported drift is the maximum over every time step of
abs(I(t)-I(0))/abs(I_exact). Initial error against the exact line integral is
recorded separately. These are sampled-time maxima, not a continuous-time bound.

Three studies vary one numerical control at a time:

- Time: L=48, N=512; dt=1/160, 1/320, 1/640, 1/1280.
- Space: L=48, dt=1/3200; N=64, 96, 128, 192, 256, 384.
- Domain: dx=1/8, dt=1/3200; L=16, 24, 32, 48 (N grows with L).

The reference run is L=48, N=512, dt=1/3200. Temporal convergence approaches
fourth order for these smooth data. Spatial and domain refinement eventually
reach the remaining error floor; a plateau is not evidence of failed refinement.
The shortest box shows finite-domain error even though the pulse does not reach
an endpoint. Boundary values and derivative jumps are recorded at t=0 and T.
For these initial/final centers (−2 and +2), the largest endpoint magnitude over
0<=t<=1 occurs at one of those two times.

## Independent analytic checks and negative control

At three distinct kappa values the program compares the travelling wave with
the rank-one reconstruction, checks twice the total diagonal derivative of K
using a fourth-order centered difference, and checks the Marchenko integral by
384-point Gauss–Legendre quadrature over a finite interval. Its omitted
exponential tail is bounded by the factor exp(−48). A separate quadrature over
20 pulse widths in each direction verifies M, P, E; their omitted tails are
exponentially small. These floating-point checks supplement the on-site exact
derivations; they are not general inverse-scattering proofs.

The wrong-speed control translates the correct shape at speed 2*kappa^2 rather
than 4*kappa^2. It preserves all three line integrals exactly but fails the KdV
equation and has relative final L2 error about 1.23. This is a deliberate failed
approximation, not a second evolution solver.

The checks assert method identities, resolved reference accuracy, strong
independent refinement and a distinguished negative control. Saved data carry
the full tested parameters, measured norms, invariant drifts and representative
exact samples. No solver runs during routine site builds.

## Reconstruction source and figure

Tuncay Aktosun, “Inverse Scattering Transform and the Theory of Solitons,”
*Encyclopedia of Complexity and Systems Science*, 4960–4971 (2009),
[DOI](https://doi.org/10.1007/978-0-387-30440-3_295),
[arXiv v1](https://arxiv.org/html/0905.4746v1), §§ VII–IX, especially (7.9) and
(8.1)–(8.3). His Schrödinger potential is q=−u in this site's positive-pulse
convention. The displayed reconstruction and its signs are derived in the
Library article; the numerical quadratures do not rely on copied source code.

The spatial-convergence figure is generated from `results.json` by
`figures-src/kdv/plot_convergence.py` from the repository root, using Matplotlib
3.9.4 in addition to the NumPy solver environment. Its black solid circles show
relative trajectory error; dashed gray squares show the maximum normalized
drift of P and E. Both curves use the recorded spatial study without smoothing.
The accompanying page table provides the semantic data alternative.