Nov 1 - Taylor Ch. 12.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ca987c1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Getting Started: The Driven Damped Pendulum\n",
    "\n",
    "# We will use Python for numerical solutions of the equation of motion for the driven damped pendulum.\n",
    "\n",
    "# We are following the example in Chapter 12 of Taylor \"Classical Mechanics\"\n",
    "\n",
    "# First, we will clear the working space before getting started.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "678abf12",
   "metadata": {},
   "outputs": [],
   "source": [
    "# The Differential Equation of Motion\n",
    "\n",
    "# Define the left and right sides of the equation separately.\n",
    "\n",
    "import sympy as sp\n",
    "\n",
    "t = sp.symbols('t')\n",
    "phi = sp.Function('phi')(t)\n",
    "beta = sp.symbols('beta')\n",
    "omega0 = sp.symbols('omega0')\n",
    "gamma = sp.symbols('gamma')\n",
    "omega = sp.symbols('omega')\n",
    "\n",
    "# Replace driving frequency, natural frequency, and damping parameters with the numbers used for our example.\n",
    "eqn_left = phi.diff(t, t) + 2 * beta * phi.diff(t) + omega0**2 * sp.sin(phi)\n",
    "force_fun = gamma * omega0**2 * sp.cos(omega * t)\n",
    "\n",
    "eqn = eqn_left.subs([(beta, omega0 / 4), (omega0, 1.5 * omega), (omega, 2 * sp.pi)]) - force_fun\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "49a4c1e8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# The Linearized Equation\n",
    "\n",
    "# It will be handy to have the linearized equation of motion.\n",
    "\n",
    "t = sp.symbols('t')\n",
    "phi = sp.Function('phi')(t)\n",
    "beta = sp.symbols('beta')\n",
    "omega0 = sp.symbols('omega0')\n",
    "gamma = sp.symbols('gamma')\n",
    "omega = sp.symbols('omega')\n",
    "A = sp.symbols('A')\n",
    "\n",
    "eqn_left_lin = phi.diff(t, t) + 2 * beta * phi.diff(t) + omega0**2 * phi\n",
    "eqn_lin = eqn_left_lin.subs([(beta, omega0 / 4), (omega0, 1.5 * omega), (omega, 2 * sp.pi)]) - force_fun\n",
    "\n",
    "soln_lin = sp.dsolve(eqn_lin, phi, ics={phi.subs(t, 0): A, phi.diff(t).subs(t, 0): 0})\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "403600aa",
   "metadata": {},
   "outputs": [
    {
     "ename": "SyntaxError",
     "evalue": "invalid syntax (1830672801.py, line 7)",
     "output_type": "error",
     "traceback": [
      "\u001b[0;36m  Cell \u001b[0;32mIn[6], line 7\u001b[0;36m\u001b[0m\n\u001b[0;31m    omega0_val = ###\u001b[0m\n\u001b[0m                 ^\u001b[0m\n\u001b[0;31mSyntaxError\u001b[0m\u001b[0;31m:\u001b[0m invalid syntax\n"
     ]
    }
   ],
   "source": [
    "# Plots for Small Driving Amplitude\n",
    "\n",
    "# Choose \\[Gamma]=0.2, a relatively small driving term. Also choose initial amplitude zero. This means that the pendulum will not be driven to large amplitudes, and the linearized equation ought to be a good approximation.\n",
    "\n",
    "# So, solve the nonlinear equation numerically, and plot it against the linearized solution.\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.integrate import solve_ivp\n",
    "\n",
    "t_max = 10\n",
    "gamma_val = 0.2\n",
    "omega0_val = ###\n",
    "beta_val = omega0_val / 4\n",
    "\n",
    "eqn_non_lin = eqn.subs(gamma, gamma_val)\n",
    "\n",
    "soln_non_lin = solve_ivp(lambda t, phi: [phi[1], -2 * beta_val * phi[1] - omega0_val**2 * np.sin(phi[0]) + gamma_val * omega0_val**2 * np.cos(omega_val * t)], [0, t_max], [0, 0], method='RK45', t_eval=np.linspace(0, t_max, 500))\n",
    "\n",
    "plt.plot(soln_non_lin.t, soln_non_lin.y[0], label='Nonlinear')\n",
    "plt.plot(soln_non_lin.t, soln_lin.subs([(A, 0), (gamma, gamma_val), (omega, 2 * np.pi)]).rhs.subs(t, soln_non_lin.t), label='Linearized')\n",
    "\n",
    "plt.xlabel('t')\n",
    "plt.ylabel('phi(t)')\n",
    "plt.title('Comparison of Nonlinear and Linearized Solutions')\n",
    "plt.legend()\n",
    "plt.grid(True)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "8f9ee1fe",
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'soln_non_lin' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[10], line 28\u001b[0m\n\u001b[1;32m     25\u001b[0m     ax\u001b[38;5;241m.\u001b[39mset_ylabel(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124my\u001b[39m\u001b[38;5;124m'\u001b[39m)\n\u001b[1;32m     26\u001b[0m     ax\u001b[38;5;241m.\u001b[39mset_title(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mPendulum Animation\u001b[39m\u001b[38;5;124m'\u001b[39m)\n\u001b[0;32m---> 28\u001b[0m ani \u001b[38;5;241m=\u001b[39m FuncAnimation(fig, update, frames\u001b[38;5;241m=\u001b[39m\u001b[38;5;28mlen\u001b[39m(\u001b[43msoln_non_lin\u001b[49m\u001b[38;5;241m.\u001b[39mt), interval\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m500\u001b[39m)\n\u001b[1;32m     30\u001b[0m plt\u001b[38;5;241m.\u001b[39mshow()\n",
      "\u001b[0;31mNameError\u001b[0m: name 'soln_non_lin' is not defined"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAGdCAYAAAAfTAk2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAf30lEQVR4nO3dfWzV5f3/8deRllPR9ohUWqoFijPcBE2khNIuFbdgKd7BZJEb7ZxxjM4oAjEC4gLBhAIzjJlyM2vdNHHAFHD8wQh1CGH2AEIAO6gkarmZ9IhFOKcTV+6u7x/8OD+PpxRw/bQ9b56P5PzR61yf0+v6BO2TTz/n4HPOOQEAABhyXXsvAAAAoLUROAAAwBwCBwAAmEPgAAAAcwgcAABgDoEDAADMIXAAAIA5BA4AADAnqb0X0B7Onz+vo0ePKjU1VT6fr72XAwAAroBzTo2NjcrKytJ117V8jeaaDJyjR48qOzu7vZcBAAB+gCNHjui2225rcc41GTipqamSLpygtLS0dl4NAAC4EpFIRNnZ2dGf4y25JgPn4q+l0tLSCBwAABLMldxewk3GAADAHAIHAACYQ+AAAABzCBwAAGAOgQMAAMwhcAAAgDkEDgAAMIfAAQAA5hA4AADAHAIHAACYQ+AAAABzCBwAAGAOgQMAAMwhcAAAgDkEDgAAMIfAAQAA5hA4AADAHAIHAACYQ+AAAABzCBwAAGAOgQMAAMwhcAAAgDkEDgAAMIfAAQAA5hA4AADAHAIHAACYQ+AAAABzCBwAAGAOgQMAAMwhcAAAgDkEDgAAMIfAAQAA5hA4AADAHAIHAACYQ+AAAABzCBwAAGAOgQMAAMwhcAAAgDkEDgAAMIfAAQAA5hA4AADAHAIHAACYQ+AAAABz2iRwli5dqpycHKWkpCg3N1dbt25tcf6WLVuUm5urlJQU9enTR8uXL7/k3JUrV8rn82n06NGtvGoAAJCoPA+cVatWacqUKZo1a5Z2796twsJCjRw5UocPH252fl1dne6//34VFhZq9+7devHFFzV58mStXr06bu6hQ4f0/PPPq7Cw0OttAACABOJzzjkvv0FeXp4GDRqkZcuWRcf69++v0aNHq6ysLG7+9OnTtW7dOtXW1kbHSktLtXfvXgWDwejYuXPnNGzYMD355JPaunWrTp48qffee++K1hSJRBQIBBQOh5WWlvbDNwcAANrM1fz89vQKzunTp7Vr1y4VFRXFjBcVFam6urrZY4LBYNz8ESNGaOfOnTpz5kx0bO7cubrlllv01FNPXXYdTU1NikQiMQ8AAGCXp4HT0NCgc+fOKSMjI2Y8IyNDoVCo2WNCoVCz88+ePauGhgZJ0ocffqjKykpVVFRc0TrKysoUCASij+zs7B+wGwAAkCja5CZjn88X87VzLm7scvMvjjc2Nurxxx9XRUWF0tPTr+j7z5w5U+FwOPo4cuTIVe4AAAAkkiQvXzw9PV2dOnWKu1pz7NixuKs0F2VmZjY7PykpSd26ddO+fft08OBBPfTQQ9Hnz58/L0lKSkrSgQMHdPvtt8cc7/f75ff7W2NLAAAgAXh6Badz587Kzc1VVVVVzHhVVZUKCgqaPSY/Pz9u/saNGzV48GAlJyerX79+qqmp0Z49e6KPhx9+WD/5yU+0Z88efv0EAAC8vYIjSdOmTVNJSYkGDx6s/Px8vfbaazp8+LBKS0slXfj10RdffKG33npL0oV3TJWXl2vatGmaOHGigsGgKisrtWLFCklSSkqKBg4cGPM9brrpJkmKGwcAANcmzwNn7NixOn78uObOnav6+noNHDhQ69evV69evSRJ9fX1MZ+Jk5OTo/Xr12vq1KlasmSJsrKy9Oqrr2rMmDFeLxUAABjh+efgdER8Dg4AAImnw3wODgAAQHsgcAAAgDkEDgAAMIfAAQAA5hA4AADAHAIHAACYQ+AAAABzCBwAAGAOgQMAAMwhcAAAgDkEDgAAMIfAAQAA5hA4AADAHAIHAACYQ+AAAABzCBwAAGAOgQMAAMwhcAAAgDkEDgAAMIfAAQAA5hA4AADAHAIHAACYQ+AAAABzCBwAAGAOgQMAAMwhcAAAgDkEDgAAMIfAAQAA5hA4AADAHAIHAACYQ+AAAABzCBwAAGAOgQMAAMwhcAAAgDkEDgAAMIfAAQAA5hA4AADAHAIHAACYQ+AAAABzCBwAAGAOgQMAAMwhcAAAgDkEDgAAMIfAAQAA5hA4AADAHAIHAACYQ+AAAABzCBwAAGAOgQMAAMwhcAAAgDkEDgAAMIfAAQAA5hA4AADAHAIHAACYQ+AAAABzCBwAAGBOmwTO0qVLlZOTo5SUFOXm5mrr1q0tzt+yZYtyc3OVkpKiPn36aPny5THPV1RUqLCwUF27dlXXrl01fPhw7dixw8stAACABOJ54KxatUpTpkzRrFmztHv3bhUWFmrkyJE6fPhws/Pr6up0//33q7CwULt379aLL76oyZMna/Xq1dE5mzdv1vjx4/XBBx8oGAyqZ8+eKioq0hdffOH1dgAAQALwOeecl98gLy9PgwYN0rJly6Jj/fv31+jRo1VWVhY3f/r06Vq3bp1qa2ujY6Wlpdq7d6+CwWCz3+PcuXPq2rWrysvL9Ytf/OKya4pEIgoEAgqHw0pLS/sBuwIAAG3tan5+e3oF5/Tp09q1a5eKiopixouKilRdXd3sMcFgMG7+iBEjtHPnTp05c6bZY06dOqUzZ87o5ptvbvb5pqYmRSKRmAcAALDL08BpaGjQuXPnlJGRETOekZGhUCjU7DGhUKjZ+WfPnlVDQ0Ozx8yYMUO33nqrhg8f3uzzZWVlCgQC0Ud2dvYP2A0AAEgUbXKTsc/ni/naORc3drn5zY1L0sKFC7VixQqtWbNGKSkpzb7ezJkzFQ6Ho48jR45c7RYAAEACSfLyxdPT09WpU6e4qzXHjh2Lu0pzUWZmZrPzk5KS1K1bt5jxV155RfPmzdP777+vu+6665Lr8Pv98vv9P3AXAAAg0Xh6Badz587Kzc1VVVVVzHhVVZUKCgqaPSY/Pz9u/saNGzV48GAlJydHx373u9/p5Zdf1oYNGzR48ODWXzwAAEhYnv+Katq0aXr99df1xhtvqLa2VlOnTtXhw4dVWloq6cKvj777zqfS0lIdOnRI06ZNU21trd544w1VVlbq+eefj85ZuHChXnrpJb3xxhvq3bu3QqGQQqGQ/vOf/3i9HQAAkAA8/RWVJI0dO1bHjx/X3LlzVV9fr4EDB2r9+vXq1auXJKm+vj7mM3FycnK0fv16TZ06VUuWLFFWVpZeffVVjRkzJjpn6dKlOn36tH7+85/HfK/Zs2drzpw5Xm8JAAB0cJ5/Dk5HxOfgAACQeDrM5+AAAAC0BwIHAACYQ+AAAABzCBwAAGAOgQMAAMwhcAAAgDkEDgAAMIfAAQAA5hA4AADAHAIHAACYQ+AAAABzCBwAAGAOgQMAAMwhcAAAgDkEDgAAMIfAAQAA5hA4AADAHAIHAACYQ+AAAABzCBwAAGAOgQMAAMwhcAAAgDkEDgAAMIfAAQAA5hA4AADAHAIHAACYQ+AAAABzCBwAAGAOgQMAAMwhcAAAgDkEDgAAMIfAAQAA5hA4AADAHAIHAACYQ+AAAABzCBwAAGAOgQMAAMwhcAAAgDkEDgAAMIfAAQAA5hA4AADAHAIHAACYQ+AAAABzCBwAAGAOgQMAAMwhcAAAgDkEDgAAMIfAAQAA5hA4AADAHAIHAACYQ+AAAABzCBwAAGAOgQMAAMwhcAAAgDkEDgAAMIfAAQAA5rRJ4CxdulQ5OTlKSUlRbm6utm7d2uL8LVu2KDc3VykpKerTp4+WL18eN2f16tUaMGCA/H6/BgwYoLVr13q1fAAAkGA8D5xVq1ZpypQpmjVrlnbv3q3CwkKNHDlShw8fbnZ+XV2d7r//fhUWFmr37t168cUXNXnyZK1evTo6JxgMauzYsSopKdHevXtVUlKiRx99VNu3b/d6OwAAIAH4nHPOy2+Ql5enQYMGadmyZdGx/v37a/To0SorK4ubP336dK1bt061tbXRsdLSUu3du1fBYFCSNHbsWEUiEf3973+PzikuLlbXrl21YsWKy64pEokoEAgoHA4rLS3tf9keAABoI1fz89vTKzinT5/Wrl27VFRUFDNeVFSk6urqZo8JBoNx80eMGKGdO3fqzJkzLc651Gs2NTUpEonEPAAAgF2eBk5DQ4POnTunjIyMmPGMjAyFQqFmjwmFQs3OP3v2rBoaGlqcc6nXLCsrUyAQiD6ys7N/6JYAAEACaJObjH0+X8zXzrm4scvN//741bzmzJkzFQ6Ho48jR45c1foBAEBiSfLyxdPT09WpU6e4KyvHjh2LuwJzUWZmZrPzk5KS1K1btxbnXOo1/X6//H7/D90GAABIMJ5ewencubNyc3NVVVUVM15VVaWCgoJmj8nPz4+bv3HjRg0ePFjJycktzrnUawIAgGuLp1dwJGnatGkqKSnR4MGDlZ+fr9dee02HDx9WaWmppAu/Pvriiy/01ltvSbrwjqny8nJNmzZNEydOVDAYVGVlZcy7o5577jndc889WrBggUaNGqW//e1vev/99/XPf/7T6+0AAIAE4HngjB07VsePH9fcuXNVX1+vgQMHav369erVq5ckqb6+PuYzcXJycrR+/XpNnTpVS5YsUVZWll599VWNGTMmOqegoEArV67USy+9pN/+9re6/fbbtWrVKuXl5Xm9HQAAkAA8/xycjojPwQEAIPF0mM/BAQAAaA8EDgAAMIfAAQAA5hA4AADAHAIHAACYQ+AAAABzCBwAAGAOgQMAAMwhcAAAgDkEDgAAMIfAAQAA5hA4AADAHAIHAACYQ+AAAABzCBwAAGAOgQMAAMwhcAAAgDkEDgAAMIfAAQAA5hA4AADAHAIHAACYQ+AAAABzCBwAAGAOgQMAAMwhcAAAgDkEDgAAMIfAAQAA5hA4AADAHAIHAACYQ+AAAABzCBwAAGAOgQMAAMwhcAAAgDkEDgAAMIfAAQAA5hA4AADAHAIHAACYQ+AAAABzCBwAAGAOgQMAAMwhcAAAgDkEDgAAMIfAAQAA5hA4AADAHAIHAACYQ+AAAABzCBwAAGAOgQMAAMwhcAAAgDkEDgAAMIfAAQAA5hA4AADAHAIHAACYQ+AAAABzCBwAAGAOgQMAAMzxNHBOnDihkpISBQIBBQIBlZSU6OTJky0e45zTnDlzlJWVpeuvv1733nuv9u3bF33+66+/1rPPPqu+ffuqS5cu6tmzpyZPnqxwOOzlVgAAQALxNHAmTJigPXv2aMOGDdqwYYP27NmjkpKSFo9ZuHChFi1apPLycn300UfKzMzUfffdp8bGRknS0aNHdfToUb3yyiuqqanRn//8Z23YsEFPPfWUl1sBAAAJxOecc168cG1trQYMGKBt27YpLy9PkrRt2zbl5+frk08+Ud++feOOcc4pKytLU6ZM0fTp0yVJTU1NysjI0IIFCzRp0qRmv9c777yjxx9/XN98842SkpIuu7ZIJKJAIKBwOKy0tLT/YZcAAKCtXM3Pb8+u4ASDQQUCgWjcSNLQoUMVCARUXV3d7DF1dXUKhUIqKiqKjvn9fg0bNuySx0iKbvRK4gYAANjnWRGEQiF17949brx79+4KhUKXPEaSMjIyYsYzMjJ06NChZo85fvy4Xn755Ute3ZEuXAVqamqKfh2JRC67fgAAkLiu+grOnDlz5PP5Wnzs3LlTkuTz+eKOd841O/5d33/+UsdEIhE98MADGjBggGbPnn3J1ysrK4ve6BwIBJSdnX0lWwUAAAnqqq/gPPPMMxo3blyLc3r37q2PP/5YX375ZdxzX331VdwVmosyMzMlXbiS06NHj+j4sWPH4o5pbGxUcXGxbrzxRq1du1bJycmXXM/MmTM1bdq06NeRSITIAQDAsKsOnPT0dKWnp192Xn5+vsLhsHbs2KEhQ4ZIkrZv365wOKyCgoJmj8nJyVFmZqaqqqp09913S5JOnz6tLVu2aMGCBdF5kUhEI0aMkN/v17p165SSktLiWvx+v/x+/5VuEQAAJDjPbjLu37+/iouLNXHiRG3btk3btm3TxIkT9eCDD8a8g6pfv35au3atpAu/mpoyZYrmzZuntWvX6l//+pd++ctfqkuXLpowYYKkC1duioqK9M0336iyslKRSEShUEihUEjnzp3zajsAACCBePq2o7fffluTJ0+Ovivq4YcfVnl5ecycAwcOxHxI3wsvvKBvv/1WTz/9tE6cOKG8vDxt3LhRqampkqRdu3Zp+/btkqQf/ehHMa9VV1en3r17e7gjAACQCDz7HJyOjM/BAQAg8XSIz8EBAABoLwQOAAAwh8ABAADmEDgAAMAcAgcAAJhD4AAAAHMIHAAAYA6BAwAAzCFwAACAOQQOAAAwh8ABAADmEDgAAMAcAgcAAJhD4AAAAHMIHAAAYA6BAwAAzCFwAACAOQQOAAAwh8ABAADmEDgAAMAcAgcAAJhD4AAAAHMIHAAAYA6BAwAAzCFwAACAOQQOAAAwh8ABAADmEDgAAMAcAgcAAJhD4AAAAHMIHAAAYA6BAwAAzCFwAACAOQQOAAAwh8ABAADmEDgAAMAcAgcAAJhD4AAAAHMIHAAAYA6BAwAAzCFwAACAOQQOAAAwh8ABAADmEDgAAMAcAgcAAJhD4AAAAHMIHAAAYA6BAwAAzCFwAACAOQQOAAAwh8ABAADmEDgAAMAcAgcAAJhD4AAAAHMIHAAAYA6BAwAAzCFwAACAOZ4GzokTJ1RSUqJAIKBAIKCSkhKdPHmyxWOcc5ozZ46ysrJ0/fXX695779W+ffsuOXfkyJHy+Xx67733Wn8DAAAgIXkaOBMmTNCePXu0YcMGbdiwQXv27FFJSUmLxyxcuFCLFi1SeXm5PvroI2VmZuq+++5TY2Nj3NzFixfL5/N5tXwAAJCgkrx64draWm3YsEHbtm1TXl6eJKmiokL5+fk6cOCA+vbtG3eMc06LFy/WrFmz9Mgjj0iS3nzzTWVkZOgvf/mLJk2aFJ27d+9eLVq0SB999JF69Ojh1TYAAEAC8uwKTjAYVCAQiMaNJA0dOlSBQEDV1dXNHlNXV6dQKKSioqLomN/v17Bhw2KOOXXqlMaPH6/y8nJlZmZedi1NTU2KRCIxDwAAYJdngRMKhdS9e/e48e7duysUCl3yGEnKyMiIGc/IyIg5ZurUqSooKNCoUaOuaC1lZWXR+4ACgYCys7OvdBsAACABXXXgzJkzRz6fr8XHzp07JanZ+2Occ5e9b+b7z3/3mHXr1mnTpk1avHjxFa955syZCofD0ceRI0eu+FgAAJB4rvoenGeeeUbjxo1rcU7v3r318ccf68svv4x77quvvoq7QnPRxV83hUKhmPtqjh07Fj1m06ZN+uyzz3TTTTfFHDtmzBgVFhZq8+bNca/r9/vl9/tbXDMAALDjqgMnPT1d6enpl52Xn5+vcDisHTt2aMiQIZKk7du3KxwOq6CgoNljcnJylJmZqaqqKt19992SpNOnT2vLli1asGCBJGnGjBn61a9+FXPcnXfeqd///vd66KGHrnY7AADAIM/eRdW/f38VFxdr4sSJ+uMf/yhJ+vWvf60HH3ww5h1U/fr1U1lZmX72s5/J5/NpypQpmjdvnu644w7dcccdmjdvnrp06aIJEyZIunCVp7kbi3v27KmcnByvtgMAABKIZ4EjSW+//bYmT54cfVfUww8/rPLy8pg5Bw4cUDgcjn79wgsv6Ntvv9XTTz+tEydOKC8vTxs3blRqaqqXSwUAAIb4nHOuvRfR1iKRiAKBgMLhsNLS0tp7OQAA4Apczc9v/i0qAABgDoEDAADMIXAAAIA5BA4AADCHwAEAAOYQOAAAwBwCBwAAmEPgAAAAcwgcAABgDoEDAADMIXAAAIA5BA4AADCHwAEAAOYQOAAAwBwCBwAAmEPgAAAAcwgcAABgDoEDAADMIXAAAIA5BA4AADCHwAEAAOYQOAAAwBwCBwAAmEPgAAAAcwgcAABgDoEDAADMIXAAAIA5BA4AADCHwAEAAOYQOAAAwBwCBwAAmEPgAAAAcwgcAABgDoEDAADMIXAAAIA5BA4AADCHwAEAAOYQOAAAwBwCBwAAmEPgAAAAcwgcAABgDoEDAADMSWrvBbQH55wkKRKJtPNKAADAlbr4c/viz/GWXJOB09jYKEnKzs5u55UAAICr1djYqEAg0OIcn7uSDDLm/PnzOnr0qFJTU+Xz+dp7Oe0uEokoOztbR44cUVpaWnsvxyzOc9vgPLcdznXb4Dz/f845NTY2KisrS9dd1/JdNtfkFZzrrrtOt912W3svo8NJS0u75v/jaQuc57bBeW47nOu2wXm+4HJXbi7iJmMAAGAOgQMAAMwhcCC/36/Zs2fL7/e391JM4zy3Dc5z2+Fctw3O8w9zTd5kDAAAbOMKDgAAMIfAAQAA5hA4AADAHAIHAACYQ+BcA06cOKGSkhIFAgEFAgGVlJTo5MmTLR7jnNOcOXOUlZWl66+/Xvfee6/27dt3ybkjR46Uz+fTe++91/obSBBenOevv/5azz77rPr27asuXbqoZ8+emjx5ssLhsMe76ViWLl2qnJwcpaSkKDc3V1u3bm1x/pYtW5Sbm6uUlBT16dNHy5cvj5uzevVqDRgwQH6/XwMGDNDatWu9Wn7CaO3zXFFRocLCQnXt2lVdu3bV8OHDtWPHDi+3kBC8+PN80cqVK+Xz+TR69OhWXnUCcjCvuLjYDRw40FVXV7vq6mo3cOBA9+CDD7Z4zPz5811qaqpbvXq1q6mpcWPHjnU9evRwkUgkbu6iRYvcyJEjnSS3du1aj3bR8Xlxnmtqatwjjzzi1q1b5z799FP3j3/8w91xxx1uzJgxbbGlDmHlypUuOTnZVVRUuP3797vnnnvO3XDDDe7QoUPNzv/8889dly5d3HPPPef279/vKioqXHJysnv33Xejc6qrq12nTp3cvHnzXG1trZs3b55LSkpy27Zta6ttdThenOcJEya4JUuWuN27d7va2lr35JNPukAg4P7973+31bY6HC/O80UHDx50t956qyssLHSjRo3yeCcdH4Fj3P79+52kmP9xB4NBJ8l98sknzR5z/vx5l5mZ6ebPnx8d++9//+sCgYBbvnx5zNw9e/a42267zdXX11/TgeP1ef6uv/71r65z587uzJkzrbeBDmzIkCGutLQ0Zqxfv35uxowZzc5/4YUXXL9+/WLGJk2a5IYOHRr9+tFHH3XFxcUxc0aMGOHGjRvXSqtOPF6c5+87e/asS01NdW+++eb/vuAE5dV5Pnv2rPvxj3/sXn/9dffEE08QOM45fkVlXDAYVCAQUF5eXnRs6NChCgQCqq6ubvaYuro6hUIhFRUVRcf8fr+GDRsWc8ypU6c0fvx4lZeXKzMz07tNJAAvz/P3hcNhpaWlKSnJ/j8ld/r0ae3atSvmHElSUVHRJc9RMBiMmz9ixAjt3LlTZ86caXFOS+fdMq/O8/edOnVKZ86c0c0339w6C08wXp7nuXPn6pZbbtFTTz3V+gtPUASOcaFQSN27d48b7969u0Kh0CWPkaSMjIyY8YyMjJhjpk6dqoKCAo0aNaoVV5yYvDzP33X8+HG9/PLLmjRp0v+44sTQ0NCgc+fOXdU5CoVCzc4/e/asGhoaWpxzqde0zqvz/H0zZszQrbfequHDh7fOwhOMV+f5ww8/VGVlpSoqKrxZeIIicBLUnDlz5PP5Wnzs3LlTkuTz+eKOd841O/5d33/+u8esW7dOmzZt0uLFi1tnQx1Ue5/n74pEInrggQc0YMAAzZ49+3/YVeK50nPU0vzvj1/ta14LvDjPFy1cuFArVqzQmjVrlJKS0gqrTVyteZ4bGxv1+OOPq6KiQunp6a2/2ARm/xq3Uc8884zGjRvX4pzevXvr448/1pdffhn33FdffRX3t4KLLv66KRQKqUePHtHxY8eORY/ZtGmTPvvsM910000xx44ZM0aFhYXavHnzVeym42rv83xRY2OjiouLdeONN2rt2rVKTk6+2q0kpPT0dHXq1Cnub7fNnaOLMjMzm52flJSkbt26tTjnUq9pnVfn+aJXXnlF8+bN0/vvv6+77rqrdRefQLw4z/v27dPBgwf10EMPRZ8/f/68JCkpKUkHDhzQ7bff3so7SRDtdO8P2sjFm1+3b98eHdu2bdsV3fy6YMGC6FhTU1PMza/19fWupqYm5iHJ/eEPf3Cff/65t5vqgLw6z845Fw6H3dChQ92wYcPcN998490mOqghQ4a43/zmNzFj/fv3b/GmzP79+8eMlZaWxt1kPHLkyJg5xcXF1/xNxq19np1zbuHChS4tLc0Fg8HWXXCCau3z/O2338b9v3jUqFHupz/9qaupqXFNTU3ebCQBEDjXgOLiYnfXXXe5YDDogsGgu/POO+Pevty3b1+3Zs2a6Nfz5893gUDArVmzxtXU1Ljx48df8m3iF+kafheVc96c50gk4vLy8tydd97pPv30U1dfXx99nD17tk33114uvq22srLS7d+/302ZMsXdcMMN7uDBg84552bMmOFKSkqi8y++rXbq1Klu//79rrKyMu5ttR9++KHr1KmTmz9/vqutrXXz58/nbeIenOcFCxa4zp07u3fffTfmz25jY2Ob76+j8OI8fx/vorqAwLkGHD9+3D322GMuNTXVpaamuscee8ydOHEiZo4k96c//Sn69fnz593s2bNdZmam8/v97p577nE1NTUtfp9rPXC8OM8ffPCBk9Tso66urm021gEsWbLE9erVy3Xu3NkNGjTIbdmyJfrcE0884YYNGxYzf/Pmze7uu+92nTt3dr1793bLli2Le8133nnH9e3b1yUnJ7t+/fq51atXe72NDq+1z3OvXr2a/bM7e/bsNthNx+XFn+fvInAu8Dn3/+5WAgAAMIJ3UQEAAHMIHAAAYA6BAwAAzCFwAACAOQQOAAAwh8ABAADmEDgAAMAcAgcAAJhD4AAAAHMIHAAAYA6BAwAAzCFwAACAOf8Ht4uZEzvoVekAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Animation of the Pendulum\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.animation import FuncAnimation\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.animation import FuncAnimation\n",
    "\n",
    "length = 1.0\n",
    "\n",
    "# Define soln_non_lin with the numerical solution\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "line, = ax.plot([], [], lw=2)\n",
    "\n",
    "def update(frame):\n",
    "    ax.cla()\n",
    "    ax.set_aspect('equal')\n",
    "    ax.plot([0, length * np.sin(soln_non_lin.t[frame])], [0, -length * np.cos(soln_non_lin.t[frame])], lw=2)\n",
    "    ax.set_xlim(-length, length)\n",
    "    ax.set_ylim(-length, length)\n",
    "    ax.set_xlabel('x')\n",
    "    ax.set_ylabel('y')\n",
    "    ax.set_title('Pendulum Animation')\n",
    "\n",
    "ani = FuncAnimation(fig, update, frames=len(soln_non_lin.t), interval=500)\n",
    "\n",
    "plt.show()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "2cb8bf90",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Same thing, but large initial amplitude\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Solve the nonlinear equation numerically\n",
    "soln_non_lin = solve_ivp(lambda t, phi: [phi[1], -2 * beta_val * phi[1] - omega0_val**2 * np.sin(phi[0]) + gamma_val * omega0_val**2 * np.cos(omega_val * t)], [0, t_max], [3 * np.pi / 4, 0], method='RK45', t_eval=np.linspace(0, t_max, 500))\n",
    "\n",
    "# Plot the solutions\n",
    "plt.plot(soln_non_lin.t, soln_non_lin.y[0], label='Nonlinear')\n",
    "plt.plot(soln_non_lin.t, soln_lin.subs([(A, 3 * np.pi / 4), (gamma, 0.2), (omega, 2 * np.pi)]).rhs.subs(t, soln_non_lin.t), label='Linearized')\n",
    "\n",
    "plt.xlabel('t')\n",
    "plt.ylabel('\\u03C6')\n",
    "plt.title('Comparison of Nonlinear and Linearized Solutions')\n",
    "plt.legend()\n",
    "plt.grid(True)\n",
    "plt.show()\n",
    "\n",
    "#They start out differently, but converge on the same driver after the transients die out."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "caa68c10",
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'beta_val' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[13], line 10\u001b[0m\n\u001b[1;32m      7\u001b[0m eqn_non_lin \u001b[38;5;241m=\u001b[39m eqn\u001b[38;5;241m.\u001b[39msubs(gamma, \u001b[38;5;241m0.9\u001b[39m)\n\u001b[1;32m      9\u001b[0m \u001b[38;5;66;03m# Solve the nonlinear equation numerically\u001b[39;00m\n\u001b[0;32m---> 10\u001b[0m soln_non_lin \u001b[38;5;241m=\u001b[39m \u001b[43msolve_ivp\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43;01mlambda\u001b[39;49;00m\u001b[43m \u001b[49m\u001b[43mt\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mphi\u001b[49m\u001b[43m:\u001b[49m\u001b[43m \u001b[49m\u001b[43m[\u001b[49m\u001b[43mphi\u001b[49m\u001b[43m[\u001b[49m\u001b[38;5;241;43m1\u001b[39;49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m-\u001b[39;49m\u001b[38;5;241;43m2\u001b[39;49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43m \u001b[49m\u001b[43mbeta_val\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43m \u001b[49m\u001b[43mphi\u001b[49m\u001b[43m[\u001b[49m\u001b[38;5;241;43m1\u001b[39;49m\u001b[43m]\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m-\u001b[39;49m\u001b[43m \u001b[49m\u001b[43momega0_val\u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m2\u001b[39;49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43m \u001b[49m\u001b[43mnp\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43msin\u001b[49m\u001b[43m(\u001b[49m\u001b[43mphi\u001b[49m\u001b[43m[\u001b[49m\u001b[38;5;241;43m0\u001b[39;49m\u001b[43m]\u001b[49m\u001b[43m)\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m+\u001b[39;49m\u001b[43m \u001b[49m\u001b[43mgamma_val\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43m \u001b[49m\u001b[43momega0_val\u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m2\u001b[39;49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43m \u001b[49m\u001b[43mnp\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mcos\u001b[49m\u001b[43m(\u001b[49m\u001b[43momega_val\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43m \u001b[49m\u001b[43mt\u001b[49m\u001b[43m)\u001b[49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43m[\u001b[49m\u001b[38;5;241;43m0\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mt_max\u001b[49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43m[\u001b[49m\u001b[38;5;241;43m0\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m0\u001b[39;49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mmethod\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43mRK45\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mt_eval\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mnp\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mlinspace\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;241;43m0\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mt_max\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m500\u001b[39;49m\u001b[43m)\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m     12\u001b[0m \u001b[38;5;66;03m# Plot the solutions\u001b[39;00m\n\u001b[1;32m     13\u001b[0m plt\u001b[38;5;241m.\u001b[39mplot(soln_non_lin\u001b[38;5;241m.\u001b[39mt, soln_non_lin\u001b[38;5;241m.\u001b[39my[\u001b[38;5;241m0\u001b[39m], label\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mNonlinear\u001b[39m\u001b[38;5;124m'\u001b[39m)\n",
      "File \u001b[0;32m/opt/conda/lib/python3.9/site-packages/scipy/integrate/_ivp/ivp.py:557\u001b[0m, in \u001b[0;36msolve_ivp\u001b[0;34m(fun, t_span, y0, method, t_eval, dense_output, events, vectorized, args, **options)\u001b[0m\n\u001b[1;32m    554\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m method \u001b[38;5;129;01min\u001b[39;00m METHODS:\n\u001b[1;32m    555\u001b[0m     method \u001b[38;5;241m=\u001b[39m METHODS[method]\n\u001b[0;32m--> 557\u001b[0m solver \u001b[38;5;241m=\u001b[39m \u001b[43mmethod\u001b[49m\u001b[43m(\u001b[49m\u001b[43mfun\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mt0\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43my0\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mtf\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mvectorized\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mvectorized\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43moptions\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    559\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m t_eval \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n\u001b[1;32m    560\u001b[0m     ts \u001b[38;5;241m=\u001b[39m [t0]\n",
      "File \u001b[0;32m/opt/conda/lib/python3.9/site-packages/scipy/integrate/_ivp/rk.py:94\u001b[0m, in \u001b[0;36mRungeKutta.__init__\u001b[0;34m(self, fun, t0, y0, t_bound, max_step, rtol, atol, vectorized, first_step, **extraneous)\u001b[0m\n\u001b[1;32m     92\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mmax_step \u001b[38;5;241m=\u001b[39m validate_max_step(max_step)\n\u001b[1;32m     93\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mrtol, \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39matol \u001b[38;5;241m=\u001b[39m validate_tol(rtol, atol, \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mn)\n\u001b[0;32m---> 94\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mf \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mfun\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mt\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43my\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m     95\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m first_step \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n\u001b[1;32m     96\u001b[0m     \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mh_abs \u001b[38;5;241m=\u001b[39m select_initial_step(\n\u001b[1;32m     97\u001b[0m         \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mfun, \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mt, \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39my, \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mf, \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mdirection,\n\u001b[1;32m     98\u001b[0m         \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39merror_estimator_order, \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mrtol, \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39matol)\n",
      "File \u001b[0;32m/opt/conda/lib/python3.9/site-packages/scipy/integrate/_ivp/base.py:138\u001b[0m, in \u001b[0;36mOdeSolver.__init__.<locals>.fun\u001b[0;34m(t, y)\u001b[0m\n\u001b[1;32m    136\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mfun\u001b[39m(t, y):\n\u001b[1;32m    137\u001b[0m     \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mnfev \u001b[38;5;241m+\u001b[39m\u001b[38;5;241m=\u001b[39m \u001b[38;5;241m1\u001b[39m\n\u001b[0;32m--> 138\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mfun_single\u001b[49m\u001b[43m(\u001b[49m\u001b[43mt\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43my\u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[0;32m/opt/conda/lib/python3.9/site-packages/scipy/integrate/_ivp/base.py:20\u001b[0m, in \u001b[0;36mcheck_arguments.<locals>.fun_wrapped\u001b[0;34m(t, y)\u001b[0m\n\u001b[1;32m     19\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mfun_wrapped\u001b[39m(t, y):\n\u001b[0;32m---> 20\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m np\u001b[38;5;241m.\u001b[39masarray(\u001b[43mfun\u001b[49m\u001b[43m(\u001b[49m\u001b[43mt\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43my\u001b[49m\u001b[43m)\u001b[49m, dtype\u001b[38;5;241m=\u001b[39mdtype)\n",
      "Cell \u001b[0;32mIn[13], line 10\u001b[0m, in \u001b[0;36m<lambda>\u001b[0;34m(t, phi)\u001b[0m\n\u001b[1;32m      7\u001b[0m eqn_non_lin \u001b[38;5;241m=\u001b[39m eqn\u001b[38;5;241m.\u001b[39msubs(gamma, \u001b[38;5;241m0.9\u001b[39m)\n\u001b[1;32m      9\u001b[0m \u001b[38;5;66;03m# Solve the nonlinear equation numerically\u001b[39;00m\n\u001b[0;32m---> 10\u001b[0m soln_non_lin \u001b[38;5;241m=\u001b[39m solve_ivp(\u001b[38;5;28;01mlambda\u001b[39;00m t, phi: [phi[\u001b[38;5;241m1\u001b[39m], \u001b[38;5;241m-\u001b[39m\u001b[38;5;241m2\u001b[39m \u001b[38;5;241m*\u001b[39m \u001b[43mbeta_val\u001b[49m \u001b[38;5;241m*\u001b[39m phi[\u001b[38;5;241m1\u001b[39m] \u001b[38;5;241m-\u001b[39m omega0_val\u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39m\u001b[38;5;241m2\u001b[39m \u001b[38;5;241m*\u001b[39m np\u001b[38;5;241m.\u001b[39msin(phi[\u001b[38;5;241m0\u001b[39m]) \u001b[38;5;241m+\u001b[39m gamma_val \u001b[38;5;241m*\u001b[39m omega0_val\u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39m\u001b[38;5;241m2\u001b[39m \u001b[38;5;241m*\u001b[39m np\u001b[38;5;241m.\u001b[39mcos(omega_val \u001b[38;5;241m*\u001b[39m t)], [\u001b[38;5;241m0\u001b[39m, t_max], [\u001b[38;5;241m0\u001b[39m, \u001b[38;5;241m0\u001b[39m], method\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mRK45\u001b[39m\u001b[38;5;124m'\u001b[39m, t_eval\u001b[38;5;241m=\u001b[39mnp\u001b[38;5;241m.\u001b[39mlinspace(\u001b[38;5;241m0\u001b[39m, t_max, \u001b[38;5;241m500\u001b[39m))\n\u001b[1;32m     12\u001b[0m \u001b[38;5;66;03m# Plot the solutions\u001b[39;00m\n\u001b[1;32m     13\u001b[0m plt\u001b[38;5;241m.\u001b[39mplot(soln_non_lin\u001b[38;5;241m.\u001b[39mt, soln_non_lin\u001b[38;5;241m.\u001b[39my[\u001b[38;5;241m0\u001b[39m], label\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mNonlinear\u001b[39m\u001b[38;5;124m'\u001b[39m)\n",
      "\u001b[0;31mNameError\u001b[0m: name 'beta_val' is not defined"
     ]
    }
   ],
   "source": [
    "# Same thing, but larger driving term\n",
    "\n",
    "# Re-set initial amplitude to zero.\n",
    "\n",
    "\n",
    "# Define the nonlinear equation with gamma = 0.9\n",
    "eqn_non_lin = eqn.subs(gamma, 0.9)\n",
    "\n",
    "# Solve the nonlinear equation numerically\n",
    "soln_non_lin = solve_ivp(lambda t, phi: [phi[1], -2 * beta_val * phi[1] - omega0_val**2 * np.sin(phi[0]) + gamma_val * omega0_val**2 * np.cos(omega_val * t)], [0, t_max], [0, 0], method='RK45', t_eval=np.linspace(0, t_max, 500))\n",
    "\n",
    "# Plot the solutions\n",
    "plt.plot(soln_non_lin.t, soln_non_lin.y[0], label='Nonlinear')\n",
    "plt.plot(soln_non_lin.t, soln_lin.subs([(A, 0), (gamma, 0.9)]).rhs.subs(t, soln_non_lin.t), label='Linearized')\n",
    "\n",
    "plt.xlabel('t')\n",
    "plt.ylabel('\\u03C6')\n",
    "plt.title('Comparison of Nonlinear and Linearized Solutions')\n",
    "plt.legend()\n",
    "plt.grid(True)\n",
    "plt.show()\n",
    "# OK, the linear approximation isn't so good anymore, but at least the response still looks sinusoidal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "b3da8167",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Animated...\n",
    "\n",
    "# import numpy as np\n",
    "# import matplotlib.pyplot as plt\n",
    "# from matplotlib.animation import FuncAnimation\n",
    "\n",
    "len_val = 1.0\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.set_xlim(-len_val, len_val)\n",
    "ax.set_ylim(-len_val, len_val)\n",
    "ax.set_aspect('equal')\n",
    "\n",
    "line, = ax.plot([], [], lw=2)\n",
    "\n",
    "def update(frame):\n",
    "    x = len_val * np.sin(soln_non_lin.y[0, frame])\n",
    "    y = -len_val * np.cos(soln_non_lin.y[0, frame])\n",
    "    line.set_data([0, x], [0, y])\n",
    "    return line,\n",
    "\n",
    "ani = FuncAnimation(fig, update, frames=len(soln_non_lin.t), interval=500)\n",
    "\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')\n",
    "plt.title('Pendulum Animation')\n",
    "plt.grid(True)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d9e73b55",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now, try driver with \\[Gamma] > 1\n",
    "\n",
    "# This means that the forcing function exceeds the weight of the pendulum bob\n",
    "\n",
    "# Define the equation of motion\n",
    "def equation(t, phi, beta, omega0, gamma, omega):\n",
    "    return phi[1], -2 * beta * phi[1] - omega0**2 * np.sin(phi[0]) + gamma * omega0**2 * np.cos(omega * t)\n",
    "\n",
    "# Parameters\n",
    "tMax = 40\n",
    "beta = 0.375\n",
    "omega0 = 1.5 * 2 * np.pi\n",
    "gamma = 1.06\n",
    "omega = 2 * np.pi\n",
    "\n",
    "# Nonlinear equation\n",
    "eqnNonLin = lambda t, phi: equation(t, phi, beta, omega0, gamma, omega)\n",
    "\n",
    "# Solve the nonlinear equation numerically\n",
    "solnNonLin = solve_ivp(eqnNonLin, [0, tMax], [0, 0], method='RK45', t_eval=np.linspace(0, tMax, 1000))\n",
    "\n",
    "# Plot the nonlinear solution\n",
    "plt.figure()\n",
    "plt.plot(solnNonLin.t, solnNonLin.y[0], label='Nonlinear')\n",
    "plt.xlabel('t')\n",
    "plt.ylabel('phi')\n",
    "plt.title('Nonlinear Solution')\n",
    "plt.legend()\n",
    "plt.show()\n",
    "\n",
    "# Linearized equation\n",
    "eqnLeftLin = lambda t, phi: equation(t, phi, beta, omega0, omega0/4, omega)\n",
    "eqnLin = lambda t, phi: eqnLeftLin(t, phi)\n",
    "\n",
    "# Solve the linearized equation analytically\n",
    "solnLin = np.array([np.pi/2 * np.sin(omega0 * t) for t in solnNonLin.t])\n",
    "\n",
    "# Plot the linearized solution\n",
    "plt.figure()\n",
    "plt.plot(solnNonLin.t, solnLin, label='Linearized')\n",
    "plt.xlabel('t')\n",
    "plt.ylabel('phi')\n",
    "plt.title('Linearized Solution')\n",
    "plt.legend()\n",
    "plt.show()\n",
    "\n",
    "# Well, all hell breaks loose. The nonlinear oscillator is a very different sort of thing when the driving amplitude gets \"large.\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "e825e020",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Animated...\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.animation import FuncAnimation\n",
    "\n",
    "# Define the length of the pendulum\n",
    "length = 1.0\n",
    "\n",
    "# Function to update the plot\n",
    "def update(time):\n",
    "    plt.clf()\n",
    "    plt.plot([0, length * np.sin(solnNonLin.y[0, time])], [0, -length * np.cos(solnNonLin.y[0, time])], 'b-')\n",
    "    plt.xlim(-length, length)\n",
    "    plt.ylim(-length, length)\n",
    "    plt.gca().set_aspect('equal')\n",
    "    plt.title('Pendulum Animation')\n",
    "\n",
    "# Create the figure and axis\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "# Animate the pendulum\n",
    "ani = FuncAnimation(fig, update, frames=len(solnNonLin.t), interval=500)\n",
    "\n",
    "# Display the animation\n",
    "plt.show()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "fe9a532f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# But the nonlinear oscillator is still periodic!\n",
    "\n",
    "# Wait long enough for the transients to die away.\n",
    "\n",
    "from scipy.integrate import solve_ivp\n",
    "\n",
    "# Define the equation of motion eqnNonLin\n",
    "def eqnNonLin(t, phi):\n",
    "    return phi[1], -2 * beta_val * phi[1] - omega0_val**2 * np.sin(phi[0]) + gamma_val * omega0_val**2 * np.cos(omega_val * t)\n",
    "\n",
    "# Define initial conditions\n",
    "phi0 = 0\n",
    "phi_dot0 = 0\n",
    "\n",
    "# Define the time range\n",
    "t_span = (0, 40)\n",
    "\n",
    "# Solve the differential equation\n",
    "solnNonLin = solve_ivp(eqnNonLin, t_span, [phi0, phi_dot0], method='RK45', dense_output=True)\n",
    "\n",
    "# Get the values at specific time points\n",
    "time_points = [30, 31, 32, 33, 34, 35, 36, 37, 38]\n",
    "values = solnNonLin.sol(time_points)[0]\n",
    "\n",
    "# Print the values\n",
    "for t, value in zip(time_points, values):\n",
    "    print(f\"phi({t}) = {value:.5f}\")\n",
    "\n",
    "\n",
    "# You can put in more digits if you want, and you'll have to wait longer to see it be periodic.\n",
    "\n",
    "# BUT, the point is that the period (=1) is still the same as the driving term."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c4cd1a68",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now make the driving a little bit larger\n",
    "\n",
    "# Don't bother plotting the linear solution anymore.\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Define the equation of motion eqnNonLin\n",
    "def eqnNonLin(t, phi):\n",
    "    return phi[1], -2 * beta_val * phi[1] - omega0_val**2 * np.sin(phi[0]) + gamma_val * omega0_val**2 * np.cos(omega_val * t)\n",
    "\n",
    "# Define initial conditions\n",
    "phi0 = 0\n",
    "phi_dot0 = 0\n",
    "\n",
    "# Define the time range\n",
    "t_span = (0, tMax)\n",
    "\n",
    "# Solve the differential equation\n",
    "solnNonLin = solve_ivp(eqnNonLin, t_span, [phi0, phi_dot0], method='RK45', dense_output=True)\n",
    "\n",
    "# Plot the solution\n",
    "t_values = np.linspace(0, tMax, 500)\n",
    "phi_values = solnNonLin.sol(t_values)[0]\n",
    "\n",
    "plt.plot(t_values, phi_values, linewidth=2)\n",
    "plt.xlabel('t')\n",
    "plt.ylabel('phi')\n",
    "plt.title('Nonlinear Pendulum Solution')\n",
    "plt.show()\n",
    "\n",
    "# Before running the code, make sure to define the values for beta_val, omega0_val, gamma_val, omega_val, and tMax according to your problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "688bcb9b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Is it still periodic at long times?\n",
    "\n",
    "# Define the equation of motion eqnNonLin\n",
    "def eqnNonLin(t, phi):\n",
    "    return phi[1], -2 * beta_val * phi[1] - omega0_val**2 * np.sin(phi[0]) + gamma_val * omega0_val**2 * np.cos(omega_val * t)\n",
    "\n",
    "# Define initial conditions\n",
    "phi0 = 0\n",
    "phi_dot0 = 0\n",
    "\n",
    "# Define the time range\n",
    "t_span = (0, 40)\n",
    "\n",
    "# Solve the differential equation\n",
    "solnNonLin = solve_ivp(eqnNonLin, t_span, [phi0, phi_dot0], method='RK45', dense_output=True)\n",
    "\n",
    "# Evaluate the solution at specific time points\n",
    "time_points = [30, 31, 32, 33, 34, 35, 36, 37, 38]\n",
    "phi_values = solnNonLin.sol(time_points)[0]\n",
    "\n",
    "# Print the phi values with number formatting\n",
    "print(np.around(phi_values, decimals=5))\n",
    "\n",
    "# Is the period now equal to two? Yes!\n",
    "\n",
    "# Plot the solution\n",
    "t_values = np.linspace(0, 40, 500)\n",
    "phi_values = solnNonLin.sol(t_values)[0]\n",
    "\n",
    "plt.plot(t_values, phi_values, linewidth=2)\n",
    "plt.xlabel('t')\n",
    "plt.ylabel('phi')\n",
    "plt.title('Nonlinear Pendulum Solution')\n",
    "plt.xlim(25, 40)\n",
    "plt.ylim(-8.7, -7.8)\n",
    "plt.show()\n",
    "\n",
    "\n",
    "# However, the result is dominated by the fundamental harmonic with period=1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a5e59649",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Once again, increase driver a little bit\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.integrate import solve_ivp\n",
    "\n",
    "# Define the equation of motion eqnNonLin\n",
    "def eqnNonLin(t, phi):\n",
    "    return phi[1], -2 * beta_val * phi[1] - omega0_val**2 * np.sin(phi[0]) + gamma_val * omega0_val**2 * np.cos(omega_val * t)\n",
    "\n",
    "# Define initial conditions\n",
    "phi0 = 0\n",
    "phi_dot0 = 0\n",
    "\n",
    "# Define the time range\n",
    "t_span = (0, 40)\n",
    "\n",
    "# Solve the differential equation\n",
    "solnNonLin = solve_ivp(eqnNonLin, t_span, [phi0, phi_dot0], method='RK45', dense_output=True)\n",
    "\n",
    "# Plot the solution\n",
    "t_values = np.linspace(0, tMax, 500)\n",
    "phi_values = solnNonLin.sol(t_values)[0]\n",
    "\n",
    "plt.plot(t_values, phi_values, linewidth=2)\n",
    "plt.xlabel('t')\n",
    "plt.ylabel('phi')\n",
    "plt.title('Nonlinear Pendulum Solution')\n",
    "plt.show()\n",
    "\n",
    "\n",
    "# Now (!) the dominant period=3, not 1. (Expand the scale if you want to see this more clearly, graphically.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "84269f87",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Evaluate to see the dominant period.\n",
    "# In Python, there is no direct equivalent to NumberForm in Mathematica. However, you can achieve similar formatting using the round function. \n",
    "\n",
    "indices = [30, 31, 32, 33, 34, 35, 36, 37, 38]\n",
    "phi_values = np.array([solnNonLin.sol(t)[0] for t in indices])\n",
    "\n",
    "formatted_phi_values = [round(phi, 6) for phi in phi_values]\n",
    "print(formatted_phi_values)\n",
    "\n",
    "# You can tabulate again, to see that the dominant period is 3, but also that the harmonics 1 and 2 are there, too.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d6ca25ad",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Something new: Change the initial conditions\n",
    "\n",
    "solnIC0 = solnNonLin\n",
    "\n",
    "tMax = 15\n",
    "\n",
    "solnNonLin = solve_ivp(lambda t, phi: [phi[1], -2 * beta_val * phi[1] - omega0_val**2 * np.sin(phi[0]) + gamma_val * omega0_val**2 * np.cos(omega_val * t)], [0, tMax], [-np.pi/2, 0], method='RK45', t_eval=np.linspace(0, tMax, 500))\n",
    "\n",
    "plt.plot(solnNonLin.t, solnNonLin.y[0], label='Solution')\n",
    "plt.plot(solnNonLin.t, solnIC0.y[0], label='Initial Condition')\n",
    "plt.xlabel('t')\n",
    "plt.ylabel('\\u03C6(t)')\n",
    "plt.title('Comparison of Solutions')\n",
    "plt.legend()\n",
    "plt.show()\n",
    "\n",
    "\n",
    "# A completely different \"attractor\"! The overall period is now=2. although harmonic period=1 still dominates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "91c04e54",
   "metadata": {},
   "outputs": [],
   "source": [
    "# A Series of Driving Strengths\n",
    "\n",
    "tMax = 40\n",
    "\n",
    "eqnNonLin = eqn.subs(gamma, 1.06)\n",
    "solnN1p06 = solve_ivp(lambda t, phi: [phi[1], -2 * beta_val * phi[1] - omega0_val**2 * np.sin(phi[0]) + gamma_val * omega0_val**2 * np.cos(omega_val * t)], [0, tMax], [-np.pi/2, 0], method='RK45', t_eval=np.linspace(0, tMax, 500))\n",
    "\n",
    "eqnNonLin = eqn.subs(gamma, 1.078)\n",
    "solnN1p078 = solve_ivp(lambda t, phi: [phi[1], -2 * beta_val * phi[1] - omega0_val**2 * np.sin(phi[0]) + gamma_val * omega0_val**2 * np.cos(omega_val * t)], [0, tMax], [-np.pi/2, 0], method='RK45', t_eval=np.linspace(0, tMax, 500))\n",
    "\n",
    "eqnNonLin = eqn.subs(gamma, 1.081)\n",
    "solnN1p081 = solve_ivp(lambda t, phi: [phi[1], -2 * beta_val * phi[1] - omega0_val**2 * np.sin(phi[0]) + gamma_val * omega0_val**2 * np.cos(omega_val * t)], [0, tMax], [-np.pi/2, 0], method='RK45', t_eval=np.linspace(0, tMax, 500))\n",
    "\n",
    "eqnNonLin = eqn.subs(gamma, 1.0826)\n",
    "solnN1p0826 = solve_ivp(lambda t, phi: [phi[1], -2 * beta_val * phi[1] - omega0_val**2 * np.sin(phi[0]) + gamma_val * omega0_val**2 * np.cos(omega_val * t)], [0, tMax], [-np.pi/2, 0], method='RK45', t_eval=np.linspace(0, tMax, 500))\n",
    "\n",
    "\n",
    "# Now, plot these (on the next slide). You see a progression from period=1, to period=2, to period=4, to period=8. This phenomenon is called a Period-Doubling Cascade."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5e4b5d39",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "fig, axs = plt.subplots(4, 2, figsize=(10, 15))\n",
    "\n",
    "axs[0, 0].plot(solnN1p06.t, solnN1p06.y[0], linewidth=2)\n",
    "axs[0, 0].set_xlim(0, 10)\n",
    "axs[0, 0].set_ylim(-2.5, -1.7)\n",
    "\n",
    "axs[0, 1].plot(solnN1p06.t, solnN1p06.y[0], linewidth=2)\n",
    "axs[0, 1].set_xlim(28, 40)\n",
    "axs[0, 1].set_ylim(-2.5, -1.7)\n",
    "\n",
    "axs[1, 0].plot(solnN1p078.t, solnN1p078.y[0], linewidth=2)\n",
    "axs[1, 0].set_xlim(0, 10)\n",
    "axs[1, 0].set_ylim(-2.5, -1.7)\n",
    "\n",
    "axs[1, 1].plot(solnN1p078.t, solnN1p078.y[0], linewidth=2)\n",
    "axs[1, 1].set_xlim(28, 40)\n",
    "axs[1, 1].set_ylim(-2.5, -1.7)\n",
    "\n",
    "axs[2, 0].plot(solnN1p081.t, solnN1p081.y[0], linewidth=2)\n",
    "axs[2, 0].set_xlim(0, 10)\n",
    "axs[2, 0].set_ylim(-2.5, -1.7)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "537ab363",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get values of the function at regular periods\n",
    "\n",
    "values_N1p06 = [solnN1p06.y[0][i] for i in [30, 31, 32, 33, 34, 35, 36, 37, 38]]\n",
    "values_N1p078 = [solnN1p078.y[0][i] for i in [30, 31, 32, 33, 34, 35, 36, 37, 38]]\n",
    "values_N1p081 = [solnN1p081.y[0][i] for i in [30, 31, 32, 33, 34, 35, 36, 37, 38]]\n",
    "values_N1p0826 = [solnN1p0826.y[0][i] for i in [30, 31, 32, 33, 34, 35, 36, 37, 38]]\n",
    "\n",
    "formatted_N1p06 = [f\"{value:.4f}\" for value in values_N1p06]\n",
    "formatted_N1p078 = [f\"{value:.4f}\" for value in values_N1p078]\n",
    "formatted_N1p081 = [f\"{value:.4f}\" for value in values_N1p081]\n",
    "formatted_N1p0826 = [f\"{value:.4f}\" for value in values_N1p0826]\n",
    "\n",
    "print(formatted_N1p06)\n",
    "print(formatted_N1p078)\n",
    "print(formatted_N1p081)\n",
    "print(formatted_N1p0826)\n",
    "\n",
    "\n",
    "# These will be useful later, when we draw a \"bifurcation diagram\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee420332",
   "metadata": {},
   "outputs": [],
   "source": [
    "# One more for now: \\[Gamma]=1.09\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "tMax = 40\n",
    "gamma = 1.09\n",
    "\n",
    "def equation(phi, t):\n",
    "    return phi''(t) + 2 * (omega0 / 4) * phi'(t) + omega0^2 * np.sin(phi(t))\n",
    "\n",
    "solnNonLin = odeint(equation, [-np.pi/2, 0], np.linspace(0, tMax, 100))\n",
    "\n",
    "t_values = np.linspace(0, tMax, 100)\n",
    "phi_values = solnNonLin[:, 0]\n",
    "\n",
    "plt.plot(t_values, phi_values, linewidth=2)\n",
    "plt.xlabel('t')\n",
    "plt.ylabel('φ(t)')\n",
    "plt.title('Nonlinear Pendulum Motion')\n",
    "plt.show()\n",
    "\n",
    "### The motion is no longer periodic!\n",
    "\n",
    "phi_values_selected = solnNonLin[[24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40], 0]\n",
    "formatted_phi_values = [f'{value:.6f}' for value in phi_values_selected]\n",
    "\n",
    "print(formatted_phi_values)\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
}