{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "###### Content under Creative Commons Attribution license CC-BY 4.0, code under MIT license (c)2014 L.A. Barba, C.D. Cooper, G.F. Forsyth." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Spreading out" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Welcome to the second lesson of the course module: \"_Spreading out: parabolic PDEs.\"_ We're studying the heat equation in one spatial dimension:\n", "\n", "$$\n", "\\begin{equation}\n", "\\frac{\\partial T}{\\partial t} = \\alpha \\frac{\\partial^2 T}{\\partial x^2}\n", "\\end{equation}\n", "$$\n", "\n", "where $\\alpha$ is the thermal diffusivity and $T$ is the temperature.\n", "\n", "In the previous lesson, we reviewed the numerical solution of the 1D diffusion equation with a forward-time, centered-space scheme: an _explicit_ scheme. What does that mean? \n", "\n", "The solution for $T$ at time step $t^{n+1}$ was calculated using different combinations of $T$ values from the *previous* time step $t^n$. We have complete knowledge of the parts that feed into the solution update at each spatial point. \n", "\n", "*Implicit* methods work differently: we will use more data from the \"future\" in the update, including several values of $T$ at $t^{n+1}$. This will make the scheme more difficult to apply, but there are several reasons why it may be worth the effort.\n", "\n", "In lesson 1, we discussed two disadvantages of explicit methods: (1) boundary effects drag behind by one time step; (2) stability requirements constrain the time-step to very small values. Both of these issues are resolved by implicit schemes." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Implicit schemes" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's move things around a bit and try combining the Euler time step with an evaluation of the spatial derivative on the *updated* solution at $t^{n+1}$. The discretized form of the equation is now as follows:\n", "\n", "$$\n", "\\begin{equation}\n", "\\frac{T_{i}^{n+1}-T_{i}^{n}}{\\Delta t}=\\alpha\\frac{T_{i+1}^{n+1}-2T_{i}^{n+1}+T_{i-1}^{n+1}}{\\Delta x^2}\n", "\\end{equation}\n", "$$\n", "\n", "The stencil for this discretization doesn't look anything like the other stencils we've used until now. Check it out." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "![stencil-implicitcentral](./figures/stencil-implicitcentral.png)\n", ".\n", "#### Figure 1. Stencil of the implicit central scheme." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "From the previous time-step, we only know $T_i^{n}$, but what about $T_i^{n+1}$, $T_{i-1}^{n+1}$ and $T_{i+1}^{n+1}$? What can we do?\n", "\n", "No need to panic! Let's start by putting what we *do know* on the right-hand side of the equation and what we *don't know* on the left. We get:\n", "\n", "$$\n", "\\begin{equation}\n", "-T_{i-1}^{n+1} + \\left( 2 + \\frac{\\Delta x^2}{\\alpha\\Delta t}\\right) T_{i}^{n+1} - T_{i+1}^{n+1} = T_{i}^{n}\\frac{\\Delta x^2}{\\alpha\\Delta t}\n", "\\end{equation}\n", "$$\n", "\n", "It looks like there are a lot of unknowns and just one equation! \n", "\n", "What does it look like with $i=1$?\n", "\n", "$$\n", "\\begin{equation}\n", "-T_{0}^{n+1} + \\left( 2 + \\frac{\\Delta x^2}{\\alpha\\Delta t}\\right) T_{1}^{n+1} - T_{2}^{n+1} = T_{1}^{n}\\frac{\\Delta x^2}{\\alpha\\Delta t}\n", "\\end{equation}\n", "$$\n", "\n", "and $i=2$?\n", "\n", "$$\n", "\\begin{equation}\n", "-T_{1}^{n+1} + \\left( 2 + \\frac{\\Delta x^2}{\\alpha\\Delta t}\\right) T_{2}^{n+1} - T_{3}^{n+1} = T_{2}^{n}\\frac{\\Delta x^2}{\\alpha\\Delta t}\n", "\\end{equation}\n", "$$\n", "\n", "What about $i=3$?\n", "\n", "$$\n", "\\begin{equation}\n", "-T_{2}^{n+1} + \\left( 2 + \\frac{\\Delta x^2}{\\alpha\\Delta t}\\right) T_{3}^{n+1} - T_{4}^{n+1} = T_{3}^{n}\\frac{\\Delta x^2}{\\alpha\\Delta t}\n", "\\end{equation}\n", "$$\n", "\n", "Can you see the common element across equations? Here's a little help:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$T_{i}^{n+1}$ also appears in the equation for $T_{i-1}^{n+1}$ and $T_{i+1}^{n+1}$. We might have enough equations if we apply this for all $i$-values at the same time, don't you think? In fact, this is a linear system of equations for the unknown values $T_{i}^{n+1}$ on the spatial grid." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### What about boundary conditions? " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's look at the boundary points of the example from the previous lesson with a Dirichlet BC at $x=0$ and a Neumann BC at $x=1$, discretizing with $N$ mesh points.\n", "\n", "The value $T_0^{n+1}$ is known at every time-step from the BC, so putting all unknown terms on the left-hand side of the equation and the known values on the right side yields the following for the $i=1$ equation:\n", "\n", "$$\n", "\\begin{equation}\n", "-T_{2}^{n+1} + \\left( 2 + \\frac{\\Delta x^2}{\\alpha\\Delta t}\\right) T_{1}^{n+1} = T_{1}^{n}\\frac{\\Delta x^2}{\\alpha\\Delta t} + T_{0}^{n+1}\n", "\\end{equation}\n", "$$\n", "\n", "That was easy!\n", "\n", "On the other hand, for $i=N-2$, the equation reads\n", "\n", "$$\n", "\\begin{equation}\n", "-T_{N-3}^{n+1} + \\left( 2 + \\frac{\\Delta x^2}{\\alpha\\Delta t}\\right) T_{N-2}^{n+1} - T_{N-1}^{n+1} = T_{N-2}^{n}\\frac{\\Delta x^2}{\\alpha\\Delta t}\n", "\\end{equation}\n", "$$\n", "\n", "The discretized Neumann boundary condition on the right side of the rod is\n", "\n", "$$\n", "\\begin{equation}\n", "\\frac{T^{n}_{N-1} - T^{n}_{N-2}}{\\Delta x} = q\n", "\\end{equation}\n", "$$\n", "\n", "But we can just as easily write that at time step $n+1$ (the boundary conditions apply at every time-step):\n", "\n", "$$\n", "\\begin{equation}\n", "\\frac{T^{n+1}_{N-1} - T^{n+1}_{N-2}}{\\Delta x} = q\n", "\\end{equation}\n", "$$\n", "\n", "Inserting the Neumann boundary condition in the equation for $i=N-2$ yields\n", "\n", "$$\n", "\\begin{equation}\n", " -T_{N-3}^{n+1} + \\left( 1 + \\frac{\\Delta x^2}{\\alpha\\Delta t} \\right) T_{N-2}^{n+1} = T_{N-2}^{n} \\frac{\\Delta x^2}{\\alpha\\Delta t} + \\Delta x q\n", "\\end{equation}\n", "$$\n", "\n", "Make sure you work this out with pen and paper: it's important to recognize where these terms come from!\n", "\n", "Now we can write the linear system of equations in matrix form as follows:\n", "\n", "$$\n", "[A][x] = [b]+[b]_{b.c.}\n", "$$\n", "\n", "where the matrix of coefficients $[A]$ is a sparse matrix—most of the matrix elements are zero—with three non-zero diagonals. We write below the system expanded out, so you can see the structure of the matrix, with $\\sigma=\\frac{\\alpha\\Delta t}{\\Delta x^2}$:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{align}\n", " \\left[\n", " \\begin{array}{cccccc}\n", " \\left( 2 + \\frac{1}{\\sigma} \\right) & -1 & 0 & \\cdots & & 0 \\\\\n", " -1 & \\left( 2 + \\frac{1}{\\sigma} \\right) & -1 & 0 & \\cdots & 0 \\\\\n", " 0 & & \\ddots & & & \\vdots \\\\\n", " \\vdots & & & & \\left( 2 + \\frac{1}{\\sigma} \\right) & \\\\\n", " 0 & \\cdots & & & -1 & \\left( 1 + \\frac{1}{\\sigma} \\right)\n", " \\end{array}\n", " \\right] \\cdot \n", " \\left[\n", " \\begin{array}{c} \n", " T_1^{n+1} \\\\\n", " T_2^{n+1} \\\\\n", " \\vdots \\\\ \\\\\n", " T_{N-2}^{n+1}\n", " \\end{array}\n", " \\right] =\n", " \\left[\n", " \\begin{array}{c} \n", " T_1^n \\frac{1}{\\sigma} \\\\\n", " T_2^{n} \\frac{1}{\\sigma} \\\\\n", " \\vdots \\\\ \\\\\n", " T_{N-2}^{n} \\frac{1}{\\sigma}\n", " \\end{array}\n", " \\right] +\n", " \\begin{bmatrix}\n", " T_0^{n+1} \\\\\n", " 0 \\\\\n", " \\vdots \\\\\n", " 0 \\\\\n", " q\\Delta x\n", " \\end{bmatrix}\n", "\\end{align}\n", "$$\n", " \n", "Notice that the Dirichlet boundary condition adds only a term to the right-hand side of the system. The Neumann boundary condition both adds a term to the right-hand side and modifies the matrix $[A]$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Problem set up" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We'll re-use the problem from lesson 1: we have a graphite rod, with [thermal diffusivity](http://en.wikipedia.org/wiki/Thermal_diffusivity) $\\alpha=1.22\\times10^{-3} {\\rm m}^2/{\\rm s}$, length $L=1{\\rm m}$, and temperature held at $T=100{\\rm C}$ on the left end, $x=0$, and $0{\\rm C}$ everywhere else initially. We'll compute the evolution of temperature on the length of the rod.\n", "\n", "Let's start like we did in the previous lesson: import your libraries and set up the discretization." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy\n", "from matplotlib import pyplot\n", "%matplotlib inline" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# Set the font family and size to use for Matplotlib figures.\n", "pyplot.rcParams['font.family'] = 'serif'\n", "pyplot.rcParams['font.size'] = 16" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# Set parameters.\n", "L = 1.0 # length of the rod\n", "nx = 51 # number of locations on the rod\n", "dx = L / (nx - 1) # distance between two consecutive locations\n", "alpha = 1.22e-3 # thermal diffusivity of the rod\n", "q = 0.0 # temperature gradient on the right side of the rod\n", "\n", "# Define the locations along the rod.\n", "x = numpy.linspace(0.0, L, num=nx)\n", "\n", "# Set the initial temperature along the rod.\n", "T0 = numpy.zeros(nx)\n", "T0[0] = 100.0" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Solving a linear system" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We need to solve the linear system of equations written above to advance the solution in time. Luckily, we can rely on our friends from SciPy who have developed some nice linear solvers, so we don't need to write our own.\n", "\n", "From `scipy.linalg`, let's import `solve`: a function to solve linear systems. Make sure to explore the documentation of [`scipy.linalg`](https://docs.scipy.org/doc/scipy/reference/linalg.html). We'll need to define our own custom functions to generate the coefficient matrix and the right-hand side of the linear system. You should carefully study the code below." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "from scipy import linalg" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "def lhs_operator(N, sigma):\n", " \"\"\"\n", " Computes and returns the implicit operator\n", " of the system for the 1D diffusion equation.\n", " We use backward Euler method, Dirichlet condition\n", " on the left side of the domain and zero-gradient\n", " Neumann condition on the right side.\n", " \n", " Parameters\n", " ----------\n", " N : integer\n", " Number of interior points.\n", " sigma : float\n", " Value of alpha * dt / dx**2.\n", " \n", " Returns\n", " -------\n", " A : numpy.ndarray\n", " The implicit operator as a 2D array of floats\n", " of size N by N.\n", " \"\"\"\n", " # Setup the diagonal of the operator.\n", " D = numpy.diag((2.0 + 1.0 / sigma) * numpy.ones(N))\n", " # Setup the Neumann condition for the last element.\n", " D[-1, -1] = 1.0 + 1.0 / sigma\n", " # Setup the upper diagonal of the operator.\n", " U = numpy.diag(-1.0 * numpy.ones(N - 1), k=1)\n", " # Setup the lower diagonal of the operator.\n", " L = numpy.diag(-1.0 * numpy.ones(N - 1), k=-1)\n", " # Assemble the operator.\n", " A = D + U + L\n", " return A" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "def rhs_vector(T, sigma, qdx):\n", " \"\"\"\n", " Computes and returns the right-hand side of the system\n", " for the 1D diffusion equation, using a Dirichlet condition\n", " on the left side and a Neumann condition on the right side.\n", " \n", " Parameters\n", " ----------\n", " T : numpy.ndarray\n", " The temperature distribution as a 1D array of floats.\n", " sigma : float\n", " Value of alpha * dt / dx**2.\n", " qdx : float\n", " Value of the temperature flux at the right side.\n", " \n", " Returns\n", " -------\n", " b : numpy.ndarray\n", " The right-hand side of the system as a 1D array of floats.\n", " \"\"\"\n", " b = T[1:-1] / sigma\n", " # Set Dirichlet condition.\n", " b[0] += T[0]\n", " # Set Neumann condition.\n", " b[-1] += qdx\n", " return b" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we'll define a function that steps in time using the implicit central-space scheme. Remember that for an implicit method, a step in time is performed by solving the entire linear system. This is a fundamental difference between implicit and explicit methods, and implies a considerable computational cost." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "def btcs_implicit(T0, nt, dt, dx, alpha, q):\n", " \"\"\"\n", " Computes and returns the temperature along the rod\n", " after a given number of time steps.\n", " \n", " The function uses Euler implicit in time,\n", " central differencing in space, a Dirichlet condition\n", " on the left side, and a Neumann condition on the\n", " right side.\n", " \n", " Parameters\n", " ----------\n", " T0 : numpy.ndarray\n", " The initial temperature distribution as a 1D array of floats.\n", " nt : integer\n", " Number of time steps to compute.\n", " dt : float\n", " Time-step size.\n", " dx : float\n", " Distance between two consecutive locations.\n", " alpha : float\n", " Thermal diffusivity of the rod.\n", " q : float\n", " Value of the temperature gradient on the right side.\n", " \n", " Returns\n", " -------\n", " T : numpy.ndarray\n", " The temperature distribution as a 1D array of floats.\n", " \"\"\"\n", " sigma = alpha * dt / dx**2\n", " # Create the implicit operator of the system.\n", " A = lhs_operator(len(T0) - 2, sigma)\n", " # Integrate in time.\n", " T = T0.copy()\n", " for n in range(nt):\n", " # Generate the right-hand side of the system.\n", " b = rhs_vector(T, sigma, q * dx)\n", " # Solve the system with scipy.linalg.solve.\n", " T[1:-1] = linalg.solve(A, b)\n", " # Apply the Neumann boundary condition.\n", " T[-1] = T[-2] + q * dx\n", " return T" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We solve the linear system for every time step, but the $A$ matrix does not change. Thus, you can generate it only once and then use it for all time steps. Let's try this out!" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "# Set the time-step size based on CFL limit.\n", "sigma = 0.5\n", "dt = sigma * dx**2 / alpha # time-step size\n", "nt = 1000 # number of time steps to compute\n", "\n", "# Compute the temperature along the rod.\n", "T = btcs_implicit(T0, nt, dt, dx, alpha, q)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now plot the solution!" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaQAAAEbCAYAAACV0PCVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xd4VGX2wPHvSSeE0BJ6IIAU6VW6BhvFigpiXSvquljWte+6Krq66q666IrYd/1ZQMGGKDaE0KT33msg1EAoKef3x70x2TFtkpnMTHI+z3Ofydz73pszLyEn975NVBVjjDEm0MICHYAxxhgDlpCMMcYECUtIxhhjgoIlJGOMMUHBEpIxxpigYAnJGGNMUAhoQhKRhiLyjYhY33NjjKniApaQRGQYMAdoWUK5SBEZIyJrRGSFiMwWkf5FlL1HRFaJyDIRWSQil/ojdmOMMb4XyDukh4DzgFkllBsLXAkMUNUOwNvAdyLSpWAhEXkI+DNwkap2Ah4EJorIEJ9HbowxxuckUDM1iEiEqmaLyLvA71RVCinTBlgN3KKqbxfYvxLYoqoXuO9rATuBf6jqYwXKTQGSVbW9fz+NMcaY8grYHZKqZpei2DBAgJ889v8InC8ice77wUBsEeXaiUjb8sRqjDHG/4K9l10nIBfY5rF/MxABtCtQLm+/Z7mCx40xxgSpiEAHUIIEIFNVczz2H3Ff6xYoB5BRQrn/ISKjgFEAYdXiuyfWq0/N6N88OaxycnNzCQsL9r9VKobVRT6ri3xWF/nWrVuXrqqJvrhWsCekopQ2axRbTlXHA+MBohu20piRL3Jux4Y8P7wTsVGhWjXlN336dFJSUgIdRlCwushndZHP6iKfiGz11bWCPcWnA7EiEu6xv4b7ur9AuYL7iypXpHqxQlx0BFOW7+ayf89m+4HMMgVsjDGmbII9IS3DiTHJY39zIBunB15eOYDkQsoVPF6k2Ajhszv70SKhOmv2ZHDxK6nM3pBe0mnGGGN8JNgT0mRAgRSP/QOBaaqa12b0DZBZRLlVqrqmNN/stHpxTL6zHwPbJHIwM4vr3v6Ft1I3Y4sYGmOM/wV1QlLVtThtPA+LSAKAiNyEM7vDowXKHQLGAHeKSAu33LnAIOBP3nzPmtUiefN3PblzYEtycpUxX63ivolLOZHl2a/CGGOMLwWs5V5EnseZqaGp+36Je+gMVT1VoOho4K/ALBHJwulJd76qLil4PVV9VkROAF+JSDaQAwxX1anexhYeJtw/qC3tGtbkTxOXMmnRTjbsPcrr13WnYc1qXn9WY4wxJQtYQlLV+0tZLgtnSqA/l6LsS8BL5QztVxd0akiLxOqM+u8Clu04zEVjU3nt2u70TK7jq29hjDHGFdSP7ILB6Q3j+eLO/vQ7rS7pR09x9Rtz+b95PuvlaIwxxmUJqRRqV4/ivRvP4Jb+zcnKUR6dvIKHJy3nVHZuoEMzxphKwxJSKUWEh/HnC9vx4pWdiY4I48NftnHVG3PZe+REoEMzxphKwRKSl4Z1bcInt/elYc0YFm49yEWvpLJk+6FAh2WMMSHPElIZdGxSky/+0J8zkuuQduQkI8bNYcKC7YEOyxhjQpolpDJKrBHN+7f04vo+zTiVk8sDnyzjr5+vICvH2pWMMaYsLCGVQ1REGE9e0oFnL+tIVHgY783ZyjVvziP96MlAh2aMMSHHEpIPjDyjKR+O6k29GtH8svkAF49NZfmOw4EOyxhjQoolJB/p3qw2X47uT9emtdh1+ARXjJvNpEU7Ah2WMcaEDEtIPlQ/PoaPRvVmZM8kTmbn8scJS3nyy1VkW7uSMcaUyBKSj0VHhPPMZR156tIORIQJb8/azHVv/cJ+a1cyxphiWULyAxHh2t7N+HBUbxLiopmzaT8XvzKLFTutXckYY4piCcmPeibX4avR/emSVIudh45zxbjZfL5kZ6DDMsaYoGQJyc8a1Izh49t6M6JHE05k5XL3R0sY85W1KxljjCdLSBUgOiKcv1/eiTGXtCciTHgrdTPXv23tSsYYU5AlpAoiIlzXJ/nXdqXZG61dyRhjCrKEVME825Uuf202kxfbeCVjjLGEFAB57Up545Xu/djGKxljjCWkAImOCOfZyzvxt2EdiQx3xitd+5bNg2eMqbosIQXY1b2a8tGo3iTWiGbuJmcevKW2vpIxpgqyhBQEujerw5TR/enerDa7Dp9g+OtzmDDf1lcyxlQtlpCCRL34GD68tTfX9m7KqexcHvh0GY9OXs6pbGtXMsZUDZaQgkhURBhPXdqR5y7vRFREGP83bxsjx88h7ciJQIdmjDF+ZwkpCI3omcTE2/rQsGYMi7Yd4sKxqczfciDQYRljjF9ZQgpSnZNq8eXo/vRuUYd9GSe5avxc3pu9BVUNdGjGGOMXlpCCWEJcNO/f3Iub+zcnO1f56xcruW/CUo6fygl0aMYY43OWkIJcRHgYf7mwHf+6qivVIsOZtHgnl702m237MwMdmjHG+JQlpBBxcedGTL6zL8l1Y1m9+wgXvZLK9LV7Ax2WMcb4jCWkENK2QTyf/6E/57Stx+HjWdz47nxe+XE9ubnWrmSMCX2WkEJMzWqRvHF9D+49tzUAL0xbx6j/LuTw8awAR2aMMeVjCSkEhYUJd5/bird/15P4mAi+X53GJa+ksmbPkUCHZowxZRb0CUlEeojIVBFZLSLLReQXERnuUSZSRMaIyBoRWSEis0Wkf6BirigD29bjq9EDOL1hPFv2ZzLsVVsi3RgTuoI6IYlIMvADkA50VNWOwNvABBG5qEDRscCVwABV7eCW+U5EulRsxBWvad1YJt3Rl8u6NeZ4Vg53f7SEx79YaVMOGWNCTlAnJGAoEA/8U1WzAVR1HHAEuBpARNoAo4BnVXWfW+ZNYBPwdCCCrmjVosL5x/DOjLm0A5Hhwruzt3DVG3NtyiFjTEgJ9oSU7b5G5O0QEcGJO9zdNQwQ4CePc38EzheROH8HGQxEhOt6N+Pj2/rQID6GhVsPcsG/Upm7aX+gQzPGmFIJ9oT0EbAG+LOIxIlIGPAIEA2Mc8t0AnKBbR7nbsZJZO0qKNag0K1pbb66qz99WtQl/ehJrnlzHq//vNGmHDLGBD0p6heViGQBZfkttltVm5Urqv+NoxHwDnAWcBQ4DNykqj+7x6cBfVS1hsd5twBvAENVdWoR1x6F87iPxMTE7hMmTPBV2AGXk6tMWp/FlM1Od/Du9cO5uUM0sZFS4rlHjx4lLq5K3FiWyOoin9VFPquLfAMHDlyoqj18ca2IYo6lA296eT0BRpQ9HI+LOe1DPwBfA3WAE+71J4nItUUlmgKxFEtVxwPjAdq0aaMpKSnljjmYnHM2XLJyD/dNXMrCtGwO5AivXduNtg3iiz1v+vTpVLa6KCuri3xWF/msLvyjuIS0R1X/4u0FReSCcsTjaQxQC7hbVY+7+z4SkSuB99y7p3QgVkTCVbXgrKN5d0xVuhHl/PYN+LJ+DW5/fyFr9mRw6auzeOayjgzr2iTQoRljzP8org3pmzJes6znFaYjsKNAMsqzDkgEmgPLcD5HkkeZ5jidIlb7MJ6QlJxQncm/78fl3ZpwIiuXez9eyp8/W87JbJs13BgTPIpMSKr6cFkuWNbzirAXaCginndyzXDatw4Ck92vUzzKDASmqWqGD+MJWdWiwnlheCf+NqwjUeFhvD93GyPGzWH7AZs13BgTHMrVy87t9eZPY3HGIT3pdvdGRAYClwEfq2q6qq7FaQd6WEQS3DI3AS2BR/0cX0gREa7u1ZRP7+hLk9rVWLrjMBeOTeXHNWmBDs0YY4pPSCLS2J2q5xcRGVNIkXNF5HsRae2P4FT1E2Aw0AdYJSIrgJdxEs0NBYqOBiYCs9wytwLnq+oSf8QV6jo2qclXo/NnDb/p3QW88O1acmzWcGNMABXXqQGcQaddgedwpuPxtBI4BkwXkZ6q6vOJ1FT1W+DbEspkAX92N1MKtWKjeOP6HoybsZEXvl3LKz9tYNG2g7w8smugQzPGVFElPXK7ABitqo+q6mbPg6q6U1UvAf4DPOCPAI3/hIUJv085jfdv6UVCXBSzN+7ngn/NZO0B6+xgjKl4JSWkxjiDS0vyJM7AVROC+rZMYMpdAzgjuQ57M07y9/kneG36Rlv4zxhToUpKSLkeY3sKpaqZONP3mBBVPz6GD27txW1ntSBX4e/frOGW/yzg4LFTgQ7NGFNF+LKXXKQPr2UCICI8jIeHnM493aKpWS2SH9fs5cKxqSzadjDQoRljqoCSEtJmESnxUZyInAls9U1IJtC61Itgyl396ZxUi52HjjNi3BzeSt1sE7QaY/yqpIT0MvCBiAwoqoC7Muv/Ac/7MjATWE1qxzLxtj7c1K852bnKmK9Wccf7izh8PCvQoRljKqliu32r6nQRGQ/8LCJLgblAGs7MCA2A3jjLPzyYN/u2qTyiIsJ47KJ29EyuzQOfLOOblXtYufswr1zVjc5JtQIdnjGmkilpHBKq+oSILAOeAG7zOLwMuFhVp/gjOBMchnRsSLtG8dz5wSJW7DzCFeNm8/CQ07mxXzLuBBrGGFNuperUoKqTVbUTTjfwfu7WRFW7WDKqGprVrc6nd/Tlhr7JZOUoT361itv+u5DDmfYIzxjjG171slPV3ao6x912+SsoE5yiI8J5/OL2jLu2GzViIpi2Ko2h/5rJYuuFZ4zxgSITkog8XpYLlvU8EzoGd2jIlNED6NykJjsPHWf4uDm8OXOT9cIzxpRLcXdIF5fxmmU9z4SQpnVjmXh731974T01ZTU3v7eAAzaQ1hhTRsV1amggIg9TiqXAPVQvRzwmhOT1wuvdog73f7KMH9fsZcjLM3h5ZFd6t6gb6PCMMSGm2IQEPF2Ga+4pYywmRJ3fvgHtG9fk7g8Xs2DrQa5+Yy53ndOK0We3IjzMeuEZY0qnuEd2kWXcmvgxXhOkGteqxkejenPnwJYo8NL367n6jbnsOXwi0KEZY0JEcUuY55Rxs0lWq6iI8DDuH9SW/97Ui4S4aOZtPsDQf820FWmNMaXi7yXITRXUv1UCU+8ewIBWCRw4doqb3l3AE1+u5GS2rbNkjCmaJSTjF4k1onnvxjN4cHBbIsKEd2ZtYdirs9mw92igQzPGBClLSMZvwsKEO1Ja8skdfWlaJ5ZVu49w0dhUJszfbmOWjDG/YQnJ+F2XpFpMuas/l3ZpxPGsHB74dBmjP1xsM4cbY/6HJSRTIWrERPLSyK78c0RnqkeF89Wy3Qx9eSYLtx4IdGjGmCDhdUISkb4i8oiIPO2+7y8iNhjWlMpl3Zrw1V0D6Ng4f9qhF79bR3aOdc40pqordUISkeoiMhVIBZ4CbnQPXQwsF5Fkn0dnKqXmCc7M4bef5YxZevmH9Yx4fQ7bD2QGOjRjTAB5c4f0LFAbuAhoAewFUNUHgMeAZ3wenam0oiLCeGhIW/7vll40iI9h0bZDDHl5JpMX7wh0aMaYAPEmIQ0FzlPVKaq6Bfj1GYuqvg+09nFspgro2zKBb+4ZwJAODTh6Mpt7P17K3R9ZhwdjqiJvElKWqmYUczy+vMGYqqlWbBT/vqYbf7+8I9Uiw/l8yS6GvjyTeZv2Bzo0Y0wF8iYhZYrIJYUdEJFBwCHfhGSqIhHhyp5NmXJXfzq56yyNfGMuz05dw6ls6/BgTFXgTUJ6GpgkIj+IyBggQUQeEpEPgc+BJ/0SoalSWiTG8ekdfRl99mkIMO7njVz66izWpRV3c26MqQxKnZBU9VPgOqAN8CjOrN5/A/oDv1PVL/0SoalyIsPDuO/8Nky8vQ9JdaqxavcRLhybyjuzNpObazM8GFNZeTUOSVU/AJKAjsBA97Wpqn7sh9hMFde9WR2m3n0mI3o04VR2Lk98uYrfvfOLLWlhTCXlzTikfSKyE2iuqitV9Wf31e9/sorI5SIyQ0QWisgmEVkgItcVOB4pImNEZI2IrBCR2SLS399xGf+Li47guSs6M+7a7tSOjWTm+nQGvTSDL5buCnRoxhgf8+YOKQzooaqb/BVMYUTkXpxHhFeranecR4brgHMKFBsLXAkMUNUOwNvAdyLSpSJjNf4zuEMDvr3nTFLaJHL4eBZ3fbiYP3ywiIPHTgU6NGOMj3iTkFaq6u6iDorIhT6Ix/OayTgDcm9T1R0AqpoF/Al4xS3TBhgFPKuq+9wybwKbKNsS7CZI1YuP4Z0bevK3YR2JdefDG/TSDH5auzfQoRljfMCbhPSmiPxJRIo6xx+97K4DDqnq/II7VXWXqi5w3w4DBPjJ49wfgfNFJM4PcZkAERGu7tWUqXcPoEez2uzNOMmN78znkcnLOXYyO9DhGWPKwZuEdB1wH7DHbaOZVnADWvohvr7AFrcNaabbRjRbRG4qUKYTzqwR2zzO3QxEAO38EJcJsGZ1q/PxbX14aEhbosLD+GDeNoa8PJNfNtvs4caEKiltnwQRyQCWFFOks6r6dLYGEVkBJAPLce6E9gKXAx8Cf1XVp91k2EdVa3icewvwBjBUVacWcf1ROI/7SExM7D5hwgRfhh+yjh49Slxc6NxYbs/IZfyyk2zPyEWA85tFcHnrKKLCpdzXDrW68Ceri3xWF/kGDhy4UFV7+OJa3iSkxaratazHy0JENuDceQ1Q1dQC+z8FBgOJwGcUnpBuBcZTTEIqqE2bNrp27Vpfhh+ypk+fTkpKSqDD8Mqp7FzG/rief0/fSE6u0iKhOs8P70z3ZrXLdd1QrAt/sbrIZ3WRT0R8lpC8eWR3UQnHB5QnkCLkDc/3vDNbDMTiPI5LB2JFJNyjTF6CsgnRqoCoCGcw7eTf96VVvTg2pR9j+LjZPDN1NSeycgIdnjGmFLyZqaGkdQFGlzOWwqxxXz3jzCmwf5n7muRRpjmQDaz2Q1wmSHVqUosvR/fnjhSnSfP1nzdx4dhUlm63qRaNCXbeDIztW9wG3OCH+PKmI+rksb8DcBxYCUwGFEjxKDMQmFbCDOWmEoqJDOfBwW359I6+tEiszoa9Rxn271k8O3WN3S0ZE8S8eWSXCswsZjvN59HBx8B84Km87tsiMgC4AnhaVY+p6lqctqKHRSTBLXMTTtvTo36IyYSIrk1r8/VdA7h1QHMUZ6LWof+aycKt1hPPmGAU4UXZjcDtHvuqA21xljF/2VdB5VHVHBEZDPwdWCkiJ4CTwB9U9Y0CRUcDfwVmiUgWTtvT+apaXK9AUwXERIbz6AXtGNKxIfdPXMrGfce4YtwcbuzbnPsHtaFalGfTozEmULxJSE+r6g+F7P9CRN4BXgIm+iasfKp6ALi1hDJZwJ/dzZjf6Na0NlPuGsC/fljP6zM28faszfywJo1nL+tEn5Z1Ax2eMQbvOjW8W8yxffy2nceYoBITGc4Dg9vy2e/70bZBDbbuz+SqN+by6OTlZJywJdONCTRvOjU0KmRrLCIdROQvQIwf4zTGZzo2qckXf+jPPee2IiJM+L952zjvnzP4blVaoEMzpkrz5pHdDpzebIU5AdxY/nCMqRhREWHcc25rBndowIOfLmfp9kPc+p8FXNCpIY9f1J7EGtGBDtGYKsebhLSd306gmoMznc88VbUBqCbktG0Qz6Q7+vLe7C08/+1apizbTer6dB694HSGd28S6PCMqVK8SUifqupbfovEmAAJDxNu6t+c89rV55HJy5m5Pp0HPlnG50t2ckmj3ECHZ0yV4c04pM8L2ykivURknIg091FMxgREUp1Y/nPTGfxzRGdqxUYya8N+Hk09zqs/beBUtiUmY/zNm4T0YhH79wCHgffLH44xgSUiXNatCT/88Swu7dKIrFx4/tu1XDQ21QbUGuNn3iSkQufyV9WtqvogYHOxm0qjblw0L43syv09YmhWN5a1aRlc/tocHpm8nMOZ1kXcGH8otg1JRC4if5bvJiIyvrBiQBMg0sexGRNw7RPC+faSAbzy4wZen7GRD+ZtY9rKNB67qB0XdWqISPnXXDLGOErq1HAaMMT9ukaBrwvKwlmdtdjZFIwJVTGR4fxpUBsu7tKIRyYtZ8HWg9z14WImLtjOExe3p0WiPRwwxheKfWSnqi+qapKqJgGr87722Fqo6jmqOquCYjYmIFrXr8GE2/rwzGUdqVktkpnr0xn80kz+MW2tzSJujA9404Z0qd+iMCZEhIUJV53RlB/vO4sRPZpwKieXsT9u4LwXf+aH1TbTgzHl4c1cdluLOy4iJS4TbkxlUTcumueu6Mwnt/ehbYMabD9wnJvfW8Ct/1nAjoOZgQ7PmJDkzcBYRKQhMBJoAXjOreKTNdWNCSU9kuvw1ej+vDt7Cy9+t47vVqUxc/0+7kw5jVvPbEFMpC1vYUxplTohiUhP4HvgFE4Hh33uoQSc5LTH59EZEwIiwsO4ZUALLuzUiKemrOKrZbv5x3frmLhwB49d2I5zTq9nvfGMKQVv2pCeBW5R1UQKdHDAWaTvUeAVfwRoTKhoUDOGV67uxge39qJ1/Ti2Hcjklv8s4KZ357M5/VigwzMm6HmTkBJU9TcL8Klqrqo+A5zju7CMCV19WyYw5a4BPHZhO2pER/DT2n0MenEGz32zhmMnswMdnjFBy5uEdLLgGxGpUeDraKCVr4IyJtRFhodxU//m/PinFIZ3d3rj/Xv6Rs7+x3QmLdpBbm5RK7kYU3V5k5AyROR34jwMXwBMEpELROQCYDKw2y8RGhPCEmtE8/zwzkz6fV86NalJ2pGT/HHCUoa9NptF2w4GOjxjgoo3CWkccC3QEhjjvn4BfAn0Ae7zeXTGVBLdmtbms9/344XhnalXI5ql2w9x2b9nc/dHi9l9+HigwzMmKJS6l53bfvRrG5KIdAIG4PSwS1XVdN+HZ0zlERYmXNG9CUM6NOC16RsZP3MTny/Zxbcr93D7WS257cyWVIuybuKm6ir1HZKI/OJuzQFU9aiqTlXVzywZGVN61aMj+NOgNvzwx7O4oGNDTmTl8tL360l54ScmLthOjrUvmSrKm0d2pwN3qupmfwVjTFWSVCeWV6/pxsejetOxsdO+dP8ny7hwbCqp6+1vPFP1eJOQlqnq/KIOikg3H8RjTJXTq0VdPr+zHy9d2YXGtaqxevcRrn1rHje88wvr0jICHZ4xFcabhPSZiIwo5vib5Q3GmKoqLEy4tGtjfrjvLB4c3JYa0RFMX7uPwS/N4OFJy0g7ciLQIRrjd97MZdcauFdEHgZWA0c9jif5LCpjqqiYyHDuSGnJiB5N+NcP63l/3jY+/GU7kxfv5MZ+zbn9rJbUrGZrYZrKyZuEdD2wF2fuugGFHK9RyD5jTBnUjYvmiUs6cH3fZF74di1TV+zhtenOirW/T2nJ7/om28StptLx5pHdqiIW6Pt1AT9/BWlMVdUyMY7Xru3OZ3f2o3eLOhw+nsUzU9cw8IXpTJi/neyc3ECHaIzPeJOQbivh+JXlCcQYU7QuSbX48NbevHtjT05vGM/uwyd44NNlDHppBl8t22VTEZlKwZsF+n7J+1pEGrkDYxGRMPf4Ot+HZ4zJIyKktKnHlNH9eXlkF5LqVGPjvmP84YPFXDA2le9XpaFqicmELm/ukBCRkSKyDtgO5K0Q+76IvCAVtOCLiMwUERWR5Ir4fsYEm7Aw4ZIujfnhjyk8PawDDeJjWL37CLf8ZwHD/j2b1PXplphMSPJmpoargPeApThz2eUNkHgQaA485PPofhvD5UD/Io7FicgrIrJWRFaJyDQRae/vmIwJlKiIMK7p1Yzp96fwlwvbUbd6FEu2H+Lat+Zx1Rtzmbdpf6BDNMYr3twhPQicq6rDVfVxIBNAVbfj9MC7wvfh5RORKOAZ4OsiikwEugJdVbUdMA+YLiKN/RmXMYEWExnOzf2bM+OBgdw/qA3xMRHM3XSAK8fPZeT4Ocy1xGRChDcJKUpVZxZ2QFWP4V0X8rK4E2fZi9/MFiEi5wGDgb+oaqa7ewwQDjzi57iMCQrVoyO4c+BpzHzwbO4+pxU13MQ0cvxcrnx9DnM2WmIywc2bhBQjIomFHRCRevhxHJKI1AHup+jkcjmQBaTm7VDVU8As95gxVUbNapHce15rUh88m3vPbU18TATzNh/gqjfmMuL1OczeYG1MJjh5NXUQMFNErhORlkCYiNQXkUHAFGCCXyJ0PAa8r6pbijjeCdjlJqGCNgP13YRpTJVSs1okd5/bitSHzuaP57WmZrVIftl8gKvfnMflr83mh9XWK88EFyntD6SIxOAkpfMBBcR9BafH3WWFJITyByhyGjADOF1VD4vI48BfgeZ5Ccrt+XdSVTt6nPsU8CjQTlV/M3BXREYBowASExO7T5jgz5waOo4ePUpcXFygwwgKlakujmcr323NYtqWLI5mOfuSaoRxQYtIzmgQTlgJHWUrU12Ul9VFvoEDBy5U1R6+uJY3C/SdAAaLyGDgXKAukA58p6rTfBFMEZ4DnlXVw2U4t9j/Yao6HhgP0KZNG01JSSnDt6h8pk+fjtWFo7LVxRDg2MlsPvxlG2/M3MT2IycZt/Qk3+yI5fazWjCsW2OiIwqfkqiy1UV5WF34h9cdEVT1G+AbP8TyGyIyAOhAybNApAONCtmf165lrbnGuKpHR3DLgBZc16cZkxbt5LXpG9myP5OHJi3nxe/XcWO/5lzdqynxMTaJq6lYXiUkEakF3AX0wUkAu4DZwFhVPeT78DgPp6fc/ALjbhu4r1+LyCmcjg7LgB4iEuXx2LA5kKaqe/0QmzEhLToinKvOaMrw7k2Ysnw3//5pI2vTMnh26hpe+XEDI3smcVP/5jSqVS3QoZoqwpuBsd2AjTgdDDrhJIrOwOPABhHp4uvgVPUxVW2pql3yNmCce3iou+9rYBIQCfQtEG+U+/5TX8dlTGUSER7GJV0a8809A3jnxp70aVGXoyezeTN1M2c+9xP3fryE1buPBDpMUwV4c4f0KvAV8Iiq7szbKSJNgL8B/6ZAQqhIqjpNRL4FxojIIHcs0qNArhubMaYEIsLANvUY2KYey3cc5vUZG/l6+W4mL97J5MU7aVc3jJz6aQxsU4+wsAqZKcxUMd4kpGRV7eO5U1V3iMgNwA6fRVUIERmKk1z+55Gde9cEMBx4Fli3roErAAAZT0lEQVQiIjluPCkFk6cxpnQ6NqnJK1d3Y/uBTN5K3czH87ezan8ON7+3gOS6sdzQN5kreiQRF+3v8fCmKvHmp2l7UQdUNVdEthbcJyLxquqz+3z30VxR0wahqhk4szkYY3wkqU4sj1/cnnvPbc3fPp5Oalo4W/Zn8viXq/jHtHWM6JnEDX2TSaoTG+hQTSXgzcDYD0XkMRH5nyQmIhEi8hfc7tMFTC9vcMaY4FAzNpIhzSP5+f4UXrumGz2Ta5NxMpu3Ujdz5vM/cfO78/lp7V5bl8mUizd3SEOB3sBdIrIBOALEA6cBOcBSd0bwPC19FqUxJihEhIcxpGNDhnRsyPIdh3l71mamLNvND2v28sOavTStE8s1vZoyokcStatHBTpcE2K8uUPqDSzBWao8C6jmvq4G1rnvC27W6mlMJdaxSU1evLILsx8+mwcGt6FxrWpsO5DJM1PX0OuZH7hvwlIWbTto0xOZUvPmDmmDqg4obWERWVyGeIwxISYhLprfp5zGbWe2ZPravfx37lZ+XrePTxft4NNFO2hTvwYjz0hiWNfG1Iq1uyZTNG8S0kVeXtvb8saYEBYeJpxzen3OOb0+W/cf44N52/hk4Q7WpmXwxJereGbqGoZ0aMDInk3p3aIOFbTItAkh3sxlV2y3bhGZpqrnl7a8Mabyala3Og8PPZ37zm/D96vT+PCXbaRuSOfzJbv4fMkukuvGMryHc9dkM0GYPN5OHXQmMBBnLJDnDIxdfRWUMaZyiIoIY2jHhgzt2JDtBzKZuGA7ExbsYMv+TJ7/di0vTFtLv5YJXN69MYPaNyA2ysY1VWWl/td3l314DDgOHCR/6Yk8flugzxgT+pLqxPLH89tw1zmtmLF+H58u3Ml3q9JI3ZBO6oZ0qketYGjHhlzevQlnJNex2SCqIG/+HLkZuFhVvyrsoHViMMaURkR4GGe3rc/ZbetzODOLL5ft4tNFO1i87RATF+5g4sIdNKwZw0WdG3Fx50a0bxRv7U1VhDcJaVtRych1fjHHjDHmN2rGRnJt72Zc27sZG/cdZdKiHXy2eBc7Dx1n/IxNjJ+xiRaJ1bnYTU4tEm1RvMrMm4T0HxEZqKo/FXH8VWCED2IyxlRBLRPjuH9QW+47rw2Lth3ki6W7mLJsN5v2HeOl79fz0vfr6dA4niEdnDap5gnVAx2y8TFvetm9LiIvi8gzwAYg06PIQJ9GZoypksLChB7JdeiRXIfHLmzHrI37+WLJLr5duYcVO4+wYucRnv92LW0b1HCTUwNa1bcm7MrA204No3E6NTTBOjUYY/wsIjyMs1onclbrRJ7O6sDM9elMXb6b71ansWZPBmv2ZPDi9+s4rV4cg9rX59zT69O5SS3rEBGivHlkdztwgapOLeygdWowxvhTTGQ457Wrz3nt6nMyO4fZG/YzdcVupq1KY8Peo2zYe5RXf9pIQlw0555ej3NOr0//0xKoFuU5QsUEK28S0taikpHLOjUYYypEdEQ4A9vWY2Dbejydk8u8TQf4fnUa361KY+eh43w0fzsfzd9OdEQYA1olkNKmHme1TrRlMoKcNwlpsoicrao/FnHcOjUYYypcZHgY/Vsl0L9VAn+9qB1r9mTw/ao0vl+dxtIdh/l+9V6+X70XgBaJ1TmrdSIpberRq3kdYiLt7imYeJOQWgKjRWQ71qnBGBOERITTG8ZzesN4Rp/TirQjJ/hxzV5mrNtH6vp0Nu07xqZ9x3hn1haiI8Lo1aIu/VrWpW/LBNo1iifc2p4CypuEdD2wF2jsbp6sU4MxJqjUj4/hqjOactUZTcnKyWXxtkP8vG4vP6/bx4qdR5ixbh8z1u0DID4mgt4t6tLvtAT6tqzLafXibEBuBfMmIa1S1SLnq7NODcaYYBYZHsYZzetwRvM63D+oLfsyTjJ7YzpzNu5n1sZ0th84zrRVaUxblQZA3epR9EiuTc/kOvRMrkO7RvFEhnuzhJzxljcJ6bYSjl9ZnkCMMaYiJdaI5pIujbmki/PAZ/uBTOZs3M/sjenM2riffRkn+XZlGt+udBJUtchwujatRc/kOoQfyqbzsVO2Kq6PeTMw9pe8r0WkEZCgqstEJExVc1V1nV8iNMaYCpBUJ5akOrGM6JmEqrJ1fybztxxg/pYDLNhykE3px5i9cT+zN+4H4J8Lv6NZ3Vi6JNWic5NadGlai3YN462jRDl4u/zESOBJnA4Oe3Dakt4XkV3A/WprFRtjKgERITmhOskJ1RneIwmAfRknWbj1APO3HGTGiq1sOwpb92eydX8mny/ZBUBEmHBavTjaNYqnXcN42jeqSbuG8dSMjQzkxwkZ3szUcBXwLvAF8AEw0j30IPAS8BDwjI/jM8aYoJBYI5rBHRoyuENDpsftpd+AM1mXlsGS7YdYuv0QS7cfZt3ejF9nkJjEzl/PbVK7Gm0bxNOqfhyt68fRql4NWibG2aBdD97cIT0InKuqMwFE5GIAVd0uItcDM7CEZIypIiLDw2jfqCbtG9Xkml7NAMg8lc2aPRms2nWElbuOsGr3EdbsPsKOg8fZcfA4369O+/V8EUiqHUurenE0d+/GkutWJzkhloY1q1XJLujeJKSovGTkSVWPiYgt9WiMqdJioyLo1rQ23ZrW/nVfdk4um9KPsS4tg3VpR9mwN4P1aUfZnH6MbQcy2XbAc0gnRIWH0bRuLM3qxNK4djUa16pGI3drUrsaiXHRlXK+Pm+SSIyIJKrqPs8DIlIPG4dkjDG/EREeRuv6NWjtMSP5qexctu4/xvq9TnLakn6Mrfsz2bz/GPsyTv46P19hIsOFejViqBcfTWJctPua/7529ShqxUZSOzaK+JgIIkKku7o3CekzYKaIPA3MBsJEpD7QBXgKmOCH+IwxplKKigijVf0ahS6dcexkNlv3Z7LtwDF2HjrBzoPH2XXoOLsOH2fnwePsP3aKnYeOs/PQ8VJ9rxoxEdSOjaJmtUiqR4dTPSqC2OgIqkeFExsVQVx0ONGR4URHhBEZnrcJURFhRISFUVH5zJuE9AjQDngPZ+kJAXa5x6YCj/k2NGOMqZqqR0c4PfUaxRd6/PipHPZmnGBfxkn2Zpx0X/PfH8zM4nDmKQ4dz+Lw8SwyTmSTcSK7gj+F97wZh3QCGCwig4FzgAQgHfhOVaf5KT5jjDEeqkWF06xudZrVLXnV3JxcJeNEFocyszh0PIvMk9kcO5VD5qlsjp10Xo+ezOZ4Vg7ZOUpWTi5ZObmcylb3NRf9zfJ3+cb78HMVm5BEJO+uZ4OqfgCgqt8A3/gwBmOMMX4SHibUio2iVqx/ZpUYf73vrlXSk8E7cB7NBaw7h4h0EZE3RGShiCwVkVUi8i8RSfQoFycir4jIWrfMNBFpH6i4jTHGeKekR3Z7VPWJComkaB8BK4Ez3e7ljYEfcB4fdlbVvFa9iUA80FVVM0VkDDBdRLqo6s7CL22MMSZYlHSHVOqpgETk4XLGUpwHVfUYgJtcngdaAUPd730eMBj4i6rmdeofA4TjdMYwxhgT5Eq6Q4oSkSRK98juGvwzU0MnVT3lsS+vd1/e6LPLgSwgNa+Aqp4SkVnusTv9EJcxxhgfKikhtQO2VEAcRSokGQG0xrl7m+G+7wTsKqTsZuBCEamnqnv9GKYxxphyKikhpQHjSnEdAUaVP5xSfCORcOAm4K0CS14kABmFFD/ivtbFWe3W81qjcONOTExk+vTpPo83FB09etTqwmV1kc/qIp/VhX/4rFODOz6pIvwFyAbuLUXZYh81qup43G70bdq00ZSUlHIHVxlMnz4dqwuH1UU+q4t8Vhf+4bNODarau5yxlEhEbgRGAENUteAkT+kUPpde3r79/o7NGGNM+ZSUkKJEJMlzzE8giMh1wH3A2YW0By0DGomI58iv5kCatR8ZY0zwKykhxQE/A89VQCxFEpFryV+PaY+770K3DQhgEhAJ9C1wTpT7/tMKDtcYY0wZFNuGpKrJFRRHkUTkGuANnLajc0V+bRYaAOwGUNVpIvItMEZEBrljkR4FcoG/VXzUxhhjvBUKi+qNBWJwBsN6KtjhYjjwLLBERHKAHUCKzdJgjDGhIegTkqrWKWW5DGwArDHGhKzQWEbQGGNMpWcJyRhjTFCwhGSMMSYoWEIyxhgTFCwhGWOMCQqWkIwxxgQFS0jGGGOCgiUkY4wxQcESkjHGmKBgCckYY0xQsIRkjDEmKFhCMsYYExQsIRljjAkKlpCMMcYEBUtIxhhjgoIlJGOMMUHBEpIxxpigYAnJGGNMULCEZIwxJihYQjLGGBMULCEZY4wJCpaQjDHGBAVLSMYYY4KCJSRjjDFBwRKSMcaYoGAJyRhjTFCwhGSMMSYoWEIyxhgTFCwhGWOMCQqWkIwxxgQFS0jGGGOCQqVKSCJST0T+T0TWutsnItIk0HEZY4wpWaVJSCISBXwHRAHtgXbAMeAnEYkLZGzGGGNKVmkSEvA7oBPwoKpmq2oO8CDQArgjoJEZY4wpUWVKSJcD21R1U94OVd0DrHKPGWOMCWKVKSF1AjYXsn8z0LGCYzHGGOOliEAH4EMJwMJC9h8BYkWkmqoeL3hAREYBo9y3J0VkhZ9jDBUJQHqggwgSVhf5rC7yWV3ka+OrC1WmhFQUKeqAqo4HxgOIyAJV7VFhUQUxq4t8Vhf5rC7yWV3kE5EFvrpWZXpklw7UKGR/DSDT8+7IGGNMcKlMCWkZkFzI/ubA8ooNxRhjjLcqU0KaBDQTkeS8HSJSHzgd+LQU54/3T1ghyeoin9VFPquLfFYX+XxWF6KqvrpWQLkDYxcAq4FrgFzgLaA/0FVVjwYwPGOMMSWoNHdIqnoKOA/IwRl7tBqIB862ZGSMMcGv0twhGVPRROQp4FHgRlV9N8DhGONTItIQeAcYpKpF9lb2pUpzh1SY8ky2KiKRIjJGRNaIyAoRmS0i/f0ds7+UtS5EpKGIPCEiv4jIYrc+JolIyA429sUkvG75P/opxApT3roQkc4i8rmILHJ/NtaKyHP+jNlfyvn7oqGIvOnWwTIRWSkij4hIpL/j9gcRGQbMAVqW8fx7RGSVWxeLROTSUp2oqpVyw5lkdSkwEWe8VTjwHrAeiCvF+eOAdUCi+/4W4DjQJdCfrSLrokA9JLnvY9zrZAIdA/3ZKvrnosB1/gN8BShwQ6A/VyDqAugL7AL6Fdh3J7Al0J+tIusC5w/7xcAKoK67r6v7++KFQH+2MtbHPKAV8K6TJrw69yGcYTgt3ffnAVnAkBLPDfQH92OF3ur+smhRYF8DnDam+0s4tw1Op4ibPPavBKYE+rNVcF2MA27x2NfSvd7YQH+2iqyLAuW7ARuBQSGekMrzcyE47bT3e+yPLM0vnmDbylkX7dxz7/XY/zmwO9CfrYz1EeG+epWQgFo4qyw86bF/CrCypPMr8yO78ky2OgznP9xPHvt/BM4PweUsylMXfwDe9ti3y32t7bMIK44vJuH9J07b0Unfh1ehylMX/YG2OHeJv1LVLFWd6utAK0B56iLbffWc+SbvTivkqGp2yaUKNRiIpfDfne1EpG1xJ1fmhFSeyVY74dwhbSvk3Aicv4hCSZnrQp2lPHI9drd2X6eXP7QKV65JeN1n4dWAj30cVyCUpy76uq813TaklW57wVMiUs2nUVaM8vwfWQd8ANyWNw5SRM7GeVQ11qdRBr9O7qtnXW72OF6oypyQEoCMQvb/OtlqCedmqrOmkue5AHV9EF9FKk9dFGYUzuPL/5Y3sAAoc124DdR/B+5T9zlEiCvPz0WS+/oh8LSqtgeuBW7AeVQVasr7f+R3wNfAehHZBXwG3KOqY3wbZtBLcF8967JUvzsrc0IqSnm6L1ZI18cK5PXncf/yuxIYoaqh/siqoNLUxR04z8FT/R1MgJWmLmLc17dU9RcAVV2Gk7DPE5Gz/BVcBSuxLkQkGucR1RlAsqo2AlKAh0XkUf+GFzJK9bumMiek8ky2mo7zV5Hn89+86+33QXwVyScTz4pIZ5zeZRer6iofxleRylQXIlILeBhnFeLKojw/F3l/AS/x2L/Yfe1ZztgqWnnq4macNrX7VXUngKouAl4AxohIF18HG8TyluTwrMtS/e6szAmpPJOtLsOpmySP/c1xGjBXlze4ClbuiWdFpBPOY4iRqjrbd6FVuLLWRW+cf/uJIrJERJYAb7rHnnT3PebTSP2vPD8Xa9xXz98hOUXsD3blqYu8Nqb1HvvX4dwZhFpyLo9l7muyx/7mHscLFWo/NN4o9WSrIlJfRArWxWScbpwpHtccCExT1cKeNQez8tRFXjL6HLgu73GVOxDwdT/H7Q9lqgtV/UZVk1S1S96GMzYN4DF335MV8gl8pzw/F1/jJB/PRuoO7ut8XwfrZ+Wpi73ua1OPazZzX0PtiUqpiUhddx7RPN/gjFFM8Sg6EFilqmsoTqD7u/uxH30UTjb+GKdnXBjONBj/M9AN6IfzH+s1j/PHAWuBBPf9TYT2wNgy1QXOX3/7gNdwGq3ztnuA6YH+bBX9c+FxrRRCexxSef+P/BPYDbRy3zd2z50W6M9WkXWB89f/EWAaUMPd1xTYgDNerVqgP1856uVdihiH5H7uE8BUj/0Pub8zWrjvz6WUA2Mr7YqxqnpKRM4DXsQZS6A4I6k9J1s9ChzG+Y9V0Gjgr8AsEcnCeWZ+vqp6PjMPeuWsiydwes7c7m4F/ey3oP3EBz8XiEg9nF8+eePRnhSRe3AGEPts9Ux/80Fd3I/TZvC1iOTgDIr9FOf/TUgpT12o6mYROQN4HJgvIqdw6uJbYIyG4OKgIvI8Trf1pu77vN97Z6gzkTU4f6AfIH9cIgCq+qyInAC+EpFsnAQ+XEsxPs0mVzXGGBMUKnMbkjHGmBBiCckYY0xQsIRkjDEmKFhCMsYYExQsIRljjAkKlpCMMcYEBUtIxhhjgoIlJGOMMUHBEpIxBhF5XEQ0b+JYEfGcWLgs10wqcD0Vkcd9EKqpxCrt1EGm6iowtU9TnGXWl+LMuhyNM8fWZ8DrHlPCICJ/B85R1R5efK9aOPP6fRaK00p5UmfSWF9dazvQBUBEbEoYUyK7QzKVjqrudX+xfuG+76KqnVW1LXA3cCmwVETaepy6l98uW1+SWjhzt1WlNW+M8QtLSKZKUWfhtHOAQziTgsYWOPYPVb0sYMEZU8VZQjJVjjtb8WM40+ffDCAir4rINretIzmvrIj0EZGfRWSxiCwVkakiMsw9NgxnXSDIX6RvifsYDxF5RER+EZGFIrJcRD4UkcYFrt3DLX9KRN4VkftFZI6I7BKR1z3WmUFEaonIayKyVUSWuduLItK8QJk6IvKGW2ad+/2HlLWu3PgOiMgWERkqIj+JyB4RmSwi8SLST0S+EZGdIjJRRGqW9XsZE/D1NmyzzV8bxa/lUg1nBdivC+y7AWfZgWT3fQ3gIHCN+16A5yiwDhTOypiFromEcxfWyf06HGdpg0VAuEe5LTjLGVzkvm/nxjaqQJkoYAEwA3d9HqAVzmPGe9z30e715wLx7r4r3GsNLKGuHi+mrt7FWXLhL+77+m69vI+zbDdAA7fMU0VcQ4HHA/0zYVtwb3aHZKokddaoSSd/Vc/CtMFpI9rsnqM4SeWTUn6b3qq6zD03B3gd6AoU1mkiTVW/dMuuwlkiPKXA8euA7sAj6nbGUNX1wHichJNXpivwZ1U94pb5BCeRlXeNojhgrHvNNCAVuAp4w923B5iJszKoMWVivexMVSYlHF8LpAGfichY4GNVXQe8Usrr1xaRz4DTcJJG3iO4FsA8j7LrPN4fwLkTyXOe+/o/S4Or6p8LvD0X505ktse1VgDXi0ikqmaVMnZP+1X1kEd8nvv24yz5bUyZ2B2SqZLczgx1ga1FlVHVDKAXMBn4E7DWbZPpV4rrdwZ+wlnCuos6vf6GuoejCzkl0+N9Ls5jvjwJQKaqnizm2ybgJqQC7VlLcO5aDuB0gS8rz/i0iH3hGFNGdodkqqrzcH55fl1cIVXdCtzmLlF+OfA0MFVEklX1QDGnXomTeJ5W1exiypVWOhArItHFJKV0nETW3X1EaExIsTskU+WISAzwJE7b0FvFlOsoIo+A0+akqu8D9+J0dkh2i+U9AhP3nO4i0pr8u6CCA0IblCPs79zX7h4xPigi9xUoE4HTKaJgma4i8no5vrcxFcISkqlSRKQH8CNOI/1Qt3NDUeoC94lIK/dcAfoCe4DVbpk04DjQxH3/MtAbmOK+/6N7bhTwYDlC/y+wEHhaRKq71+yAM0vENI8y/xCRGm6ZOjidEdaW43sbUzEC3c3PNtt8vQH1gCU47Sbqfr0Ep+daKk57UA2Pc17FmaVBgVXALThtMv8AlrvnrwS+Ajp6nHsbTlvUCpz2phh3/yhgA04y+MH9vup+n2eBlu51T7mxTnbPmw0cdbclQJS7vxYwzv1eS9zPcrZHLHlltuFMmbQAuLMUdfY4hXT7xmkHO+DGuAQnSU8uxb6WHtexbt+2lbiJqk0xZUxV5058+ldVLannYVmvr8ATqvq4P65vKgd7ZGeMAeduLM1fs33jPNo8WtI5pmqzOyRjjDFBwe6QjDHGBAVLSMYYY4KCJSRjjDFBwRKSMcaYoGAJyRhjTFCwhGSMMSYoWEIyxhgTFP4fg+e29HIPcPwAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Plot the temperature along the rod.\n", "pyplot.figure(figsize=(6.0, 4.0))\n", "pyplot.xlabel('Distance [m]')\n", "pyplot.ylabel('Temperature [C]')\n", "pyplot.grid()\n", "pyplot.plot(x, T, color='C0', linestyle='-', linewidth=2)\n", "pyplot.xlim(0.0, L)\n", "pyplot.ylim(0.0, 100.0);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Not too impressive, this looks just like the result from *explicit* forward in time, centered in space for $\\alpha\\frac{\\Delta t}{\\Delta x^2} = \\frac{1}{2}$. \n", "\n", "But try $\\alpha\\frac{\\Delta t}{\\Delta x^2} = 5$, which violates the stability condition of the *explicit* scheme:" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaQAAAEbCAYAAACV0PCVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xd4VGX2wPHvSSeFmkIn1NCLgEoPIorYdW1rL4sVy/7sfUXXtmtZdFexrK69ooiAoNIREaRX6WBooYcSUs7vj3tj4pg2yUxmJjmf57nPZO59782Zl5CTe98mqooxxhgTaGGBDsAYY4wBS0jGGGOChCUkY4wxQcESkjHGmKBgCckYY0xQsIRkjDEmKAQ0IYlIIxGZJCLW99wYY2q4gCUkETkX+AFoXUa5SBEZJSKrRGSZiMwRkf4llL1dRFaIyBIR+VlEzvFH7MYYY3wvkHdI9wJDgdlllBsNXAQMUNXOwJvAFBHpXrSQiNwLPAicqapdgXuAT0TkNJ9HbowxxuckUDM1iEiEquaKyFvAlaoqxZRJA1YC16nqm0X2Lwc2qurp7vu6wK/AP1X14SLlvgZSVbWTfz+NMcaYygrYHZKq5paj2LmAAFM99n8PnCIi8e77YUBsCeU6ikj7ysRqjDHG/4K9l11XIB/Y7LF/AxABdCxSrmC/Z7mix40xxgSpiEAHUIZE4LCq5nnsP+C+NihSDuBgGeV+R0RGACMAwmrV7pmUnEKd6D88Oaxx8vPzCQsL9r9VqobVRSGri0JWF4XWrFmTqapJvrhWsCekkpQ3a5RaTlXHAGMAohu11ZiLn2do10Y886euxEaFatVU3rRp00hPTw90GEHB6qKQ1UUhq4tCIrLJV9cK9hSfCcSKSLjH/gT3dXeRckX3l1SuRMmxQnx0BOOXbOP8//zAlj2HKxSwMcaYign2hLQEJ8ZmHvtbArk4PfAKygGkFlOu6PESxUYIY2/qS2qDWFZuO8DZL8/mh3Vl5jFjjDE+EuwJaSygQLrH/sHAZFUtaDOaBBwuodwKVV1Vnm/WNiWBL2/uz6B2Sew5dIzL3viRt+dsxBYxNMYY/wvqhKSqq3HaeO4TkUQAEbkGZ3aHB4qU2weMAm4WkVZuuZOBU4E7vfmedWIjefOq3lw/qBV5+coj45Zzz2dLyM717FdhjDHGlwLWci8iz+LM1NDcfb/IPXS8qh4rUnQk8AgwW0RycHrSnaKqi4peT1WfEpGjwHgRyQXygAtUdaK3sYWHCfed1oGOjWpzz2dL+Hj+Vn7ZmcWrl/UkuXaM15/VGGNM2QKWkFT1rnKWy8GZEujBcpR9AXihkqH95uzuTWidFM/17yxg4eZ9nDF6Fq9e3pMezev56lsYY4xxBfUju2DQuUkdvrylH8e3rM/Og9lc9OpcPp6/JdBhGWNMtWMJqRwS46N577oTuKJPC47l5XP3p0t4dNxycvLyAx2aMcZUG5aQyikyPIzHzu7M0+d3ITJceGvORq54Yx57Dh0r+2RjjDFlsoTkpYt6N+fDEX1ISojmh/W7OXP0LJZn7A90WMYYE/IsIVVAzxb1+OqW/nRrVpdf9x3h/P/MYdzijECHZYwxIc0SUgU1rBPDRyNO5IKeTTmak8+tHyzkyYkrycu3QbTGGFMRlpAqISYynGf+1JW/ndWJ8DDh1enrufqtn9h/OCfQoRljTMixhFRJIsKVfVN599oTqB8XxYw1uzjr5Vms2eG5EoYxxpjSWELykT6tGzDuln50alybTbsPc+7Ls5m0bHugwzLGmJBhCcmHmtaL5dMb+nJWt8YcOpbHDe8u4LnJq8m3diVjjCmTJSQfqxUVzosXd+eB4R0IE/jX92v5y//mc+CotSsZY0xpLCH5gYjwl4GtePua46lTK5LvVu3knJdns3ZnVqBDM8aYoGUJyY8GtE3iq1v6075hAut3HeKcl2fz7YodgQ7LGGOCkiUkP2veIJbPbuzL8C4NycrO5br/zefFb3+xdiVjjPFgCakKxEVH8PKfj+PuYWmIwPPfruGGdxdw0NqVjDHmN5aQqoiIcFN6G968qje1YyKYvGIH57w8m3W7rF3JGGPAElKVG5yWzLhb+tMuJZ51uw5xzkuz+W6ltSsZY4wlpABITYxj7E39OK1zQw5m53Lt29auZIwxlpACJC46gn9fehx3nVrYrnS9tSsZY2owS0gBJCLcPLiwXWmK265k45WMMTWRJaQgUNCulJaS4LQrvTybycttHjxjTM1iCSlIpCbG8flNfTm9SyOysnMZ8Y7Ng2eMqVksIQWRuOgIXvpzD+49rf1v8+Bd97/57D9i7UrGmOrPElKQERFuGNSat685nrqxkXy/aidnv2TrKxljqj9LSEGqYB68Do1qs3H3Yc55eTZfL9kW6LCMMcZvLCEFsWb1Y/n8xr6c070xh4/lcfP7P/PkhJXk5uUHOjRjjPE5S0hBrlZUOM9f1J1HzuxIeJjw6oz1XPHmPHZnZQc6NGOM8SlLSCFARLi6X0vev+4EEuOjmLNuN2eOnsWSrfsCHZoxxviMJaQQckKrBowfOYAezeuSsf8of3rlBz7+aUugwzLGGJ+whBRiGtaJ4cMRJ3LpCc05lpvP3Z8t4b7Pl5Kdmxfo0IwxplIsIYWg6Ihwnji3C8+c35WoiDA+mLeZC1+dS8a+I4EOzRhjKizoE5KI9BKRiSKyUkSWisg8EbnAo0ykiIwSkVUiskxE5ohI/0DFXFUu7N2MT2/oQ5O6tVi8ZR9njJ7FnLWZgQ7LGGMqJKgTkoikAt8BmUAXVe0CvAl8LCJnFik6GrgIGKCqnd0yU0Ske9VGXPW6Nq3LVyP7M6BtInsOHeOyN37klenrULUph4wxoSWoExIwHKgNPKequQCq+gpwAPgzgIikASOAp1R1l1vmdWA98EQggq5q9eOieOvq47llcBvyFZ6auIob3/3ZlrIwxoSUYE9Iue5rRMEOERGcuMPdXecCAkz1OPd74BQRifd3kMEgPEy489Q0XruiFwnREUxavp2zX57NLzblkDEmRAR7QvoQWAU8KCLxIhIG3A9EA6+4ZboC+cBmj3M34CSyjlUUa1AY2jGFcSOdpSzW7zrE2S/P5qvFGYEOyxhjyiQltTWISA5QkYaIbaraolJR/T6OxsB/gUFAFrAfuEZVp7vHJwN9VDXB47zrgNeA4ao6sYRrj8B53EdSUlLPjz/+2FdhB1x2rvLf5dnM3eZ0Bx/aIoKL0qKICJMyz83KyiI+vkbcWJbJ6qKQ1UUhq4tCgwcPXqCqvXxxrYhSjmUCr3t5PQEurHg4Hhdz2oe+AyYA9YGj7vU/F5HLSko0RWIplaqOAcYApKWlaXp6eqVjDianDFHembuJUeNXMGVTLntI4N+XHkdK7ZhSz5s2bRrVrS4qyuqikNVFIasL/ygtIW1X1Ye8vaCInF6JeDyNAuoCt6lqwSCbD0XkIuBt9+4pE4gVkXBVLTo6tOCOabcP4wkpIsIVfVLp1LgON7/3Mws27eX0f81k9CXH0ad1g0CHZ4wxv1NaG9KkCl6zoucVpwuwtUgyKrAGSAJaAktwPkczjzItcTpFrPRhPCGpZ4t6jL+1P31bNyAzy+ka/qp1DTfGBJkSE5Kq3leRC1b0vBLsBBqJiOedXAuc9q29wFj363SPMoOByapq3cyAxPho3rn2BG5Kb01evvLkxFVc/84CW43WGBM0KtXLzu315k+jccYhPeZ290ZEBgPnAR+paqaqrsZpB7pPRBLdMtcArYEH/BxfSAkPE+4e1t7pGh4TweQVOzjrpVksz9gf6NCMMab0hCQiTdypeuaJyKhiipwsIt+KSDt/BKeqnwLDgD7AChFZBryIk2iuKlJ0JPAJMNst8xfgFFVd5I+4Qt3Qjil8PXIAnRrXZtPuw5z77zl89JNnr3ljjKlaZd3hnAv0AKbgTMfjaTlwCJgmIk18HBsAqvqNqg5W1Q6q2llVu6rqs6qaXaRMjqo+qKppbpk+qjrTH/FUF80bxPLZjX255PhmHMvN557PlnLXJ4s5csxmDTfGBEZZCel0YKSqPqCqGzwPquqvqno28D/gbn8EaPwnJjKcJ8/ryj8u6EZMZBifLNjKuf+ezfZDtkS6MabqlZWQmuAMLi3LYzgDV00I+lPPpoy9qR8tE+NYtf0gj845woSl2wIdljGmhikrIeV7jO0plqoexpm+x4SoDo1qM+6Wfgzv0pCjeXDTez/z6LjlHMu1f1ZjTNXwZS+5SB9eywRAQkwkL//5OC7tEEVkuPDWnI1c8OoPbN17ONChGWNqgLIS0gYRKfNRnIgMBDb5JiQTSCLC0BaRfHx94cJ/p/9rFt+t3BHo0Iwx1VxZCelF4H0RGVBSAXdl1veAZ30ZmAmsHs3r8fWt/RnSPpn9R3K49u35PDVxFbl59gjPGOMfpc1lh6pOE5ExwHQRWQzMBXbgzIzQEDgRZ/mHewpm3zbVR93YKF67ohdjZq7n2W9W88r0dSzYtId/XdKDRnVqBTo8Y0w1U2Ybkqr+DTgfZ0G864GHgUfcrwU4S1X/4c8gTeCEhQk3DGrNB385kZTa0fy0cS/DX5zJ1NU7Ax2aMaaaKVenBlUdq6pdcbqB93O3pqraXVW/9meAJjgc37I+E24dwMB2Sew9nMPV//2JpyauIsce4RljfMSrXnaquk1Vf3A3W4a0hmkQH81bV/Xm7mFphIcJr0xfxyVj5pKxz3MydmOM8V6JCUlEHq3IBSt6ngkNYWHCTelt+HDEiTSsHcN8d42lqavsEZ4xpnJKu0M6q4LXrOh5JoT0Tq3PhNsGkJ7mPsJ76yee+HqFDaQ1xlRYab3sGorIfZRjKXAPcZWIx4SQ+nFRvHllb16dsZ5/TF7NazM3MG/jXkZf3IPmDWIDHZ4xJsSUmpCAJypwze0VjMWEoLAw4cb01hzfsj63frDQHUg7k6fO78rpXRsFOjxjTAgp7ZFdZAW3pn6M1wSpni3qMeHWAZzaKYWD2bnc/P7P3D92KUdzbDkLY0z5lLaEeV4FN2tEqKHqxEbyymU9eezsTkSFh/H+j5s5+6XZ/LLDVpE3xpTN30uQmxpGRLiiTypjb+5Lq8Q4Vu84yJkvzeLDeZtR1UCHZ4wJYpaQjF90alyHcSP7c95xTTiak8+9ny/llvcXsv9ITqBDM8YEKUtIxm/ioyN47sLuPH9RN+Kiwvl66TaGvziTBZv2BDo0Y0wQsoRk/O7cHk2ZcNsAujWtw6/7jnDhq3MZ/d0v5OXbIzxjTCFLSKZKtGgQxyc39OWGQa3Jy1f+OWUNf35tLtv227RDxhiH1wlJRPqKyP0i8oT7vr+I2GBYU6aoiDDuPa0971x7PEkJ0fy4YQ/DXpjJxKXbAh2aMSYIlDshiUiciEwEZgGPA1e7h84ClopIqs+jM9XSgLZJTLxtACe5i//d+N7P3P3pYg5l5wY6NGNMAHlzh/QUUA84E2gF7ARQ1btx1kh60ufRmWorMT6aN67sxWNndyI6IoyP52/ljNGzWLxlX6BDM8YEiDcJaTgwVFW/VtWNwG8DYFX1XaCdj2Mz1VzBmKWvRvanfcMENmQe4vz/zOHf09ZahwdjaiBvElKOqpY25L52ZYMxNVO7lAS+uLkfV/dLJTdfeWbSai59fS6/2jpLxtQo3iSkwyJydnEHRORUwJ61mAqLiQznkTM78dbVvUmMj2bu+j0Me2EGXy76NdChGWOqiDcJ6QngcxH5TkRGAYkicq+IfAB8CTzmlwhNjZKelsyk2wdwcocUDh7N5bYPF3HrBzbDgzE1QbkTkqp+BlwOpAEP4Mzq/XegP3Clqn7llwhNjZMYH81rV/TkyfO6UCsynHGLMzjthRnMWZcZ6NCMMX7k1TgkVX0faAZ0AQa7r81V9SM/xGZqMBHhkuObOzM8NKtLxv6jXPr6j/x9wkqyc21JC2OqI2/GIe0SkV+Blqq6XFWnu69+7w4lIueLyAwRWSAi60VkvohcXuR4pIiMEpFVIrJMROaISH9/x2X8r2ViHJ/e0IfbhrQlTIQxM9Zz9kuzWZFxINChGWN8zJs7pDCgl6qu91cwxRGRO3AeEf5ZVXviPDJcAwwpUmw0cBEwQFU7A28CU0Ske1XGavwjMjyMO4a245Mb+tCiQSyrth/k7JdnWfdwY6oZbxLSclUtcY4XETnDB/F4XjMVZ0Du9aq6FUBVc4A7gZfcMmnACOApVd3llnkdWE/FlmA3Qeq45s6qtJed2JycPKd7+IWv/sDGzEOBDs0Y4wPeJKTXReROESnpHH/0srsc2KeqPxXdqaoZqjrffXsuIMBUj3O/B04RkXg/xGUCJC46gsfP6cJbV/cmpXY0Czbt5bQXZ/Lu3E22AKAxIc6bhHQ58H/AdreNZnLRDWjth/j6AhvdNqSZbhvRHBG5pkiZrjizRmz2OHcDEAF09ENcJsDS05L55vaBnNWtMUdy8njwi2Vc9d+f2L7/aKBDM8ZUkJT3r0oROQgsKqVIN1X16WwNIrIMSAWW4twJ7QTOBz4AHlHVJ9xk2EdVEzzOvQ54DRiuqhNLuP4InMd9JCUl9fz44499GX7IysrKIj4+dG4sf9yWy/9WZHMoB2Ij4NIOUfRtHIGIVPraoVYX/mR1UcjqotDgwYMXqGovX1zLm4S0UFV7VPR4RYjIWpw7rwGqOqvI/s+AYUAS8AXFJ6S/AGMoJSEVlZaWpqtXr/Zl+CFr2rRppKenBzoMr+w4cJR7P1vC1NW7ABjaMYUnzu1MckJMpa4binXhL1YXhawuComIzxKSN4/szizj+IDKBFKCgrnzPO/MFgKxOI/jMoFYEQn3KFOQoHb7IS4TZFJqx/DmVb155k9dSYiOYMqKHZzy/AzGLc6wtiVjQoQ3MzVsLaPIyErGUpxV7qtnnHlF9i9xX5t5lGkJ5AIr/RCXCUIiwoW9mvHNHQMZ0DaRfYdzuPWDhdz8/s/szsoOdHjGmDJ4MzC2b2kbcJUf4iuYjqirx/7OwBFgOTAWUCDdo8xgYHIZM5Sbaqhx3Vr875rjefK8LsRFhTNh6XZOeX4G45fY3ZIxwSzCi7KzcH7xV6WPgNuBx0XkDFXNEpEBwJ+Ax1T1ELBaRMYA94nIeFXNdHvhtQYuq+J4TZAomHqof5tE7vlsCXPW7eaW9xcyvtM2HjunU6XblowxvudNQloH3OCxLw5oj7OM+Yu+CqqAquaJyDDgaWC5iBwFsoFbVPW1IkVHAo8As0UkB6ft6RRVLa1XoKkBmtWP5d1rT+CDnzbz5IRVTFq+nbkbdvPImR05p3sTn/TEM8b4hjcJ6QlV/a6Y/eNE5L/AC8AnvgmrkKruAf5SRpkc4EF3M+Z3wsKES09oQXpaMvd9vpQZa3Zxx0eLGb94G0+c24WGdexuyZhg4E2nhrdKObaLP7bzGBNUmtStxdtXuz3xYiL4btVOhj43nQ/nbba2JWOCgDedGhoXszURkc4i8hBgf2aaoFfQE2/KHYMY0j6Zg9m53Pv5Uv782o82J54xAebNI7utlNyp4ShwdeXDMaZqNKwTw+tX9mLc4gz+9tUKfli/m1NfmMEdQ9txXf+WRIR7tVSYMcYHvElIW/jjBKp5ONP5/KiqNgDVhBQR4ezuTRjQNolR41cwduGvPDVxFV8tzuDp87vSuUmdQIdoTI3iTUL6TFXf8FskxgRI/bgonr+oO2d3b8wDY5exPOMAZ788m+sGtKRnlLUtGVNVvHku8WVxO0XkBBF5RURa+igmYwIiPS2ZyXcM5Jp+LclX5dXp63lg1hFm/rIr0KEZUyN4k5CeL2H/dmA/8G7lwzEmsOKiI3j4zI58fmNf2jdMYNcR5fI35nHbhwvJtOmHjPErbxJSsSMIVXWTqt4D2Fzsptro0bweX43szwXtIomJDOPLRRkM+afTRTzflk03xi9KbUMSkTMpnOW7qTtFzx+KAU2BSB/HZkxARYaHcXqrKEaefTwPfrmMGWt2ce/nS/ns5638/dwutE1JKPsixphyK6tTQxvgNPfrhCJfF5WDszprqbMpGBOqmjeI5e2re/PVkm089tUKftq4l+H/msl1A1ox8qQ2xEZ50zfIGFOSUh/ZqerzqtpMVZsBKwu+9thaqeoQVZ1dRTEbU+VEhLO6Nea7vw7ikuObk5On/GfaOoY+N4PJy7cHOjxjqgVv2pDO8VsUxoSIOrGRPHleFz6/qS8dG9Xm131HGPHOAq57+ye27Dkc6PCMCWnezGW3qbTjIlLmMuHGVBfHNa/HuFv68ciZHUmIjuDblTsZ+vx0Xp66luzcvLIvYIz5A68efotII+BioBUQ7XHYJ2uqGxMqIsLDuLpfS07v0ognJqzky0UZPPvNaj5bsJVHzurEoHZJgQ7RmJBS7oQkIr2Bb4FjOB0cCkYLJuIkJ3uQbmqk5NoxvHhxDy7q1YwHv1zG+l2HuPLNeQztmMLDZ3SkWf3YQIdoTEjwpg3pKeA6VU2iSAcHnEX6HgBe8keAxoSKvm0SmXTbQO47rT1xUeFMWbGDIc9N57kpazhyzB7jGVMWbxJSoqr+YQE+Vc1X1SeBIb4Ly5jQFBURxvWDWvP9nemc070xx3Lz+dd3v3Dyc9OZtGybrbtkTCm8SUi/mzdFRBKKfB0NtPVVUMaEupTaMbxwcQ8+uaEPHdzeeDe8+zOXvfEjq7YfCHR4xgQlbxLSQRG5UkQEmA98LiKni8jpwFhgm18iNCaE9U6tz/iR/Rl1difq1Ipk9trdDH9xJg+MXcpumxvPmN/xJiG9AlwGtAZGua/jgK+APsD/+Tw6Y6qB8DDh8j6pTLsznav6piIivPfjZtL/MY3XZ67nWG5+oEM0Jih4Mw7pE1UdqqprVXUz0BU4AzgPaKuqs/wVpDHVQb24KB49qxOTbhvAwHZJHDyay+Nfr2TYCzP4buUOa18yNV65E5KIzHO3lgCqmqWqE1X1C1XN9F+IxlQvbVMSePvq3rx5VS9aJcaxPvMQ1749n8vfmMfyjP2BDs+YgPHmkV0H4GZV3eCvYIypKUSEk9qnMOn2gTx0Rkdqx0Qwa20mZ4yexZ2fLGb7/qOBDtGYKudNQlqiqj+VdFBEjvNBPMbUKFERYVzbvyUz7h7Mtf1bEhEmfLpgK+n/mMo/J68mKzs30CEaU2W8SUhfiMiFpRx/vbLBGFNT1Y2N4qEzOvLtXwcxvEtDjubkM/r7taQ/O433ftxEbp51fDDVnzdz2bUD7hCR+4CVQJbH8WY+i8qYGqpFgzj+fWlPFmzaw+Nfr2Th5n08MHYZb8zcwJ2npnFa54Y4Iy+MqX68SUhXADtx5q4bUMxxWz7TGB/p2aI+n9/YlwlLt/PsN6tYn3mIm977mW5N63DPsPb0bZMY6BCN8TlvEtIKVe1R0kERWeiDeIwxLhHh9K6NOKVTCh/9tIUXv/uFxVv38+fXf2RguyTuPjWNzk3qBDpMY3zGmzak68s4flFlAjHGFC8yPIzLTmzB9LvSufOUdiRERzBjzS7OGD2LkR8sZN0uz6fnxoQmbwbGziv4WkQai0hX9+sw9/ga34dnjCkQGxXBLSe1Zcbdg7muf0uiwsP4anEGQ5+bzl2fLLYVa03I8+YOCRG5WETWAFuAghVi3xWRf0gVtbSKyEwRURFJrYrvZ0ywqRcXxYNndGTaXelccnwzRIRPFmzlpH9O48EvltoYJhOyvJmp4RLgbWAxzlx2B91D9wAtgXt9Ht0fYzgf6F/CsXgReUlEVovIChGZLCKd/B2TMYHSuG4tnjyvK9/9dRDn9mhCbr7y7tzNDHp2Ko+PX0GmTd5qQow3d0j3ACer6gWq+ihwGEBVt+D0wPuT78MrJCJRwJPAhBKKfAL0AHqoakfgR2CaiDTxZ1zGBFpqYhzPX9SdybcPZHiXhmTn5vP6rA30f/p7Hh+/gp0H7Y7JhAZvElKUqs4s7oCqHsK7HnsVcTPOshd/mC1CRIYCw4CHVLXgQfooIBy4389xGRMU2qYk8O9LezJ+ZH9O7pDM0RwnMQ18ZiqjLDGZEOBNQooRkaTiDohIMn4chyQi9YG7KDm5nA/kAL/NOK6qx4DZ7jFjaozOTerw+pW9GT+yP0M7pnA0J583Zm1gwNNT+dtXy9l5wBKTCU5eTR0EzBSRy0WkNRAmIikicirwNfCxXyJ0PAy8q6obSzjeFchwk1BRG4AUN2EaU6N0blKH167oxde39ufUTilk5+bz39kb6f/MVB78Yqn1yjNBR8q7BouIxOAkpVMABcR9BafH3XnFJITKByjSBpgBdFDV/SLyKPAI0LIgQbk9/7JVtYvHuY8DDwAdVXVlMdceAYwASEpK6vnxx/7MqaEjKyuL+Pj4QIcRFKpTXWw+kMe4dTnM35EHQJjAiY0iOL1VJE3iy/7btDrVRWVZXRQaPHjwAlXt5YtrlbvdR1WPAsNEZBhwMtAAyASmqOpkXwRTgmeAp1S1IgvFlNoVXVXHAGMA0tLSND09vQLfovqZNm0aVheO6lYXVwC/7DjIf6av48tFGczJyGVORi6ndkrh5sFt6Nq0bonnVre6qAyrC//wuiOCqk4CJvkhlj8QkQFAZ8qeBSITaFzM/oJ2rd2+jMuYUNY2JYHnLuzOHSe3Y8yM9Xw0fwvfLN/BN8t30L9NIiMGtmJA20SbxNVUOa8SkojUBW4F+uAkgAxgDjBaVff5PjyG4vSU+6nIf46G7usEETmG09FhCdBLRKI8Hhu2BHao6k4/xGZMSGtWP5ZR53Rm5JA2vDFrA+/+sIlZazOZtTaT9g0TGDGwFWd0bUxUhFfj542pMG8Gxh4HrMPpYNAVJ1F0Ax4F1opId18Hp6oPq2prVe1esAGvuIeHu/smAJ8DkUDfIvFGue8/83VcxlQnyQkx3HdaB+bcO4S7h6WRlBDNqu0H+evHixn4zFTGzFjHwaM5gQ7T1ADe/OnzMjAeaKGqTVS1s6o2BlrgDFb9tz8CLA+3DesbYJSIxLq7HwDygb8HKi5jQkmd2EhuSm/DrHsG88z5XWmTHM/2A0f5+4RV9H3yez5YmW0984xfefPILlVV+3juVNWtInIVsNVnURVDRIbjJJffPbJz75oALgCeAhaJSJ4bT7qq/ury8adZAAAZRklEQVTPuIypbqIjwrmwdzP+1LMp09bs5NXp6/lxwx6+2QRTnp3KyR1SuKZ/S05oWd/amYxPeZOQtpR0QFXzRWRT0X0iUltVD1Q4sj9+jwmUPG0QqnoQZzYHY4wPhIUJJ7VP4aT2KSzdup8nP5/LTzvymLxiB5NX7KBDo9pc3S+Vs7o1JiYyPNDhmmrAm0d2H4jIwyLyuyQmIhEi8hBu9+kiplU2OGNMcOjStA4jukYz+96TuHVIWxLjo1i57QB3f7qEfk99zzOTVrF1rz3OM5XjzR3ScOBE4FYRWQscAGoDbYA8YLE7I3iB1j6L0hgTFJITYvjr0HbclN6arxZn8N/ZG1mx7QD/nraO/0xfx5D2yVx2YgsGtk0iLMwe5xnveJOQTgQWFXlfC2f+uJVF3hdlP43GVFMxkeFc0MtpZ1qwaS/vzN3ExKXb+XblTr5duZPm9WO57MTmXNCzGfXiogIdrgkR3iSktao6oLyFRWRhBeIxxoQQEaFXan16pdbnoTOy+Xj+Ft6bu5nNew7z9wmr+Mc3azi1c0Mu7t2MPq0a2F2TKZU3CelML6/tbXljTAhLjI/mpvQ2XD+wNdNW7+SduZuYvmYXXy3O4KvFGTSvH8tFbu+9lNoxgQ7XBCFv5rIrtVu3iExW1VPKW94YUz2FhwlDOqQwpEMKW/ce5pP5W/lk/hY27znMs9+s5rkpaxiclswFvZoyOC3ZZoIwv/F26qCBwGCcsUCe/Tx7+CooY0z10LReLHcMbcetQ9oy45ddfDRvC9+u3PHbVj8uirO6Neb845rSuUltG9dUw5U7IbnLPjwMHAH2Urj0RAG/LdBnjAlt4WHC4LRkBqcls+tgNp//vJXPft7Kmh1ZvDVnI2/N2Ui7lHjOP64p5/ZoQrI90quRvLlDuhY4S1XHF3fQOjEYY8ojKSGa6we1ZsTAVizPOMCnC7YybnEGa3Zk8eTEVTw9aRV9WjfgrG6NGdapEXViIwMdsqki3iSkzSUlI9cppRwzxpjfERE6N6lD5yZ1uH94B6av2cVnC7by/aqdzF67m9lrd/PQF8sZlJbEWd0ac3KHFGpF2YwQ1Zk3Cel/IjJYVaeWcPxl4EIfxGSMqWGiIsIY2jGFoR1T2H84h2+Wb+fLxb/yw7rdTFmxgykrdhAbFc7JHVIY3qUhg9olW3KqhrzpZfeqiLwoIk8CawHPeUIG+zQyY0yNVCc2kgt7N+PC3s3YefAoXy/ZxrjFGSzcvI9xizMYtziDWpHhDG6fxLDOjTipfTLx0V6vNWqCkLedGkbidGpoinVqMMb4WXJCDFf3a8nV/VqyZc9hJizdxoRl21m8ZR8Tlm5nwtLtREWEMbBtEqd2SuGk9sk0iI8OdNimgrz5s+IG4HRVnVjcQevUYIzxp2b1Y7l+UGuuH9SaX/cdYdKy7Uxato35m/b+1o1cBI5rXo+TO6QwtGMyrZPirSt5CPEmIW0qKRm5rFODMaZKNKlbi2v7t+Ta/i3ZeeAo3yzfzpSVO/lhXSYLNu1lwaa9PD1pFakNYhnSIYXBacn0Sq1ny2QEOW8S0lgROUlVvy/huHVqMMZUueTaMVzeJ5XL+6SSlZ3LzDW7mLJyB1NX7WTj7sO8MWsDb8zaQK3IcPq0bsCgdkkMapdEamJcoEM3HrxJSK2BkSKyBevUYIwJQvHREZzWpRGndWlEbl4+P2/ex3erdjBjTSYrtx3g+1U7+X7VTgCa149lYLtE+rVO5MRWDWxW8iDgTUK6AtgJNHE3T9apwRgTNCLCwzi+ZX2Ob1mf+06DHQeOMmPNLqav2cXMXzLZvOcw787dzLtzNyMCHRrWpm/rBvRt04DeqfVJiLEBuVXNm4S0QlVLnK/OOjUYY4JZSu0YLujVjAt6NSMvX1m0ZR9z1mYyZ91uFmzey4ptB1ix7QCvz9pAeJjQuXFtertLa/RKrUei9d7zO28S0vVlHL+oMoEYY0xVCQ8TeraoR88W9Rg5pC1Hc/L4edNe5qzbzZx1mSzeuv+37fVZGwBolRRH7xZOcsrJyic/X219Jx/zZmDsvIKvRaQxkKiqS0QkTFXzVXWNXyI0xhg/i4kMp2+bRPq2SQTSyMrOZdHmfczbuIf5G/ewcPM+1u86xPpdh/ho/hYAnvxpMl2b1aFb07p0b+ZsNils5Xi7/MTFwGM4HRy247QlvSsiGcBdquo5WNYYY0JOfHQE/dsm0r9tIgA5efmsyDjATxv3MH/jXuat28Geo7m/zblXoGHtGDo1rk3HxrXp2Kg2nRrXoVn9WjYWqpy8manhEuAtYBzwPnCxe+ge4AXgXuBJH8dnjDEBFxkeRrdmdenWrC7XDYBp06bR4bgTWbxlH4u27GPx1n0s2bKf7QeOsv3AUb5ze/IBJERH0KFRbdo3SqBtcjxtU5xXm1Hij7y5Q7oHOFlVZwKIyFkAqrpFRK4AZmAJyRhTQ6TUjuGUTg05pVNDAPLzlQ27D7Fy2wFWZBxgeYbTSWLXwWzmbdzDvI17fnd+/bgo2iTH0zY5npaJcbRMjKNFgzia1a9FdETNHMDrTUKKKkhGnlT1kIjY7IbGmBorLExonRRP66R4zuja+Lf9Ow8eZUXGAX7ZkcUvOw+yZkcWa3dmsefQMeZt2MO8Db9PVGECjevWIrVBHM0bxNKkbi2a1qtF47rOlpIQTUR49Vz23ZskEiMiSaq6y/OAiCRj45CMMeYPkhNiSE6LIT0t+bd9qsr2A0fdJJXFxsxDbNx9iE27D7N172G27j3C1r1HnCkIPISHCQ1rx9CwTgzJCdEkJUSTFO++ulu92CjqxUURFxUeUu1X3iSkL4CZIvIEMAcIE5EUoDvwOPCxH+IzxphqR0RoVKcWjerUYmC7pN8dO5abz9a9h9m0+zCbdh8iY/9Rft17hF/3HSFj3xF2Hszm133O+7JEhAl1YyOpUyuSerFR1KkVSVx0BHHR4cRGRRAXFU5sdARx0RHERIQRFRFGVHgYkeFhREaEERkuRIaHEVZFSc2bhHQ/0BF4G2fpCQEy3GMTgYd9G5oxxtQ8URFhtEqKp1VSfLHHs3Pz2LbvKDsOHGVXVja7Dmaz86DzWvD1/sPH2Hs4hyM5eWRmHSMz6xhwqGo/SAV4Mw7pKDBMRIYBQ4BEIBOYoqqT/RSfMcaYIqIjwklNjCvX5LBHc/I4cCSHfUdy2Hc4h32Hj3EkJ49D2XkcPpZLVnYuh4/lcSg7lyM5eeTmKTl5+eTk5ZOd67wey83/w+J3RW3y3UcrPSGJSMFdz1pVfR9AVScBk3wYgzHGGD+IiQwnJjLcrwN25WbfXausrho34jyaC1irmIh0F5HXRGSBiCwWkRUi8i8RSfIoFy8iL4nIarfMZBHpFKi4jTHGeKesR3bbVfVvVRJJyT4ElgMD3e7lTYDvcB4fdlPVgpa9T4DaQA9VPSwio4BpItJdVX8NTOjGGGPKq6w7pHJPBSQi91UyltLco6qHANzk8izQFhjufu+hwDDgIVUtWKdpFBCO0xnDGGNMkCvrDilKRJpRvkd2l+KfmRq6quoxj30Fvfvqua/nAznArIICqnpMRGa7x3z4lNMYY4w/lJWQOgIbqyCOEhWTjADa4dy9zXDfdwUyiim7AThDRJJVdSfGGGOCVlkJaQfwSjmuI8CIyodTjm8kEg5cA7xRZMmLROBgMcUPuK8NcFa79bzWCNy4k5KSmDZtms/jDUVZWVlWFy6ri0JWF4WsLvzDZ50a3PFJVeEhIBe4oxxlS33UqKpjgDEAaWlpmp6eXungqoNp06ZhdeGwuihkdVHI6sI/fNapQVVPrGQsZRKRq4ELgdNUNavIoUyKn0uvYN/uYo4ZY4wJImUlpCgRaeY55icQRORy4P+Ak4ppD1oCNBaRKI/9LYEd1n5kjDHBr6yEFA9MB56pglhKJCKXUbge03Z33xluGxDA50Ak0LfIOVHu+8+qOFxjjDEVUGobkqqmVlEcJRKRS4HXcNqOTi4ylfoAYBuAqk4WkW+AUSJyqjsW6QEgH/h71UdtjDHGW6GwqN5oIAZnMKynoh0uLgCeAhaJSB6wFUi3WRqMMSY0BH1CUtX65Sx3EBsAa4wxIat6roNrjDEm5FhCMsYYExQsIRljjAkKlpCMMcYEBUtIxhhjgoIlJGOMMUHBEpIxxpigYAnJGGNMULCEZIwxJihYQjLGGBMULCEZY4wJCpaQjDHGBAVLSMYYY4KCJSRjjDFBwRKSMcaYoGAJyRhjTFCwhGSMMSYoWEIyxhgTFCwhGWOMCQqWkIwxxgQFS0jGGGOCgiUkY4wxQcESkjHGmKBgCckYY0xQsIRkjDEmKFhCMsYYExQsIRljjAkKlpCMMcYEBUtIxhhjgoIlJGOMMUHBEpIxxpigUK0Skogki8h7IrLa3T4VkaaBjssYY0zZqk1CEpEoYAoQBXQCOgKHgKkiEh/I2IwxxpSt2iQk4EqgK3CPquaqah5wD9AKuDGgkRljjClTdUpI5wObVXV9wQ5V3Q6scI8ZY4wJYtUpIXUFNhSzfwPQpYpjMcYY46WIQAfgQ4nAgmL2HwBiRaSWqh4pekBERgAj3LfZIrLMzzGGikQgM9BBBAmri0JWF4WsLgql+epC1SkhlURKOqCqY4AxACIyX1V7VVlUQczqopDVRSGri0JWF4VEZL6vrlWdHtllAgnF7E8ADnveHRljjAku1SkhLQFSi9nfElhataEYY4zxVnVKSJ8DLUQktWCHiKQAHYDPynH+GP+EFZKsLgpZXRSyuihkdVHIZ3UhquqrawWUOzB2PrASuBTIB94A+gM9VDUrgOEZY4wpQ7W5Q1LVY8BQIA9n7NFKoDZwkiUjY4wJftXmDsmYqiYijwMPAFer6lsBDscYnxKRRsB/gVNVtcTeyr5Ube6QilOZyVZFJFJERonIKhFZJiJzRKS/v2P2l4rWhYg0EpG/icg8EVno1sfnIhKyg419MQmvW/6vfgqxylS2LkSkm4h8KSI/uz8bq0XkGX/G7C+V/H3RSERed+tgiYgsF5H7RSTS33H7g4icC/wAtK7g+beLyAq3Ln4WkXPKdaKqVssNZ5LVxcAnOOOtwoG3gV+A+HKc/wqwBkhy318HHAG6B/qzVWVdFKmHZu77GPc6h4Eugf5sVf1zUeQ6/wPGAwpcFejPFYi6APoCGUC/IvtuBjYG+rNVZV3g/GG/EFgGNHD39XB/X/wj0J+tgvXxI9AWeMtJE16dey/OMJzW7vuhQA5wWpnnBvqD+7FC/+L+smhVZF9DnDamu8o4Nw2nU8Q1HvuXA18H+rNVcV28Alznsa+1e73Rgf5sVVkXRcofB6wDTg3xhFSZnwvBaae9y2N/ZHl+8QTbVsm66Oiee4fH/i+BbYH+bBWsjwj31auEBNTFWWXhMY/9XwPLyzq/Oj+yq8xkq+fi/Ieb6rH/e+CUEFzOojJ1cQvwpse+DPe1ns8irDq+mIT3OZy2o2zfh1elKlMX/YH2OHeJv1HVHFWd6OtAq0Bl6iLXffWc+abgTivkqGpu2aWKNQyIpfjfnR1FpH1pJ1fnhFSZyVa74twhbS7m3Aicv4hCSYXrQp2lPPI9drdzX6dVPrQqV6lJeN1n4bWAj3wcVyBUpi76uq913Dak5W57weMiUsunUVaNyvwfWQO8D1xfMA5SRE7CeVQ12qdRBr+u7qtnXW7wOF6s6pyQEoGDxez/bbLVMs49rM6aSp7nAjTwQXxVqTJ1UZwROI8v36lsYAFQ4bpwG6ifBv5P3ecQIa4yPxfN3NcPgCdUtRNwGXAVzqOqUFPZ/yNXAhOAX0QkA/gCuF1VR/k2zKCX6L561mW5fndW54RUksp0X6ySro9VyOvP4/7ldxFwoaqG+iOrospTFzfiPAef5e9gAqw8dRHjvr6hqvMAVHUJTsIeKiKD/BVcFSuzLkQkGucR1fFAqqo2BtKB+0TkAf+GFzLK9bumOiekyky2monzV5Hn89+C6+32QXxVyScTz4pIN5zeZWep6gofxleVKlQXIlIXuA9nFeLqojI/FwV/AS/y2L/Qfe1dydiqWmXq4lqcNrW7VPVXAFX9GfgHMEpEuvs62CBWsCSHZ12W63dndU5IlZlsdQlO3TTz2N8SpwFzZWWDq2KVnnhWRLriPIa4WFXn+C60KlfRujgR59/+ExFZJCKLgNfdY4+5+x72aaT+V5mfi1Xuq+fvkLwS9ge7ytRFQRvTLx771+DcGYRacq6MJe5rqsf+lh7HixVqPzTeKPdkqyKSIiJF62IsTjfOdI9rDgYmq2pxz5qDWWXqoiAZfQlcXvC4yh0I+Kqf4/aHCtWFqk5S1Waq2r1gwxmbBvCwu++xKvkEvlOZn4sJOMnHs5G6s/v6k6+D9bPK1MVO97W5xzVbuK+h9kSl3ESkgTuPaIFJOGMU0z2KDgZWqOoqShPo/u5+7EcfhZONP8LpGReGMw3G7wa6Af1w/mP9x+P8V4DVQKL7/hpCe2BsheoC56+/XcB/cBqtC7bbgWmB/mxV/XPhca10QnscUmX/jzwHbAPauu+buOdODvRnq8q6wPnr/wAwGUhw9zUH1uKMV6sV6M9XiXp5ixLGIbmf+ygw0WP/ve7vjFbu+5Mp58DYartirKoeE5GhwPM4YwkUZyS152SrWcB+nP9YRY0EHgFmi0gOzjPzU1TV85l50KtkXfwNp+fMDe5W1HS/Be0nPvi5QESScX75FIxHe0xEbscZQOyz1TP9zQd1cRdOm8EEEcnDGRT7Gc7/m5BSmbpQ1Q0icjzwKPCTiBzDqYtvgFEagouDisizON3Wm7vvC37vHa/ORNbg/IG+h8JxiQCo6lMichQYLyK5OAn8Ai3H+DSbXNUYY0xQqM5tSMYYY0KIJSRjjDFBwRKSMcaYoGAJyRhjTFCwhGSMMSYoWEIyxhgTFCwhGWOMCQqWkIwxxgQFS0jGGETkURHRgoljRcRzYuGKXLNZkeupiDzqg1BNNVZtpw4yNVeRqX2a4yyzvhhn1uVonDm2vgBe9ZgSBhF5Ghiiqr28+F51ceb1+yIUp5XypM6ksb661hagO4CI2JQwpkx2h2SqHVXd6f5iHee+766q3VS1PXAbcA6wWETae5y6kz8uW1+Wujhzt9WkNW+M8QtLSKZGUWfhtCHAPpxJQWOLHPunqp4XsOCMqeEsIZkax52t+GGc6fOvBRCRl0Vks9vWkVpQVkT6iMh0EVkoIotFZKKInOseOxdnXSAoXKRvkfsYDxG5X0TmicgCEVkqIh+ISJMi1+7llj8mIm+JyF0i8oOIZIjIqx7rzCAidUXkPyKySUSWuNvzItKySJn6IvKaW2aN+/1Pq2hdufHtEZGNIjJcRKaKyHYRGSsitUWkn4hMEpFfReQTEalT0e9lTMDX27DNNn9tlL6WSy2cFWAnFNl3Fc6yA6nu+wRgL3Cp+16AZyiyDhTOypjFromEcxfW1f06HGdpg5+BcI9yG3GWMzjTfd/RjW1EkTJRwHxgBu76PEBbnMeMt7vvo93rzwVqu/v+5F5rcBl19WgpdfUWzpILD7nvU9x6eRdn2W6Ahm6Zx0u4hgKPBvpnwrbg3uwOydRI6qxRk0nhqp7FScNpI9rgnqM4SeXTcn6bE1V1iXtuHvAq0AMortPEDlX9yi27AmeJ8PQixy8HegL3q9sZQ1V/AcbgJJyCMj2AB1X1gFvmU5xEVtk1iuKB0e41dwCzgEuA19x924GZOCuDGlMh1svO1GRSxvHVwA7gCxEZDXykqmuAl8p5/Xoi8gXQBidpFDyCawX86FF2jcf7PTh3IgWGuq+/WxpcVR8s8vZknDuROR7XWgZcISKRqppTztg97VbVfR7xee7bjbPktzEVYndIpkZyOzM0ADaVVEZVDwInAGOBO4HVbptMv3JcvxswFWcJ6+7q9Pob7h6OLuaUwx7v83Ee8xVIBA6ranYp3zYRNyEVac9ahHPXsgenC3xFecanJewLx5gKsjskU1MNxfnlOaG0Qqq6CbjeXaL8fOAJYKKIpKrqnlJOvQgn8TyhqrmllCuvTCBWRKJLSUqZOImsp/uI0JiQYndIpsYRkRjgMZy2oTdKKddFRO4Hp81JVd8F7sDp7JDqFit4BCbuOT1FpB2Fd0FFB4Q2rETYU9zXnh4x3iMi/1ekTAROp4iiZXqIyKuV+N7GVAlLSKZGEZFewPc4jfTD3c4NJWkA/J+ItHXPFaAvsB1Y6ZbZARwBmrrvXwROBL523//VPTcKuKcSob8DLACeEJE495qdcWaJmOxR5p8ikuCWqY/TGWF1Jb63MVUj0N38bLPN1xuQDCzCaTdR9+tFOD3XZuG0ByV4nPMyziwNCqwArsNpk/knsNQ9fzkwHujice71OG1Ry3Dam2Lc/SOAtTjJ4Dv3+6r7fZ4CWrvXPebGOtY9bw6Q5W6LgCh3f13gFfd7LXI/y0kesRSU2YwzZdJ84OZy1NmjFNPtG6cdbI8b4yKcJD22HPtae1zHun3bVuYmqjbFlDE1nTvx6SOqWlbPw4peX4G/qeqj/ri+qR7skZ0xBpy7sR3+mu0b59FmVlnnmJrN7pCMMcYEBbtDMsYYExQsIRljjAkKlpCMMcYEBUtIxhhjgoIlJGOMMUHBEpIxxpigYAnJGGNMUPh/JVbOh2q88C8AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Increase the CFL number.\n", "sigma = 5.0\n", "dt = sigma * dx**2 / alpha # time-step size\n", "nt = 100 # number of time steps to compute\n", "\n", "# Compute the temperature along the rod.\n", "T = btcs_implicit(T0, nt, dt, dx, alpha, q)\n", "\n", "# Plot the temperature along the rod.\n", "pyplot.figure(figsize=(6.0, 4.0))\n", "pyplot.xlabel('Distance [m]')\n", "pyplot.ylabel('Temperature [C]')\n", "pyplot.grid()\n", "pyplot.plot(x, T, color='C0', linestyle='-', linewidth=2)\n", "pyplot.xlim(0.0, L)\n", "pyplot.ylim(0.0, 100.0);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**It didn't blow up!**\n", "\n", "We were not able to use such a large time step with the explicit scheme. You can try out other values of `sigma` and you'll get a stable solution. In fact, this is an *unconditionally stable* scheme—the most valuable feature of implicit methods is that they give stable solutions without a constraint on the choice of time step. \n", "\n", "Using the implicit scheme, we can always advance in time using larger time steps. But each time step requires the solution of a linear system, which is computationally expensive. This is the trade-off between explicit and implicit methods. \n", "To experiment further, set different values of the Neumann boundary flux and see if the solution behaves as you expect." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### A word of warning" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Implicit methods allow you to use significantly larger time steps, because they are not subject to stability constraints. But that doesn't mean you can use just _any_ large time step! Remember that Euler's method is a first-order method, so the _accuracy_ gets worse as you increase the time step, in direct proportion. In fact, you can lose the ability to capture the correct physics if your time step is too large. Numerical stability does not imply accuracy!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Dig deeper" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You see how matrix `A` is mostly full of zeros? We call such a matrix *sparse*, and there are many ways to make more efficient calculations taking advantage of their particular structure. First of all, you can optimize the memory usage. Check out SciPy's [sparse-matrix storage formats](https://docs.scipy.org/doc/scipy/reference/sparse.html): you don't need too store $(N-2)^2$ elements! For example, a `coo_matrix` format stores only $3*N_\\text{nonzero}$, where $N_\\text{nonzero}$ is the number of non-zero elements in `A`. Make sure to explore this topic a little more. It's an important topic in numerical PDEs." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "###### The cell below loads the style of the notebook" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from IPython.core.display import HTML\n", "css_file = '../../styles/numericalmoocstyle.css'\n", "HTML(open(css_file, 'r').read())" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (MOOC)", "language": "python", "name": "py36-mooc" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.6" } }, "nbformat": 4, "nbformat_minor": 1 }