Numerics and validation of KAM tori

Home / KAM course

Working directory: run this guide’s commands from toWebSite/codes/KAMExample/, unless stated otherwise. Setup and complete workflow · Original Markdown.

A standalone KAM example

Build and run with a C++17 compiler and GNU make:

make
./exampleKAM.x
./exampleKAM.x --help
make check

No TorKam, FFTW, GMP, MPFR, or other numerical library is required. The only non-standard C++ interface is POSIX/GNU getopt_long for command-line parsing (available on the course’s Linux environment).

The numerical routines now live in the local common library ../lib/kam.h and ../lib/kam.cpp, shared with ../continuation. The Makefile compiles that source automatically. Keep codes/lib/ alongside this directory when copying the example; it is self-contained within the course repository.

Hamiltonian and notation

The two-degree-of-freedom example is

\[ H_\varepsilon(\varphi,I)=\tfrac12(I_1^2+I_2^2) +\varepsilon[\cos\varphi_1+\cos\varphi_2+\cos(\varphi_1-\varphi_2)]. \]

At zero coupling, I=omega is an invariant torus. The mixed cosine couples the two rotators. We prescribe omega=(1,(1+sqrt(5))/2) by default and find

\[ K(\theta)=(\theta+u(\theta),\omega+v(\theta)),\qquad X_H(K)-DK\,\omega=0. \]

Both u and v are periodic; the nonperiodic lift theta is never sent through the FFT. Angles have period 2 pi, and the ordering is (phi1,phi2,I1,I2). The initial guess is u=v=0.

Connection to TorKam and the slides

The geometric step follows the separately BSD-licensed src_UTILS/KAMLagrangian.cpp in /home/figueras/work/UTILS/lib_TorKam/torkam-dev, especially compute_N_LAG, compute_invP_LAG, compute_torsion_LAG, compute_etaK, solve_xi_of_K_LAG, and KAM_LAG_step. The standalone implementation replaces its multiprecision grids and linear tools with std::vector, std::complex, a radix-2 FFT and checked 2-by-2 matrix operations. The common source credits the original routines and retains their copyright/license notice.

Each correction follows the slide sequence:

  1. Compute L=DK and E=X_H(K)-L*omega using spectral differentiation.
  2. Construct N=-J*L*(L^T*L)^(-1) and store P=(L,N) as its two blocks.
  3. Use Q=-J*P^T*J as the approximate inverse and compute eta=-Q*E=(N^T*J*E,-L^T*J*E).
  4. Compute T=-N^T*J*(Lie_omega N+DX_H(K)*N).
  5. Solve the nonzero Fourier modes of Lie_omega xiN=etaN, where Lie_omega=-omega.dot(partial_theta) and the divisor is -i*k.dot(omega).
  6. Solve <T>*c=<etaL-T*tildeXiN> to fix <xiN>=c.
  7. Solve Lie_omega xiL=etaL-T*xiN, setting <xiL>=0 to fix the phase.
  8. Update the periodic part by Delta=L*xiL+N*xiN and recompute the residual.

The normal equation’s average vanishes geometrically; the program prints its numerical size before projecting it out. Singular/ill-conditioned Gram or average torsion matrices and numerically unresolved divisors cause failure. Custom frequencies are checked only on the retained modes; this finite check is not a Diophantine proof.

Fourier resolution and output

-n N retains modes with |k1|,|k2| < N/2; Nyquist lines are excluded to avoid ambiguities in real-valued interpolation. Nonlinear operations are evaluated on a 2N by 2N zero-padded grid. This avoids aliasing into the retained band for quadratic products of retained series and reduces aliasing for the nonpolynomial operations. The latter are checked through the final residual.

Convergence requires the unprojected residual to be below --tolerance on the work grid and a 4N by 4N check grid. The program also reports the energy oscillation and the largest embedding coefficient in the outer quarter of the Fourier band. These are floating-point numerical checks, not rigorous bounds or an existence validation. A stalled iteration may need a larger -n, a smaller perturbation, or a tolerance compatible with double precision.

./exampleKAM.x -e 0.01 -n 32 -t 1e-11 -o torus.csv
./exampleKAM.x -e 0.03 -n 64 -t 1e-10 -o torus-finer.csv

The CSV contains the base grid, the lifted embedding, its periodic parts, the four residual components, and the energy. Its first comment line records epsilon and omega. phi1 and phi2 are unwrapped; reduce modulo 2*pi for plots if desired. Output is written only on convergence. Exit status is zero on success, one on numerical/I/O failure, and two for invalid options. make clean removes the executable and preserves computed CSV files.

make check runs analytic checks of Fourier differentiation, the cohomological inverse, dealiased multiplication and Hamiltonian derivatives, followed by integrable and coupled end-to-end computations. The latter use /dev/null instead of creating output files.

Rigorous existence validation

The separate, self-contained C++ program in ../validation/ checks the modified Figueras–Haro a posteriori theorem for this Hamiltonian and the exact golden-ratio frequency. It requires MPFI/MPFR/GMP and does not link TorKam or the numerical common library.

From this directory, for example:

make -C ../validation
../validation/validateKAM.x --input results/torus-eps0p01.csv \
  --epsilon 1/100 --precision 256 --refine 2 --certificate validation.json

The optional refinement improves the seed at multiprecision; the validator then independently checks a frozen Fourier polynomial with interval arithmetic and analytic tail bounds. The exact frequency and action offset are specified algebraically, rather than treating the rounded CSV frequency as Diophantine. See the validator’s README for a successful nonzero-coupling certificate, rechecking instructions, precise input semantics, and scope.

Reproduce the numerical slides

The optional plotting stage requires Python 3, NumPy and Matplotlib:

python3 -m pip install numpy matplotlib
make figures

plot_results.py runs the solver at epsilon=0.01 and 0.03, both with N=64 and tolerance 1e-12. It saves fresh torus CSV files, solver logs, convergence histories, diagnostic tables, and summary.json under results/. The four figures (PDF and PNG) and a small TeX file of measured values are written to ../../slides/images/kam-*. The standalone solver itself still needs only C++.

The plots show:

  • Sampled residual versus quasi-Newton correction count, on the 2N work grid.
  • The two action displacements over the parameter torus for epsilon=0.03.
  • Fourier coefficient envelopes and the residual of truncations of the same converged torus to N=8,16,32,64, each checked on a 4N grid. These are not independent Newton runs at the smaller resolutions.
  • K(theta0+omega*t) against independent RK4 integration on 0 <= t <= 20, with time-step halving to distinguish integration error from torus error.

The script independently recomputes the final fine-grid residual using NumPy and checks the trajectory agreement before completing. Its figure labels and slide captions refer to sampled numerical errors, not rigorous sup-norm bounds. It leaves the default torus.csv from manual runs untouched.