{ "cells": [ { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## BVPs, BCs, nonlinear, nonuniform grids\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Boundary conditions\n", "* Dirichlet conditions are straightforward\n", " * Unknowns are interior points. Points near the boundaries are written in terms of the boundaries.\n", "* Neumann and Robin conditions:\n", " * $y^{\\prime} = \\alpha$ is given at the boundary (or $y^{\\prime} + \\beta y = \\alpha$ for Robin).\n", " * Interior cells need $y$ on the boundary, not $y^{\\prime}$.\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "#### Two Approaches\n", "**Approach 1. Ghost Cell Method** Include the boundary point in the list of unknowns (like interior points).\n", "* Discretize the boundary point exactly as if it was an interior point.\n", "* For central differences, this will reference one point past the boundary.\n", " ``` \n", " \\\\|\n", " \\\\|\n", " \\\\|\n", " * * * *\n", " i-1 \\\\|i i+1 \n", " -1 \\\\|0 1 2\n", " \\\\|\n", " \\\\|\n", " ```\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "* Use the BC $y^{\\prime}=\\alpha$ as an equation for the new unknown $y_{-1}$.\n", " $$y_0^{\\prime} = \\alpha \\rightarrow \\frac{y_1-y_{-1}}{2\\Delta x} = \\alpha \\rightarrow$$\n", " \n", " $$y_{-1}=y_{1} - 2\\Delta x\\alpha.$$\n", " \n", "* Substitute this equation in for $y_{i-1} = y_{-1}$ appearing in the finite difference equation (FDE) at point $i=0$. The first two equations below are the FDE at $i=0$. The last equation has the substitution for $y_{-1}$:\n", " $$l_iy_{i-1} + a_iy_i + u_{i}y_{i+1} = F_i$$\n", " $$\\mbox{or,}$$\n", " $$l_0y_{-1} + a_0y_0 + u_0y_{1} = F_0, $$\n", " $$\\rightarrow a_0y_0 + (l_0+u_0)y_1 = F_0 + 2\\Delta x\\alpha l_0. $$\n", "* That is, the unknown $y_{-1}$ is in terms of $y_1$, so when we substitute it into the FDE at $i=0$, the coefficient of $y_1$ and the RHS are modified in the $i=0$ equation. \n", "* This is called the **Ghost Cell Method** since a false or *ghost* point arises, which is then handled with the boundary condition equation.\n", " \n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "\n", "* **Advantages**\n", " * Solves directly for the unknown boundary value.\n", " * Uses a uniform stencil (central difference everywhere, even at the boundaries.\n", "\n", "* **Disadvantage**\n", " * Higher order $y^{\\prime\\prime\\prime}$ or higher order FDA can lead to multiple outside points, which can be awkward.\n", " * Remedy this using one-sided differences." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "**Method 2** Don't include the boundary point in the list of unkowns. \n", "\n", " ``` \n", " \\\\|\n", " \\\\|\n", " \\\\|\n", " * * * * \n", " \\\\|i-1 i i+1\n", " \\\\|-1 0 1 2 \n", " \\\\|\n", " ```\n", "* The first unknown point is $i=0$.\n", "* The FDE at this point is in terms of point $i-1=-1$\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "* Use the Boundary condition at $i-1$ with a one-sided difference to get another equation to write $y_{-1}$ in terms of the interior points we are solving for ($y_0$, $y_{1}$, etc.)\n", "* At the left side for the $i=0$ FDE (which has a $y_{-1}$ term):\n", " * Write a one sided difference for $y^{\\prime}_{-1}$.\n", " $$y^{\\prime}_{-1} = \\alpha = \\frac{-\\frac{3}{2}y_{-1} + 2y_0 - \\frac{1}{2}y_1}{\\Delta x}.$$\n", " * Solve for $y_{-1}$:\n", " \n", " $$ y_{-1} = \\frac{\\alpha\\Delta x - 2y_0 + \\frac{1}{2}y_{1}}{-3/2}.$$\n", " \n", " * The FDE at point $i=0$ is \n", " $$l_0y_{-1} + a_0y_0 + u_0y_1 = F_0.$$\n", " * Insert the above BC equation for $y_{-1}$ (blue) into this FDE to get the final FDE at the $i=0$ point.\n", " $$ (a_0 + \\frac{4}{3}l_0)y_0 + (u_0 - \\frac{1}{3}l_0)y_1 = F_0 + \\frac{2}{3}\\alpha l_0\\Delta x.$$\n", "* A similar procedure is done at the right side of the domain.\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "* **Advantage**\n", " * Extends easily to higher order\n", " * On the homework, if we exclude the BC point, then we don't have to divide by $r=0$ at the cylindrical centerline.\n", "* **Disadvantage**\n", " * Two different stencils are needed. One for the fully interior points $i=1,\\,2,\\ldots$, and one for the point $i=0$ next to the boundary. (Similarly for the right side of the domain.)\n", " * The boundary point remains unknown.\n", " * Once the solution to the interior points $i=0,\\,1,\\ldots$ is found, we can find $y_{-1}$, (the boundary point), using the above blue equation." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Nonlinear relaxation methods\n", "* The usual approach is to **iterate**\n", " * Use a linearized form of the equation.\n", " * Guess a solution.\n", " * Iterate to improve it.\n", "* Apply the FDA to $y^{\\prime\\prime}$, $y^{\\prime}$ as before.\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "\n", "* But consier a term like $y^\\prime y$\n", " - This results in nonlinear terms:\n", " $$y^{\\prime}y \\rightarrow \\left(\\frac{y_{i+1}-y_{i-1}}{2\\Delta x}\\right)\\cdot y_i.$$\n", "* We can linearize these terms by splitting the product and *lagging* part of it. For example\n", " $$y_i^2\\rightarrow y_i\\cdot y_i \\rightarrow y_i^{new}\\cdot y_i^{old},$$\n", " where $y_i^{old}$ is the value from the previous iteration, which is known.\n", "* For $$y^{\\prime\\prime} + P(x,y)y^{\\prime} + Q(x,y)y = F(x),$$\n", "use $$y^{\\prime\\prime} + P(x,y^{old})y^{\\prime} + Q(x,y^{old})y = F(x),$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "\n", "* Example\n", " $$\\nabla\\cdot\\vec{q} = F \\rightarrow -\\nabla\\cdot(k\\nabla T) = F.$$\n", " $$-\\frac{d}{dx}\\left(k\\frac{dT}{dx}\\right) = F.$$\n", " $$\\frac{dT}{dx}\\frac{dk}{dx} + k\\frac{d^2T}{dx^2} = -F.$$\n", " * $k = k(T)$ so the above equation is nonlinear. So lag $T$ when evaluating $k$. Use $k(T^{old})$.\n", " * This is now linear, but we have to iterate.\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "**Question** What are other approaches that could (or should?) be used?\n", "- Newton's method?\n", "- What about a Taylor Series linearization instead of the linearization shown?\n", " - Consider $yy$\n", " - above linearization: $yy\\approx y_0y$\n", " - Taylor series: $yy\\approx y_0^2 + 2y_0(y-y_0)$\n", " - Plot these for some arbitrary $y_0$" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEOCAYAAACn00H/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABDuUlEQVR4nO3dd1xUV/r48c9B6YKiCNi7YsPeNRqNPZaosTBo3GRTNn1TNSbZ1N3NxjXlu3E3yS/FyIiNWGPUGDX2ngiKvTeaCgrSOb8/LiDIoIIwBZ736zUv4J47c5+xzMO5zylKa40QQghRXE62DkAIIYRjkgQihBCiRCSBCCGEKBFJIEIIIUpEEogQQogSqWzrAKzJ19dXN2zY0NZhCCGEw9i7d2+81rqmpbYKlUAaNmzInj17bB2GEEI4DKXUmaLa5BaWEEKIEpEEIoQQokQkgQghhCgRSSBCCCFKRBKIEEKUYykpKWRnZ5fJa1s1gSil7lNKLVdKXVBKaaXU1FvalVLqHaXURaVUilJqo1Kq9S3n+Cil5iqlEnMec5VS1az5PoQQwhFkZmYSFhbGvHnzSE1NLfXXt3YPpApwAHgBSLHQ/hrwMvAc0AWIBX5RSnnlO2ce0BEYkvPoCMwtw5iFEMLhaK1Zvnw5586d48SJE3zzzTdcvXq1VK9h1QSitV6ltX5Da70YKNCnUkop4EXgn1rrcK31AeARwAsIzjmnJUbSeEJrvV1rvR14EnhQKdXCim9FCCHs2ubNm4mMjMz7OT4+nq+//pqzZ8+W2jXsqQbSCAgA1uYe0FqnAJuAnjmHegBJwLZ8z9sKJOc7pwCl1BNKqT1KqT1xcXFlEbcQQtiVgwcPsmHDhkLH09LSKM09oOwpgQTkfI255XhMvrYAIE7n+xPI+T423zkFaK2/0lp31lp3rlnT4mx8IYQoNy5cuMDSpUstto0YMYIGDRqU2rXsKYEIIYS4B4mJiYSFhZGZmVmorXfv3rRv375Ur2dPCSQ656v/Lcf987VFAzVz6iVAXu3EL985QghR4aSlpTFv3jySk5MLtbVs2ZL+/fuX+jXtKYGcwkgCA3MPKKXcgD7crHlsxxjJ1SPf83oAnhSsiwghRIWRnZ3N4sWLiY2NLdRWq1YtHnroIfL93l1qrD0PpIpSqr1Sqn3Otevn/Fw/p5bxKfC6UmqMUqoN8D1G0XwegNb6ELAa+FIp1UMp1QP4EliptT5izfciDA0bNmTmzJm2DkOICktrzerVqzl+/HihNi8vLyZNmoSzs3OZXNvaPZDOwO85D3fg3Zzv38tp/xfwCfAFsAeoBQzSWl/P9xrBwH5gTc5jPzDZGsELIYS92blzJ7t37y503NnZmeDgYLy8vCw8q3So0hzSZe86d+6s72Y/kHfffdcK0dze3/72t7s6T2vNxx9/zJdffsnFixdp2rQpr7/+OiEhITz22GPs2LGDPXv24O7uTlZWFv369aNq1aqsXLkSgGnTprFkyRLOnj2Lv78/48eP57333sPNzS3vGqtWreLdd98lIiICDw8PevbsyaJFixgyZAi//fZboXiEENZx+PBhFixYYLFt4sSJtGhx79PjlFJ7tdadLbVVqA2lyqM333yTxYsX88UXX9CiRQu2b9/O448/jo+PD59//jkdOnTglVde4YsvvuDDDz/k2LFjRERE5D3f09OTb7/9ljp16hAVFcVTTz2Fq6sr77//PgCrV69m5MiRTJs2je+++47MzEzWrl1LdnY2P/74I+3atePRRx/lL3/5i63+CISokC5evMiPP/5osW3w4MGlkjzuRBKIA0tOTmbWrFmsXbuWPn36ANCoUSN27drFF198wfDhwzGbzfTq1YsaNWrwj3/8g+XLl+Pn55f3Gm+99Vbe9w0bNuSNN95g5syZeQnk/fffZ9y4cXzwwQd55wUFBQHg4eFBpUqV8PLyIiDA4jQcIUQZSEhIICwsjIyMjEJtnTt3plu3blaJQxKIA4uKiiI1NZUhQ4YUGGGRkZFB7t7vXbp0YcaMGbzzzjs8/fTTDB06tMBrLF68mE8//ZTjx4+TlJREVlYWWVlZee2///47U6dOtcbbEULchdTUVObNm0dSUlKhtmbNmjF06NAyGXFliSQQB5a7RPOKFSuoX79+gbbcURdaa7Zs2UKlSpU4ceIEWuu8f1w7duxg4sSJ/O1vf+OTTz6hWrVqLF++nFdeecW6b0QIcVeysrJYuHAhlpZl8vf3Z+zYsTg5WW9slCQQC+62gG1rrVq1wtXVlTNnzhQ5SWjWrFns27ePTZs2MWzYMP7v//6P559/HoCtW7dSp06dArexzpw5U+D5HTp04Ndff+Xxxx+3+PouLi4FeixCiLKhtWbFihWcOnWqUJuXlxfBwcG4urpaNSZJIA7My8uLV155hVdeeQWtNffddx9JSUns2LEDJycnunXrxowZM5g3bx49e/Zk9uzZPPbYYwwYMIDWrVvTvHlzLly4gNlspkePHqxZs4awsLAC15gxYwYjRoygadOmBAcHo7Vm7dq1PPnkk3h4eNCwYUM2b95MSEgIrq6u+Pr62uhPQ4jybePGjezfv7/QcRcXF4KDg/H29rZ6TDKM18FprfnPf/7Df//7X06cOIG3tzft27fn+eefZ/r06XTt2pVvv/027/zJkycTERHBrl27cHV1Zfr06fy///f/SElJYdCgQQwcOJCnn366wHDc5cuX8+6773LgwAG8vLzo2bMnCxcuxM3NjR07dvDkk09y5MiRUl/pUwhh2LdvHytWrCh0XClFcHAwTZs2LbNr324YryQQIYSwY8ePH2fevHkWfzkbMWIEHTt2LNPr3y6B2NNaWEIIIfK5dOkSixYtspg8evfuXebJ404kgQghhB26evUqZrOZ9PT0Qm1BQUFlsrpucUkCEUIIO3Pjxg3MZrPFpdkbNWrEyJEjrTbX43YkgQghhB3JyMhg/vz5XL58uVCbn58f48ePp1KlSjaIrDBJIEIIYSeys7MJDw/n3Llzhdq8vb0xmUwFFjq1NUkgQghhB7TW/PTTTxw5UnhrI1dXV0wmk03metyOJBAhhLADv/32G/v27St0vFKlSkycOLHAIqj2QhKIEELY2N69ewvtrZProYceylsc1d5IAhFCCBs6dOgQP/30k8W2IUOG0Lp1aytHdPckgQiL3nnnHdq0aWPrMIQo106fPk14eLjFiYK9evWy2r4eJSUJxEEppW77kD08hLBv0dHRzJ8/3+Jq1u3atWPAgAE2iKp4ZDVeB3Xp0qW871euXMnjjz9e4Ji7u7stwiogIyMjb18SIcRNV65cITQ0lLS0tEJtzZo1Y8SIEXYxUfBOpAfioAICAvIe1apVK3AsOTmZKVOmEBAQgKenJx07dmTlypV5z33vvfcs3p7q1atX3l4ht8rOzub999+nXr16uLq60rZtW5YtW5bXfvr0aZRShIWF0b9/f9zd3fnyyy9L900LUQ5cv36d0NBQi7PM69aty7hx4+xmouCdSA+kCC+++CJ//PGHVa/Zvn17Pv3003t+naSkJIYOHcoHH3yAu7s7CxYsYMyYMURERBAYGMijjz7Ke++9x65du+jatSsAR44cYdu2bcyePdvia3722Wd8/PHH/O9//6Nz586EhoYyZswY9u7dS/v27fPOmz59OjNnzuSbb76R3ocQt0hNTcVsNnP16tVCbb6+vkyaNAkXFxcbRFYy0gMph9q1a8dTTz1F27Ztadq0KTNmzKBjx44sXrwYMH7LGTJkSIF9Qr799ls6depEu3btLL7mzJkzeeWVVwgODqZ58+a899579OnTh5kzZxY477nnnmPcuHE0atSIunXrlt2bFMLBZGRkEBYWRkxMTKE2b29vQkJC8PDwsEFkJSc9kCKURk/AVpKTk3n33XdZuXIlly5dIiMjg9TUVIKCgvLOefzxx3nkkUf45JNPcHFxYe7cuQW2ts3v2rVrXLx4kV69ehU43rt3b1atWlXgWOfOFrcNEKJCy8rKYtGiRZw9e7ZQm7u7O5MnT6Zq1ao2iOzeSAIph1555RVWr17NzJkzadasGR4eHkyZMqXAstDDhw/Hw8OD8PBwqlatSkJCAsHBwcW+1q2FPk9Pz3uOX4jyRGvNsmXLOHbsWKE2Z2dnTCaTw24FbVe3sJRSlZRS7yulTimlUnO+fqCUqpzvHKWUekcpdVEplaKU2qiUst+ZNjawZcsWpkyZwtixYwkKCqJu3bqcOHGiwDmVK1dm6tSpfPvtt3z77beMGTOmyN+AvL29qV27Nlu3bi10nVatWpXZ+xDC0WmtWbVqFZGRkYXanJycmDBhAnXq1LFBZKXD3nogrwPPAI8AkUAQMAdIA97POec14GVgKnAEeBv4RSnVQmt93doB26PmzZuzZMkSRo0ahbOzM++++y6pqamFzvvzn//MRx99hJOTE2vXrr3ta7766qu8/fbbNGvWjE6dOhEaGsrmzZstrt0jhDCsX78eS9toK6UYO3YsTZo0sUFUpcfeEkhPYIXWOnf3+NNKqeVANzB6H8CLwD+11uE5xx4BYoFgQMaNArNmzeKxxx6jT58++Pj48OKLL1pMII0bN6Zv376cOXOGfv363fY1n3/+ea5fv85rr71GTEwMLVq0IDw8vMiiuxAV3datW9myZYvFtgcffLBc9N6VpSn0tqKUmgY8DQzSWh9WSrUC1gD/0FrPVko1Bk4AXbXWu/M97ycgXmv9yO1ev3PnztrSbwMVWatWrTCZTMyYMcPWoQhRbuzZs6fI9a0eeOCBQgNS7JlSaq/W2uLoGHvrgXwEeAFRSqksjPg+1FrnTk4IyPl66zi4GMDijUSl1BPAEwD169cv9YAdVVxcHIsXL+b06dM8+eSTtg5HiHIjIiKiyOTRu3dvh0oed2JvCWQCMAXjdtRBoD3wmVLqlNb6m5K8oNb6K+ArMHogpRSnw/Pz88PX15cvv/zSYUeACGFvDh06xNKlSy22de7cmf79+1s3oDJmbwnkY2Cm1np+zs+RSqkGwHTgGyA657g/kH9AtX++NnEX7OnWpRDlwYkTJ4pcWTcoKIhhw4Y5xPpWxWFXw3gBD+DWpSmzuBnnKYxEMTC3USnlBvQBtlkjQCGEuNXp06eLXFk3MDCQkSNHlrvkAfbXA1kBTFNKncK4hdUBeAn4AUBrrZVSnwJvKKUOA0eBN4EkYJ5NIhZCVGjnz58nLCyMzMzMQm2NGzdm7NixDrM4YnHZWwJ5DmO+x2zAD7gEfA28l++cfwHuwBeAD7ATY9SWzAERQljVpUuXCA0NLbDKQ6769eszYcIEKle2t4/Z0mNXw3jLmgzjFUKUltjYWL7//ntSUlIKtdWqVYspU6bg5uZmg8hK1+2G8dpbDUQIIexefHw8P/zwg8Xk4efnR0hISLlIHnciCUQIIYrhypUrzJkzx+KGUDVq1GDy5MkOtyx7SUkCKac2btyIUor4+Pgyv1a/fv149tlny/w6Qtja1atXmTNnDklJSYXafHx8mDJlClWqVLFBZLYhCaSc6tmzJ5cuXaJGjRq2DkWIciEhIYE5c+Zw7dq1Qm1Vq1ZlypQpeHt72yAy2ym/wwMqOBcXFwICAu58ohDijhITE5kzZw6JiYmF2ry8vJgyZQrVqlWzfmA2Jj0QB7dp0ya6d+9OlSpVqFq1Kl27duXAgQOFbmF9//33VKlShZ9//pnAwEA8PDwYOXIkiYmJLF68mGbNmlG1alUmT55coDDYr18/nnrqKV544QV8fHzw8fHh1VdfJTs7u8iY0tPTef3116lbty4eHh506dKFNWvWlPmfhRBl4dq1a8yZM4eEhIRCbZ6enkyZMoXq1atbPzA7ID2Q27C0wvn48fD003DjBgwbVrh96lTjER8P48YVbv/LX2DCBDh3DiZPLti2cWPx4svMzGTUqFE89thjmM1mMjIy2LdvX5GTltLS0vj3v/+N2WwmPT2dsWPHMnbsWNzd3QkPD+fy5cuMGTOG2bNn8/LLL+c9z2w2M3XqVLZv305ERASPP/44tWrV4qWXXrJ4nT/96U+cOHGCefPmUbduXVatWsWIESPYvXu3LP8uHEpu8rh69WqhNk9PTx555JEKvZacJBAHdu3aNRISEhgxYkTexjSBgYEAxMTcumCxkXC++OILWrRoAUBwcDCffPIJMTExef8JRo0axYYNGwokkFq1avH555+jlCIwMJCjR48ya9YsiwnkxIkThIWFcfr06bzVj5999lnWrVvHl19+yezZsws9Rwh7lJs8rly5Uqgtd5vomjVr2iAy+yEJ5DZu1yPw8Lh9u6/v7dvr1St+j+NW1atXZ+rUqQwePJgBAwYwYMAAxo0bV+Sy9a6urnnJA8Df35+AgIACv0H5+/sTFRVV4Hndu3cvsI5Pjx49eOutt7h27VqhouG+ffvQWhfaLCctLa3crUQqyq9r167xww8/WEwe7u7uTJ48GT8/PxtEZl8kgTi47777jhdffJHVq1ezfPlyZsyYwdKlS3F1dS107q1LKiilcHZ2LnTsdvWNO8nOzkYpxe7duwu9tru7e4lfVwhruV3Pw83NjcmTJ8sAlRySQMqBdu3a0a5dO15//XWGDh3KnDlzeOKJJ0rt9Xfu3InWOq8XsmPHDmrXrm1xyGKHDh3QWhMdHc39999fajEIYQ13Sh5TpkyhVq1aNojMPskoLAd26tQppk2bxrZt2zhz5gwbNmwgIiKi1PdavnjxIi+++CJHjhxh8eLFfPzxx/z1r3+1eG7z5s0xmUxMnTqVxYsXc/LkSfbs2cPMmTP58ccfSzUuIUpTYmIi33//vSSPYpAeiAPz8PDg6NGjPPzww8THx+Pv74/JZOL1119n69atpXYdk8lEVlYW3bp1QynFY489VmQCAeO22ocffshrr73G+fPnqV69Ol27dpUeibBbuZMELQ3VleRRNFmNV9xWv379aNOmDf/5z39sHYoQZSJ3eRJLkwRzC+YVOXncbjVe6YEIISqsy5cv88MPP1hcnsTd3Z0pU6ZIwfw2JIEIISqkuLg4fvjhB4sLI+bO8/D397dBZI5DEoi4rY33OllFCDsUHR3N3LlzuXHjRqG23OVJZJ7HnUkCEUJUKBcuXCA0NJTU1NRCbVWqVKnwy5MUhyQQIUSFcfbs2by14G7l7e3NlClTZAuEYpAEIoSoEE6ePMn8+fPJyMgo1FatWjWmTJmCj4+PDSJzXJJAhBDl3tGjR1m4cCFZWVmF2qpXr86UKVOoWrWqDSJzbJJAhBDl2oEDB1iyZInFNd5q1qzJ5MmT8fLyskFkjk8SiBCi3Nq7dy8rV6602BYQEMDkyZPx8PCwclTlhyQQIUS5tH37dtauXWuxrW7duphMJtzc3KwcVfkiCUQIUa5ordmwYQObN2+22N6oUSMmTpyIi4uLlSMrf+xuNV6lVC2l1BylVJxSKlUpFaWU6puvXSml3lFKXVRKpSilNiqlWtsyZiGEfdBas2rVqiKTR/PmzQkODpbkUUrsKoEopaoBWwEFDAdaAs8BsflOew14Oed4l5y2X5RSUgUTogLLyspiyZIlFLVgaps2bRg/fnyhjdVEydnbn+RrwCWt9ZR8x07lfqOMHY1eBP6ptQ7POfYIRhIJBr60XqhCCHuRnp7OokWLOH78uMX2Tp06MWzYMJyc7Op3Zodnb3+ao4GdSqkFSqlYpdQfSqln1c0NuRsBAUBeZUxrnQJsAnpaPVohhM2lpKQwd+7cIpNH7969GT58uCSPMmBvf6KNgaeBk8Bg4DPgn8AzOe256yrH3PK8mHxtBSilnlBK7VFK7YmLiyv9iIUQNnPt2jW+++47zp8/b7H9gQceYMCAAdz8HVSUJnu7heUE7NFaT8/5+XelVDOMBFKiHY201l8BX4GxoVSpRCmEsLm4uDhCQ0Mt7uWhlOLBBx+kY8eONois4rC3HsglIOqWY4eA+jnfR+d8vXWRfv98bUKIcu78+fN89913FpNHpUqVePjhhyV5WIG9JZCtQItbjjUHzuR8fwojUQzMbVRKuQF9gG3WCFAIYVtHjx5lzpw5pKSkFGpzcXHBZDLRsmVLG0RW8djbLaxPgG1KqRnAAqAD8DzwBoDWWiulPgXeUEodBo4CbwJJwDybRCyEsJrff/+dFStWoHXhu9Genp6YTKYKvX+5tdlVAtFa71ZKjQb+DrwFnM35Ojvfaf8C3IEvAB9gJzBIa33dutEKIaxFa83mzZvZsGGDxXYfHx9CQkKoXr26lSOr2OwqgQBorX8CfrpNuwbeyXkIIcq57OxsfvrpJ/bt22exPSAgAJPJRJUqVawcmbC7BCKEELnS09MJDw/n6NGjFtsbNWrEhAkTcHV1tXJkAiSBCCHsVFJSEmFhYVy8eNFie5s2bRg9ejSVKlWycmQilyQQIYTdiY+Px2w2k5CQYLG9R48eDBw4UCYI2pgkECGEXTlz5gzz588nNTXVYvuQIUPo1q2blaMSlkgCEULYjcjISJYtW2Zx7/JKlSoxZswYWrVqZYPIhCWSQIQQNqe1ZtOmTWzcuNFiu7u7OxMnTqR+/foW24VtSAIRQthUZmYmK1euZP/+/RbbfXx8MJlM1KhRw8qRiTuRBCKEsJkbN26wYMECzp49a7G9bt26TJw4EU9PTytHJu6GJBAhhE3Ex8cTFhbGlStXLLa3atWK0aNH4+zsbOXIyo/z588TFhbGqVOnmD179p2fUEySQIQQVnfy5EkWLVpU5EirXr16yT4eJZSQkEB4eDhms5mNGzeitaZXr16kp6eX+l7wkkCEEFa1Z88eVq1aZXFBRCcnJ4YPHy5LsRdTWloaq1atwmw2s3LlStLS0mjWrBnvvPMOwcHBNG3atEyuKwlECGEV2dnZrFmzhl27dllsd3NzY/z48TRq1MjKkTmm7OxstmzZQmhoKIsWLSIhIYGaNWvy5JNPEhISQufOncu8BycJRAhR5lJSUli8eDEnT5602O7j40NwcDC+vr5WjszxHDx4kNDQUObNm8fZs2fx8PDgoYceIiQkhAceeIDKla33sX7XV1JKVdZaZ5ZlMEKI8udOxfKGDRvy8MMP4+HhYeXIHMeFCxcICwsjNDSU/fv3U6lSJQYNGsTf//53Ro0aZbOViIuTqi4ppeYA32itD5VVQEKI8uPYsWOEh4eTlpZmsb1Dhw4MHz5cFkS0IDExkfDwcEJDQ/OK4V27duXzzz9nwoQJ+Pn52TrEYiWQN4A/AX9VSu0C/h+wQGudVCaRCSEcltaarVu38uuvv1psV0oxcOBAunfvLiOt8klPT+fnn38mNDSUFStWkJaWRtOmTXn77bcJDg6mefPmtg6xgLtOIFrrr4GvlVItgUeBD4BPlVKLMHolW8soRiGEA8nIyGD58uUcOHDAYrurqyvjxo0rs5FBjiY7O5utW7diNptZuHAhV69epWbNmjzxxBOYTCa6du1qt0m22NWWnNtXryqlpgFPAx8DjyiljgGfAl9prbNLNUohhEO4evUqCxYsICYmxmJ79erVmTRpkhTLgaioqLxi+JkzZ/Dw8GD06NF5xXBHmEBZ7ASilHIBxmD0QvoDW4BvgNoY+5f3AyaWXohCCEdw8uRJFi9eTEpKisX2Jk2aMHbsWNzd3a0cmf24ePEiYWFhmM1mfv/9d5ycnBg4cCAffvihTYvhJVWcUVgdMZLGJCAD+AF4Vmt9NN85K4E9pR2kEMJ+aa3Ztm0bv/76q8XJgQA9e/ZkwIABODk5WTk627t27Ro//vgjoaGhrF+/Hq01Xbp04bPPPmPChAn4+/vbOsQSK04PZDewFngCWFbEkN7TwPxSiEsI4QDS0tJYvnw5UVFRFtsrV67MiBEjCAoKsnJktpWens7q1avziuGpqak0btyYt956C5PJZHfF8JIqTgJZiTHy6qeiahxa62SMkVpCiHIuPj6ehQsXEhcXZ7G9atWqTJgwgVq1alk5MtvI7YmZzWYWLFjAlStX8PX15c9//jMmk4lu3brZbTG8pIqTQK5j9C4SlVLfA99prY+VSVRCCLsWFRXFsmXLSE9Pt9jesGFDxo0bVyGWYT906BBmsxmz2czp06dxd3fPK4YPHDjQIYrhJVWcYbwhSilvwITRy5imlNqC0StZpLW2XDkTQpQbWVlZ/Prrr2zfvr3Ic3r06MEDDzxQrusdly5dypsZnr8Y/t577zF69Gi8vLxsHaJVFGsUltb6GvBf4L9KqdbAn4Evgc+VUguAT2WWuhDl0/Xr11m8eHGRmz85OzszcuRI2rRpY+XIrOP69esFiuHZ2dl07tyZTz75hIkTJxIQEGDrEK2uRKtuKaVqA6OAB4FMIByoB0QopaZrrWeWRnBKqenA34EvtNbP5hxTwN8wivk+wE7gGa31wdK4phCisJMnT/Ljjz+SnJxssb169ep2s7xGaUpPT2fNmjWYzWaWLVuWVwyfMWMGJpOJFi1a2DpEmyrOMF5njKTxKDAQ+B34FxCWu5yJUmokxvDee04gSqnuGEki4pam14CXganAEeBt4BelVAut9fV7va4Q4qbs7Gw2b97Mb7/9VuQQ3ZYtWzJq1ChcXV2tHF3Z0Fqzffv2vGL45cuXqVGjBo899hgmk0mWX8mnWIspAgqYB0zTWt/6wQ6wCbh6r0EppaoCZoxk9bd8xxXwIvBPrXV4zrFHgFggGON2mhCiFCQnJ/Pjjz8WuQS7UooBAwbQs2fPcvGBevjwYcxmM/PmzePkyZO4u7szatQoTCYTgwcPLtfF8JIqTgL5K0ax3PIelIDWOgEojd1gvgIWa603KKX+lu94IyAAYz5K7jVTlFKbgJ5IAhGiVJw+fZrw8HCSkiyvlVqlShXGjRtHgwYNrBxZ6YqOjmb+/PmEhoayd+9enJycGDBgAG+//TZjxoypMMXwkirOKKy5ZRlILqXU40BTIMRCc26V6taFdmKAOkW83hMYt8KoX79+KUUpRPmUnZ3Npk2b2LRpU5G3rBo2bMjYsWMdbtmNXNevX2fJkiWYzWbWrVtHdnY2nTp1YtasWUycOLHCzFspDXa1I6FSqgVG0by31jqjNF5Ta/0VRo+Gzp07W/4fIYTIG2V0+vTpIs/p06cP/fr1c7ghuhkZGaxduxaz2czSpUtJSUmhYcOGvPHGG5hMJgIDA20dokOyqwQC9AB8gYP57qlWAu5TSj0FtM455g/kH0voD0RbK0ghypujR4+ybNkybty4YbE9d9tUR1qCXWvNzp07CQ0NZcGCBcTHx1O9enWmTp1KSEgIPXr0KBe1G1uytwSylMKLMX4HHMPomRzFSBQDMdbmQinlBvQBXrValEKUE5mZmaxbt46dO3cWeU79+vUZO3Ys3t7eVoys5I4ePZo3M/zEiRO4ubkxcuRIQkJCGDx4MC4uLrYOsdywqwSSU4RPyH9MKZUMXNFaH8j5+VPgDaXUYYyE8iaQhDE6TAhxl+Li4ggPDy9y7w6A++67j759+9r9LauYmBjmz5+P2Wxm9+7deSPE3nzzTcaMGeMwyc/R2FUCuUv/AtyBL7g5kXCQzAER4u5ordmzZw9r164lM9PSotrGKKsxY8bQqFFpDKosG0lJSSxdupTQ0FB++eUXsrOz6dChA//+97+ZOHEitWvXtnWI5Z7dJxCtdb9bftbAOzkPIUQxJCcns2LFCo4cOVLkOc2aNWPUqFF2uRBi/mJ4bs2mQYMGTJs2jZCQEFq2bGnrECsUu08gQojScezYMZYtW1bkciS5CwLa27LjWmt27dpFaGgo8+fPJz4+Hh8fHyZPnozJZKJXr152f4utvJIEIkQ5l56ezi+//MKePUVvFurr68uYMWPsag6EpWL4iBEjMJlMDB06VIrhdkASiBDl2Pnz51myZAlXrlwp8pxOnTrZzVIdMTExLFiwgNDQ0LxieP/+/ZkxYwZjxoyhatWqtg5R5CMJRIhyKCsri02bNrF58+YiZ5S7u7szYsQIm9cNkpOTWbp0KWazmbVr15KVlUX79u35+OOPmTRpEnXqWFxkQtgBSSBClDMxMTEsXbqU6Oii59Y2adKEUaNG2Wytp9z5J6GhoSxZsiSvGP7aa69hMplo3br1nV9E2JwkECHKiezsbLZt28bGjRvJysqyeE7lypV54IEH6Nq1q9UL5Vprdu/ejdlsZv78+cTGxuLj40NISAghISFSDHdAkkCEKAfi4uJYtmwZFy5cKPKc2rVr89BDD+Hr62vFyOD48eN5xfBjx47h6upaoBheXvYRqYgkgQjhwO6m16GU4r777qNPnz5UqlTJKnHFxsayYMECzGYzO3fuRClFv379mDZtGmPGjKFatWpWiUOULUkgQjiomJgYli1bxqVLl4o8p2bNmowePdoqs7KTk5NZtmwZZrOZNWvWkJWVRVBQEB9//DETJ06kbt26ZR6DsC5JIEI4mKysLDZv3szmzZvJzs4u8rwePXrQv39/Klcuu//mucVws9nMkiVLSE5Opl69erz66quYTCbatGlTZtcWticJRAgHcv78eZYvX05cXFyR59SoUYNRo0ZRr169Mokhdy0ts9lMWFgYsbGxVKtWjeDgYEJCQujdu7cUw21Ma7h0CdLSoCyXM5MEIoQDSE9PZ/369bdddh2ge/fu9O/fv0wmBZ44cYJ58+YRGhrK0aNHcXFxySuGDxs2TIrhNjZ3LuzdCxERxuPyZZg0CeaV4TrlkkCEsHNHjhxh1apVXLt2rchzfH19GTVqVKnXGeLi4li4cCGhoaHs2LEDpRR9+/bl1VdfZdy4cVIMt6LsbDh9+maCiIgAT0+YM8do//e/4dgxaNsWxowxvvboUbYxSQIRwk5dv36d1atXExUVVeQ5Sil69epF3759S63WcePGjQLF8MzMTIKCgvjoo4+YNGlSmd0aEzclJEBkJJw6BVOmGMfGjoWlS43vlYImTaBPn5vP+eUXqFEDrHn3UBKIEHYmOzub3bt3s379etLT04s8r1atWowcOZKAgIB7vmZmZibr16/PmxmelJRE3bp1eemllzCZTAQFBd3zNURhGRlQubKREJYsgW++MXoW584Z7UoZicPTE/70Jxg+3OhZtGljHMuvZk3rxy8JRAg7cvHiRVauXHnbobmVK1fm/vvvp3v37vdUrNZas2/fvrxl0qOjo6latSoTJkwgJCSE++67T4rhpSghAXbtMhJEZKTxNSoKDh2Cxo0hJgbOnjV6FUFBRqJo1w48PIznjxxp0/AtkgQihB1ISUlh/fr1t11yHYw1rIYPH46Pj0+Jr3Xy5Mm8meFHjhzBxcWF4cOHExISwrBhw3BzcyvxawtISTESQ26d4k9/MhLCr7/CuHHGObVqGclh0CDIHe/w1FPGw5FIAhHChrTW/PHHH6xbt44bN24UeZ6HhweDBw+mbdu2JVrDKj4+Pq8Yvn37dsDY7/zll19m3Lhx95SQKqrsbDhzBlxcoE4do4A9YoTxNXd6jrs7dO9uJJB+/WDDBqNnUaOGTUMvNZJAhLCRixcv8vPPP3P+/PnbntexY0ceeOAB3N3di/X6N27cYPny5ZjNZlavXk1mZiatW7fmH//4B8HBwdSvX/9ewq9wMjLg669v9iwOHIDr1+GNN+DDDyEgAFq2hAkTjIQRFGTcmspdPaZGDSOJlCeqqL0CyqPOnTvrO90iEKKsJScns379evbt23fb8/z8/Bg+fHixPuizsrJYv349ZrOZ8PBwkpKSqFOnDsHBwXnFcHvartbeZGbC0aM3axQREdCqFXz0kTE5z8fHGOXUtu3NGkXv3kbiKK+UUnu11p0ttUkPRAgrycrKYvfu3WzcuJG0tLQiz3NxcaFfv3507dr1rhY/zC2G584Mj46Oxtvbm/Hjx+cVw621iKIjiYkxEkRi4s3aRKdOxjEwRkcFBkL79sbPShm3p3x9je+FJBAhrOL48eOsWbOG+Pj4257XunVrBg0ahLe39x1f89SpU3kzww8fPoyzszPDhw/HZDLx4IMPSjE8R3q6UacA+PJLWLTISBK5q8H4+99MIK++avQ02rY1ehW3Tq63xVBZeyYJRIgyFBcXx9q1azl+/Phtz6tZsyZDhw6l0R0WLrp8+TILFy7EbDazdetWwCiG//Wvf2XcuHFUr1691GJ3RDExN4fK5j5OnzaG0Lq6GsNkExONYnfbtjeHy+YKCbFV5I5JEogQZSA5OZmNGzeyd+/eIvckB3B1daVv3763vV2VkpLCihUrCA0N5eeffyYzM5NWrVrxj3/8g0mTJtGgQYOyeht269q1m3WKyEh46y1jaOzcuUYvAoxFBIOCjN5FWpqRQD780HiI0iEJRIhSlJGRwY4dO9iyZcttZ5EDdOjQgQEDBuB565RijHrJxo0bCQ0NJTw8nOvXr1O7dm1eeOEFQkJCaNeuXYUohmdmGnWHmjWN2sNvv8HUqUavIlfVqjB5spFAJkyAXr2Mmdo22u69QpEEIkQpyM7OZv/+/WzYsIHr16/f9tx69eoxZMiQQps85c4JyZ0ZfvHiRby9vRk3bhwhISH07du33BfDExLgu+9u3n46eNDoPXz9Nfz5z0aS6N4dnnji5kio+vVvFrXr1TMewjrsahivUmo6MAZoAaQBO4DpWusD+c5RwN+AJwAfYCfwjNb64J1eX4bxitKmtebo0aP8+uuvt92jA6Bq1aoMHDiQVq1aFeg9nD59Oq8YfujQIZydnRk6dCghISE8+OCDxZ7/Ye9SU43lO/LXKR58EF54Aa5cMeZLBATcrFHkTsKTaSu24UjDePsBs4HdgALeA9YppVppra/knPMa8DIwFTgCvA38opRqobW+/a9+QpSis2fPsm7dOs7lrnxXBBcXF/r06UO3bt3y9um4fPkyixYtwmw2s2XLFgB69+7N//73P8aNG0eNcjBVWWujaB0RYQyJHTrUmKHt72/UMADc3Ix5FrnLeVSvboyO8vW1Xdzi7tlVD+RWSqkqQCIwWmu9Iqf3cRH4j9b6w5xz3IFY4BWt9Ze3ez3pgYjScOnSJdavX3/HkVVKKTp16kS/fv3w9PQkJSWFlStX5hXDMzIyaNmyJSEhIQQHB9OwYUPrvIEykFukBnjvPWNp8chIY8QTGLedclZQ4b//NXoZ7dpB06Y3Z2oL++RIPZBbeQFOwNWcnxsBAcDa3BO01ilKqU1AT6BQAlFKPYFxu0uWbhD3JDY2lt9+++22+3PkCgwMZMCAAfj4+PDbb78RGhrK4sWLuX79OrVq1eK5554jJCSE9u3bO1wx/PRp2L274KqyuZsdAZw4YXwNDr55Cyr/1uh/+Yu1IxZlxd4TyGfAH0DO7y7kbnwQc8t5MUAdSy+gtf4K+AqMHkjphyjKu8uXL/Pbb78RGRl5x3Pr1avHgAEDuHr1Kh999BFhYWFcuHABLy8vxo4di8lk4v7773eIYnhc3M0EceAA/O9/xq2ojz+G2bONJT1atIAuXYzeRHa2cSx3hzxR/tltAlFKzQJ6A7211lm2jkdUPJcvX2bTpk1ERkbedi4HGBMBAwMD2bFjB8OHD+fgwYNUrlyZYcOGMWvWLEaMGGG3xfC0NKOo3ayZsUnRvHnw8ssQHX3zHH9/uHjRKGS/8AI89phRu5DJ7hWbXSYQpdQnwETgfq31yXxNuf+k/YGz+Y7752sT4p7ExcWxZcuWu0ocLi4u3LhxgwULFrBp0yYAevXqxX//+18efvhhuyyGnzsHZvPNEVCHD0NWFqxbBwMGGMNgBw++OUs7KMhIILmaN7dd7MK+2F0CUUp9BkzASB6Hb2k+hZEoBmKM1EIp5Qb0AV61Zpyi/ImOjmbz5s13rHFkZGRw7tw5zpw5w/bt20lPTycwMJAPPviA4ODgOy5HYg1JScZtp9wkERkJzz1nzMqOjYXp043eRLt2MHr0zUQBxo54+ffaFqIodpVAlFJfAJOB0cBVpVRuzSNJa52ktdZKqU+BN5RSh4GjwJtAEjDPBiGLcuDcuXNs3ryZY8eOFXlOdnY2Z86cISoqikOHDpGUlERAQADPPPMMJpOJjh072qQYnpVlFK0jI8HPz/jgj4szvs/l5WUkiNzw2rWDq1ehWjWrhyvKGbtKIMDTOV9/veX4u8A7Od//C3AHvuDmRMJBMgdEFIfWmmPHjrFt2zbOnDlT5HnR0dFERERw8OBBEhMTqVKlCmPGjMFkMjFgwACrFsNTU2/WHJ55xhgJdeCAsYUqwKRJRgLx9TX2rwgMNHoVDRoUXH68cmVJHqJ02PU8kNIm80BEVlYWBw4cYNu2bcTGxlo8JyEhgQMHDhAREUFsbCyVKlVi8ODBTJ48mZEjR+Lh4VHmcR46BHv3Fhwq26wZbNxotD/wgDFRL3eYbFCQUdS20zq9cGCOPA9EiFKRmprK3r172blzp8W1qlJSUoiKiiIiIiKvR9K2bVtmzJhBcHAwvmUwNVpruHDhZoKIjoZPPjHaXn4Zfv7Z2MeiZUsjYfTocfO569aVejhCFJv0QES5duXKFXbu3Mnvv/9ORkZGgbaMjAyOHTtGREQEx44dIysrCz8/Px566CFeeuklmpficKOkJGNhwC5djLkS//oX/POfRi0iV8OGxnaqzs6wf7/xtVmzm8t8CGEL0gMRFYrWmlOnTrFz506OHj1aoC23GB4ZGcnBgwdJS0vD09OT++67j8cff5zx48eXSl3j4MGbO99FRhqFbq2Nr40bG4/x4wtuapS/LtGu3T2HIESZkwQiyo20tDQiIiLYvXt3oZVxY2JiiIiIIDIykmvXruHi4kLLli0ZMGAAjz32GC1btiz2KKrLlwtuahQRAZ99Zqz7dOiQsSZUs2bGntqTJxtJIXdL1HHjbm6jKoSjkgQiHF5sbCx79uxh//79BTZxSkxMJDIyksjISGJiYlBK0bRpUwYNGsTIkSPp168f9e5i84j0dDhyxEgQrVsbCWH3buja9eY5uYsDZmYaPz/4oHHbygr1diFsRhKIcEiZmZlERUWxd+9ezp69uShBSkoKhw4dIiIigtM5q/vVrVuXYcOG0b59e3r16kX37t0t7h2utTFU1t0dkpPhySdvztTOLZ+88YaRQFq2NNaEyt3UqFatgkNlZYkPURFIAhEOJTY2ln379hEREUFKzgSIzMzMvGL40aNHycrKokaNGvTr14+2bdvSuHFjunbtSseOHXHL98m+Zw/88cfN208RETBypLEjnocH/P67UdgePvxmoggMNJ5bpQq88or1378Q9kQSiLB7qampHDx4kD/++IPz588DRjH87NmzREREEBUVRWpqKp6ennTu3JmgoCBq165No0aN6Ny5Ky4uLThwwImPPjJGQL31lvG6kycbvQtPTyM5jBtnDJcFozdx8I57XApRsUkCEXYpOzub06dPs3//fqKiosjMKS7ExMTk1TUSExNxdnamZcuWtG3bllq1WpGcXIthw2rSpUsXPvjAn2eeMW5HgZEUeva8mUB++MHYAa9RIyOxCCGKR+aBCLuSmyAiIiLyJvwlJibmzQzPXwyvXXss6elDiI+vRVxcbRISqlC5siYpSeHqauxZcfiwUdzOnant6WnjNyiEg5F5IMKu5S4dEhkZmbe8iHHbKorff7/E+fM+QDvc3Z/Ey6sDwcFmatXKZvv2Xmzd2odmzTIZPtw5Z0kPlVfMfvrpoq8phLh3kkCETVy/fp2oqCgOHDjA+fPnSU+vTHR0DSIj/bl69UdOnVpPVtYYYFfec5ydE/HziyUgoClDhzbhmWfa4u3thIuLi+3eiBAVmCQQYTWJiYkcPHiIAweOcunSKRISvFmzZhAXLtTg2jU/wChEuLjsp1OnRBo3bkRCwir8/WNp0CCRbt2a0aFDBwICJjrcPuJClEeSQESZycqCVasS2LDhMnv2pHPqVBViYzvQs2cKrVrtZO/eXRw5Mors7K04OR2kbt2rtG/vRFCQH5UrD8PJyYmmTRMJCupCixYtqFxZ/rkKYU/kf6S4ZxkZxiKAkZGwf7/Gw+MK3brt5fDhY7z88hNkZlbDzS0FX98L+PuvZv9+Mxs3LkIpRZMmgwkKCqJFixa4ulYDoE6dOrRt25Y2bdrgKVVvIeyWJBBx17SGmBi4dAk6dDCOjRwJa9Zo0tONW0pOTtm0aXOe7OztAEyY8A1xcbs5enQDp0+fAowE0b37UFq3bk2VKlUAqFmzJq1bt6Zt27YWZ4kLIeyPJBBxWz//DGvX3pytHRcHDRpoNm48w8mTJ1HKm86d0/D3j8HfPxZf3zggncOHjxMREcGRI0fIysrCx8eHvn37EhQURI0aNQCoUaMGrVq1ok2bNvjl34NVCOEQJIFUcFrDmTM3l/LIXftp3z5j69MVK+D77zXNm2fSvfsVqlU7h4vLYebMOQFAx47G62RnZ3Pu3DlWrzaWSU9JScHDw4NOnTrRtm1b6tati1IKPz8/AgMDadWqFX5+flIMF8KBSQKpQBITjT20IyJg4kTw8TH2zp4+/eY5jRtDUJDmyJEYrl07RZcuF6hf/xRpaTcsvmZsbGzexL/ExEQqV65MYGAgQUFBNGnShEqVKlGvXj2aN29Oy5Yt83ofQgjHJwmkHMrMNEZAuboaPYl33jGSRs5OrQC0aAH9+8OwYeDllYG/fwyenqeIjz/N+fPnWbw4vcjXv3btWt7M8OjoaJRSNG7cmP79+xMYGEiVKlVo0qQJzZo1o3nz5nl1DiFE+SIJxMGlpMCWLQVXlI2Kgq+/NhYLVMrYBa9nT2N58tats6hdOx6tz7Fs2QUuXrxIfHws8fG3v05qaiqHDh0iMjKSU6dOobWmdu3aDB48mDZt2tCoUSOaNm1K06ZNadCggQy5FaICkP/lDiI11UgMuYmiWzdjS9SEBBg0yDinVi1jzacBA4z9KtLT0/H1jeH776OJjo7m0qVLRETE8vvvWXd1zczMTI4fP05kZCRHjhwhMzMTHx8f+vTpQ5cuXejWrRuNGzemSZMmVMu/H6sQokKQBGJntIazZ43d7Fq3Nn7u0MGoXWTlfO67uRmP8eMhIADWrs2kVq3LaB1LXFwcsbGxbNsWy08/XS3B9TXnzp0jIiIirxju7u5Op06dGDx4MAMHDqRx48YEBATgJEvYClGhSQKxA/PmGbehcvfWvnYN+vWDDRuMW1B9+8KDD2bTqFESdepcxts7loSEy4SGXiE+Pp7ExMR7jiEuLi6vGJ6QkICzszOdOnVi9OjRPPzwwzRo0IBKlSrd+5sVQpQbkkCsIDMTjh+/WaOIjDT2qFi3zmj/4Ydstm9XtGiRzvDhydSrl0i9ejEsXnyexMRE6tdPICkpifPnIWc/pVJx/fr1vGL4pUuXcHJyokuXLphMJh555BG8vb1L72JCiHLHYROIUupp4FWgFnAQeFFrvdnacWitycrKIjMzk4yMDC5dymT/fs3Bg4pJk+LJzEznrbf8WLKkJmDM1K5d+zp168bz+ec/k5KSTIcO2XTvnl5gT+3Ll41HaUtLS+PkyZNERUVx8OBBsrOz6dChA6+99hoTJ04kICCg9C8qhCiXHDKBKKUmAJ8BTwNbcr7+rJRqpbU+W1rXCQ8P5/Tp0+RuuqW1RmtNdnY2aWkKrTNQKovjx5uwbVsPYmL8SU72yXt+QsIP+PgkUK1afUaP9sHfPwZf3zicnY1ixtWcEkW+bbpLnY+PD76+vpw5c4YtW7awYcMGUlJSaNSoEW+88QYmk4nA3I2+hRCiGBwygQAvAd9rrb/O+fk5pdQQ4C/A9KKfVjypqakkJSWRmurKmTP1iYnxJzbWj5gYf+LjfZk69XsaNDhHZmZlUlLcadbsGP7+sfj5Gct6VKli7KXaoMFZGjQotbxmkbe3N35+fvj6+lKzZk1q1qzJyZMnWbhwIe+//z7x8fHUqFGDP/3pT5hMJnr06CGzwIUQ98ThEohSygXoBMy8pWkt0LMsrhkbW5OwsGAAqlW7ip9fLC1bHqZKlSQAAgOPEBh4pCwuXYCXlxc+Pj5Ur14976uvry/Vq1fP21TpyJEjmM1mzGYzJ0+exM3NjVGjRhESEsKgQYNk8yUhRKlxuAQC+AKVgJhbjscAD9x6slLqCeAJgPr16xfrQrm/odeqFc2jj36Dn18sbm5Fz9C+F56ennh5eeHl5UWVKlWoWrUq3t7eVK1alWrVquHt7V3k5Lzo6Gjmz5+P2Wxmz549ODk50b9/f95++20eeughKYYLIcqEIyaQYtFafwV8BdC5c2ddnOfmJhBn50zq1y96+FOlSpWoXLkylStXxtnZGRcXF5ydnXF1dc17uLm54e7ujru7Ox4eHnkPT09PPDw8ij2n4vr16yxdupTQ0FDWrVtHdnY2HTt2ZNasWUyYMIHatWsX6/WEEKK4HDGBxANZgP8tx/2B6NK80KhRo8jKyspLJEoplFI4OTnlPSpVqmS1WkJGRga//PILoaGhLF26lJSUFBo2bMi0adMwmUy0atXKKnEIIQQ4YALRWqcrpfYCA4FF+ZoGAuGleS0PD4/SfLkS0Vqzc+dOzGYzCxYsIC4ujurVq/PII48QEhJCz549pRguhLAJh0sgOWYBc5VSu4CtwFNAbeB/No2qFB09ejSvGH7ixAnc3NwYMWIEkydPZvDgwVIMF0LYnEMmEK31AqVUDeBNjImEB4BhWuszt3+mfYuJickrhu/evRulFP379+fNN99kzJgxUgwXQtgVh0wgAFrr2cBsW8dxr5KSkli6dClms5lffvmFrKws2rdvz8yZM5k4cSJ16tSxdYhCCGGRwyYQR5aZmVmgGH7jxg0aNGjA66+/LsVwIYTDkARiJVprdu/eTWhoKPPnzycuLg4fHx8mT56cVwyX5dGFEI5EEkgZO378OGazmdDQUI4fP46rqysjR47EZDIxZMgQXF1dbR2iEEKUiCSQMhAbG8uCBQsIDQ1l165dKKW4//77mT59OmPHjqVq1aq2DlEIIe6ZJJBSkpyczLJlywgNDWXt2rVkZWXRrl07/vWvfzFp0iTq1q1r6xCFEKJUSQK5B5mZmaxbtw6z2cySJUtITk6mfv36vPrqq5hMJtq0aWPrEIUQosxIAimm3GK42Wxm/vz5xMbGUq1aNUwmEyaTid69e0sxXAhRIUgCuUu5xXCz2cyxY8dwcXFhxIgRhISEMHToUCmGCyEqHEkgd5CcnMyAAQPYuXMnSin69u3L66+/ztixY6lWrZqtwxNCCJuRBHIHnp6eNGvWjLFjx0oxXAgh8pEEchfmzp1r6xCEEMLuSLVXCCFEiUgCEUIIUSKSQIQQQpSIJBAhhBAlIglECCFEiUgCEUIIUSKSQIQQQpSIJBAhhBAlorTWto7BapRSccCZEj7dF4gvxXAcgbzn8q+ivV+Q91xcDbTWNS01VKgEci+UUnu01p1tHYc1yXsu/yra+wV5z6VJbmEJIYQoEUkgQgghSkQSyN37ytYB2IC85/Kvor1fkPdcaqQGIoQQokSkByKEEKJEJIEIIYQoEUkgQgghSkQSyB0opZ5WSp1SSqUqpfYqpfrYOqayopSarpTarZS6ppSKU0qtUEq1sXVc1pTzZ6CVUv+xdSxlSSlVSyk1J+fvOVUpFaWU6mvruMqKUqqSUur9fP+XTymlPlBKlZtdWZVS9ymlliulLuT8G556S7tSSr2jlLqolEpRSm1USrW+l2tKArkNpdQE4DPg70AHYBvws1Kqvk0DKzv9gNlAT6A/kAmsU0pVt2VQ1qKU6g48AUTYOpaypJSqBmwFFDAcaAk8B8TaMKyy9jrwDPA8EAi8kPPzdFsGVcqqAAcw3luKhfbXgJcx/q67YPx9/6KU8irpBWUU1m0opXYCEVrrx/MdOwYs1lqXp394FimlqgCJwGit9Qpbx1OWlFJVgX3An4G/AQe01s/aNqqyoZT6O9BXa93L1rFYi1JqJXBZa/1IvmNzgBpa6wdtF1nZUEolAc9qrb/P+VkBF4H/aK0/zDnmjpFEXtFaf1mS60gPpAhKKRegE7D2lqa1GL+hVwReGP9Grto6ECv4CuMXgw22DsQKRgM7lVILlFKxSqk/lFLP5nzIlFdbgPuVUoEASqlWGL3sVTaNynoaAQHk+zzTWqcAm7iHz7Nyc/+vDPgClYCYW47HAA9YPxyb+Az4A9hu4zjKlFLqcaApEGLrWKykMfA08AnwT6A98H85beW19vMRxi9EUUqpLIzPvg+11rNtG5bVBOR8tfR5VqekLyoJRFiklJoF9AZ6a62zbB1PWVFKtcCocfXWWmfYOh4rcQL25LsN+7tSqhlGTaC8JpAJwBQgGDiIkTQ/U0qd0lp/Y8vAHJncwipaPJAF+N9y3B+Itn441qOU+gSYBPTXWp+0dTxlrAdGb/OgUipTKZUJ9AWezvnZ1bbhlYlLQNQtxw4B5XVwCMDHwEyt9XytdaTWei4wi/JVRL+d3M+sUv08kwRSBK11OrAXGHhL00CM0VjlklLqM24mj8O2jscKlgJtMX4jzX3sAebnfJ9uk6jK1lagxS3HmlPyvXIcgQfGL4T5ZVFxPgNPYSSKvM8zpZQb0Id7+DyTW1i3NwuYq5TahfGf7imgNvA/m0ZVRpRSXwCTMYqsV5VSufdNk7TWSTYLrAxprROAhPzHlFLJwBWt9QFbxGQFnwDblFIzgAUYQ9SfB96waVRlawUwTSl1CuMWVgfgJeAHm0ZVinJGTTbN+dEJqK+Uao/xb/msUupT4A2l1GHgKPAmkATMK/FFtdbyuM0Do9h4GkjD6JHcZ+uYyvC96iIe79g6Niv/OWzEGO5o81jK8D0OB/YDqTkfJs+TM6y/PD4wCuifYvSyUoCTGLUvN1vHVorvsV8R/3+/z2lXwDsYtzBTgd+ANvdyTZkHIoQQokQqyv0/IYQQpUwSiBBCiBKRBCKEEKJEJIEIIYQoEUkgQgghSkQSiBBCiBKRBCKEEKJEJIEIIYQoEUkgQgghSkQSiBA2oJSaopS6fOtqv0ops1Jqua3iEqI4JIEIYRuLMP7/jco9kLOt7kOA7E8hHIIkECFsQBvbiZqBR/MdDgauAT/ZJCghikkSiBC28zUwUClVN+fnR4E5WutMG8YkxF2T1XiFsCGl1G5gGcbGVpFAoNb6iE2DEuIuyYZSQtjW18BrGNvqbpXkIRyJ9ECEsCGllBfGBj/OwFNa6+9sHJIQd01qIELYkNb6OrAQY8fLhTYOR4hikQQihO3VAhZorZNtHYgQxSE1ECFsRCnlA/QBBgHtbByOEMUmCUQI2/kdqA68obU+YOtghCguKaILIYQoEamBCCGEKBFJIEIIIUpEEogQQogSkQQihBCiRCSBCCGEKJH/D196XAwQwUERAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "y = np.linspace(0,10,100)\n", "\n", "y0 = 3\n", "y2 = y**2\n", "y0y = y0 * y\n", "yta = y0**2 + 2*y0*(y-y0)\n", "\n", "plt.rc('font', size=14)\n", "plt.plot(y,y2, color='grey', lw=5)\n", "plt.plot(y,yta, 'k-')\n", "plt.plot(y,y0y, 'b--')\n", "plt.xlabel(\"y\")\n", "plt.ylabel(\"yy\")\n", "plt.legend(['exact', 'Taylor', 'simple'], frameon=False);" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Nonuniform grids\n", "#### Method 1\n", "* Set an arbitrary grid with arbitrary $\\Delta x$ spacing between points.\n", "\n", "```\n", " * * * \n", "i-1 i i+1\n", "```\n", "\n", "* Let $\\Delta x_{i-1} = x_i - x_{i-1}$, and $\\Delta x_{i} = x_{i+1}-x_i$.\n", "* Then a central difference approximation is \n", "$$f^{\\prime}(x_i) = f_i^{\\prime} \\approx \\frac{f_{i+1}-f_{i-1}}{\\Delta x_{i-1} + \\Delta x_i}.$$\n", " * A Taylor Series gives second order when $\\Delta x_{i-1}=\\Delta x_i$, but only first order when $\\Delta x_{i-1}\\ne \\Delta x_i$.\n", " * These are second and first order **asymptotically** as $\\Delta x\\rightarrow 0$. \n", " * In practice, accuracy is not severely compromized if $\\Delta x_{i-1}$ is not too different from $\\Delta x_i$.\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "\n", "* For a second derivative, we have\n", "$$f^{\\prime\\prime}_i\\approx \\frac{\\frac{f_{i+1}-f_i}{\\Delta x_i}-\\frac{f_i-f_{i-1}}{\\Delta x_{i-1}}}{\\frac{\\Delta x_{i-1}+\\Delta x_i}{2}}.$$\n", " " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "#### Method 2\n", "* Stretch the grid **analytically**\n", "* Let $x$ be the nonuniform grid and\n", "* Let $\\eta$ be a corresponding uniform grid.\n", "* For example\n", "$$\\eta = \\ln(x+1) \\rightarrow x = e^{\\eta} - 1.$$\n", "\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "\n", "* Then \n", " $$ \\frac{dy}{dx} = \\frac{dy}{d\\eta}\\frac{d\\eta}{dx} = \\frac{dy}{d\\eta}\\left(\\frac{1}{x+1}\\right).$$\n", " * This allows us to transform all the $dy/dx$ derivatives on the nonuniform grid to $dy/d\\eta$ derivatives on a uniform grid times some known factor $1/(x+1)$ at any given point. \n", " * So, if we are evaluating $dy/dx$ at point $i$, we would have \n", " $$y^{\\prime}_i = \\frac{y_{i+1}-y_{i-1}}{\\Delta \\eta}\\left(\\frac{1}{x_i+1}\\right).$$\n", "* This works great, and is easy to use. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Questions\n", "- What would you do to have an arbitrary grid spacing that is not defined by some known function?\n", "- Could you develop a second order approximation to a central derivative on a nonuniform grid?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "celltoolbar": "Slideshow", "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.9.7" }, "varInspector": { "cols": { "lenName": 16, "lenType": 16, "lenVar": 40 }, "kernels_config": { "python": { "delete_cmd_postfix": "", "delete_cmd_prefix": "del ", "library": "var_list.py", "varRefreshCmd": "print(var_dic_list())" }, "r": { "delete_cmd_postfix": ") ", "delete_cmd_prefix": "rm(", "library": "var_list.r", "varRefreshCmd": "cat(var_dic_list()) " } }, "types_to_exclude": [ "module", "function", "builtin_function_or_method", "instance", "_Feature" ], "window_display": false } }, "nbformat": 4, "nbformat_minor": 4 }