Verify a two-soliton collision

Can a computation reproduce a soliton collision, including its position shifts? This laboratory compares periodic numerical evolution with the full exact two-soliton field on a long interval. Separate refinements expose different errors, and endpoint measurements distinguish numerical inaccuracy from the finite separation of the pulses. The calculation tests this specified finite-time problem; it does not prove a general inverse-scattering theorem.

Required background. The collision lesson derives the interacting field and its signed shifts. The one-soliton numerical laboratory develops Fourier differentiation, dealiasing and integrating-factor RK4. Python and NumPy are needed to execute the experiment; its tables can be interpreted without running code.

Helpful background. The conservation lesson explains why preserved integrals alone do not establish a small solution error.

The exact field through the collision

Use the dimensionless, positive-field equation

$$u_t+6uu_x+u_{xxx}=0,\qquad x\in\mathbb R,$$

with smooth decaying fields and no background. For $\kappa_1\gt\kappa_2\gt0$, the target is

$$\begin{aligned} \tau&=1+e^{\eta_1}+e^{\eta_2}+A e^{\eta_1+\eta_2},\\ \eta_j&=2\kappa_j(x-4\kappa_j^2t-x_j),\\ A&=\left(\frac{\kappa_1-\kappa_2}{\kappa_1+\kappa_2}\right)^2, \qquad u=2\partial_x^2\log\tau. \end{aligned}$$

Here $x_j$ are phase parameters, not both incoming center intercepts. The Library derivation obtains this field from the two-bound-state Marchenko system, using Aktosun 2009, § IX, equations (9.1)–(9.3) with the opposite field sign translated.

The experiment fixes

$$\kappa_1=1,\qquad \kappa_2=\frac12,\qquad A=\frac19,\qquad x_j=\frac{\log A}{4\kappa_j}, \qquad -4\le t\le4.$$

These phases give $u(-x,-t)=u(x,t)$. The pulses are well separated at both endpoints and overlap near zero. See the three exact profiles and their endpoint markers; those curves evaluate the exact formula, independently of numerical time stepping.

Evaluate the field without overflow

Large positive $\eta_j$ make direct exponentiation unsafe even when $u$ remains small. Write the four logarithmic terms as

$$z=(0,\eta_1,\eta_2,\eta_1+\eta_2+\log A), \qquad w_s=\frac{e^{z_s-\max z}}{\sum_r e^{z_r-\max z}}.$$

Their $x$-slopes are $r=(0,2\kappa_1,2\kappa_2,2\kappa_1+2\kappa_2)$. Since $(\log\tau)_x=\sum_s w_sr_s$, differentiation gives the stable formula

$$u=2\sum_s w_s(r_s-\overline r)^2, \qquad \overline r=\sum_s w_sr_s.$$

This centered variance avoids subtracting nearly equal large moments. Exponent rescaling avoids overflow; sufficiently small tails can still underflow to zero in binary64.

The program tests this evaluator in several independent ways. An exact rational expansion makes all coefficients of the KdV bilinear identity vanish. Analytic derivatives give a sampled PDE residual below $4.5\times10^{-15}$. A direct two-by-two Marchenko matrix solve, at moderate positions and times for two unequal parameter pairs, agrees within $1.7\times10^{-14}$. These checks test signs, interaction coefficient and normalization before accepting the target. The computation notes define the samples and calculations.

Evolve a periodic approximation

The code samples the full exact field at $t=-4$, projects it onto retained Fourier modes, and evolves a periodic field on $[-\ell/2,\ell/2)$. This is a different boundary problem from the decaying line. Its usefulness depends on keeping boundary mismatch small throughout the selected time interval.

For $N$ equally spaced points, retain integer modes $|m|\lt N/3$. This strict two-thirds cutoff prevents aliases from the quadratic nonlinearity from entering the retained modes. With spatial wave numbers $k_m=2\pi m/\ell$ and Fourier coefficients $v_m$, evolve

$$\dot v_m=ik_m^3v_m-3ik_m\widehat{u^2}_m.$$

The existing one-soliton solver supplies the integrating-factor RK4 step. It applies classical RK4 to the transformed variable and resets the integrating factor at every step. This is the method described in Kassam and Trefethen 2005, pp. 1215–1216, equations (1.5)–(1.9), PDF. Independent stage, Fourier-convolution and linear-wave checks run again; the old one-soliton inputs and results are not changed.

Define the errors before comparing them

At each of the nine checkpoint times $t=-4,-3,\ldots,4$, compare every spatial sample with the exact field, without shifting the numerical profile to align its peaks. If $e_j=u_j^{\rm num}-u_j^{\rm exact}$, use

$$\varepsilon_2(t)= \frac{\sqrt{\Delta x\sum_j|e_j|^2}} {\sqrt{\Delta x\sum_j|u_j^{\rm exact}|^2}}, \qquad \varepsilon_\infty(t)= \frac{\max_j|e_j|}{\max_j|u_j^{\rm exact}|}.$$

The tables report the maximum over these nine checkpoints, not a bound at every intermediate time. The saved initial row also records the error introduced by projection.

The exact line integrals are

$$M=4\sum_j\kappa_j=6,\qquad P=\frac{16}{3}\sum_j\kappa_j^3=6,\qquad E=\frac{32}{5}\sum_j\kappa_j^5=6.6.$$

Their definitions are $M=\int u\,dx$, $P=\int u^2\,dx$ and $E=\int(u^3-u_x^2/2)\,dx$. Numerical drift is measured at every integration step, from the computed initial value, and divided by the corresponding exact line integral. Initial integral error is stored separately. This prevents a poorly resolved initial field from appearing accurate merely because its incorrect integral stays constant.

Refine time, space and interval separately

Time steps

Fix $\ell=80$ and $N=768$. The duration is $8$, so the first row uses $6{,}400$ steps.

Step $h$	Maximum $\varepsilon_2$	Maximum $\varepsilon_\infty$
$1/800$	$6.45\times10^{-7}$	$5.92\times10^{-7}$
$1/1600$	$2.68\times10^{-8}$	$2.46\times10^{-8}$
$1/3200$	$1.26\times10^{-9}$	$1.17\times10^{-9}$
$1/6400$	$6.08\times10^{-11}$	$5.76\times10^{-11}$

The observed orders $\log_2(\varepsilon_2(h)/\varepsilon_2(h/2))$ are $4.59$, $4.41$ and $4.37$. Report that finite sequence as measured: it is not an independent proof of an exact asymptotic order. The construction has classical fourth order for the fixed smooth semidiscrete problem; error constants, other discretization errors and floating-point effects affect a fitted order.

Spatial resolution

Fix $\ell=80$ and $h=1/1600$. The last column is the largest relative drift among $M,P,E$.

Points $N$	Maximum $\varepsilon_2$	Maximum $\varepsilon_\infty$	Largest invariant drift
128	$3.22\times10^{-1}$	$3.35\times10^{-1}$	$5.45\times10^{-11}$
192	$1.25\times10^{-2}$	$1.35\times10^{-2}$	$7.03\times10^{-10}$
256	$3.70\times10^{-4}$	$3.74\times10^{-4}$	$1.90\times10^{-9}$
384	$2.49\times10^{-6}$	$2.10\times10^{-6}$	$2.64\times10^{-9}$
512	$3.15\times10^{-8}$	$3.47\times10^{-8}$	$2.66\times10^{-9}$
768	$2.68\times10^{-8}$	$2.46\times10^{-8}$	$2.66\times10^{-9}$

The coarse run preserves its integrals very well while representing the field poorly. Increasing $N$ eventually exposes the time-error floor at this fixed step. The finer time study lowers that floor.

The artificial periodic interval

Fix $\Delta x=1/8$ and $h=1/1600$, increasing $N$ along with $\ell$. The final column is the largest absolute exact-field value at either boundary, sampled at $161$ times.

Length $\ell$	Points $N$	Maximum $\varepsilon_2$	Sampled boundary magnitude
40	320	$2.34\times10^{-3}$	$8.04\times10^{-3}$
48	384	$7.48\times10^{-7}$	$2.71\times10^{-6}$
64	512	$2.68\times10^{-8}$	$4.45\times10^{-12}$
80	640	$2.68\times10^{-8}$	$1.39\times10^{-15}$

The saved results also include endpoint value and derivative mismatches. The last two rows are limited by other errors; the table does not establish that the boundary error is exactly zero or bound it between sample times.

The final reference combines $\ell=80$, $N=768$ and $h=1/6400$. It has maximum checkpoint errors $6.08\times10^{-11}$ and $5.76\times10^{-11}$ in the two norms. Its largest relative invariant drift is $2.25\times10^{-12}$.

Separate two meanings of phase error

For an isolated endpoint maximum, locate the zero of $u_x$ near the predicted center. The code brackets each maximum within $1/\kappa_j$ of its asymptotic position and bisects; it uses the derivative of the Fourier interpolant for numerical data. It does not assign two peak trajectories during the overlap.

For endpoints $-T$ and $T$, a measured shift is

$$\widehat{\Delta x_j}(T) =x_j^{\rm peak}(T)-x_j^{\rm peak}(-T)-8\kappa_j^2T.$$

Even the exact field has finite-separation corrections to this estimate. At $T=4$:

Pulse	Asymptotic shift	Exact finite-time estimate	Numerical estimate
Fast	$1.098612288668$	$1.098611106194$	$1.098611106122$
Slow	$-2.197224577336$	$-2.197224576956$	$-2.197224576964$

The numerical errors relative to exact finite-time estimates are about $-7.18\times10^{-11}$ and $-8.55\times10^{-12}$. In contrast, the finite-time biases relative to asymptotic shifts are $-1.18\times10^{-6}$ and $3.80\times10^{-10}$. Refining the time step cannot remove those latter biases. The broad slow-pulse tail perturbs the fast maximum more strongly.

Using the exact formula alone, the fast-pulse bias decreases from approximately $-4.80\times10^{-4}$ at $T=2$ to $-1.18\times10^{-6}$ at $T=4$ and $-2.93\times10^{-9}$ at $T=6$. These are progressively better separated endpoints, not additional PDE solver runs.

Reject an attractive but incorrect approximation

Add two independent one-soliton fields $u_1+u_2$, using their correct incoming intercepts. Each term solves KdV, but their sum has residual

$$(u_1+u_2)_t+6(u_1+u_2)(u_1+u_2)_x+(u_1+u_2)_{xxx} =6\partial_x(u_1u_2).$$

The incoming sum starts close to the exact field: relative $L^2$ error is $3.73\times10^{-6}$ at $t=-4$. At $t=0$, its maximum sampled PDE residual is $3.21$; by $t=4$, its relative field error is $0.860$. It lacks the outgoing position shifts. Good incoming agreement therefore cannot justify replacing the interacting solution by a sum throughout the collision.

Reproduce the saved calculation

Download the complete collision experiment ZIP. It includes the shared solver in the correct sibling folder. Extract it and follow the root README; from integrable-kdv-collision/kdv-collision, the command is python3 experiment.py --check with the supplied NumPy dependency installed.

For inspection or individual downloads: collision experiment, inputs, saved results, requirements and notes, plus the shared one-soliton solver. Preserve sibling folders named kdv-collision and kdv; both programs retain their filename experiment.py.

With NumPy installed, a checkout uses:

python3 public/computations/kdv-collision/experiment.py --check

For separate downloads, run python3 computations/kdv-collision/experiment.py --check from above those folders. The notes give virtual-environment instructions. The recorded calculation uses Python 3.9.6, NumPy 2.0.2 and binary64 arithmetic; last digits can vary across platforms.

The check rejects changed inputs, reruns the independent checks and all refinements, and compares scientific outputs with the stated tolerances. It writes no files. Regenerating results requires the explicit --write option; routine site builds do not execute the solver.

Exercises

Guided: check the field at maximum overlap

For the stated parameters at $t=0$, show that $\tau=(1+e^x)^3$ and derive $u(x,0)$. Does that profile imply a single travelling soliton of inverse width $1/2$?

Hint

Here $e^{\eta_1}=3e^{2x}$ and $e^{\eta_2}=3e^x$. Differentiate $3\log(1+e^x)$ twice. Compare the result with the amplitude–width relation for one soliton.

Solution

The mixed term is $e^{3x}$, so $\tau=1+3e^x+3e^{2x}+e^{3x}=(1+e^x)^3$. Hence

$$u(x,0)=\frac{6e^x}{(1+e^x)^2} =\frac32\operatorname{sech}^2(x/2).$$

A single soliton with $\kappa=1/2$ would have height $2\kappa^2=1/2$, not $3/2$. This is an instantaneous collision profile. Its later evolution cannot be inferred by translating this shape rigidly.

Independent: diagnose the coarse calculation

A colleague accepts the $N=128$ run because its largest invariant drift is below $6\times10^{-11}$. Assess that decision using the recorded field errors. Which refinement addresses the dominant problem? Would reducing $h$ alone establish an accurate continuum solution?

Hint

Compare the first and last spatial rows, and distinguish conservation of the projected field from its agreement with the exact field.

Solution

The relative field errors reach $0.322$ and $0.335$, despite the small drift. The initial field and its evolution lack sufficient spatial resolution. Increasing $N$ at fixed $h$ reduces those errors to about $2.7\times10^{-8}$, where time error becomes important. Reducing $h$ on the coarse grid mainly converges to its inaccurate finite-dimensional evolution. A continuum claim also needs spatial and boundary checks, with the actual norm and finite-time scope stated.

Transfer: request a more accurate phase shift

You need the fast soliton’s asymptotic shift to absolute accuracy $10^{-8}$. Is the reference calculation at $T=4$ enough? Propose a changed experiment using the exact endpoint-time evidence, and explain why simply extending the time on an unchanged small periodic box can fail.

Hint

Compare the fast numerical error with its finite-time bias. The fast pulse travels at speed $4$.

Solution

No. Its numerical error is about $7.2\times10^{-11}$, but the endpoint bias is about $1.18\times10^{-6}$. The exact $T=6$ bias is about $2.93\times10^{-9}$, so those endpoints are a promising starting point. Recheck full-field errors, peak errors and the boundary mismatch for a numerical run over $[-6,6]$; the existing run does not certify that new experiment.

At $t=6$, the fast asymptotic center is near $24.55$. A box of length $48$ ends at $24$, so this pulse has reached the periodic seam. Increase the interval as needed while retaining spatial resolution, and refine the time step separately. A smaller time step alone removes neither finite-separation bias nor the wrong boundary problem.

References

• Aktosun, Tuncay. “Inverse Scattering Transform and the Theory of Solitons.” In Encyclopedia of Complexity and Systems Science, Springer, 2009, pp. 4960–4971. DOI. Author version, arXiv:0905.4746v1, § IX, equations (9.1)–(9.3).
• Kassam, Aly-Khan, and Lloyd N. Trefethen. “Fourth-order time-stepping for stiff PDEs.” SIAM Journal on Scientific Computing 26(4), 2005, pp. 1214–1233. DOI. Author PDF, pp. 1215–1216, equations (1.5)–(1.9).