Home / KAM course

Mathematical contract and theorem-to-code map

1. Scope and conventions

The target is Theorem 2.18 of Figueras–Haro (2024), specialized using Remarks 2.20–2.21 to d=n=2, no additional first integrals, canonical geometry, and Diophantine exponent tau=1. The angle-action/affine-lift adaptation is the one described in the paper’s Section 2.2.

Ambient coordinates are unchanged from KAMExample: z=(phi1,phi2,I1,I2), with the physical angles of period 2 pi. For the paper’s conventions choose

\[ \Omega=J_p=\begin{pmatrix}0&-I_2\\ I_2&0\end{pmatrix},\quad G=I_4,\quad \alpha=I\cdot d\varphi,\quad \Omega=d\alpha,\quad X_H=\Omega^{-1}\nabla H. \]

Here the displayed matrix represents the two-form. Thus the vector field is (partial_I H,-partial_phi H), as in the numerical example. All constant geometry norms are 1, their derivative bounds are 0, all additional-integral bounds are 0, and c_DPhi=1.

Only the parameter is changed to the paper’s unit torus:

\[ x=\theta/(2\pi),\qquad K(x)=(2\pi x+u(2\pi x),\omega_r+v(2\pi x)),\qquad \omega=\omega_r/(2\pi). \]

Consequently

\[ L=D_xK=2\pi D_\theta K,\quad \rho=\rho_r/(2\pi),\quad \rho_\infty=\rho_{\infty,r}/(2\pi),\quad \delta=(\rho-\rho_\infty)/6. \]

The substep cuts have ratio 1/2; an iteration spends 3 delta_j, so the limiting width is exactly rho_infinity. The constants use a=2. Ambient coordinates, domain distances, the Hamiltonian and its derivatives are not rescaled. Averages have normalization 1/M^2, with no additional 2 pi factor.

2. Exact candidate and arithmetic

All retained Fourier coefficients are exact dyadics; both Nyquist lines are zero and conjugate symmetry holds exactly. validateCandidate() enforces this contract, including for direct calls to the bounds API. The action offset is the exact algebraic frequency. The candidate is therefore a specified real analytic embedding once injectivity has been checked.

MPFI encloses real quantities; complex numbers are pairs of real intervals. Every proof operation uses these intervals. upper() deliberately replaces an interval evaluation of a bound by an exact dyadic majorant. lower() is used for distances and positive gaps. No ambiguous interval order is used: strict tests compare an upper endpoint with a lower endpoint. Division by an interval containing zero, empty intervals, NaNs and infinities cannot pass a proof test.

The numerical proposal stage may take midpoints and freeze coefficients. Its output is a new candidate, not a claimed enclosure of an unknown exact torus. The subsequent certification is independent of the numerical proposal.

3. Diophantine and Rüssmann bounds

Put phi=(1+sqrt(5))/2 and psi=(1-sqrt(5))/2. For nonzero integer (p,q),

\[ (p+q\phi)(p+q\psi)=p^2+pq-q^2\in\mathbb Z\setminus\{0\},\qquad |p+q\psi|\le |p|+|q|. \]

Thus omega_r=(1,phi) satisfies |k.omega_r| >= 1/|k|_1 for every nonzero integer mode. In unit-parameter coordinates,

\[ \gamma=1/(2\pi),\qquad\tau=1. \]

No finite scan of divisors is used as a replacement for this argument.

Lemma 2.14 and Corollary 2.15 specialize to the uniform constants

\[ \widehat c_R =\sqrt{\frac{2^3\zeta(2,2)\Gamma(3)}{(2\pi)^4}} =\frac{\sqrt{\pi^2/6-1}}{\pi^2},\qquad \widehat c_R^1=4\widehat c_R. \]

These hold at every strip loss used in the infinite iteration. Both are evaluated with outward rounding, rather than approximated by finite divisor sums. c_R and c_R1 in the report are their computed upper bounds.

4. Norms, Gram inverse, and embedding

Scalar functions use strip sup norms. Matrix bounds follow the paper’s norm

\[ \|A\|_\rho=\max_i\sum_j\sup_{|\Im\theta_\ell|\le\rho}|A_{ij}(\theta)|. \]

Each finite polynomial is bounded by its weighted Fourier sum. Derivative coefficients are exactly ik_j c_k in radians. Matrix transpose norms are bounded separately; a vector and its transpose do not have the same norm.

Writing U=D_theta u and V=D_theta v,

\[ L^TL=(2\pi)^2\left(I+U+U^T+U^TU+V^TV\right). \]

At both the initial and wider strips the code checks

\[ r_B=b_U+b_{U^T}+b_{U^T}b_U+b_{V^T}b_V<1. \]

The Banach-algebra Neumann series gives a holomorphic inverse throughout the strip, with

\[ b_B=\frac1{(2\pi)^2(1-r_B)},\qquad b_N\le b_L b_B,\qquad b_{N^T}\le b_B b_{L^T},\qquad N=J_pLB. \]

This is a transpose Gram matrix, not a Hermitian matrix. Its inverse is symmetric, so B and B^T share the same norm bound.

The additional sufficient test ||D_theta u|| < 1 makes the angular map theta+u(theta) injective modulo the period lattice: apply the Lipschitz estimate on the convex universal-cover strip, translating one point by a lattice element if necessary. This proves embedding, not just immersion. At the limit the code checks

\[ b_{D_\theta u}+r_{D_xK}/(2\pi)<1. \]

5. Explicit domains and Hamiltonian bounds (H1)

The complex domain U0 is defined by |Im phi_i| < R0_i and |I_i| < R0_(i+2). Its radii are upward-rounded bounds on the corresponding candidate components, enlarged by domainMargin. The outer extension domain U has another domainMargin of reserve. The record contains all eight radii and a lower bound for dist(K(T_rho), boundary U0).

For the literal functional-independence assumption in the theorem’s real setting, use the real slice of U0 intersected with I1>0. The code checks

\[ I_1\ge1-\|v_1\|_\rho>0, \qquad 1-\|v_1\|_\rho-r_K>0. \]

Since partial H/partial I1 = I1, the Hamiltonian is regular there. The total displacement bound keeps the real iteration and its limit in this regular domain. The normal intermediate displacement is also dominated by the full one-step majorant. The larger complex analytic-extension domains are allowed to contain critical points.

Let S=max(cosh R1,cosh R2), M=cosh(R1+R2) for the outer angular radii. Then the implementation uses

\[ \begin{aligned} c_X&=\max(R_{I_1},R_{I_2},|\epsilon|(S+M)),\\ c_{DX}=c_{DX^T}&=\max(1,|\epsilon|(S+2M)),\\ c_{D^2X}&=|\epsilon|(S+4M),\\ c_{T_H}&=1+|\epsilon|(S+2M),\qquad c_{DT_H}=c_{D^2X}. \end{aligned} \]

These are componentwise derivative-sum bounds as required by H1. The equal transpose bound for DX follows from the symmetry of the force Jacobian.

6. The theorem’s torsion and its average (H2)

For

\[ A=DX_H=\begin{pmatrix}0&I\\ \epsilon C&0\end{pmatrix},\quad C=\begin{pmatrix} \cos\varphi_1+\cos(\varphi_1-\varphi_2)&-\cos(\varphi_1-\varphi_2)\\ -\cos(\varphi_1-\varphi_2)&\cos\varphi_2+\cos(\varphi_1-\varphi_2) \end{pmatrix}, \]

Eq. (2.9) is

\[ T_H=\Omega A-A\Omega =\operatorname{diag}(-(I+\epsilon C),I+\epsilon C),\qquad T=N^TT_H(K)N. \]

This formula is used even off invariance. The derivative-based torsion in the double-precision course solver is a different off-torus quantity.

Grid values of L, B, N and T are enclosed with MPFI, including each 2-by-2 inverse. The sampled mean is enlarged by the analytic mean-alias bound below in every entry. A direct interval inverse of this enlarged average then bounds ||<T>^-1||. Invertibility on a grid alone is not used to certify strip-wide invertibility.

The strict budgets are

\[ \begin{aligned} \sigma_L=\sigma_{L^T}&=f\max(b_L,b_{L^T}),\\ \sigma_N=\sigma_{N^T}&=f\max(b_N,b_{N^T}),\\ \sigma_B&=f b_B,\qquad \sigma_{\langle T\rangle^{-1}}=f b_{\langle T\rangle^{-1}}, \end{aligned} \]

rounded upward, with f=sigmaFactor>1. Every gap between a budget and an initial bound is required to have a strictly positive lower bound.

7. Rigorous FFT, interpolation and mean-alias bounds

The FFT implements

\[ \widetilde f_k=M^{-2}\sum_{i,j=0}^{M-1} f(2\pi i/M,2\pi j/M)e^{-2\pi\mathrm i(k_1i+k_2j)/M}. \]

It encloses the discrete coefficients. For a function bounded by F on a wider strip w, its true coefficients satisfy |fhat_k| <= F exp(-w |k|_1).

Let q=exp(-(w-rho)) and use the half-open FFT rectangle. The exact alias identity is ftilde_k = sum_m fhat_(k+Mm). Each omitted mode folds to a representative with no larger l1 index. Consequently

\[ \begin{aligned} \|f-I_Mf\|_\rho &\le2F\sum_{k\text{ outside}}q^{|k|_1}\\ &\le 8F\,q^{M/2}\frac{1+q}{(1-q)^2}. \end{aligned} \]

The union bound includes both axes and the Nyquist modes; its overcounting is conservative. It is a simpler, less sharp bound derived from the same alias identity as Figueras–Haro–Luque, Section 3. Thus

\[ b_E=\max_i\sum_{k\text{ in grid}}|\widetilde E_{i,k}|e^{\rho_r|k|_1} +C_{\mathrm{grid}}M_E \]

bounds the full continuous residual. All weights and sums are enclosed. The independent wider-strip bound M_E comes from

\[ E_\varphi=v-D_\theta u\,\omega_r,\qquad E_I=\epsilon(\sin\varphi_1+\sin(\varphi_1-\varphi_2), \sin\varphi_2-\sin(\varphi_1-\varphi_2))-D_\theta v\,\omega_r. \]

The exact constant action offset has been canceled symbolically.

For the mean, let q=exp(-w M). Only frequencies in M Z^2 alias into zero, so

\[ |\langle f\rangle-M^{-2}\sum f(\theta_{ij})| \le F\left[\left(\frac{1+q}{1-q}\right)^2-1\right] =\frac{4Fq}{(1-q)^2}. \]

The wider-strip torsion bound is M_T <= b_Nt(w) c_Th(w) b_N(w). It bounds every matrix entry, so the above error may safely inflate each sampled mean entry before inversion. These computations require only the finite candidate, not an assumed analyticity width of the unknown exact torus.

8. Appendix B implementation map

The code uses g=gamma*delta,

\[ e_* = \max(b_{\eta^L},b_{\eta^N}/g),\quad e_*'=e_*/\delta,\quad \chi=e_*/(\gamma\delta^2). \]

The projected residual majorants are conservatively obtained as b_etaL=b_Nt*b_E, b_etaN=b_DKt*b_E. Their exact zero-average compatibility is the geometric identity for an exact symplectic Hamiltonian, as in the paper; no numerical zero-mode scan is used to prove it.

Terms multiplying first-integral or nonconstant-geometry derivatives vanish in this specialization. In particular Q_etaL4=0.

There are two different pair conventions in the tables. The E, E^T, DE, (DE)^T and redNN coefficients combine as C_L+g C_N; the already-weighted normal/tangent correction pairs and redLN combine as C_L+C_N. The code preserves this distinction.

Conservative common-transpose bounds

For the intermediate normal correction the code uses the common bound

\[ C_{\Delta\bar{DK},\mathrm{common}} =4\sigma_N(C_{\widehat\xi^N_0}^L+C_{\widehat\xi^N_0}^N+c_R^1). \]

The factor 4 dominates the factors d=2 and n=2 in Table 2. The final common derivative bound uses 2n=4, covering both (3.64) and (3.65). Since L=DK, the same bound applies to changes in L and L^T.

For constant canonical geometry, directly transposing N=J_p L B gives

\[ C_{\Delta N^T}\le\sigma_B C_{\Delta L^T} +C_{\Delta B}\sigma_{L^T}. \]

The common budgets implement this bound and its non-transposed counterpart. This avoids identifying row and column norms or relying on the untransposed symbols printed in Table 3’s row (3.71). Subsequent constants are enlarged consistently. The numerical certificate therefore uses conservative majorants, not a claim of reproducing the paper’s sharpest constants.

Final inequalities and distances

Set

\[ Q_\eta=\max(Q_{\eta^L},2Q_{\eta^N}),\qquad \kappa=4Q_\eta\chi, \]

and compute C_Delta from all six positive geometry-budget gaps. Then

\[ \mathcal C=\max\{\gamma\delta C_{sym},\ C_{\xi^L},\ \delta C_{dist}+2Q_\eta,\ C_\Delta+4Q_\eta\}. \]

The code requires nu<1, C*chi<1, kappa<1, and checks each component of the maximum explicitly. It then computes

\[ r_K=\frac{2}{2-\kappa}\frac{C_{\Delta K}}{\gamma\delta}e_*,\qquad r_{DK}=\frac{C_{\Delta DK,\mathrm{common}}}{1-\kappa}\chi, \]

and the analogous inverse/frame distances. The additional limiting embedding and regular real-domain reserves must also pass. All comparisons use strict interval endpoint inequalities.

References checked locally

• Figueras_Haro_A_modified_parameterization_method_for_invariant_Lagrangian_tori_for_partially_integrable_Hamiltonian_systems_Physica_D_462_2024.pdf — Theorem 2.18, Remarks 2.20–2.24, Lemma 2.14, Corollary 2.15, Lemma 3.6, Section 3.3, Appendix B.
• Figueras_Haro_Luque_Rigorous_Computer-Assisted_Application_of_KAM_Theory:_A_Modern_Approach_Foundations_of_Computational_Mathematics_1_71_2016.pdf — Section 3 (alias identity and interpolation bounds) and Section 5 (validation methodology). This local filename is the online-first version; the published volume is Found. Comput. Math. 17 (2017), 1123–1193.