inigial commit

eigenvectors_and_nodes.py

@@ -0,0 +1,47 @@
+# %% IMPORTS
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.linalg import eig
+from numpy import real
+# set global font size for plots:
+plt.rcParams.update({'font.size': 14})
+# %% EIGENVECTORS AND NODES
+# define systems:
+A_saddle   = np.array([[2, 0], [0, -2]])  # saddle point
+A_instable = np.array([[2, 0], [0, 2]])   # unstable node
+A_stable   = np.array([[-2, 0], [0, -2]]) # stable knot
+
+# calculate eigenvalues and eigenvectors:
+eigvals_saddle, eigvecs_saddle = eig(A_saddle)
+eigvals_instable, eigvecs_instable = eig(A_instable)
+eigvals_stable, eigvecs_stable = eig(A_stable)
+
+# plot the phase portraits for each system:
+def plot_phase_portrait(A, eigvecs, eigvals, title):
+    fig, ax = plt.subplots(figsize=(6, 6))
+    x, y = np.meshgrid(np.linspace(-5, 5, 20), np.linspace(-5, 5, 20))
+    u = A[0, 0] * x + A[0, 1] * y
+    v = A[1, 0] * x + A[1, 1] * y
+    ax.streamplot(x, y, u, v, color='b')
+    for vec in eigvecs.T:
+        # indicate the direction of the eigenvectors and consider the sign of the real part:
+        ax.quiver(0, 0, 3*real(vec[0]), 3*real(vec[1]), scale=1, 
+                  scale_units='xy', angles='xy', color='r')
+    ax.set_title(title+f"\neigenvalues (real part): {real(eigvals)}")
+    ax.axis('equal')
+    plt.tight_layout()
+    plt.savefig(f'figures/eigenvectors_and_nodes_{title}.png', dpi=120)
+    plt.show()
+
+# saddle point:
+plot_phase_portrait(A_saddle, eigvecs_saddle, eigvals_saddle, 'saddle point')
+
+# unstable node:
+plot_phase_portrait(A_instable, eigvecs_instable, eigvals_instable, 'unstable node')
+
+# stable knot:
+plot_phase_portrait(A_stable, eigvecs_stable, eigvals_stable, 'stable knot')
+
+
+
+# %%

phae_plane_analysis_pendulum.py

@@ -0,0 +1,79 @@
+# %% IMPORT
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.integrate import solve_ivp
+# set global font size for plots:
+plt.rcParams.update({'font.size': 12})
+# create a folder "figures" to save the plots (if it does not exist):
+import os
+if not os.path.exists('figures'):
+    os.makedirs('figures')
+# %% EXAMPLE: PENDULUM
+g = 9.81  # acceleration due to gravity in m/s^2
+l = 1.0   # length of the pendulum rope in meters
+
+# define the function for the phase portrait:
+def pendulum_phase_portrait(t, z):
+    theta, v = z  # theta is the angle, v is the angular velocity
+    dtheta_dt = v
+    dv_dt = -(g / l) * np.sin(theta)
+    return [dtheta_dt, dv_dt]
+
+# angles and speeds for phase portrait:
+theta_vals = np.linspace(-2 * np.pi, 2 * np.pi, 100)
+v_vals = np.linspace(-10, 10, 100)
+Theta, V = np.meshgrid(theta_vals, v_vals)
+U = Theta, V = np.meshgrid(theta_vals, v_vals)
+dTheta_dt, dV_dt = pendulum_phase_portrait(None, [Theta, V])
+
+# calculate the trajectory for the pendulum:
+Z0 = [1.0, 0] # initial conditions
+t_span = [0, 10] # time span
+t_eval = np.linspace(*t_span, 1000) # time points for the solution
+sol = solve_ivp(pendulum_phase_portrait, t_span, Z0, t_eval=t_eval) # solve the ODE
+
+# potential energy:
+theta_energy = np.linspace(-2*np.pi, 2*np.pi, 400)
+U_energy = l * (1 - np.cos(theta_energy))
+
+# plot potential energy:
+fig = plt.figure(1, figsize=(6, 4))
+plt.plot(theta_energy, U_energy, 'r-')
+plt.title('Potential energy of the pendulum')
+plt.xlabel('Theta (rad)')
+plt.ylabel('Potential energy (joules)')
+plt.tight_layout()
+plt.savefig('figures/pendulum_phase_energy.png', dpi=120)
+plt.show()
+
+# plot phase portrait:
+fig = plt.figure(1, figsize=(6, 4))
+speed = np.sqrt(dTheta_dt**2 + dV_dt**2)
+#ax[1].streamplot(Theta, V, dTheta_dt, dV_dt, density=1.5)
+plt.streamplot(Theta, V, dTheta_dt, dV_dt, color=speed, cmap='cool', density=1.5)
+# plot the trajectory:
+plt.plot(sol.y[0], sol.y[1], '-', label=f'Trajectory for $(\Theta_0,v_0)$=({Z0[0]},{Z0[1]})', c='r', lw=2)
+# plot nullclines:
+# the x-nullcline is just the y-axis since dt x=0 -> y=dt Theta=0:
+plt.axhline(0, color='red', linestyle='--', label='x-nullcline ($\dot{x}=0$, $y=0$)')
+# the y-nullcline is the set of points where dt y=0, i.e., x=n*pi for n in Z:
+for n in range(-2, 3):  # for n from -2 to 2 to cover the relevant points
+    plt.axvline(n * np.pi, color='blue', linestyle='--', label='y-nullcline ($\dot{y}=0$, $x=n\pi$)' if n == -2 else "")
+
+# calculate and plot the separatrix for points near the saddle points:
+epsilon = 1e-8  # small push
+for theta0 in [ -3.0*np.pi + epsilon, 3.0*np.pi - epsilon]:
+    sol_separatrix = solve_ivp(pendulum_phase_portrait, [-10, 10], [theta0, 0], 
+                               t_eval=np.linspace(-10, 10, 1000))
+    plt.plot(sol_separatrix.y[0], sol_separatrix.y[1], 'k--', lw=1.5, 
+                 label='Separatrix' if theta0 == -3.0*np.pi + epsilon else "")
+
+plt.xlim(-2*np.pi, 2*np.pi)
+plt.legend(frameon=True, loc='upper right', fontsize=10)
+plt.title('Phase portrait of the pendulum')
+plt.xlabel('Theta (rad)')
+plt.ylabel('Angular velocity (rad/s)')
+plt.tight_layout()
+plt.savefig(f'figures/pendulum_phase_portrait_z{Z0[0]}_{Z0[1]}.png', dpi=120)
+plt.show()
+# %% END