{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "48974d97",
   "metadata": {},
   "source": [
    "# Eigenvalue problem: Coupled oscillators \n",
    "- by Börge Göbel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "ce43e36c",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import integrate\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "421f538a",
   "metadata": {},
   "source": [
    "![Coupled_oscillators](figure_08_coupled_oscillators.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef352f47",
   "metadata": {},
   "source": [
    "\\begin{align}\n",
    "m\\ddot{r}_1 &= -k(r_1-r_l-a) -k(r_1-r_2-a)\\\\\n",
    "m\\ddot{r}_2 &= -k(r_2-r_1-a) -k(r_2-r_3-a)\\\\\n",
    "m\\ddot{r}_3 &= -k(r_3-r_2-a) -k(r_3-r_r-a)\n",
    "\\end{align}\n",
    "\n",
    "The constants can be dropped, since we can transform \\\\( (r_1 - a) \\rightarrow r_1, (r_2 - 2a) \\rightarrow r_2, (r_3 - 3a) \\rightarrow r_3, \\\\) and use \\\\(r_l=0, r_r=4a\\\\)\n",
    "\n",
    "\\begin{align}\n",
    "m\\ddot{r}_1 &= -k(r_1) -k(r_1-r_2)\\\\\n",
    "m\\ddot{r}_2 &= -k(r_2-r_1) -k(r_2-r_3)\\\\\n",
    "m\\ddot{r}_3 &= -k(r_3-r_2) -k(r_3)\n",
    "\\end{align}\n",
    "\n",
    "This is identical to\n",
    "\n",
    "\\begin{align}\n",
    "\\begin{pmatrix}\\ddot{r}_1\\\\\\ddot{r}_2\\\\\\ddot{r}_3\\end{pmatrix}=-\\frac{k}{m}\n",
    "\\begin{pmatrix}2&-1&0\\\\-1&2&-1\\\\0&-1&2\\end{pmatrix}\n",
    "\\begin{pmatrix}r_1\\\\ r_2\\\\ r_3\\end{pmatrix}\n",
    "\\end{align}\n",
    "\n",
    "We can, of course, solve this system of equations numerically. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "86a7ad2c",
   "metadata": {},
   "source": [
    "## 1. Solving the coupled differential equations numerically "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b46f20b",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "6a70a97b",
   "metadata": {},
   "source": [
    "## 2. The eigenvalue problem"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "462014cf",
   "metadata": {},
   "source": [
    "### 2.1 Why is it an eigenvalue problem?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6392a25b",
   "metadata": {},
   "source": [
    "The equation for 3 oscillators was \n",
    "\n",
    "\\begin{align}\n",
    "\\begin{pmatrix}\\ddot{r}_1\\\\\\ddot{r}_2\\\\\\ddot{r}_3\\end{pmatrix}=-\\frac{k}{m}\n",
    "\\begin{pmatrix}2&-1&0\\\\-1&2&-1\\\\0&-1&2\\end{pmatrix}\n",
    "\\begin{pmatrix}r_1\\\\ r_2\\\\ r_3\\end{pmatrix}\n",
    "\\end{align}\n",
    "\n",
    "Solving the system numerically is nice. However, it would be much easier if the matrix would be just a diagonal matrix or even a scalar:\n",
    "\n",
    "\\begin{align}\n",
    "\\begin{pmatrix}\\ddot{q}_1\\\\\\ddot{q}_2\\\\\\ddot{q}_3\\end{pmatrix}=-\\frac{k}{m}\n",
    "\\lambda\n",
    "\\begin{pmatrix}q_1\\\\ q_2\\\\ q_3\\end{pmatrix}\n",
    "\\end{align}\n",
    "\n",
    "In this case the solution would be harmonic oscillators with a frequency \\\\( \\omega = \\sqrt{\\frac{k}{m}\\lambda}\\\\). \n",
    "\n",
    "We must find a unitary matrix \\\\(\\underline{u}\\\\) with \\\\(\\underline{u}^{-1}\\underline{u} = \\underline{u}\\,\\underline{u}^{-1} = 1\\\\) and \\\\( \\underline{u}^{-1}\\underline{A}\\,\\underline{u} = \\lambda\\\\) so that\n",
    "\n",
    "\\begin{align} \n",
    "\\ddot{\\vec{r}}&=-\\frac{k}{m}\\underline{A}\\vec{r}\\\\\n",
    "\\ddot{\\vec{r}}&=-\\frac{k}{m}\\left(\\underline{u}\\,\\underline{u}^{-1}\\right)\\,\\underline{A}\\,\\left(\\underline{u}\\,\\underline{u}^{-1}\\right)\\vec{r}\\\\\n",
    "\\left(\\underline{u}^{-1}\\ddot{\\vec{r}}\\right)&=-\\frac{k}{m}\\left(\\underline{u}^{-1}\\,\\underline{A}\\,\\underline{u}\\right)\\,\\left(\\underline{u}^{-1}\\vec{r}\\right)\\\\\n",
    "\\ddot{\\vec{q}}&=-\\frac{k}{m}\\lambda \\vec{q}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1a38f224",
   "metadata": {},
   "source": [
    "This means we must find:\n",
    "\n",
    "\\begin{align}\n",
    "\\underline{u}^{-1}\\,\\underline{A}\\,\\underline{u} &= \\lambda\\\\\n",
    "\\underline{A}\\,\\underline{u} &= \\lambda\\underline{u}\n",
    "\\end{align}\n",
    "\n",
    "or\n",
    "\n",
    "\\begin{align}\n",
    "\\underline{u}^{-1}\\,\\underline{A}\\,\\underline{u}\\vec{r} &= \\lambda\\vec{r}\\\\\n",
    "\\underline{A}\\,\\underline{u}\\vec{r} &= \\lambda\\underline{u}\\vec{r}\\\\\n",
    "\\end{align}\n",
    "\n",
    "In other words, we must find the eigenvalues and eigenvectors of \\\\( \\underline{A} \\\\)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3751672b",
   "metadata": {},
   "source": [
    "### a) Calculate the eigenvalues using a numpy routine"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "428ab94f",
   "metadata": {},
   "source": [
    "We want to determine the eigenvalues of the matrix\n",
    "\n",
    "\\begin{align}\n",
    "\\underline{A} = \\begin{pmatrix}2&-1&0\\\\-1&2&-1\\\\0&-1&2\\end{pmatrix}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "218c128a",
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.array([[2,-1,0],[-1,2,-1],[0,-1,2]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "408a3296",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "5959358c",
   "metadata": {},
   "source": [
    "### b) Exercise: Program a routine ourselves"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5e3f01ba",
   "metadata": {},
   "source": [
    "We have to solve the following equation\n",
    "\n",
    "\\begin{align}\n",
    "0 = \\det(\\underline{A}-\\lambda\\underline{1})\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5d4856c1",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "2d0916fd",
   "metadata": {},
   "source": [
    "\\begin{align}\n",
    "0 &= \\det\\begin{pmatrix}2-\\lambda&-1&0\\\\-1&2-\\lambda&-1\\\\0&-1&2-\\lambda\\end{pmatrix}\\\\\n",
    "0 &= (2-\\lambda)\\det\\begin{pmatrix}2-\\lambda&-1\\\\-1&2-\\lambda\\end{pmatrix}-(-1)\\det\\begin{pmatrix}-1&-1\\\\0&2-\\lambda\\end{pmatrix}\\\\\n",
    "0 &= (2-\\lambda)^3-(2-\\lambda)-(2-\\lambda)\\\\\n",
    "0 &= -\\lambda^3 +6\\lambda^2 - 12\\lambda + 8 -4 +2\\lambda\n",
    "\\end{align}\n",
    "\n",
    "We have to find the roots of a cubic function (characteristic polynomial)\n",
    "\n",
    "\\begin{align}\n",
    "0 &= -\\lambda^3 +6\\lambda^2 - 10\\lambda + 4\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d799f5d4",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "1bc60590",
   "metadata": {},
   "source": [
    "### c) Analyzing the eigensystem"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a59c622a",
   "metadata": {},
   "source": [
    "The analytical solutions are:\n",
    "\\\\( \\quad \\lambda_1 = 2+\\sqrt{2}, \\lambda_2 = 2, \\quad \\lambda_3 = 2-\\sqrt{2}\\\\)\n",
    "\n",
    "This leads us to the eigenfrequencies: \\\\( \\omega = \\sqrt{\\frac{k}{m}\\lambda}\\\\)\n",
    "\n",
    "\\begin{align} \n",
    "\\omega_1 = \\sqrt{\\frac{k}{m}}\\sqrt{2+\\sqrt{2}}, \\quad \\omega_2 = \\sqrt{\\frac{k}{m}}\\sqrt{2}, \\quad \\omega_3 = \\sqrt{\\frac{k}{m}}\\sqrt{2-\\sqrt{2}}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "78cc94f8",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "6fefcf01",
   "metadata": {},
   "source": [
    "The eigenvectors are \n",
    "\\begin{align}\n",
    "\\begin{pmatrix}-\\frac{1}{2}\\\\-\\frac{1}{\\sqrt{2}}\\\\\\frac{1}{2}\\end{pmatrix}\\quad\\quad\n",
    "\\begin{pmatrix}\\frac{1}{\\sqrt{2}}\\\\0\\\\\\frac{1}{\\sqrt{2}}\\end{pmatrix}\\quad\\quad\n",
    "\\begin{pmatrix}-\\frac{1}{2}\\\\\\frac{1}{\\sqrt{2}}\\\\\\frac{1}{2}\\end{pmatrix}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2971acc0",
   "metadata": {},
   "source": [
    "## 3. Fourier transform"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "734a9800",
   "metadata": {},
   "source": [
    "\\\\(\n",
    "\\tilde{y}(\\omega) = \\frac{1}{\\sqrt{2\\pi}} \\int_{-\\infty}^{\\infty} y(t) e^{i\\omega t}\\mathrm{d}t\n",
    "\\\\)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "3bd63c85",
   "metadata": {},
   "outputs": [],
   "source": [
    "def integralTrapezoidal(data):\n",
    "    a = 0\n",
    "    for i in range( len(data[0]) -1 ):\n",
    "        a = a + ( data[1,i+1] + data[1,i] ) / 2 * (data[0,i+1] - data[0,i])\n",
    "    return a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "81b082f5",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "a5999d83",
   "metadata": {},
   "source": [
    "Eigenvalues:\n",
    "\n",
    "\\begin{align} \n",
    "\\omega_1 = \\sqrt{\\frac{k}{m}}\\sqrt{2+\\sqrt{2}}, \\quad \\omega_2 = \\sqrt{\\frac{k}{m}}\\sqrt{2}, \\quad \\omega_3 = \\sqrt{\\frac{k}{m}}\\sqrt{2-\\sqrt{2}}\n",
    "\\end{align}\n",
    "\n",
    "Eigenvectors:\n",
    "\n",
    "\\begin{align}\n",
    "\\begin{pmatrix}-\\frac{1}{2}\\\\-\\frac{\\sqrt{3}}{2}\\\\\\frac{1}{2}\\end{pmatrix},\\quad\\quad\n",
    "\\begin{pmatrix}\\frac{\\sqrt{3}}{2}\\\\0\\\\\\frac{\\sqrt{3}}{2}\\end{pmatrix},\\quad\\quad\n",
    "\\begin{pmatrix}-\\frac{1}{2}\\\\\\frac{\\sqrt{3}}{2}\\\\\\frac{1}{2}\\end{pmatrix}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "75377957",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "d4b025b0",
   "metadata": {},
   "source": [
    "## 4. Fit harmonic functions with the eigenfrequency to the numerical solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8924520",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "17ef80b2",
   "metadata": {},
   "source": [
    "We consider the model function \n",
    "\n",
    "\\\\( f(t) = A_1\\cos(\\omega_1t+\\phi_1) + A_2\\cos(\\omega_2t+\\phi_2) + A_3\\cos(\\omega_3t+\\phi_3)\\\\)\n",
    "\n",
    "The fitting parameters are\n",
    "\n",
    "\\\\( \\vec{a} = (A_1, \\phi_1, A_2, \\phi_2, A_3, \\phi_3)\\\\)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "a65fab92",
   "metadata": {},
   "outputs": [],
   "source": [
    "def functionModel(t, a):\n",
    "    return a[0]*np.cos(omega1*t+a[1]) + a[2]*np.cos(omega2*t+a[3]) + a[4]*np.cos(omega3*t+a[5])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee572618",
   "metadata": {},
   "source": [
    "There are many reasonable definitions of an error function but a very common choice is: \\\\( \\Delta = \\sum_{i=1}^n \\left(y_i-f(x_i)\\right)^2\\\\)\n",
    "\n",
    "\\\\( f \\\\) is the fit function that is determined by the coefficients \\\\( a_i \\\\) in our case.\n",
    "\n",
    "\\\\( (x_i, y_i) \\\\) are the data points that we try to fit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "562af681",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "a74ea56d",
   "metadata": {},
   "source": [
    "We can use several different methods to minimize the error, e. g. a Monte-Carlo algorithm. Here, we will use the gradient descent method. The coefficients \\\\( a_i \\\\) will be updated along the gradient direction of the error function \\\\( \\nabla_{\\vec{a}} \\Delta\\\\). The gradient consists of elements \\\\( \\frac{\\partial}{\\partial a_k} \\Delta = -2 \\sum_{i=1}^n \\left(y_i-f(x_i)\\right) \\frac{\\partial}{\\partial a_k}f(x_i) \\\\)\n",
    "\n",
    "\\\\( \\nabla_{\\vec{a}} \\Delta = \\begin{pmatrix}\n",
    "\\cos(\\omega_1t+\\phi_1)\\\\\n",
    "-A_1\\sin(\\omega_1t+\\phi_1)\\\\\n",
    "\\cos(\\omega_2t+\\phi_2)\\\\\n",
    "-A_2\\sin(\\omega_2t+\\phi_2)\\\\\n",
    "\\cos(\\omega_3t+\\phi_3)\\\\\n",
    "-A_3\\sin(\\omega_3t+\\phi_3)\n",
    "\\end{pmatrix}\\\\)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "be90b0f2",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "c622375b",
   "metadata": {},
   "source": [
    "## 5. Generalization to n oscillators"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "424c4358",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "6b0e121f",
   "metadata": {},
   "source": [
    "### 5.1 Numerical solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e1d635a9",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "7ec3537d",
   "metadata": {},
   "source": [
    "### 5.2 Eigenfrequencies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1993074f",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "e7fb7035",
   "metadata": {},
   "source": [
    "### 5.3 Fourier transform"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c7a7bf9d",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "e33a2996",
   "metadata": {},
   "source": [
    "## 6. Alternative geometry: Periodic boundary conditions"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "3d3109d6",
   "metadata": {},
   "source": [
    "![Coupled_oscillators_circle](figure_08_coupled_oscillators_circle.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38f94aae",
   "metadata": {},
   "source": [
    "\\\\( \n",
    "\\begin{pmatrix}\\ddot{r}_1\\\\\\ddot{r}_2\\\\\\ddot{r}_3\\\\\\ddot{r}_4\\\\\\ddot{r}_5\\\\\\ddot{r}_6\\end{pmatrix}=-\\frac{k}{m}\n",
    "\\begin{pmatrix}2&-1&0&0&0&-1\\\\-1&2&-1&0&0&0\\\\0&-1&2&-1&0&0\\\\0&0&-1&2&-1&0\\\\0&0&0&-1&2&-1\\\\-1&0&0&0&-1&2&\\end{pmatrix}\n",
    "\\begin{pmatrix}r_1\\\\ r_2\\\\ r_3\\\\r_4\\\\ r_5\\\\ r_6\\end{pmatrix}\n",
    "\\\\)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3e42f950",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "9d476ab4",
   "metadata": {},
   "source": [
    "### 6.1 Numerical solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b9e965f",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "24971161",
   "metadata": {},
   "source": [
    "### 6.2 Eigenfrequencies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "feac9741",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "9a48308a",
   "metadata": {},
   "source": [
    "The first eigenvalue is zero. The corresponding mode is the translation that is allowed due to the periodic boundary conditions. It was forbidden before, since the system had fixed edges.\n",
    "\n",
    "Numerically, this value will be a very small number. Here it is negative which gives an error when calculating the square root. The correct value for the square root (and the frequency) is zero: The eigenmode is not an oscillation but a translation."
   ]
  }
 ],
 "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}