## The FPU problem – part 3

June 17, 2013

### 2.3. Approximation with Symplectic Algorithms

The reason why we are focusing on symplectic algorithms is because for the theoretical result reported in 1994.
This result states that for any given symplectic map ${\Psi_\tau}$, analytic and ${\tau}$-close to the identity, there exists an analytic autonomous system ${H_\tau}$ whose time-1 flow ${\Phi_{H_\tau}}$ is close to ${\Psi_\tau}$, more precisely it differs from ${\Psi_\tau}$ for a quantity exponentially small in ${1/\tau}$.
The assumption that the map ${\Psi_\tau}$ is near to the identity allows one to consider a formal (ie, possibly nonconverging) solution,

$\displaystyle H_\tau = \tau \; h_1 + \tau^2 \; h_2 + \cdots,$

s.t. its formal time-one flow agrees order by order with ${\Psi_\tau}$. But whereas the formal solution can be found directly, by constructing order by order the Hamiltonian vector field of the interpolating Hamiltonian ${H_\tau}$, the solution found in 1994 has been constructed indirectly, following a Nekhoroshev-like approach, in order to prove the exponential estimate

$\displaystyle \| \Psi_\tau - \Phi_{H_\tau} \| \leq C \tau \; e^{-\tau^*/\tau},$

for some suitable ${C,\tau^*>0}$. The concise statement of the main result is the following:

Theorem 11 (Benettin, Giorgilli) Let ${\Psi_\tau}$ be a symplectic map, which is analytic and close to the identity. Then ${\exists \tau^{*}>0}$ (which depends on the constant of analicity of ${\Psi_\tau}$), and an interpolating Hamiltonian ${H_\tau}$ s.t. its time-1 flow ${\Phi_{H_\tau}}$ satisfies

$\displaystyle \| \Psi_\tau - \Phi_{H_\tau} \| = O(\tau \; e^{-\tau^*/\tau}).$

Remark 4 If the above estimate holds, then we can deduce some important consequencies:
1st property the symplectic map ${\Psi_\tau}$ admits an almost-integral of motion, up to an exponentially large number of iterations, ie ${\exists c>0}$ s.t.

$\displaystyle | H_\tau(\Psi^n_\tau(q,p) - H_\tau(q,p) | \leq c \; n \; \tau \; e^{-\tau^*/\tau}.$

\indent In particular, for any ${k \geq 1}$ we get

$\displaystyle | H_\tau(\Psi^n_\tau(q,p) - H_\tau(q,p) | \leq c \; \tau^k,$

for ${n \leq \tau^{k-1} e^{\tau^*/\tau}}$;

2nd property the distance between the n-th iterate of ${\Psi_\tau}$ and the flow ${\Phi_{H_\tau}}$,

$\displaystyle d_n := \| \Psi_\tau - \Phi_{H_\tau} \|,$

grows slowly with ${n}$, precisely ${\exists C,\lambda>0}$ s.t.

$\displaystyle d_n < C \frac{e^{k\lambda\tau}-1}{\lambda} \; e^{-\tau^*/\tau}.$

Exploiting the previous inequality we have that ${d_n = O(e^{-\tau^*/\tau})}$ for ${n=O(\tau^{-1})}$, and only on much larger ${n = O(\tau^{-2})}$ divergence possibly occures. Note that this estimate is valid in general, even if the map ${\Psi_\tau}$ admits exponential divergence of nearby orbits.

Corollary 12 The above exponential estimate and the previous remark have important consequencies from the numerical point of view.
Assume that one has to study a system described by a Hamiltonian H. Integrating numerically its equation of motion means that one approximate the time-${\tau}$ flow ${\Phi^\tau_H}$ by a mapping ${\Psi_\tau}$ close to the flow. For example, ${\Psi_\tau}$ is an algorithm of order ${k}$ we have that

$\displaystyle \| \Phi^\tau_H - \Psi_\tau \| = O(\tau^{k+1}).$

Now, if the algorithm ${\Psi_\tau}$ is symplectic, then by previous theorem we know that there exists an interpolating Hamiltonian ${H_\tau}$, and we denote ${K_\tau=\tau^{-1}H_\tau}$. For any ${t}$ one has ${\Phi^t_{H_\tau} = \Phi^{\tau t}_{K_\tau}}$, hence ${K_\tau}$ is a time-${\tau}$ interpolating Hamiltonian; this leads to

$\displaystyle \| \Psi_\tau - \Phi^\tau_{K_\tau} \| \leq C \; \tau \; e^{-\tau^*/\tau}$
Besides, one has to remember that also ${K_\tau}$ is close to the original Hamiltonian ${K}$, precisely

$\displaystyle \|H-K_\tau\| \leq C'\tau^k.$

Finally, by the above remark we have that the algorithm ${\Psi_\tau}$ follows the orbits of the interpolating Hamiltonian ${K_\tau}$ for a number of iterations ${n=O(\tau^{-1})}$ with exponentially small error (ie, for common values of ${\tau}$ the error is much smaller than the computer roundoff error).

Remark 5 Now we are able to understand what actually happens when we use a symplectic integrator of order ${k \geq 1}$ for the numerical integration of the equations of motion of an Hamiltonian system:
step 1. first we replace the original Hamiltonian ${H}$ by another autonomous Hamiltonian ${K_\tau}$, which is ${\tau^k}$-close to it;
step 2. then we compute exactly (ie, within the roundoff error) the orbits of ${K_\tau}$ through its time-1 flow (and the approximation is reliable up to an exponentially large number of iterations).
Note that this explains the good properties of energy conservation showed by symplectic algorithms. Indeed, if we set ${H(t):=H(q(t),p(t))}$, then ${\exists C''>0}$ s.t.

$\displaystyle |H(t)-H(0)| \leq |H(t)-K_\tau(t)| + |K_\tau(t)-K_\tau(0)| + |K_\tau(0)-H(0)| \leq C''\tau^k, \\$
almost uniformly in ${t}$ (ie, up to times exponentially long in ${1/\tau}$).

Another reason that explains why splitting algorithms are used to integrate numerically the FPU model is indeed the fact that these algorithms preserve the resonance relations of the model, which is a crucial property for the dynamics, since it determines the interactions between normal modes.
Remark 6 FPU model Resonance Relations
In many numerous experiments, including Fermi’s ones, initial data are taken with all energy in few modes of low order; in this case, one has to remember that for large ${N}$ (and this is what happens in the thermodynamic limit) the low frequencies are almost resonant, being integer multiples of the lowest one, so for ${1 \leq j \ll N}$ we get

$\displaystyle \omega^2_j \sim j \frac{\pi}{N+1} - \frac{1}{24} \big( \frac{j \pi}{N+1} \big)^3;$

$\displaystyle |\omega_j-j\omega_1| = O \big(\frac{j^3}{N^3}\big).$

This situation may recall the typical presence of small denominators in perturbation theory. If we assume that the integration step is small (for example ${\tau \ll 1/4}$), then also the new frequencies will satisfy the same resonance relation,

$\displaystyle |\omega_j-j\omega_1| = \frac{j(j^2-1)}{6} \big( \frac{\pi}{N+1} \big)^3 \big( \frac{1}{4}-\tau^2 \big) + \cdots = O(j^3/N^3).$