Update notebook.ipynb

18383da1 · tegenolf · GitHub · 6a25d63b · 18383da1
Unverified Commit 18383da1 authored 8 months ago by tegenolf Committed by GitHub 8 months ago
--- a/week1/notebook.ipynb
+++ b/week1/notebook.ipynb
@@ -47,7 +47,7 @@
   "source": [
    "m = 1 # point mass\n",
    "g = 1 # magnitude of the gravitational field\n",
-    "l = 1 # length of the rod"
+    "L = 1 # length of the rod"
   ]
  },
  {

 %% Cell type:markdown id:consistent-cornell tags:

 <h1 style="text-align: center; vertical-align: middle;">Numerical Methods in Accelerator Physics</h1>
 <h2 style="text-align: center; vertical-align: middle;">Python examples -- Week 1</h2>

 %% Cell type:markdown id:terminal-diabetes tags:

 <h2>Run this first!</h2>

 Imports and modules:

 %% Cell type:code id:helpful-slovakia tags:

 ``` python
 from config import *
 %matplotlib inline
 ```

 %% Cell type:markdown id:atlantic-finance tags:

 <h2>Phase space plots based on Lagragian (slide 21)</h2>

 %% Cell type:code id:excited-technical tags:

 ``` python
 m = 1 # point mass
 g = 1 # magnitude of the gravitational field
-l = 1 # length of the rod
+L = 1 # length of the rod
 ```

 %% Cell type:code id:expected-blood tags:

 ``` python
 def plot_arrow(theta, p, dt=1):
    # equations of motion:
    dtheta = p / (m * L*L) * dt
    dp = - m * g * L * np.sin(theta) * dt

    # plotting
    plt.scatter([theta], [p], s=50, c='r', marker='*', zorder=10)
    plt.annotate('',
                 xytext=(theta, p),
                 xy=(theta + dtheta, p + dp),
                 arrowprops=dict(facecolor='black', shrink=0.03))
 ```

 %% Cell type:code id:together-birthday tags:

 ``` python
 plot_arrow(1, 0)
 plot_arrow(0, -1)
 plot_arrow(-1, 0)
 plot_arrow(0, 1)
 #plot_arrow(1/1.41, 1/1.41)

 # increase momentum:
 #for i in np.arange(0.4, 2.5, 0.2):
 #    plot_arrow(0, i)

 # increase angle:
 #for i in np.arange(1, np.pi * 1.1, 0.2):
 #    plot_arrow(i, 0)

 set_axes()
 ```

 %% Cell type:markdown id:sunset-dependence tags:

 <h2>Hamiltonian of the pendulum (slide 25)</h2>

 %% Cell type:code id:practical-keeping tags:

 ``` python
 def hamiltonian(theta, p):
    T = p * p / (2 * m * L*L)
    U = m * g * L * (1 - np.cos(theta))
    return T + U
 ```

 %% Cell type:code id:accessory-perfume tags:

 ``` python
 TH, PP = np.meshgrid(np.linspace(-np.pi * 1.1, np.pi * 1.1, 100),
                     np.linspace(-3, 3, 100))

 HH = hamiltonian(TH, PP)

 plt.contour(TH, PP, HH, levels=10)
 plt.colorbar(label=r'$\mathcal{H}(\theta, p)$')

 set_axes()
 ```

 %% Cell type:code id:random-honor tags:

 ``` python
 TH, PP = np.meshgrid(np.linspace(-np.pi * 1.1, np.pi * 1.1, 100),
                     np.linspace(-3, 3, 100))

 HH = hamiltonian(TH, PP)

 plt.contourf(TH, PP, HH, cmap=plt.get_cmap('hot_r'))
 plt.colorbar(label=r'$\mathcal{H}(\theta, p)$')

 plot_arrow(2, 0)

 set_axes()
 ```

 %% Cell type:markdown id:raised-herald tags:

 <h2>Explicit Euler Method (slide 28)</h2>

 %% Cell type:code id:brilliant-scale tags:

 ``` python
 def solve_euler(theta, p, dt=0.1):
    theta_next = theta + dt * p / (m * L*L)
    p_next = p - dt * m * g * L * np.sin(theta)
    return (theta_next, p_next)
 ```

 %% Cell type:markdown id:cordless-submission tags:

 <h3>First example</h3>

 %% Cell type:code id:generous-carroll tags:

 ``` python
 theta_ini = -1.1
 p_ini = 0
 n_steps = 100
 ```

 %% Cell type:code id:infrared-emergency tags:

 ``` python
 results_euler = np.zeros((n_steps, 2), dtype=np.float32)
 results_euler[0] = (theta_ini, p_ini)

 for k in range(1, n_steps):
    results_euler[k] = solve_euler(*results_euler[k - 1])
 ```

 %% Cell type:code id:prescribed-digit tags:

 ``` python
 plt.contour(TH, PP, HH, levels=10, linestyles='--', linewidths=1)

 plt.plot(results_euler[:, 0], results_euler[:, 1], c='b')

 set_axes()
 ```

 %% Cell type:code id:ethical-property tags:

 ``` python
 plt.plot(
    hamiltonian(results_euler[:, 0], results_euler[:, 1]),
    c='b', label='Euler integrator')
 plt.axhline(hamiltonian(theta_ini, p_ini), c='r', label='theory')

 plt.xlabel('Steps $k$')
 plt.ylabel(r'$\mathcal{H}(\theta, p)$')
 plt.legend();
 ```

 %% Cell type:markdown id:spanish-option tags:

 <h3>Second example: even worse!</h3>

 %% Cell type:code id:martial-ebony tags:

 ``` python
 theta_ini2 = -3.
 p_ini2 = 0
 n_steps = 100
 ```

 %% Cell type:code id:elementary-float tags:

 ``` python
 results_euler2 = np.zeros((n_steps, 2), dtype=np.float32)
 results_euler2[0] = (theta_ini2, p_ini2)

 for k in range(1, n_steps):
    results_euler2[k] = solve_euler(*results_euler2[k - 1])
 ```

 %% Cell type:code id:educational-group tags:

 ``` python
 plt.contour(TH, PP, HH, levels=10, linestyles='--', linewidths=1)
 plt.plot(results_euler2[:, 0], results_euler2[:, 1], c='b')
 set_axes()
 ```

 %% Cell type:code id:structured-discount tags:

 ``` python
 h_sep = hamiltonian(-np.pi, 0) # Value of Hamiltonian at separatrix
 h_sep
 ```

 %% Cell type:code id:decreased-deficit tags:

 ``` python
 plt.plot(
    hamiltonian(results_euler2[:, 0], results_euler2[:, 1]),
    c='b', label='Euler integrator')
 plt.axhline(hamiltonian(theta_ini2, p_ini2), c='r', label='theory')
 plt.axhline(h_sep, c='k', label='$\mathcal{H}_{sep}$')

 plt.xlabel('Steps $k$')
 plt.ylabel(r'$\mathcal{H}(\theta, p)$')
 plt.legend(bbox_to_anchor=(1.05, 1))
 plt.ylim(top=2.1);
 ```

 %% Cell type:code id:identified-opportunity tags:

 ``` python
 ```