{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tutorial 6: BVP, PDE, and 3D\n",
    "_Boundary Value Problems, Partial Differential Equations, and how to plot our problems in 3D._\n",
    "\n",
    "That may sound like a lot, however, it is very similar to what you have already learned.\n",
    "If a partiticular section is difficult to understand, try looking at the code line by line."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3D plotting in python\n",
    "Our standard matplotlib cannot handle 3D graphing. Just as before we will `import` the tools we need.\n",
    "\n",
    "```python\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "```\n",
    "We set our Z axes equal to certain data, perhaps a function $z = f(x,y)$\n",
    "```python\n",
    "Z = data\n",
    "```\n",
    "\n",
    "We can use `numpy.meshgrid` to create our meshgrid\n",
    "```python\n",
    "X, Y = np.meshgrid(np.arange(3), np.arange(3))\n",
    "```\n",
    "\n",
    "We create a figure\n",
    "```python\n",
    "fig = plt.figure()\n",
    "```\n",
    "\n",
    "Add our plot to the figure\n",
    "```python\n",
    "ax = fig.add_subplot(1,1,1, projection='3d')\n",
    "```\n",
    "\n",
    "Choose the arrangement\n",
    "```python\n",
    "ax.plot_surface(X, Y, Z)\n",
    "```\n",
    "\n",
    "And that's it!\n",
    "```python\n",
    "plt.show()\n",
    "```\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "\n",
    "data = [[1, 2, 3], [1, 2, 3], [1, 2, 3]]\n",
    "T_array = [1, 2, 3]\n",
    "R_array = [3, 4,5]\n",
    "data = np.array(data)\n",
    "Z = data\n",
    "X, Y = np.meshgrid(np.arange(3), np.arange(3))\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(1,1,1, projection='3d')\n",
    "ax.plot_surface(X, Y, Z)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The toolkit does not strictly prevent you from swapping the axes (or their signs)\n",
    "\n",
    "```python\n",
    "ax.plot_surface(Y, Z, X)\n",
    "ax.plot_surface(-X, -Z, Y)\n",
    "```\n",
    "\n",
    "However, if you are inconsistent, you risk confusing data and graphs. Be especially careful when plotting multiple things onto the same figure.\n",
    "\n",
    "\n",
    "Here two planes will appear on the same subplot, and the X of one is the Y of another. You are not prevented from doing this. \n",
    "```python\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(1,1,1, projection='3d')\n",
    "ax.plot_surface(Y, Z, X)\n",
    "ax.plot_surface(X, Z, Y)\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "\n",
    "data = [[1, 2, 3], [1, 2, 3], [1, 2, 3]]\n",
    "T_array = [1, 2, 3]\n",
    "R_array = [3, 4,5]\n",
    "data = np.array(data)\n",
    "Z = data\n",
    "X, Y = np.meshgrid(np.arange(3), np.arange(3))\n",
    "\n",
    "fig1 = plt.figure()\n",
    "ax1 = fig1.add_subplot(1,1,1, projection='3d')\n",
    "ax1.plot_surface(Y, Z, X)\n",
    "\n",
    "fig2 =  plt.figure()\n",
    "ax2 = fig2.add_subplot(1,1,1, projection='3d')\n",
    "ax2.plot_surface(-X, -Z, Y)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Optionally, just as in 2D, we can add labels:\n",
    "```python\n",
    "ax.set_xlabel('$x$')\n",
    "ax.set_ylabel('$y$')\n",
    "ax.set_zlabel('$z$')\n",
    "```\n",
    "\n",
    "Unfortunately, the $z$ label does not automatically fit within the plot.\n",
    "\n",
    "Furthermore, the labels assume you have chosen the default:\n",
    "```python\n",
    "ax.plot_surface(X, Y, Z)\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "\n",
    "data = [[1, 2, 3], [1, 2, 3], [1, 2, 3]]\n",
    "T_array = [1, 2, 3]\n",
    "R_array = [3, 4,5]\n",
    "data = np.array(data)\n",
    "Z = data\n",
    "X, Y = np.meshgrid(np.arange(3), np.arange(3))\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(1,1,1, projection='3d')\n",
    "ax.plot_surface(X, Y, Z)\n",
    "ax.set_xlabel('$x$')\n",
    "ax.set_ylabel('$y$')\n",
    "ax.set_zlabel('$z$')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Saddle Point Problem:\n",
    "\n",
    "We have a function, $z = f(x,y) = x^2 - y^2$, where $-L < x < L$ and $-L < y < L$, provided $L = 10$, give the 3D plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "# https://en.wikipedia.org/wiki/Saddle_point\n",
    "L = 10.\n",
    "x = np.linspace(-L, L, 100)\n",
    "y = np.linspace(-L, L, 100)\n",
    "X, Y = np.meshgrid(x, y)\n",
    "Z = (X**2 - Y**2)\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(1,1,1, projection='3d')\n",
    "ax.plot_surface(X, Y, Z);\n",
    "# plt.show();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Of course the number of intervals affects the resolution and accuracy:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "# https://en.wikipedia.org/wiki/Saddle_point\n",
    "L = 10.\n",
    "x = np.linspace(-L, L, 8)\n",
    "y = np.linspace(-L, L, 8)\n",
    "X, Y = np.meshgrid(x, y)\n",
    "Z = (X**2 - Y**2)\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(1,1,1, projection='3d')\n",
    "ax.plot_surface(X, Y, Z);\n",
    "# plt.show();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Partial Differential Equations\n",
    "\n",
    "Before we step into the territory of code, let’s find out what `Partial Differential Equation` (PDE) means, and how it differs from `Ordinary Differential Equations` (ODE). \n",
    "\n",
    "Wikipedia explains this well, whereas an ODE “is a differential equation containing one or more functions of one independent variable and its derivatives”, a PDE “may be with respect to _more than_ one independent variable”. \n",
    "\n",
    "In essence, we are doing what we did with ODEs, but simply with _more variables_.\n",
    "\n",
    "As you may very well suspect, we solve PDEs in a similar way to ODEs.\n",
    "\n",
    "```\n",
    "\n",
    "```\n",
    "\n",
    "There is a specific case of differential equations, where we have an Initial Value Problem (IVP); in such a problem we may have the general solution, but we still need a particular solution in order to solve our problem. It is with the Initial Condition(s) (IC) that we do so.\n",
    "\n",
    "In layman’s terms, an IVP is a problem where we have a sense of what the function is (either by having calculated it ourselves or by having been provided the function), and our job is to specify the terms of the function, so it can satisfy some conditions.\n",
    "\n",
    "```\n",
    "\n",
    "```\n",
    "\n",
    "Now let's start coding a problem involving PDE and ICs. This particular problem is called a `Heat Equation`, `Diffusion Equation`, or `Heat Diffusion Equation`, it describes the diffusion of heat over an object as time progresses.\n",
    " \n",
    "\n",
    "We define a function, in this case:\n",
    "```python\n",
    "def calc_heat(nu):\n",
    "    nx = 41\n",
    "    dx = 2 / (nx - 1)\n",
    "    nt = 20    #the number of timesteps we want to calculate\n",
    "    # nu = 0.3   #the value of viscosity\n",
    "    sigma = .2 #sigma is a parameter, we'll learn more about it later\n",
    "    dt = dx**2 #dt is defined using sigma ... more later!\n",
    "\n",
    "\n",
    "    u = np.ones(nx)      #a numpy array with nx elements all equal to 1.\n",
    "    u[int(.5 / dx):int(1 / dx + 1)] = 2  #setting u = 2 between 0.5 and 1 as per our I.C.s\n",
    "    \n",
    "    plt.plot(np.linspace(0, 2, nx), u,'k--');\n",
    "```\n",
    "\n",
    "Our heat diffusion has a range of possible forms it can take on depending on the viscosity\n",
    "```python\n",
    "    un = np.ones(nx) #our placeholder array, un, to advance the solution in time\n",
    "\n",
    "    for n in range(nt):  #iterate through time\n",
    "        un = u.copy() #copy the existing values of u into un\n",
    "        for i in range(1, nx - 1):\n",
    "            u[i] = un[i] + nu * dt / dx**2 * (un[i+1] - 2 * un[i] + un[i-1])\n",
    "            u[0] = 1\n",
    "            u[-1] = 1\n",
    "    plt.plot(np.linspace(0, 2, nx), u,'r--');\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np               #loading our favorite library\n",
    "import matplotlib.pyplot as plt   #and the useful plotting library\n",
    "plt.style.use('ggplot')\n",
    "%matplotlib inline\n",
    "from ipywidgets import interact\n",
    "from IPython.display import clear_output, display, HTML\n",
    "def calc_heat(nu):\n",
    "    nx = 41\n",
    "    dx = 2 / (nx - 1)\n",
    "    nt = 20    #the number of timesteps we want to calculate\n",
    "    # nu = 0.3   #the value of viscosity\n",
    "    sigma = .2 #sigma is a parameter, we'll learn more about it later\n",
    "#     dt = sigma * dx**2 / nu #dt is defined using sigma ... more later!\n",
    "    dt = dx**2 #dt is defined using sigma ... more later!\n",
    "\n",
    "\n",
    "    u = np.ones(nx)      #a numpy array with nx elements all equal to 1.\n",
    "    u[int(.5 / dx):int(1 / dx + 1)] = 2  #setting u = 2 between 0.5 and 1 as per our I.C.s\n",
    "    \n",
    "    plt.plot(np.linspace(0, 2, nx), u,'k--');\n",
    "\n",
    "    un = np.ones(nx) #our placeholder array, un, to advance the solution in time\n",
    "\n",
    "    for n in range(nt):  #iterate through time\n",
    "        un = u.copy() #copy the existing values of u into un\n",
    "        for i in range(1, nx - 1):\n",
    "            u[i] = un[i] + nu * dt / dx**2 * (un[i+1] - 2 * un[i] + un[i-1])\n",
    "            u[0] = 1\n",
    "            u[-1] = 1\n",
    "    plt.plot(np.linspace(0, 2, nx), u,'r--');\n",
    "    plt.xlabel('x')\n",
    "    plt.ylabel('T')\n",
    "\n",
    "    plt.xlim(0,2)\n",
    "    plt.ylim(1,2)\n",
    "    plt.show()\n",
    "    return None\n",
    "interact(calc_heat, nu=(0., .5, 0.02));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is our heat diffusion specifically when $nu = 0.3$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy                 #loading our favorite library\n",
    "from matplotlib import pyplot    #and the useful plotting library\n",
    "%matplotlib inline\n",
    "\n",
    "nx = 41\n",
    "dx = 2 / (nx - 1)\n",
    "nt = 20    #the number of timesteps we want to calculate\n",
    "nu = 0.3   #the value of viscosity\n",
    "sigma = .2 #sigma is a parameter, we'll learn more about it later\n",
    "dt = sigma * dx**2 / nu #dt is defined using sigma ... more later!\n",
    "\n",
    "\n",
    "u = numpy.ones(nx)      #a numpy array with nx elements all equal to 1.\n",
    "u[int(.5 / dx):int(1 / dx + 1)] = 2  #setting u = 2 between 0.5 and 1 as per our I.C.s\n",
    "\n",
    "un = numpy.ones(nx) #our placeholder array, un, to advance the solution in time\n",
    "\n",
    "for n in range(nt):  #iterate through time\n",
    "    un = u.copy() #copy the existing values of u into un\n",
    "    for i in range(1, nx - 1):\n",
    "        u[i] = un[i] + nu * dt / dx**2 * (un[i+1] - 2 * un[i] + un[i-1])\n",
    "        \n",
    "pyplot.plot(numpy.linspace(0, 2, nx), u);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.full_like(x, 673)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PDE and 3D Plotting\n",
    "Let's continue the previous example in 3D.\n",
    "\n",
    "Starting with the setup of our tools:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy\n",
    "from matplotlib import pyplot, cm\n",
    "from mpl_toolkits.mplot3d import Axes3D ##library for 3d projection plots\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we declare our variables:\n",
    "\n",
    "$nt$ our timesteps, $nu$ our viscosity, $\\sigma$, $dx$, and $dt$\n",
    "\n",
    "Using $nx$ and $ny$, the end of our ranges, we define two linspaces:\n",
    "\n",
    "```python\n",
    "x = numpy.linspace(0, 2, nx)\n",
    "y = numpy.linspace(0, 2, ny)\n",
    "```\n",
    "And also our initial conditions (ICs)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "###variable declarations\n",
    "nx = 31\n",
    "ny = 31\n",
    "nt = 17\n",
    "nu = .01\n",
    "dx = 2 / (nx - 1)\n",
    "dy = 2 / (ny - 1)\n",
    "sigma = .25\n",
    "dt = sigma * dx * dy / nu\n",
    "\n",
    "x = numpy.linspace(0, 2, nx)\n",
    "y = numpy.linspace(0, 2, ny)\n",
    "\n",
    "u = numpy.ones((ny, nx))  # create a 1xn vector of 1's\n",
    "un = numpy.ones((ny, nx))\n",
    "\n",
    "###Assign initial conditions\n",
    "# set hat function I.C. : u(.5<=x<=1 && .5<=y<=1 ) is 2\n",
    "u[int(.5 / dy):int(1 / dy + 1),int(.5 / dx):int(1 / dx + 1)] = 2  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With our linspaces we define a meshgrid:\n",
    "\n",
    "```python\n",
    "x = numpy.linspace(0, 2, nx)\n",
    "y = numpy.linspace(0, 2, ny)\n",
    "\n",
    "X, Y = numpy.meshgrid(x, y)\n",
    "```\n",
    "\n",
    "And we graph as before,\n",
    "\n",
    "starting with a figure\n",
    "```python\n",
    "fig = pyplot.figure()\n",
    "```\n",
    "\n",
    "Add our plot to the figure\n",
    "```python\n",
    "ax = fig.gca(projection='3d')\n",
    "```\n",
    "\n",
    "Choose the arrangement\n",
    "```python\n",
    "surf = ax.plot_surface(X, Y, u, rstride=1, cstride=1, cmap=cm.viridis, linewidth=0, antialiased=False)\n",
    "\n",
    "ax.set_xlim(0, 2)\n",
    "ax.set_ylim(0, 2)\n",
    "ax.set_zlim(1, 2.5)\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = pyplot.figure()\n",
    "ax = fig.gca(projection='3d')\n",
    "X, Y = numpy.meshgrid(x, y)\n",
    "surf = ax.plot_surface(X, Y, u, rstride=1, cstride=1, cmap=cm.viridis,\n",
    "        linewidth=0, antialiased=False)\n",
    "\n",
    "ax.set_xlim(0, 2)\n",
    "ax.set_ylim(0, 2)\n",
    "ax.set_zlim(1, 2.5)\n",
    "\n",
    "ax.set_xlabel('$x$')\n",
    "ax.set_ylabel('$y$');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Just like in the 2D example, we are going to simulate heat diffusion, there is a lot here, so let's go through this step by step.\n",
    "\n",
    "We define our diffuse function, which takes in the timestep as a parameter:\n",
    "```python\n",
    "def diffuse(nt):\n",
    "```\n",
    "We define a 2D array $u$ which is defined by the derivatives of $y$ and $x$\n",
    "```python\n",
    "    u[int(.5 / dy):int(1 / dy + 1),int(.5 / dx):int(1 / dx + 1)] = 2  \n",
    "```\n",
    "\n",
    "We iterate through time, and copy the existing values of $u$ into $un$, just like before. However, this time we have an extra dimension, so we iterate over every _$i$_ and _$j$_ within $u$.\n",
    "```python\n",
    "    for n in range(nt + 1):\n",
    "        un = u.copy()\n",
    "        for i in range(1, u.shape[0]-1):\n",
    "            for j in range(1, u.shape[1] - 1):\n",
    "                u[i,j] = un[i, j] + nu * dt / dx**2 * \\\n",
    "                (un[i+1, j] - 2*(un[i,j]) + un[i-1,j]) + \\\n",
    "                nu * dt / dy**2 * \\\n",
    "                (un[i,j+1] - 2*un[i, j] + un[i, j-1])\n",
    "```\n",
    "                \n",
    "Here we are setting the outer square of our 2D array to 1, convince yourself that this is correct:\n",
    "```python\n",
    "    u[0, :] = 1\n",
    "    u[-1, :] = 1\n",
    "    u[:, 0] = 1\n",
    "    u[:, -1] = 1\n",
    "```\n",
    "\n",
    "And now we plot the 3D graph.\n",
    "\n",
    "    -Starting with a figure\n",
    "    \n",
    "    -Add our plot to the figure\n",
    "    \n",
    "    -Choose the arrangement\n",
    "\n",
    "```python\n",
    "    fig = pyplot.figure()\n",
    "    ax = fig.gca(projection='3d')\n",
    "    surf = ax.plot_surface(X, Y, u[:], rstride=1, cstride=1, cmap=cm.viridis,\n",
    "        linewidth=0, antialiased=True)\n",
    "    ax.set_zlim(1, 2.5)\n",
    "    ax.set_xlabel('$x$')\n",
    "    ax.set_ylabel('$y$');\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "###Run through nt timesteps\n",
    "def diffuse(nt):\n",
    "    u[int(.5 / dy):int(1 / dy + 1),int(.5 / dx):int(1 / dx + 1)] = 2  \n",
    "    \n",
    "#     un+1i,j=uni,j+νΔtΔx2(uni+1,j−2uni,j+uni−1,j)+νΔtΔy2(uni,j+1−2uni,j+uni,j−1)\n",
    "#     u[1:-1, 1:-1] = (un[1:-1,1:-1] + \n",
    "#                 nu * dt / dx**2 * \n",
    "#                 (un[1:-1, 2:] - 2 * un[1:-1, 1:-1] + un[1:-1, 0:-2]) +\n",
    "#                 nu * dt / dy**2 * \n",
    "#                 (un[2:,1: -1] - 2 * un[1:-1, 1:-1] + un[0:-2, 1:-1]))\n",
    "\n",
    "    for n in range(nt + 1):\n",
    "        un = u.copy()\n",
    "        for i in range(1, u.shape[0]-1):\n",
    "            for j in range(1, u.shape[1] - 1):\n",
    "                u[i,j] = un[i, j] + nu * dt / dx**2 * \\\n",
    "                (un[i+1, j] - 2*(un[i,j]) + un[i-1,j]) + \\\n",
    "                nu * dt / dy**2 * \\\n",
    "                (un[i,j+1] - 2*un[i, j] + un[i, j-1])\n",
    "        u[0, :] = 1\n",
    "        u[-1, :] = 1\n",
    "        u[:, 0] = 1\n",
    "        u[:, -1] = 1\n",
    "\n",
    "    \n",
    "    fig = pyplot.figure()\n",
    "    ax = fig.gca(projection='3d')\n",
    "    surf = ax.plot_surface(X, Y, u[:], rstride=1, cstride=1, cmap=cm.viridis,\n",
    "        linewidth=0, antialiased=True)\n",
    "    ax.set_zlim(1, 2.5)\n",
    "    ax.set_xlabel('$x$')\n",
    "    ax.set_ylabel('$y$');\n",
    "    \n",
    "# http://nbviewer.jupyter.org/github/barbagroup/CFDPython/blob/master/lessons/09_Step_7.ipynb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we plug in 10 timesteps into the diffuse function we defined before:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "diffuse(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Boundary value problems\n",
    "\n",
    "Let’s take another moment before we start coding again. What is a Boundary Value Problem (BVP)? Well, it is a type of problem where we must satisfy certain boundary conditions.\n",
    "If this sounds familiar, it should be. Our old friend Wikipedia comes to save the day again, “a [BVP] has conditions specified at the extremes (‘boundaries’) of the independent variable in the equation whereas an [IVP] has all of the conditions specified at the same value of the independent variable”.\n",
    "\n",
    "```\n",
    "\n",
    "```\n",
    "\n",
    "Imagine we have some problem involving the height of a skydiver\n",
    "\n",
    "an IVP would define a function $h(t)$,\n",
    "\n",
    "perhaps involving $h'(t) = k \\cdot h(t)$ and $h(0) = C$, where k and C are specified\n",
    "\n",
    "Here the we only ever use the value $t_i = 0$ in $h(0)$\n",
    "\n",
    "```\n",
    "\n",
    "```\n",
    "\n",
    "How about a BVP? Well, here there has been defined two `boundaries`, such as the initial height $h(0) = h_1$ and the final height $h(t_f) = h_2$. \n",
    "\n",
    "We would need the differential equation as well. However, you can see here we use $t_i$ and $t_f$. These were the aforementioned \"extremes\", before the skydiver jumps (the highest possible height) and after they land (the lowest possible height).\n",
    "\n",
    "For our `heat equations`, these boundaries could be end points on an object. There are multitude of problems that need boundary conditions in order to solve.\n",
    "```\n",
    "\n",
    "```\n",
    "\n",
    "\n",
    "We have a few new tools here as well:\n",
    "```python\n",
    "from scipy.optimize import least_squares\n",
    "from scipy.integrate import ode\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib\n",
    "# import snips as snp\n",
    "# snp.prettyplot(matplotlib)\n",
    "from scipy.optimize import least_squares\n",
    "from scipy.integrate import ode\n",
    "from scipy.integrate import odeint\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The method we will be using is solve by shooting.\n",
    "\n",
    "We will test our ODE function, \"odefun\", and see if overshoots our initial boundary condition or undershoots it.  \n",
    "\n",
    "Provided that we know a limit $u_1(0)$, our miss-shot can give us an $u_2(0)$ using fsolve.\n",
    "\n",
    "In this case we will guess $u_2(0) = 1$\n",
    "\n",
    "Our odefun is:\n",
    "\n",
    "```python\n",
    "def odefun(U, y):\n",
    "    u1, u2 = U\n",
    "    mu = 1\n",
    "    Pdrop = -100\n",
    "    du1dy = u2\n",
    "    du2dy = 1.0 / mu * Pdrop\n",
    "    return [du1dy, du2dy]\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We need u_1(0) and u_2(0), but we only have u_1(0). We need to guess a value for u_2(0) \n",
    "# and see if the solution goes through the u_2(d)=0 boundary value.\n",
    "\n",
    "d = 0.1 # plate thickness\n",
    "\n",
    "def odefun(U, y):\n",
    "    u1, u2 = U\n",
    "    mu = 1\n",
    "    Pdrop = -100\n",
    "    du1dy = u2\n",
    "    du2dy = 1.0 / mu * Pdrop\n",
    "    return [du1dy, du2dy]\n",
    "\n",
    "u1_0 = 0 # known\n",
    "u2_0 = 1 # guessed\n",
    "\n",
    "dspan = np.linspace(0, d)\n",
    "\n",
    "U = odeint(odefun, [u1_0, u2_0], dspan)\n",
    "\n",
    "plt.plot(dspan, U[:,0])\n",
    "plt.plot([d],[0], 'ro')\n",
    "plt.xlabel('d')\n",
    "plt.ylabel('$u_1$');\n",
    "# plt.savefig('images/bvp-shooting-1.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And so we see here that we have undershot, we want to overshoot for our range, so we try again. This time with $u_2(0) = 7$.\n",
    "\n",
    "You can interact with the $u_2(0)$ here to see how it affects our odefun"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ipywidgets import interact, interactive, fixed, interact_manual\n",
    "\n",
    "def calc_poiss(u2_0):\n",
    "    u1_0 = 0 # known\n",
    "#   u2_0 = 7 # guessed\n",
    "    dspan = np.linspace(0, d)\n",
    "    U = odeint(odefun, [u1_0, u2_0], dspan)\n",
    "    plt.plot(dspan, U[:,0])\n",
    "    plt.plot([d],[0], 'ro')\n",
    "    plt.xlabel('d')\n",
    "    plt.ylabel('$u_1$')\n",
    "    \n",
    "interact(calc_poiss, u2_0=(0., 7, .2));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the previous slider, you can see that we satisfy our boundary condition when $u_2(0) \\simeq 5$, but we want to know how to calculate this number, and we want something more precise.\n",
    "\n",
    "as mentionned earlier, since we know $1 < u_2(0) < 7$ we can use fsolve to figure this out for us. Luckily that comes with scipy\n",
    "\n",
    "```python\n",
    "from scipy.optimize import fsolve\n",
    "```\n",
    "And we have our known initial condition, and the size of d spans from 0 to 0.1\n",
    "```python\n",
    "u1_0 = 0 # known\n",
    "dspan = np.linspace(0, d)\n",
    "```\n",
    "\n",
    "Now we need to give fsolve an function to solve. This is our objective function.\n",
    "```python\n",
    "def objective(u2_0):\n",
    "    U = odeint(odefun, [u1_0, u2_0], dspan)\n",
    "    u1 = U[:,0]\n",
    "    return u1[-1]\n",
    "```\n",
    "\n",
    "Now we set $u_2(0)$ equal to what we determined it should be.\n",
    "```python\n",
    "u2_0, = fsolve(objective, 1.0)\n",
    "```\n",
    "\n",
    "And using odeint as before, we get our Numerical solution:\n",
    "```python\n",
    "U = odeint(odefun, [u1_0, u2_0], dspan)\n",
    "\n",
    "plt.plot(dspan, U[:,0], label='Numerical solution')\n",
    "plt.plot([d],[0], 'ro')\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.integrate import odeint\n",
    "from scipy.optimize import fsolve\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "d = 0.1 # plate thickness\n",
    "Pdrop = -100\n",
    "mu = 1\n",
    "\n",
    "def odefun(U, y):\n",
    "    u1, u2 = U\n",
    "    du1dy = u2\n",
    "    du2dy = 1.0 / mu * Pdrop\n",
    "    return [du1dy, du2dy]\n",
    "\n",
    "u1_0 = 0 # known\n",
    "dspan = np.linspace(0, d)\n",
    "\n",
    "def objective(u2_0):\n",
    "    U = odeint(odefun, [u1_0, u2_0], dspan)\n",
    "    u1 = U[:,0]\n",
    "    return u1[-1]\n",
    "\n",
    "u2_0, = fsolve(objective, 1.0)\n",
    "\n",
    "# now solve with optimal u2_0\n",
    "U = odeint(odefun, [u1_0, u2_0], dspan)\n",
    "\n",
    "plt.plot(dspan, U[:,0], label='Numerical solution')\n",
    "plt.plot([d],[0], 'ro')\n",
    "\n",
    "# plot an analytical solution\n",
    "u = -(Pdrop) * d**2 / 2 / mu * (dspan / d - (dspan / d)**2)\n",
    "plt.plot(dspan, u, 'r--', label='Analytical solution')\n",
    "print(u2_0)\n",
    "\n",
    "plt.xlabel('d')\n",
    "plt.ylabel('$u_1$')\n",
    "plt.legend(loc='best');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bibliography:\n",
    "\n",
    "https://en.wikipedia.org/wiki/Ordinary_differential_equation\n",
    "\n",
    "https://en.wikipedia.org/wiki/Partial_differential_equation\n",
    "\n",
    "https://en.wikipedia.org/wiki/Boundary_value_problem\n",
    "\n",
    "http://nbviewer.jupyter.org/github/barbagroup/CFDPython/blob/master/lessons/09_Step_7.ipynb\n",
    "\n",
    "http://kitchingroup.cheme.cmu.edu/blog/2013/02/15/Plane-Poiseuille-flow-BVP-solve-by-shooting-method/"
   ]
  }
 ],
 "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}