{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "plt.style.use('default')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Coupled Oscillators\n",
    "## Lecture 10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mass-spring system\n",
    "\n",
    "Equation of motion:\n",
    "\n",
    "$$ x'' + \\omega_0^2 x = 0 $$\n",
    "\n",
    "where $\\omega_0 = \\sqrt{k/m}$ is the natural angular frequency. \n",
    "\n",
    "The solution is\n",
    "\n",
    "$$ x(t) = A \\cos(\\omega_0 t + \\phi) $$\n",
    "\n",
    "where the coefficients $A$ and $\\phi$ are determined from the \n",
    " initial conditions for $x(0)$ and $x'(0)$.\n",
    " \n",
    "*See handwritten notes for lecture.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Coupled mass-spring systems\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For a system of coupled harmonic oscillators, the motion can be understood as a superposition of independent normal modes.\n",
    "\n",
    "The trick is finding the normal modes coordinates\n",
    "\n",
    "$$ q_1, q_2, \\ldots $$\n",
    "\n",
    "that decouple the equations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Matrix Approach"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By rewriting the equations of motion as matrix equation\n",
    "\n",
    "$$ x'' = A x $$\n",
    "\n",
    "the normal modes are the *eigenvectors* of $A$ and the *eigenvalues* are the frequencies of the normal modes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Finding Eigenvectors and Eigenvalues with numpy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can create array using the np.array( ) function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.array([[-2, 1],[1,-2]])\n",
    "print(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "NumPy comes with many mathematical tools that are useful, including a sub-package for methods from linear algebra."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy.linalg as LA\n",
    "help(LA)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to calculate the eigenvectors and eigenvalues of the matrix $A$.  This can be done with the `LA.eig()` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "help(LA.eig)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "LA.eig(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w, v = LA.eig(A)\n",
    "print(w[0], v[:,0])\n",
    "print(w[1], v[:,1]) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can recognize these at the eigenvectors we found earlier\n",
    "\n",
    "$$ q_1 = \\frac{1}{\\sqrt{2}} \\Big( \\begin{matrix}\n",
    "      1\\\\\n",
    "      1\\\\\n",
    "   \\end{matrix}\\Big) $$\n",
    "\n",
    "$$ q_2 = \\frac{1}{\\sqrt{2}} \\Big( \\begin{matrix}\n",
    "      -1\\\\\n",
    "      1\\\\\n",
    "   \\end{matrix}\\Big) $$\n",
    "      \n",
    "with the corresponding eigenvalues -1 and -3."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Any general solution for the coupled oscillator is a linear combination of these two normal modes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = np.arange(0, 30, 0.01)\n",
    "q1 = v[:,0].reshape(2,1) * np.cos(np.sqrt(-w[0])* t)\n",
    "q2 = v[:,1].reshape(2,1) * np.cos(np.sqrt(-w[1])* t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note: the `.reshape(2,1)` is needed because we need to the eigenvector to be written as a column vector and not row vector."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With these two normal modes, we can generate arbitrary solutions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(2, 1)\n",
    "\n",
    "C1 = 1\n",
    "C2 = 1\n",
    "\n",
    "x = C1*q1 + C2*q2\n",
    "\n",
    "plt.sca(axes[0])\n",
    "plt.plot(t, x[0])\n",
    "plt.xlabel('$t$')\n",
    "plt.ylabel('$x_1$')\n",
    "\n",
    "plt.sca(axes[1])\n",
    "plt.plot(t, x[1])\n",
    "plt.xlabel('$t$')\n",
    "plt.ylabel('$x_2$')\n",
    "\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To visualize the different modes of the coupled mass-spring system, try this:\n",
    "\n",
    "https://phet.colorado.edu/sims/normal-modes/normal-modes_en.html\n",
    "\n",
    "(two masses, horizontal polarization)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Four mass, five spring system"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can extend this analysis to many more masses and many more springs:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.array([[-2, 1, 0, 0],\n",
    "              [ 1,-2, 1, 0],\n",
    "              [ 0, 1,-2, 1],\n",
    "              [ 0, 0, 1,-2]])\n",
    "print(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w, v = LA.eig(A)\n",
    "print(w[0], v[:,0])\n",
    "print(w[1], v[:,1]) \n",
    "print(w[2], v[:,2])\n",
    "print(w[3], v[:,3])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, the eigenvalues are not in order.  Let's sort them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "idx = w.argsort()[::-1]   \n",
    "w = w[idx]\n",
    "v = v[:,idx]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(4):\n",
    "    print(w[i], v[:,i])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can define our eigenmodes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = np.arange(0, 30, 0.01)\n",
    "q1 = v[:,0].reshape(4,1) * np.cos(np.sqrt(-w[0])* t)\n",
    "q2 = v[:,1].reshape(4,1) * np.cos(np.sqrt(-w[1])* t)\n",
    "q3 = v[:,2].reshape(4,1) * np.cos(np.sqrt(-w[2])* t)\n",
    "q4 = v[:,3].reshape(4,1) * np.cos(np.sqrt(-w[3])* t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And create arbitrary solutions!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(4,1, figsize=(8,8))\n",
    "\n",
    "C1 = 1\n",
    "C2 = 1\n",
    "C3 = 1\n",
    "C4 = 1\n",
    "\n",
    "x = C1*q1 + C2*q2 + C3*q3 + C4*q4\n",
    "\n",
    "for i in range(4):\n",
    "    plt.sca(axes[i])\n",
    "    plt.plot(t, x[i])  \n",
    "    plt.ylabel('$x_{}$'.format(i))\n",
    "plt.xlabel('$t$')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, compare with the animation\n",
    "\n",
    "https://phet.colorado.edu/sims/normal-modes/normal-modes_en.html"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}