Numerical Methods in Accelerator Physics
Lecture Series by Dr. Adrian Oeftiger
Lecture 3
Run this notebook online!
Interact and run this jupyter notebook online:
Also find this lecture rendered as HTML slides on github $\nearrow$ along with the source repository $\nearrow$.
Run this first!
Imports and modules:
from config import (np, plt, hamiltonian,
plot_hamiltonian, solve_leapfrog, dt)
import PyNAFF
%matplotlib inline
Refresher!
- coordinate system: $\zeta=(x,x',y,y',z,\delta)$
- rms emittance (phase-space area, thermal energy)
- non-linearities and filamentation: microscopic level (Liouville theorem) vs. macroscopic level (emittance growth)
- emittance preservation and injection mismatch
- modelling error (vs. discretisation error)
- discrete frequency analysis:
- FFT: Fast Fourier Transform
- NAFF: Numerical Analysis of Fundamental Freqencies
Today!
- Chaos and Early Indicators
- Numerical artefacts: rounding and truncation
Non-linearities
Non-linearities in particle dynamics can originate from e.g.:
- stray fields in drift spaces
- higher-order multipole components in dipole magnets (steering) and quadrupole magnets (focusing)
- higher-order multipole magnets (sextupoles, octupoles) used to control various properties of the beam
- effects of beam fields on individual particles within the beam (space-charge forces, beam-beam effects when colliding two beams)
Effects can be varied and quite dramatic! $\implies$ require understanding of nonlinear dynamics to optimise design and operation of many accelerator systems!
Effects of Non-linearities
Design accelerator based on linear quasi-periodic oscillations $\implies$ focusing around reference particle, phase-space trajectory should be an ellipse!
Non-linearities in beam lines can have several effects:
- shape of phase-space ellipse can become distorted
- phase advance per period (frequency) can depend on oscillation amplitude
- motion can be stable for small amplitude but unstable at large amplitude: chaotic motion
- features such as "phase-space islands" can appear in phase-space portrait
Chaotic Dynamics
Tentative definition:
- sensitive to initial conditions: butterfly effect
- qualitatively recurring
but chaos is not:
- (strictly) periodic or regular
- "escalating"
... what is the problem?
$\implies$ motion not predictable despite of deterministic dynamics! (correct initial condition? numerical precision? ...)
$\implies$ cannot exclude sudden changes of amplitude on long-term time scales (synchrotron storage times) from short-term simulations in case of chaotic motion!
Deterministic Chaos
Necessary condition: non-linearity
"The exponential divergence or convergence of nearby trajectories (Lyapunov exponents) is conceptually the most basic indicator of deterministic chaos."
M. Sano and Y. Sawada (1985)
Phys. Rev. Lett. 55, 1082
One possible conceptual approach:
(Maximum) Lyapunov Exponent
Consider two nearby trajectories $\zeta_1, \zeta_2$ evolving over path length $s$:
- infinitesimal distance $|\zeta_1-\zeta_2|\rightarrow 0$
- infinite evolution $s\rightarrow \infty$
Chaotic behaviour if distance grows or shrinks exponentially:
$$|\zeta_1(s_0 + s) - \zeta_2(s_0 + s)| = \underbrace{|\zeta_1(s_0) - \zeta_2(s_0)|}\limits_{\Delta\zeta(s_0)} \cdot e^{\lambda_1 s}$$$\implies$ $\lambda_1$: maximum Lyapunov exponent
Lyapunov Spectrum
When chaotic: perturbation aligns along direction of strongest expansion (or weakest contraction), associated with maximum Lyupunov exponent
$\implies$ orthogonal directions for further Lyapunov exponents $\lambda_i$
$\implies$ as many $\lambda_i$ as phase-space coordinates
Example: Lorenz Attractor
Studied by E. Lorenz in 1963 for atmospheric convection:
$$\begin{align}
\frac{dx}{dt}&=\sigma(y-x) \\
\frac{dy}{dt}&=x(\rho - z) - y \\
\frac{dz}{dt}&=xy - \beta z
\end{align}$$Lorenz used $\sigma=10, \beta=\frac{8}{3}, \rho=28$ and investigated chaotic motion.
Chaotic Motion in Lorenz Attractor
Two identical Lorenz oscillators with initial conditions. Second oscillator is slightly perturbed ($10^{-14}$) at $t = 30$:
Returning to Discrete Pendulum
N = 11
thetas = np.linspace(0, 0.99 * np.pi, 11)
ps = np.zeros_like(thetas)
plt.scatter(thetas, ps, c='b', marker='.')
plot_hamiltonian()
Time evolution
Numerical integration of equations of motion for distribution of pendulums, using leapfrog (drift+kick+drift):
n_steps = 1000
results_thetas = np.zeros((n_steps, N), dtype=np.float32)
results_thetas[0] = thetas
results_ps = np.zeros((n_steps, N), dtype=np.float32)
results_ps[0] = ps
for k in range(1, n_steps):
results_thetas[k], results_ps[k] = solve_leapfrog(results_thetas[k - 1], results_ps[k - 1])
plt.scatter(results_thetas, results_ps, c='b', marker='.', s=1)
plot_hamiltonian()
Chaos near Unstable Fixed Point
Let us investigate the unbounded, continuously rotating pendulum, with an energy just above the separatrix value:
N = 2
thetas = np.pi * np.ones(N, dtype=np.float32)
ps = np.linspace(0.01, 0.05, N)
results_thetas2 = np.zeros((n_steps, N), dtype=np.float32)
results_thetas2[0] = thetas
results_ps2 = np.zeros((n_steps, N), dtype=np.float32)
results_ps2[0] = ps
for k in range(1, n_steps):
results_thetas2[k], results_ps2[k] = solve_leapfrog(results_thetas2[k - 1], results_ps2[k - 1])
We plot subsequent passages with the ever increasing $\theta$ projected to the interval $\bigl[0, 2\pi\bigr)$ :
plt.scatter(results_thetas2 % (2 * np.pi), results_ps2, c='b', marker='.', s=1)
plt.scatter([np.pi], [0], c='r', marker='o')
plt.xlim(np.pi - 0.1, np.pi + 0.1)
plt.ylim(-0.01, 0.1)
plt.xlabel(r'$\theta$')
plt.ylabel(r'$p$');
$\implies$ Observation near the red unstable fixed point:
The phase-space trajectory becomes non-regular!
$\implies$ The discrete pendulum dynamics features chaotic motion.
(The discrete pendulum equations are also known as the Chirikov standard map.)
Qualitative Investigation of Local Lyapunov Exponent
First investigate around stable fixed point $\theta=0$, $p=0$.
(Later come back here and investigate the around unstable fixed point $\theta=\pi$, $p=0$.)
dist = 1e-10
thetas = np.pi * np.ones(N, dtype=np.float64) # change this line
ps = np.array([0.001, 0.001 + dist], dtype=np.float64)
n_steps = 100000
results_thetas3 = np.zeros((n_steps, N), dtype=np.float64)
results_thetas3[0] = thetas
results_ps3 = np.zeros((n_steps, N), dtype=np.float64)
results_ps3[0] = ps
for k in range(1, n_steps):
results_thetas3[k], results_ps3[k] = solve_leapfrog(results_thetas3[k - 1], results_ps3[k - 1])
plt.plot(results_thetas3[:, 0], label='$p_{ini}$')
plt.plot(results_thetas3[:, 1], label='$p_{ini} + \Delta p_{ini}$')
# plt.xlim(100000-1000, 100000)
plt.xlabel('Steps $k$')
plt.ylabel(r'$\theta$')
plt.legend();
Euclidean Distance $|(\theta_2,p_2)-(\theta_1,p_1)|$
results_dist = np.sqrt(
(results_thetas3[:, 1] - results_thetas3[:, 0])**2 +
(results_ps3[:, 1] - results_ps3[:, 0])**2
)
plt.plot(results_dist)
plt.yscale('log')
plt.xlabel('Steps $k$')
plt.ylabel('Phase-space distance');
$\implies$ how do you interpret this picture?
Exercise: Go back and investigate the distance evolution for two trajectories near to the unstable fixed point!
Evaluating Local Maximum Lyapunov Exponent
Nearly exponential increase over two periods (first 4000 steps), then bounded by system size.
Local Maximum Lyapunov exponent estimated by simple linear regression:
$$\lambda_1 \approx \mathrm{slope}\left(\frac{1}{k\Delta t} \ln\left(\frac{|(\theta_2,p_2)-(\theta_1,p_1)|}{10^{-10}}\right)\right)$$B, A = np.polyfit(
x=np.arange(n_steps)[:4000],
y=np.log(results_dist[:4000] / dist),
deg=1
)
B
0.000934171219815854
plt.plot(results_dist[:4000] / dist)
plt.plot(np.exp(B * np.arange(n_steps)[:4000] + A))
plt.yscale('log');
Frequency Diffusion
Use NAFF algorithm as early indicator of chaotic motion in periodic systems:
vs.
Further reading: "Introduction to Frequency Map Analysis" by J. Laskar in Hamiltonian Systems with Three or More Degrees of Freedom, Springer (1999), pp. 134-150
Example: Earth-Moon System
Precession of Earth is stabilised by presence of Moon:
Example: CERN LHC
Frequency Map Analysis (FMA) to identify diffusion due to non-linear resonances:
$\implies$ always seek to maximise dynamic aperture in accelerators
$\implies$ use of "early chaos indicators" such as Lyapunov exponents or FMA: identification and mitigation/correction of sources
Try Concept on Discrete Pendulum
Can we detect frequency diffusion for the discrete pendulum?
n_steps = 100000
theta_ini = np.pi - 0.001
p_ini = 0
results_thetas4 = np.zeros(n_steps, dtype=np.float64)
results_thetas4[0] = theta_ini
results_ps4 = np.zeros(n_steps, dtype=np.float64)
results_ps4[0] = p_ini
for k in range(1, n_steps):
results_thetas4[k], results_ps4[k] = solve_leapfrog(results_thetas4[k - 1], results_ps4[k - 1])
plt.plot(results_thetas4)
plt.xlim(0, 1000)
plt.xlabel('Steps $k$')
plt.ylabel(r'$\theta$');
window_length = 1000
freqs_naff = []
for signal in results_thetas4.reshape((n_steps // window_length, window_length)):
freq_naff = PyNAFF.naff(signal - np.mean(signal), turns=window_length, nterms=1)[0, 1]
freqs_naff += [freq_naff]
plt.plot(np.arange(0, n_steps, window_length), freqs_naff, ls='none', marker='.')
plt.xlabel('Steps $k$')
plt.ylabel('NAFF determined frequency');
plt.figure(figsize=(10, 5))
plt.hist(freqs_naff);
Summing up Concepts
- sources for non-linearities in accelerators
- deterministic chaos: (bounded) deterministic dynamics with a positive Lyapunov exponent
- sensitivity to initial conditions
- prevents prediction for long-term time scales
- early indicators of chaos:
- (maximum) Lyapunov exponent
- Frequency Map Analysis (FMA)
Floating Point Representation
Represent real numbers $r\in\mathbb{R}$ on computers by "floats":
- finite and discrete set of numbers
- floating-point form: $r=\underbrace{6.022}\limits_\text{significand}\times 10{\underbrace{{}^{23}}\limits_\text{exponent}}$
Todays standard: IEEE 754
$r=\pm a \times 2^{e}$
e.g. a "single-precision" float (FP32):
- 1 sign bit
- $t=23$ bits for significand $a$
- 8 bits for exponents in range $e\in[-126, 127]$
- $\pm\infty$ and NaN
Numerical Artefacts
Floating point representation of numbers on computers $\leadsto$ numerical artefacts:
- rounding: storing the least significant bit depending on arithmetic operation
- truncating: e.g. π cannot be stored exactly but is truncated
Smallest "machine" precision in resolving real numbers: $2^{-t}$
- FP32: $2^{-23}\approx 10^{-7}$
- FP64: $2^{-52}\approx 10^{-15}$
Think about what this means for simulating long-term dynamics, e.g. for LHC storage times...
(Drastic) Example of Truncation Error
The python library mpmath works with arbitrary numerical precision (workdps for decimal precision):
import mpmath as mp
with mp.workdps(7):
x = mp.mpf(10) / 9 # == 1.11...
for _ in range(30):
print (x)
x = (x - 1) * 10
1.111111 1.111111 1.11111 1.111104 1.111045 1.110449 1.104488 1.044884 0.4488373 -5.511627 -65.11627 -661.1627 -6621.627 -66226.27 -662272.7 -6622737.0 -6.622738e+7 -6.622738e+8 -6.622738e+9 -6.622738e+10 -6.622738e+11 -6.622738e+12 -6.622738e+13 -6.622738e+14 -6.622738e+15 -6.622738e+16 -6.622738e+17 -6.622738e+18 -6.622738e+19 -6.622738e+20
Overview: 3 Sources of Simulation Errors
- Discretisation error
- Modelling error
- Numerical artefacts
- (Input error...)
Summary
- sources for non-linearities in accelerators
- deterministic chaos
- early indicators:
- (maximum) Lyapunov exponent
- Frequency Map Analysis
- numerical artefacts: rounding and truncation, machine precision
Literature
- chapter 2.1 in Justin Solomon, Numerical Algorithms