Nov 11 -  Nonlinear Motion in "State Space".ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "defbb6ee",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Nonlinear Motion in \"State Space\"\n",
    "\n",
    "# Continue with the damped, driven pendulum. Use the following block to set the problem each time:\n",
    "\n",
    "import numpy as np\n",
    "from scipy.integrate import solve_ivp\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.animation import FuncAnimation\n",
    "\n",
    "# Define the equation\n",
    "def eqn(t, phi, phi_dot):\n",
    "    beta = omega_0 / 4\n",
    "    omega = 2 * np.pi\n",
    "    gamma = 0.6\n",
    "    return phi_dot, -2 * beta * phi_dot - omega_0**2 * np.sin(phi) + gamma * omega_0**2 * np.cos(omega * t)\n",
    "\n",
    "# Set the parameters\n",
    "omega_0 = 1.5 * omega\n",
    "t_max = 40\n",
    "len = 1\n",
    "\n",
    "# Solve the equation\n",
    "sol = solve_ivp(eqn, [0, t_max], [-np.pi/2, 0], method='RK45', dense_output=True)\n",
    "\n",
    "# Calculate the position of the pendulum bob\n",
    "t = np.linspace(0, t_max, 500)\n",
    "phi = sol.sol(t)[0]\n",
    "x_bob = len * np.sin(phi)\n",
    "y_bob = -len * np.cos(phi)\n",
    "\n",
    "# Create an animation\n",
    "fig, ax = plt.subplots()\n",
    "line, = ax.plot([], [], lw=2)\n",
    "ax.set_xlim(-len, len)\n",
    "ax.set_ylim(-len, 0)\n",
    "ax.set_aspect('equal')\n",
    "ax.set_xlabel('x')\n",
    "ax.set_ylabel('y')\n",
    "\n",
    "def animate(frame):\n",
    "    line.set_data([0, x_bob[frame]], [0, y_bob[frame]])\n",
    "    return line,\n",
    "\n",
    "animation = FuncAnimation(fig, animate, frames=len(t), blit=True)\n",
    "\n",
    "# Uncomment the following line to save the animation as a video file\n",
    "# animation.save('pendulum_animation.mp4', writer='ffmpeg')\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1a7f3d80",
   "metadata": {},
   "outputs": [],
   "source": [
    "# The solution: Animated\n",
    "\n",
    "\n",
    "# Set the parameters\n",
    "len = 1\n",
    "t_max = 40\n",
    "\n",
    "# Create a figure and axes\n",
    "fig, ax = plt.subplots()\n",
    "ax.set_aspect('equal')\n",
    "ax.set_xlim(-len, len)\n",
    "ax.set_ylim(-len, 0)\n",
    "line, = ax.plot([], [], lw=2)\n",
    "dot, = ax.plot([], [], 'ko')\n",
    "\n",
    "# Animation update function\n",
    "def update(time):\n",
    "    x = len * np.sin(sol.sol(time)[0])\n",
    "    y = -len * np.cos(sol.sol(time)[0])\n",
    "    line.set_data([0, x], [0, y])\n",
    "    dot.set_data(x, y)\n",
    "    return line, dot\n",
    "\n",
    "# Create the animation\n",
    "animation = FuncAnimation(fig, update, frames=np.linspace(0, t_max, 500),\n",
    "                          interval=1000/30, blit=True, repeat=False)\n",
    "\n",
    "# Uncomment the following line to save the animation as a video file\n",
    "# animation.save('pendulum_animation.mp4', writer='ffmpeg')\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c51ea888",
   "metadata": {},
   "outputs": [],
   "source": [
    "# The solution: Plot the angle\n",
    "\n",
    "\n",
    "\n",
    "# Set the parameters\n",
    "t_max = 40\n",
    "\n",
    "# Generate values for t\n",
    "t = np.linspace(0, t_max, 1000)\n",
    "\n",
    "# Evaluate \\[Phi][t] using the solution\n",
    "phi = sol.sol(t)[0]\n",
    "\n",
    "# Create the plot\n",
    "plt.plot(t, phi, lw=2)\n",
    "\n",
    "# Set plot labels and title\n",
    "plt.xlabel('t')\n",
    "plt.ylabel('\\[Phi]')\n",
    "plt.title('Plot of \\[Phi][t]')\n",
    "\n",
    "# Show the plot\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f9e0174d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# The solution: Plot the angular velocity\n",
    "\n",
    "\n",
    "# Set the parameters\n",
    "t_max = 40\n",
    "\n",
    "# Generate values for t\n",
    "t = np.linspace(0, t_max, 1000)\n",
    "\n",
    "# Evaluate \\[Phi]'[t] using the solution\n",
    "phi_prime = sol.sol(t)[1]\n",
    "\n",
    "# Create the plot\n",
    "plt.plot(t, phi_prime, lw=2)\n",
    "\n",
    "# Set plot labels and title\n",
    "plt.xlabel('t')\n",
    "plt.ylabel('\\[Phi]\\'')\n",
    "plt.title('Plot of \\[Phi]\\'[t]')\n",
    "\n",
    "# Show the plot\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20d8eb7e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# The solution: Plot the State Space\n",
    "\n",
    "\n",
    "# Set the parameters\n",
    "t_max = 40\n",
    "omega = 2 * np.pi\n",
    "\n",
    "# Generate values for t\n",
    "t = np.linspace(5, t_max, 1000)\n",
    "\n",
    "# Evaluate \\[Phi][t] and \\[Phi]'[t] using the solution\n",
    "phi = sol.sol(t)[0]\n",
    "phi_prime = sol.sol(t)[1]\n",
    "\n",
    "# Create the parametric plot\n",
    "plt.plot(phi, phi_prime, lw=2)\n",
    "\n",
    "# Set plot limits and aspect ratio\n",
    "plt.xlim(-1.1 * np.pi, 1.1 * np.pi)\n",
    "plt.ylim(-1.2 * omega * np.pi, 1.2 * omega * np.pi)\n",
    "plt.gca().set_aspect(1)\n",
    "\n",
    "# Set plot labels and title\n",
    "plt.xlabel('\\[Phi]')\n",
    "plt.ylabel('\\[Phi]\\'')\n",
    "plt.title('Parametric Plot of [\\[Phi][t], \\[Phi]\\'[t]]')\n",
    "\n",
    "# Show the plot\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "155e6a38",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "These are various combinations to try\n",
    "\n",
    "(1) \\[Gamma]=\\[Beta]=0 and small \\[Phi][0], say \\[Pi]/10: Elliptical state space orbits!\n",
    "(2) Same but with large \\[Phi][0], say 9\\[Pi]/10\n",
    "(3) Add damping  \\[Beta]=\\[Omega]0/4\n",
    "(4) Turn on driving with \\[Gamma]=0.2; reset \\[Phi][0]=0\n",
    "(5) Drive more strongly with \\[Gamma]=0.9\n",
    "(5a) Adjust tMin to 5 on state space plot: Period=1!\n",
    "\n",
    "Now run through our \"standard\" period doubling cascade\n",
    "(6a) \\[Gamma]=1.06 and \\[Phi][0]=-\\[Pi]/2\n",
    "(6a) \\[Gamma]=1.078 and \\[Phi][0]=-\\[Pi]/2\n",
    "(6a) \\[Gamma]=1.081 and \\[Phi][0]=-\\[Pi]/2\n",
    "(6a) \\[Gamma]=1.0826 and \\[Phi][0]=-\\[Pi]/2 (run longer, blow up scale)\n",
    "\n",
    "Now chaotic driving strength\n",
    "(7) \\[Gamma]=1.09 and \\[Phi][0]=-\\[Pi]/2 (run longer, blow up scale)\n",
    "(8) \\[Gamma]=1.105 (Taylor Fig.12.24; run w/tMax=200)\n",
    "(9) \\[Gamma]=1.5 and \\[Beta]=\\[Omega]0/8 (Taylor Fig.12.27; run w/tMax=200)\n",
    "  Note: Try Mod[\\[Phi],2\\[Pi]] for \"rolling\" motion\n",
    "\"\"\""
   ]
  }
 ],
 "metadata": {
  "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.9.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}