diff --git a/week3/config3.py b/week3/config3.py index 8b137891791fe96927ad78e64b0aad7bded08bdc..21a4848bc2e6b21b3c3494db9d91cb28e80932d2 100644 --- a/week3/config3.py +++ b/week3/config3.py @@ -1 +1,133 @@ +#!/usr/bin/env python +m = 1 +g = 1 +L = 1 + +dt = 0.1 + +def hamiltonian(theta, p): + T = p * p / (2 * m * L*L) + U = m * g * L * (1 - np.cos(theta)) + return T + U + +def solve_euler(theta, p, dt=0.1): # for comparison + theta_next = theta + dt * p / (m * L*L) + p_next = p - dt * m * g * L * np.sin(theta) + return (theta_next, p_next) + +def solve_leapfrog(theta, p, dt=dt): + theta_half = theta + dt / 2 * p / (m * L*L) + p_next = p - dt * m * g * L * np.sin(theta_half) + theta_next = theta_half + dt / 2 * p_next / (m * L*L) + return (theta_next, p_next) + +def emittance(theta, p): + N = len(theta) + + # subtract centroids + theta = theta - 1/N * np.sum(theta) + p = p - 1/N * np.sum(p) + + # compute Σ matrix entries + theta_sq = 1/N * np.sum(theta * theta) + p_sq = 1/N * np.sum(p * p) + crossterm = 1/N * np.sum(theta * p) + + # determinant of Σ matrix + epsilon = np.sqrt(theta_sq * p_sq - crossterm * crossterm) + return epsilon + +# more stuff: + +# for TUDa jupyter hub (https://tu-jupyter-i.ca.hrz.tu-darmstadt.de/) +# --> install dependencies via +# $ pip install -r requirements_noversions.txt --prefix=`pwd`/requirements +import sys +sys.path.append('./requirements/lib/python3.8/site-packages/') + +import warnings +warnings.filterwarnings('ignore') + +import numpy as np +np.random.seed(0) + +import matplotlib +import matplotlib.pyplot as plt + +import seaborn as sns +sns.set_context('talk', font_scale=1.2, rc={'lines.linewidth': 3}) +sns.set_style('ticks', + {'grid.linestyle': 'none', 'axes.edgecolor': '0', + 'axes.linewidth': 1.2, 'legend.frameon': True, + 'xtick.direction': 'out', 'ytick.direction': 'out', + 'xtick.top': True, 'ytick.right': True, + }) + +def set_axes(xlim=(-np.pi * 1.1, np.pi * 1.1), ylim=(-3, 3)): + plt.xlim(xlim) + plt.ylim(ylim) + plt.xlabel(r'$\theta$') + plt.ylabel(r'$p$') + plt.xticks([-np.pi, -np.pi/2, 0, np.pi/2, np.pi], [r"$-\pi$", r"$-\pi/2$", "0", r"$\pi/2$", r"$\pi$"]) + +def plot_hamiltonian(xlim=(-np.pi * 1.1, np.pi * 1.1), ylim=(-3, 3)): + TH, PP = np.meshgrid(np.linspace(*xlim, num=100), + np.linspace(*ylim, num=100)) + HH = hamiltonian(TH, PP) + plt.contourf(TH, PP, HH, cmap=plt.get_cmap('hot_r'), levels=12, zorder=0) + plt.colorbar(label=r'$\mathcal{H}(\theta, p)$') + set_axes(xlim, ylim) + +def plot_macro_evolution(results_thetas, results_ps): + n_steps, N = results_thetas.shape + + centroids_theta = 1/N * np.sum(results_thetas, axis=1) + centroids_p = 1/N * np.sum(results_ps, axis=1) + + var_theta = 1/N * np.sum(results_thetas * results_thetas, axis=1) + var_p = 1/N * np.sum(results_ps * results_ps, axis=1) + + results_emit = np.zeros(n_steps, dtype=np.float32) + for k in range(n_steps): + results_emit[k] = emittance(results_thetas[k], results_ps[k]) + + fig, ax = plt.subplots(1, 3, figsize=(12, 4)) + + plt.sca(ax[0]) + plt.plot(centroids_theta, label=r'$\langle\theta\rangle$') + plt.plot(centroids_p, label=r'$\langle p\rangle$') + plt.scatter([0, 0], [centroids_theta[0], centroids_p[0]], marker='*', c='k') + plt.xlabel('Steps $k$') + plt.ylabel('Centroid amplitude') + plt.legend() + + plt.sca(ax[1]) + plt.plot(var_theta, label=r'$\langle\theta^2\rangle$') + plt.plot(var_p, label=r'$\langle p^2\rangle$') + plt.scatter([0, 0], [var_theta[0], var_p[0]], marker='*', c='k') + plt.xlabel('Steps $k$') + plt.ylabel('Variance') + plt.legend() + + plt.sca(ax[2]) + plt.plot(results_emit) + plt.scatter([0], [results_emit[0]], marker='*', c='k') + plt.xlabel('Steps $k$') + plt.ylabel('RMS emittance $\epsilon$') + + plt.tight_layout() + + return ax + +def get_boundary_ids(thetas, ps): + psf = ps.flatten(); tsf = thetas.flatten() + ps_min = ps.min(); ps_max = ps.max(); thetas_min = thetas.min(); thetas_max = thetas.max() + + seq = np.arange(len(psf)) + i_right = seq[(psf == 0) * (tsf == thetas_max)] + i_bottom = seq[psf == ps_min] + i_left = seq[(psf == 0) * (tsf == thetas_min)] + i_top = seq[psf == ps_max] + + return np.concatenate((i_right, i_top, i_left, i_bottom, i_right)) diff --git a/week3/nmap3.ipynb b/week3/nmap3.ipynb new file mode 100644 index 0000000000000000000000000000000000000000..af0bfe29954571c0276cfb1764d77a52f502be87 --- /dev/null +++ b/week3/nmap3.ipynb @@ -0,0 +1,1283 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "incomplete-medline", + "metadata": {}, + "source": [ + "<h1 style=\"text-align: center; vertical-align: middle;\">Numerical Methods in Accelerator Physics</h1>\n", + "<h2 style=\"text-align: center; vertical-align: middle;\">Python examples -- Week 3</h2>" + ] + }, + { + "cell_type": "markdown", + "id": "atlantic-france", + "metadata": {}, + "source": [ + "<h2>Run this first!</h2>\n", + "\n", + "Imports and modules:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "weighted-binary", + "metadata": {}, + "outputs": [], + "source": [ + "from config3 import *\n", + "%matplotlib inline\n", + "import PyNAFF" + ] + }, + { + "cell_type": "markdown", + "id": "normal-gasoline", + "metadata": {}, + "source": [ + "<h3>Pendulum parameters, Hamiltonian</h3>" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "supported-aggregate", + "metadata": {}, + "outputs": [], + "source": [ + "m = 1 # point mass\n", + "g = 1 # magnitude of the gravitational field\n", + "L = 1 # length of the rod" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "unlimited-suspension", + "metadata": {}, + "outputs": [], + "source": [ + "### Values of hamiltonian for comparison with theory\n", + "TH, PP = np.meshgrid(np.linspace(-np.pi * 1.1, np.pi * 1.1, 100), \n", + " np.linspace(-3, 3, 100))\n", + "\n", + "HH = hamiltonian(TH, PP)" + ] + }, + { + "cell_type": "markdown", + "id": "other-avatar", + "metadata": {}, + "source": [ + "<h2>Statistics</h2>" + ] + }, + { + "cell_type": "markdown", + "id": "overhead-listing", + "metadata": {}, + "source": [ + "<h3>Collection of pendulums (slides 14 & 17)</h3>" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "enormous-section", + "metadata": {}, + "outputs": [], + "source": [ + "theta_grid = np.linspace(-0.5 * np.pi, 0.5 * np.pi, 21)\n", + "p_grid = np.linspace(-0.3, 0.3, 5)\n", + "\n", + "thetas, ps = np.meshgrid(theta_grid, p_grid)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "improved-failure", + "metadata": {}, + "outputs": [], + "source": [ + "N = len(theta_grid) * len(p_grid)\n", + "N" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "incredible-turkish", + "metadata": {}, + "outputs": [], + "source": [ + "plt.scatter(thetas, ps, c='b', marker='.')\n", + "\n", + "plot_hamiltonian()" + ] + }, + { + "cell_type": "markdown", + "id": "prompt-clearance", + "metadata": {}, + "source": [ + "<h3>Time evolution</h3>" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "signed-horizon", + "metadata": {}, + "outputs": [], + "source": [ + "n_steps = 100" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "japanese-taiwan", + "metadata": {}, + "outputs": [], + "source": [ + "results_thetas = np.zeros((n_steps, N), dtype=np.float32)\n", + "results_thetas[0] = thetas.flatten()\n", + "\n", + "results_ps = np.zeros((n_steps, N), dtype=np.float32)\n", + "results_ps[0] = ps.flatten()\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": "markdown", + "id": "commercial-chinese", + "metadata": {}, + "source": [ + "<h3>Observations of centroids (slide 19)</h3>" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "labeled-necklace", + "metadata": {}, + "outputs": [], + "source": [ + "centroids_theta = 1/N * np.sum(results_thetas, axis=1)\n", + "centroids_p = 1/N * np.sum(results_ps, axis=1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "broadband-genome", + "metadata": {}, + "outputs": [], + "source": [ + "plt.plot(centroids_theta, label=r'$\\langle\\theta\\rangle$')\n", + "plt.plot(centroids_p, label=r'$\\langle p\\rangle$')\n", + "\n", + "plt.xlabel('Steps $k$')\n", + "plt.ylabel('Centroid amplitude')\n", + "plt.legend(bbox_to_anchor=(1.05, 1));" + ] + }, + { + "cell_type": "markdown", + "id": "surface-grammar", + "metadata": {}, + "source": [ + "<h3>Observations of RMS size</h3>" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "micro-syntax", + "metadata": {}, + "outputs": [], + "source": [ + "var_theta = 1/N * np.sum(results_thetas * results_thetas, axis=1)\n", + "var_p = 1/N * np.sum(results_ps * results_ps, axis=1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "empirical-tiffany", + "metadata": {}, + "outputs": [], + "source": [ + "plt.plot(var_theta, label=r'$\\langle\\theta^2\\rangle$')\n", + "plt.plot(var_p, label=r'$\\langle p^2\\rangle$')\n", + "\n", + "plt.xlabel('Steps $k$')\n", + "plt.ylabel('Variance')\n", + "plt.legend(bbox_to_anchor=(1.05, 1));" + ] + }, + { + "cell_type": "markdown", + "id": "medical-hazard", + "metadata": {}, + "source": [ + "<h3>Distribution evolution over time</h3> " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "patient-viking", + "metadata": {}, + "outputs": [], + "source": [ + "k = 18\n", + "\n", + "plt.scatter(results_thetas[k], results_ps[k], c='b', marker='.')\n", + "\n", + "plot_hamiltonian()" + ] + }, + { + "cell_type": "markdown", + "id": "auburn-prompt", + "metadata": {}, + "source": [ + "<h3>RMS emittance evolution (slide 20) $$\\epsilon = \\sqrt{\\langle \\theta^2\\rangle \\langle p^2\\rangle - \\langle \\theta\\,p\\rangle^2}$$</h3>" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "painted-patrol", + "metadata": {}, + "outputs": [], + "source": [ + "def emittance(theta, p):\n", + " N = len(theta)\n", + " \n", + " # subtract centroids\n", + " theta = theta - 1/N * np.sum(theta)\n", + " p = p - 1/N * np.sum(p)\n", + " \n", + " # compute Σ matrix entries\n", + " theta_sq = 1/N * np.sum(theta * theta)\n", + " p_sq = 1/N * np.sum(p * p)\n", + " crossterm = 1/N * np.sum(theta * p)\n", + " \n", + " # determinant of Σ matrix\n", + " epsilon = np.sqrt(theta_sq * p_sq - crossterm * crossterm)\n", + " return epsilon" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "theoretical-maximum", + "metadata": {}, + "outputs": [], + "source": [ + "results_emit = np.zeros(n_steps, dtype=np.float32)\n", + "\n", + "for k in range(n_steps):\n", + " results_emit[k] = emittance(results_thetas[k], results_ps[k])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "portable-hudson", + "metadata": {}, + "outputs": [], + "source": [ + "plt.plot(results_emit)\n", + "\n", + "plt.xlabel('Steps $k$')\n", + "plt.ylabel('RMS emittance $\\epsilon$');" + ] + }, + { + "cell_type": "markdown", + "id": "sticky-insured", + "metadata": {}, + "source": [ + "<h3>Exercise: Linearized pendulum motion</h3>" + ] + }, + { + "cell_type": "markdown", + "id": "gentle-anthony", + "metadata": {}, + "source": [ + "$$\\left.\\begin{matrix}\\,\n", + " \\theta_{k+1/2} &=& \\theta_k + {\\cfrac{\\Delta t}{2}}\\cdot\\cfrac{p_{k}}{mL^2} \\\\\n", + " p_{k+1} &=& p_k - \\Delta t\\cdot m g L\\sin\\theta_{k+1/2} \\\\\n", + " \\theta_{k+1} &=& \\theta_{k+1/2} + {\\cfrac{\\Delta t}{2}}\\cdot\\cfrac{p_{k+1}}{mL^2}\n", + "\\end{matrix}\\right\\} \\quad \\stackrel{\\text{Taylor}}{\\Longrightarrow} \\quad \\left\\{\\begin{matrix}\\,\n", + " \\theta_{k+1/2} &=& \\theta_k + {\\cfrac{\\Delta t}{2}}\\cdot\\cfrac{p_{k}}{mL^2} \\\\\n", + " p_{k+1} &=& p_k - \\Delta t\\cdot m g L\\color{red}{\\left(\\theta_{k+1/2}-\\frac{\\theta_{k+1/2}^3}{3!}+\\frac{\\theta_{k+1/2}^5}{5!}-\\frac{\\theta_{k+1/2}^7}{7!}+...\\right)} \\\\\n", + " \\theta_{k+1} &=& \\theta_{k+1/2} + {\\cfrac{\\Delta t}{2}}\\cdot\\cfrac{p_{k+1}}{mL^2}\n", + "\\end{matrix}\\right.$$ " + ] + }, + { + "cell_type": "markdown", + "id": "challenging-withdrawal", + "metadata": {}, + "source": [ + "$$\\left.\\begin{matrix}\\,\n", + " \\theta_{k+1/2} &=& \\theta_k + {\\cfrac{\\Delta t}{2}}\\cdot\\cfrac{p_{k}}{mL^2} \\\\\n", + " p_{k+1} &=& p_k - \\Delta t\\cdot m g L\\sin\\theta_{k+1/2} \\\\\n", + " \\theta_{k+1} &=& \\theta_{k+1/2} + {\\cfrac{\\Delta t}{2}}\\cdot\\cfrac{p_{k+1}}{mL^2}\n", + "\\end{matrix}\\right\\} \\quad \\stackrel{\\text{First order}}{\\Longrightarrow} \\quad \\left\\{\\begin{matrix}\\,\n", + " \\theta_{k+1/2} &=& \\theta_k + {\\cfrac{\\Delta t}{2}}\\cdot\\cfrac{p_{k}}{mL^2} \\\\\n", + " p_{k+1} &=& p_k - \\Delta t\\cdot m g L\\cdot\\color{red}{\\theta_{k+1/2}} \\\\\n", + " \\theta_{k+1} &=& \\theta_{k+1/2} + {\\cfrac{\\Delta t}{2}}\\cdot\\cfrac{p_{k+1}}{mL^2}\n", + "\\end{matrix}\\right.$$ " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "polish-runner", + "metadata": {}, + "outputs": [], + "source": [ + "def solve_leapfrog_linear(theta, p, dt=dt):\n", + " theta_half = theta + dt / 2 * p / (m * L*L)\n", + " p_next = p - dt * m * g * L * theta_half\n", + " theta_next = theta_half + dt / 2 * p_next / (m * L*L)\n", + " return (theta_next, p_next)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "north-verification", + "metadata": {}, + "outputs": [], + "source": [ + "results_thetas_linear = np.zeros((n_steps, N), dtype=np.float32)\n", + "results_thetas_linear[0] = thetas.flatten()\n", + "\n", + "results_ps_linear = np.zeros((n_steps, N), dtype=np.float32)\n", + "results_ps_linear[0] = ps.flatten()\n", + "\n", + "for k in range(1, n_steps):\n", + " results_thetas_linear[k], results_ps_linear[k] = solve_leapfrog_linear(results_thetas_linear[k - 1], results_ps_linear[k - 1])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "abandoned-gibson", + "metadata": {}, + "outputs": [], + "source": [ + "var_theta_linear = 1/N * np.sum(results_thetas_linear * results_thetas_linear, axis=1)\n", + "var_p_linear = 1/N * np.sum(results_ps_linear * results_ps_linear, axis=1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "informative-draft", + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "plt.plot(var_theta_linear, label=r'$\\langle\\theta^2\\rangle$')\n", + "plt.plot(var_p_linear, label=r'$\\langle p^2\\rangle$')\n", + "\n", + "plt.xlabel('Steps $k$')\n", + "plt.ylabel('Variance')\n", + "plt.legend(bbox_to_anchor=(1.05, 1));" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "positive-covering", + "metadata": {}, + "outputs": [], + "source": [ + "results_emit_linear = np.zeros(n_steps, dtype=np.float32)\n", + "\n", + "for k in range(n_steps):\n", + " results_emit_linear[k] = emittance(results_thetas_linear[k], results_ps_linear[k])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "furnished-birmingham", + "metadata": {}, + "outputs": [], + "source": [ + "plt.plot(results_emit_linear)\n", + "\n", + "plt.xlabel('Steps $k$')\n", + "plt.ylabel('RMS emittance $\\epsilon$');" + ] + }, + { + "cell_type": "markdown", + "id": "ongoing-springfield", + "metadata": {}, + "source": [ + "<h3>Centroid offset in (original) pendulum</h3>" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "analyzed-mechanism", + "metadata": {}, + "outputs": [], + "source": [ + "results_thetas2 = np.zeros((n_steps, N), dtype=np.float32)\n", + "results_thetas2[0] = thetas.flatten() + 0.1 * np.pi\n", + "\n", + "results_ps2 = np.zeros((n_steps, N), dtype=np.float32)\n", + "results_ps2[0] = ps.flatten()\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": "desirable-former", + "metadata": {}, + "outputs": [], + "source": [ + "centroids_theta2 = 1/N * np.sum(results_thetas2, axis=1)\n", + "centroids_p2 = 1/N * np.sum(results_ps2, axis=1)\n", + "\n", + "plt.plot(centroids_theta, label=r'no offset: $\\langle\\theta\\rangle$')\n", + "plt.plot(centroids_theta2, label=r'with offset: $\\langle\\theta\\rangle$')\n", + "plt.plot(centroids_p2, label=r'with offset: $\\langle p\\rangle$')\n", + "\n", + "plt.xlabel('Steps $k$')\n", + "plt.ylabel('Centroid amplitude')\n", + "plt.legend(bbox_to_anchor=(1.05, 1));" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "urban-business", + "metadata": {}, + "outputs": [], + "source": [ + "results_emit2 = np.zeros(n_steps, dtype=np.float32)\n", + "\n", + "for k in range(n_steps):\n", + " results_emit2[k] = emittance(results_thetas2[k], results_ps2[k])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "complex-welcome", + "metadata": {}, + "outputs": [], + "source": [ + "plt.plot(results_emit, label='no offset')\n", + "plt.plot(results_emit2, label='with offset')\n", + "\n", + "plt.xlabel('Steps $k$')\n", + "plt.ylabel('RMS emittance $\\epsilon$')\n", + "plt.legend();" + ] + }, + { + "cell_type": "markdown", + "id": "subtle-april", + "metadata": {}, + "source": [ + "<h3>Long-term behaviour</h3>" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "forty-lotus", + "metadata": {}, + "outputs": [], + "source": [ + "n_steps = 3000" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "linear-roman", + "metadata": {}, + "outputs": [], + "source": [ + "results_thetas3 = np.zeros((n_steps, N), dtype=np.float32)\n", + "results_thetas3[0] = thetas.flatten() + 0.1 * np.pi\n", + "\n", + "results_ps3 = np.zeros((n_steps, N), dtype=np.float32)\n", + "results_ps3[0] = ps.flatten()\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": "entire-transportation", + "metadata": {}, + "outputs": [], + "source": [ + "plot_macro_evolution(results_thetas3, results_ps3);" + ] + }, + { + "cell_type": "markdown", + "id": "completed-chapel", + "metadata": {}, + "source": [ + "<h2>Oscillation frequencies</h2>" + ] + }, + { + "cell_type": "markdown", + "id": "downtown-fleece", + "metadata": {}, + "source": [ + "<h3>Employ FFT on pendulum motion (slide 31)</h3>" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "weighted-native", + "metadata": {}, + "outputs": [], + "source": [ + "## Observe two particles:\n", + "i1 = N // 2 - 10\n", + "i2 = N // 2\n", + "\n", + "plt.plot(results_thetas3[:, i1])\n", + "plt.plot(results_thetas3[:, i2])\n", + "plt.xlabel('Steps $k$')\n", + "plt.ylabel(r'$\\theta$');" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "upset-mixture", + "metadata": {}, + "outputs": [], + "source": [ + "## Obtain the FFT spectra for both:\n", + "spec1 = np.fft.rfft(results_thetas3[:, i1])\n", + "spec2 = np.fft.rfft(results_thetas3[:, i2])\n", + "\n", + "freq = np.fft.rfftfreq(n_steps)\n", + "\n", + "## Frequency of linearized system\n", + "freq_theory = np.sqrt(1 / 1) * 0.1 / (2 * np.pi)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "passing-locking", + "metadata": {}, + "outputs": [], + "source": [ + "plt.plot(freq, np.abs(spec1))\n", + "plt.plot(freq, np.abs(spec2))\n", + "plt.axvline(freq_theory, c='k', zorder=0, ls='--')\n", + "plt.xlim(0, 0.05)\n", + "plt.xlabel('Phase advance per $\\Delta t$ [$2\\pi$]')\n", + "plt.ylabel('FFT spectrum')\n", + "plt.yscale('log');" + ] + }, + { + "cell_type": "markdown", + "id": "motivated-traffic", + "metadata": {}, + "source": [ + "<h3>Frequency vs Amplitude</h3>" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "refined-memory", + "metadata": {}, + "outputs": [], + "source": [ + "specs = np.abs(np.fft.rfft(results_thetas3.T))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "boxed-nation", + "metadata": {}, + "outputs": [], + "source": [ + "max_ids = np.argmax(specs, axis=1)\n", + "amplitudes = np.max(specs, axis=1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "framed-former", + "metadata": {}, + "outputs": [], + "source": [ + "i_range = np.where(ps.flatten() == 0)[0]\n", + "\n", + "plt.scatter(results_thetas3[0][i_range], freq[max_ids][i_range])\n", + "\n", + "plt.axhline(freq_theory, c='k', ls='--', zorder=0)\n", + "plt.xticks([-np.pi/2, 0, np.pi/2, np.pi], [r\"$-\\pi/2$\", \"0\", r\"$\\pi/2$\", r\"$\\pi$\"])\n", + "plt.xlabel(r'Initial $\\theta$ ($p=0$)')\n", + "plt.ylabel('Phase advance / \\n' + r'integration step $\\Delta t$ [$2\\pi$]');\n", + "plt.ylim(0.012, 0.0165);" + ] + }, + { + "cell_type": "markdown", + "id": "major-explosion", + "metadata": {}, + "source": [ + "<h3>NAFF Algorithm (slide 32)</h3>" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "forbidden-killing", + "metadata": {}, + "outputs": [], + "source": [ + "freqs_naff = []\n", + "\n", + "for signal in results_thetas3.T[i_range]:\n", + " freq_naff = PyNAFF.naff(signal, turns=n_steps - 1, nterms=1)\n", + " try:\n", + " freq_naff = freq_naff[0, 1]\n", + " except IndexError:\n", + " freq_naff = 0\n", + " freqs_naff += [freq_naff]\n", + "\n", + "freqs_naff = np.array(freqs_naff)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "secret-ghost", + "metadata": {}, + "outputs": [], + "source": [ + "plt.scatter(results_thetas3[0][i_range], freq[max_ids][i_range])\n", + "plt.plot(results_thetas3[0][i_range], freqs_naff, c='r', marker='.')\n", + "\n", + "plt.axhline(freq_theory, c='k', ls='--', zorder=0)\n", + "plt.xticks([-np.pi/2, 0, np.pi/2, np.pi], [r\"$-\\pi/2$\", \"0\", r\"$\\pi/2$\", r\"$\\pi$\"])\n", + "plt.xlabel(r'Initial $\\theta$ ($p=0$)')\n", + "plt.ylabel('Phase advance / \\n' + r'integration step $\\Delta t$ [$2\\pi$]');" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "revised-median", + "metadata": {}, + "outputs": [], + "source": [ + "plt.plot(results_thetas3[0][i_range], 100 * (freq[max_ids][i_range] - freqs_naff) / freqs_naff)\n", + "plt.xlabel(r'Initial $\\theta$ ($p=0$)')\n", + "plt.ylabel('Error (FFT - NAFF) / NAFF ' + r'[%]');" + ] + }, + { + "cell_type": "markdown", + "id": "curious-handbook", + "metadata": {}, + "source": [ + "<h2>Deterministic chaos</h2>" + ] + }, + { + "cell_type": "markdown", + "id": "changing-likelihood", + "metadata": {}, + "source": [ + "<h3>Discrete pendulum (slide 41)</h3>" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "demanding-aaron", + "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": "advance-combining", + "metadata": {}, + "outputs": [], + "source": [ + "plt.scatter(thetas, ps, c='b', marker='.')\n", + "\n", + "plot_hamiltonian()" + ] + }, + { + "cell_type": "markdown", + "id": "cathedral-feature", + "metadata": {}, + "source": [ + "<h3>Time evolution (using Leapfrog)</h3>" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "gothic-abuse", + "metadata": {}, + "outputs": [], + "source": [ + "n_steps = 1000" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "enhanced-simpson", + "metadata": {}, + "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": "moral-script", + "metadata": {}, + "outputs": [], + "source": [ + "plt.scatter(results_thetas, results_ps, c='b', marker='.', s=1)\n", + "\n", + "plot_hamiltonian()" + ] + }, + { + "cell_type": "markdown", + "id": "passing-rebate", + "metadata": {}, + "source": [ + "<h3>Chaos near unstable fixed point</h3>\n", + "\n", + "Let us investigate the unbounded, continuously rotating pendulum, with an energy just above the separatrix value:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "designing-papua", + "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": "regulated-wholesale", + "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": "mechanical-collection", + "metadata": {}, + "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": "optional-garbage", + "metadata": {}, + "source": [ + "<h3>Qualitative Investigation of Local Lyapunov Exponent</h3>\n", + "\n", + "First investigate around stable fixed point $\\theta=0$, $p=0$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "residential-naples", + "metadata": {}, + "outputs": [], + "source": [ + "dist = 1e-10\n", + "\n", + "thetas = 0 * 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": "running-trial", + "metadata": {}, + "outputs": [], + "source": [ + "n_steps = 100000" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "accepting-occurrence", + "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": "transsexual-insertion", + "metadata": {}, + "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.xlim(100000-1000, 100000)\n", + "\n", + "plt.xlabel('Steps $k$')\n", + "plt.ylabel(r'$\\theta$')\n", + "plt.legend();" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "silver-liquid", + "metadata": {}, + "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.xlim(100000-1000, 100000)\n", + "\n", + "plt.xlabel('Steps $k$')\n", + "plt.ylabel(r'$\\theta$')\n", + "plt.legend();" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "assigned-brighton", + "metadata": {}, + "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": "binary-subcommittee", + "metadata": {}, + "outputs": [], + "source": [ + "plt.plot(results_dist)\n", + "\n", + "#plt.xlim(100000-1000, 100000)\n", + "\n", + "plt.yscale('log')\n", + "plt.xlabel('Steps $k$')\n", + "plt.ylabel('Phase-space distance');" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "colonial-assignment", + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "plt.plot(results_dist)\n", + "\n", + "plt.xlim(100000-1000, 100000)\n", + "\n", + "plt.yscale('log')\n", + "plt.xlabel('Steps $k$')\n", + "plt.ylabel('Phase-space distance');" + ] + }, + { + "cell_type": "markdown", + "id": "molecular-flesh", + "metadata": {}, + "source": [ + "Let us now investigate around unstable fixed point $\\theta=\\pi$, $p=0$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "unknown-sentence", + "metadata": {}, + "outputs": [], + "source": [ + "thetas2 = np.pi * np.ones(N, dtype=np.float64) # change this line" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "opponent-fourth", + "metadata": {}, + "outputs": [], + "source": [ + "results_thetas4 = np.zeros((n_steps, N), dtype=np.float64)\n", + "results_thetas4[0] = thetas2\n", + "\n", + "results_ps4 = np.zeros((n_steps, N), dtype=np.float64)\n", + "results_ps4[0] = ps\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": "abandoned-immunology", + "metadata": {}, + "outputs": [], + "source": [ + "plt.plot(results_thetas4[:, 0], label='$p_{ini}$')\n", + "plt.plot(results_thetas4[:, 1], label='$p_{ini} + \\Delta p_{ini}$')\n", + "\n", + "#plt.xlim(100000-1000, 100000)\n", + "\n", + "plt.xlabel('Steps $k$')\n", + "plt.ylabel(r'$\\theta$')\n", + "plt.legend();" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fixed-restriction", + "metadata": {}, + "outputs": [], + "source": [ + "results_dist2 = np.sqrt(\n", + " (results_thetas4[:, 1] - results_thetas4[:, 0])**2 + \n", + " (results_ps4[:, 1] - results_ps4[:, 0])**2\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "express-candy", + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "plt.plot(results_dist2)\n", + "\n", + "#plt.xlim(100000-1000, 100000)\n", + "\n", + "plt.yscale('log')\n", + "plt.xlabel('Steps $k$')\n", + "plt.ylabel('Phase-space distance');" + ] + }, + { + "cell_type": "markdown", + "id": "permanent-production", + "metadata": {}, + "source": [ + "<h3>Evaluating local maximum Lyapunov exponent</h3>\n", + "\n", + "Nearly exponential increase over two periods (first 4000 steps), then bounded by system size.\n", + "\n", + "Local Maximum Lyapunov exponent estimated by simple linear regression:\n", + "\n", + "$$\\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)$$" + ] + }, + { + "cell_type": "markdown", + "id": "confirmed-battlefield", + "metadata": {}, + "source": [ + "Near stable fixed point:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "lightweight-hindu", + "metadata": {}, + "outputs": [], + "source": [ + "B, A = np.polyfit(\n", + " x=np.arange(n_steps)[:4000],\n", + " y=np.log(results_dist[:4000] / dist),\n", + " deg=1\n", + ")\n", + "\n", + "B" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "incorporated-gnome", + "metadata": {}, + "outputs": [], + "source": [ + "plt.plot(results_dist[:4000] / dist)\n", + "plt.plot(np.exp(B * np.arange(n_steps)[:4000] + A))\n", + "\n", + "plt.yscale('log');" + ] + }, + { + "cell_type": "markdown", + "id": "maritime-buffer", + "metadata": {}, + "source": [ + "Near unstable fixed point:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "surprised-traffic", + "metadata": {}, + "outputs": [], + "source": [ + "B, A = np.polyfit(\n", + " x=np.arange(n_steps)[:4000],\n", + " y=np.log(results_dist2[:4000] / dist),\n", + " deg=1\n", + ")\n", + "\n", + "B" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "earned-evaluation", + "metadata": {}, + "outputs": [], + "source": [ + "plt.plot(results_dist2[:4000] / dist)\n", + "plt.plot(np.exp(B * np.arange(n_steps)[:4000] + A))\n", + "\n", + "plt.yscale('log');" + ] + }, + { + "cell_type": "markdown", + "id": "casual-institution", + "metadata": {}, + "source": [ + "<h3>Frequency diffusion (discrete pendulum - slide 42)</h3>" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "pleasant-provider", + "metadata": {}, + "outputs": [], + "source": [ + "n_steps = 100000" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "exempt-trouble", + "metadata": {}, + "outputs": [], + "source": [ + "theta_ini = np.pi - 0.001\n", + "p_ini = 0" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "restricted-nicholas", + "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": "cooked-europe", + "metadata": {}, + "outputs": [], + "source": [ + "plt.plot(results_thetas4)\n", + "\n", + "plt.xlim(0, 1000)\n", + "\n", + "plt.xlabel('Steps $k$')\n", + "plt.ylabel(r'$\\theta$');" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "demonstrated-humanity", + "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": "vulnerable-donna", + "metadata": {}, + "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');" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "dominican-wales", + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure(figsize=(10, 5))\n", + "plt.hist(freqs_naff);" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "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.9.18" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}