02_PDEs.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "<table>\n",
    " <tr align=left><td><img align=left src=\"https://i.creativecommons.org/l/by/4.0/88x31.png\">\n",
    " <td>Text provided under a Creative Commons Attribution license, CC-BY. All code is made available under the FSF-approved MIT license. (c) Kyle T. Mandli</td>\n",
    "</table>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Review: Partial Differential Equations\n",
    "\n",
    "The purpose of this lecture is to review the bare-essential background for our discussion of numerical methods for  partial different equations (PDEs).  Most of what we will use is general knowledge of PDE classification, solutions techniques for each type of PDE, and some knowledge of the solution properties.  This also includes questions related to well-posedness, boundary conditions, and other issues."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "First some notation to get us started, we will use the notation\n",
    "$$\n",
    "    \\frac{\\partial u}{\\partial x} = u_x\n",
    "$$\n",
    "to denote a partial derivative of the function $u$ with respect to the variable $x$.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Higher order or mixed derivatives are expressed similarly as\n",
    "$$\n",
    "    \\frac{\\partial^2 u}{\\partial x^2} = u_{xx} \\quad \\quad \\frac{\\partial^2 u}{\\partial y \\partial x} = u_{xy}\n",
    "$$\n",
    "noting that in general that the order of the derivatives can matter.  \n",
    "\n",
    "When can we reorder the derivatives in these cases?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Here we have singled out variables usually associated with cartesian space but it will also be useful to consider polar-coordinate derivatives such as\n",
    "$$\n",
    "    \\frac{\\partial u}{\\partial \\theta} = u_\\theta\n",
    "$$\n",
    "or in time\n",
    "$$\n",
    "    \\frac{\\partial u}{\\partial t} = u_t.\n",
    "$$    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## 1. PDE Classification\n",
    "\n",
    "Classical classification of PDEs is based on a general form of a linear, second order PDE looking like\n",
    "$$\n",
    "    a u_{xx} + 2b u_{xy} + c u_{yy} + d u_x + e u_y + f u = g.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Note that in the general framework for classifications that $x$ or $y$ could also be $t$ or any other variable.  Formally taking a Fourier transform of \n",
    "$$\n",
    "    a u_{xx} + 2b u_{xy} + c u_{yy} + d u_x + e u_y + f u = g.\n",
    "$$\n",
    "allows us to transform this general equation into\n",
    "$$\n",
    "    P(X,Y) = a X^2 + 2b X \\cdot Y + c Y^2 + \\cdots\n",
    "$$\n",
    "which looks like a general, second order polynomial.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The terminology for PDE classification comes from the fact that the equation above describes general [conic sections](https://en.wikipedia.org/wiki/Conic_section) which are often classified as elliptic, parabolic, or hyperbolic based on their discriminant\n",
    "$$\n",
    "    b^2 - a c.\n",
    "$$\n",
    "We then classify the PDE according to the sign of the discriminant so that\n",
    " - if $b^2 - a c < 0$ the PDE is elliptic,\n",
    " - if $b^2 - a c = 0$ the PDE is parabolic, and\n",
    " - if $b^2 - a c > 0$ the PDE is hyperbolic."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The reason why we care about classifications is that often there are a class of numerical methods that work well for one classification but not for others.  We will investigate more on why and how this works as we introduce each type of PDE but for now lets examine some of the basic examples of PDEs from each classification."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Examples\n",
    "\n",
    "**Laplace's Equation:** In the $n$ dimensional case Laplace's equation is\n",
    "$$\n",
    "    \\Delta u = \\nabla^2 u = \\sum^n_{i=1} u_{x_i x_i} = 0\n",
    "$$\n",
    "and when $n = 2$ using $x$ and $y$ as the dimensions we have\n",
    "$$\n",
    "    u_{xx} + u_{yy} = 0.\n",
    "$$\n",
    "What type of PDE is this?\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Identifying the coefficients from above as $a = 1$, $b = 0$, and $c = 1$ then the discriminant tells us that Laplace's equations is elliptic in character since\n",
    "$$\n",
    "    b^2 - a c = -1.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "**Heat Equation:**  The heat equation is defined in one spatial dimension as\n",
    "$$\n",
    "    u_t = u_{xx}.\n",
    "$$\n",
    "\n",
    "What type of PDE is the heat equation?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The heat equation is parabolic.  This can be seen by replacing the $y$-derivatives above with $t$ and therefore finding the coefficients $a=-1$, $b=0$, and $c=0$ meaning the discriminant is zero."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "What about the heat equation in multiple space dimensions?\n",
    "$$\n",
    "    u_t = u_{xx} + u_{yy}.\n",
    "$$\n",
    "Unfortunately the identification breaks apart here as we have three variables."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "**Wave Equation:** The wave equation in one spatial dimension is defined as\n",
    "$$\n",
    "    u_{tt} - u_{xx} = 0.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "By process of elimination you probably guessed that this equation is hyperbolic.  Again replacing the $y$-derivatives with $t$ we find $a = -1$, $b = 0$, $c = 1$ so that\n",
    "\n",
    "$$\n",
    "    b^2 - a c = -(-1)(1) = 1.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "**Navier-Stokes Equations**:  Models incompressible fluid flow\n",
    "$$\\begin{aligned}\n",
    "    u_t + (u^2)_x + (uv)_y = - \\frac{1}{\\rho} p_x + \\nu (u_{xx} + v_{yy})\\\\\n",
    "    v_t + (uv)_x + (v^2)_y = - \\frac{1}{\\rho} p_y + \\nu (u_{xx} + v_{yy})\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Even ignoring the fact that these equations are unclassifiable via our nomenclature they contain both something that looks like a hyperbolic piece (the left-hand side) and an elliptic piece (the right-hand side) or parabolic given the time derivative.  Because of this the Navier-Stokes equations are often considered parabolic however approaches to solving them often use techniques from each of the classifications."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "One of the more useful ways to think about PDE classification that is more general is to examine how information propagates through the domain.  \n",
    "\n",
    " - **Hyperbolic PDEs:** information travels at finite speeds (group velocity is real and finite).  \n",
    " - **Parabolic PDEs:** describe information movement that is infinite (think about what happens in the case of diffusion in a vacuum).  \n",
    " - **Elliptic PDEs:** are limiting case between hyperbolic and elliptic equations in terms of information propagation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## 2. Fourier Analysis\n",
    "\n",
    "Fourier analysis is one of the basic tools for analyzing and solving linear PDEs.  For finite domains we usually use Fourier series to represent the solution.  The use of the Fourier series is equivalent to considering the eigenfunctions of the differential operators as a basis of the solution.  This linear algebra interpretation of PDEs is immensely powerful and we will come back to this again and again."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Fourier Series\n",
    "\n",
    "The Fourier series of a function comes in a number of different flavors depending on the type of function we are considering.  Recall that $\\sin x$ is an odd function and $\\cos x$ is an even function.  On finite intervals then for odd functions we can use the Fourier $\\sin$ series and for even functions we can use the Fourier $\\cos$ series.  If the function is neither odd or even then we can use the full series, often written as\n",
    "$$\n",
    "    f(x) = a_0 + \\sum^\\infty_{n=1} a_n \\cos \\left( \\frac{n \\pi x}{L} \\right) + \\sum^\\infty_{n=1} b_n \\sin \\left( \\frac{n \\pi x}{L} \\right),\n",
    "$$\n",
    "where\n",
    "$$\\begin{aligned}\n",
    "    a_0 &= \\frac{1}{L} \\int^L_0 f(x) dx, \\\\\n",
    "    a_n &= \\frac{2}{L} \\int^L_0 f(x) \\cos \\frac{n \\pi x}{L} dx, \\\\\n",
    "    b_n &= \\frac{2}{L} \\int^L_0 f(x) \\sin \\frac{n \\pi x}{L} dx.\n",
    "\\end{aligned}$$\n",
    "These coefficients are projections of the function $f(x)$ onto the Fourier basis of $\\cos$ and $\\sin$ which form an orthogonal basis.\n",
    "\n",
    "How does this relate to linear algebra and finite dimesional vector spaces?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Fourier Transform\n",
    "\n",
    "The Fourier transform of a function $v(x)$ is defined by\n",
    "$$\n",
    "    \\widehat{v}(\\xi) = \\frac{1}{\\sqrt{2 \\pi}} \\int^\\infty_{-\\infty} v(x) e^{-i \\xi x} dx.\n",
    "$$\n",
    "In order for the transform to exist we require that the function $v \\in L^2$, that is\n",
    "$$\n",
    "    ||v||_2 = \\left ( \\int^\\infty_{-\\infty} \\left|v(x) \\right|^2 dx \\right )^{1/2} < \\infty.\n",
    "$$\n",
    "The function $\\widehat{v}(\\xi) \\in L^2$ since from Perseval's identity we have\n",
    "$$\n",
    "    ||v(x)||_2 = ||\\widehat{v}(\\xi)||_2.\n",
    "$$\n",
    "We can also use the inverse transform to get back the original function:\n",
    "$$\n",
    "    v(x) = \\frac{1}{\\sqrt{2 \\pi}} \\int^\\infty_{-\\infty} \\widehat{v}(\\xi) e^{i \\xi x} d\\xi.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "One of the more important properties of the Fourier Transform we have already used above is that the Fourier transform turns spatial derivatives into algebraic relations:\n",
    "$$\n",
    "    \\widehat{v_x} = i \\xi \\widehat{v}.\n",
    "$$\n",
    "This is a result of the functions \n",
    "$$\n",
    "    e^{i\\xi x}\n",
    "$$ \n",
    "being eigenfunctions of the operator\n",
    "$$\n",
    "    \\frac{\\partial}{\\partial x}\n",
    "$$\n",
    "This translates into the discrete case as well where we will find that the discrete analogs of the eigenfunctions will be eigenvectors of the equivalent finite difference operators."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "**Exercise:** Try taking the Fourier transform of the advection equation:\n",
    "$$\n",
    "    u_t + c u_x = 0.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Since the time derivative is not effected by the Fourier transform we end up with\n",
    "$$\n",
    "    \\widehat{u}_t = -i \\xi c \\widehat{u}\n",
    "$$\n",
    "which is an infinite system of ODEs.  Now try solving this system to find the solution to the advection equation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We know the solution to the ODEs are of the form\n",
    "$$\n",
    "    \\widehat{u}(\\xi, t) = \\widehat{u}(\\xi, 0) e^{-i \\xi c t}.\n",
    "$$\n",
    "Using the inverse transform we find\n",
    "$$\\begin{aligned}\n",
    "    u(x,t) &= \\frac{1}{\\sqrt{2 \\pi}} \\int^\\infty_{-\\infty} \\widehat{u}(\\xi, 0)  e^{-i \\xi c t} e^{i \\xi x} d\\xi, \\\\\n",
    "    &= \\frac{1}{\\sqrt{2 \\pi}} \\int^\\infty_{-\\infty} \\widehat{u}(\\xi, 0)  e^{i \\xi (x - c t)} d\\xi, \\\\\n",
    "    &= u(x - c t, 0).\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## 3. Elliptic Equations\n",
    "\n",
    "Elliptic equations are most often boundary value problems on some domain.  Above we mentioned Laplace's equation which is a special case of another classic elliptic equation, the Poisson problem\n",
    "$$\n",
    "    \\nabla^2 u = f.\n",
    "$$\n",
    "These types of problems are often found when we want to find the steady-state of some time dependent problem such as the heat equation.  These also can arise in time-dependent problems where there are some dynamics that is much faster than the time-scales of the rest of the problem, one example is the incompressibility condition used in the Navier-Stokes equation (the Poisson problem for pressure)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The essential behavior of elliptic boundary value problems is that the solution in the interior of the domain $\\Omega$ is determined completely by the boundary conditions and the function $f$.  This is important when considering numerical methods for these types of problems as we often must solve large, sparse system of linear equations which couple every point to every other point."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "In general a linear operator $L$ of the form\n",
    "$$\n",
    "    L = \\sum^n_{i,j=1} A_{ij} \\frac{\\partial^2}{\\partial x_i \\partial x_j} + \\sum^n_{i=1} B_i \\frac{\\partial}{\\partial x_i} + C\n",
    "$$\n",
    "can be termed elliptic if the matrix $A$ is positive or negative definite (i.e. $v^T A v$ has the same sign for $v \\neq 0$ and $v \\in \\mathbb{R}^n$)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Solving Elliptic Equations\n",
    "\n",
    "As an example of solving an elliptic PDE let's consider Laplace's equation on a rectangular domain\n",
    "$$\\begin{aligned}\n",
    "    &u_{xx} + u_{yy} = 0 \\quad \\Omega = (0, L) \\times (0, H) \\\\\n",
    "    &u(x,0) = f_1(x) \\\\\n",
    "    &u(x,H) = f_2(x) \\\\\n",
    "    &u(0,y) = g_1(y) \\\\\n",
    "    &u(L,y) = g_2(y).\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Since we have non-homogeneous boundary conditions on each side of the rectangle we cannot directly use separation of variables and instead break down the problem into 4 sub-parts so that\n",
    "$$\n",
    "    u(x,y) = u_1(x,y) + u_2(x,y) + u_3(x,y) + u_4(x,y)\n",
    "$$\n",
    "where each $u_i(x,y)$ has one of the non-homogeneous sides of the rectangle retained and the others are kept at 0.  For instance for $u_1(x,y)$ we have\n",
    "$$\\begin{aligned}\n",
    "    \\nabla^2 u_1 = 0 \\quad & \\quad \\Omega = (0, L) \\times (0, H) \\\\\n",
    "    u_1(x,0) &= f_1(x) \\\\\n",
    "    u_1(x,H) &= 0 \\\\\n",
    "    u_1(0,y) &= 0 \\\\\n",
    "    u_1(L,y) &= 0.\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Using separation of variables we assume a solution of the form $u_4(x,y) = h(x) \\phi(y)$ which substituted into the PDE leads to\n",
    "$$\\begin{aligned}\n",
    "    &\\phi \\frac{\\text{d}^2 h}{\\text{d} x^2} + h \\frac{\\text{d}^2 \\phi}{\\text{d} y^2} = 0 \\\\\n",
    "    &\\frac{1}{h} \\frac{\\text{d}^2 h}{\\text{d} x^2} = - \\frac{1}{\\phi} \\frac{\\text{d}^2 \\phi}{\\text{d} y^2} = -\\lambda \\Rightarrow \\\\\n",
    "    &h''(x) + \\lambda h(x) = 0 \\quad\\quad \\phi''(y) - \\lambda \\phi(y) = 0.\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The boundary conditions for the two functions are\n",
    "$$\n",
    "    h(0) = h(L) = 0 \\quad \\phi(0) = f_1(x) \\quad \\phi(H) = 0.\n",
    "$$\n",
    "First solving the easier homogeneous BVP for $h(x)$ we find\n",
    "$$\n",
    "    h(x) = A \\cos (\\sqrt{\\lambda} x) + B \\sin(\\sqrt{\\lambda} x)\n",
    "$$\n",
    "which using the BCs leads us to\n",
    "$$\\begin{aligned}\n",
    "    &h(0) = 0 = A \\\\\n",
    "    &h(L) = 0 = B \\sin(\\sqrt{\\lambda} L).\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The latter we have two choices to satisfy the BC, we could assume $B=0$ leading to the trivial solution or\n",
    "$$\n",
    "    \\sqrt{\\lambda} L = n \\pi.\n",
    "$$\n",
    "Using the latter condition leads to the eigenvalues \n",
    "$$\n",
    "    \\lambda_n = \\left(\\frac{n \\pi}{L}\\right)^2, \\quad n = 1,2,\\ldots\n",
    "$$\n",
    "so that we have the solutions\n",
    "$$\n",
    "    h_n(x) = B_n \\sin\\left(\\frac{n \\pi x}{L}\\right).\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "For the $\\phi$ equation we have solutions of the form\n",
    "$$\n",
    "    \\phi(y) = A e^{\\sqrt{\\lambda} y} + B e^{-\\sqrt{\\lambda} y}\n",
    "$$\n",
    "or\n",
    "$$\n",
    "    \\phi(y) = A \\cosh(\\sqrt{\\lambda} y) + B \\sinh(\\sqrt{\\lambda} y).\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "To satisfy the homogeneous boundary condition we need to modify these solutions slightly so that they are shifted\n",
    "$$\n",
    "    \\phi(y) = A \\cosh \\left(\\sqrt{\\lambda} (y - H)\\right) + B \\sinh \\left(\\sqrt{\\lambda} (y - H)\\right)\n",
    "$$\n",
    "which allows us to satisfy the homogeneous BC by\n",
    "$$\\begin{aligned}\n",
    "    \\phi(H) = 0 &= A \\cosh \\left(\\sqrt{\\lambda} (H - H)\\right) + B \\sinh \\left(\\sqrt{\\lambda} (H - H)\\right) \\\\\n",
    "    &=A \\cosh (0) + B \\sinh (0) \\\\\n",
    "    &=A\n",
    "\\end{aligned}$$\n",
    "and we are then left with\n",
    "$$\n",
    "    \\phi(y) = B \\sinh \\left(\\frac{n \\pi}{L} (y - H)\\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Combining the two solutions we then have\n",
    "$$\n",
    "    u_{1,n}(x,y) = B_n \\sin\\left(\\frac{n \\pi x}{L}\\right) \\sinh \\left(\\frac{n \\pi}{L} (y - H)\\right)\n",
    "$$\n",
    "and by superposition we have\n",
    "$$\n",
    "    u_1(x,y) = \\sum^\\infty_{i=1} B_n \\sin\\left(\\frac{n \\pi x}{L}\\right) \\sinh \\left(\\frac{n \\pi}{L} (y - H)\\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Now, to satisfy the last boundary condition and find the coefficients $B_n$ we have\n",
    "$$\n",
    "    u_1(x,0) = f_1(x) = \\sum^\\infty_{i=1} -B_n \\sin\\left(\\frac{n \\pi x}{L}\\right) \\sinh \\left(\\frac{n \\pi H} {L}\\right)\n",
    "$$\n",
    "and using the orthogonality of the $\\sin$ functions we have\n",
    "$$\\begin{aligned}\n",
    "    \\int^L_0 f_1(x) \\sin\\left(\\frac{m \\pi x}{L}\\right) dx &= \\sum^\\infty_{i=1} -B_n \\sinh \\left(\\frac{n \\pi H}{L} \\right ) \\int^L_0 \\sin\\left(\\frac{n \\pi x}{L}\\right) \\sin\\left(\\frac{m \\pi x}{L}\\right) dx \\\\\n",
    "    \\int^L_0 f_1(x) \\sin\\left(\\frac{m \\pi x}{L}\\right) dx &= -B_m \\frac{L}{2} \\sinh \\left(\\frac{n \\pi H}{L} \\right ) \\\\\n",
    "    B_m &= \\frac{-2}{L \\sinh \\left(\\frac{m \\pi H}{L} \\right )} \\int^L_0 f_1(x) \\sin\\left(\\frac{m \\pi x}{L}\\right) dx.\n",
    "\\end{aligned}$$\n",
    "We then have found the solution to the $u_1$ component of the general problem above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## 4. Parabolic Equations\n",
    "\n",
    "Parabolic equations are similar to elliptic equations however often also include time-dependence (this is not always the case but will be for our considerations).  We can for instance take any elliptic linear operator as defined above and create a parabolic equation by writing\n",
    "$$\n",
    "    u_t - Lu = f.\n",
    "$$\n",
    "\n",
    "The classic example of a parabolic equation is the heat (or diffusion) equation defined as\n",
    "$$\n",
    "    u_t = \\nabla^2 u.\n",
    "$$\n",
    "These PDEs require both an initial condition $u(x, 0)$ as well as boundary conditions on $\\Omega$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Solving Parabolic Equations\n",
    "\n",
    "Similar to elliptic problems we again will use Fourier series to solve parabolic equations.  Let us again consider our model problem, the heat equation, in one spatial dimension:\n",
    "$$\\begin{aligned}\n",
    "    &u_t = \\kappa u_{xx}, \\quad \\Omega = (0, T) \\times (0, L), \\\\\n",
    "    &u(x, 0) = f(x), \\\\\n",
    "    &u(0, t) = 0, \\quad \\text{and} \\\\\n",
    "    &u(L, t) = 0.\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Using separation of variables we assume $u(x,t) = G(t) \\phi(x)$ and can then write the two ODEs\n",
    "$$\\begin{aligned}\n",
    "    &\\frac{\\text{d} G}{\\text{d} t} = -\\kappa \\lambda G, \\\\\n",
    "    &\\phi''(x) + \\lambda \\phi(x) = 0.\n",
    "\\end{aligned}$$\n",
    "The time dependent ODE has solutions of the form\n",
    "$$\n",
    "    G(t) = C e^{-\\lambda \\kappa t}.\n",
    "$$\n",
    "The BVP for $\\phi$ has solutions of the form\n",
    "$$\n",
    "    \\phi(x) = A \\cos \\left( \\sqrt{\\lambda} x \\right) + B \\sin \\left( \\sqrt{\\lambda} x \\right).\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Using the boundary conditions we find that $A = 0$ and for non-trivial solutions that\n",
    "$$\n",
    "    \\lambda = \\left( \\frac{n \\pi}{L} \\right )^2,\n",
    "$$\n",
    "leading to solutions of the form\n",
    "$$\n",
    "    u_n(x,t) = B_n e^{-\\left(\\frac{n \\pi}{L}\\right)^2 \\kappa t} \\sin \\left( \\frac{n \\pi x}{L} \\right),\n",
    "$$\n",
    "and by superposition we have\n",
    "$$\n",
    "    u(x,t) = \\sum^\\infty_{i=1} B_n e^{-\\left(\\frac{n \\pi}{L}\\right)^2 \\kappa t} \\sin \\left( \\frac{n \\pi x}{L} \\right).\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Evaluating the initial condition we have\n",
    "$$\n",
    "    u(x,0) = f(x) = \\sum^\\infty_{i=1} B_n \\sin \\left( \\frac{n \\pi x}{L} \\right),\n",
    "$$\n",
    "which gives us the $B_n$ by again using the orthogonality of the $\\sin$ functions\n",
    "$$\\begin{aligned}\n",
    "    \\int^L_0 f(x) \\sin \\left( \\frac{m \\pi x}{L} \\right) dx &= \\sum^\\infty_{i=1} B_n \\int^L_0\\sin \\left( \\frac{n \\pi x}{L} \\right) \\sin \\left( \\frac{m \\pi x}{L} \\right) dx, \\\\\n",
    "    \\int^L_0 f(x) \\sin \\left( \\frac{m \\pi x}{L} \\right) dx &= B_m \\frac{L}{2}, \\\\\n",
    "    B_m &= \\frac{2}{L} \\int^L_0 f(x) \\sin \\left( \\frac{m \\pi x}{L} \\right) dx,\n",
    "\\end{aligned}$$\n",
    "therefore solving the PDE."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## 5. Hyperbolic Equations\n",
    "\n",
    "Hyperbolic equations often arise as conservation or balance laws.  For instance the conservation of momentum (Newton's second law) is an example of a balance law.  Similarly energy conservation can be written as a hyperbolic equation.  The classic hyperbolic equation is of course the wave equation\n",
    "$$\n",
    "    u_{tt} - c^2 u_{xx} = 0.\n",
    "$$\n",
    "However often when considering numerical methods for hyperbolic equations we decompose second order PDEs into a system of first order PDEs such as\n",
    "$$\n",
    "    u_t + A u_x = 0,\n",
    "$$\n",
    "where now $A$ is a matrix.  In this case the PDE is hyperbolic if the matrix $A$'s eigenvalues $\\lambda$ are real and the matrix is diagonalizable. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Solving Hyperbolic Equations\n",
    "\n",
    "The most common approach to solving hyperbolic PDEs is to find the characteristics of the given problem and drawing on the ($x$,$t$) plane places where something interesting might be happening.  If the equation is non-linear this could lead to shocks (discontinuities) even though the initial condition was smooth.  This will be problematic for us numerically as we shall see later in the course."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The basic approach to the method of characteristics employs the assumption that our solution is of the form\n",
    "$$\n",
    "    u(x,t) = u(x(t), t),\n",
    "$$\n",
    "in other words our solution can be described by a single characteristic variable, in this case $t$.  Let us suppose that we are trying to find the solution to\n",
    "$$\n",
    "    u_t + c(x,t,u) u_x = 0.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "If we take the full (or total) derivative of the new solution we find\n",
    "$$\n",
    "    \\frac{\\text{d} u}{\\text{d} t} = \\frac{\\partial u}{\\partial t} + \\frac{\\text{d} x}{\\text{d} t} \\frac{\\partial u}{\\partial x},\n",
    "$$\n",
    "that, after transforming the original equation to use $\\xi$, leads us to the conclusion that our original PDE can be rewritten to be two ODEs\n",
    "$$\\begin{aligned}\n",
    "    &\\frac{\\text{d} u}{\\text{d} t} = 0, \\quad \\text{and}\\\\\n",
    "    &\\frac{\\text{d} x}{\\text{d} t} = c(x,t,u).\n",
    "\\end{aligned}$$\n",
    "Now we know that $u(x(t),t)$ is constant along the curves $x(t)$ and in fact\n",
    "$$\n",
    "    u(x,t) = u(x(t), 0).\n",
    "$$\n",
    "\n",
    "The solution to the PDE then turns to solving for the characteristic curves given by $x(t)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "As a simple example lets solve the linear advection equation on an infinite domain\n",
    "$$\\begin{aligned}\n",
    "    &u_t + c u_x = 0, \\quad \\Omega = (-\\infty, \\infty) \\times (0, T) \\\\\n",
    "    &u(x,0) = f(x), \\quad \\text{and} \\\\\n",
    "    &u_x(\\pm \\infty, t) \\rightarrow 0.\n",
    "\\end{aligned}$$\n",
    "The characteristic equation is\n",
    "$$\n",
    "    \\frac{\\text{d} x}{\\text{d} t} = c,\n",
    "$$\n",
    "whose solution is $x(t) = c t + x_0$.  We can draw the characteristics then as\n",
    "\n",
    "![Constant Wave Speed](./images/characteristics.png)\n",
    "\n",
    "where the equations for the lines are $x(t) = c t + x_0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "For nonlinear equations, i.e. ones where the speed may be dependent on the solution $c(u)$, the initial condition will determine the characteristic slopes leading to two cases.  Let's consider Burger's equation\n",
    "$$\n",
    "    u_t + u u_x = 0,\n",
    "$$\n",
    "with the same setup as the linear PDE above.  Try to sse if you can figure out the characteristic curves for the initial conditions\n",
    "$$\n",
    "    u(x,0) = \\left \\{ \\begin{aligned}\n",
    "        4, \\;\\; & x<0, \\\\\n",
    "        3, \\;\\; & x>0.\n",
    "    \\end{aligned} \\right. \n",
    "$$\n",
    "What are the characteristics for this problem?  What ends up being the problem in this case?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The characteristic curves will be\n",
    "$$\n",
    "    x(t) = \\left \\{ \\begin{aligned}\n",
    "        4 t + x_0, \\;\\; & x<0, \\\\\n",
    "        3 t + x_0, \\;\\; & x>0,\n",
    "    \\end{aligned} \\right. \n",
    "$$\n",
    "which lead to characteristic curve collision and therefore a shock.\n",
    "\n",
    "![Shock](./images/characteristics_shock.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We can find the equation for the location of the shock by using the Rankine-Hugoniot condition\n",
    "$$\n",
    "    \\frac{\\text{d} x_s}{\\text{d} t} = \\frac{f(x_+, t) - f(x_-, t)}{u_+ - u_-},\n",
    "$$\n",
    "where $x_s(t)$ is the location of the shock, $f$ is the flux function defined from $u_t + f(u)_x = 0$, and the notation $u_\\pm$ denotes values to the left and right of the shock depending on the sign."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "What if we had the initial condition\n",
    "$$\n",
    "    u(x,0) = \\left \\{ \\begin{aligned}\n",
    "        3, \\;\\; & x<0, \\\\\n",
    "        4, \\;\\; & x>0,\n",
    "    \\end{aligned} \\right.\n",
    "$$\n",
    "instead?  What type of solutions do we now expect?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The characteristic curves are now\n",
    "$$\n",
    "    x(t) = \\left \\{ \\begin{aligned}\n",
    "        3 t + x_0, \\;\\; & x<0, \\\\\n",
    "        4 t + x_0, \\;\\; & x>0,\n",
    "    \\end{aligned} \\right. \n",
    "$$\n",
    "leaving a section of the $(x,t)$ plane without and characteristics. Instead we assume a smoothly varying solution, called a rarefaction, and fill in the missing solution.\n",
    "\n",
    "![Rarefaction](./images/characteristics_rarefaction.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "In this case we want to find smoothly varying solutions inside of this wedge which we find by assuming that the solution is self-similar.  This implies that $u$ is of the form\n",
    "$$\n",
    "    u(x, t) = u(\\xi),\n",
    "$$\n",
    "where for first order hyperbolic PDEs the self-similar variable is\n",
    "$$\n",
    "    \\xi = x / t.\n",
    "$$\n",
    "Try plugging this into the original equation then will give us a solution in terms of $u(\\xi)$ for the interior of the wedge."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "$$\\begin{aligned}\n",
    "    &u_t + u u_x = 0 \\\\\n",
    "    &u' \\frac{\\partial \\xi}{\\partial t} + u u' \\frac{\\partial \\xi}{\\partial x} = 0 \\\\\n",
    "    &u' \\left ( -\\frac{x}{t^2} + u \\frac{1}{t} \\right ) = 0 \\\\\n",
    "    &u' \\left ( -\\xi + u\\right ) = 0\n",
    "\\end{aligned}$$\n",
    "This last equation can now be solved."
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "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.7.6"
  },
  "latex_envs": {
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 0
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}