Examples_of_PDES.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Try this notebook on Binder:\n",
    "[![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/echoi/compgeodyn/master?labpath=Examples_of_PDES.ipynb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load_ext watermark\n",
    "# %watermark -v -m -p numpy,matplotlib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# CERI 7315/8315 Computational Methods for Geodynamics: Examples of PDES"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Potential equation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let $\\Omega$ be a domain in $\\mathbb{R}^{2}$. Find a function in $\\Omega$ with\n",
    "\\begin{equation}\n",
    "u_{xx} + u_{yy} = 0.\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### One solution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Identify $\\mathbb{R}^{2}$ with the complex plane if $w(z) = u(x,y) + i v(x,y)$ is a holomorphic function in $\\Omega$. Then, $u$ and $v$ satisfy the potential equation.\n",
    "\n",
    "Roughly speaking, a _holomorphic_ function is a differentiable function of complex variable.\n",
    "That means the real and complex part of $w(z)$ should satisfy the following relations:\n",
    "\n",
    "\\begin{equation}\n",
    "\\begin{split}\n",
    "\\frac{\\partial u}{\\partial x} &= \\frac{\\partial v}{\\partial y} \\\\\n",
    "\\frac{\\partial u}{\\partial y} &= -\\frac{\\partial v}{\\partial x}\n",
    "\\end{split}\n",
    "\\end{equation}\n",
    "\n",
    "One can easily verify that $u$ and $v$ with these properties are the solutions to the potential equation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When $\\Omega$ is a disk, i.e., $\\Omega = \\left\\{(x,y) \\in \\mathbb{R}^{2} | x^{2} + y^{2} < 1 \\right\\}$, a complex solution exists.\n",
    "\n",
    "Since $z^{k} = (r e^{i\\phi})^{k}$ is holomorphic, $r^{k}\\cos k\\phi$ and $r^{k}\\sin k\\phi$ satisfy the potential equation.\n",
    "\n",
    "When the boundary values are given as a Fourier expansion using these functions,\n",
    "\n",
    "\\begin{equation}\n",
    " u( \\cos \\phi, \\sin \\phi ) = a_{0} + \\sum_{k=1}^{\\infty} \\left( a_{k} \\cos k\\phi + b_{k} \\sin k\\phi \\right),\n",
    "\\end{equation}\n",
    "\n",
    "the solution in the interior is \n",
    "\n",
    "\\begin{equation}\n",
    " u( r, \\phi ) = a_{0} + \\sum_{k=1}^{\\infty} r^{k} \\left( a_{k} \\cos k\\phi + b_{k} \\sin k\\phi \\right),\n",
    "\\end{equation}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Poisson equation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above potential equation is the special case of the Poisson equation:\n",
    "\n",
    "Let $\\Omega$ be a domain in $\\mathbb{R}^{d}$ ($d$ = 2 or 3), $f: \\Omega \\rightarrow \\mathbb{R}$ be a prescribed charge density in $\\Omega$. Then $u$ satisfying\n",
    "\n",
    "\\begin{equation}\n",
    "  - \\Delta u = f \\text{ in } \\Omega\n",
    "\\end{equation}\n",
    "\n",
    "describe the potential in $\\Omega$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Wave equation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{equation}\n",
    "\\begin{split}\n",
    "  u_{tt} &= u_{xx} \\\\\n",
    "  u_{tt} &= c^{2} \\Delta u.\n",
    "\\end{split}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To get a solution for $c=1$, we transform the variables\n",
    "\\begin{align}\n",
    "  \\xi &= x + t \\\\\n",
    "  \\eta &= x - t.\n",
    "\\end{align}\n",
    "\n",
    "Then the partial derivaties are given as\n",
    "\\begin{align}\n",
    " u_{x} &= u_{\\xi} + u_{\\eta} \\\\\n",
    " u_{t} &= u_{\\xi} - u_{\\eta},\n",
    "\\end{align}\n",
    "\n",
    "\\begin{align}\n",
    " u_{xx} &= u_{\\xi\\xi} + 2 u_{\\xi \\eta} + u_{\\eta\\eta} \\\\\n",
    " u_{tt} &= u_{\\xi\\xi} - 2 u_{\\xi \\eta} + u_{\\eta\\eta}.\n",
    "\\end{align}\n",
    "\n",
    "This result implies that $4 u_{\\xi\\eta} = 0$ to satisfy the wave equation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The general solution is\n",
    "\n",
    "\\begin{equation}\n",
    "  u(x,t) = \\phi(\\xi) + \\psi(\\eta) = \\phi(x+t) + \\psi(x-t),\n",
    "\\end{equation}\n",
    "where $\\phi$ and $\\psi$ are functions to be determined from the initial conditions.\n",
    "\n",
    "When $u(x,0) = f(x) = \\phi(x) + \\psi(x)$ and $u_{t}(x,0) = g(x) = \\phi^{\\prime} (x) + \\psi^{\\prime} (x)$, we get\n",
    "\n",
    "\\begin{align}\n",
    "  \\phi^{\\prime} = \\frac{1}{2} ( f^{\\prime} + g ) \\\\\n",
    "  \\psi^{\\prime} = \\frac{1}{2} ( f^{\\prime} - g ).\n",
    "\\end{align}\n",
    "\n",
    "From this, we can express $\\phi$ and $\\psi$ as follows:\n",
    "\n",
    "\\begin{align}\n",
    "  \\phi(\\xi) = \\frac{1}{2} f(\\xi) + \\frac{1}{2} \\int_{x_{0}}^{\\xi} g(s)\\,ds \\\\\n",
    "  \\psi(\\eta) = \\frac{1}{2} f(\\eta) - \\frac{1}{2} \\int_{x_{0}}^{\\eta} g(s)\\,ds.\n",
    "\\end{align}\n",
    "\n",
    "Finally, the solution $u(x,t)$ is given as\n",
    "\\begin{equation}\n",
    "  u(x,t) = \\frac{1}{2} \\left[ f(x+t) + f(x-t) \\right] + \\frac{1}{2} \\int_{x-t}^{x+t} g(s)\\,ds.\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If $f(x)$ is a triangle function,\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\n",
    "  f(x) =\n",
    "  \\begin{cases}\n",
    "   x+1 & \\text{ for } -1 < x \\le 0 \\\\\n",
    "   1-x & \\text{ for } 0 < x \\le 1 \\\\\n",
    "   0   & \\text{otherwise}.\n",
    "  \\end{cases}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "num_points = 97\n",
    "x = np.linspace(-3,3,num_points)\n",
    "\n",
    "def f(x):\n",
    "    val = 0.0\n",
    "    if x > -1.0 and x <= 0.0:\n",
    "        val = x + 1.0\n",
    "    elif x > 0.0 and x <= 1.0:\n",
    "        val = 1.0 - x\n",
    "    return val\n",
    "\n",
    "def eval_u(u, x, t):\n",
    "    for j in range(len(t)):\n",
    "        for i in range(num_points):\n",
    "            u[j,i] = 0.5 * ( f( x[i] - t[j] ) + f( x[i] + t[j] ) )\n",
    "\n",
    "t=np.array([0.0, 0.3, 0.7, 1.5])\n",
    "u = np.zeros((len(t), num_points))\n",
    "eval_u(u, x, t)\n",
    "\n",
    "plt.plot( x, u[0,:], linewidth=2, label=\"t=0\")\n",
    "plt.plot( x, u[1,:], linewidth=2, label=\"t=0.3\")\n",
    "plt.plot( x, u[2,:], linewidth=2, label=\"t=0.7\")\n",
    "plt.plot( x, u[3,:], linewidth=2, label=\"t=1.5\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Heat equation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{align}\n",
    " u_{t} = \\sigma \\, u_{xx} \\\\\n",
    " u_{t} = \\sigma \\, \\Delta u \\\\\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We review the solution for a rod with $\\sigma=1$ and length of $\\pi$ as well as boundary conditions $u(0,t) = u(\\pi,t) = 0$.\n",
    "\n",
    "We assume that the initial condition is given as a Fourier expansion:\n",
    "\\begin{equation}\n",
    "  u(x,0) = \\sum_{k=1}^{\\infty} a_{k} \\sin kx, \\quad 0 < x < \\pi.\n",
    "\\end{equation}\n",
    "\n",
    "Noting that $e^{-k^{2}t} \\sin kx$ satisfy the heat equation, we get the solution\n",
    "\\begin{equation}\n",
    "  u(x,t) = \\sum_{k=1}^{\\infty} a_{k} e^{-k^{2}t} \\sin kx, \\quad t \\ge 0.\n",
    "\\end{equation}\n",
    "\n",
    "Note that high-frequency components decay fast and for that reason, the diffusion process cannot be reversed: If $t < 0$, those high frequency components grow fast, leading to the divergence of the solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "numx = 101\n",
    "x = np.linspace(0, np.pi, numx)\n",
    "numt = 21\n",
    "t = np.linspace(0, 2.0, numt)\n",
    "u = np.zeros((numt, numx))\n",
    "for i in range(numt):\n",
    "    for j in range(numx):\n",
    "        u[i,j] = np.exp(-t[i]) * np.sin(x[j]) + np.exp(-4.0*t[i])*np.sin(2.0*x[j])\n",
    "plt.plot(x,u[0,:])\n",
    "plt.plot(x,u[5,:])\n",
    "plt.plot(x,u[10,:])\n",
    "plt.plot(x,u[15,:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "nt = np.linspace(0, -2.0, numt)\n",
    "nu = np.zeros((numt, numx))\n",
    "for i in range(numt):\n",
    "    for j in range(numx):\n",
    "        nu[i,j] = np.exp(-nt[i]) * np.sin(x[j]) + np.exp(-4.0*nt[i])*np.sin(2.0*x[j])\n",
    "plt.plot(x,nu[0,:])\n",
    "plt.plot(x,nu[3,:])\n",
    "plt.plot(x,nu[6,:])\n",
    "plt.plot(x,nu[9,:])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.8.8 ('base')",
   "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.8.8"
  },
  "vscode": {
   "interpreter": {
    "hash": "ac42d88c93b784c771ddc161da243913e0e94454e081953e68028d558236f8c0"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}