{ "cells": [ { "cell_type": "markdown", "id": "91890c18", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "

Numerical Methods of Accelerator Physics

\n", "

Lecture by Dr. Adrian Oeftiger

\n", "\n", "\n", "\n", "

Part 3: 04.11.2022

" ] }, { "cell_type": "markdown", "id": "9ce0e12b", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "

Run this notebook online!

\n", "\n", "Interact and run this jupyter notebook online:\n", "\n", "
\n", "1. via the public mybinder.org service:
\n", "\n", "

\n", "\n", "

\n", "
\n", "\n", "
\n", "2. on the local TU Darmstadt jupyterhub $\\nearrow$ (using your TU ID)\n", "\n", "$\\implies$ make sure you installed all the required python packages (see the [README](./README.md))!\n", "
\n", "\n", "Finally, also find this lecture rendered [as HTML slides on github $\\nearrow$](https://aoeftiger.github.io/TUDa-NMAP-03/) along with the [source repository $\\nearrow$](https://github.com/aoeftiger/TUDa-NMAP-03)." ] }, { "cell_type": "markdown", "id": "e1649e74", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "

Announcements

\n", "\n", "Starting flipped classroom next week:\n", "- first video material: lecture part 4\n", " - available on moodle by Monday night, 7.11.2022\n", " - please watch it until Friday, 11.11.2022\n", " - note your questions\n", " \n", "- Friday, 11.11.2022:\n", " - lecture time: discuss open questions & solve problems together" ] }, { "cell_type": "markdown", "id": "abaa09ce", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "

Refresher!

\n", "\n", "- coordinate system: $\\zeta=(x,x',y,y',z,\\delta)$\n", "- rms emittance (phase-space area, thermal energy)\n", "- non-linearities and filamentation: microscopic level (Liouville theorem) vs. macroscopic level (emittance growth)\n", "- emittance preservation and injection mismatch\n", "- modelling error (vs. discretisation error)\n", "- discrete frequency analysis:\n", " - FFT: Fast Fourier Transform\n", " - NAFF: Numerical Analysis of Fundamental Freqencies" ] }, { "cell_type": "markdown", "id": "e7a5ad4a", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "

Today!

\n", "\n", "1. Chaos and Early Indicators\n", "2. Numerical artefacts: rounding and truncation" ] }, { "cell_type": "markdown", "id": "f7a570fa", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "

Non-linearities

\n", "\n", "Non-linearities in particle dynamics can originate from e.g.:\n", "\n", "- stray fields in drift spaces\n", "- higher-order multipole components in dipole magnets (steering) and quadrupole magnets (focusing)\n", "- higher-order multipole magnets (sextupoles, octupoles) used to control various properties of the beam\n", "- effects of beam fields on individual particles within the beam (space-charge forces, beam-beam effects when colliding two beams)\n", "\n", "Effects can be varied and quite dramatic! $\\implies$ require understanding of nonlinear dynamics to optimise design and operation of many accelerator systems!" ] }, { "cell_type": "markdown", "id": "102fe6c4", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "

Effects of Non-linearities

\n", "\n", "Design accelerator based on linear quasi-periodic oscillations $\\implies$ focusing around reference particle, phase-space trajectory should be an ellipse!\n", "\n", "Non-linearities in beam lines can have several effects:\n", "- shape of phase-space ellipse can become distorted\n", "- phase advance per period (frequency) can depend on oscillation amplitude\n", "- motion can be stable for small amplitude but unstable at large amplitude: chaotic motion\n", "- features such as \"phase-space islands\" can appear in phase-space portrait" ] }, { "cell_type": "markdown", "id": "542d4bad", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "

Chaotic Dynamics

\n", "\n", "\"Double\n", "\n", "Tentative definition:\n", " \n", "- sensitive to initial conditions: butterfly effect\n", "- qualitatively recurring\n", "\n", "but chaos is not:\n", "\n", "- (strictly) periodic or regular\n", "- \"escalating\"\n", "\n", "

image by Ari Rubinsztejn

" ] }, { "cell_type": "markdown", "id": "425525be", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "
\n", "

... what is the problem?

\n", "\n", "$\\implies$ motion not predictable despite of deterministic dynamics! (correct initial condition? numerical precision? ...)\n", "\n", "$\\implies$ cannot exclude sudden changes of amplitude on long-term time scales (synchrotron storage times) from short-term simulations in case of chaotic motion!\n", "
" ] }, { "cell_type": "markdown", "id": "a89b0da5", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "

Deterministic Chaos

\n", "\n", "Necessary condition: non-linearity\n", "\n", "\"The exponential divergence or convergence of nearby trajectories (Lyapunov exponents) is conceptually the most basic indicator of deterministic chaos.\"\n", "

M. Sano and Y. Sawada (1985)
Phys. Rev. Lett. 55, 1082

\n", "\n", "One possible conceptual approach: \n", "
deterministic chaos $\\doteq$ \"bounded, deterministic dynamics with a positive Lyapunov exponent.\"
\n" ] }, { "cell_type": "markdown", "id": "179822c0", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "

(Maximum) Lyapunov Exponent

\n", "\n", "Consider two nearby trajectories $\\zeta_1, \\zeta_2$ evolving over path length $s$:\n", "\n", "- infinitesimal distance $|\\zeta_1-\\zeta_2|\\rightarrow 0$\n", "- infinite evolution $s\\rightarrow \\infty$\n", "\n", "Chaotic behaviour if distance grows or shrinks exponentially:\n", "\n", "$$|\\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}$$\n", "\n", "$\\implies$ $\\lambda_1$: maximum Lyapunov exponent\n", "\n", "
$$ \\lambda_1 \\doteq \\lim\\limits_{s\\rightarrow \\infty}\\,\\lim\\limits_{\\Delta\\zeta(s_0)\\rightarrow 0}\\,\\frac{1}{s}\\,\\ln\\left(\\frac{\\Delta\\zeta(s_0+s)}{\\Delta\\zeta(s_0)}\\right)$$
" ] }, { "cell_type": "markdown", "id": "04563e27", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "

Lyapunov Spectrum

\n", "\n", "When chaotic: perturbation aligns along direction of strongest expansion (or weakest contraction), associated with maximum Lyupunov exponent \n", "\n", "$\\implies$ orthogonal directions for further Lyapunov exponents $\\lambda_i$\n", "\n", "$\\implies$ as many $\\lambda_i$ as phase-space coordinates" ] }, { "cell_type": "markdown", "id": "90124625", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "

Example: Lorenz Attractor

\n", "\n", "Studied by E. Lorenz in 1963 for atmospheric convection:\n", "\n", "\"Lorenz\n", "\n", "\n", "$$\\begin{align}\n", " \\frac{dx}{dt}&=\\sigma(y-x) \\\\\n", " \\frac{dy}{dt}&=x(\\rho - z) - y \\\\\n", " \\frac{dz}{dt}&=xy - \\beta z\n", "\\end{align}$$\n", " \n", "Lorenz used $\\sigma=10, \\beta=\\frac{8}{3}, \\rho=28$ and investigated chaotic motion.\n", "\n", "

image by Dan Quinn

" ] }, { "cell_type": "markdown", "id": "477918a0", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "

Chaotic Motion in Lorenz Attractor

\n", "\n", "Two identical Lorenz oscillators with initial conditions.\n", "Second oscillator is slightly perturbed ($10^{-14}$) at $t = 30$:\n", " \n", "\"Trajectories\n", "\n", "

image by G. Ansmann

" ] }, { "cell_type": "code", "execution_count": null, "id": "2c1f9d9c", "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [], "source": [ "from config import (np, plt, hamiltonian, \n", " plot_hamiltonian, solve_leapfrog, dt)\n", "%matplotlib inline" ] }, { "cell_type": "code", "execution_count": null, "id": "24d57556", "metadata": {}, "outputs": [], "source": [ "N = 11\n", "\n", "thetas = np.linspace(0, 0.99 * np.pi, 11)\n", "ps = np.zeros_like(thetas)" ] }, { "cell_type": "code", "execution_count": null, "id": "6245ca65", "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [], "source": [ "plt.scatter(thetas, ps, c='b', marker='.')\n", "\n", "plot_hamiltonian()" ] }, { "cell_type": "markdown", "id": "32d2cb24", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "

Time evolution

\n", "\n", "Numerical integration of equations of motion for distribution of pendulums, using leapfrog (drift+kick+drift):" ] }, { "cell_type": "code", "execution_count": null, "id": "14014c24", "metadata": {}, "outputs": [], "source": [ "n_steps = 1000" ] }, { "cell_type": "code", "execution_count": null, "id": "eb251771", "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "results_thetas = np.zeros((n_steps, N), dtype=np.float32)\n", "results_thetas[0] = thetas\n", "\n", "results_ps = np.zeros((n_steps, N), dtype=np.float32)\n", "results_ps[0] = ps\n", "\n", "for k in range(1, n_steps):\n", " results_thetas[k], results_ps[k] = solve_leapfrog(results_thetas[k - 1], results_ps[k - 1])" ] }, { "cell_type": "code", "execution_count": null, "id": "a3968683", "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [], "source": [ "plt.scatter(results_thetas, results_ps, c='b', marker='.', s=1)\n", "\n", "plot_hamiltonian()" ] }, { "cell_type": "markdown", "id": "77a9290a", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "

Chaos near Unstable Fixed Point

" ] }, { "cell_type": "code", "execution_count": null, "id": "a8c2e6e3", "metadata": {}, "outputs": [], "source": [ "N = 2\n", "\n", "thetas = np.pi * np.ones(N, dtype=np.float32)\n", "ps = np.linspace(0.01, 0.05, N)" ] }, { "cell_type": "code", "execution_count": null, "id": "07057a0c", "metadata": {}, "outputs": [], "source": [ "results_thetas2 = np.zeros((n_steps, N), dtype=np.float32)\n", "results_thetas2[0] = thetas\n", "\n", "results_ps2 = np.zeros((n_steps, N), dtype=np.float32)\n", "results_ps2[0] = ps\n", "\n", "for k in range(1, n_steps):\n", " results_thetas2[k], results_ps2[k] = solve_leapfrog(results_thetas2[k - 1], results_ps2[k - 1])" ] }, { "cell_type": "code", "execution_count": null, "id": "f5f60253", "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [], "source": [ "plt.scatter(results_thetas2 % (2*np.pi), results_ps2, c='b', marker='.', s=1)\n", "plt.scatter([np.pi], [0], c='r', marker='o')\n", "\n", "plt.xlim(np.pi - 0.1, np.pi + 0.1)\n", "plt.ylim(-0.01, 0.1)\n", "plt.xlabel(r'$\\theta$')\n", "plt.ylabel(r'$p$')" ] }, { "cell_type": "markdown", "id": "3ef66e0c", "metadata": {}, "source": [ "$\\implies$ observation near red unstable fixed point:\n", "- phase-space trajectory becomes non-regular, obvious during subsequent passages when wrapping $\\theta$ to interval $\\bigl[0, 2\\pi\\bigr)$" ] }, { "cell_type": "markdown", "id": "86c12fe5", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "

Qualitative Investigation of Local Lyapunov Exponent

\n", "\n", "First around stable fixed point $\\theta=0$, $p=0$ -- then around unstable fixed point $\\theta=\\pi$, $p=0$:" ] }, { "cell_type": "code", "execution_count": null, "id": "c6eaf7e6", "metadata": {}, "outputs": [], "source": [ "dist = 1e-10\n", "\n", "thetas = np.pi * np.ones(N, dtype=np.float64) # change this line\n", "ps = np.array([0.001, 0.001 + dist], dtype=np.float64)" ] }, { "cell_type": "code", "execution_count": null, "id": "7ad87d3a", "metadata": {}, "outputs": [], "source": [ "n_steps = 100000" ] }, { "cell_type": "code", "execution_count": null, "id": "2194797a", "metadata": {}, "outputs": [], "source": [ "results_thetas3 = np.zeros((n_steps, N), dtype=np.float64)\n", "results_thetas3[0] = thetas\n", "\n", "results_ps3 = np.zeros((n_steps, N), dtype=np.float64)\n", "results_ps3[0] = ps\n", "\n", "for k in range(1, n_steps):\n", " results_thetas3[k], results_ps3[k] = solve_leapfrog(results_thetas3[k - 1], results_ps3[k - 1])" ] }, { "cell_type": "code", "execution_count": null, "id": "adc3f6b3", "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [], "source": [ "plt.plot(results_thetas3[:, 0], label='$p_{ini}$')\n", "plt.plot(results_thetas3[:, 1], label='$p_{ini} + \\Delta p_{ini}$')\n", "\n", "plt.xlabel('Steps $k$')\n", "plt.ylabel(r'$\\theta$')\n", "plt.legend();" ] }, { "cell_type": "markdown", "id": "18217739", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "

Distance $|(\\theta_2,p_2)-(\\theta_1,p_1)|$

" ] }, { "cell_type": "code", "execution_count": null, "id": "39edde36", "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "results_dist = np.sqrt(\n", " (results_thetas3[:, 1] - results_thetas3[:, 0])**2 + \n", " (results_ps3[:, 1] - results_ps3[:, 0])**2\n", ")" ] }, { "cell_type": "code", "execution_count": null, "id": "7ce0195a", "metadata": { "scrolled": false }, "outputs": [], "source": [ "plt.plot(results_dist)\n", "\n", "plt.yscale('log')\n", "plt.xlabel('Steps $k$')\n", "plt.ylabel('Phase-space distance');" ] }, { "cell_type": "markdown", "id": "d4933779", "metadata": { "slideshow": { "slide_type": "-" } }, "source": [ "" ] }, { "cell_type": "markdown", "id": "96665c03", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "

Frequency Diffusion

\n", "\n", "Use NAFF algorithm as early indicator of chaotic motion in periodic systems:\n", "\n", "
\n", "Regular orbits $\\Longleftrightarrow$ fixed frequency up to numerical accuracy!
\n", "\n", "

vs.

\n", "\n", "
\n", "Non-regular orbits $\\Longleftrightarrow$ frequency evolves in time indicating chaotic diffusion of orbit!
\n", "\n", "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" ] }, { "cell_type": "markdown", "id": "87366e15", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "

Example: Earth-Moon System

\n", " \n", "Precession of Earth is stabilised by presence of Moon:\n", "\n", "\"Earth\n", "\n", "

image by J. Laskar

" ] }, { "cell_type": "markdown", "id": "96ecbb1e", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "

Example: CERN LHC

\n", "\n", "\"Frequency\n", "\n", "Frequency Map Analysis (FMA) to identify diffusion due to non-linear resonances:\n", "\n", "
Concept of dynamic aperture: smallest amplitude where particles survive (before chaos and amplitude increase can lead to loss!)!
\n", "\n", "$\\implies$ always seek to maximise dynamic aperture in accelerators\n", "\n", "$\\implies$ use of \"early chaos indicators\" such as Lyapunov exponents or FMA: identification and mitigation/correction of sources\n", "\n", "

image by Y. Papaphilippou

" ] }, { "cell_type": "markdown", "id": "4187490b", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "

Try Concept on Discrete Pendulum

\n", "\n", "Can we detect frequency diffusion in discrete pendulum?" ] }, { "cell_type": "code", "execution_count": null, "id": "0c266200", "metadata": {}, "outputs": [], "source": [ "n_steps = 100000" ] }, { "cell_type": "code", "execution_count": null, "id": "c52150c2", "metadata": {}, "outputs": [], "source": [ "theta_ini = np.pi - 0.001\n", "p_ini = 0" ] }, { "cell_type": "code", "execution_count": null, "id": "2e0c52b1", "metadata": {}, "outputs": [], "source": [ "results_thetas4 = np.zeros(n_steps, dtype=np.float64)\n", "results_thetas4[0] = theta_ini\n", "\n", "results_ps4 = np.zeros(n_steps, dtype=np.float64)\n", "results_ps4[0] = p_ini\n", "\n", "for k in range(1, n_steps):\n", " results_thetas4[k], results_ps4[k] = solve_leapfrog(results_thetas4[k - 1], results_ps4[k - 1])" ] }, { "cell_type": "code", "execution_count": null, "id": "1e16f9f6", "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [], "source": [ "plt.plot(results_thetas4)\n", "\n", "plt.xlabel('Steps $k$')\n", "plt.ylabel(r'$\\theta$');" ] }, { "cell_type": "code", "execution_count": null, "id": "c9809515", "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [], "source": [ "import PyNAFF" ] }, { "cell_type": "code", "execution_count": null, "id": "0047646d", "metadata": {}, "outputs": [], "source": [ "window_length = 1000\n", "\n", "freqs_naff = []\n", "\n", "for signal in results_thetas4.reshape((n_steps // window_length, window_length)):\n", " freq_naff = PyNAFF.naff(signal - np.mean(signal), turns=window_length, nterms=1)[0, 1]\n", " freqs_naff += [freq_naff]" ] }, { "cell_type": "code", "execution_count": null, "id": "55aaebf6", "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "plt.plot(np.arange(0, n_steps, window_length), freqs_naff, ls='none', marker='.')\n", "\n", "plt.xlabel('Steps $k$')\n", "plt.ylabel('NAFF determined frequency');\n", "\n", "# plt.ylim(1.59e-2, 1.6e-2)" ] }, { "cell_type": "code", "execution_count": null, "id": "180034d5", "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [], "source": [ "plt.figure(figsize=(10, 5))\n", "plt.hist(freqs_naff);" ] }, { "cell_type": "markdown", "id": "a2bedc61", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "

Summing up Concepts

\n", "\n", "- sources for non-linearities in accelerators\n", "- deterministic chaos: (bounded) deterministic dynamics with a positive Lyapunov exponent\n", " - sensitivity to initial conditions\n", " - prevents prediction for long-term time scales\n", "- early indicators of chaos:\n", " - (maximum) Lyapunov exponent\n", " - Frequency Map Analysis (FMA)" ] }, { "cell_type": "markdown", "id": "6f181746", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "

Floating Point Representation

\n", "\n", "Represent real numbers $r\\in\\mathbb{R}$ on computers by \"floats\":\n", "- finite and discrete set of numbers\n", "- floating-point form: $r=\\underbrace{6.022}\\limits_\\text{significand}\\times 10{\\underbrace{{}^{23}}\\limits_\\text{exponent}}$\n", "\n", "Todays standard: IEEE 754\n", "\n", "$r=\\pm a \\times 2^{e}$\n", "\n", "e.g. a \"single-precision\" float (FP32):\n", "\n", "- 1 sign bit\n", "- $t=23$ bits for significand $a$\n", "- 8 bits for exponents in range $e\\in[-126, 127]$\n", "- $\\pm\\infty$ and NaN" ] }, { "cell_type": "markdown", "id": "e5d9d0fe", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "

Numerical Artefacts

\n", "\n", "Floating point representation of numbers on computers $\\leadsto$ numerical artefacts:\n", "- rounding: storing the least significant bit depending on arithmetic operation\n", "- truncating: e.g. π cannot be stored exactly but is truncated\n", "\n", "Smallest \"machine\" precision in resolving real numbers: $2^{-t}$\n", "- FP32: $2^{-23}\\approx 10^{-7}$\n", "- FP64: $2^{-52}\\approx 10^{-15}$\n", "\n", "
Accumulation of errors $\\implies$ gradual loss of significance
" ] }, { "cell_type": "markdown", "id": "e3f79c47", "metadata": { "slideshow": { "slide_type": "-" } }, "source": [ "" ] }, { "cell_type": "markdown", "id": "7b65193c", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "

(Drastic) Example of Truncation Error

\n", "\n", "The python library `mpmath` works with arbitrary numerical precision (`workdps` for decimal precision):" ] }, { "cell_type": "code", "execution_count": null, "id": "a1e6b892", "metadata": {}, "outputs": [], "source": [ "import mpmath as mp\n", "\n", "with mp.workdps(7):\n", " x = mp.mpf(10) / 9 # == 1.11...\n", " \n", " for _ in range(30):\n", " print (x)\n", " x = (x - 1) * 10" ] }, { "cell_type": "markdown", "id": "e38a356e", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "

Overview: Sources of Simulation Errors

\n", " \n", "1. Discretisation error\n", "2. Modelling error\n", "3. Numerical artefacts\n", "4. Input error..." ] }, { "cell_type": "markdown", "id": "0b0325dc", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "

Summary

\n", "\n", "- sources for non-linearities in accelerators\n", "- deterministic chaos\n", "- early indicators:\n", " - (maximum) Lyapunov exponent\n", " - Frequency Map Analysis\n", "- numerical artefacts: rounding and truncation, machine precision" ] }, { "cell_type": "markdown", "id": "0f42e186", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "

Literature

\n", "\n", "- chapter 2.1 in Justin Solomon, [Numerical Algorithms](https://people.csail.mit.edu/jsolomon/share/book/numerical_book.pdf)" ] } ], "metadata": { "celltoolbar": "Slideshow", "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.0" }, "rise": { "enable_chalkboard": true, "footer": "

Fachbereich Elektrotechnik und Informationstechnik (etit) | Institut für Teilchenbeschleunigung und Elektromagnetische Felder (TEMF) | Dr. Adrian Oeftiger

", "header": "", "scroll": true, "theme": "simple", "transition": "none" } }, "nbformat": 4, "nbformat_minor": 5 }