02_Keep_it_steady.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Keep it Steady\n",
    "\n",
    "_Solving steady-state one-dimensional differential equations with the finite difference approximations._"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Derivatives in Engineering equations\n",
    "\n",
    "This Module is focused on solving engineering problems where the governing differential equations are across space and time. The general form of these differential equations are [Partial Differential Equations (check out the 3Blue1Brown video)](https://www.youtube.com/watch?v=ly4S0oi3Yz8). A partial differential equation (PDE) relates the derivatives in space and time. The main difference between _ordinary_ and _partial_ differential equations (ODEs and PDEs) is that ODEs only have one indendent variable e.g. all derivatives are $\\frac{d^i}{dx^i}$. In PDEs, there are multiple derivatives, sometimes time and space. These equations come up in a number of engineering applications e.g. \n",
    "\n",
    "* [Conductive heat transfer](https://en.wikipedia.org/wiki/Heat_equation) $q(x,y,t)+\\rho C_{p}\\frac{\\partial T}{\\partial t} = -\\kappa \\nabla^2T$\n",
    "* [Kirchoff plate mechanics](https://en.wikipedia.org/wiki/Kirchhoff%E2%80%93Love_plate_theory) $D\\nabla^2 \\nabla^2 w(x,y) = -q(x,y,t)-2\\rho h \\frac{\\partial w^2}{\\partial t^2}$\n",
    "* [Euler Beam mechanics](https://en.wikipedia.org/wiki/Euler%E2%80%93Bernoulli_beam_theory) $EI \\frac{\\partial^4 w}{\\partial x^4} = -\\mu \\frac{\\partial^2 w}{\\partial t^2}+q(x,t)$\n",
    "* [Elastica Buckling equation](https://www.continuummechanics.org/columnbuckling.html) $EI\\frac{d^2y}{dx^2}=-Py$\n",
    "* ...\n",
    "\n",
    "## Discussion\n",
    "\n",
    "What other PDE equations describe engineering problems?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "# One dimensional, steady-state differential equations\n",
    "\n",
    "We can start our discussion of these differential equations with the steady-state solutions along one axis. These analyses include static deflection and buckling of beams and steady-state temperature through a pipe or heat fin. \n",
    "\n",
    "1. Steady-state Temperature: $-\\kappa A \\frac{d^2T}{dx^2}=Q(x)$\n",
    "\n",
    "2. Static beam deflection: $EI \\frac{\\partial^4 w}{\\partial x^4} = q(x)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Temperature in a fin meant to remove heat to enviroment\n",
    "\n",
    "<img src=\"../images/thermal_fin.png\" style=\"width: 500px;\"/> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the heat fin shown above, a hot substrate is at temperature $T_{base}$ and we want to remove heat with a cylindrical fin that is connected to the substrate. We know the base of the fin is the same temperature as the substrate, $T(x=0)=T_{0}=T_{base}$, but we do not know the heat flux into the fin, $\\frac{dT}{dx}(0)=??$. We also do not know the temperature at the end of the fin, $T(L)=??$. So we could use our first [shooting method](https://github.uconn.edu/rcc02007/CompMech03-IVPs/blob/master/notebooks/04_Getting_to_the_root.ipynb), but we also have a new approach. We use our new knowledge of [Linear Algebra](https://github.uconn.edu/rcc02007/CompMech04-LinearAlgebra) and [finite differences](./01_Revisiting_derivatives.ipynb) to create a linear system of equations. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is what we do know:\n",
    "\n",
    "__Equations:__\n",
    "\n",
    "$\\frac{d^2T}{dx^2} \\approx \\frac{T_{i-1}-2T_{i}+T_{i+1}}{\\Delta x^2}$\n",
    "\n",
    "$-\\kappa \\pi R^2 \\frac{T_{i-1}-2T_{i}+T_{i+1}}{\\Delta x^2} = h2\\pi R (T_{i}-T_{\\infty})$\n",
    "\n",
    "where , R is the radius of the cylinderical fin, $\\Delta x$ is the step size, $\\kappa$ is the conductive heat transfer of the fin, $h$ is the convective heat transfer between the fin and the environment. \n",
    "\n",
    "__Boundary Conditions:__\n",
    "\n",
    "The base of the fin is the same temperature as the base, $T(x=0) = T_{base}=T_{0}.$\n",
    "\n",
    "The tip of the fin has the same convective heat transfer as the rest of the fin, so the energy conducted to the tip has to be removed through convection, $hA(T(L)-T_{\\infty})=-\\kappa A\\frac{dT}{dx} \\rightarrow h(T_{6}-T_{\\infty})=-\\kappa \\frac{T_{7}-T_{5}}{\\Delta x}.$ \n",
    "\n",
    "1. $T_{0} = T_{base}$\n",
    "\n",
    "2. $T_7=\\frac{-h\\Delta x}{k}(T_{6}-T_{\\infty})+T_5$\n",
    "\n",
    "If we divide the fin into 6 equally spaced sections, we have 6 equations with 6 unknown temperatures based upon finite difference equations as such,\n",
    "\n",
    "1. $T_0-2T_1+T_2+\\Delta x^2 h'(T_{\\infty}-T_1) = 0$\n",
    "\n",
    "2. $T_1-2T_2+T_3+\\Delta x^2 h'(T_{\\infty}-T_2) = 0$\n",
    "\n",
    "3. $T_2-2T_3+T_3+\\Delta x^2 h'(T_{\\infty}-T_3) = 0$\n",
    "\n",
    "4. $T_3-2T_4+T_4+\\Delta x^2 h'(T_{\\infty}-T_4) = 0$\n",
    "\n",
    "5. $T_4-2T_5+T_5+\\Delta x^2 h'(T_{\\infty}-T_5) = 0$\n",
    "\n",
    "6. $T_5-2T_6+T_7+\\Delta x^2 h'(T_{\\infty}-T_6) = 0 \\leftarrow T_7 = \\frac{-h\\Delta x}{k}(T_{6}-T_{\\infty})+T_5$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where $h'=\\frac{2h}{\\kappa R}$ is the modified convective heat transfer for the fin. And our boundary conditions give us values for $T_{0}~and~T_{7}.$ We can plug in constants for forced air convection, $h=100~W/m^2K$, aluminum fin, $\\kappa=200~W/mK$, and 60-mm-long and 1-mm-radius fin, the air is room temperature, $T_{\\infty}=20^oC$, and the base is $T_{base}=T_{0}=100^oC$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams.update({'font.size': 22})\n",
    "plt.rcParams['lines.linewidth'] = 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "h=100 # W/m/m/K\n",
    "k=200 # W/m/K\n",
    "R=1E-3# radius in m\n",
    "L=60E-3# length in m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\left[\\begin{array}{cccccc}\n",
    "2.1 & -1 & 0 & 0 & 0 & 0 \\\\\n",
    "-1 & 2.1 & -1 & 0 & 0 & 0 \\\\\n",
    "0 & -1 & 2.1 & -1 & 0 & 0 \\\\\n",
    "0 & 0 & -1 & 2.1 & -1 & 0 \\\\\n",
    "0 & 0 & 0 & -1 & 2.1 & -1 \\\\\n",
    "0 & 0 & 0 & 0 & -2 & 2.105 \\\\\n",
    "\\end{array}\\right]\n",
    "\\left[\\begin{array}{c}\n",
    "T_{1}\\\\\n",
    "T_{2}\\\\\n",
    "T_{3}\\\\\n",
    "T_{4}\\\\\n",
    "T_{5}\\\\\n",
    "T_{6}\\end{array}\\right]=\n",
    "\\left[\\begin{array}{c}\n",
    "10+T_0\\\\\n",
    "10\\\\\n",
    "10\\\\\n",
    "10\\\\\n",
    "10\\\\\n",
    "10+h\\Delta x/k T_{\\infty}\\\\\n",
    "\\end{array}\\right]$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "h' = 1000.0, and step size dx= 0.01\n",
      "2.1\n",
      "finite difference A:\n",
      "------------------\n",
      "[[ 2.1   -1.     0.     0.     0.     0.   ]\n",
      " [-1.     2.1   -1.     0.     0.     0.   ]\n",
      " [ 0.    -1.     2.1   -1.     0.     0.   ]\n",
      " [ 0.     0.    -1.     2.1   -1.     0.   ]\n",
      " [ 0.     0.     0.    -1.     2.1   -1.   ]\n",
      " [ 0.     0.     0.     0.    -2.     2.105]]\n",
      "\n",
      "finite difference b:\n",
      "------------------\n",
      "[102.    2.    2.    2.    2.    2.1]\n",
      "\n",
      "finite difference solution T(x):\n",
      "------------------\n",
      "[79.5141801  64.97977822 54.94335415 48.4012655  44.6993034  43.46727164]\n",
      "\n",
      "finite difference solution at x (mm)=\n",
      "------------------\n",
      "[10. 20. 30. 40. 50. 60.]\n"
     ]
    }
   ],
   "source": [
    "hp = 2*h/k/R\n",
    "N=6\n",
    "dx=L/N\n",
    "\n",
    "print('h\\' = {}, and step size dx= {}'.format(hp,dx))\n",
    "diag_factor=2+hp*dx**2 # diagonal multiplying factor\n",
    "print(diag_factor)\n",
    "Tinfty=20\n",
    "T0 = 100\n",
    "\n",
    "A = np.diag(np.ones(N)*diag_factor)-np.diag(np.ones(N-1),-1)-np.diag(np.ones(N-1),1)\n",
    "A[-1,-2]+= -1\n",
    "A[-1,-1]+= h/k*dx\n",
    "b = np.ones(N)*hp*Tinfty*dx**2\n",
    "b[0]+=T0\n",
    "b[-1]+=h*dx/k*(Tinfty)\n",
    "print('finite difference A:\\n------------------')\n",
    "print(A)\n",
    "print('\\nfinite difference b:\\n------------------')\n",
    "print(b)\n",
    "T=np.linalg.solve(A,b)\n",
    "print('\\nfinite difference solution T(x):\\n------------------')\n",
    "print(T)\n",
    "print('\\nfinite difference solution at x (mm)=\\n------------------')\n",
    "print(np.arange(1,7)*dx*1000)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Observing our set of equations\n",
    "\n",
    "The equations and solution are shown above for the finite difference set of six coupled linear equations. We solved for the solutions of temperature at $x=[10,~20,~30,~40,~50,~60]~mm$. We didn't include a solution for $x=0~mm$ because it was included in our boundary conditions. If we specified the temperature at \n",
    "$x=60~mm$, then we would exclude that point as well. \n",
    "\n",
    "When we approximate a derivative as $\\frac{df}{dx}\\approx \\frac{\\Delta f}{\\Delta x}$, what we're doing is changing a differential equation into a set of coupled equations. \n",
    "\n",
    "Below, lets compare our solution to the [analytical solution](https://en.wikipedia.org/wiki/Fin_(extended_surface)#Solutions) for a fin with the same boundary conditions\n",
    "\n",
    "$T(x)=20+80\\frac{\\cosh(s(L-x))+\\frac{h}{sk}\\sinh(s(L-x))}{\\cosh(sL)+\\frac{h}{sk}\\sinh(sL)}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAscAAAEaCAYAAAD9vvpCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd3xT1fsH8M/TpDvde9GWFsooFNqyUWQKFFQoSwQFnOzx9Su4+ArqV0F/IiiCbAfIRmQIomUjdDCkjEIZLauli+6V5vz+SNJvSNORkjYtPO/X677S3HvOvU9SrE9OnnsOCSHAGGOMMcYYA0yMHQBjjDHGGGMNBSfHjDHGGGOMqXByzBhjjDHGmAonx4wxxhhjjKlwcswYY4wxxpiK1NgBNHbOzs7Cz8/P2GEwxlijEhcXly6EcDF2HIwxpo2T40fk5+eH2NhYY4fBGGONChElGTsGxhjThcsqGGOMMcYYU2kwyTERBRHRdCL6mYguE5GCiAQRDatB39FEdJSIsokoj4hiiWgyEVX5+oioPxH9QUSZRFRARPFE9D4RmRvulTHGGGOMscaiIZVVTAQwXd9ORLQUwCQARQD+AlAKoDeAbwH0JqLhQogyHf3eAbAAQBmAQwCyAPQA8AmAQUTUWwhRULuXwhhjjDHGGqMGM3IMIB7AFwBGAggEcLi6DkQUCWVinAKgrRBikBBiCIBmAC4BGAJgio5+4QA+B1AAoJsQoo8QYjiApgCOAOgM4FNDvCjGGGOMMdZ4NJjkWAixSgjxjhBisxDiWg27vat6nC2EuKpxrlQoR6IBYI6O8oo5AAjAAiHEKY1+eQDGA1AAmERE9rV5LYwxxhhjrHFqMMmxvojIG0AYgBIAW7SPCyEOA7gDwB3KkWB1PzMAA1RP1+vodx3A3wDMAAw0eOAA4O4OEFXc3N3r5HKMMcYYY6xmGm1yDKC96vGCEKKwkjYxWm0BIAiAFYDMKkaodfUznNRU/fYzxhhjjLF60ZiTY3/VY1VzZSZrtdX8ORmV09WPMcYYY4w95hpzcixTPeZX0SZP9WhjgH7liOgN1XRxsWlpadUGyhhjjDHGGofGnByT6lHUU79yQogVQohwIUS4iwuvfsoYY4wx9rhozMlxrupRVkUb9bFcjX217ccYY4wxxh5zjTk5vql69K2ijY9WW82fm+jZz3Dc3HTuzrJxRIlcUSeXZIwxxhhj1WtIK+Tp64zqsTURWVYyY0UHrbYAcBlAIQBHIgqoZMaKjjr6GU5KSvmPF+/m4IWlx1FSpkyKJ/15Be/0b1Enl2WMMcYYY1VrtCPHQohbAE5DOR/xcO3jRNQDgDeUq+f9rdGvBMDvqqcv6ejXFEAXKOdP3mPwwLW08rTFO/2Dyp8vO3wNp65n1PVlGWOMMcaYDo02OVb5TPW4gIgC1TuJyBXAd6qnnwshtGsVPofyhrzZRNRRo58MwBoo35fvhBAP6ixyDRO6+aN7oDMAQAhg5qazyC4srY9LM8YYY4wxDQ0mOSaiUCI6qd4AhKoO/VdrfzkhxFYAy6BcBe88Ee0iou0ArgJoBeBXAN9qX0sIEQPlEtJWAE4Q0R9EtBnANQA9AJwC8H7dvNKKTEwIXw4Pgb2VKQDgbnYRPvg1HkLUekINxhhjjDFWCw0mOQZgC6CTxqaeY7iZ1v6HCCEmQVkecRrKxPZZAIkApgCIFEKU6bqYEGIhlMtIH4SyNnkwgHQAHwDoIYQoMNQLqwl3Owt8PrRN+fNd5+7i17N36jMExhhjjLEnHvHo5KMJDw8XsbGxBjvf7K3/YFPsLQCAjbkUe6c/BR9HK4OdnzHGGgIiihNChBs7DsYY09aQRo4ZgLmDW8HPSZkM5xbLMXPTWcjLeHo3xhhjjLH6oHdyTEQyIupCRC8Q0Tgiel713LouAnzSWJtL8fWo9pCYKBfyi03KwrJDumabY4wxxhhjhlaj5JiInIloNhGdAJAJ4BiAbQBWA9iuep5JRCeI6B0icq6ziJ8A7XzsMaN3s/LnX/91FWeSs4wYEWOMMcbYk6HK5JiIvIhoDYBbAP4LoDOAUgDnARwGsEv1GA9Arjr+OYBbRLSaiLzqMPbH2qSegejg5wAAKFMIzNh0FvnFciNHxRhjjDH2eKs0OSaiuQASAIwDcB3Ae1DO6mArhGgnhOglhHhB9RgCwA7KleXeh3LZ5fEAElTnYXqSmBC+GtEONubKRQyTMgowb9cFI0fFGGOMMfZ4q2rk+CMARwF0FUK0FkIsEELEVTE1mlwIESuE+EwI0RJAdyjLLf5j8KifED6OVvj4heDy55tjb+P38/eMGBFjjDHG2OOtquT4aSHEACHEySraVEoIcUII0R/KuYdZLb3Q3gvPhXiWP5+z/TxSsouMGBFjjDHG2OOr0uRYCHHMEBcw1HmeZB+/EAwve0sAQHZhKf615SwUCp6fmjHGGGPM0Hie40bAztIUX40IASlnd8PxxAysPnbDuEExxhhjjD2GqputYhwRzSeialcxIqIOqrYvGy48ptapqRMmPRNQ/vyL/Qm4cDfbiBExxhhjjD1+qpqtohmAlQCeA3CuBuc6B2AwgJVE5GeI4NjDZvRpjrbedgCAkjIFpm88i6JSnfdHMsYYY4yxWqhq5PgV1fH3hBCl1Z1ICFEC4F0AplBO/8YMzFRigq9HtoOlqQQAkHg/D5/tvWTkqBhjjDHGHh9VJce9AKQLIfbW9GRCiH0A7gPo86iBMd2ausgwd3Cr8uc//J2Eg5fvGzEixhhjjLHHR1XJcXMAsbU4ZxyAFrULh9XEqA4+6NvKrfz5v7eeQ3pesREjYowxxhh7PFSVHNsCyKzFObMA2NQuHFYTRIQFkW3hYmMOAEjPK8E7W/+BEDy9G2OMMcbYo6gqOc4B4FCLc9oDyKtdOKymHK3N8H/DQ8qfR12+j59PJRsxIsYYY4yxxq+q5DgJQIdanLMDgJu1iobp5enmLpjQzb/8+Se7LyLxfq4RI2KMMcYYa9yqSo4PAnAmouE1PRkRjQDgourL6sE7/YPQwl1ZxVIsV2DaL2dRLOfp3RhjjDHGaqOq5Hg5AAWA74iodXUnIqJgAN8BkAP43jDhsepYmErw9ah2MJMqf5UX7+Xgqz+uGDkqxhhjjLHGqdLkWAiRCGAhACcA0arV7wK12xFRIBF9DOAklDXK/yeEuFpXAbOKWrjbYk7//00QsuLodZxITDdiRIwxxhhjjVOVy0cD+ADKUWBLAO8DSCCiTCJKUG2ZABIAvAfACsAqIcS7dRox02lcVz883dwFACAEMGvzOTwoKDFyVIwxxhhjjUuVybFQmghgKIAzAAjK2SiaqTZ71b4zACKFEG/WbbisMiYmhC+HtYWjtRkAICWnCO/tOM/TuzHGGGOM6aG6kWMAgBDiVyFEOABvAM8DeE21PQ/AWwgRLoTYUXdhsppwtbXA50PblD/fez4FW+NuGzEixhhjjLHGRapPYyHEXQB36ygWZgD9WrvjxY5N8Eu0cs7jj367gI7+jvB1sjZyZIwxxhhjDV+NRo5Z4/LhoJZo6qxMhvNLyjBj01nIyxRGjooxxhhjrOGrNDkmotGGuIChzsNqzspMisWj2kNqQgCAM8kP8E1UopGjYowxxhhr+KoaOf6ZiM4SUSQR6VV+QURSIhpGROcA/PRoIbLaaONth1n9mpc//ybqKuKSMo0YEWOMMcZYw1dVcvw6AHcAmwHcJaKviWgwEbnoakxErkT0HBEtgbIueRMAVyhv3GNG8ObTAejo7wgAUAhgxqazyC0qNXJUjDHGGGMNV1WLgKyGcrq2Bap20wD8CiCFiNKI6DIRRase0wDcA7ADwBTVKf4LoLkQYm2dvgJWKYkJYdHIdrCxUA7838osxEe/XTRyVIwxxhhjDVd18xznCiHeg3IKt1cB7AbwAMpV85oDCFc9OgHIhDI5fhnK6d0+FELk1mHsrAa87C3x6ZD/Te+27fRt7P6HJxxhjDHGGNOlRrXEQogiAGsBrCUiAuALZcmEHZTJ8n0hRFKdRckeyXMhnjh0+T62n7kDAHhv+3mENnGAp72lkSNjjDHGGGtY9J7KTbVq3k0hRLQQ4oAQIoYT44Zv3vOt4e2gTIZziuSYtfksyhS8eh5jjDHGmCae5/gJYWNhiq9HtoNqdjecvJ6JlUevGzcoxhhjjLEGhpPjJ0i4nyOm9Awsf/5/fyQg/k62ESNijDHGGGtYODl+wkzt3QwhPvYAgNIygWkbz6CwpMzIUTHGGGOMNQycHD9hTCUmWDyyHazMJACA62n5+GQPT+/GGGOMMQZwcvxE8nO2xkeDW5c/X38qGX9eTDViRIwxxhhjDQMnx0+o4eHeGBDsXv78nW3/4H5ukREjYowxxhgzPk6On1BEhP8OaQM3W3MAQGZ+Cf695R8IwdO7McYYY+zJVaNFQNjjycHaDF+NaIeXVp0CABy+koYfTtzEuG7+Ro6MMcYeFhcXRwB6SqXS4UT0lBBCZuyYGGONBxHlCSGOyuXyLQAOhoWFVToaSPqOFBKRJYARALoAcAFwWAixRHXMD4AngFghREntwm9cwsPDRWxsrLHDeCSf7rmIlUdvAADMpCbYPbU7mrvZGDkqxtjjjIjihBDhNWkbFxdHJiYmsy0tLV9zcXERtra2eVKptEy5YCtjjFVNCAG5XC7JycmRpaWlUWFh4SqFQrGgsgRZr7IKIuoB4DqANQDeAPACgHYaTboDOApgcO3CZ8bw9rNBaOlhCwAokSsw7ZczKJbz9G6MsQajp6Wl5WuBgYE5Tk5O2aamppwYM8ZqjIhgampa5uTklB0YGJhjaWn5GoCelbWvcXJMREEAdgNwBfAjgFcBaP91+hVAMZRJM2skzKUSLBnVDuZS5T+Hyym5+GJfgpGjYowxJalUOtzFxUVIpVKFsWNhjDVuUqlU4eLiAqlUOryyNvqMHL8PwArAS0KI8UKItdoNhBB5AC4BaK93tMyomrnZ4P2IluXPVx27gaNX04wYEWOMKRHRU7a2tnnGjoMx9niwtbXNI6KnKjuuT3LcC8A/QoiN1bS7BcBDj/OyBmJsZ1/0DHIpf/6vzeeQlf9ElI4zxhowIYRMKpVyrRdjzCCkUqm8qpt69UmOXQDU5Lt2OZQjzKyRISIsHBYCJ2szAMD93GLM2c7TuzHGjI9rjBljhlLd3xN9kuMHALxq0K4pAP4+vpFysTHHwmFty5/vv5CKzbG3jBgRY4wxxlj90Sc5jgUQrpquTSciCgEQAuDEo4XFjKl3SzeM7exb/vyj3y7iRnq+ESNijDHGGKsf+iTH3wMwB7CZiJpoHyQiDwCrNdqyRuy9gS0R4GINACgsLcOMjWdQWsY3ijPGGGPs8Vbj5FgI8RuU8xuHA0gkIvXo8DNE9BeAawBCAawQQhw0eKSsXlmaSbB4VHuYSpR1OeduZ2Pxn1eNHBVjjLHGyMvLqw0RhSUkJJjVx/WIKIyIwurjWtoSEhLMiCjMy8urjTGuzx6dXouACCFeA/AOgFwAnVW7/aCcSFkO4AMhxERDBsiMJ9jLDm/3Cyp//t2hRETfyDRiRIwxxp50kZGRfkQUtmTJEidjx8IeT1J9OwghviSiJVAmx00BSKCcvu2YEKLAwPExI3v9qaY4fCUNJ65lQCGAmZvO4vcZT8HWwtTYoTHGGGM6nT59+oKxY2CNlz4r5PUioqcBQAhRIoQ4IoRYJ4RYLYT4gxPjx5OJCeH/RoTAzlKZDN95UIi5v8YbOSrGGGOscu3bty9q3759kbHjYI2TPmUVfwKYW1eBPAoi8iaib4gogYgKiaiIiK4S0XIialpFv9FEdJSIsokoj4hiiWgyEelVbvK487CzxH+H/K906tezd7Hz7B0jRsQYY0xTVFSU9ZtvvukdHBzc0snJKcTU1DTU1dW1bf/+/Zv+9ddf1trtZ82a5UlEYbNmzfK8deuWdPTo0b5ubm5tzczMQr28vNpMmjTJq6CgoMJksFlZWSZffvmlc58+fQKaNGkSbGlp2d7Kyqp9y5YtW82ePds9Ly+vRhNSKxQK+Pr6BhNRmK741Pr16xdARGGff/65i7qWd/v27U4AMH36dD91bbF2mUVVNcfFxcX05ZdfOnfq1Km5nZ1dOzMzs1APD482PXv2DFy2bJmjZtsrV66Yvfvuu+6dOnVq7u7u3tbMzCzUzs6uXadOnZovX77cUdf5WeOnT1lFJoDUugqktoioPYAoAPYAbgPYrzoUDuBNAC8R0bNCiBNa/ZYCmASgCMBfAEoB9AbwLYDeRDRcCMErMqlEtPXAwQRvbI27DQD4YEc8Qps4wMeR13thjDFj++CDD7yio6NtAgICCkNCQvLNzMwU169ft9i/f7/Dn3/+6bBixYrrEyZMyNLud/v2bdPw8PBWQgiEhYXl5ebmSuLi4mTLli1zv3z5smVUVFSiZvvo6Girf//7376Ojo5yf3//orZt2xZkZmZK/vnnH9nChQu9fv/9d/uTJ08mWFlZVbl6lImJCV599dX7//nPf3y+/fZbl969e1eYL/TGjRumUVFRdtbW1oo333wzo6ioyGTo0KEZMTExslu3bpmHhobm+fn5FavbBwUFFWufQ1taWpqkX79+zc6ePWttZmYmQkND85ydnUtTUlLM4uLiZFeuXLGcOHFi+c01q1atcvriiy88fXx8ips2bVoYHh5edu/ePdO4uDib6Ohom1OnTlmvXbuWFwN4zOiTHMcCaFFXgTyCpVAmxisBTBZClAIAEZkCWA5gAoBlUM6/DNWxSCgT4xQATwshrqr2uwE4CGAIgCkAFtffy2j4PnquNaJvZCI5swC5xXL8a/M5/PJGZ0hMeOUqxpjx+M3ZY5RZCQzh5ucRcYY4z6xZs1I6dep03cfHR665f8OGDXavvPJKwKxZs3yHDx+ebWNj89CcnFu2bHEeOXJk+rp165ItLCwEAJw+fdriqaeeannw4EG7P/74w7pfv37liWtgYGDxzp07r0RERORKJJLy86Snp0uGDh3a9OjRo7affvqp26effppSXcyTJ0/OWLBggdfevXsd7927d9vDw+Oh2BcvXuxSVlZGQ4cOTXdwcFAAUGzbtu1mZGSk361bt8xfeeWV9GnTpmXo8z6NGjXK7+zZs9bt2rXL37FjxzU/P79S9bGCggLavXu3jWb7iIiI7BEjRmSFh4c/VKJx/vx58379+jVft26d69ixYzN79erFiwE8RvQpH/g/AO2IaFRdBaMvIrIA0EX1dK46MQYA1c8fqp62JSLNIc53VY+z1Ymxqk8qAPVsG3O4vOJhMnMpFo1sV54MR9/MxPLD14wcFWOMsWHDhuVoJ8YAMHr06OwBAwZkZWdnS/bs2WOjfdzd3b1k1apV5YkxAISGhhYNGTIkAwD++OMPW832AQEBpc8999xDiTEAODs7ly1dujQZAHbu3OlQk5idnJzKhgwZklFSUkJLly59aOaJ4uJiWr9+vTMATJ8+/X5NzledEydOWEZFRdlbWVkp9uzZk6iZGAOAlZWVGDFiRI7mvh49ehRoJ8YA0KZNm+K33377HgBs2rSpRq+XNR76jBxnAPgOwHrVyOsOAEkACnU1FkKcfvTwqlUG5RRyUgC6hi/V/7HnQxUnEXkDCANQAmBLhQ5CHCaiO1Auld0ZvNrfQ8J8HTC1VyC+Vs15vOjAFXQPdEaIj72RI2OMsSfbvXv3pFu2bLGLj4+3zM7OlsjlcgKAhIQES9WjuXafrl275spksgolEC1atCgCgLt371aYmkihUODAgQOygwcPym7fvm1WVFRkIoSAEMrTJCUlVbhOZWbOnHl//fr1LuvWrXOdN29eqjrp/vHHH+3T09NNO3bsmBsWFmaQG+t2795tBwB9+vR54OnpWeGDRGUKCgpo+/btdtHR0Vbp6emmxcXFBACpqammAJCYmFjj18saB33LKgSUSehQ1VYZoee5a0UIUapagORZAPOISLus4hNV09VC/V8t0F71eEEIoTOxBxADZXLcHpwcVzClZyCOXEnD6eQHGHg+Cm5tJ0Bk3wc1aQJ8+inw0kvGDpEx9gQxVGlCY/bFF184z50716eoqKjSbzxzcnIk2vt8fHxKdLW1tbUtA4Di4uKHznfr1i3p888/H3jmzJlKb6LLy8urcJ3KhIWFFXXp0iXn77//tt26davdyJEjswFgxYoVrgDw1ltvpdX0XNVJSkoyA4CgoKAaJ9t//vmn9ZgxYwLUibAu+rxe1jjoUzZwWrXFafxc2XbGsGFWaRKAqwBeB3CdiHYQ0Q4ANwCMgrJu+G2N9v6qx6Qqzpms1ZZpkEpM8PXI9hhx5Qg+3/ct3B+kgoQAkpKAN94A1q83doiMMfbEOHLkiNXs2bN95XI5ffjhh7fPnj0bn52dfaasrCxOCBE3efLkFAAQQlT4htXERL/qwVdeecXvzJkz1qGhoXk7duy4cufOnXNFRUWnhRBxhYWFtfrGePLkyfcBYPny5S4AEBsbaxEbGytzcXEpHTNmTIWbCOtLbm6uyahRowJTU1NNR4wYkX7kyJFL6enpZ+VyeZwQIm7btm1XAd3vK2vcajy6K4QIr8tAaksIcZ2IugL4EcAAAN4ah2MBHNGsRQYgUz1WVTyfp3qsUJ8FAET0BoA3AKBJkya1CbvRa+Jkhf/8vR5Wcq2bgwsKgPff59FjxhirJxs3bnQQQmD8+PH358+fX2FWqevXrxvka/+cnByTw4cP20kkEuzfvz/R2dn5oRmdLly4UKvrjBo1Knv27NklR44csUtISDD7+uuvXQFg7NixaaamhltwytfXtwQArly5YlGT9vv375dlZGRIW7duXbBp06YKA2pXrlzhcorHVKO/4UyVGMcDCATwPABnAC4AXgDgAGAbEWnOz6z+hFflNDNVEUKsEEKECyHCXVxcanuaRs8q9a7uA8nJuvczxhgzuKysLCmgu0Ti7t270mPHjtlW7KW/zMxMiUKhgJWVVZl2YgwA69atq9VyzhKJBBMmTLivUCjwySefuP/6669OEolETJs2LV1XezMzMwEA6prqmoqIiMgGgD///NP+3r171Q4OpqenSwHAw8NDZ+nJli1beJ7jx1SjTo6JyB7Ar1CO8PYXQvwmhMgQQqQLIXYC6A/ljXgfElEzVbdc1aOs4hnLqY/lVtHmiUeVjJoXeXjVcySMMfbkUtfQbty40Sk7O7v8/+tZWVkmY8aM8cvNzTVITay3t3epra1tWW5urkR7AYytW7farly50q225546dWq6hYWF4ueff3bJz8836dev3wNfX99SXW09PT1LAODSpUs1GgFW69atW2HPnj2z8/PzTQYNGhSQlJT00LB0QUEBbd68ufyDRHBwcBEAnDx50ubMmTPl1yorK8Pbb7/tcfr06aryCNaI6bN8dKg+W10GrSECylHik0KI69oHhRCJAE5BWT7yjGr3TdWjbxXn9dFqy3T59FPA6uFFQAqk5ni/44s4kajzAz9jjDEDmzx5crq7u3vJxYsXrfz9/dv069cvoG/fvgH+/v5tz58/bzV8+HCD/EGWSqWYMWPGPQCYOHGif/v27VsMHjzYv23bti2GDx/e7PXXX6/1QmEuLi5lQ4YMKV98Q12HrEtkZOQDExMTrFmzxq179+7NRowY4Tty5EjfAwcOVHqToNovv/xyIzg4uOD06dOyoKCgNl27dm02ePBg/44dOwa5u7uHzJw5szw36N69e0HPnj2z8/LyJJ07d2719NNPNxs0aFBTPz+/4K+//trjrbfeqnYuZ9Y46TNyHAvlLA412aING2al1EOX2VW0eaB6VH/KVd8s2JqILCvp00GrLdPlpZeAFSsAX18IItyzd8Wc/lOwLagHJvwQwwkyY4zVAxcXl7LY2NhLL774YrqVlZXi0KFDdufPn7fu379/Vmxs7CVvb2+dI7C1MW/evNS1a9deCwkJyU9MTLQ4ePCgvUQiwXfffXfjm2++ufMo5+7bt28OAAQGBhZFRETkVdaua9euhatWrboeHBycf+bMGdmWLVucN2/e7FyTkWQ3N7ey6Ojoy59++mlyq1atCs6fP2/9xx9/ONy+fdssPDw87z//+c9tzfZ79+699t57791p0qRJcXR0tM3x48dtmjVrVvTHH38kRERE5FR2Hda40f9mOKumIZF6KjdtJlCOwjqqjp8DUCaE6KCjrUER0SsA1kE5u0Sg1o136uncrkE5EjxSCLFZtT8OQCiAV4QQP2r16QHgEJSr53kJIR5aTUhbeHi4iI2NNcjraexupudj1IqTSMlRzpJjYWqCNa90QNdAZyNHxhhraIgorqY3ep87d+5mSEgIf9p+zPXt2zfgzz//tP/888+TZ8+ebbAp3BjT5dy5c84hISF+uo7VeORYdQNaBx1bmBDCGcqV6i4AuA2gk0Eir97vAAqgHEFeRETld46qfl4CZWKcBWC/Rr/PVI8LiChQo48rlAudAMDn1SXG7GF+ztbY+EZnuNsqP7wXlSp4BJkxxli1jh49ahUVFWVvb28vnzhxol5LQjNmaAa7IU8IcQrKGSJ64+F5heuMEOI+lPMclwGYDOU8x78R0S4o5zl+C0AxgAlCiGyNflsBLAPgDuA8Ee0iou1QzpfcCsqb/L6tj9fwuOEEmTHGWE2NHDnSd9CgQU0HDBgQpFAoMHv27Lu2trY8MMWMyqCzVahuiosBMM6Q563mmj8A6AjgJyiXhO4HoC+Us1SsBhAqhPhVR79JAF6CctGSHlCuspcIYAqASCFEhWlqWM1wgswYY6wmNm/e7Pz777872Nvbyz/88MPbc+bM4XIKZnR1scRzOpTJar0RQpwG8HIt+m0AsMHwETF1gqyuQVYnyFyDzBhjTE0I8cQv/c0aHoOOHBORFZT1xpXeZcqeHDyCzBhjjLHGRp95jh2r2JoQUT8AewB4AvirziJmjQonyIwxxhhrTPQZOU6rYrsB5cwRPQDcB/CuYcNkjRknyIwxxhhrLPRJjrMAZFaypQCIA7AAQIgQ4qZhw2SNHSfIjDHGGGsM9Jnn2FkI4VLJ5iWE6CiEeFc1vRpjFVSaIF/jBJkxxhhjDYNBb8hjrDo6E+R1nCAzxhhjrGHQ54a83y9f/L0AACAASURBVIhoWg3aTSGi3x4tLPY483O2xi+cIDPGGGOsAdJn5HgQgHY1aBcCIKJ24bAnhT8nyIwxxhhrgOqirMIMAC/9yKrFCTJjjDHGGpq6SI7bA8iog/OyxxAnyIwxxhhrSKpMjlV1xr9p1BD31tynte0lomsAWgM4VueRs8cGJ8iMMfbofv75Z/vQ0NAWMpmsPRGFEVHYiRMnLBMSEsyIKMzLy6uNoa7l5eXVhojCEhISzAx1zroSGRnpR0RhS5YscdLcP2vWLE8iCps1a5andp+CggKaNGmSV5MmTYLNzMxCiSisRYsWrdTHExMTTZ977jl/V1fXtlKpNIyIwiZMmOBTH6+H1T1pNccHafwsAPiotqokAJj9KEGxJ486QX5xxUmk5BSVJ8hrxnVA1wBnY4fHGGMN2vHjxy3HjRvXFAA6d+6c6+bmVgoALi4uZQpF/VU6RkZG+m3fvt1p8eLFN6dNm9Zov0WeMWOG18qVK92cnJzkffr0eWBpaanw8fEpAQCFQoEhQ4YExsfHWwUEBBR17tw519TUVHTs2DHf2HEzw6guOR6seiQAvwH4E8DiStqWALgjhLhooNjYE4YTZMYYq52tW7c6lJWV0ZQpU1K++eabO5rHiouL6fTp0xfMzMyEoa73xx9/XCkpKSE/P79SQ52zvv373/++P3bs2Ex3d3e59rHdu3c7AMDBgwcvt2nTpljz2JUrV8zi4+OtPDw8Si5dunTB1NS0vkJm9aTK5FgIsUf9MxHFATiquY8xQ+MEmTHG9Hf79m0zAGjWrFmR9jFzc3PRvn37CvsfRevWrYurb9WweXh4yD08PCokxgCQkpJiBgDaiTEA3LhxwwwAvL29Szgxfjzps0JeByHEx3UZDGMA1yAzxhqA5csd4enZBiYmYfD0bIPlyx2NHZIu6rrZrVu3OgHA9OnT/dT1xpGRkX4AUFXNsbotAKxcudKhXbt2LaysrNpbW1u379KlS/P9+/fLdF1Xu+ZYfY3t27dXiENXvW9KSopk2rRpns2bN29lZWXV3tLSsn2rVq1azps3z7W4uJj0fR9ycnJMpk6d6uXj4xNsZmYW6u7u3nbMmDFNUlJSJNW9d5o1x+rXJYR46P1RvwYiChs4cGAQAMTExMg0j2ueu7i4mBYuXOgSFhYWZGtr287c3DzU19c3+LXXXvO+e/duhYFJ9bkjIyP9UlJSJOPGjfPx8vJqY2pqGtqnT58AzbaJiYmm48eP9/Hz8wu2sLAIlclk7UNDQ1ssWbLESVcJTceOHYOIKGz37t02R48eterVq1egvb19OwsLi9CgoKBWixYtqnTkSaFQYNWqVQ5PP/10M0dHxxBTU9NQV1fXtl26dGn+3//+10VXn23bttn26tUr0MnJKcTU1DTUxcWl7eDBg/2jo6MtK7tOQ1NdWQVjRsEjyIwxo1m+3BEzZ/qiqEg5gHTvnhlmzvQFALz1VqYxQ9PWvn37gqFDh2bExMTIbt26ZR4aGprn5+dXDADdunXLq+l5ZsyY4fnNN994hIaG5vXs2TP70qVLlidPnrQZPHhw87179yb06dOnynpaW1tbRWVxAEBQUFD5z9HR0ZaDBg1qlpaWZurm5lbaqVOnXIVCgXPnzsk++ugjn/3799tHRUVdtbCwqFEZSE5Ojkm3bt2C4uPjrWQyWdnTTz+dLZFIsGvXLscjR47YNmvWrLCm70NERERWRkaGVJ3kDx06tLxuOigoqHjo0KEZ9+/fNz127Jitk5OTvEePHtna58jMzDTp27dvs9OnT8tkMllZcHBwga2tbVl8fLzV6tWr3fbu3etw8ODBhKCgoBIdfaXh4eGt8vLyJOHh4blt27YVDg4O5aPbu3btshk9enRAXl6epEmTJsVPPfVUdn5+vsm5c+dk06dP9zt48KDNjh07bup6bXv37rVduXKlm7+/f9FTTz2VfefOHfMzZ85Yz5o1y/fBgweSefPmpWq2LyoqooiIiKZRUVH2EokEISEheZ6eniXp6emmV65csXz//fdt3nvvvTTNPuPHj/dZt26dq0QiEW3atCnw8PAouXnzpvnu3bsdDxw44PDDDz9cGzlyZIX3rKHROzkmIkcoF/loDsAWynpkbUIIMf0RY2NPOE6QGWNGMX++V3lirFZUZIL5870aWnI8duzYB2PHjn0QGRnpd+vWLfNXXnklvTY3wq1bt8710KFDl5566qkCACgrK8OYMWN8N27c6Dx37lzPPn36XK2qv4eHh3zbtm03q4sjLy+Phg4dGpiWlmY6Z86cO/Pnz09RlyakpqZKhgwZ0vTvv/+2fe+99zy++uqruzWJ/V//+pdnfHy8VbNmzQoPHjx4xcvLSw4A6enpkn79+jWLioqyr+n7sGLFitsAQEROALBt27abmsefffbZvN27d9scO3bMtmnTpkXaxwHg5Zdf9jt9+rSsf//+WT/++GOSi4tLGQDI5XJMnTrVa/ny5e5jx471j46OTtDue+jQIbtu3brl7Nq165qDg8NDw8BJSUmmY8aMCSgsLJQsWbLk5uTJkzNMTJT/TBMTE00HDx7c7Ndff3VasmRJrq73ftmyZe6LFi26OWPGjPJj3333nePkyZP9v/rqK4+33347zcbGpvyaEydO9I6KirL39fUt3rFjR6JmaY5cLsemTZvsNM+/cOFCl3Xr1rkGBgYWbd68+Zpm+59++sl+/PjxTV9//XX/Xr16nVe/Jw2VXvMcE9HrAJIBrAPwHoApGttk1aZ+ztgjq6zE4u9rjfYmaMZYQ6eqN63x/sfAO++8c0edGAOARCLBl19+eQcA4uLibGpT6qDL0qVLne/cuWM2cODArM8++yxFs2bXzc2tbMOGDTelUqlYu3atS01m2cjLy6MNGza4AMBXX311S50YA4Czs3PZsmXLkogMEnqNxMXFWezZs8fB09OzZMuWLTc0k0CpVIpvv/32TvPmzQtjYmJkusoMpFKpWL16dZJ2YgwAn3/+uWtOTo7k9ddfT5k6dWp5YgwAgYGBpd9///1NAPj+++9ddcX27LPPZmkmxgAwadKkzKZNmxbl5eVJjh07ZqXef+fOHenPP//sYmJigq1btyZq16xLpVK89NJL5SPAcrkcX375pQcAbNy48Zp2+7Fjxz4YPXp0em5urmTFihUPldg0RDVOjomoJ4DlqqffAIhR/fwvAMsA3FM9XwpgmqECZEydILvZmgNQJsjj10VzgswYqxvu7hW+7q5y/2MgMjKywlfdXl5ecltb27KSkhJKTU2ttHZXH/v377cDgGHDhmXpOu7n51fq6+tb/ODBA2l8fLx5dec7fvy4dUFBgYmrq2vpoEGDcrWPd+rUqbB58+Y1Lqt4VL/99psdAPTu3TtbJpNVKAuRSCTo2LFjHgAcOXLEWvt4q1atCnSVWwDAX3/9ZQcAL774os73rnv37gVWVlaKy5cvWxUUFFT4RDBw4ECd5QwBAQFFAHDr1q3yTyp79uyxkcvl1K5du7zw8PBqb+b8+++/rdLS0kwDAwOLwsLCdLZ/5plncgHg5MmTFV53Q6NPWcVM1WN/IcQxIloLoIMQYhEAENG/AawA8BKUq+QxZjD+ztbY+EYXjFrxN1JzissT5LXjOqJLQIP/EMoYa0zmzr3zUM0xAFhYKDB37p0qejVqgYGBOhMymUxWlpOTIyksLDTIirrJycnmADBhwoSmEyZMqLJtSkqKtG3btlXOipGUlGQKAN7e3pW28/b2Lk5ISKiXm8GuX79uDgA//fSTy08//aTzhjW1tLS0CjmYt7d3pR/Abt26ZQ4APXr0aFldHKmpqVJ/f/+Hptnz8/PTeW4bG5syACjS+PeelJRkDgCBgYE1muXk6tWr5gCQmJhooX1zoraMjIwGf7+bPgF2BHBGCKFz9TshRAERvQrgJoB5AMY9cnSMaeAEmTFWL9R1xfPneyElxQzu7iWYO/dOQ6s3NiSJxCADw9UqK1NWGTzzzDPZjo6OOqdRU2vodam6qF9f69atC4KCgqocsQ4ODq6QeFpYWFRaS6JQKAhQ3jRobm5eZc2JrpsZNcswDE0uV/4qXV1dS7t3755TVdugoCCDTitYF/RJjh0AHNJ4XgoARGQlhCgAACFEMREdA9DbYBEypoETZMZYvXjrrczHORk2FtXsBRZvvvlm2qhRox551oImTZqUAsCdO3cqLcG4fft2teUZhqJeRa9bt26533///W1Dntvd3b0kOTnZ/KOPPrpbk1KHR+Hr61sMKEeCa9JePSrt4uJSqusmxcZGn48RGQA05ztU/9Hw1WonBcBZCqsz6gSZa5AZY6xhUa/CJ5fLdd4F169fvxwA2LJli4MhrtetW7cCS0tLRWpqqunvv/9eYU7mmJgYiytXrtTb/LqDBw/OAYB9+/bZl5YadvHAnj17ZgPAhg0b6nzO7YiIiFypVCrOnj0rO336dLUJco8ePQrs7e3lly9ftqpJrXhDp09ynASgicbzf6Ccxi1SvYOIHAD0AHDLINExVglOkBljrOHx9PQsAYBLly7pTKhmzpyZ5u7uXrJ9+3anmTNneubm5lbIQ2JiYiwWL15co0E2GxsbxahRo9IBYNasWU00F9jIyMiQTJw40Ve9oEd96N69e0GfPn0eJCcnm0dERARcu3atwhJ6SUlJpvPnz3fVN3n+4IMPUmQyWdk333zj/tlnn7no6v/XX39Zr1mz5pE/eHh5eclfeumlNIVCgeHDhwf8888/DyW8crkcGzZsKJ/KzdzcXMyaNeteWVkZhgwZEnjw4EEr7XPm5OSYfP/99441SbaNTZ+yiigAs4nISwhxB8AeANkAPiIifwC3AYwCYAdgrcEjZUwLl1gwxljDEhkZ+eDrr7/2XLNmjdulS5csPT09S4gIr732Wnrfvn3z7ezsFDt37kx84YUXAr/++muPtWvXugYFBRW4urqWpqWlmd66dcv87t27Zm3bts2fPn16jUY7Fi1adOfUqVOyixcvWgUFBQV37tw5VyKRiJMnT9ra2NjIe/Xq9UCfuY4f1aZNm27079+/2YEDB+xbt25tFxQUVODt7V2Sm5sruXfvntn169ctFAoF3n777TRTU9MaZ+6BgYGlGzZsuDZ27NiA9957r8miRYs8AgMDCx0dHctSUlJMk5OTzdPS0kwjIiKyJkyYoHNGC30sW7bs9s2bN80PHz5sFxoa2rpdu3b5Hh4eJRkZGaYJCQmWmZmZ0tGjR8ep23/44Yf3k5KSzFavXu3Wq1evls2bNy/09fUtVigUUL/uoqIiky1btlwNDQ1t0HXH+owcbwKwBcrFPyCEyAYwEUAZgPEAPgTQDMBlKG/IY6zO8QgyY4w1HF27di1ctWrV9eDg4PwzZ87ItmzZ4rx582ZnzZHkjh07Fp4/f/7inDlz7vj6+hZdvHjRav/+/Q43btywcHFxKZ0+ffq9FStWJNX0mnZ2dooTJ04kTJo0KcXW1rbs8OHDdmfOnJENHDgwKzo6+rK9vX293tjn6OioOHHiRMK33357Izw8PDc5Odl8//799vHx8VYSiUSMHj06bdu2bVetrKz0HtIePHhw7j///BM/ZcqUFEdHR/m5c+dkBw4csL93756Zn59f8bvvvntnwYIFBplVxdLSUvz111+JS5cuvdGhQ4e8q1evWu7bt8/h2rVrFkFBQQWfffZZsnafVatW3d67d2/CoEGDMnNyciSHDh2yO3XqlE1hYaFJ7969s5ctW3ajX79+NV650VjoUb9uIKIAAC8AcIQyMd4shKhy6pXHSXh4uIiNjTV2GE+8G+n55SPIAGBhasIjyIw1YEQUJ4QIr0nbc+fO3QwJCUmv65gYY0+Oc+fOOYeEhPjpOqbPIiAmpGOZGSHENSHE/wkh3hdC/PQkJcas4eARZMYYY4wZgj5lFXL8b1U8xhocTpAZY4wx9qj0SY7zoCybYKzB0pUgT1gXwwkyY4wxxmpEn+Q4AYBHXQXCmKFoJ8iFpWWcIDPGGGOsRvRJjtcA6E5EwXUVDGOG4u9sjV9e78wJMmOMMcb0UuPkWAixDMB6AH8R0WQi8qy7sBh7dE1dZJwgM8YYY0wv+sxWkQNgGABnAEsA3CKiIiLK0bE98nrpjBkCJ8iMMcYY04c+ZRUy1UYam5nGfs3NxrBhMlZ7nCAzxhhjrKb0SY5t9NwYazA4QWaMMcZYTehTc5yvz1aXQTNWG5wgM8YYY6w6+owcM9bocYLMGGOMsaronRwTkTcR/YeI9hFRHBHN1zgWSkSjiYjLKliDxQkyY4wxxiqjV3JMRKOgXCVvLoB+ANoB8NZo4gvgJwAvGCpAxuqCrgR529sLUOTlA5iYAH5+wPr1xg2SMcYYY/VOn6ncOgD4EcpZKuYD6K36WdNeKJeZft5QATJWVzQT5OcuHMT83Utgcfc2IASQlAS88QYnyIwxxtgTRp+R4zmq9hFCiHlCiIPaDYQQxVCOLLc2UHyM1Sl1gvzusZ9gJS9++GBBAfD++8YJjDHG9EBEYUQUZuw4GHsc6JMcdwcQI4Q4VE272wA8ah0RY/WsqYsM7tlpOo+J5OR6joYxxlhVlixZ4kREYZGRkX7GjuVJ0rFjxyAiCtu9e/djf1+ZPsmxPYCkGp7TrHbhMGYc1KSJzv0pdi44djW9nqNhjDHGmLHokxynA/CrQbvmAO7VKhrGjOXTTwErq4d2FUjN8Vn3sRiz+hQ++u0CCkvKjBQcY4wxxuqLPsnxCQBhRNS2sgZE1ANASwCHHzUwxurVSy8BK1YAvr4AEQo9vPHJ8zPwW+ueAIB1J24i4pujOHfrgZEDZYzVOWfnEBCFVdicnUOMHVpNfPnll84tW7ZsZWlp2d7e3r5dv379AmJiYix0tY2KirJ+8803vYODg1s6OTmFmJqahrq6urbt379/07/++staVx+5XI6FCxe6tG/fvoWNjU07U1PTUCcnp5BWrVq1fP31173v3r0r1e6Tk5Nj8sEHH7gFBwe3lMlk7S0sLEIDAwNbz5o1yzM7O7vGuYiXl1eb6dOn+wHA9u3bndS11rrKLIqLi2nhwoUuYWFhQba2tu3Mzc1DfX19g1977TWdMWqWa6SlpUnGjRvn4+Hh0cbCwiK0adOmrRcuXOiibhsbG2sxcODApk5OTiEWFhahbdq0ablt2zZbXTFr1oPr87upzfs2a9YsTyIKmzVrlueVK1fMhg0b5ufm5tZWKpWGTZgwwUf9vixdutRx8ODB/n5+fsHW1tbtLS0t2wcEBLSeOHGiV2pqqkTznLt377YhorCYmBgZAAwePLi55vuuLrOortxFfZ6OHTsGVbY/NzfXZNq0aZ7+/v6tLSwsQlu0aNFKs21KSopk2rRpns2bN29lZWXV3tLSsn2rVq1azps3z7W4uFh7kohaq/CPowpLAEQC2E5ELwshTmgeJKJQAGsAKAAsNVSAjNWbl15SbgAsAczMLcb97f/gz0v3AQDX0/IxdNkJTOkZiCm9AmEq4TV0GHssZWTo/n9jZfsbkFdffdVn3bp1rmFhYXl9+vR5cP78eesDBw7YHz161PbXX3+9+uyzz+Zptv/ggw+8oqOjbQICAgpDQkLyzczMFNevX7fYv3+/w59//umwYsWK6xMmTMjS7DNy5Ei/7du3O1lYWChCQ0PzHB0d5ZmZmdLk5GTzVatWuY0aNSrL09NTrm5/7do102effbb5tWvXLBwcHOTt2rXLMzc3V5w/f9560aJFHnv27LE/duxYgouLS7Vfz0VERGTFxcVZnz59Wubj41PcoUOH8tfTrVu38p8zMzNN+vbt2+z06dMymUxWFhwcXGBra1sWHx9vtXr1are9e/c6HDx4MCEoKKhE+xrZ2dmSjh07tsjPz5eEh4fnZWVlSWNiYmSzZ89ukp2dLenZs2fu888/39zT07OkS5cuuTdu3DCPj4+3GjlyZLNdu3YlDBgwIE/7nLX53TzK+5aYmGjeoUOHVubm5oqwsLA8uVxO9vb2ZQBw+/Zt6ZQpU/xtbW3LmjZtWtS6deuC3Nxcyfnz562WL1/uvnv3bofo6OjLHh4ecgDw8vIqHTp0aMbhw4ftMjIypN27d89xdXUtVV/Ly8urVPv6tVFcXExdu3YNun79ukWHDh1yW7VqVVhSUlKe8EZHR1sOGjSoWVpamqmbm1tpp06dchUKBc6dOyf76KOPfPbv328fFRV11cLCQjxyMEKIGm9Qzm+sAFAG4K7q8T6Aa6qfFQDe1+ecjX0LCwsT7PGlUCjExugk0erD34Xv7N3l23PfHBVXU3ONHR5jjRaAWFHDv7Nnz569KYSIrbdNOaGj7q0+49BjAyAACAsLi7K9e/deVu8vKyuLnTRp0j0Awt3dvTg/Pz9Os9+WLVuuJCcnn9U+3/r1669KpVKFnZ2dPCcn57R6f0JCwj/qc+nqd/z48Qu3b98+q3n9du3a5QEQL7/8cqrmuXJzc+Oef/75DABi6NCh6TV9rYsXL75RXZ+IiIhMAKJ///6Z9+/fP6PeX1paGvvWW2/dAyA6dOiQq+u86n6a79WmTZuuABBWVlZlnp6exXPnzr2l2feNN95IASA6d+6cY4jfTW3ft5kzZ95VXy8yMjK9sLAwTjuezMzM0+vXr79aVFT00LHc3Ny4YcOGpQMQo0ePvq/dr0OHDrkAxK5duxJq83vZtWtXgq73Xb0fgGjRokWBrn9Xubm5cV5eXsUAxJw5c26XlJSUH0tJSTnTpUuXbABi5syZd2v670j1d0Xn3xy9hr6EEPMBDIdyujZ3KOc5dgbgD+AGgDFCiE/1OSdjDRkRYWSHJvh9+tPo4OdQvv/c7WxELDmKdcdvQKF49A+pjDFmCGPHjk3THLk0MTHB4sWL73h7exenpKSY/fDDDw6a7YcNG5bj4+Mj1z7P6NGjswcMGJCVnZ0t2bNnT/nsBPfu3ZMCQOvWrQt09evatWuhl5dX+f6tW7fanj171jokJCR/zZo1t2xsbBTqYzKZTPz4449Jjo6O8p07dzqmpaVJtM9XG3FxcRZ79uxx8PT0LNmyZcsNzZFVqVSKb7/99k7z5s0LY2JiZNHR0Zba/a2trRVr1qxJtrKyKv/jPmLEiJygoKDCgoICE1dX19J58+alavaZP3/+PdW1ZZV9va/P7+ZR3zc7O7uylStXJusaRXVwcFCMHj0629zc/KFjMplMrFmzJlkikYjff//dQbtffViyZEmSrn9XS5cudb5z547ZwIEDsz777LMUU1PT8mNubm5lGzZsuCmVSsXatWtdFAqFdne96f29sBBimxCiNZQ35/UC0BdACyFEoBBiwyNHxFgD1MTJChvf6II5A1rATFVOUSxX4KNdF/Hymmjcyy40coSMMQaMGzcuQ3ufVCrFkCFDMgHg8OHDFabhunfvnnTJkiVOb7zxhvfIkSN9IyMj/SIjI/0SEhIsASAhIcFc3bZt27ZF1tbWikOHDtnNmTPH/cqVK1XOTrVnzx47AHjuueeyJJKKua+tra2iTZs2+WVlZXT06FGdNc76+u233+wAoHfv3tkymaxCciiRSNCxY8c8ADhy5EiFawYHB+erSwo0+fn5FanPq33Mzc2tzN7eXl5aWkraNbtq+vxuHvV96969e46Dg0OVWeLx48ct586d6/byyy83GTZsmF9kZKTfuHHjmpiamoqsrCypoT6s1JSTk5O8b9+++bqO7d+/3w4Ahg0blqXruJ+fX6mvr2/xgwcPpPHx8ea62uij1vVTQohkADwJLHtiSEwIb/UIQI/mLpi56Swup+QCAI4lpqPfoiP4+PlgPN/OE0QGuyeAMcb0oquGFgD8/PxKAODu3bsPJbNffPGF89y5c32KiooqHSzLyckpT5IcHBwU33zzzc2pU6f6LViwwGvBggVerq6upaGhoXkDBgzIfu211zI1R1yTkpLMAeDjjz/2/vjjj72rij01NdUgNd3Xr183B4CffvrJ5aeffnKpqm1aWlqFa3p4eOh8D62trRUA4O3trfO4lZWV4sGDBygsLNT5Xurzu3nU983Hx6dYV1sAyM7ONhk6dKh/VFSUfVXnzcrKktSkDtxQPD09K405OTnZHAAmTJjQdMKECVWeJyUlRdq2bdtKz1UTtf6HSES2ALygrBO5J4So8EmKscdRSw9b7JzSDYsOXMX3R65BCCC3SI4Zm87iwMVUfPJCMByseapvxhotJye5zpvvnJwqjCY2NkRUnrgeOXLEavbs2b4SiUR8+OGHtyMjIx/4+/uXymQyhYmJCaZMmeK1dOlSdyHEQ5/4x48fnzV48OCcX375xf7o0aOymJgY2b59+xz27dvnsGDBAs/Dhw9fDgwMLAWAsrIyAoAOHTrkVZWwAUDTpk11Jo/6KitT5nOtW7cuCAoKqvJrveDg4CLtfSYmVX+pXt3x2tL83Tzq+2ZpaVlpvd+0adO8oqKi7AMCAormz59/u1u3bgXu7u5ydZmFq6tr27S0NFMhDFsyWF25Q1U30ql/p88880y2o6Njlf8dGiKh1zs5JqIxAGYCCIGy5hgABBGdA7BYCPHjowZVG0RkCWAqlDXRzaBciCQVQCyAr4UQx7XamwCYCGA8gBZQ3lD4D4DvhBC/1GPorBEyl0owZ0AL9G7pilmbz+JWpvLv757z9xBzMxMLhrVFzyBXI0fJGKuV9PRzxg6htq5cuWLWpUuXCgnhzZs3zQDA3d29fGaBjRs3OgghMH78+Pvz589P1e6jHoHVxdnZuWzq1KkZU6dOzQCACxcumL/66qu+p06dspk5c6b3rl27bgCAl5dXCQAMGTIk891339W9FKmB+fj4lABAt27dcr///vvb9XHNmtDnd1OX79uePXscAeCXX3651qFDh4c+HOTk5Jikp6eb6u5ZNTMzMwEA+fn5Oj893Lhxo9ajRh4eHiU3b960ePPNN9NGjRpV54OxNf74Q0o/AfgBQHtVZGn2cwAAIABJREFU3xzVZqLat5aI1lM9f69MRP5QJrYLADSBcp7l3QDSADwPoKdWewmAHQC+hTKR/gPAMQAdAGwgoiX1Fjxr1Dr4OeL36U/jxY4+5fvu5xZj/NoYvL/jPPKLG/1AE2OsEfnhhx+ctPfJ5XLs3LnTEQB69OiRq96flZUlBf6XTGq6e/eu9NixYzrn7dWldevWxe++++49ALh06VL5ikoDBgzIBoAdO3Y46vM6qqJOwuRyuc5cY/DgwTkAsG/fPvvSUoPMMmYQ+vxu6uJ9U8vOzpYAQNOmTSu8OStXrnSsbMTY1NRUoY5ZF/W/o2vXrlW4yREA9u3bZ1e7iIF+/frlAMCWLVvq5UZBfb4bmATgJShXypsGwEEI4SCEcIByaempUE7rNgrAZEMHWhkisgZwAEAggI8BeAshnhdCjBBCdATgAWCzVrcZAJ4DcBFAcyHEUCFEBIA2UI42TyWi5+vrNbDGTWYuxWdD22L1K+Fwlv1voGX9qWQMXHIUcUk67x9gjDGD+/HHH132798vUz9XKBSYNWuWZ3Jysrmrq2vpyy+/XP4HKSgoqAgANm7c6KS5oERWVpbJmDFj/HJzcyvckHX8+HHLlStXOuTl5VVITHfu3GkPAKoptwAAY8aMedC6deuCmJgY2ejRo5voulnt4sWLZp999lmVtcGa1ElYYmKizsUzunfvXtCnT58HycnJ5hEREQHXrl2rMBKalJRkOn/+fNf6TJ71+d3Uxfum5u/vXwQAX3zxxUN9jxw5YvXJJ594VdbPw8OjFAAuXLigM/l9+umnC6ytrRWJiYkW33///UNJ/eeff+6yb9++Wie2M2fOTHN3dy/Zvn2708yZMz1zc3Mr5K8xMTEWixcvrvABpDb0Kat4A0AxgGeEEJc0DwghcgAsJaIoAGdUbb81RIA18AGAAAA/CiHmah8UQmQAKL9DVDVq/I7q6UQhRKpG26tENBvAOgDvA9hZh3Gzx0zvlm7YP8Me7++Ix74LKQCApIwCDF9+AhOfCcD03s1hJuWFQxhjdefFF19MHzhwYFB4eHiuq6traXx8vNXNmzctLCwsFGvWrLn+/+3deXxU1f3/8dcnK4SdJKgIQkBRUcQFVFBxb/WrrVqrrdpa97r8qt9al1r7ta1WRavW1rqBtdZWrah1q637AsqiYKUqYtkhbCbsYc3y+f1x7sAwmSyTdRLez8djHpe559w7585JyGfOnPs58dkbrrjiitKHH36414wZM/KKioqGDBs2rMzd+eijj7pkZ2dXnXHGGaXPPPNMQfz5586dm3vJJZcMuPLKK6sGDx68oXfv3lvKy8vt888/zysuLs7t1KlT1S233LIkVj8zM5MXX3xx9oknnrjHU089VfjSSy/l77nnnht69+69ZcWKFVlLlizJXbBgQW5+fn5FfacPHHPMMesLCgrKZ8yYkbfvvvvuPWjQoI3Z2dk+cuTIsquuumoFwNNPPz3vhBNO2OONN97ovs8++3Tbc889N/Tp02fLunXrMpcuXZozd+7cDlVVVVxzzTUl2dnZLZKPM5W+aY73LeaGG25YesEFFwy44447dn3hhRd67rHHHhuXLVuW8/HHH3c+6aSTVk6bNq1z4o2bAKeddtqq5557Lv/mm2/u89Zbb3UtLCwsB/jZz362bOjQoZu7dOlSdfXVVy+55ZZb+lx22WVFY8aMKSwsLCyfOXNm3uLFi3MuvfTSZQ899NDODXnvunXrVvXiiy/OPvXUU3e/9957d/nTn/7Ua88999zQq1ev8pKSkuxFixblLlmyJGe//fZbH/sZaIxU/lIPAt5NDIzjRWXvEEZxm52Z5QAXR09H1/OwEUAvoNjdxycpfwYoB4abWY2foESSye+cy4PfO5C7zxhKl9zw2bPK4f535nDaAx/w3+Xr6jiDiEjDjR07dtGtt966cPXq1Vlvvvlm95UrV2Yfd9xxq997772ZJ5100nYrsBUWFlZOnTr1i7POOqs0Ly+v6t133+326aefdjrhhBNWTZ069Ys+ffpUG1Y98sgjy2644YbFw4cPL1u2bFnOm2++2f2DDz7o2rFjx6qLL754+ccff/z5qFGjNsQfM3DgwPJPPvnki9tuu23h3nvvvWH27NkdX3311R6zZs3q2Llz58pLLrlk+VNPPTW7vtfYsWNHf+mll2YdddRRa4qLi3NffPHF/HHjxhWMHz9+ayq0nj17Vk2cOPHLP/zhD/OGDRu2buHChbmvvfZa988++ywvMzPTzz777JLnnntuVnxmjeaWSt9A079vMeeff/6ql19++b+HHHLIumXLluW89dZb3cvKyjJvvvnmRc8///y8mo4755xz1owePXphUVHRpkmTJnUZN25cwbhx4woWLVq0NZC++eabl//2t7+dP2jQoI2fffZZp4kTJ3bt16/fptdff33mSSedtDbVtsY7+OCDN3766aczfvrTny7u16/fphkzZuS99tprPebNm9ehsLCw/Kqrrlo6ZsyYBY15jRir792IZrYceMvdz66j3lPAMe6+UxO0r642jQAmAovcfTczGwmcDOQDy4BX3X1SwjE/IiyF/by7f6uG8/4b2B842d1fqa0Nw4YN86lTpzb+YqTdKV61gWuf+Q+T5m77EJuTlcG1X9uTCw8vIiNDKd9kx2Vm09x9WH3qTp8+ff7QoUNLm7tNIs3BzA4CcPdprd0W2Wb69OkFQ4cO7Z+sLJWR43eBkWZW41SMqGwE4Ya4ljAk2s4ys8eAD4AbCNM6bgImmtmzUSaLmKJoW9uni1j+5qJa6ojUqk+PPJ646BD+7+TBW6dTbKmo4tZ/fsFZYydTvGpDHWcQERGRlpZKcPxzwo13Y82s2go7ZtYZGBPV+XnTNK9OsQnfo4BzgbsIUzp6ELJULAZOB+6POyY2GT7pKiyR2Ncb1a4TwMwuMbOpZja1pKRFMtNIG5WRYVx4eBGv/Ohw9t11243fU+at5IR7J/DM1EU0dS5JERERabhUbsj7BvAscAFwipn9A4jNTekflXcDHgVOTszm5u73NLaxScSC+yzgEXe/Nq7sJTNbAnwI/MDMfu3uc4nLzdzQF3X3MYQPAgwbNkyRjdRpj5268PfLDuO+t2dx/zuzqXIo21zBtc/+hzdmLOf2bw0hv3OjV7wUERGRRkolOL6LbQFld+B7cc/jI+HEdf0sqtccwXH83U1jEwvdfaqZTQOGAUcBc+OO6ZxYP06sTHdPSZPJycrgJ1/bk6P36sVPxk1nXmn48uL1Gcv5eOEqbv/Wfhw/uNmn6ouISAvSXOO2J5Xg+B4aMdraTObH/bumOyznEYLjWPqQ2DH9ajlvbEWH+bXUEWmQA3frwStXHs7t/5zJXyaHqe+lZVu4+PGpnDmsD/938mC6dGjQAkUiIiLSSPUOjt39muZsSAN9HPfvfMKKeIliORpj84hjxwxPdkIzywP2jZ7+u7ENFEkmLyeLW07dl+MG78S1z0znq3UhZ/64qcVMnLOCe87cn4OLmnxhJBEREalDm16RwN0XA1Oip8cmlptZD+DA6Gks39okwkp+fcxsVJLTngFkAx9F5xdpNkcOKuT1H4/i5P122bqveNVGvjNmErf/8ws2V1S2YutE0oduXBWRplLX/ydtOjiO3BptbzKz/WM7zawD8CDhJsFphKAYd68EfhNVe9DMesUdswfbFhOJnVekWXXPy+EPZx/I7767P107hC9z3OHh8XP55n0fMGNJo/Kmi7R5ZlZWUVFRbflcEZGGqKioyDKzaguvxKQy5xgAMxsKHA30BpKuaw64u1+V6rkbwt1fNrO7gGuAKWY2hbBc9MFRGxcDZ/n2HxN+S0j/9g1CjuS3CKPFxxGu6T5319LR0qJO2X9XDinK59pnpzNhVljv4Mvl6zjl/vf58fGD+OGogWRq4RDZAbn7hLVr1x6Xn5+/prXbIiJt39q1azu7+xs1ldc7OI5GYv8KnBbbVUt1B1okOAZw92vNbCLwI+AAII+wkMc9wGh3L0moX2lmpwKXA+cDXwcqCSPMD7j7ky3VdpF4O3frwOMXHMxfJi/gtn9+wabyKsornTtf/ZK3v/iKu88cSr/8Tq3dTJEWVVFR8UxJScnx3bp1y8jKyqpq7faISNtVUVGRUVJSQkVFxTM11Ull+eh7gP8l3Ng2DpjFtpvcqnH3+2sqa0+0fLQ0lzklZVw9bjrTF63eui8vJ5OfnzSYsw7uS2IucZG2JJXlo6dNm2YZGRnXd+zY8aLCwkK6du1alpWVVaHfARGpD3enoqIia+3atZ1LSkrYuHHjI1VVVXccdNBBSYPgVILjxYQR2YOixTQEBcfSvCoqq7j/nTn8/u1ZVFZt+109Zq9ejD59CL261DSzSSS9pRIcQwiQgaOzsrLOMLMj3L22XPUiItsxszJ3nxCNGL9TU2AMqQXHG4A33P2UJmpnu6DgWFrCf4pX8+OnP2FOybZVz3vkZXPbaUM4ccgutRwpkp5SDY5FRFpKKtkq5pF+i4CI7BD269OdV648gvMP679136oN5Vz2xMdc/fQnrN1U3nqNExERaUdSCY4fB44ys/zmaoyI1KxDdia/+MY+PHHRIezSbdt0ir//ezEn/HY8E2eXtmLrRERE2odUguO7gInAm2Z2SDO1R0TqcNjuBbz6v6M47YBdt+5bsmYTZz8yhV+9/DmbyrVwiIiISEPVOziOFs/4LpAJTDSz1Wb2mZn9J8ljerO1WETo1jGb335nfx4450B65GVv3f+nD+Zz8n3vs/C+R6B/f8jICNsnnmi1toqIiLQlqeQ53hV4Dygi5DjuCgyuobrmJou0gP8ZsgvD+vXg+uf+wztfhnTeg995mYJX/wAVm0OlBQvgkkvCv885p5VaKiIi0jakkq3iScLI8YfAfcBsas9z/HlTNDDdKVuFpAN356kPF/HrV2bw+u/Opc/akuqV+vWD+fNbvG0iyShbhYikq1SC46+AjcBe7r6xWVvVhig4lnSyYMV6+hZ0ISPJlzduxoaNW+iUm/Kq8SJNTsGxiKSrVG7I6wBMVmAskr765XfCduubtGxxlwJGjn6bO1+dyVdrN7Vwy0RERNqGVILjTwGlcRNJc3bbbZCXt92+DVm53DnqXNZsLOeBd+dw2B1vc80z0/ly2bpWaqWIiEh6SiU4/i1wpJkd2FyNEZEmcM45MGZMmGNsRtVuuzH1xtFMH3XS1irllc6z04r5+r3jOffRD5kwq4T6TrESERFpz1KZc9wTuA74ITAaeA0oBqqS1Xf3lU3UxrSmOcfSVlRWOW/MWMbYCfOYtmBVtfK9du7CxUcM4BtDe5OTlcrnZpHUac6xiKSrVILj2MoCRt2p2tzdd4i7fhQcS1s0bcEqHpkwl9c+X0ZVwm/zTl1zOW9kEWcfshvdOmYnP4FIIyk4FpF0lUpwXEoK+YvdvbChjWpLFBxLW7ZgxXoefX8e46YWszFhZb28nEy+M7wvFxxWRN+eeTWcQaRhFByLSLqqd3AsySk4lvZg9YYtPDFlIY9NnE/Jus3blWUYnDhkFy4+YgD79+3eSi2U9kbBsYikKwXHjaTgWNqTzRWVvPjJEh6ZMJf/Lq++xs/w/j24+IgBHLf3TmRkWCu0UNoLBccikq4UHDeSgmNpj9yd8bNKGTt+Lu/PLq1WXlTQiQsPL+L0A/vQMSezFVoobZ2CYxFJVykHx2Y2DLgCGAEUAk+7++VR2dHR/rHunmT92vZHwbG0dzOWrOWRCXN5afoSKhLu3uuRl833R/Tn3BH9KOic20otlLZIwbGIpKuUgmMzuxa4DYgNFTnwZ3e/ICr/GvAv4HJ3f7iJ25qWFBzLjmLpmo08NnE+T05ZyLpNFduV5WRlcPqBu3Lh4QPYvVfnVmqhtCUKjkUkXaWSreJrwKvAcuCnwHhgDvBYXHCcAXxFWGb65GZpcZpRcCw7mrLNFTz90SIefX8ei1dXX03+mL16cfERAzh0QE/MNC9ZklNwLCLpKpVcxFcDW4AT3H06UO0Pn7tXmdl/gUFN1kIRSSudc7O48PAifjCiH69+voyx4+cyvXjN1vK3Z37F2zO/Yt9du3LxEQP4nyG7kJ2pRUVERKRtSOUv1nBgSiwwrkUxsEvDmyQibUFWZgYn79ebF644jHE/HMHxg3ci/vPyZ4vXctXfPuHIO9/hkQlzWbepvPUaKyIiUk+pjBx3ApbVo14eYRU9EdkBmBkHF/Xk4KKezC0p44/vz+PZacVsrggryy9Zs4lfv/IFv3tzFt89uC/nH1ZE7+4dW7nVIiIiyaUy53g+sMrdD4jbV0XcnONo32yg3N33buK2piXNORapbkXZZv46eSGPT5rPivVbtivLyjBO2i8sKrLvrt1ap4HS6jTnWETSVSrTKt4D9jOzUTVVMLNTgQHAW41tmIi0Xfmdc7nquD344KfHcPu3hjCgsNPWsooq58VPlnDyfe9z1pjJvD1zOVVVyrcuIiLpIZWR4yHANGAd8CPgBaAMeAy4BPgm8DDQGdjP3Wc1Q3vTjkaORepWVeW88+VXjJ0wl8lzV1Yr371XZy46vIhTD9iVDtlaVGRHoJFjEUlXqeY5voAQAGcAFYQ5y5sJc4xzgCrgInf/c9M3NT0pOBZJzafFaxg7YS6vfLqUyoQR44LOOZw7oj/fO7QfPTvltFILpSUoOBaRdJVSfiV3fxQ4nJDvuIoQFHcgLAryLnDMjhQYi0jqhvTpxu/POoDx1x3NRYcX0Tl3233BpWVbuOeN/zJy9Fv8/IVPmVe6Hp54Avr3h4yMsH3iiVZru4iItH8pLx+99UCzHKA3ITBe6u4bmrJhbYVGjkUaZ+2mcv724UL+9MF8lq7ZtF3ZKTPe4c7X7id3S9z+vDwYMwbOOaeFWypNSSPHIpKuagyOzexR4P1otFhqoOBYpGmUV1bxyn+WMnbCXD5fshaA9x88nz5rS6rVrezbl8yFC1u6idKEFByLSLqqLc/xedFWwbGINLvszAxOPWBXTtm/N5PmrGDshLn0XluatK4tKubYu99l5MACRgzM59AB+ZqjLCIiTSKVRUBERJqdmTFy9wJG7l5A+U19yCheVK3Okq4FzClZz5yS9fxl8gIA9t6lKyMG5DNyYD4HD+hJ1w7ZLd10ERFpBxQci0jayh59O1xyCWzYdkvD5pxc7jnqvGp1v1i6li+WruXRD+aRYTBk126MiEaWh/fvQV6O/rsTEZG66a+FiKSv2E13N94ICxfCbruRe+ut3Hbmdzl9wSomzVnBxDmlTC9es11auCqH6cVrmF68hofem0N2prF/3+6MGJDPiIEFHLBbd+VTFhGRpGq7Ia/a0tBSnW7IE2l9ZZsr+Gj+SibNWcGkOSv4bMkaakvEk5uVwUH9ejByYAiW9+vTjezMlDJbSiPphjwRSVd1BcdlQPI7Ymrn7j6wMQ1rKxQci6SfNRvKmTxvxdZg+cvl62qt3yknk+FFPaM5ywUM7t2VzAxrodbumBQci0i6qis4bih39x3iO0sFxyLpr7RsM5PnrmDinBVMnrOCuaXra63ftUMWh0Q3940cWMCgnTpjpmC5KSk4FpF0VVdw/CpwR0NO7O7vNaJdbYaCY5G2Z+majVtHlSfOWcHi1RtrrZ/fKYdDB4ZgecSAfIoKOilYbiQFxyKSrjTnuJEUHIu0be7OopUbmTS3lIlRsFyybnOtx+zctQMjB+ZvDZj79Mhroda2HwqORSRdKThuJAXHIu2LuzOnZD2T5oRgefLcFazaUF7rMX17dmTkgAJG7h5Glnt17VC90hNPbJd1g1tv3aGXwFZwLCLpSsFxIyk4FmnfqqqcmcvWMWnuCibNKWXK3JWs21xR6zEDCzttv3rfC89Uy9dMXh6MGbPDBsgKjkUkXSk4biQFxyI7lorKKj5fsjaaglHK1Pmr2FheWesxU8ZcyE6rllcv6NcP5s9vnoamOQXHIpKutAiIiEgKsjIzGNq3O0P7dueyowaypaKK6cWrmTh7BZPmlvLxgtVsqdw+2U/hqq+SnqtqwUIufXwqBV1yKeiUQ0GXXPI75VLQOYf8zrkUds6la8cs3fwnItKCagyO3V0Z8UVE6pCTlcHw/j0Z3r8nV7EHm8ormZawet+SrgX0WVtS7dglXQt4fUaSEeU42ZlGfqdc8jvnUNA5bAs7xz8PwXRB51x6dsrRYiYiIo2kkWMRkSbUITuTw3Yv4LDdC4A9KdtcwYKev6LXz68mZ/OmrfU2ZOVy56hz6zxfeaWzbO0mlq3dVGddgO552SFojkaiCzptH0THRqQLuuSQl9PAPwG6uVBE2jEFxyIizahzbhb7XHMZ7NIVbrwRX7iQqj59WXvDLzj766fw9bItlJZtZkXZZkrKtrCibHN4vn4Lpes2s35L7fOZE63eUM7qDeXMrkfdjtmZFHTJiaZybBuBrj5KnUv3jtlkZFgIjONvLlywIDwHBcgi0i7UeEOe1I9uyBOR5rRxS+V2wfKK9ZspjQLq0vhgumwLKzdsobn+S8/MMHp2yuEfd3+vSW4u1A15IpKuNHIsIpLGOuZk0rdnHn171r3QSEVlFas2lG8NlkujwLk0yYh0admWajcO1qayyilZt7nGmwtZuLDe5xIRSWcKjkVE2omszAwKu+RS2CW3zrruzrrNFVuD6GrTOrbu30JJ2WbWbQq5nWu6uZDddmvqyxERaRUKjkVEdkBmRtcO2XTtkE1RQac662+uqGRF2RYqdr2NyuuuJHPTxm2FeXnhpjwRkXZAOX9ERKROuVmZ9O7ekd1+dBGZj4wNc4zNwnYHXulPRNofjRyLiEhqzjlHwbCItFsaORYRERERiSg4FhERERGJKDgWEREREYkoOBYRERERiSg4FhERERGJKDgWEREREYmYu7d2G9o0MysBFjTiFAVAaRM1R5qG+iQ9qV/ST2P6pJ+7FzZlY0REmoKC41ZmZlPdfVhrt0O2UZ+kJ/VL+lGfiEh7pGkVIiIiIiIRBcciIiIiIhEFx61vTGs3QKpRn6Qn9Uv6UZ+ISLujOcciIiIiIhGNHIuIiIiIRBQci4iIiIhEFBy3AjM728wmmNkaMyszs6lmdoWZqT+agZllm9mxZna3mU02s6VmtsXMFpvZs2Z2VB3Hq79aiJndZmYePa6ppZ76pJmZWUczu87MPjKz1Wa2wczmmdkzZnZYkvoZUR9MjfpkTdRHZ7VG+0VEGkpzjluYmd0PXA5sAt4CyoFjgS7A88AZ7l7Zei1sf8zsOOCN6OkyYBqwHhgM7Bvtv8Xdb0pyrPqrhZjZcGAS4UO7Ade6+11J6qlPmpmZFQGvA7sDXwGTgc1Af2B/4GZ3/3Vc/Uzg78A3gbWEfskl9EsucJ+7X9mClyAi0nDurkcLPYDTAQeWAnvE7d8JmBGVXdXa7WxvD+AY4FngiCRl3wEqovf+aPVXq/VRLvA5sJgQ4DpwTZJ66pPm74tOwOzovbwZyE4ozwcGJez7SVT/c2CnuP17ED6QOnBKa1+bHnrooUd9HvoKsmXdEG2vd/dZsZ3uvhy4LHr6U3013LTc/W13/7a7T0hS9jTwWPT0ewnF6q+WczNhJP9SYE0t9dQnze/nwEDgcXe/yd3L4wvdfYW7/zf2PBo1vi56elnUF7G6s4Dro6c3Nm+zRUSahv6AtBAz6wMcBGwBnkksd/f3CKNmOwOHtmzrdnj/jrZ9YjvUXy3HzA4hjDw+6e4v11JPfdLMzCwHuDh6Orqeh40AegHF7j4+SfkzhKkvw81s18a3UkSkeSk4bjkHRNvP3X1jDXU+SqgrLWOPaLs0bp/6qwWYWQfgz8BK4Ko6qqtPmt9BhGkTi9z9CzMbGd0k+bCZ/crMRiQ5JvZef5SkDHffQJhuAWG+sohIWstq7QbsQIqi7YJa6ixMqCvNzMx2Bs6Lnj4XV6T+ahm3AnsC33X30jrqqk+a35BoO8vMHgN+kFB+k5k9B3w/7gNKfftlf9QvItIGaOS45XSOtutrqVMWbbs0c1sEMLMs4K9AN+CthK/01V/NzMxGAv8LvBDN/a6L+qT59Yy2o4BzgbsIGSt6AKcQpq2cDtwfd4z6RUTaFQXHLceirXLnpY+HCKmmFlH9Zjz1VzMys47Anwhpvy6v72HRVn3SfGJ/E7KAP7r7te4+x91Xu/tLwKmE9/8HZjYgqqt+EZF2RcFxy1kXbTvXUidWtq6WOtIEzOx3wIWENFPHuvuyhCrqr+Z1GzAIuNrdl9ZVOaI+aX7x79vYxEJ3n0rIE54BHJVwjPpFRNoFzTluOfOjbb9a6vRNqCvNwMzuBq4ESgiB8awk1eZHW/VX8zgNqCKMQCbOa90r2l5mZicDs939ItQnLWF+3L/n1VBnHjCMkBUk/hj1i4i0CwqOW04sXdg+ZtaxhrvthyfUlSZmZncCVwMrgOPdfUYNVdVfzS8DOLKW8gHRo3v0XH3S/D6O+3c+4QNkooJoG5tHHDtmeJK6mFke21aiVL+ISNrTtIoW4u6LCH9EcoAzEsvN7EhCnt1lhCV0pYmZ2WjgWmAVITCeXlNd9Vfzcvf+7m7JHoTUbhCWjzZ33z86Rn3SzNx9MTAlenpsYrmZ9QAOjJ5OjbaTCEtM9zGzUUlOewaQDXwUnV9EJK0pOG5Zt0fbO8xs99hOM+sFPBA9He3uVS3esnbOzG4hrNS1mhAY12cES/2VftQnze/WaHuTmW3NSxzlpH6QkN1lGtEHEHevBH4TVXsw6ovYMXuwbTGR2HlFRNKauesG45ZkZg8QlrndBLxJWDnqWKAr8ALw7eiPjTQRM/sm8GL0dCrbFiRINNPdt1sVTP3V8uLy617r7nclKVefNDMz+w1wDWE1wimEaUgHA70J6dyOjp+rHy0h/TzwDUIGkrcIo8XHAR2A+9z9ypa8BhGRhlJw3ArM7GzgCkLC/UxgJvAo8KBGvJqemZ1HSBtWl/fc/agkx6sEZ1gjAAAM3klEQVS/WlBdwXFUR33SzMzsNOBHhBXw8ggLebxEGJmvNhfZzDIIafnOJ9xUWQn8B3jA3Z9sqXaLiDSWgmMRERERkYjmHIuIiIiIRBQci4iIiIhEFByLiIiIiEQUHIuIiIiIRBQci4iIiIhEFByLiIiIiEQUHIuIiIiIRBQc76DMzM2sWpJrM5sflfVv+VZJovbcHxb81Mw+N7NN0XWujsrOi54/1sJtuj563RNa8nWbipnlmdlSM/vIzKy12yMi0hYpOJYmZ2b9owBjfmu3RdLaFcDtwK7AK8CfgVZbSc3MdgFuBMa7+6ut1Y7GcPcNwK3AMODcVm6OiEibpBXydlCxUWN3t4T9A4FsYI67lzfw3P2BecACd+/fqIbu4KIPGP2AInef37qtaVpm9h4wCviau7+RUNYN2AVY4+5LW6g9Y4CLgWPd/e2WeM3mYGY5wAKgivBzs6WVmyQi0qZo5Fi24+5z3H1mQwNjkRT0jbazEgvcfU30c9hSgXE+8H1gLvBOS7xmc4mC4b8CvYEzWrk5IiJtjoLjdszMhpjZ82a20szWm9nHZnZRHcckneNqZt3N7LZofugGM9toZsVm9q6Z3RBX7zHCqDFAv9jc5sRpFmZWaGZXmdmrZjYvmnO6xswmm9kVZpaZpG1bp2tE81UvN7NPovasMrMXzWzfWq4t38xuNrN/m9na6D2ZZWaPmdnIJPU7mdl10fzNtdE1f25mvzSzzrW9j0nO1cXMLjGzF8xsdtTmsqgtN5pZx1TOF9e+G81senSu9dH78TMzy0tS/6jo/XvXzLKjY2dG7/1XZvZXM9utltc73cwmRq+zysxeN7Mj4s9bz3a/G31zURTtmhf3M3JeVCfpnOPGXkMtLgA6AI97kq/TYm2OXv+w6Od2VfQz+5qZ7R9X99zoZ6Ys+t37q5ntnOScW6/RzHqY2e/NbGH0c/aFmV0aV3cfMxtnZsuj8g/N7Ou1XM+fo+3lDXgvRER2aFmt3QBpHmZ2JPAvoCPwJfBvwtfUD5vZ4BTPlQd8AAwGvgLeBNZH5xsMHEqYOwrwPtAZOD2q82zcqUrj/v114F6gGJgNTAF2AkYAhwDHm9lpyQKVyGPAd4DxhJHH4cA3gaPM7AB3n5twDQcQ5rXuAqwE3gU2EaYsnBVVmxhXvw/wWnR9JcCkqP5w4BfAaWZ2lLuvqqF9iYYCDxPevy+BqUB+dK2/Br5pZke6+6b6nMzMCoC3gSHAKuANwIGjCXNOzzSzY9x9ZZLDswk/G4cA7wFfEN73c4BRZrafu69OeL2fRed1wvu0ENiHMMr6+3q+BzGvAvOBbwOdgOeAsqhsdj3PkfI11OHUaPtmHfW+AVwFTCP8fAwFvgYcambDgB8CV0Zteg04LGrT/mZ2YA1THLoTfr66En5/8gnTTR60ML1kPPA6YarEO8AehJ/DV6I+Hp94Qnf/zMyWAyPMrNDdS+r3NoiICO6uRzt7EALiYkIgcxvR3PKo7EhC0Oqh+6sdOz8q6x+379xo3z+ArIT6mcAxCfv6R/Xn19LGvYFDkuzfhRDIO/CdGs7rhK+/B8aV5RKCXwfGJhzXBVgUlT0IdEwoLwQOj3tuhADQgfuAvIT39i9R2WMp9Ekf4BggI2F/d0KQ58D19emPaP+4aP94oHvc/h6EDzIOPJVwzFFx799HQK+4sm6EgM+BGxOOOwioBLYAJySUXRl3zndT/DlNem1R2XnJ3uOGXkMd7ciLrm0L0KGGOu9G560Cvh23P4NwE6EDnwLLgMFx5T0JH4Yc+H4N1+jAM/GvDZwY7V8XvU8/STj2N1H5W7Vc1/NRnTNT6Rc99NBDjx39oWkV7dO3CRkA5gD/5+5bR1/d/T3goRTPt1O0fdPdK+IL3L3SG3Dzkrt/4e5TkuxfClwXPf12Lae40t3nxB23GfhV9PTYhLoXEoLTycDl7r4x4TVL3P39uF0nEEYhJwNXecgAEKu7EbiUMAJ8jpn1qKWN8a9R7O5vu3tVwv7VhAATar/ercysX1S3CrjE40ZIPYxkXxyVnWlmfZOcwoEL3P2ruOPWAHdETxPfvysIQeDjnpDFwd1/Txj1b2mpXkNt9iGMRM/zukfu/+buW78NifrzzujpvsBN7j4jrnwl237fjq7hnOuAy+Jf293/BUwnfAuzxN3vTjgm9k3N4WaWXcN5Y+04oPZLEhGReJpW0T4dGW3/5u6VScr/Alydwvk+jLbXm1kp8A9P7SvrpMwsizCaOgLYmTDn0wgjvQCDaji0gvDVfKKZ0bZ3wv5Yzto/xn9QqMX/RNvnEoNZAHdfb2ZTo3rDCV9518nMjPA1+yhCsN6RcL2xjCE1XW+iI6JjJrn7zMRCd59hZh8SpruMAp5IqLLQ3T9Nct6a3r/Yz1NNadaeIkxvaEmpXkNtekXbFfWom+znbnYd5bEbDmtq01R3L02yfzZh2ka1c7r7SjNbQZiCkU8YsU4Um1KzU5IyERGpgYLj9qlPtJ1XQ/n8VE7m7u+Z2Z3ANURTCsxsJmF+5HPu/lqqDTSzQcALhOkVNelaw/6liSPYUTvXhviT3ISiftG2WiBZgwHR9jdm9ps66hbW54RmthPwd6DajX9xarreRLtG25r6F8K3BofG1Y23sIZj1kbbDjW83oIajqtpf3NK9Rpq0y3h2NoUJ+5w9zLbtt5GtXK2zaeuqU3Jjok/rrby/FrOG7ue7jWUi4hIEgqOd0wpJ7d29+vN7CHgFOBwwgjoxcDFZvY6cFKygLUWzxIC45cIX0t/QchpWxkFzl+ybUQ1UbXR3CYWy5TxHnV/kKhvYPgIITD+APgl4Svz1e5ebiEv7eYU2hd7X2rrx9pWR2vo+1fT6zV3fzT3a8a+BanPh5NaXzfZNw2NPWc9ymsSu5763jQqIiIoOG6vFkfb/jWUF9Wwv1buPo+QYeJeADM7nPCV+tcIqbDG1Oc8ZrYXIcvCV8C3kkz92L0h7avFAmAvYE/CaHddFkXbZ9z9/sa+uJl1IkzBqAROTjIlJdXrjY0kDqilTqyPF9dSp76WROfrR7gRMlH/JniN1hSbt5zfqq1oerHr+arWWiIish3dkNc+vRdtv2tJ8gUTUks1WnQT22PR06FxRbF0VTV9+OoZbZfUMCe6SdoXJzbt4wKL+/67Fv+Ktk21gEI3wu/auhrmaqd6vRMIo7iHRqPs2zGzvQlzgKsI2SwaK3aOs2oo/24TvEZr+pwwcl9kDcg3ncZiKRs/btVWiIi0MQqO26dngaWEEclfxgeE0WjvZamczMxOM7NRZpaRsL8jcFz0NH56QQkhQN6phmwOswiB275mNirhnOdTcxDWUI8QRj9HAveZ2XZzNC0sSHJ43K4XCCnBjjSzh8ysJwnMbICZXVHP119O+Gq7u5mdnXCeE0jt5kjcfQEhN3AGIW91bM4sZtadkE85Axjn7ouSnyUl9xP66wdmdnx8gZldTrihss2KMpBMIWSsOKiVm9OUDiVKsdfK7RARaVMUHLdDUeqx7xEWrfg5MMPMnjSzdwijyvWa/hDnyOi4ZdFqYH81s5cJX+8fSrjR7eG41y8n5BzOAv5tZk+Y2SNmNjoqLwEeiMrfMbO3o/Z9CjwKjG7wxSfh7usIc6W/IqQlK7awmt7TZjY5uo6L4upXERaF+JSwqMN8M5tgZk+Z2Rtm9iVRmrx6vn4lYQENgCcsrDL3pJlNIYxS39OAy7oM+IyQ93eumT1nZs8Rpj0cQZjTXN/gvVbu/hFhnnQO8Fr0XjxhZp8Q8kD/LqqabIGLtuKFaHtcrbXaCDMbQshSMcm1AIiISEoUHLdTUe7hQwk3vO1MCPZ6AFe4e0ojlYSpE3cA/yXkcj0DOJiQaurHwMFRjtl4FwN/JNzcdiYh13D81+9XAZcQgriDCYseLI+2qQbvdXL3qYR5zqMJaa+OB04iTPF4koTcz+5eHLXr/xEWJdmHsOrfvoS8tHcB30rh9e8m5CaeHJ3rZMIc5O+5+40NuJ5Swojt/xHmFZ8YPRYBNwKHefLV8RrE3W8h9OMU4EDCe1dKyCc8NaqWLB1ZW/EYsBE4t55Tb9LdD6LtA63aChGRNsjql/ZVRCQ5M/sj4YbMa5IsVtFmRNlYfggc25CFbdJFlP1kAWEqTJEnX7JaRERqoJFjEamTmQ2K5jPH77Nojvj5hBvanmqVxjWdXxK+FfhFK7ejsS4hfFv0MwXGIiKp08ixiNTJzH4NXEuYYrKIsLrfYEKKtyrgUncf23otbBpmdh1hCtGJiUtltwVmlkeYD19MmO6k/+BFRFKk4FhE6mRmI4EfEVLEFRJWISwBJgL3uvsHrdg8ERGRJqPgWEREREQkojnHIiIiIiIRBcciIiIiIhEFxyIiIiIiEQXHIiIiIiIRBcciIiIiIpH/D7Gktq2SDBSlAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "L=60e-3\n",
    "s=np.sqrt(hp)\n",
    "F=lambda x: 20+80*(np.cosh(s*L-s*x)+h/s/k*np.sinh(s*L-s*x))/(np.cosh(s*L)+h/s/k*np.sinh(s*L))\n",
    "x=np.arange(0,N+1)*dx\n",
    "plt.plot(x*1000,F(x),label='analytical')#a*np.cosh(s*x)+b*np.sinh(s*x))\n",
    "plt.plot(x[1:]*1000,T,'ro',label='finite difference')\n",
    "plt.plot(x[0],100,'rs',label='base temperature')\n",
    "plt.xlabel('distance along fin (mm)')\n",
    "plt.ylabel('Temperature (C)')\n",
    "plt.legend(bbox_to_anchor=(1,0.5),loc='center left');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise\n",
    "\n",
    "Decrease $\\Delta x$ from $10~mm$ to $5~mm$. \n",
    "\n",
    "a. Make a plot of the Temperature along the fin. \n",
    "\n",
    "b. Plot the heat flux through the fin $-\\kappa \\frac{dT}{dx}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Static beam deflections\n",
    "\n",
    "In the [Module 04 project](https://github.uconn.edu/rcc02007/CompMech04-LinearAlgebra), we solved a finite element analysis problem with stretching beams. A beam is very stiff in tension and compression, but it is usually much less stiff if bent. The [Euler-Bernoulli](https://en.wikipedia.org/wiki/Euler%E2%80%93Bernoulli_beam_theory) beam theory relates the static deflection of a beam, $w(x)$, to the bending stiffness, $EI$,  and applied load, $q$ as such\n",
    "\n",
    "$EI \\frac{d^4w}{dx^4} = q.$\n",
    "\n",
    "If $q=q(x)=constant=q_0$, then the equation can be integrated $\\times4$ to create a polynomial of order 4\n",
    "\n",
    "$w'''(x)=q_0x+A$\n",
    "\n",
    "$w''(x)=\\frac{q_0x^2}{2}+Ax+B$\n",
    "\n",
    "$w'(x)=\\frac{q_0x^3}{6}+\\frac{Ax^2}{2}+Bx+C$\n",
    "\n",
    "$w(x)=\\frac{q_0x^4}{24}+\\frac{Ax^3}{6}+\\frac{Bx^2}{2}+Cx+D$\n",
    "\n",
    "where $A,~B,~C,~and~D$ are integration constants. We have four unknowns, so we need four boundary conditions. \n",
    "\n",
    "Let's consider the [simply supported beam](https://www.efunda.com/formulae/solid_mechanics/beams/casestudy_display.cfm?case=simple_uniformload) with a uniform load. We are considering linear beam mechanics, so we have an underlying assumption that the tension in the beam is negligible. The loads and boundary constraints are shown at the top of the figure below. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "![Boundary conditions for 4 different beam supports](../images/beams_BCs.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "solving for A, B, C, and D there results\n",
    "\n",
    "$w(x) = -\\frac{qx\\left(L^3-2x^2L+x^3\\right)}{24EI},$\n",
    "\n",
    "but we can use the same finite difference method we used above to solve for the deflection $w(x_i)$ at different locations, $x_i$. We begin creating the numerical solution with the same process we did with the heat transfer problem, we estimate the _differential equation_ as a _finite difference_ equation and use the following material and geometry (1-m steel rod 1-cm-by-1-cm) with 100 N/m load applied\n",
    "\n",
    "* $L=1~m$\n",
    "* $E=200e9~Pa$\n",
    "* $I=\\frac{0.01^4}{12}~m^4$\n",
    "* $q=100~N/m$\n",
    "\n",
    "__Finite difference equation dividing bar length $L$ into 6 segments:__\n",
    "\n",
    "* $\\frac{d^4w}{dx^4} \\approx \\frac{w(x_{i+2})−4w(x_{i+1})+6w(x_i)−4w(x_{i-1})+w(x_{i-2})}{h^4}=\\frac{q}{EI}.$\n",
    "\n",
    "1. $w_{-1} - 4w_0 +6w_1-4w_2+w_3 =\\frac{q}{EI}h^4$\n",
    "\n",
    "2. $w_{0} - 4w_1 +6w_2-4w_3+w_4 =\\frac{q}{EI}h^4$\n",
    "\n",
    "3. $w_{1} - 4w_2 +6w_3-4w_4+w_5 =\\frac{q}{EI}h^4$\n",
    "\n",
    "4. $w_{2} - 4w_3 +6w_4-4w_5+w_6 =\\frac{q}{EI}h^4$\n",
    "\n",
    "5. $w_{3} - 4w_4 +6w_5-4w_6+w_7 =\\frac{q}{EI}h^4$\n",
    "\n",
    "__Boundary Conditions:__\n",
    "\n",
    "* $w(0)=0=w_0$\n",
    "* $w(L)=0=w_6$\n",
    "* $w''(0)=0=\\frac{w_{-1}-2w_{0}+w_{1}}{h^2} \\rightarrow w_{-1}=-w_{1}$\n",
    "* $w''(L)=0=\\frac{w_{5}-2w_{6}+w_{7}}{h^2} \\rightarrow w_{7}=-w_{5}$\n",
    "\n",
    "__Final linear equation set:__\n",
    "\n",
    "1. $5w_1-4w_2+w_3 =\\frac{q}{EI}h^4$\n",
    "\n",
    "2. $-4w_1 +6w_2-4w_3+w_4 =\\frac{q}{EI}h^4$\n",
    "\n",
    "3. $w_{1} - 4w_2 +6w_3-4w_4+w_5 =\\frac{q}{EI}h^4$\n",
    "\n",
    "4. $w_{2} - 4w_3 +6w_4-4w_5+w_6 =\\frac{q}{EI}h^4$\n",
    "\n",
    "5. $w_{3} - 4w_4 +5w_5-4w_6 =\\frac{q}{EI}h^4$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\left[\\begin{array}{cccccc}\n",
    "5 & -4 & 1 & 0 & 0  \\\\\n",
    "-4 & 6 & -4 & 1 & 0 \\\\\n",
    "1 & -4 & 6 & -4 & 1 \\\\\n",
    "0 & 1 & -4 & 6 & -4 \\\\\n",
    "0 & 0 & 1 & -4 & 5 \n",
    "\\end{array}\\right]\n",
    "\\left[\\begin{array}{c}\n",
    "w_{1}\\\\\n",
    "w_{2}\\\\\n",
    "w_{3}\\\\\n",
    "w_{4}\\\\\n",
    "w_{5}\n",
    "\\end{array}\\right]=\n",
    "\\frac{qh^4}{EI}\\left[\\begin{array}{c}\n",
    "1\\\\\n",
    "1\\\\\n",
    "1\\\\\n",
    "1\\\\\n",
    "1\n",
    "\\end{array}\\right]$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "finite difference A:\n",
      "------------------\n",
      "[[ 5. -4.  1.  0.  0.]\n",
      " [-4.  6. -4.  1.  0.]\n",
      " [ 1. -4.  6. -4.  1.]\n",
      " [ 0.  1. -4.  6. -4.]\n",
      " [ 0.  0.  1. -4.  5.]]\n",
      "\n",
      "finite difference b:\n",
      "------------------\n",
      "[-0.00046296 -0.00046296 -0.00046296 -0.00046296 -0.00046296]\n",
      "\n",
      "deflection of beam (mm)\n",
      "-------------\n",
      " [-4.05092593 -6.94444444 -7.98611111 -6.94444444 -4.05092593]\n",
      "at position (m) \n",
      "-------------\n",
      " [0.16666667 0.33333333 0.5        0.66666667 0.83333333]\n"
     ]
    }
   ],
   "source": [
    "L=1\n",
    "h=L/6\n",
    "E=200e9\n",
    "I=0.01**4/12\n",
    "q=100\n",
    "\n",
    "A=np.diag(np.ones(5)*6)\\\n",
    "+np.diag(np.ones(4)*-4,-1)\\\n",
    "+np.diag(np.ones(4)*-4,1)\\\n",
    "+np.diag(np.ones(3),-2)\\\n",
    "+np.diag(np.ones(3),2)\n",
    "A[0,0]+=-1\n",
    "A[-1,-1]+=-1\n",
    "\n",
    "b=-np.ones(5)*q/E/I*h**4\n",
    "\n",
    "w=np.linalg.solve(A,b)\n",
    "xnum=np.arange(0,L+h/2,h)\n",
    "print('finite difference A:\\n------------------')\n",
    "print(A)\n",
    "print('\\nfinite difference b:\\n------------------')\n",
    "print(b)\n",
    "print('\\ndeflection of beam (mm)\\n-------------\\n',w*1000)\n",
    "print('at position (m) \\n-------------\\n',xnum[1:-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f40a40c6510>]"
      ]
     },
     "execution_count": 65,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi40LCBodHRwOi8vbWF0cGxvdGxpYi5vcmcv7US4rQAAIABJREFUeJzt3Xl4VNX9x/H3NwkBFGQxiMqqICqgogZ3BATUuu+Ke2tL3Vrrvu/aWrW2+nOlLlXrvoKIC6gIIghRWRUEEWWzBNlkDUnO748TFDGQCXNnzsydz+t55iEzGe793EzmmzPnnnuOOecQEZH4yAsdQEREoqXCLiISMyrsIiIxo8IuIhIzKuwiIjGjwi4iEjMq7CIiMaPCLiISMyrsIiIxUxBip0VFRa5t27Yhdi0ikrU+/fTTBc65ZjU9L0hhb9u2LSUlJSF2LSKStczs20Sep64YEZGYUWEXEYkZFXYRkZhRYRcRiRkVdhGRmIlkVIyZHQrcC+QDjzrn7ohiu2sV3zaEBcvKfvV4UYNCSq7rE+WuREQil+4alnSL3czygQeA3wAdgb5m1jHZ7a6ruh/Ixh4XEckk6a5hUXTF7AVMd87NcM6VAc8DR0ew3V/ZN28y5+UPTMWmRURSqi5l3FjwJC0oTfm+oijsLYBZ69yfXfXYL5hZPzMrMbOS0tJNO7AeeeO4rOAF2tq8TUsqIhJI3/z3+W3BO2xrP6R8X1EUdqvmsV+tkO2c6++cK3bOFTdrVuMVsdX6d/kRrKGAC/IHbNL/FxEJoS5lnFvwBqMqOjLW7ZTy/UVR2GcDrda53xKYG8F2f2UBjXimojfH5n9EK/tfKnYhIhK5k/KHsbUt4r6KY9OyvygK+1hgBzPbzswKgVOASDvCixoU/vT1I+VHUEE+5+cP/MXjIiIZqXw1F9R5gzGVOzKq8pfjSlJVw5Ie7uicKzezC4F38MMdH3fOTU462Tp+NRxo8Hj6ljxO3wvvi3I3IiLRG/cMW/MDW5/1KDPbHZSWXUZygZJzbrBzroNzrp1z7vYotrlR+/8FLA8++mfKdyUissnKy2DEPdCyK2zfM227zc4rTxu1gN1Ph8+ehiWzQ6cREanehOdhySzofiVYdeNMUiM7CzvAARcDDkbeGzqJiMivVayB4XfDtrtD+95p3XX2FvbGraHLqfDpk7BU49pFJMNMfAkWf5v21jpkc2EHOOASqCyHj3USVUQySEW5b61vvQt0ODTtu8/uwt50O9jtFCh5HH7UuHYRyRCTX4WFXwdprUO2F3aAbpdCRRmM+r/QSUREoLICht8FW3WCHQ8PEiH7C/uW7WCXE2HsY7B8Qeg0IpLrvngdFnwF3S+HvDAlNvsLO0C3y2DNSvhYrXYRCaiyEj68C4p2hJ1TMsltQuJR2Jt1gM7HwZh/w/LUz5wmIlKtLwdC6ZfQ/YpgrXWIS2EHOPAKWLMCRj8QOomI5KLKSvjwTijqAJ3SM9nXhsSnsG+1E3Q6Bj7pDysWhk4jIrlm6pswfzIceDnk5QeNEp/CDv4HWvYjjH4odBIRySXOwYd/h6btoNNxodPErLA37wQ7HwWfPAwrF4VOIyK5Yupg+H6ib1zmJz1pbtLiVdjBn7RYvRQ+eSR0EhHJBWtb602280OvM0D8CvvWu8BOR8DoB2HVktBpRCTuvnoH5o2HAy/LiNY6JFnYzexEM5tsZpVmVhxVqKR1v8IXdbXaRSSVnIMP74DGbWDXk0On+UmyLfZJwHHA8AiyRGeb3aDDb2DUA7Bqaeg0IhJX04fC3M/91Cb5dUKn+UlShd0596VzbmpUYSLV/QpYtRjG9A+dRETiyDkYdgc0ag279Q2d5hfi18e+Vos9YIeDYdT9sPrH0GlEJG6+fg/mlEC3i6EgNYtSb6oaC7uZDTWzSdXcajURgpn1M7MSMyspLS3d9MS10f0qP+xxzL/Tsz8RyQ3O+atMt2gJXU4LneZXajyF65yLZE0n51x/oD9AcXGxi2KbNWq5J7Tv4ycH26sf1G2Qlt2KSMzN+ABmfQKH3wMFdUOn+ZX4dsWs1eMqWLkQxqrVLiIRWNu3vkVL2P300Gmqlexwx2PNbDawL/Cmmb0TTawItSz2rfaR98HqZaHTiEi2W9ta73ZJRrbWIflRMa8551o65+o655o75w6JKlik1GoXkShkQWsdcqErBtRqF5FoZEFrHXKlsINa7SKSnJ9a6y0yurUOuVTYWxZD+95qtYvIpsmS1jrkUmGHqnHtarWLSC39orV+Rug0Ncqtwt6qq1rtIlJ7WdRah1wr7KBWu4jUTpa11iEXC7ta7SJSG1nWWodcLOygVruIJCYLW+uQq4X9F612zfwoIhvw9ftZ11qHXC3sAD2u8a12rbIkItVxDj74KzRqlVWtdcjlwt5yT+hwqJ/5UWujisj6pg3x860feFlWtdYhlws7QI+r/SpLox8OnUREMolzMOyvfi3TDJxvvSa5Xdi37QI7Hu7XRl25OHQaEckUX73t1zI98PKMWss0Ubld2MHPIbN6CYx+MHQSEckEa/vWm2wHu50SOs0mUWHfZlfY+SgY9SCsWBg6jYiENuVN+H4CdL8yK1vroMLu9bgaypb5ha9FJHdVVsKwv8GW7WGXE0On2WTJrqB0l5lNMbMJZvaamTWOKlhaNe8InY7xQx+X/xA6jYiE8uVA+N+kqtZ6jUtCZ6xkW+xDgM7OuV2Br4Crk48USPeroGw5fHxf6CQiEkJlpb/KtKgDdD4+dJqkJLs03rvOufKqu6OBlslHCmSrnWCXE2BMf1hWGjqNiKTbF69B6Zd+QEVefug0SYmyj/13wFsRbi/9ul8J5atg5L9CJxGRdKqsgGF/h2Y7Q8djQ6dJWo2F3cyGmtmkam5Hr/Oca4Fy4JmNbKefmZWYWUlpaYa2iIt2gF1PhrGPwo/fh04jIuky8WVYMLWqtZ79Y0rMOZfcBszOAs4FejnnViTyf4qLi11JSUlS+02ZhTPg/q5Q/Ds47K7QaUQk1SrW+Pd83YbQ78OMLuxm9qlzrrim5yU7KuZQ4ErgqESLesZrur1fqLbkCVj8Xeg0IpJq456BRd/AQddldFGvjWSP4n6gITDEzMaZWTwmXTnwcjCDD+8MnUREUmnNKv8+b9kVdjg4dJrIJDVQ0znXPqogGaVRSyg+x4+QOeBi2LJd6EQikgqfPQlL58AxD/rGXEzE43NHKhxwsZ+qc9gdoZOISCqUrYDhd0PbbrBd99BpIqXCviENm8Ne/WDiSzD/y9BpRCRqY/rD8vm+bz1GrXVQYd+4/S+CwgZ+pjcRiY9VS/31Ku37QOt9QqeJnAr7xmzWFPa9wM8fMXdc6DQiEpXRD8HKRXDQtaGTpIQKe032PR/qN4EPbg+dRESisGKhn8l1pyNg291Dp0kJFfaa1Gvku2SmvQvffRI6jYgk6+P7YPWP0POa0ElSRoU9EXv1g82bwfu3hk4iIslYNt9Pz935eGjeKXSalFFhT0Th5tDtMpg5Ar7+IHQaEdlUw++G8tV+cZ0YU2FPVPFvoVEreO8WvyaiiGSXRd9CyeN+ypCieF5buZYKe6IK6vq/8nM/gy/fCJ1GRGpr2B1geX567phTYa+N3U6Boh19X3tFec3PF5HMMP9LGP8c7N0PGrUInSblVNhrIy8fel0PC76CCc+HTiMiiXr/Nj8t7wGXhE6SFirstbXTEdBiT/jgb35mOBHJbLPGwpRBsN+f/UWHOUCFvbbMoNcNsHS2PxEjIpnLOXjvZj9ceZ/zQqdJGxX2TbF9D38bcbe/0EFEMtOMD/ww5W6XQd0GodOkjQr7pup1A6z4AUY9EDqJiFTHORh6MzRq7Ycr55Bkl8a71cwmVK2e9K6ZbRtVsIzXYk/Y+Sj4+H5YviB0GhFZ3xcDYN446Hm1H66cQ5Jtsd/lnNvVOdcFGATcEEGm7HHQdbBmOYy4J3QSEVlXRbkfCdNsJ9j15NBp0i6pwu6cW7rO3c2B3Loks9mOsNupMPZRWDwrdBoRWWv8s/DDNDjoej9MOcck3cduZreb2SzgNDbSYjezfmZWYmYlpaWlye42c/S4yv+rxThEMkPZCv9+bNkVdjo8dJogaizsZjbUzCZVczsawDl3rXOuFfAMcOGGtuOc6++cK3bOFTdr1iy6IwitcSvY+4/+qrbvJ4ZOIyKjH4Qf50GfW2O35F2iaizszrnezrnO1dwGrPfUZ4HjUxMzw3W7xM/bPvSm0ElEctvyH2DkvbDjYdBm39Bpgkl2VMwO69w9CpiSXJwsVb8JHHgZTB8KM4aFTiOSu4bfBWXLoPdNoZMElWwf+x1V3TITgIOBiyLIlJ26/sFP6zvkBqisDJ1GJPcs/MYPZNj9DD+wIYclOyrm+KpumV2dc0c65+ZEFSzr1Knnhz/OGw+TXw2dRiT3vH8r5BXEfhGNROjK0yjtchI038XPTVG+OnQakdwx51OY9ArsdyFssU3oNMGpsEcpLw/63AyLv4Oxj4VOI5IbnIMhN8JmW/oZHEWFPXLte/kJwobfBSsXh04jEn/Th/qJvrpfCfW2CJ0mI6iwp0KfW2DlQhj5r9BJROKtssIPWGiyHeyZWxN9bYwKeypss5vvbx/9ECzJ3fPJIik3/nmY/4WfbbWgMHSajKHCnioHXQeu0k9EJCLRK1vu31/b7gGdjg2dJqOosKdKkzZ+xZbxz8Lcz0OnEYmfj++HH+fCIX/N2akDNkSFPZW6XerP1L9znT9zLyLRWDrPn8PqeHROTx2wISrsqVSvEfS8Br79CKa8GTqNSHy8fxtUluf81AEbosKeanuc7Sf7H3I9lJeFTiOS/eaNh3HP+FlVm24fOk1GUmFPtfwCOPh2WDgDxv47dBqR7OYcvHOtn3iv22Wh02QsFfZ02KE3tOsFH/4dViwMnUYke00d7C9G6nkN1G8cOk3GUmFPl4Nvg9U/+uIuIrVXXgbvXg9FHWDPs0OnyWgq7OnSvCPscZafVnTBtNBpRLJPyWOw8GvftZlfJ3SajKbCnk49r4GC+v4SaBFJ3IqFMOwO2L4n7NAndJqMp8KeTg22ggMv9f2EMz4MnUYkewy/C1YvhUNu18VICYiksJvZZWbmzKwoiu3F2t7nQePW8PbVUFEeOo1I5iudCmP6w+6nQ/NOodNkhaQLu5m1AvoA3yUfJwfUqedPpM6fDJ8+ETqNSGZzDt6+CupsDgepCzNRBRFs45/AFcCACLYVe8W3DWHBsnyeqdOJTm/eSI9XG7KYhhQ1KKTkOvUdiqxVfNsQuqwYxaOF73PLmjN4/LYxAHqvJCCpFruZHQXMcc6NT+C5/cysxMxKSktLk9ltVluwrAwwbi4/kwas5NKCl9Z5XETW+nHZMq4veJpplS14quLnQq73Ss1qLOxmNtTMJlVzOxq4Fkjo85Fzrr9zrtg5V9ysWbNkc2e9r1wrnq7ow6n579HRZoaOI5JxzskfTJu8+dxcfiblkXQu5I4aC7tzrrdzrvP6N2AGsB0w3sxmAi2Bz8xs69RGjo9/lh/PYhpwY52nAM3+KPKTJXO4sGAAb1d05aPKXUKnyTqb3BXjnJvonNvKOdfWOdcWmA3s4Zz7PrJ0MbeUBtxVfjJ7503hyLxRoeOIZI4hN5BPJbeVnxY6SVbSOPbAXqzowcTKtlxd51m/IoxIrvt2FEx6mYcrjmC22yp0mqwUWWGvarkviGp7cVXU4JfrMlaSx01rzmJbWwgf/TNQKpEMUVkBb10OW7TkxbonVPuU9d9D8ms6I5Fm1Q/TOhxenQIj74Mup0HT7dKeSyQjfPYkfD8RTniCkZ2PCJ0ma6krJlP0vhnyCvxc0yK5aMVCeO9WaHOAFqdOkgp7pthiG+h+BUx9E6a+HTqNSPoNvQlWLYHf/F3zwSRJhT2T7HO+X0bvrcuhbEXoNCLpM2us74bZ5zzYunPoNFlPhT2TFBTC4ffA4u9gxN2h04ikR0U5DLoYGm4LPa4KnSYWVNgzTdv9Ybe+/kRq6dTQaURSb0x/+N9EOPRvULdh6DSxoMKeifrcCoWbwZuX+tntROJq6Vz44HZo3xs6Hh06TWyosGeiBs2g141+0d6JL4VOI5I671wDleVw2F06YRohFfZMtedvocWe/hd/5eLQaUSiN/09mPwadLsUmm4fOk2sqLBnqrw8fyJ1xQ/w/q2h04hEa80qGHwZNG0H+18UOk3sqLBnsm27QNc/wNjHYM5nodOIRGfkv2DhDDj8H1BQN3Sa2FFhz3QHXesXwR50sZ9HQyTb/fA1jLgHOh8P7XqGThNLKuyZrl4jPwxs3jgY/VDoNCLJcQ7euMi30g/5a+g0saXCng06HQcdDoX3b4OF34ROI7LpPnvKj/bqcws01Jo8qaLCng3MfF9kXgEM+ovGtkt2WjoP3r3eT/K1x1mh08RasotZ32Rmc8xsXNXtsKiCyXoatYTeN8KMYTDu2dBpRGpv8GVQsRqOus+P+pKUieKn+0/nXJeq2+AIticbUnwOtN7Xj21fNj90GpHEfTEApgzyc8Fs2S50mtjTn81skpcHR94Ha1bA4MtDpxFJzMpF/vd1611h3z+FTpMToijsF5rZBDN73MyaRLA92ZhmHfy87V+8DlPeDJ1GpGbvXg/LF8BR/wf5WrQtHWos7GY21MwmVXM7GngIaAd0AeYB/9jIdvqZWYmZlZSWlkZ2ADlp/79A885+krBVS0KnEdmwGR/C50/Dfhf6C+4kLcxFNMLCzNoCg5xzNc6SX1xc7EpKSiLZb86a8xk82gv2OBOOvDd0GpFfK1sBD+0LlgfnfQx16odOlPXM7FPnXHFNz0t2VMw269w9FpiUzPakFlrs4Vdc+vQ/vlUkkmk+uB0WzfTnhVTU0yrZPvY7zWyimU0AegIXR5BJEtXzWj+J0oALYNXS0GlEfjZzJIx6AIp/B9t1C50m5yRV2J1zZzjndnHO7eqcO8o5Ny+qYJKAws3g2Edg6Rx45+rQaUS81T/C6+dBkzZ+0RhJOw13zHatuvqTqZ//F6a+HTqNCLx7nV+395iHoW6D0Glykgp7HPS4yo+SeePPsGJh6DSSy6YN9ed99vsTtNk3dJqcpcIeBwV14diHfVF/89LQaSRXrVwEAy+EZjv78z8SjAp7XGy9i2+5T34VJr0SOo3kosGXw/JS38ioUy90mpymwh4n+/8FWhT7VvuP34dOI7lk8ut+4fUDr9CFSBlAhT1O8gt8a2nNShj4Z03vK+mxbL5f4Wvb3aHbJaHTCCrs8VO0A/S+Caa9A589GTqNxJ1zvhFRttyPgsmvEzqRoMIeT3v9EbbvAW9dBaVTQ6eROBvzb/jqLehzM2y1U+g0UkWFPY7y8vyFS4Wbw8u/gzWrQieSOPp+oh+zvsMhsPe5odPIOlTY46rh1nDMQ/C/STDk+tBpJG7KlvtGQ/0mcMyDfvlGyRgq7HHW4WDY5wIY0x+maHEridDbV8GCaXDcI7B5Ueg0sh4V9rjrfaNfuWbA+bB0bug0EgeTXoXPnoIDLvbnciTjqLDHXUFdOOEJKC+DV/tBZUXoRJLNFn0Lb1RdL9HzmtBpZANU2HNBUXs4/G6YOQJG3BM6jWSrijXwyu8BByc8pqGNGUyFPVfs1hd2ORGG/Q2+Gx06jWSjYXfA7DFw5L+gSdvQaWQjVNhzhRkcfg80buVHMyzTurNSC9OGwIh/QJfTofPxodNIDZIu7Gb2JzObamaTzezOKEJJitTbAk56Clb8AC//FirKQyeSbLDwG98F07wzHHZX6DSSgGTXPO0JHA3s6pzrBNwdSSpJnW12gyP+5fvb378ldBrJdGtWwotnAA5Ofsqv2iUZL9kW+3nAHc651QDOufnJR5KU69IXis+BkffCFwNCp5FM5RwMusRfYXrcv6Hp9qETSYKSLewdgG5m9omZfWhmXaMIJWlw6N/8kLXXz4fSr0KnkUxU8jiMfxa6XwUdDgmdRmqhxsJuZkPNbFI1t6OBAqAJsA9wOfCiWfXXFptZPzMrMbOS0lKduAuuoK7vby+oBy+c5hcgFllr1lh460po3we6Xxk6jdRSjYXdOdfbOde5mtsAYDbwqvPGAJVAtdcXO+f6O+eKnXPFzZo1i/YoZNM0agEnPgE/TIcBF2j+dvGWlcKLZ8IW28Jx/f2kcpJVkn3FXgcOAjCzDkAhsCDZUJJG2x3o52//YgB8fF/oNBJaRbkfMbVyIZz8NGzWNHQi2QQFSf7/x4HHzWwSUAac5ZyafVlnvz/DnE9h6E1QtCPseGjoRBKCc/DW5X7E1DEP+xFUkpWSKuzOuTLg9IiySChmforfRTP9xUvnvOMXx5bcMvohf8J0/4v8yCnJWuo8E69wc+j7AtRrBM+erMWwc83Ut+Cda2DnI6HXTaHTSJJU2OVnW2wDp74AKxf74l62PHQiSYd5E+Dlc2DbLnCsTpbGgV5B+aVtdoUTHofvJ1RN81sZOpGk0tJ5/o94/cbQ93ldWRoTKuzyazseCof8FaYMgvduCp1GUqVsOTx3Cqxe6j+pNdw6dCKJSLKjYiSu9j7XL3028l5o2g72PCt0IolSZaX/RPb9BN9S18nyWFFhl+qZwW/uhEXfwJuX+NacLiuPB+f8mqVTBsGhf9frGkPqipENyy+AE5/0rbkXz4SZI0Mnkih88FcY8wjseyHs/cfQaSQFVNhl4+ptAae9Ao3b+JNscz8PnUiSMeoBGH4n7H4GHHyb/2QmsaPCLjXbfEs44zWo3wT+e7xmg8xWnz3tx6p3PAaOvFdFPcZU2CUxjVrAma+D5cPTx8Di70InktqY/Dq88Wdo18vPrZ6XHzqRpJAKuyRuy3a+5V62DJ46GpZpXZWsMH2oX9qu5V5+Yq+CwtCJJMVU2KV2tu4Mp73spxx4+lhYsTB0ItmY70bDC2dAs538WPXCzUMnkjRQYZfaa7UXnPIsLPgKnjzKz98tmeebEfD0cX5e9TNe9VeXSk5QYZdN066nbwH+MB3+c5i/NF0yx7Sh8MwJ0LgVnP0mNNgqdCJJIxV22XTtDoLTX4Glc+GJ3+iEaqb4cpCfKqBoBzh7sKYKyEEWYl2M4uJiV1JSkvb9SrSKbxvCgmVldLHpPFl4B8uoz6ll17J889aUXNcndLycdM0tN3JLxX1MdNtzVtkVLKUBAEUNCvWaxICZfeqcK67peUm12M3sBTMbV3WbaWbjktmeZJcFy8oAGOfa07fsOupRxouFt9B4+YzAyXLU5//ltop7KXE7cnrZ1T8Vdfj5tZLckFRhd86d7Jzr4pzrArwCvBpNLMk2X7i2nFx2PQa8UHgrzNXf+LT6pD8MuICPKjtzdtkVLKd+6EQSUCR97GZmwEnAc1FsT7LTdNeSE8tuYCV1fZ/7lDdDR4q/ygp4+2q/VumOh/P7NZexirqhU0lgUZ087Qb8zzk3LaLtSZb61m3Nsatvga12hudPg5H3+dkEJXqrf4TnT4XRD8Le58HJT1NGndCpJAPUOG2vmQ0Fqjutfq1zbkDV132pobVuZv2AfgCtW7euZUzJJqU09kPsXjsXhlwPP0yDw++BfBWdyCyZDc+eAvO/gMP/AV1/HzqRZJAaC7tzrvfGvm9mBcBxwJ41bKc/0B/8qJhaZJQMVdSgsNqTckUNCqFOfTjhCfigPYy4GxbNhJOe8hOJSXLmfArP9YU1K+G0F6H9z2/Rjb4mkjOSHu5oZocCVzvnuif6fzTcMceMew4G/gmatIFTX/Rzzsimmfy6/yTUoJn/WW61c+hEkkZpGe5Y5RR00lQ2pktfOGugn1emfw+Y9EroRNmnfDW8dRW8dJZf+OT376uoywYlXdidc2c75x6OIozEWJv9oN8wPxnVy7+DARf6xZSlZgumw6O94ZOHYK8/wllv+Ba7yAZoSgFJnyZt4LeD4YBL4PP/Qv+e8L/JoVNltnHPwSMHwpJZcMpzcNidUKde6FSS4VTYJb3y60DvG/2iHasW++I+9lENiVzf6h/h1X7w+rmwbRc4dyTsdFjoVJIlVNgljO17+GK13YHw5qV+zPvSuaFTZYaZH/lW+sSXoMc1vuulUYvQqSSLqLBLOGtHdhx8O3z9Htzf1S+2XFEeOlkYy0r9iJf/HA6V5XDWIOhxpZaxk1pTYZew8vJgvwvh/NH+BOs71/iRM7PGhE6WPpWVUPI43F8ME1+GbpfC+Z9A2/1DJ5MspcIumaHpdr71ftLTsHIhPNYHBv45/kvvzZvgj3XQxX4Y43kjodcNULhZ6GSSxWq88lQkbcyg41F+daZhd8Doh+DLN2D/i/wl83Ub1LyNbLFwBoy4B8Y9A/WbwrGPwK4n+5+BSJK00IZkru8n+blmvn7fF7/9LoSuf4B6W4ROtul++BqG3w0TXoC8Aij+LfS4SlMtSEISvfJUhV0y36yxMPxOmPYu1GsM+14Ie/eDeo1CJ0vcgmkw/C4/0iW/EIp/5z+JaNk6qQUVdomfOZ/Bh3fCV2/5or7rydDlVNimS2Z2YZSX+T9G456FqYP9xGhdz4F9/wQNm4dOJ1lIhV3ia954GHmvX7S5YjU029nPR7PLSbDFNmGzOQfzxvkrRie+5E8Eb74V7H467HO+pgKQpKiwS/ytXASTX/NFdPYYsDzYvifsfCS0PQC2bJ+elnx5Gcz9HL4Z7ic4K/0S8uv6K0V3OxXaHQT5GqcgyVNhl9yyYDqMf86flFwyyz/WoLkv8G32h7bdoGiHaAp9eRnM/QxmjvBXic4aA2tW+O+17Oq7hzodqxOiEjkVdslNzvmhhGuL7syP4Md5/nuFDaBx659vjVr5fxs096399VWUwdI5sPg7WPxt1b/f+dWLKquujm3e+ec/Hm32h823TN+xSs5JtLDr86HEi5lfyGPLdrDn2esU+o+gdMrPRfq7UbBqSeLbbbiN/0PQohg6HQct9vRXym7WNGWHIrKpVNgl3tYt9Otbudh32yybX/3/zSuARi1hixaaKleySlKF3cy6AA8D9YBy4HznXA5N8iFZrX5jqN+Y4tuGbHCd0JLrtIyfZJ9k54q5E7jZOdcFuKHqvkhWqa6ob+xxkUyXbGF3wNrruxsBmlBbRCSwZPvY/wK8Y2Z34/9I7LehJ5pZP6AfQOtxPtMmAAAEh0lEQVTWrZPcrYiIbEiNhd3MhgLVTWhxLdALuNg594qZnQQ8BvSubjvOuf5Af/DDHTc5sYiIbFSNhd05V22hBjCzp4CLqu6+BDwaUS4REdlEyfaxzwW6V319EDAtye2JpF1Rg8JaPS6S6ZLtY/8DcK+ZFQCrqOpDF8kmJdf1CR1BJFJJFXbn3EfAnhFlERGRCGjNUxGRmFFhFxGJGRV2EZGYUWEXEYmZIPOxm1kp8O0m/vciYEGEcbKBjjk36JhzQzLH3MY5V+P6ikEKezLMrCSRiebjRMecG3TMuSEdx6yuGBGRmFFhFxGJmWws7P1DBwhAx5wbdMy5IeXHnHV97CIisnHZ2GIXEZGNyNjCbmaHmtlUM5tuZldV8/26ZvZC1fc/MbO26U8ZrQSO+RIz+8LMJpjZe2bWJkTOKNV0zOs87wQzc2aW1SMoEjleMzup6nWebGbPpjtj1BL4vW5tZh+Y2edVv9uHhcgZJTN73Mzmm9mkDXzfzOy+qp/JBDPbI9IAzrmMuwH5wNfA9kAhMB7ouN5zzgcervr6FOCF0LnTcMw9gc2qvj4vF4656nkNgeHAaKA4dO4Uv8Y7AJ8DTarubxU6dxqOuT9wXtXXHYGZoXNHcNwHAnsAkzbw/cOAtwAD9gE+iXL/mdpi3wuY7pyb4ZwrA54Hjl7vOUcDT1Z9/TLQy8wsjRmjVuMxO+c+cM6tqLo7GmiZ5oxRS+R1BrgVv1D6qnSGS4FEjvcPwAPOuUUAzrn5ac4YtUSOOXZrJzvnhgMLN/KUo4GnnDcaaGxm20S1/0wt7C2AWevcn131WLXPcc6VA0uALdOSLjUSOeZ1nYP/i5/NajxmM9sdaOWcG5TOYCmSyGvcAehgZiPNbLSZHZq2dKmRyDHfBJxuZrOBwcCf0hMtqNq+32sl2YU2UqW6lvf6w3cSeU42Sfh4zOx0oJifV6/KVhs9ZjPLA/4JnJ2uQCmWyGtcgO+O6YH/RDbCzDo75xanOFuqJHLMfYH/OOf+YWb7Ak9XHXNl6uMFk9L6lakt9tlAq3Xut+TXH89+ek7VCk6N2PhHn0yXyDFjZr3xC4kf5ZxbnaZsqVLTMTcEOgPDzGwmvi9yYBafQE3093qAc26Nc+4bYCq+0GerRI75HOBFAOfcKKAefj6VOEvo/b6pMrWwjwV2MLPtzKwQf3J04HrPGQicVfX1CcD7ruqsRJaq8ZiruiUewRf1bO97hRqO2Tm3xDlX5Jxr65xriz+vcJRzriRM3KQl8nv9Ov4kOWZWhO+amZHWlNFK5Ji/A3oBmNnO+MJemtaU6TcQOLNqdMw+wBLn3LzIth767PFGziofBnyFP6N+bdVjt+Df2OBf/JeA6cAYYPvQmdNwzEOB/wHjqm4DQ2dO9TGv99xhZPGomARfYwPuAb4AJgKnhM6chmPuCIzEj5gZBxwcOnMEx/wcMA9Yg2+dnwOcC5y7zuv8QNXPZGLUv9e68lREJGYytStGREQ2kQq7iEjMqLCLiMSMCruISMyosIuIxIwKu4hIzKiwi4jEjAq7iEjM/D+Qm6s2WoOEEQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x=np.linspace(0,L)\n",
    "w_an=-q*x*(L**3-2*x**2*L+x**3)/24/E/I\n",
    "\n",
    "plt.plot(xnum,np.block([0,w*1000,0]),'s')\n",
    "plt.plot(x,w_an*1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise \n",
    "\n",
    "Divide the simply-supported beam into 12 sections and plot the deflection as a function of distance along the beam with a uniform load. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Discussion\n",
    "\n",
    "What is the convergence rate for this central difference method? _Hint: check the magnitude of round-off error for [central difference methods](./01_Revisiting_derivatives.ipynb)._"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# What We've Learned\n",
    "\n",
    "* The difference between a _PDE_ and an _ODE_\n",
    "* How to approximate differential equations with boundary conditions\n",
    "* Solve steady-state heat transfer problem\n",
    "* Solve static deflection of elastic beam\n",
    "* Demonstrate convergence of finite difference solutions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problems\n",
    "\n",
    "![Thermal fin connected to a constant temperature heat sink](../images/thermal_connect.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Consider the thermal fin shown above connected to a heat sink with constant temperature. $h'=\\frac{2h}{\\kappa R}$ is the modified convective heat transfer for the fin. And our boundary conditions give us values for $T_{0}~and~T_{7}.$ We can plug in constants for forced air convection, $h=100~W/m^2K$, aluminum fin, $\\kappa=200~W/mK$, and 60-mm-long and 1-mm-radius fin, the air is room temperature, $T_{\\infty}=20^oC$, the base is $T_{base}=T_{0}=100^oC$, and the sink is $T_{sink}=25^oC$. Use the following finite difference equation to solve for the temperature along the fin and the heat flux through the fin given, \n",
    "\n",
    "$T(x=0)=100^oC,~and$\n",
    "\n",
    "$T(x=60~mm)=25^oC.$\n",
    "\n",
    "$\\frac{T_{i-1}-2T_i+T_{i+1}}{\\Delta x^2}+ h'(T_{\\infty}-T_i) = 0$\n",
    "\n",
    "a. Set up and solve the finite difference equations for $\\Delta x=10~mm$, plot the resulting temperature $T(x)$. \n",
    "\n",
    "b. Set up and solve the finite difference equations for $\\Delta x=5~mm$, plot the resulting temperature $T(x)$. \n",
    "\n",
    "c. Set up and solve the finite difference equations for $\\Delta x=1~mm$, plot the resulting temperature $T(x)$. \n",
    "\n",
    "d. Plot the heat flux through the fin, $-\\kappa \\frac{dT}{dx}$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. Consider the encastre beam shown in the __Static Beam deflections__ section. Use the following material and geometry (1-m steel rod 1-cm-by-1-cm) with 100 N/m load applied\n",
    "\n",
    "$EI \\frac{d^4w}{dx^4} = q.$\n",
    "\n",
    "We can approximate the function as a finite difference approximation as such,\n",
    "\n",
    "$\\frac{d^4w}{dx^4} \\approx \\frac{w(x_{i+2})−4w(x_{i+1})+6w(x_i)−4w(x_{i-1})+w(x_{i-2})}{h^4}=\\frac{q}{EI}.$\n",
    "\n",
    "* $L=1~m$\n",
    "* $E=200e9~Pa$\n",
    "* $I=\\frac{0.01^4}{12}~m^4$\n",
    "* $q=100~N/m$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "a. Solve for the four integration constants using the boundary conditions shown in the __Static Beam deflections__ section. $w(x)=\\frac{q_0x^4}{24}+\\frac{Ax^3}{6}+\\frac{Bx^2}{2}+Cx+D$\n",
    "\n",
    "b. Create a finite difference approximation with 10, 20, 30, and 40 segments. \n",
    "\n",
    "c. Plot the error between the maximum predicted numerical deflection (b) and the analytical deflection (a). What is the convergence rate of the finite difference approximation?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}