notebooks/Vibration_1d.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Vibrations of an elastic rod"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We consider the eigenvalue problem associated to the longitudinal vibrations of an elastic rod\n",
    "\n",
    "\\begin{align}\n",
    "  -\\partial_{xx} u = \\omega^2 u, \\quad \\Omega, \\quad u_{|\\partial \\Omega} = 0\n",
    "\\end{align}\n",
    "\n",
    "where $\\Omega = (0, 1)$.\n",
    "\n",
    "The weak formulation, for $u \\in \\mathcal{V}_h$, is given by \n",
    "\n",
    "\\begin{align}\n",
    "  a(u,v) = \\omega^2 b(u,v), \\quad \\forall v \\in \\mathcal{V}_h\n",
    "  \\label{eq:vibration-wf}\n",
    "\\end{align}\n",
    "\n",
    "where $a(u,v) = \\int_{\\Omega} \\nabla v \\cdot \\nabla u~d\\Omega$ and $b(u,v) = \\int_{\\Omega} v u~d\\Omega$.\n",
    "\n",
    "The bilinear forms $a$ and $b$ give the Stiffness and Mass matrices after discretization. Their associated GLT symbols are respectivaly $\\mathfrak{s}_p$ and $\\mathfrak{m}_p$. Therefor, the eigenvalues will be approximated by a uniform sampling of $\\frac{\\mathfrak{s}_p}{\\mathfrak{m}_p}$. \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Formal model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympde.calculus import grad, dot\n",
    "from sympde.topology import ScalarFunctionSpace\n",
    "from sympde.topology import Domain\n",
    "from sympde.topology import elements_of\n",
    "from sympde.expr     import BilinearForm\n",
    "from sympde.expr     import integral\n",
    "\n",
    "DIM = 1\n",
    "domain = Domain('\\Omega', dim=DIM)\n",
    "\n",
    "V = ScalarFunctionSpace('V', domain)\n",
    "\n",
    "u,v = elements_of(V, names='u,v')\n",
    "\n",
    "laplace = BilinearForm((u,v), \n",
    "                       integral(domain, dot(grad(v), grad(u))))\n",
    "mass    = BilinearForm((u,v), \n",
    "                       integral(domain, u*v))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Symbolic expressions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# imports from gelato\n",
    "from gelato.expr     import gelatize\n",
    "from gelato.printing import latex"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We shall need to defined the following symbols"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# imports from sympy\n",
    "from sympy import Symbol\n",
    "from sympy import ratsimp\n",
    "from sympy import expand\n",
    "\n",
    "tx = Symbol(\"tx\")\n",
    "nx = Symbol(\"nx\", integer=True)\n",
    "N  = Symbol(\"N\", integer=True)\n",
    "wl = Symbol('\\omega_l')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We compute the GLT expression by calling the function **gelatize** for a given degree $p$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = 3\n",
    "\n",
    "sl = gelatize(laplace, p)\n",
    "sm = gelatize(mass, p)\n",
    "\n",
    "expr = sl/sm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Print the latex expression using the **latex** function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import Math"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - \\frac{- 46872 \\cos{\\left (\\theta_x \\right )} + 168 \\cos{\\left (3 \\theta_x \\right )} - 59136}{5955 \\cos{\\left (\\theta_x \\right )} + 600 \\cos{\\left (2 \\theta_x \\right )} + 5 \\cos{\\left (3 \\theta_x \\right )} + 6040} - \\frac{42}{5}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Math(latex(ratsimp(expr/nx**2)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "expr = expr.subs({tx: wl / N})\n",
    "expr = expr.subs({nx: N-p}) #N = nx + p\n",
    "expr = expand(expr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us print again our expression in Latex,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - \\frac{N^{2} \\cos{\\left (\\frac{\\omega_l}{N} \\right )}}{4 \\left(\\frac{397 \\cos{\\left (\\frac{\\omega_l}{N} \\right )}}{840} + \\frac{\\cos{\\left (\\frac{2 \\omega_l}{N} \\right )}}{21} + \\frac{\\cos{\\left (\\frac{3 \\omega_l}{N} \\right )}}{2520} + \\frac{151}{315}\\right)} - \\frac{2 N^{2} \\cos{\\left (\\frac{2 \\omega_l}{N} \\right )}}{5 \\left(\\frac{397 \\cos{\\left (\\frac{\\omega_l}{N} \\right )}}{840} + \\frac{\\cos{\\left (\\frac{2 \\omega_l}{N} \\right )}}{21} + \\frac{\\cos{\\left (\\frac{3 \\omega_l}{N} \\right )}}{2520} + \\frac{151}{315}\\right)} - \\frac{N^{2} \\cos{\\left (\\frac{3 \\omega_l}{N} \\right )}}{60 \\left(\\frac{397 \\cos{\\left (\\frac{\\omega_l}{N} \\right )}}{840} + \\frac{\\cos{\\left (\\frac{2 \\omega_l}{N} \\right )}}{21} + \\frac{\\cos{\\left (\\frac{3 \\omega_l}{N} \\right )}}{2520} + \\frac{151}{315}\\right)} + \\frac{2 N^{2}}{3 \\left(\\frac{397 \\cos{\\left (\\frac{\\omega_l}{N} \\right )}}{840} + \\frac{\\cos{\\left (\\frac{2 \\omega_l}{N} \\right )}}{21} + \\frac{\\cos{\\left (\\frac{3 \\omega_l}{N} \\right )}}{2520} + \\frac{151}{315}\\right)} + \\frac{3 N \\cos{\\left (\\frac{\\omega_l}{N} \\right )}}{2 \\left(\\frac{397 \\cos{\\left (\\frac{\\omega_l}{N} \\right )}}{840} + \\frac{\\cos{\\left (\\frac{2 \\omega_l}{N} \\right )}}{21} + \\frac{\\cos{\\left (\\frac{3 \\omega_l}{N} \\right )}}{2520} + \\frac{151}{315}\\right)} + \\frac{12 N \\cos{\\left (\\frac{2 \\omega_l}{N} \\right )}}{5 \\left(\\frac{397 \\cos{\\left (\\frac{\\omega_l}{N} \\right )}}{840} + \\frac{\\cos{\\left (\\frac{2 \\omega_l}{N} \\right )}}{21} + \\frac{\\cos{\\left (\\frac{3 \\omega_l}{N} \\right )}}{2520} + \\frac{151}{315}\\right)} + \\frac{N \\cos{\\left (\\frac{3 \\omega_l}{N} \\right )}}{10 \\left(\\frac{397 \\cos{\\left (\\frac{\\omega_l}{N} \\right )}}{840} + \\frac{\\cos{\\left (\\frac{2 \\omega_l}{N} \\right )}}{21} + \\frac{\\cos{\\left (\\frac{3 \\omega_l}{N} \\right )}}{2520} + \\frac{151}{315}\\right)} - \\frac{4 N}{\\frac{397 \\cos{\\left (\\frac{\\omega_l}{N} \\right )}}{840} + \\frac{\\cos{\\left (\\frac{2 \\omega_l}{N} \\right )}}{21} + \\frac{\\cos{\\left (\\frac{3 \\omega_l}{N} \\right )}}{2520} + \\frac{151}{315}} - \\frac{9 \\cos{\\left (\\frac{\\omega_l}{N} \\right )}}{4 \\left(\\frac{397 \\cos{\\left (\\frac{\\omega_l}{N} \\right )}}{840} + \\frac{\\cos{\\left (\\frac{2 \\omega_l}{N} \\right )}}{21} + \\frac{\\cos{\\left (\\frac{3 \\omega_l}{N} \\right )}}{2520} + \\frac{151}{315}\\right)} - \\frac{18 \\cos{\\left (\\frac{2 \\omega_l}{N} \\right )}}{5 \\left(\\frac{397 \\cos{\\left (\\frac{\\omega_l}{N} \\right )}}{840} + \\frac{\\cos{\\left (\\frac{2 \\omega_l}{N} \\right )}}{21} + \\frac{\\cos{\\left (\\frac{3 \\omega_l}{N} \\right )}}{2520} + \\frac{151}{315}\\right)} - \\frac{3 \\cos{\\left (\\frac{3 \\omega_l}{N} \\right )}}{20 \\left(\\frac{397 \\cos{\\left (\\frac{\\omega_l}{N} \\right )}}{840} + \\frac{\\cos{\\left (\\frac{2 \\omega_l}{N} \\right )}}{21} + \\frac{\\cos{\\left (\\frac{3 \\omega_l}{N} \\right )}}{2520} + \\frac{151}{315}\\right)} + \\frac{6}{\\frac{397 \\cos{\\left (\\frac{\\omega_l}{N} \\right )}}{840} + \\frac{\\cos{\\left (\\frac{2 \\omega_l}{N} \\right )}}{21} + \\frac{\\cos{\\left (\\frac{3 \\omega_l}{N} \\right )}}{2520} + \\frac{151}{315}}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Math(latex(expr))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Well this can lead to a complicated expression. Let us take the expansion of the eigenvalues with respect to $\\frac{1}{N}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# here we consider an expansion up to the forth order\n",
    "order = 4\n",
    "expr = expr.series(1/N, 0, order)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let us check again the eigenvalue expression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{- 6 N \\omega_l^{2} + 9 \\omega_l^{2}}{N^{2}} + \\omega_l^{2} + O\\left(\\frac{1}{N^{4}}; N\\rightarrow \\infty\\right)$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Math(latex(expr))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<link href='http://fonts.googleapis.com/css?family=Fenix' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' type='text/css'>\n",
       "<style>\n",
       "    @font-face {\n",
       "        font-family: \"Computer Modern\";\n",
       "        src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
       "    }\n",
       "    div.cell{\n",
       "        width:600px;\n",
       "        margin-left:16% !important;\n",
       "        margin-right:auto;\n",
       "    }\n",
       "    h1 {\n",
       "        font-family: 'Alegreya Sans', sans-serif;\n",
       "    }\n",
       "    h2 {\n",
       "        font-family: 'Fenix', serif;\n",
       "    }\n",
       "    h3{\n",
       "\t\tfont-family: 'Fenix', serif;\n",
       "        margin-top:12px;\n",
       "        margin-bottom: 3px;\n",
       "       }\n",
       "\th4{\n",
       "\t\tfont-family: 'Fenix', serif;\n",
       "       }\n",
       "    h5 {\n",
       "        font-family: 'Alegreya Sans', sans-serif;\n",
       "    }\t   \n",
       "    div.text_cell_render{\n",
       "        font-family: 'Alegreya Sans',Computer Modern, \"Helvetica Neue\", Arial, Helvetica, Geneva, sans-serif;\n",
       "        line-height: 135%;\n",
       "        font-size: 120%;\n",
       "        width:600px;\n",
       "        margin-left:auto;\n",
       "        margin-right:auto;\n",
       "    }\n",
       "    .CodeMirror{\n",
       "            font-family: \"Source Code Pro\";\n",
       "\t\t\tfont-size: 90%;\n",
       "    }\n",
       "/*    .prompt{\n",
       "        display: None;\n",
       "    }*/\n",
       "    .text_cell_render h1 {\n",
       "        font-weight: 200;\n",
       "        font-size: 50pt;\n",
       "\t\tline-height: 100%;\n",
       "        color:#054BCD;\n",
       "        margin-bottom: 0.5em;\n",
       "        margin-top: 0.5em;\n",
       "        display: block;\n",
       "    }\t\n",
       "    .text_cell_render h5 {\n",
       "        font-weight: 300;\n",
       "        font-size: 16pt;\n",
       "        color: #054BCD;\n",
       "        font-style: italic;\n",
       "        margin-bottom: .5em;\n",
       "        margin-top: 0.5em;\n",
       "        display: block;\n",
       "    }\n",
       "    \n",
       "    .warning{\n",
       "        color: rgb( 240, 20, 20 )\n",
       "        }  \n",
       "</style>\n",
       "<script>\n",
       "    MathJax.Hub.Config({\n",
       "                        TeX: {\n",
       "                           extensions: [\"AMSmath.js\"]\n",
       "                           },\n",
       "                tex2jax: {\n",
       "                    inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
       "                    displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
       "                },\n",
       "                displayAlign: 'center', // Change this to 'center' to center equations.\n",
       "                \"HTML-CSS\": {\n",
       "                    styles: {'.MathJax_Display': {\"margin\": 4}}\n",
       "                }\n",
       "        });\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# css style\n",
    "from IPython.core.display import HTML\n",
    "def css_styling():\n",
    "    styles = open(\"../styles/custom.css\", \"r\").read()\n",
    "    return HTML(styles)\n",
    "css_styling()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}