{ "metadata": { "name": "", "signature": "sha256:be31e716b6db4cbb9e3b2c8b7b2641cd5e4730e3aa1a1514e5f95c4fab355502" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "Diffusion\n", "==\n", "\n", "$$\\frac{\\partial u}{\\partial t} = \\frac{\\partial}{\\partial x}D\\frac{\\partial u}{\\partial x}$$" ] }, { "cell_type": "code", "collapsed": false, "input": [ "%matplotlib inline\n", "from pylab import *\n", "from numpy import *" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 11 }, { "cell_type": "code", "collapsed": false, "input": [ "#dimensions of the computational domain:\n", "maxx = 10.\n", "maxt = 1.\n", "\n", "#difffusivity\n", "D = 0.5" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 12 }, { "cell_type": "markdown", "metadata": {}, "source": [ "We discretize the simulation domain:\n", "$$t_n = t_0 + n\\Delta t$$\n", "$$x_i = x_0 + i\\Delta x$$\n", "and seek for the numerical approximation $U^n_i$ to the actual solution $U(t_n, x_i)$ on the grid points" ] }, { "cell_type": "code", "collapsed": false, "input": [ "#discretization parameters\n", "nx = 100 # number of unknown grid points in spatial direction\n", "SC = 1.0 # stability criterion SC = 2*D*dt/dx^2, for FTCS should be SC <= 1\n", "\n", "def diffusion_init(maxx, maxt, D, nx, SC):\n", " #choose time step according to CFL condition\n", " dx = maxx/(nx+1)\n", " dt = SC*dx**2/(2*D)\n", " nt = int(maxt/dt)+1\n", " \n", " #define array for storing the solution\n", " U = zeros((nt, nx+2))\n", " \n", " x = arange(nx+2)*dx\n", " t = arange(nt)*dt\n", " return U, dx, dt, x, t" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 13 }, { "cell_type": "markdown", "metadata": {}, "source": [ "To obtain an explicit scheme for propagating the solution in time, we may replace the derivatives with forward differences in time and central differences in space (FTCS), assuming D=const.\n", "$$\\frac{U^{n+1}_i-U^n_i}{\\Delta t} = D\\frac{U^n_{i+1}-2U^n_i+U^n_{i-1}}{\\Delta x^2}$$\n", "This leads to a simple explicit formula for $U^{n+1}_i$\n", "$$U^{n+1}_i = U^n_i + \\frac{D\\Delta t}{\\Delta x^2}(U^n_{i+1}-2U^n_i+U^n_{i-1})$$" ] }, { "cell_type": "code", "collapsed": false, "input": [ "def diffusion_solve(U, dx, dt, nx, nt, D):\n", " xint = arange(1, nx+1)\n", " for it in range(0,nt-1):\n", " U[it+1,xint] = U[it,xint] + D*dt/dx**2*(U[it,xint+1] - 2*U[it, xint] + U[it, xint-1])" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 14 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's try to propagate a square initial condition" ] }, { "cell_type": "code", "collapsed": false, "input": [ "U, dx, dt, x, t = diffusion_init(maxx, maxt, D, nx, SC=0.9)\n", "U[0,:] = 0.\n", "U[0, (xmaxx*0.2)] = 1.\n", "plot(U[0,:])\n", "ylim([0,1.1])" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 15, "text": [ "(0, 1.1)" ] }, { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAD7CAYAAABgzo9kAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEgNJREFUeJzt3W+MZXddx/H3Z7ctESvUpglqu6YEimkbqFQtFVGuQmRp\nhBpJKBUIFJE+cBF9IKVIZDdGAwEikAqutTSEBPoACC5JS6nSG4kiUGxLobuwK1R3W1v+IxIIO92v\nD+5tuQy7c2dnf7Pzu2ffr2Sz9/yZM79fd/KZT8+955xUFZKk4dm00QOQJK0PA16SBsqAl6SBMuAl\naaAMeEkaKANekgbqpOP1jZL4eUxJWoOqylq+7rg2+Koa7J/Xv/71Gz4G5+fcnN/w/hwLT9FI0kAZ\n8JI0UAZ8I6PRaKOHsK6GPL8hzw2c34ksx3qOZ9XfKKnj9b0kaSiSUIvwJqsk6fgx4CVpoAx4SRoo\nA16SBsqAl6SBMuAlaaAMeEkaKANekgbKgJekgZob8EneleSBJHetsM/bk+xNcmeSJ7cdoiRpLVbT\n4K8Hth5pY5JLgMdX1TnAK4B3NhqbJOkYzA34qvo48M0Vdnku8O7pvp8ETkvymDbDkyStVYtz8GcC\n+2eWDwBnNTiuJOkYtHpk3/I7nZ2wt4183evg/vs3ehRtnXoqvOUtsHnzRo9E0tFoEfD3Altmls+a\nrvsx27dvf/j1aDQa5H2c3/QmeNvb4KTj9rTb9ffKV8KOHfDoR2/0SKThG4/HjMfjJsda1f3gk5wN\nfLiqnniYbZcA26rqkiQXA2+tqosPs98JcT/4k06C739/WAF/+umwb9/kb0nH17HcD35uDCV5H/B0\n4Iwk+4HXAycDVNXOqroxySVJ9gHfBa5Yy0CGoAoefHB4pzJOOgmWljZ6FJKO1tyAr6rLV7HPtjbD\nWWyHDsGmTZA1/a7t1+bNk19ckhaLV7I2tLQ0rFMzD7HBS4vJgG9oiKdnwAYvLSoDviEbvKSeGPAN\nLS0Nt8Eb8NLiMeAbevDB4TZ4T9FIi8eAb8gGL6knBnxDNnhJPTHgG7LBS+qJAd+QDV5STwz4hmzw\nknpiwDdkg5fUEwO+IRu8pJ4Y8A3Z4CX1xIBvyAYvqScGfEM2eEk9MeAbssFL6okB35ANXlJPDPiG\nvF2wpJ4Y8A35wA9JPTHgG7LBS+qJAd+QDV5STwz4hmzwknpiwDdkg5fUEwO+IRu8pJ4Y8A3Z4CX1\nxIBvyAYvqScGfEM2eEk9MeAbssFL6okB35ANXlJPDPiGbPCSemLAN+TtgiX1xIBvyNsFS+rJ3IBP\nsjXJniR7k1x1mO1nJPlIkjuSfC7JS9dlpAvABi+pJysGfJLNwDXAVuA84PIk5y7bbRtwe1X9IjAC\n3pJkgD12Phu8pJ7Ma/AXAfuq6p6qOgjcAFy6bJ//AR41ff0o4OtVdUL2PRu8pJ7M65tnAvtnlg8A\nT1m2z7XAx5LcB/wU8Px2w1ssNnhJPZkXR7WKY7wWuKOqRkkeB9yS5IKq+s7yHbdv3/7w69FoxGg0\nOoqh9m9pCR7xiI0eRXs2eOn4GY/HjMfjJseaF/D3AltmlrcwafGzngr8FUBV/WeSLwO/ANy2/GCz\nAT9ENnhJx2p5+d2xY8eajzXvHPxtwDlJzk5yCnAZsGvZPnuAZwIkeQyTcP/Smke0wDwHL6knK/bN\nqlpKsg24GdgMXFdVu5NcOd2+E/hr4PokdzL5hfHqqvrGOo+7SzZ4ST2ZG0dVdRNw07J1O2defw14\nTvuhLR4bvKSeeCVrQzZ4ST0x4BuywUvqiQHfkA1eUk8M+IZs8JJ6YsA3ZIOX1BMDviEf+CGpJwZ8\nQz6yT1JPDPiGbPCSemLAN2SDl9QTA74hG7yknhjwDdngJfXEgG/IBi+pJwZ8QzZ4ST0x4BuywUvq\niQHfkA1eUk8M+IZs8JJ6YsA3ZIOX1BMDviEbvKSeGPANebtgST0x4BvydsGSemLAN2SDl9QTA74h\nG7yknhjwDdngJfXEgG/IBi+pJwZ8Q0Nt8Js2QRUcOrTRI5F0NAz4hoba4MEWLy0iA76hoTZ48Dy8\ntIgM+IZs8JJ6YsA3ZIOX1BMDviEbvKSeGPCNPPQpk00D/S9qg5cWz0Dj6Ph76FbByUaPZH3Y4KXF\nMzfgk2xNsifJ3iRXHWGfUZLbk3wuybj5KBfAkM+/gw1eWkQrnjFOshm4BngmcC/w6SS7qmr3zD6n\nAX8LPKuqDiQ5Yz0H3Kshn38HG7y0iOY1+IuAfVV1T1UdBG4ALl22z+8DH6iqAwBV9bX2w+zfUB/2\n8RAf+iEtnnkBfyawf2b5wHTdrHOA05PcmuS2JC9uOcBFMdTH9T3Ex/ZJi2de56xVHONk4ELgGcAj\ngU8k+feq2rt8x+3btz/8ejQaMRqNVj3Q3tngJbUwHo8Zj8dNjjUvku4Ftswsb2HS4mftB75WVd8D\nvpfkX4ALgBUDfmhs8JJaWF5+d+zYseZjzTtFcxtwTpKzk5wCXAbsWrbPPwJPS7I5ySOBpwB3r3lE\nC8oGL6k3K0ZSVS0l2QbcDGwGrquq3UmunG7fWVV7knwE+CxwCLi2qk64gLfBS+rN3M5ZVTcBNy1b\nt3PZ8puBN7cd2mKxwUvqjVeyNmKDl9QbA74RG7yk3hjwjdjgJfXGgG/EBi+pNwZ8IzZ4Sb0x4Bux\nwUvqjQHfiLcLltQbA74RbxcsqTcGfCM2eEm9MeAbscFL6o0B34gNXlJvDPhGbPCSemPAN2KDl9Qb\nA74RG7yk3hjwjdjgJfXGgG/EBi+pNwZ8IzZ4Sb0x4BuxwUvqjQHfiA1eUm8M+EZs8JJ6Y8A3YoOX\n1BsDvhEbvKTeGPCN+MAPSb0x4BvxkX2SemPAN2KDl9QbA74RG7yk3hjwjdjgJfXGgG/EBi+pNwZ8\nIzZ4Sb0x4BuxwUvqjQHfiA1eUm/mBnySrUn2JNmb5KoV9vuVJEtJfq/tEBeDDV5Sb1YM+CSbgWuA\nrcB5wOVJzj3Cfm8EPgJkHcbZPRu8pN7Ma/AXAfuq6p6qOgjcAFx6mP1eCbwf+Grj8S0MG7yk3swL\n+DOB/TPLB6brHpbkTCah/87pqmo2ugVig5fUm3kBv5qwfivwmqoqJqdnTthTNENv8Aa8tFjmdc57\ngS0zy1uYtPhZvwTckATgDODZSQ5W1a7lB9u+ffvDr0ejEaPR6OhH3ClvFyyphfF4zHg8bnKsTIr3\nETYmJwFfAJ4B3Ad8Cri8qnYfYf/rgQ9X1QcPs61W+l6L7vnPh+c9Dy67bKNHsj5uvRV27IBGP3eS\nVikJVbWmMyMrds6qWkqyDbgZ2AxcV1W7k1w53b5zLd90iGzwknozN5Kq6ibgpmXrDhvsVXVFo3Et\nHM/BS+qNV7I2YoOX1BsDvhEbvKTeGPCN2OAl9caAb8QGL6k3BnwjNnhJvTHgG7HBS+qNAd+IDV5S\nbwz4RmzwknpjwDdig5fUGwO+EW8XLKk3BnwjPvBDUm8M+EZs8JJ6Y8A3YoOX1BsDvhEbvKTeGPCN\n2OAl9caAb2ToDd7PwUuLx4BvZOgNftMmSODQoY0eiaTVMuAbGXqDB8/DS4vGgG9k6A0ePA8vLRoD\nvhEbvKTeGPCN2OAl9caAb8QGL6k3BnwDhw5B1eSTJkNmg5cWy8Aj6fgY+q2CH2KDlxaLAd/A0B/2\n8RAvdpIWiwHfwInU4D1FIy0OA74BG7ykHhnwDdjgJfXIgG/ABi+pRwZ8AzZ4ST0y4BuwwUvqkQHf\ngA1eUo9WFfBJtibZk2RvkqsOs/2FSe5M8tkk/5rkSe2H2i8bvKQezQ34JJuBa4CtwHnA5UnOXbbb\nl4DfqKonAX8J/H3rgfbMBi+pR6tp8BcB+6rqnqo6CNwAXDq7Q1V9oqq+PV38JHBW22H2zQYvqUer\nCfgzgf0zywem647kD4Abj2VQi8YGL6lHq4mlWu3Bkvwm8DLg1w63ffv27Q+/Ho1GjEaj1R66azZ4\nSa2Mx2PG43GTY60m4O8Ftswsb2HS4n/E9I3Va4GtVfXNwx1oNuCHxAYvqZXl5XfHjh1rPtZqTtHc\nBpyT5OwkpwCXAbtmd0jy88AHgRdV1b41j2ZBnQgP+wBvFywtmrmxVFVLSbYBNwObgeuqaneSK6fb\ndwJ/Afw08M4kAAer6qL1G3ZfToTH9YEP/JAWzap6Z1XdBNy0bN3OmdcvB17edmiLwwYvqUdeydqA\nDV5Sjwz4BmzwknpkwDdgg5fUIwO+ARu8pB4Z8A3Y4CX1yIBvwAYvqUcGfAM2eEk9MuAbsMFL6pEB\n34ANXlKPDPgGbPCSemTAN+DtgiX1yIBvwNsFS+qRAd+ADV5Sjwz4BmzwknpkwDdgg5fUIwO+ARu8\npB4Z8A3Y4CX1yIBvwAYvqUcGfAM2eEk9MuAbsMFL6pEB34ANXlKPDPgGbPCSemTAN2CDl9QjA74B\nG7ykHhnwDdjgJfXIgG/ABi+pRwZ8Az7wQ1KPDPgGfGSfpB4Z8A3Y4CX1yIBvwAYvqUcGfAM2eEk9\nmhvwSbYm2ZNkb5KrjrDP26fb70zy5PbD7JsNXlKPVgz4JJuBa4CtwHnA5UnOXbbPJcDjq+oc4BXA\nO9dprN1aWoK77x5v9DDW1Xg8HmyDH4/HGz2EdeX8TlzzGvxFwL6quqeqDgI3AJcu2+e5wLsBquqT\nwGlJHtN8pB178EH4/OfHGz2MdTUejwfb4IceEM7vxDUv4M8E9s8sH5ium7fPWcc+tMWxtASbToB3\nM4ba4KWhmvfWYK3yOFnN1z3nOas82oK54w446wT4lXbyyfDlLw/v3/ELX4DPfGajR7F+TqT5veEN\ncP75GzuenqTqyBme5GJge1VtnS5fDRyqqjfO7PN3wLiqbpgu7wGeXlUPLDvWan9ZSJJmVNXyEr0q\n8xr8bcA5Sc4G7gMuAy5fts8uYBtww/QXwreWh/uxDFCStDYrBnxVLSXZBtwMbAauq6rdSa6cbt9Z\nVTcmuSTJPuC7wBXrPmpJ0lwrnqKRJC2udf/sx2oulFokSbYkuTXJ55N8LskfT9efnuSWJF9M8tEk\np230WI9Fks1Jbk/y4enyYOaX5LQk70+yO8ndSZ4ylPkluXr6s3lXkvcmecQizy3Ju5I8kOSumXVH\nnM90/nunmfPbGzPq1TvC/N40/dm8M8kHkzx6ZttRzW9dA341F0otoIPAn1bV+cDFwB9N5/Qa4Jaq\negLwz9PlRfYq4G5++ImoIc3vbcCNVXUu8CRgDwOY3/S9sj8ELqyqJzI5rfoCFntu1zPJj1mHnU+S\n85i8T3je9GvekaT3DzAfbn4fBc6vqguALwJXw9rmt96TX82FUgulqu6vqjumr/8P2M3kWoCHL/ia\n/v27GzPCY5fkLOAS4B/44UdgBzG/aRv69ap6F0zeZ6qqbzOM+f0vkwLyyCQnAY9k8uGIhZ1bVX0c\n+Oay1Ueaz6XA+6rqYFXdA+xjkkHdOtz8quqWqjo0XfwkP7yu6Kjnt94Bv5oLpRbWtDE9mck/wmNm\nPj30ALDIV/P+DfBnwKGZdUOZ32OBrya5Psl/JLk2yU8ygPlV1TeAtwD/zSTYv1VVtzCAuS1zpPn8\nHJOMecgQ8uZlwI3T10c9v/UO+MG+g5vkVOADwKuq6juz22ryzvVCzj3J7wBfqarb+fEL2IDFnh+T\nT45dCLyjqi5k8smvHzllsajzS/I44E+As5mEwalJXjS7z6LO7UhWMZ+FnWuSPwd+UFXvXWG3Fee3\n3gF/L7BlZnkLP/obaCElOZlJuL+nqj40Xf1Akp+Zbv9Z4CsbNb5j9FTguUm+DLwP+K0k72E48zsA\nHKiqT0+X388k8O8fwPx+Gfi3qvp6VS0BHwR+lWHMbdaRfhaX581Z03ULJ8lLmZwmfeHM6qOe33oH\n/MMXSiU5hckbBLvW+XuuqyQBrgPurqq3zmzaBbxk+volwIeWf+0iqKrXVtWWqnoskzfoPlZVL2Y4\n87sf2J/kCdNVzwQ+D3yYxZ/fHuDiJD8x/Tl9JpM3yocwt1lH+lncBbwgySlJHgucA3xqA8Z3TJJs\nZXKK9NKq+v7MpqOfX1Wt6x/g2cAXmLwhcPV6f7/jMJ+nMTk3fQdw+/TPVuB04J+YvOv9UeC0jR5r\ng7k+Hdg1fT2Y+QEXAJ8G7mTSch89lPkBr2byC+suJm9AnrzIc2Pyf5H3AT9g8n7eFSvNB3jtNGv2\nAM/a6PGvYX4vA/YC/zWTL+9Y6/y80EmSBqr3z4hKktbIgJekgTLgJWmgDHhJGigDXpIGyoCXpIEy\n4CVpoAx4SRqo/wexmrTFvAjGMwAAAABJRU5ErkJggg==\n", "text": [ "" ] } ], "prompt_number": 15 }, { "cell_type": "code", "collapsed": false, "input": [ "diffusion_solve(U, dx, dt, nx, len(t), D)\n", "pcolormesh(U, rasterized=True, vmin=-1, vmax=1)" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 16, "text": [ "" ] }, { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAEACAYAAABS29YJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XuMZNV9J/DvN90zY4anBisDDGQGJUYBK/YGZW3nxZQ2\nyJlYEeD8gbHW0dhG+cebXTt/eMPYkj0TKch2lNhZrfgja4MmVmCXOA4CCSnTIWlwlI0fMRDWAwsk\nOwSMpicLNu8wXTW//FG3qk51n9P3de6t+/h+pFLfvs9zblWf/tU5555DM4OIiHTbjyw6ASIiUj0V\n9iIiPaDCXkSkB1TYi4j0gAp7EZEeUGEvItIDWxb2JG8nuUbyMWfd75F8nOSjJL9O8nxn2yGST5F8\nguR7q0y4iIhklxbZ3wHgwIZ1xwC83czeCeBJAIcAgORVAD4A4KrkmNtI6puDiEgDbFkYm9k3APxg\nw7oVMzuT/PpNAJcmy9cDuMvM1s3sBICnAbwrbnJFRKSIspH3RwHcnyxfAuA5Z9tzAPaUPL+IiERQ\nuLAn+WkAp83szi1201gMIiINsFzkIJIfBvA+AL/krP4+gMuc3y9N1m08Vv8AREQKMDMWPTZ3ZE/y\nAIBPArjezP7V2XQvgJtIbid5OYC3AfiW7xxm1tnXZz/72YWnQflT/vqYvy7nzax8jLxlZE/yLgD7\nAbyV5LMAPotx75vtAFZIAsD/NrOPmdlxkncDOA5gCOBjFiOFIiJS2paFvZl90LP69i32vxXArWUT\nJSIicakffGSDwWDRSaiU8tduXc5fl/MWA+uuaSGp2h0RkZxIwupsoBURkfZRYS8i0gMq7EVEekCF\nvYhIDxR6grasz+DTAIAljKbrljAEAOzAaWfdbPsOvDld3p7ss91Z5x7nrt+JN5zl1wEA5+KV6bqz\nknXj9a/Olkezfc57MTn3S04m3OUXnOWXPetDx73mWf9mYPvQWXb3SbPDWXbf7bM928/3bN+4frJ8\nobPuPGf5Qs++zvLLu7ZPV72ydO5sGedMl9/AzmTdbPvryToAOI3t3vVvJutPO5ly9x1iabo8cm7G\nKFk/cra7+7pGnj+ZUWBfaY7fwe8uOgkLp8heRKQHVNiLiPTAQqpxtuJ+JXarcUL71MGSyxXu4Bri\ny97Qsy60b95ruO/25Do7AvtGZqrpEFkoRfYiIj3QuMg+j/nGstPB/TYKNb61ku+bQKvf1XR1v3++\nRlmRtlFkLyLSA60JWYYZ6vKlGzr1zUukIRTZi4j0QGsi+zSdfLAl1DMntN633dcDB5jvhdNDo7kH\nrDr42RHZQJG9iEgPdCayb4y0qDskrRkiT5QfiuZ9AWzR5o+i+ayQes2IhOmvo25phWtdbc9NSUdJ\nviqYULVM2j8DNQxLl6kaR0SkBzof2edpiKusoS4UJRetCinaQFv0vBVF+XnutxpRRcpRZC8i0gON\njuznB0XLHgbnGUah0Y16bkQd41tAo7PavnHim54+EZciexGRHmhwrLdYoyVP1NbEQM5Xn96UdHrS\n4b2vOeX5NqYeNiJjiuxFRHqgl5F9I+vphxt+du16OTTy/RFpOUX2IiI9sGVhT/J2kmskH3PW7SK5\nQvJJksdIXuBsO0TyKZJPkHxvVYkeYWnTq1dGzqvIdhHpnbTI/g4ABzasuwXAipldAeCB5HeQvArA\nBwBclRxzG0l9cxARaYAtC2Mz+waAH2xYfR2Ao8nyUQA3JMvXA7jLzNbN7ASApwG8K15StzbE0vTV\nGEOUqxMfBl5p++Q9R5E0ddwIy9OXSBcU+STvNrO1ZHkNwO5k+RIAf+fs9xyAPSXS1k3DwHITNDlt\nAUONSy+SSalqFjMzALbVLmXOLyIicRSJ7NdIXmRmJ0leDOBUsv77AC5z9rs0WbfJg4cfAgAQhn2D\nH8Plgx9Lvaj7dXo5R8ujG+0VreIZLU+uW4MsEXWegdBCWa45ch+VvHkaNE36ZnV1Faurq9HOV+RP\n8F4ABwF8Pvl5j7P+TpJ/gHH1zdsAfMt3gv2HrwGgicNFREIGgwEGg8H09yNHjpQ635aFPcm7AOwH\n8FaSzwL4DIDPAbib5M0ATgC4EQDM7DjJuwEcxzhu/FhSzdNojYkCfZG2uy7PMMnuuxoaTG2HZ13D\nVNnY3pj3XaQmWxb2ZvbBwKZrA/vfCuDWsokSEZG4FtqvrOgQxnnOW8dxqRocPWdSUfrzTCwTOq7o\n9UT6Rg89iYj0QGOeGJn0tsnT06b4tVoe4VU1LWFD1P3+5GkbaP1nR3pLkb2ISA+0MO7rEF+EnveL\nTdHJS3zHtb1twUPDHYiMKbIXEemBxoU9w7keOluHuVXWn9beMyct0o4x4bgvS3n670egOm+RxWhc\nYZ/F5B+C+89g/uv66U37Sje577v+kYiEqRpHRKQHWhnZF9WZxrrQEAg+Lc6y+37pG5pIOYrsRUR6\noNFx3yhHY23ouJj71iI0+FmMBtrlwPoGiPk+FB7Kutl/DiKlKLIXEekBhTKO2qP8GJF7WVmGUa6I\nBjQTqY8iexGRHuhlZF97ZBg7as/zsFWMrNY9haEid5HoFNmLiPRAayL7qiY6yXTtyYTjee9WkTrw\nPMMplElD0VuYJx3O/So74fgi6ZuGdIEiexGRHmhxvFWtyp7YzBMZV1nXvyO411hFPXMW+SSsInTp\ns1YX9lnmMO3kQFl5uku2+B2ef3/9GUn7DITecw2/IH2jahwRkR5ocdwXx0IjvDxdKGMMl+Cruokx\nZn4EirRFqqXIXkSkB1oZ2U/qb5cDFdZulLjde3yEKLKqO5cl0k6rp3e3hwY/qyqKj3Bf0t4ffQsQ\nyU+RvYhID7Qysq9KKyNGX4Tekne1qvutoYpFNlNkLyLSAwqB2qKq3jiibwLSC4Uje5KHSH6P5GMk\n7yS5g+QukisknyR5jOQFMRMrIiLFFCrsSe4D8BsArjazn8J4IN2bANwCYMXMrgDwQPL7Qo2wNH01\nwtB5uUbY3MtmFFgfOkfWa4xS1mW5hoi0StHI/mUA6wB2klwGsBPA8wCuA3A02ecogBtKp1BEREor\nVFlpZi+S/H0A/wzgDQB/YWYrJHeb2Vqy2xqA3ZHSmSpL5J5vIvLZrRkujf8n7sCZ9AN9kfBWEfjm\nC+dTtDeO7zp50pkz4p/cwzz141nGPqpKY74JikRSqLAn+eMAPgFgH4CXAPwpyQ+5+5iZkTTf8Q8e\nfmi6vHewF/sGe73XGc6NYb+5dErbnkVrGudij2ffMEXfB7dQbmXXWZGA1dVVrK6uRjtf0ZLuZwD8\nrZm9AAAkvw7gZwGcJHmRmZ0keTGAU76D9x++puBlRUT6YTAYYDAYTH8/cuRIqfMVrbN/AsB7SJ5F\nkgCuBXAcwH0ADib7HARwT6nUOfI0tI6wPH35DLE0feU5rlK+xtBQI+kwxyvtfAtsiA3db997k+U4\nEQkrWmf/KMk/BvAdAGcAfBfAHwE4F8DdJG8GcALAjZHSKSIiJRQOjczsCwC+sGH1ixhH+YW50Vrd\nc82mqiOQrCvLdVynAYG3GlpFxjRcgohIDzQg9spvEq01LvLPq475aBsw72xTqLeO9JkiexGRHmhl\nZJ9H4/ph54nQszzklGfC8RjXrtGw4ENVqqcX2UyRvYhID3Q+sk8TigKjRIdFIum6e+PkvV7J9FV6\nv0UkSJG9iEgPdCayzzNoVt4nL0fLk/1P501WPIF6+nVPpL0tdNwCJy+Z3cMM+6a8PzEGSNPTt9I3\niuxFRHqg1eGNG50tR+gkXlm9cVNGrGxKOqanjXu/y0brajeQLmtMYT97UKrjT/a48oxnH+OhqrTr\nNKC7pYhUQ9U4IiI90JjIvirzD+ZEzG7Rb/x5ZoMKcBtlhynHbUtLZ4T0zIlYEzI3W5iqWERKUWQv\nItIDnY/sfQp33QvdraruYuw69Krq5CPcl0XONyvSB4rsRUR6oDWRfYzJxfNYWHSZofeMW0+/nnIr\ntoXe4QX1vFHULrIYiuxFRHqgNZG9y9cnfz7yL/oIfWCQrqWKotE8QxVXxb1uhdG+7x4WH+pg62Gr\nGzestUgDKLIXEemBRkf2sScfb119sVs3XzD77nHbWvaEbOveL5EGU2QvItIDjY7s04zm6un9Yev8\nU7ObhyhuXPSYIfpeD/TM8QmOLNywKD/GsNRFz9G4z4BIBVpd2BdV+x930YbYRT5UVVPDrf/SKnxF\nYlM1johID/Qysg8Z1v3Ifo7oeRgY/Cyt4XboZGObb6aqmrt/qlukyGIoshcR6YFORvZVReXmnJaV\nXAH1PVxVw3WsosBddfoi+RWO7EleQPJrJB8neZzku0nuIrlC8kmSx0heEDOxIiJSTJlqnD8EcL+Z\nXQngHQCeAHALgBUzuwLAA8nvtRhhefoK77O0KSocYmn6ynNcLUaz17r7Gs5eQ2x+udvd49zz1Z+V\nre+h732YHLPVcVnedxEpWNiTPB/AL5rZ7QBgZkMzewnAdQCOJrsdBXBDlFSKiEgpRcOhywH8C8k7\nALwTwN8D+ASA3Wa2luyzBmB33hPPPyhVPgT1RYVReoHEDiQX9ZBT7OtGuC9pg5uJSH5F/zSXAVwN\n4DfN7Nskv4QNVTZmZiTNd/CDhx+aLu8d7MW+wd6CyRAR6abV1VWsrq5GO1/Rwv45AM+Z2beT378G\n4BCAkyQvMrOTJC8GcMp38P7D1xS87Ngk8os9iUn0Ca59UfMw+/bQxOLuYetpSQiczzsoWlratlqf\nUWUTwAeuIdJWg8EAg8Fg+vuRI0dKna9Qnb2ZnQTwLMkrklXXAvgegPsAHEzWHQRwT6nUiYhIFGVC\nq/8M4E9IbgfwjwA+AmAJwN0kbwZwAsCNpVOY8A13HLt+X9otT72+2gCkbwoX9mb2KIB/79l0bfHk\nxOF+jd8eWD+RpSqhlmqBgv+rfNU422q4bl5p99D3PoSOyfN+qEumyJiGSxAR6YHOhD2hse19X9eL\nRnsj57BcZ8jTqJlhwLO0Rll3+3LgfLlmrSrYKDsq+OnyvT/B+YFVHSOSiSJ7EZEe6ExkH8Oohm6B\ncxo2W1Qd6XHvq6JykfooshcR6QFF9nVLe1jJrZsfbd4MzNfJpwXjbs8c93xnea7nTY+IdIIiexGR\nHuhVZJ+nl0facXOrenUXPTLci7Q2kJi9pkRkM0X2IiI90OjQqehwCPM9PjYf18Y+26F6+jx97psm\nz/sQo6dUk99fkaopshcR6YFGR/ZphhUOhOaNLos+QRviG1LY7YETuUfM3Pkm18kyrHFBvidoq4yu\nNbSxSFgrC/vJ1/jQsAj5qnzaUUCsB5bTyue2vMFF3wc17Ipko2ocEZEe6FUINIkC80aRpaP/UWDZ\nYz0w+FkMvnMHh0POkeYsiow135ZvXSJtoMheRKQHehXZp5nv3lfzhCUp0XOoHTV1DtqIaYih9nss\nIgAU2YuI9EInI3s3YgxNUejbnma4NPvfuGPpTKG05RF7LLJaxjZzbqd7v9KkvQ9DfSMQKUWRvYhI\nD7Qmsi/aj953jkofvil5R0MPUsXuZz9M642TpsJPzjBibxx9CxAZU2QvItIDrYns6zAM1PVHjQ5D\n3Wo8IXrsQcy858uRnjJ836o0vIFIfRTZi4j0QKsj+9AYKEWHQy4sYoC6Hkh67H72k+uclS1ZW4uQ\n/6Lvw/w3sFZ/nEUq1fm/jvkC4HS5cy3HO9cmBatxfIV5qNE1tRonsvn7VfJc3f+oilRK1TgiIj3Q\nyXAp1E1zNsCWP9vzM1y5y556iix3zjdmvMuzPm/Xy7R907pepqVn0/q0GrLgHLSb733a+9DGGcVE\nmqpUZE9yieTDJO9Lft9FcoXkkySPkbwgTjJFRKSMstU4HwdwHIAlv98CYMXMrgDwQPJ7o42wNH3V\nYui8Iuy67nlFuHTOncur/X0Q6ZnChT3JSwG8D8CXATBZfR2Ao8nyUQA3lEqdiIhEUabO/osAPgng\nPGfdbjNbS5bXAOwucf6gSfRXx7yzc9uXKow6UyLoPM875Xo2qsreOCn3q6ooXt8ORDYrVNiT/FUA\np8zsYZID3z5mZiTNt+3Bww9Nl/cO9mLfYG+RZIiIdNbq6ipWV1ejna9oZP9zAK4j+T4AbwFwHsmv\nAlgjeZGZnSR5MYBTvoP3H76m4GXz8w13nCXyCw2THJXni8nc1IEZTlG0n713ysMKJy8pcu81rLH0\n2WAwwGAwmP5+5MiRUucrVGdvZp8ys8vM7HIANwH4KzP7dQD3AjiY7HYQwD2lUiciIlHE6mc/qa75\nHIC7Sd4M4ASAGyOdvzKFI8bYTyikTUQeWE7bNzWZsaP5gvdFkbtItUoXWWb2IIAHk+UXAVxb9pwi\nIhJXZ56gjTMomr+OOO2Jzlx3seDE4qF98gyDnKuXTp6IP+Wp2fHy5ntYNJrX4Gci+fXyL6X2KoMM\npba3wTTDKYpe2nfdbTmGdaiSqnRE4tNAaCIiPdDqyD7LvLS+qgR3eOKmzZyUpSE2T2OtO159nuPq\nkHbvUwej27T/4t8/kaZSZC8i0gOtjuyLyhsl+h7yKXznCtZ/F43EC0fwRevpnfvieyhKEbrIYiiy\nFxHpgV5G9q4sk2lUxomeJ5OJrPs3z8myz4RbZ+/rsjnM0hunIgu99yI9o8heRKQHOhlOuVHismda\nwixyDcJVw12sqp99FBny7233yHjM5vWd/NiKVEqRvYhID3Q+RBp6+uIXHTrXjSjNOYyefTckwr/s\nPBqwPtq8OdSTpuhwCb66/nUnDWeNPDtsXE7h3pc8EXjakNJNeAZCpM0U2YuI9EBnIvssT9PmOYcb\nldYZVWZ5ataVNnlJ3vNVZTZ5Sb6nYtOoT75INp0p7LMoWzBEL1g8XS8Dm6OMZ+/7xzAsWF2TRePu\nt0iPqRpHRKQHehXZT4SqEtLGs587h7NqbmtqKO1P02So4byNr759Qm+qt4E277DGk/WhMeyDY9tv\nPZ592v0WkXIU2YuI9IBCqBwKz1QVai/OUUeep2q9cDV8aOe09u4MM1WJyGIpshcR6YHOR/ahoRNm\n2/0P86TNlRqss09PkNe0Dt2zbitpXS9D+6571s0p1nt1izr7zUMc53m4TfX4IuUoshcR6YFOhkuh\nB6wmkeRShgG2/NPkVdfPfotVAOL0s0+9Tg397IcZet2kDZqmtgCR/BTZi4j0QCcj+zyyDIXgjVCX\nZv8ndyydSbuIf9nhi9ar7Gefeo08XXqc2+PeF1daNK6BzkSqpcheRKQHehnZhyfF8NcnF5l4I2NC\nptY3/Ny4HKOfve/cc5F9wR444XT4JhxPHwhNdfIi8SmyFxHpgUKRPcnLAPwxgB8FYAD+yMz+G8ld\nAP4XgL0ATgC40cx+GCmtlUjrERLad7TsHnc6+wWdUDs4Ls2CuOkpOvm4e1+8YwppqkGRhSj6F7YO\n4LfM7BGS5wD4e5IrAD4CYMXMvkDytwHckrwWxle4+B6uCu0LZKi+SbuLGQrO4YafG5ezVO/4hMZl\n810vNXFZLhI8Rf6urLHHvhfps0LVOGZ20sweSZZfBfA4gD0ArgNwNNntKIAbYiRSRETKKf3dmeQ+\nAD8N4JsAdpvZWrJpDcDusuevwvy8tNkba+fWLQUiTXe1LyJ2vlS4E4d4G0wr5Luem55taYO3BbLv\n3pe0ezh3XMoctCJSTqkG2qQK588AfNzMXnG3mZlhXJ8vIiILVjiyJ7kN44L+q2Z2T7J6jeRFZnaS\n5MUATvmOffDwQ9PlvYO92DfYWzQZUaUN0hV7iON1N8rf8BPI1vUy7ZtA2vnmtjvpOcs9SYQhjn3d\nV1UPLxK2urqK1dXVaOcr2huHAL4C4LiZfcnZdC+AgwA+n/y8x3M49h++pshlRUR6YzAYYDAYTH8/\ncuRIqfMVjex/HsCHAPwDyYeTdYcAfA7A3SRvRtL1slTqapClrjhqVBqYZNwXocfumek7X6jOvujF\n8wxhrChfpD6FCnsz+xuE6/uvLZ4cERGpQq+eZMnX5973QFBgCOS0oDQQMfseqgrVwccY4viswD7e\n9BQdCC1tovYMH7m0njsikp+GSxAR6YFeRfZpQk9s+p/+LHjrckxLmKXPvS/ozjLEcWXTEgauPvT0\naNIQCSL1UWQvItIDvQ+tQk/T+iLQ09juP0nBsXHecJZ9kXaEzjGp53OjfTc9wZP4BPLv3i/vU8h6\nalakNr0s7OfnqPWXZKFCaSJTwe879b/6D/MVviF5xrPfFtgntRonkM7pARkKeFfqg2mefUUkHlXj\niIj0QC8je9d8o+ysVdL3wM+beaP5yfpAY2daNU7exlqfUDfMXNU4bvonJ/Hlc4M356pxtn4wTY21\nItVSZC8i0gMKpxyhRsI3sQPAfPQZjPLfMneg72RTaVF1HcMl+L5dAPCn3V3/ltDmzY2y4/U7POlR\n3bxIXRTZi4j0QO8j+1DPnNNzkeiryTp/NG/OrnTD58mpAz1bfNF8qN48RtdLt2dOnmvPpX+SJ+fE\ntjloB+C/X6c9ET6gHjgiVVNkLyLSA72P7F3zPUJOO0ubI9Q3sNN/Ejdw9dR7v/ya/zBfVB17ikL3\nfJNB0ULRvJvO89wNkzwFovnQffHdQ/XAEamPInsRkR5QaBXg9irZmdQnvx6IWl8/e/Y/8+w3z8w2\nTKLfQDTvTto7ibqD9eaRTa7jRvuv+HYE5tM/yZPzrcXNv8u9X6nPKohIpRTZi4j0gCJ7R6hnjq+P\n+Cs4d7q8PHKi+bOdnSYRsXOXX8zRM6fKfvaTnjmhbxJuOs9zPyWTPDn5dPPv3heX7x6qB45IfVTY\nB/gaD92CbMkZQ+CVnedMl3e89ursgEmB+NJsldv98QXPdWM3yob4ruOm50p3g5P+aZ6cqh03/3P3\nxVPwq1FWZDFUjSMi0gMKszJ4PemoeK7ThPkCLpwubx/NumniR50DTyU/nRqM5wLXCDaOVsy97i5n\n2U3nHrcGZpInJ59u/l9YuhA+r8/NgCsidVNkLyLSA4rsA3yNtT/EBdN1l+D56fLzSxdPl8978ZnZ\nSSZB7rOzVW6d/SlsFrtRNmRynVB63ubu7KQfv5z8fHG2ys2/W2fv3q8JNcqKLIYiexGRHlBkn8Fk\n8K6dTkfF53HJdPnf4ZHZzpc7B96R/Dx/tuopp8uLe/PrephqI/e6bt+Zp5zlq53045+Snx92j5v1\nQHLvixvlhwZAE5F6KLIXEekBRfY5uHXQF+L/T5ePO73S9zz1N7MDkorvx/90tsodcOGJ2Aksya2z\nv9RZfvyfZstX3pQsOKH/8Z+c5X+7M47CD3FZ1PSJSHGK7EVEeiB6ZE/yAIAvYTzNxZfN7POxr9EE\n/4ifmC6/A49Nl197rzMo2u+MhxHY40zhd3dguISm+a6zvN+dgvDY+Mdrn5nl863Os7ffwDXVJkxE\nCoka2ZNcAvDfARwAcBWAD5K8cuujuuWhB23RSajU6suLTkG1Tqw+k75Ti3U5f13OWwyxI/t3AXja\nzE4AAMn/CeB6AI9Hvk6j/DneP13+ka/cil95LSnw3zf+8cXvLCBREX3R+TZiFwGDAXD2/bPBz/78\n/e/ffFBLPbP6DPYN9i46GZXpcv66nLcYYhf2ezD/CM5zAN4d+RqN9tdX/gLOvH9clbH0a7cuODXx\nrT4EHHkIGNmnFp0UEckhdgNtt+swRERaimbxymeS7wFw2MwOJL8fAnDGbaQlqX8IIiIFmBmLHhu7\nsF8G8H8B/BKA5wF8C8AHzazTdfYiIk0Xtc7ezIYkfxPAX2Dc9fIrKuhFRBYvamQvIiLNVOsTtCQP\nkHyC5FMkf7vOa1eB5GUk/5rk90j+H5L/JVm/i+QKySdJHiO5eazfliC5RPJhkvclv3cpbxeQ/BrJ\nx0keJ/nujuXvUPLZfIzknSR3tDl/JG8nuUbyMWddMD9J/p9Kypz3LibV2QXy93vJ5/NRkl8neb6z\nLVf+aivsO/rA1TqA3zKztwN4D4D/lOTpFgArZnYFgAeS39vq4wCOY9bTqkt5+0MA95vZlQDegfFw\nRZ3IH8l9AH4DwNVm9lMYV6vehHbn7w6Myw+XNz8krwLwAYzLmgMAbiPZ9OFhfPk7BuDtZvZOAE8C\nOAQUy1+dmZ8+cGVm6wAmD1y1lpmdNLNHkuVXMX54bA+A6wAcTXY7CuCGxaSwHJKXYvxo2JcBTHoB\ndCVv5wP4RTO7HRi3N5nZS+hI/gC8jHEwsjPpOLET404Trc2fmX0DwA82rA7l53oAd5nZevKQ59MY\nl0GN5cufma2Y2eQJxm9iNkZh7vzVWdj7HrjaU+P1K5VEUj+N8Ruy28zWkk1rAHYvKFllfRHAJwGc\ncdZ1JW+XA/gXkneQ/C7J/0HybHQkf2b2IoDfB/DPGBfyPzSzFXQkf45Qfi7B/FTKXShvPgrg/mQ5\nd/7qLOw72xJM8hwAfwbg42Y2N3e4jVvAW5d3kr8K4JSZPYxZVD+nrXlLLAO4GsBtZnY1gNewoUqj\nzfkj+eMAPgFgH8YFwzkkP+Tu0+b8+WTIT2vzSvLTAE6b2Z1b7LZl/uos7L8PzA1wfhnm/zO1Eslt\nGBf0XzWze5LVayQvSrZfDP90s033cwCuI/n/ANwF4D+Q/Cq6kTdg/Nl7zsy+nfz+NYwL/5Mdyd/P\nAPhbM3vBzIYAvg7gZ9Gd/E2EPo8by5tLk3WtQ/LDGFen/kdnde781VnYfwfA20juI7kd48aFe2u8\nfnQkCeArAI6b2ZecTfcCOJgsHwRwz8Zjm87MPmVml5nZ5Rg37P2Vmf06OpA3YNzeAuBZklckq64F\n8D0A96ED+cO4sfk9JM9KPqfXYtzQ3pX8TYQ+j/cCuInkdpKXYzyV0LcWkL5SkiHjPwngejNzB0jP\nnz8zq+0F4FcwfsL2aQCH6rx2Rfn5BYzrsx8B8HDyOgBgF4C/xLj1/BiACxad1pL53A/g3mS5M3kD\n8E4A3wbwKMaR7/kdy99/xfgf2GMYN15ua3P+MP6G+TyA0xi3/31kq/wA+FRS1jwB4JcXnf4C+fso\nxnPCPeOUL7cVzZ8eqhIR6YGm9zsVEZEIVNiLiPSACnsRkR5QYS8i0gMq7EVEekCFvYhID6iwFxHp\nARX2IiK05x+cAAAAB0lEQVQ98G8dYOeiTVEFvgAAAABJRU5ErkJggg==\n", "text": [ "" ] } ], "prompt_number": 16 }, { "cell_type": "markdown", "metadata": {}, "source": [ "The time step limit seems to be to strict, let's try increasing the time step a bit" ] }, { "cell_type": "code", "collapsed": false, "input": [ "U, dx, dt, x, t = diffusion_init(maxx, maxt, D, nx, SC*1.03)\n", "U[0,:] = 0.\n", "U[0, (xmaxx*0.2)] = 1." ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 18 }, { "cell_type": "code", "collapsed": false, "input": [ "diffusion_solve(U, dx, dt, nx, len(t), D)\n", "pcolormesh(U, rasterized=True, vmin=-1, vmax=1)" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 19, "text": [ "" ] }, { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAEACAYAAABS29YJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXu8VdP6/z/jhAjpJFKkHQkRcjnlVjsnSU6RUISK+pFO\noZOKjvZakZMQ6iRO0Q5RyKW+QkVbLkdIrlFy2k5uuR1Crhm/P+YY4/nM1lrtS7u9197reb9eXuvp\nWWvOOeaY+zV9xjOe8QxjrYWiKIpSs/lDVTdAURRF2froy15RFCUH0Je9oihKDqAve0VRlBxAX/aK\noig5gL7sFUVRcoDNvuyNMXcZY9YZY94iX31jzEJjzCpjzAJjTD367kpjzPvGmPeMMZ22ZsMVRVGU\n0lOSsp8OoPMmvpEAFlprWwB42v0bxpiWAHoCaOmOuc0YoyMHRVGULGCzL2Nr7XMA/reJuxuAGc6e\nAeA0Z58K4H5r7a/W2mIAqwH8qeKaqiiKopSX8ijvhtbadc5eB6ChsxsD+Ih+9xGAPbegbYqiKEoF\nsUVhFhvVWthcvQWtxaAoipIFbFOOY9YZY/aw1n5mjGkE4HPn/xhAE/rdXs4Xwxij/wNQFEUpB9Za\nU95jy/OynwugD4Dr3eej5L/PGDMBUfhmPwAvpzvB1fYqAMA1n4wOPjtwewCAmTsr+BLoRXYixU7n\nK429NY9rjwQ6bMH1SvXbYfL/y8SNqc8+0/cV0xftAXSolL6v/Of7EoCpAAbA7t02fG+Olf60naQ/\nzYGR//vDagXfobXfDPZduCDY4zACANAb9wXfInQM9mFYHuwP0DzYjfEJAOBL7Bp8O+P7YP+IHYK9\nHX4J9kZEbaqFjTHfksSzaJdon+L3eH86X0X8dnNt29LrLU48j3aJ9il+ABiDsajuGFPu9zyAEl72\nxpj7AbQH0MAYsxbAaADjADxgjLkQQDGAswDAWrvCGPMAgBUAfgNwiS2hpOZP9baXa+0Y/fReyA2d\nm3Uvg9L8tj2eda/7rfaCoxf4yB+Qyo6bf8GXt22A/5/ZszXsZX8tAGAZ2uIOABfhTphd6AXfjF7w\nO5B/98jfYNuPgu8ZnBDsS3FrsC/BbQCAR0M+A3A8ngv2W2gV7DwUB3sddgcA1MM3wfc9dg52uhc8\nkPnF+Af8jlrYWKEv4op+aZf3ev7elPRs9mVvrT07w1cd0zmttdcBuG5LG6UoiqJULJoHX+HkVXUD\ntip5Vd2ArcwRVd2ArUzT/KZV3YStRk2+t4pgs8p+a9MAXwTbNo6GxAZfBV+C4pTVJ4yzZceV6rdz\nKDq2Y5qYPX2f6FHRMfvKO66y2nY9/g4AOIKSx+x+FLrZkfyNyF/7JwDAh3/YO/hOx8PBvtadFwCm\noT8AoAvmB99SWoayP1YFey3lOeyKLwHEQzc74Mdg/4Ltgl2akEeT/H2wMcNv+ffVMWbv721TvxKh\nyl5RFCUHqFJlv6ZOs2CbOpF6+owmaPfIEuWXdWqW1HrBf5HK3ltvgjbr+qICjhuBNwAAdp8Mar4x\nq3nxf79rpDAPhWTgzECfYI/FVcH2WThlycABJAunLBk4QOWq7myZoFU1v3lU2SuKouQAVarsW5Oy\nsQ1dzL4BxZu/rL7x5vIeV6rffvxTsLGnpK+m+z5B39fIvij3cROD/QYOBQCY+qTmm5Car5uabglI\nymWmdMvBdI05OANA2dItAUm5LEu6JfsrQ3VnS8xeVf7mUWWvKIqSA1Spsl+CdsE2O0fqKZbt8KVk\nLiTQhexE7LOsdjYfl/n7l8QmtV4QdhogYmqejkNbssvetso+bmte4yYMCfah27i/vb0zqPkGqRk4\ngGThlJSBA0gWTlkycABR9GXJwGG/xuwVjyp7RVGUHKBKlf0peDzYdufUmP0cyszpUc0UenmPy/y9\nqPL77BvBhjkUm8LfJ+j7mtMXW37c37Am2L4cQixm35DU/E6pGTiAZOGUlIEDSBZOWTJwAMnCKUsG\nDvs1Zq94VNkriqLkAPqyVxRFyQGqNIwzi0oY+yqDsckwGq4maFKqOoQgyntcxu93ppRUKnVa8CxS\n4dANH/ddaiprteyLch93b7BXgBb0udAhp1XGJmj/mJpuCUjKZUnploCkXJYl3RLQCVqdoK04VNkr\niqLkAFWq7C/AXcG2O6ZO0H61jSiqXX+bHuwE+rnPBPlKb2fzcRm/J1W+1M4J9gTTA5vC3/MooMb0\nRTmPuw3nBrvlzqnFzcyuJadbfvQH2Vb5VDwGoOR0SwB4EcdE18WK4Csp3RIQRV8WNc9+naBVPKrs\nFUVRcoAqVfb/xF+DHRZVUXzUNBGldf0a8Y+oBqq0vMdl/L57+pj90LlI4TtS+7HjHsntmP0lVD7b\n7kF/Zw3d31799OmWP9UT/0F4O9g+5bKkdEtAUi7Lkm4JSMplWdIt2a8xe8Wjyl5RFCUHqFJlfwXG\nBzssquK46S6cmfNZsBPYw30myFd6O5uPi/uuFZtU+ef2pmBPM3/DpvD38Zg9nc/FmatPX5T3OImb\nryL1bPZKzbaJLaqqSxk4238U7GfRPth/xT8BlJyBA0gWTlkycACN2WvMvuJQZa8oipID6MteURQl\nB6jSMM4YFAS79o7fAtgkdEMTtMteF/8RiNIwfQpmZCdKbVef0IWk9CVuTT9B2/9+pJCk0E7suEtz\nb4J2OlVLbdEoQ2qln6ClidhaO8sk6X//IBtZ+3RLAPgHRgIAbsOg4OuKecH26ZaApFyWJd0SkPBN\nNk6C6gRt9UKVvaIoSg5Qpcp+nFNGALB++10AAOaPpL52IvVF/sT/In82q9LyHhf3SWVGVuW17WXB\nnmluwabw9/EJWjqfKxdQffqifMf1o5RFVu5mj9TJ2Fi65c6bT7cEJOXyfNwTfE/hpGAfiVeDvRL7\nAyhbuiUgKZfZPAmazW1TBFX2iqIoOUCVKvshlLK2d50PAQB2B1Jfe5L6Yv//1gIAEhT/zDZVWt7j\n4j4p1pV4KX3MvvftSCFJaj92XNtcidnL7lxrKGXR5NHfEyl3s7uL2VO6ZUP39wjE95j16ZaApFyW\nlG4JSMplWdItAY3Za8y+4lBlryiKkgNUqbK/AxcF+w0cAmCThS11SH3tL/7Fn0T+DpgVfAkql5wN\nqrS8x8V8TdKr8uZWyiE8bKToWbrvYzF7Pt/a7J/3KPs1bgMA3Eu7ejVrnmEOiEeN20f+WAYONp+B\nA0gWTkkZOIBk4ZQlAweoHnHxbG6bIqiyVxRFyQH0Za8oipIDVGkYpy+kRv1xeAGADKmBTYba26Wm\nYfoUTCD7QhDlPS4BqmuzVu5v0saPgj3Y7BXsAvl5IEmhHT4uUYv76yb3SQuwsq4vynqNSwAA5yLD\nxH7TDCFCt4cCp1seAtm0ndMtx+DqYPuUy5LSLQFJuSxLuiWQ3ZOg2dw2JRVV9oqiKDlAlSr7+9A7\n2E+iMwDA/PHX4LOkRM0+osp+eSnyb0f1yROkmLJBlZb3uJjS7kcTqtQX7ezRwX7C/Bubwt+DRgGx\n802viRO00aKxz5B+Yj82OmycqvJLk255OW4O9mz0BFByuiUgKZdlSbcEqsckaDa3TRFU2SuKouQA\nVarsz8BDwT4TDwIA7LbbBh8XQrNUCM0cHvkXLhbfiVgY7AROJDuRYme3mqW6/dPl/p6yzwQbRlRn\nwRikkCS1z8fFSyd85j73KHXb4u0s/W/Le1zpfitbdc1xi9D2aE1/NyZDzJ72N65dPyrC9yH2Dr6/\n4PFgp0u3BIDT8AgA4Dm0C76S9pgti5oHsjsuns1tU1JRZa8oipIDlFvZG2OuBHAugN8BvAWgH4Ad\nAcwG0BRAMYCzrLXfZDrHPHQN9l24AABQb9dPg89y7HU/UmXLXeGqBhSD/jI10ySyU7NNslvZk9K+\nM32JhDPsvsFeZD7ApvD3PAqIne/C6hyzv57sbsHu4fcx/o3+bg7IELOnQmi/OpV/UBkycADJwikp\nAweQLJyyZOAA1SMuns1tU4RyKXtjTB6AAQAOt9a2AlALQC8AIwEstNa2APC0+7eiKIpSxZQ3jLMe\nwK8A6hhjtgFQB8AnALoBmOF+MwPAaVvcQkVRFGWLKVcYx1r7tTHmJgD/BfAjgKestQuNMQ2ttevc\nz9YBaLi583TCU8EehMkAgP/WkpoksU2h36Yh+CGRfwNN2tahYXCCJruqKgRR3uMSzVNDLQCwyt4Z\n7IfMhcEuGIEUkhTa4eNiE7TuOonV2bswLfP3I8jeEOyvfnThvcPo72Yj/d3sTf4/UIXLulHK5bM0\n0ToQUk40XbolAOSjCADwOg4LvpI2FNcJWp2grSrK9bI3xuwL4DIAeQC+BfCgMeZc/o211hpjbJrD\n8WxiCQDgTXyNA/J3x4H5u6f7maIoSs5SVFSEoqKiCjtfuV72AI4E8KK19isAMMY8DOBoAJ8ZY/aw\n1n5mjGkE4PN0B7dPRAqqGHn4BFH85wZcAQA4qNY74Xf2d1JlXL3wXafgjhbfwkWchikTbQkcSnYi\n9llWe+setyz6ZKX9TvoJ2kF2p2AXGZnwS/f9ZBoFxM53kJ+gXUZtOKKC72nLjsv8vTzf+agT7F2P\nTjNBy4uqaMK/9u7fBtunXHK6pf97BICJGBJsn24JSMplK7wZfMW0B0FDrAv2N6gHoGxqnv3ZPAma\nzW2rzuTn5yM/Pz/8O5lMbtH5yhuzfw9AW2PMDsYYA6AjgBUA5gEhjaEPgEe3qHWKoihKhVAuZW+t\nfcMYczeAVxGlXr4G4F8AdgbwgDHmQrjUy82dx8c8AaAA0f+1XqJa5LUbi/qy75FaO9ApOI7jNyPV\nuobj0HPI7uE+E+Qrvb11j4tUdaJ/qvoGgF+spP9NNtcEu+CvSCFJap+Pi8Xs3XUS06pLzJ6fo4zW\nuvDCu+/ciO9g8v1KfyN7paZbApJyeS+V7/B/j0DJe8yWlG4JSMplWdIt2Z+NcfFsbpuSSnnDOLDW\njgcwfhP314hUvqIoipJFlPtlXxHwrj4jMA4A0BlPBN/PG3cJNhdCs+87BXeU+DbcRZk5NC+csNWk\nrO/OqUr7DvufYF9k9gn2aKn/hhdEPKb9fgyNAvh8XuX76wJA4rvSq/zKUfa8OI5335I2f/sFKfe/\nuBHfDxkW49EcUOO60hfPoAOAeAbOMNwY7PtwTrB5NOqzcDgD5xM0Dna6Xak0Zq8x+6pCyyUoiqLk\nAFWq7HmZ+SSX8fAgzgy+hrtJyVn7H1JrBzkF9y75OtEetY+IvwMk7p/ALu4zQb7S2xV/HGUNOVWd\nmJE+A6eh7RvsMaYw2AX9kUKS1D4fx6MDf51EH1bzm89iYrtylD2PwOQ5LqTR2i4nknL/xo342pPv\nJxrxNf462P+lomed8SSAsmXgAJKFU1IGDiCKvixqnv3ZGBfP5rYpqaiyVxRFyQH0Za8oipIDVGkY\nZzktM++PaQCAPrg7+N7fuF+wTZOfgh0maNvScJ1qv5sjxX/9q+If4WqNJ9Ay+Co/jHMb2RQq6Zoa\nVimyMlmdb04OdsH7wcRS6aK03ycp5MPnCxO0XSlsNI9DOtzOS0q4p9LbZfuthJWuh0zWn8j16mmf\nXuN24rLfk++o34P9+2+ibQ7Ca8GegfMBxNMt+6Ew2PPRJdjpKlw2wdrg4xIJ6SpcliXdkv3ZPAma\nzW1TBFX2iqIoOUCVKvtD8Faw73bq6p+0E9DhtUR92Z+3D7Zp6RTcY6TgaKJuw0SalGtEyvVTXyKg\n9LtasV0xapZU8s6pqjqxMv0E7ZH2z8FOmqeDXXAeUkiS2ufjQKMDf53E/qTmy5CGuXWV/UL3KSOf\nEfQcN6yi5z6QVP46N+LrRr4NlG65j6RbLqGJ1gGYCqBs6ZYAsC+ignPpdqQC0hc90wlanaCtKlTZ\nK4qi5ABVquxXUOy8J2YDAIbjhuB7BrLLUoOmHwXbunIIpg0puFmk9s4S//MUyz/OLcaJL7RKLacQ\n2YkUu2LUrKx4iqnnm1KV9lr7z2DDSF2EgpfE/ZpUl0j7PY8C+HwhZn8Tqfm/cb9QO2nZ/9ZT9rPI\nPtG1kZ7jpzRaO4+e+8f03Hu5Ed+34tu51RfBXkMpkqekKXrG6ZbdKd1ysVt0BaQvelZSuiUgij4b\nVbDG7HMDVfaKoig5QJUq+xZYGexH3aZWSRQE36l4LNiffb9XsM0BUWaOXUCqLp/U3hTynyz+iU9E\n/iEgNQtWs/eSfS7ZidhnWe14GWFSyd1TVfU0K6k0/Y0E37vYg4OdNG8Hu+B0pJAktc/HzafRgb9O\nrDgat+cR7pfUMsgVo+y5v3uRHbVjIo3AjqPnGFtMN4D8rnSCOfFX8f2wbbBb1pE5oJlU9GwUrgWQ\nOQOnDZYGm0ejPgsn3SYlQDxm77Nwsl0Fa8y+5qLKXlEUJQfQl72iKEoOUKVhnA/QPNhdMB8AMI72\nF72XQikH1JXQhf0kSsM0h9MQfg4N7Skc8ct48W/nFlslXi3NRGRqSKfsoQufQkg7QDXPECp51rWN\nwiqw0hfzzfXBLpgvP3lNog1pv+eQD5/Ph4j8dQEg0T51j1oAm+xT6++p9CmrbGcOlfFzcOE2Wjz1\ny4v0fIfTc/+Q/D3cBO168TXZW8JinG55Ae4Ktk+55HTLP0MmtpfR8/PploCkXJaUbgnoBK1O0FY9\nquwVRVFygCpV9lwH3O8ANASTgu9iTAk272DVoEWUhmlp39lYpUOqCW8uEP+ymyP/ERl3tfqM7D3I\nnuM+OTVzKtkDyE5dsJXYIYNKHpeqqjOVSOhnZRIwaWRr34I0yj5JPj5uOo0O/HViE7TcnpHk5/b/\nWNLCtEz9kq4Pub9pVOWez0vLaVR2OT1f3qGsb+qiqnp//jT41kAqffrqlkD6CpdlSbcEJOWypHRL\noPqoYJ2grbmoslcURckBqlTZFyMv2H5J+jT0D74bMDzYJ5KS/OKbKA2z1qFSaMr+H6m9U0nt3Zha\nKGsqLbQakKacApBe5SfwPPmOIzs1NTGyXRz+Rzpvv/Tq+cYfosU/w8xuwdfESgx5urkv2AX3BBNv\npimXwN/zKIDP50cN/roAkNgxQzunc7+4e4qlrGa6/9T+yjh6oucw1Y222vLiqefoOQ4jP43Map0W\n/T1s/Han4Gu128vBnklzBCPxj2D7Inzz0C34Skq3BCTlsqR0SyC7VbDG7HMDVfaKoig5QJUqe1ZJ\nfj/aczAz+LjkLO9gtV/9KMNk4xui4MzxpPZmpKp5ALCj3aKbTIuH/kj+/7Fy/dZ97kK+tWQ3ITt1\nwVaiYwaVzJkwTlVnKpEw8gdxJ3cUu4DqnIXvSe3zceNodOCvwyOJjJk53P5FPmafaWFapn5J04d/\nTD+qGuCej51Pz3EMPcc3yd+b/J9F/rxD3w2+xVRyYwDNJ3DRs3sQdVhZMnAAycIpKQMHqD4qWGP2\nNRdV9oqiKDmAvuwVRVFygCoN43yCxsH2E2JzcEbwjcS4YJ9PO1i9gqMAAA0PpQ3Jn6Ch/Sk0tP8n\n+YdG/jdGie/QNCEKIH1IJ2Pogqozxipq+kVcfN470odKHrX/BgCcZo4OvjY2P9jjTFGwCyjS86ZE\netJ+zyEfPp8PEfnrRveRoZ0XpbmnV0tx/+lCWplCZfQc3nChNTOWnuPz9Bz/mn5RVYOTopTc/2yk\ndMta6dMtb8blwfYpl2VJtwQk5bKk0A2Q3SEPnaDNDVTZK4qi5ABVquwb45Ngv+omxLpiXvBxffHJ\ntIPVCXgGALDu06bBV+fYr4Nt7ycVeC6pwGucYrxSfPf+Q357bgdSnYtJdTZxapb2O000oN9+map8\nAVG/iavTq+Trv5U2jzD1XSOlpMFSLpEgmywhKcIVBccihSSp/dhxNDrw1+GRBLcnsQvdE7ff9WG6\n+wQ20y/p+pD6+14a/Rzqno+lPXENtcG+RbXtz5Y2b1gX9WGrRunTLf/uqlsCQF9MD7ZPuSxLuiUg\nKZclpVsC1UcF6wRtzUWVvaIoSg5QpcqeVVJrvA5AyiYAstgFiMfvH8OpAIDmjaTI14Zl9YNtTqZ6\n5reROrzMKcZB5KNl+ImbS0g3bEU+UpeJo8n/b/IPdmr2GvK9RL8l9fyWjQqEtTKiRLvbFsFOmlXB\nLrgumHjzKqTA3/MogM/3iBs1+OsCNLrYtJ1t09zTpFLcf5r+StevAHAul0P4lxuBTSbfEnpml9Dz\nfV/q1Tc/Ivp74B3OuODZCPobKkS/YPuUy7KkWwKi6DVmn71tUwRV9oqiKDlAlSr7BlTW9i20AgDk\nY3Hw3Q1ZHTQaY4Ld25XJfePnQ4Kv8RESnLbzRO2Zc0gdjnOKcbz4VlwoirHlYFKdrFy7OjVLMeQE\n7ZyUoEygxMDUcyQeS6+SH7ISWz7D/AlAvKTBI1wi4dVgInmk2AVkh+9J7ceOo9FBuA6NJLg9scwc\nbv+p7p74PqeUol/S9SH19woaVZk73QjsMfKNpOf4hvgbd5Xn/v7PUcnsE2svCr4JGBrsGzEs2Gfg\noWAvQkcAQGssD77VVH67pD1mS1Lz7M92Fawx+5qLKntFUZQcoEqV/ZdU1tbvR/scbTDRE7ODfT1t\nanI7BgIA2tdeEnyfrJXgdL2TpMStnZqan23/Rj6KC8+hWH6PNLH8BBXmStyTXqHGRgQznJo9VXw3\n/yxtu9w0CvY31hXmMlcG30VWSgskzbfBLiDl/ibF59N9z6MAPt8dbtQQrgsZXWzazkTtNPfUpxT3\nn6a/EtSvc0jNt+T4vHtm5qb0hdDqDZC2fbNW+rB1k2jNAG96c2WagmcA8Ai6B/sYvAhARpdAyRk4\ngMbsq0PbFEGVvaIoSg6gL3tFUZQcIGsmaP1+tLyw5VEaag/GxGAPxQQAkoIJAC2bvBbsb5bK0L72\n2RL+sOPTLMM/m0I6t4r/+kvFP8It8knQAqwELcyK+bnMgAt1JD74XXy15f+vy+wDwT7CnAUgXtLg\nDi6R8FQwMVayUzHqMKQwlkI7fByHgvx1llLYiNvDIaZY+/d195SpnEKmfnH+68nXg/rbJuk53O/C\nbXPFV3u4tP3n5RKOatlGnvtCt2NWfH/ZG4LN6ZadIB3zMtoAkFAiEE+3zFSvXidos7dtSirlVvbG\nmHrGmIeMMe8aY1YYY9oYY+obYxYaY1YZYxYYY+qVfCZFURRla7Mlyv5WAPOttWcYY7YBsCOAUQAW\nWmvHG2NGABjp/ksLp7H5/WhfpYUtXTA/2LyDVRIFAIBeuD/4lm84XM7VRuqZ//yIqEC/GMcOIxU5\nXdTlqlPE34ImDBNu4jbjHq2ZVO4yNyLYN7XgGRAvVbCzvQQAsNTcFnyjRMwiKbeBAiqH8CYVPQvH\n0fdJHgXQ+ca6UYO/LiCji03bmTCiCcI9HVGK+0/TXyOoX1fx4rbHSeXf7kZgN9LiqWclnTavuzzf\nlRsODHbHOgsAxNMtr6c/P063XIz8YB+CtwCULd0SEEVf3VVwTZmgVTZPuZS9MWYXAMdba+8CAGvt\nb9babwF0AzDD/WwGgNMqpJWKoijKFlFeZd8MwBfGmOkADgWwDMBlABpaa70kWgeg4eZOshO+C7aP\nkfqyCYAsdgGA3pAFRj4NcyoGBF+7Os8Gu3iNqL2G3akM8sRIHZprSUX2Sa8u55PK75Iu3XBGamwe\nABIvkN+p32n2/eDrb/YL9ld2fLBhov12u1lp+1gjCrZABjGYcLbYQ+Xn8j2pfT6ORwf+OnNpJMHt\n4VEHt98vtord57Gl6Bfnn0++FqzmqaSEcb+1j4qv4RB5juvWSAG8Ns3kuc/A+QCA4RSnvwh3BPsh\nKp/t0y0BSbksS7oloDH7bGibUnrKG7PfBsDhAG6z1h4O4AdsEq6x1lqACporiqIoVUZ5lf1HAD6y\n1r7i/v0QgCsBfGaM2cNa+5kxphGAz9Md/GwiWgy1EbWwd34e9s5vFsod88KW4yGLpliVDXGZOZdS\nhs5jrkwtALRqJsv+170gKrDexdFiHHtV6tJ8ALAdyP+U+CeeFPmHpCkbACC+dyup3Bt/+AJAfJ/X\n161I7cOMSPT9bRTxmmseDb4CEdRIyoAABX3FfrMQKQyl75M0Coidz40a/HUBYKUbXWzaTh6N+Hvy\ne+YCyLx3LfXXRNdfXahf7XnU34tTC6HVu44WT70g2UGtjpXn6zNwAKCP2+CGi+ZNjxU8kzIKPgMH\nkCwczsDhTDGN2WdX23KFoqIiFBUVVdj5yvWydy/ztcaYFtbaVQA6AnjH/dcHwPXu89F0x7dPRKtk\nefirKIqiCPn5+cjPzw//TiaTW3S+LcnGGQxgpjFmOwAfAOgHoBaAB4wxFwIoBnBW5sMVRVGUyqLc\nL3tr7RuA2ww2Tsc0vrTwBK2fEGuJFcH3Io4JNu9gdTsuAgBci1HBdzZmBXv5eqnxsv+xbwT7m1lR\nKKDOaNrVqn/6UMKag8TfzIUpYiGKDPXer/tK8huv2jGaEX3KPhN8JxmptV7bXhbsleYWAEB/K/WC\nkkZCCbyv7CRKrRxM9erD94Vix/ajpVCQv840ChtxezjExO33ISm+z8SuJffLENeHa6gPzTsUuknQ\n7lPTouez4WGpr79/L3mO76w/NNgn1n0i2D7lchzVUToDc4K9CH8ONicCrMT+AOI7p31FdZtKqn2T\nzSGPmjRBq2wZWi5BURQlB6jScgmcxtbQzeV6lQUAR2JZsHkHq/NxDwDgOlL2sTTMuqJEV74nKjCv\nVzQpuWGCKMbat1M5hZNIda4U1blw/8h/4hukWg+V34757Ptgj951p2D7+vBcTZLTG3elCdEONvrN\nNCOTjwUvBRPJtmIXnCn2aw8ihcH0fWw/Wj6fGzX46wLAYje62LSdPBrx93SVkTxOvv/EHqTyqb8W\nuv5qRv1qe1M5hKeoHMK06Nx5QyX1tPg9yTE95gB5vne7dEsAuMKlXHJ1S57YPx7PBft1SJ0Jv6Dv\nE8gkcEnploBO0OoEbfVClb2iKEoOUKXKnlWSj9nvi9XBx+qrA+1gNRs9AQCX4+bgG4Jbg/04Tgn2\nkQc8H+ziZyJ12GSo5CD+fJkoVFaX9kBS+R9Gy/anNhXfRazm9xA1f4+VfXHPMwcDAFbYwuBrafoG\nO8/K/PUMZOc7AAAgAElEQVRiExUhG/lDcCG5o9i8r+wUqlc/sDFSmEJqP7YfLY0O/HXG0UiC2wMa\ndXD7/SiF73O0kftnlT+VVP6Jrg/tUaTm3yU1P0aeQ5Nbouez9hmZZDjyBHmOnG7pdy0DgJGI9tXl\n0honUcEzngNKV/SsLOmWQPVQwRqzVzyq7BVFUXKAKlX2rJK8qipGs+Br5QpUAcASHB/s0/AIAGAy\npIjXOEip3l6UmfPqF8cFu/UJUXGvtfeKYmx8i+xh+nMPUZc7r/ki2Gt2jMostPxWsnhG7yJqNlM5\nhBftYwCAY4yUYuaMl2KKkfe00eKvcUbKAhTMDSaSsmYMBV3Efk1qxQUG0vex/Wj5fDvGrwsAs42U\nOOZ28mjE35MftQDx+2eVP4T6a80uUR/u/IP068/nSX83niPP4ZN7oxSj1udKMbblX0j5ho67/V+w\nuejZWEQ3yzuc8VxPugwcQLJwOAOHM8V+RJ1gVxcVXFNi9krFocpeURQlB9CXvaIoSg5QpWGc7fBL\nsP2EGFceXIGWwW4DmUic7yZg+6Ew+MZgdLA5DbPDbrLoZvnbUSig1blyrk/+IauSOJTw3VFSz6bB\nxo+i69aSlM3zMoRuFtkFwe5oOkX3Rpt616OdoXhXqtmuvnyBNAFJWjBVIOuEMO16sftLcUb5nkI7\nfFwsFOSuk6SwEbdnKYWYuP0+JMX3yffPIZ35RvrrSNeH37WRfm38CoVu6Dm0ujJ6Pm+9LWmhHQ6W\n5zgDfYLtdy0DJOUyU3XLdOmWgKRcliXdEsi+kEdNmaBVtg6q7BVFUXKArJmg9RNivLCF0zCX0Q5W\nHV31wvtwTvBxGuZfITUCHtsok6PHHxyp0bcWiGJseaXsYfrJQFGXea/Igp4vG+4V+az4xpGanWil\npnoX0z7YxXZKdJwZGHxNrLR5qZEa/YNsNLGZNJK6WHBPMJE8T+wCuQSWyqUD/en7JI0CYufbJ35d\nAJhMe95yO0GjEX9PftQCAPPp/rlf+to0fbhOfNzfLafIc1ixINp17PhOMnp4fONfgn1uLUm3HAG5\nwTtcGQ3e4ew5mtjnUhxrkBdsnxxQlnRLIDtUcDZcT+vSVw9U2SuKouQAWROz9zHShlQCP1Mapt8/\nlPcUvY3SMP/BaZi1pC77c2sjNXpMJ1luv+JO2bu21RSJ5RefIsvz918XFeEqMuI7w8oin7OMpHe+\nbEV1tjLnAoinMa6lWDjvSjXZ1ZfPWCKByh5Mo+JmfWj3qfA9qX0+LjY6cNfhkQS3Zy6NOrj9fpTy\nFt1nDxrN3Ez9wv21v436kPu11ePS3/wcjrkwej4vrpUyDV2aSEGzGyALvq6D5Jb6lMvHIbmnXHKD\n54DSFT0rS7ol+2uKQi/v9bTsQfVAlb2iKEoOkHUxe977M1Nmjs+wmEu7U3HxqzG4OthT8f+CfVKT\naEHQi8tEMXoVCQBvFUgsv/XjsqBnZeuomNqRpFpvJTU/1orqbGN6BPtHm4gMk5DzWol1zzUSkx7t\nVuonZV0PCqQCBJKXkp8U/1IaCXj60/dJ3o+Wz9c2fl0AGEN73nI7l9NoxN+TH7UAQBHdP/cLj358\nH7ZeLv3K/X1MUp6Dfz4nHfFY8N2FC4I9FDcF+2LaY3YmegOI73C2DEcGO10GDiBZOGXJwAE0Zq9l\nD6oXquwVRVFygCpV9vx/e6+quBgVq6/9qXCVV2tc5IpL3Q7DjcG+lAqkzd0QZeZ0PEKW2784j1Q+\nqcvlA2R5/jHLI/+rDUS1nmRFdV5A5RCesaKCG5so1s2ZLcspFt7PyihmjInmKgooRz5JZQ8KKN5e\nSFk1PalYWvie1D4fFxsduOvwSILbM51GHdz+tW6U8gndZx8azYyhfuH+OubLqA9j/TqV1Dw9h45d\no+fz+AbJwOldhzNwZI/Z2zAo2D4Lp8jN6QAlZ+AAkoVTFjXP/pqi0Mt7PVXz1QNV9oqiKDmAvuwV\nRVFygCoN4/CQzw+heWELp2GuRvNg+zRMHq5zGuYkDAn2OIwMdq86URrmotUSHvAhAwB4cbKEEjpM\nleX5L3Z2E4ZfSojiDgrdjLb/CvZRRiaEf7KDAQBrzaTgO8nKkv3pRqowFrgqA7xPbMEYsZOjyd9K\n7BckIzXQl75PUsgndr4u8esCQNJIf3M7n6LQUx13TwcbScd8nO6f+4VDXb4POzxJ/Ur93XGQPAf/\nfLo3l7TZG3BFsLk0Ble4nO9SLg+j6pYlpVsCkhxQlnRL9teUcEx5r6ehm+qBKntFUZQcIGsmaP2E\nGC9s4TTMPKwJtldrXORqLroGm9MwC5AMti+Qxopx0RJR+ScNEiW6eOTJwe7yZJRa+FRLUa297Z3B\nvsRcGOz5VtI+tzPXAIinMT5FE5+X/yaTgMltogVmBdI0JM8Wu+B0sQsfFrtbmidYSGqfj4uNDtx1\neCTB7bmZRh3xNMxolPIN3edFNJoZRf3C/dVlRdSH3K8njZP+5ufQvV3UOC5oxxPtF8XSLWXyuJ3b\nY/bVWLql/N34HamA+ATtdy45gBf56QStlj2oaaiyVxRFyQGyJmbvVRUvbGH1xarMp2EuRZvgK0sa\n5sNfi2T2KhIAnpolSrTbONm1af75UWphzxWFwXdvvb7BHk4lgI+nomHfuJRFTrfknaFu5l2pfPkC\nLpHwN7HHyjoijCI1/gLF3D196fuxNArg8/lRQ6w8wzaibOM7WKWmYdZzoxYAeIjuv7CejHJ6flMY\nbN+H3e6WfuX+7t5LnoN/PufWl1ECP8fJsXTLx4NdhA4A4umWq7FvsHkOiGP2fr6oLGqe/TVFoW/J\n9ZTsR5W9oihKDpDVMXvOzGmMT4Pt9w/lrIun0THYnJkzkTJzfIG0XvULg++Rt0Xln9FLFu7MHX9W\nsHveHf1+dn7f4Bv0jUjtK4xI5oetFGTbwdwGAOhgpSzAbCPFvwpEaCLppicKqLxBkoqYFchWqiiU\nQQxORCqFpPZH0XFJGh3468RGEtweGnVw+99xo5Rius9LaTRzmZWLcH/1LCoEEO/XM4ZLf/Nz6Hlw\n9FsueMZzL5yBM4/mao5wRc/egqQjccmNz2kOKF3RM43Z60Kpmowqe0VRlBxAX/aKoig5QNZN0PLC\nFt4T9PM0aZg8XD/epd0B8aE9p+n5UMBdP0sFxd4HyyTgQwukkmPv4eKffVlfAED/IomxTGshMZak\nlY1eOxrZOWmdjdIMXzKykOgiKwXok+ZbaZuvVcP1cC4We/ztYg/fW+yi/yKFvvT9eAr58Pl8iChW\ni4f2s+V23k2hp+bunhq6EBUAzKD7n9xCQlr9V0l/+T7sfUuG/u4kfv98/lp7srSHnuM9kII//Nz9\nhD2nW36KxsHmv6fvKBHA/+3pBK0ulKrJqLJXFEXJAbJugpYnyXiCdtc0aZhcCfNFHBNsTsOcjr7B\n9ul7rBhnrpVUwfM7ybL/mXeKv/8tkUKd1kPU/MhViWAX7CT2/VQhsp6bzOxsWwTfHWaVHCfl45F0\nGzgVXEc+2YQJBceKXfiC2PlIpZDU/nA6LkmjA3+d2EiC20Ojju7U/tVulLKS7nP4TjKaGfl9Itjc\nX/3nRH3I/Xr+hdTf9Bz6NYn2ueV0y4kYHOxTaI/ZRTQx71Muy5JuCcjfXjbWfs+GEYFSM1BlryiK\nkgNkTczeK4l0u1cBcVXm0zBZwR2BV4PNaq8Xpel5dfgPKo7mVSQA3L1Ulv33v5Di824j18vnyOKh\nm4+SdMMJ30saYleKZRfbqJ77y0Z2bBr5QzCRpFr0oXwBl0ggewKVURhKNegX0U5Tnr70/QQaBfD5\n/KghVp5BtoeNtfMeGo00dffUihaKTaE0zHFHJYJ9+SvSX74P+9PWWdzf/drIc/DPh9Mt+Tk+gu7B\n5ufu53A4TbekdEtAY/Yan88NVNkriqLkAFkTs/fqgmP2JWXmNKM9RV9H62BzhgarQF8gLVYcbf3A\nYF/URnZfmvaIxJuHJMcDACYOkEU+174iWSd/31MWEj1oJRNoNzMPANCJdoCaTGWEY6UK3OKmAorT\nJyl+XyAVh1Eoa8moYAR9T2p/KB2XJBXvrxMbSVB7JtCo4wxq/3/cKOV1us+he8po5tqPpV9uHiCj\nnyFToz7kfr2ou/Q3P4dBdSM/F7S7h0pg8PPluZrmWA0gXlqjpAwcQGP2quZzA1X2iqIoOcAWveyN\nMbWMMcuNiSSsMaa+MWahMWaVMWaBMaZeSedQFEVRtj5bGsa5FMAKIIyPRwJYaK0db4wZ4f49MtPB\nJU3Q7oANwU6XhskbSHOlw0xpmIXo5xopqYKX1pUQzB3vye5LQ7qPD/bEiVH4ZvRUibGM6Swxlqkf\nyyKfHtvK1lAf2IMBAK+at4NvqGyWhPGy3gcFLi0ySQufCrqJPWGu2AMpxFJEE6ke3oR8AoV8+Hw+\nRFRA6ZjjqU7OcGrnwxR6auLu6ZBt5wXf1F/l/v/eWfpz9JPSX74PhwyRfuX+HnSAHOefz40YFnxd\nSki3BICViFJEM6VbZpqgzbbNubXGjbI1KLeyN8bsBaALgGkAjHN3AzDD2TMAnLZFrVMURVEqhC1R\n9jcDuAJAXfI1tNauc/Y6AA03d4KyTNCmS8PkioZe1QHAkRnSMM/Bfa7hlwffPzbSJOIBpOafkcnY\nkUMSAIBxwxLBN+FJSTcceohMUD75a36wG5siAEAHu0PwFRpZzDOcShWM9XvCUqXLsVQBk2vYc1XL\nw5HKbFL7Q7m2PY0O/HXG0khiFLWnkEYdp1L717pRyss2P/gGHCKjmQlvSr+MGSajn5E3JgDE+3XI\nCdLf/Bz+XutaAPK8AOAhnBFsfr7LITPQ/u+BSySUlG4J5OYErZJ7lEvZG2P+AuBza+1yiKqPYa21\nAOwWtE1RFEWpIMqr7I8B0M0Y0wXA9gDqGmPuAbDOGLOHtfYzY0wjgIKnxLOJJQCA3/EHNM1viqb5\neWlj9iWlYbKCa44Pgs1q73gsCbZXh1xUy6tIAJj4qajOYSfITkw3Fl4NABh3o8SYR54laYOz3pRo\nVa8Gjwb7Pbfb0zKqDd+XShJMokVMo9w0wtgR5Psz/fZpsXvTU3vhN6TA+9JO4tr2dD4/ahh1Pf2W\nSicMpnYuotHIXu6eOjcoCr5ZX8r9Dz1LRjnjHpD+GleYAAAM6yv9yv09pJGofP98pqF/8PFzjKdb\nynMvdnM4ZUm3BHInZq9UL4qKilBUVFRh5yvXy95aexWAqwDAGNMewDBr7XnGmPEA+gC43n0+mu74\n9ol2AOJ/qIqiKIqQn5+P/Pz88O9kMpn5x6WgohZV+XDNOAAPGGMuBFAM4KyMR6BkJZIpM8ertd1p\n4FAcy8wRWcoq0O9XyopxBGXmDGtEan7p1cFO9o3kdsE18tvbH+gb7IuPLgz2ki9lV6d9XWngfKkp\nhtlSORiDHxR7wpnR5yg5LcbLaWNljbnQGQ0OAnNJ7Q/mcsc0OvDXmUAjiaHUntl04h7U/nVulLKE\ndq/qdbT8P/32f/cN9shrZPSTvDq6EPfrsDbS3/wcfBYO7y/7FGTLLX6+q2iuxv89xAueyd9Nugwc\noHJVt5Y9UKqKLX7ZW2ufBfCss78GaEZUURRFyQqyrhBaug1NgLgq29llWHCRq0aQxPBVscycZcH2\n6vAczAy+m0Fx5fWiOke1EXvsI5ECnXA1ZeAMkNj03H93Cna3pguC/V9XZmAZ5an3pJIE0yivfejo\n6HP8GPENp+8n0XGni4nXkArvSzuJyx3T+fyoYfhoas+ZYven6y2l0Ugjd0/tmsqGJnM/lPu/eEBh\nsCdMlf4qeCRS7qO6S79yf4+sKyrfP5/Z6Bl8/BxfpzkZfu5+Dkdj9oqSipZLUBRFyQH0Za8oipID\nZF3Vy0wTtDwE90Nz3r2K0zD3pXQ8HvJ3QBGAeHjgYt6jti6Fbt6TsMK47lGoZ+REmXCcMVXmnvt0\nfiDYr3x4cLD33yFagLQ3lR54ghYr9Zeijpjk5oyHU4xmwsNicw37aVTVMt0E7RNkD+ba9hSa8deZ\nRGGjwdSeJyjkczK1f/0OUUjqlR/lPrt1ltDVjCelX4ZOpDTMIVEfcr+OOkD6m5/D7bgIgDwvAHgO\nxwebn+8HtKeB/3vg0holhW6AmjNBqyibQ5W9oihKDpDVE7QlpWGyguM0TFZ7B1KhLK8Ou0GKeHEa\n5lUbZXl/8gDJSSxYEE0uTh4i+6QOuuzOYD/5ZH6wO7csCvaXP+4EAHjXfB98J0tdNsyUbEIMvjT6\nnHCr+Ia2EnvSW2J3FhPkDrQnexLXtqfz+VHD0EupPdIV6E3tfJdGI7vZ6J6OainF3Z5ckR/sPpfJ\nKGfyLdJfIxdEo6JkJ+nXERsl3fK6WlI0zT+f+ZBVXvwc30XLYO+epuhZWdIt2V8dJ2gVpbSoslcU\nRckBsjpmX1Ia5s4Z9qjlAmmsAv+EpQDiivE8SBGvG2tJSd0xa0Xl39Ap2rv2ijsnBd+sW6hEQg9Z\nVPTOChlVtKgXxZYbUMmCIipM1ptU/DSnsIeS2p9A6prLGs+kQmd5SOVZsvtzuWMaBvjrTKM29Ce7\niNqRT+3/tV40SnnnG7nPzj2Kgj1rjvTLoDtl9HPDhVEfFqwVNT+6iah5fg73ICqZ7J8XALxMe3Lx\n8/2E5mr830NZ0i2B6hezV5TyoMpeURQlB6iWMXuv1ljB8UIaVnstsCrYXh3+GYuC7z6cE+zBmBjs\na5vQHrNLo401pl4om3QMKBAlumSOlA5od5QsNvphXfT/0o+3/z348inD5mHKvOnv4uWTKCNm6D5i\nT/qP2ByTp3plgT+RPY3LHfP53KhhMMXpH6b4/enUzmIajez2U3RPBx0lGTFLXqHSCQUyypmalP66\nYmk0Krq2jfTrJZgs7cGQYPvn8zQtxubnyHMy/Nz9HE5Z1Dz7szlmryhbiip7RVGUHEBf9oqiKDlA\ntZmgTRfSybRH7e5YF2we8h+GaEPWInQIvtPwSLDvgGzbNPprmUi8tc3/AwBcOutfwfdw8uRgn36+\nLGNa80qjYDfb81MAwI60D+xSWeOF02UOGIVurnLwseKb9ILYXD6UIixpJ2jpsNhxHAry1ymksFFf\nas9SCjG1ofZjnygkteZjuc9250vo6uG7pV8GzJJQ1629oj4c9bXsNTumvqRh8nN4FN0ByPMCgLcg\neaP8fD+nzdD830NZ0i3Zn40TtIpSUaiyVxRFyQGqzQRtujTMdJUwgbja4zQ9rw6PwYvBNw/dgt0P\n04N9Y/3Bwb7i7WhycUYvKpEwXhYPvXS3yPW2+aJGrUt1/JxKFrQRsYt5Mn+JvmdHn1PuF99gKeqJ\nSbTnF0/ArkYqvC/tA2Tz+aY4+T/wbGqPZEKiK7VzHY1GdneLtJrlfxp8LxXJD04fL6OcGcOlvy59\nOxoV3XCw9Ov5uDvY09Ev2P758F4E/BzXokmw+bn7CfvqPkGrKFsDVfaKoig5QI2J2XMaJhdIYxXY\nAisBAEtpgU5HSsOcSWmYl/98c7AnHxwt+x80TxYJPT5cNnQ9ZaBsAfVxkcj4PfeN2tHw+eDCa8eJ\n3ZV2iZrppggGkoqeRLHyrmLGCp3thVS4xj0fx6ODwe46M2kk0Zva8xqNOg6n9uOo6OPjD+Q+2w6U\nhj4+RfqlzzwqndA16sNLfpZ0y5trXx5sfg6LXMqlf14AsAr7B5ufb7pdqapjzF5Rtjaq7BVFUXKA\nGhOzZ5XPaq9xbAerSB1ylsdiysw5Aw8F+7bag4J9xZooZj+za4/g6z15TrCXT5FCw6070zInp4jX\n02Kmw28X+wlJ/kFvN3UwZa74LqYdom6nfWCpnhmKkQqXPZ5HNp9viuuCgTJlgSckAQknUzvX00aT\ndV1Gz56dRV0vf1KueMpkGeXMHCT9NWhNNCq6oZnE7LvTfvQP4Yxgt3Ll3bg8NT9HXjSXruhZdYnZ\nK0plospeURQlB9CXvaIoSg5QrSdoM21InqlOTp4LevACneOxJNiPQCo2xnZOatY38j1TGHxPDsoP\ndudhRcH+8smdgt2gZVQhsu5jwYV3qZrkyVSLZrarODmQYjCTKCJEkRSaygRtuS5wjXs+jkNBg911\nZlPYqCe1510KMR1I1Tfh5l+/XCH32XqYNPTJG/OD3fsZCXXdfkJfAPH0Vr8jFQAcQ0vBnkM7APK8\nAKCYlo+lq4cDyN9DWUI37NfNwJWajCp7RVGUHKBaTtB6f7rdq4C42kuXhrk/pfS9CKlP0AWPB/tu\nnB/soV9Ee6nOOoFq2N8pk4tv3tgi2IecKtUZ/el+pBnVA6lm/CJS0j2dYp4m85voTU+n8Dc6h5j4\nGKk0J5vTNPvS+aY5Md5fMiWxiNrWkewfqXTCDm7Y0OBU2X3rzcfk/jvfWRTsWRfSSOmLQgDAhN0u\nCb6TIEOG+Tgl2P75rCxDuiUgo7vKUui6UEqpTqiyVxRFyQGqZcze+zPF7HkJPat8n77HirEN7YbE\n9dPPwX3BnrpbtMJowFKpIbD4wqOD3WHUv4O9/mEZbdQ9MmrfDrRw6X1KdexIteRnu4Jk/fcW36T/\nip0vZqzQWX2kwiUU+LiZNDoY7K4zm0YSPak979OoYz+K68Op/PWvyn0eMkpGM4vHSr/0Wkq17dtE\nfcjprbyXwJF4Ndh+0VumdMt0JRKA7IvZK0o2ocpeURQlB6gxMXtW+SXtYLUvad9XcWSweck+L/IZ\ntCGK2T/YRooPnDlLliu9N7ZpsA8460O5Qa/oqUTCflRGuIgKj/VsG31Oe0l8FCrHbLJp4yh8hFTy\nyC4iuyfZ09yooX9b+i2VO86ndlKNsrBQrO7Z0t/vPSD332GWjHIe7CX9df6GaFQ0uY7E7I/Hc8Fe\nRKMq/3w+oNmHTKO1dAuoqjJmryjZiip7RVGUHKBax+wz7VGbuUDalwDiirEVZaU/h+ODzZtp3F3H\nxezfpph9L4rZXyNq9od75P+fOx7n9p6l0gPvnyl2PpUXnu1GAf2pHPIkSUCJlS3mQmd1kcr7ZPNx\nvOnJYHed2TSS6EnteZ9GHfs9SAe6Amk/PC/3ecA1MppZfLX0y5lvy+hn6sHRgZzx5DcpAYCWWBFs\nvw7CPy8A+AoNgp0uAweo3LIHm/oVJdtRZa8oipID6MteURQlBzDW2sq9oDH2ahvFCCqyNjhP1LG/\nDg35N7ghP+96xCl9vGR/BVoGu9/GQgDAU7Wk1sGZ8yREsbqrVJVvfj5Nmf7NfXYSFyi9sYgmQfNd\ngcdpVMOetqONLY7ak2wqUR/gEgq86Opksv2d9qf6+UV07XwK44AmbrHAfcpWslh9N93/PLn/B7vK\nBO1JG6MFVNNr9Q0+Dt3w4jafcsl7EaR7jkD6UJ6WPVA2ZQzGVnUTthhjDKy1przHq7JXFEXJAapU\n2TMVWRs80zl8Gian7nHpBFaSvJTfp2cOWCMTtEuayU6w7ca/HOzfqIDYNl2cQb5i2gEqjxZY+YJk\nJ28vvhk/ic0lELjQ2c5I5TuyufY9L7bq467zBF2jJ7WnmBZS5dF+tH6y+bf54tqGJqCXDKd+WSP9\nMrVZdOO8eOopyEiJR1t+0VumgmeZnm86n5Y9UABV9oAqe0VRlJygXMreGNMEwN2IwsMWwL+stRON\nMfURrf9pimgTpbOstd9scmyFxexLw3b4Odi/oDYAoCHWBd+XlNLHqnM19g12T7ekyZfeBYAzF1DM\nvhPFrAdQzN7H50kxg5R9EcW9891KqWmUN8lpk5nKGn+NVLiEAsf0udyxT9/sTyu0iuja+X+jH7Oy\n94qf4virp9L9L6CYfSeJ2ftS0rNpaVdzfBBsXtzWwKVcrkPD4Ev3HEuDlj1QAFX2QPmV/a8ALrfW\nHgSgLYBBxpgDAYwEsNBa2wLA0+7fiqIoShVTITF7Y8yjAP7p/mtvrV1njNkDQJG19oBNfltlMXu/\n5J4XWvEGGZ+TZs6nQgN+kc+ATylm34hi0xMlNm3PlXYYX++A1HwxZePk0UYms5/y1xVmkp1HNsfe\nd0AqP5LNsf5isnu7zyLy9aT2FNOGJXmcjeO6wNIKLXOv2EuGUL98SjH7RlEn8CK2Irrb3WkM4jcq\nyVTwrDJi9krNQpV9BcTsjTF5AFoDWAqgobXWx0jWATQOVxRFUaqMLSqXYIzZCcAcAJdaa78zRv6n\nY621xpi0w4ZnE7IVYNP8psjLb5ruZ4qiKDlLUVERioqKKux85Q7jGGO2BfB/AJ6w1t7ifO8ByLfW\nfmaMaQRgcWnDOJ6KnqBNdw5O6ePwAC/y4TTM7q5ODk8inrlEJmjXtGsU7GYDP5WL+5DOWdQgKmVZ\n9E+x893arilSwj22I9XLZPMEbEkTtPz9n8j2u8YOlDVlKKJr5/+VfsxFdR5wnxS6WTOF7n+J3P+D\n7WSC1k9+P0L1cDjdkhex+fBNWdItM6ELpRRAwzhAOcM4JpLwdwJY4V/0jrkA+ji7D4BHNz1WURRF\nqXzKm3p5HIAlAN5ElHoJAFciEqAPANgbpUi9LImKVl8+fY9T93g3JFaSvIPVajfN2W+9TNC+WJcm\nIu+gCVqqahkmaMlXTIo5j+ohzHb1Czjdch7ZPPlRTHZJE7R5ZK8j22turqDZk9pTTNth5dEIBK4C\nZmyClqpiLrlI+uWY9dIv0+tGE7TNaXrZ70gFpN93oLzplplQNZ+7qLIvZ8zeWvs8Mo8KOmbwK4qi\nKFVE1pRL2Fqki9NySh8X1WLVyQt6fOkEjiufuYxi9kdQzPoyitn7xVSUeslF0YoKxT52l+hz2rfi\ny6PDMtWwX49UMn3Po4Zi99l/F/G9QNfO70s/XkC2H9xQOYU1t9D9L6OY/RESs/fzIVwigRe3raYk\nUV/0rCzploqyOVTZa7kERVGUnKDGK3sm3bJ43g2JlWRrSL1fn5lzzsb7gu/1WlIbuF0h5cqcQhf0\nOz91EVcxlSHIo/LCM93lOAOHSyRwhg3H3rdFKr+SzbF+zszxsbZ3ydeb2lNM5Y7zqKwDfAG0+8kn\nmxI5P4AAAAspSURBVE9hSV+J2R+2UU5yX61zAMQzcJZDLsijLb8rlap5paJQZa/KXlEUJSfQl72i\nKEoOkFNhHA/vevQzpfRxiOEryM7f+VgMIL5Rea/3ZAnBmgNognIYTdD6ydh+dHFKb1xEKYuHu7yo\nmb+Jj3ek4nALp1ty7XoP17jnNEwOEfkdrHpTPtZrdO2OlC5KG3gB090nTdquuZHu/z25/1kHnBbs\nfd3kdxE6BN+ukB3VeRFbbZdyyZPnirIlaBhHlb2iKEpOkFPKPt0ELU8M/uz2MAWA/bEq2L4aZg88\nFHzvUhpmh4f+LRfhjWP7pfreHy32flRLvtDVks+jw7lEAqt5nmhNt1CCBHpsYpdVvp9GLSZfX2rP\n+1Tbfr8x9COv8qen8QFYfMbRwT6Qyk/MwRkA4tUtV6JFsGvTHsJ+olwnaJWKQpW9KntFUZScIKeU\nvSeTYuRFPpyGeQxeBBCPK/deMyfYq5vRTk0jaaeq9u6zL128rZhP0MIkH0/nmmOsyj8km9U8q/iS\nvufaon50QHXZYvMCJ/PuWi+RXeg+nxXX6nF0/2vk/mc26xFsPx/yIo4JPh5V8SK2dCMwRdkSVNmr\nslcURckJclLZM5kyc/KwJti+QNopYUURsIai6yfOe15OKFWQJWZ/hLjevU7sA2kz2WkulM2LoN4k\nu6QMnExkysw5xH3yAq3+1J53afPaA/lxLXOfHLOXrXuxsOtxwW5GMwKPu5VlXPCsGM2CXZuKnmkW\njlLRqLJXZa8oipIT6MteURQlB8j5MA7DIYYNFDjxdXJ4ErHf11In54P6NEE5iiZo/WQsbTLOK5se\npknOPPfJ9XA4BENRlXJDUZoQCuJ61MVkn96e/sEzt37zcZq0XT1W7n/fr+X+p9c/J9h+8pvr4dSh\nwBLvJaAoFY2GcVTZK4qi5AQ5r+wzpWHyDlZ+wrCDK5sAyG5KAHDKgqflhK3p5L7qJVWTfI0qSDbf\nXuyZP0WfXIt+NdklpVtmItNxvvAD17vvTe1Z/ZPYh1OlzlAMlKteLhfz8U5/Djb34WJXJoEnxLkP\ndQGVsjVRZa/KXlEUJSfIeWXPZErD3B8rAcTjyuf9fHewi2tLCuFBBR/ICZ2i/5UKoW1LJQlmUsqi\nnw2QnW/Ln26ZiXRpmG3Ix2mYvSmF9FcqnbCtT7mkevfvJPcNdt7PkrJ6T+3zg+3nQ1Zi/+DTdEul\nslBlr8peURQlJ1BlnwHOzPnFFUjzZRMA4EsqgXz6C0/IgVLbS2L2pOZfuF1sLmHsyySwmq+IDJxM\n+MwcXmjFpRM+JvvYi+kfXuVzzF5qxuHhY08OdgMqYezLJGxHBc80A0epLFTZq7JXFEXJCVTZE5ky\nQnwRL869PwNSCO1zymBvfQ0lpbuc+vXniauuJKBg2n/E9vH096g95c3AyUS68x1APp4X6L+P2Osl\nqQZ173EG3ebyq2XxAJcwfghSCM3n1HMxOc3AUSoLVfaq7BVFUXICfdkriqLkABrGycB2lBa40QVA\nWtPqof/R5GKfNx4I9m+ShYltujiDQiLz7hGbyxf4Mgnbko8nTysaH5D6lXxcOoEnh7tSGAou9PSb\nFADFNpJtiRmHnhXsP9Ik93K32qwWBaR+ofRWRdmaaBhHlb2iKEpOoMq+FPgdlXgSsQvVtucUwnb/\noJ1jncovPltcDXcUe8oPYntFz/vLVga8Gxar/IHUznXUzjyfcklqfsmVfwo2p6zOR5dg+8nY72JL\nuxSlclBlr8peURQlJ1BlXwZ8CiYQT8O8dO2/gr2+8XbBrnukW0DE6ZYU62aNyymXVUXGNMwu9A+X\nhrn+VbrPT2Sh1K1N/l+wuYQxp1wqSmWjyl6VvaIoSk6gyr6ctKGSZVzE68zx8+RHTtEX8aIqOgf9\nMuvoSjaXQc732US00OrB4fJrLia3NFZmTVGqDlX2quwVRVFyAn3ZK4qi5AAaxqkALsHkYH9P067N\nG7r9WHeR3yapNnx1oYCqduLb6GP1Otl3dieazr0NgyqpVYpSejSMsxWUvTGmszHmPWPM+8aYERV9\nfkVRFKXsVKiyN8bUArAS0cr7jwG8AuBsa+279Jsap+yZEwZfh3yvhN3CpGT/KmtOhdN+GJB/AABa\naDV6SM15nsVFHyIvv2lVN2OrUZPvb3P3psq+4pX9nwCsttYWW2t/BTALwKkVfI2spmh1yb+pzhSt\nrOoWbF0+LPqwqpuwVanJ91eT760i2Kbkn5SJPQFaeQR8BORW/t3iNsfh9yHtAAC1zHVV3JqKp2ge\nkJwHbKzBozNFqYlUtLKv3NleRVEUpVRUdMy+LYCEtbaz+/eVAH631l5Pv9H/ISiKopSDLYnZV/TL\nfhtEE7R/RrTG8mVsMkGrKIqiVD4VGrO31v5mjPkrgKcA1AJwp77oFUVRqp5KX1SlKIqiVD6VWi6h\npi24MsY0McYsNsa8Y4x52xgzxPnrG2MWGmNWGWMWGGPqlXSubMUYU8sYs9wYM8/9uybdWz1jzEPG\nmHeNMSuMMW1q2P1d6f423zLG3GeMqV2d788Yc5cxZp0x5i3yZbwfd//vu3dOp6ppdenJcH83uL/P\nN4wxDxtjdqHvynR/lfaydwuu/gmgM4CWAM42xhxYWdffSvwK4HJr7UEA2gIY5O5pJICF1toWAJ52\n/66uXApgBSTTqibd260A5ltrDwRwCKJtBWrE/Rlj8gAMAHC4tbYVorBqL1Tv+5uO6P3BpL0fY0xL\nAD0RvWs6A7jNGJPttcDS3d8CAAdZaw8FsArAlUD57q8yb77GLbiy1n5mrX3d2d8DeBfRWoNuAGa4\nn80AcFrVtHDLMMbsBaALgGkAfBZATbm3XQAcb629C4jmm6y136KG3B+iytS/AqjjEifqIEqaqLb3\nZ619DsD/NnFnup9TAdxvrf3VWlsMYDWid1DWku7+rLULrbW/u38uBeCLUpX5/irzZZ9uwdWelXj9\nrYpTUq0RPZCG1tp17qt1ABpWUbO2lJsBXAHgd/LVlHtrBuALY8x0Y8xrxpipxpgdUUPuz1r7NYCb\nAPwX0Uv+G2vtQtSQ+yMy3U9jRO8YT01431wAhM2vy3x/lfmyr7EzwcaYnQDMAXCptZZ39IONZsCr\n3b0bY/4C4HNr7XKIqo9RXe/NsQ2AwwHcZq09HFG1n1hIozrfnzFmXwCXAchD9GLYyRhzLv+mOt9f\nOkpxP9X2Xo0xowD8Yq29bzM/2+z9VebL/mMgthFpE8T/z1QtMcZsi+hFf4+19lHnXmeM2cN93wjA\n51XVvi3gGADdjDFrANwP4ARjzD2oGfcGRH97H1lrX3H/fgjRy/+zGnJ/RwJ40Vr7lbX2NwAPAzga\nNef+PJn+Hjd93+zlfNUOY0xfROHU3uQu8/1V5sv+VQD7GWPyjDHbIZpcmFuJ169wjDEGwJ0AVlhr\nb6Gv5gLo4+w+AB7d9Nhsx1p7lbW2ibW2GaKJvWesteehBtwbEM23AFhrjGnhXB0BvINot8hqf3+I\nJpvbGmN2cH+nHRFNtNeU+/Nk+nucC6CXMWY7Y0wzAPshWuRZrTDGdEYUSj3VWvsTfVX2+7PWVtp/\nAE5GtMJ2NYArK/PaW+l+jkMUz34dwHL3X2cA9QEsQjR7vgBAvapu6xbeZ3sAc51dY+4NwKGIynC/\ngUj57lLD7m84ov+BvYVo8nLb6nx/iEaYnwD4BdH8X7/N3Q+Aq9y75j0AJ1V1+8txfxcAeB/Ah/R+\nua2896eLqhRFUXKAbM87VRRFUSoAfdkriqLkAPqyVxRFyQH0Za8oipID6MteURQlB9CXvaIoSg6g\nL3tFUZQcQF/2iqIoOcD/B/+mtHszmu0KAAAAAElFTkSuQmCC\n", "text": [ "" ] } ], "prompt_number": 19 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Obviously, the stability condition is there for a reason :). As usual, implicit methods come to the rescue. Implicit differencing leads to\n", "$$\\frac{U^{n+1}_i-U^n_i}{\\Delta t} = D\\frac{U^{n+1}_{i+1}-2U^{n+1}_i+U^{n+1}_{i-1}}{\\Delta x^2}$$\n", "Substituting $\\alpha=\\frac{D\\Delta t}{\\Delta x^2}$ and collecting the advanced values on left side yields\n", "$$-\\alpha U^{n+1}_{i+1}+(2\\alpha + 1)U^{n+1}_i-\\alpha U^{n+1}_{i-1} = U^n_i$$\n", "Or expressing in vectorized form, where $K$ is the matrix of second differences:\n", "$$U^{n+1}_i-U^n_i = \\alpha K U^{n+1}$$\n", "$$({1}-\\alpha K) U^{n+1}=U^n$$" ] }, { "cell_type": "code", "collapsed": false, "input": [ "from scipy.sparse import dia_matrix, eye\n", "from scipy.sparse.linalg import factorized\n", "def d2matrix(nelem):\n", " elements = ones((3,nelem))\n", " elements[1,:] *= -2\n", " return dia_matrix((elements, [-1,0,1]), shape=(nelem,nelem)).tocsc()\n", "\n", "def diffusion_solve_implicit(U, dx, dt, nx, nt, D):\n", " alpha = D*dt/dx**2\n", " d2x = eye(nx)-d2matrix(nx)*alpha\n", " xint = arange(1, nx+1)\n", " LU = factorized(d2x.tocsc())\n", " for it in range(0,nt-1):\n", " U[it+1,xint] = LU(U[it,xint])\n" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 20 }, { "cell_type": "code", "collapsed": false, "input": [ "U, dx, dt, x, t = diffusion_init(maxx, maxt, D, nx, SC*10.)\n", "U[0,:] = 0.\n", "U[0, (xmaxx*0.2)] = 1.\n", "#U[0,:] = sin(x*5)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 25 }, { "cell_type": "code", "collapsed": false, "input": [ "diffusion_solve_implicit(U, dx, dt, nx, len(t), D)\n", "pcolormesh(U, rasterized=True, vmin=-1, vmax=1)" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 26, "text": [ "" ] }, { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEACAYAAACj0I2EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEp9JREFUeJzt3W2MXNV9x/Hfr7t4sXkwdaSSAK6NKEQJSlFQ2lLawDSh\nrYsiyItKBSUtASmvmobkBS0PavBGStUkikLUihdNAiJpoFIJQiAhFedhDBUSJY2hBJsa2kAwCJOa\nBxts1l773xdzZ+fsesa7c2eG8f79/UirOXPm3DvnzAw/ru88/B0RAgAsf78y7gkAAIaDQAeAJAh0\nAEiCQAeAJAh0AEiCQAeAJI4Y6LZvs73T9pNF31dtb7P9hO17bK8e/TQBAItZ7Aj9dkkbFvQ9KOnc\niDhP0nZJN4xiYgCA/hwx0CPiYUmvLejbFBGHqquPSjpjRHMDAPRh0HPo10h6YBgTAQAMpnag275J\n0v6IuHOI8wEA1DRZZyPbn5J0qaSPHmEMPxIDADVEhOts1/cRuu0Nkq6TdHlEvL3IpNL+3XzzzWOf\nA+tjbawv398gFvvY4l2SHpH0Xtsv2L5G0j9IOlHSJttbbN860AwAAENxxFMuEXFll+7bRjQXAMAA\n+KZoTY1GY9xTGKnM68u8Non1Hcs86Dmbnju2Y1T7BoCsbCveqTdFAQBHJwIdAJIg0AEgCQIdAJIg\n0AEgCQIdAJIg0AEgCQIdAJIg0AEgCQIdAJIg0AEgCQIdAJKoVbFoqb6gmzShg3PXJzRb3HGnf4X2\nHzZmqugrb5/SzGH9q7R3rq9sr9S+ufZJ2nNYe17f7s7YyV3FIl4t2m906dvV5faF7d1d+t4q2mWZ\nkJkufUtxfHU51aVPkk4o2quL9sld+sr2u4r2mi5jir7ZYuyek1d22jpp3uXC9j51xu7Vqrn2jFZU\nt5d9nQUe1MRce381VpJmq/6Dxcu7HFu2Z4t25/bF/7M42GU7jNcX9aVxT2HsOEIHgCQIdABIgkAH\ngCQIdABIgkAHgCQIdABIgkAHgCQIdABIgkAHgCQIdABIgkAHgCQIdABIgkAHgCSOGOi2b7O90/aT\nRd8a25tsb7f9oO1TRj9NAMBiFjtCv13ShgV910vaFBHnSPphdR0AMGZHDPSIeFjSawu6L5N0R9W+\nQ9LHRzAvAECf6pxDPzUidlbtnZJOHeJ8AAA1DfSmaESEpBjSXAAAA6hTgm6n7XdHxMu23yPplV4D\nN298SK7yfn3j1/UbjdO6jutWzutgUaKuV8mwiblSY91v71WirN2eV5aseCQmy0dlsUdotmgfLNoz\nRfvtBZdSfyXoyv2WJrqM6TWfcuxUl3Y5ttxHL5MLLjX/MexW8m0pZeDmP2dTi9x+eNm51pjJJd9f\nabHSc5Sdwyg0m001m82h7KtOoN8n6SpJX64u7+018OKNF82rKbq0pACAY0ej0VCj0Zi7Pj09XXtf\ni31s8S5Jj0h6r+0XbF8t6e8l/aHt7ZI+Ul0HAIzZEY/QI+LKHjddMoK5AAAGwDdFASAJAh0AkiDQ\nASAJAh0AkiDQASAJAh0AkiDQASAJAh0AkiDQASAJAh0AkiDQASAJAh0AkiDQASAJAh0AkiDQASCJ\nOhWL+lK3bNf8MmEreo6TNK8qUtmeLNorippw7X3PFLXYVk7snWtPTR5abHLddSslJ3XKzZVl597s\ncnu5Xa9ScqXyIWpvV5aXm+0xtnzWj68uy7mf3GMf3RT7mp3oHB+Uj+1M9fzNdCkD2GpPdd2u/drZ\n36WvNbXDy86VYyg7h2MNR+gAkASBDgBJEOgAkASBDgBJEOgAkASBDgBJEOgAkASBDgBJEOgAkASB\nDgBJEOgAkASBDgBJEOgAkETtQLd9g+2nbD9p+07bU4tvBQAYlVqBbnu9pE9LOj8iPqDWj7NeMbxp\nAQD6Vff30HdLOiBple2DklZJenFoswIA9K3WEXpEvCrpa5J+IeklSa9HxA+GOTEAQH9qHaHbPkvS\n5yStl/SGpH+1/YmI+F45bvPGh+ba6xrrdFbjjK7761ZNplclmYkuVWrK6jcTRYmdiXkVi/bPtfdX\n7YPa19nXZHF/k52xixasKasJldV9ulUvKisTLdbut2JRe/on9BjbrUqR1JnnUu6v230X+y0fw7LS\nT7viUK/KRL0qGbXHlPsqx3arUlS2+61SRHUijEOz2VSz2RzKvuqecvmQpEciYpck2b5H0oWS5gX6\nxRsvWrDZUpICAI4djUZDjUZj7vr09HTtfdX9lMvTki6wvdK2JV0iaWvtWQAABlb3HPoTkr4j6SeS\n/qvq/qdhTQoA0L+6p1wUEV+R9JUhzgUAMAC+KQoASRDoAJAEgQ4ASRDoAJAEgQ4ASRDoAJAEgQ4A\nSRDoAJAEgQ4ASRDoAJAEgQ4ASRDoAJAEgQ4ASRDoAJBE7Z/PraNXia+ybFynbFin1NiKYuyMjqws\nOzdZtPdqZbG/mcP6Vk3s7cxhqlOCrqxM1/XRKku3zfRov7XgcmH7jU5zX9W/r9j+QHkfheOK+ays\nKrqt7DG2Zwm6dsm6cr5L2UfVjuIp3T/ReabKx7ZdNq7s2z+v7Fz30nTtMf2UnWtNf+Kwsb1Qdg6Z\ncIQOAEkQ6ACQBIEOAEkQ6ACQBIEOAEkQ6ACQBIEOAEkQ6ACQBIEOAEkQ6ACQBIEOAEkQ6ACQBIEO\nAEnUDnTbp9i+2/Y221ttXzDMiQEA+jPIz+d+Q9IDEfGntifV+SFWAMAY1Ap026slfTgirpKkiJjV\nvF/1BgC80+qecjlT0i9t3277p7a/aXvVMCcGAOhP3VMuk5LOl/SZiHjM9i2Srpf0hXLQ5o0PzbXX\nNdZpfWNd3XkCQErNZlPNZnMo+6ob6Dsk7YiIx6rrd6sV6PNcvPGiJe2sLBXWLkc3O6802Iri9nol\nw8rSdPvU+sfElDql5vaq8w+MlVOdcnSTU4fKnRxutkf77aL9ZnVZnpQq2rt3d9p73m7Pp/tuS5PF\nDauq9kmdZerkHtsVVd46g3qto1Suv9rH3hM6/8grH8OyrFz78S77yrHzt+s81+3Sc+Xro7x9tsdr\noVvpOUrN4WjVaDTUaDTmrk9PT9feV61TLhHxsqQXbJ9TdV0i6anaswAADGyQT7n8laTv2V4h6X8k\nXT2cKQEA6qgd6BHxhKTfGuJcAAAD4JuiAJAEgQ4ASRDoAJAEgQ4ASRDoAJAEgQ4ASRDoAJAEgQ4A\nSRDoAJAEgQ4ASRDoAJAEgQ4ASRDoAJAEgQ4ASQzye+gDGUYFmf3zSu/0r6xitEIzc+1VU50KOicc\n/2Zng26PVlEhqNiF9FbRblckKqoUvfpqp72rqBC0p7rcV2x+oMvdStJxRXtldbm3qJQ0W+x3Tbnh\n8UV7dXVZzr1cU6lcf7WPfVNl5aGVc+09Oumwdtm3WJWiVv9UNZ3Oa6WsUtStMtHC8cCxhCN0AEiC\nQAeAJAh0AEiCQAeAJAh0AEiCQAeAJAh0AEiCQAeAJAh0AEiCQAeAJAh0AEiCQAeAJAh0AEhioEC3\nPWF7i+37hzUhAEA9gx6hXytpq6QYwlwAAAOoHei2z5B0qaRvSfLQZgQAqGWQI/SvS7pO0qEhzQUA\nMIBaFYtsf0zSKxGxxXaj17jNGx+aa69rrNP6xro6dwcAaTWbTTWbzaHsq24JugslXWb7UrWKkZ1s\n+zsR8RfloIs3XtT3jtvlw8rycL1KjbXLkR3sUZasn1Jk5f1Naf9c+6TVnRJ0U+3SbeV0ijJv88rO\nFSXmtKvqKvp2zHYfum/BpdRfCbrdRd++cm7Fnawp539Cdfmuoq/crkvZOUmaqUrXlWXlXtevFu1T\n5trtMWXZuX1FubqZopRgWY5utsvzR9k5ZNNoNNRoNOauT09P195XrVMuEXFjRKyNiDMlXSHpRwvD\nHADwzhrW59D5lAsAjFndUy5zImKzpM1DmAsAYAB8UxQAkiDQASAJAh0AkiDQASAJAh0AkiDQASAJ\nAh0AkiDQASAJAh0AkiDQASAJAh0AkiDQASAJAh0AkiDQASCJgX8+d1R6VaAp+9tVhsoqNrPzqhd1\n2t2q4iylys2KVZ3qRaevqUoPTRUD3i7au4r2K53mzpdaly8WN+8s2mWVoT3VZVmlqCwgVCqfvHb1\nopOKvrLq0WyxkwMvddqntqsQrSkGry3a5VqLMf+3qlXiaKdOnevbVZQ9Wqxi0fznZkXRf/hLkmpE\nwNJwhA4ASRDoAJAEgQ4ASRDoAJAEgQ4ASRDoAJAEgQ4ASRDoAJAEgQ4ASRDoAJAEgQ4ASRDoAJAE\ngQ4ASdQKdNtrbf/Y9lO2f2b7s8OeGACgP3V/PveApM9HxOO2T5T0n7Y3RcS2Ic4NANCHWkfoEfFy\nRDxetd+UtE3SacOcGACgPwOfQ7e9XtIHJT066L4AAPUNFOjV6Za7JV1bHakDAMakdgk628dJ+r6k\nf46Ie7uN2bzxobn2usY6rW+sq3t3XXUrTXZwXgm6zvImikJu7RJ0ZUm0PUXxtrK9Tys7O1+7XZJ0\n+ktFrbmi1Jye7jSf+d9O+/nqslfZuXml4qrLsgTdUrRL0JX7fXUJ7d3VPM8ud1aePPu1TvPFtZ0S\nc9t1jqT5JejKsnPlY9t+vGd7PDeUmMOxrNlsqtlsDmVftQLdtiV9W9LWiLil17iLN15Ud14AcExo\nNBpqNBpz16enp2vvq+4pl9+T9ElJf2B7S/W3ofYsAAADq3WEHhH/Lr6UBABHFUIZAJIg0AEgCQId\nAJIg0AEgCQIdAJIg0AEgCQIdAJIg0AEgCQIdAJIg0AEgCQIdAJIg0AEgCQIdAJIg0AEgCQIdAJJw\nRIxmx3b8bdw4kn0fTSb8d+OewkgdPAaeQ+TwRX1p3FMYCtuKCNfZliN0AEiCQAeAJAh0AEiCQAeA\nJAh0AEiCQAeAJAh0AEiCQAeAJAh0AEiCQAeAJAh0AEiCQAeAJGoHuu0Ntp+2/YztvxnmpAAA/asV\n6LYnJP2jpA2S3i/pStvvG+bEjnbPNZ8f9xRG6ufjnsAIZX/uWN+xq+4R+m9LejYinouIA5L+RdLl\nw5vW0e/55C+q58Y9gRHK/tyxvmNX3UA/XdILxfUdVR8AYEzqBvpoqmIAAGqrVbHI9gWSNkbEhur6\nDZIORcSXizGEPgDUULdiUd1An5T035I+KuklSf8h6cqI2FZnEgCAwU3W2SgiZm1/RtK/SZqQ9G3C\nHADGa2RFogEA76yRfFM005eObK+1/WPbT9n+me3PVv1rbG+yvd32g7ZPGfdcB2F7wvYW2/dX19Os\nz/Yptu+2vc32Vtu/k2x9N1Svzydt32l7armuz/ZttnfafrLo67mWau3PVHnzR+OZ9dL1WN9Xq9fm\nE7bvsb26uK2v9Q090BN+6eiApM9HxLmSLpD0l9V6rpe0KSLOkfTD6vpydq2krep8ginT+r4h6YGI\neJ+k35T0tJKsz/Z6SZ+WdH5EfECtU6BXaPmu73a1sqPUdS223y/pz9TKmQ2SbrV9tP+cSbf1PSjp\n3Ig4T9J2STdI9dY3isWn+tJRRLwcEY9X7TclbVPrM/eXSbqjGnaHpI+PZ4aDs32GpEslfUtS+931\nFOurjnY+HBG3Sa33fyLiDSVZn6Tdah10rKo+rLBKrQ8qLMv1RcTDkl5b0N1rLZdLuisiDkTEc5Ke\nVSt/jlrd1hcRmyLiUHX1UUlnVO2+1zeKQE/7paPqaOiDaj3op0bEzuqmnZJOHdO0huHrkq6TdKjo\ny7K+MyX90vbttn9q+5u2T1CS9UXEq5K+JukXagX56xGxSUnWV+m1ltPUype2DFlzjaQHqnbf6xtF\noKd8l9X2iZK+L+naiNhT3hatd5aX5bptf0zSKxGxRZ2j83mW8/rU+iTX+ZJujYjzJb2lBacflvP6\nbJ8l6XOS1qsVACfa/mQ5Zjmvb6ElrGXZrtP2TZL2R8SdRxh2xPWNItBflLS2uL5W8/8vs+zYPk6t\nMP9uRNxbde+0/e7q9vdIemVc8xvQhZIus/1zSXdJ+ojt7yrP+nZI2hERj1XX71Yr4F9Osr4PSXok\nInZFxKykeyT9rvKsT+r9WlyYNWdUfcuO7U+pddrzE0V33+sbRaD/RNLZttfbXqHWSf37RnA/7wjb\nlvRtSVsj4pbipvskXVW1r5J078Jtl4OIuDEi1kbEmWq9mfajiPhz5Vnfy5JesH1O1XWJpKck3a8E\n61PrDd4LbK+sXquXqPXmdpb1Sb1fi/dJusL2CttnSjpbrS85Liu2N6h1yvPyiHi7uKn/9UXE0P8k\n/Yla3yR9VtINo7iPd+pP0u+rdW75cUlbqr8NktZI+oFa70o/KOmUcc91CGu9WNJ9VTvN+iSdJ+kx\nSU+odQS7Otn6/lqt/0k9qdabhsct1/Wp9a/ElyTtV+u9uKuPtBZJN1Y587SkPx73/Gus7xpJz0h6\nvsiXW+uujy8WAUASR/tnNgEAS0SgA0ASBDoAJEGgA0ASBDoAJEGgA0ASBDoAJEGgA0AS/w/4FVmf\n7fohlwAAAABJRU5ErkJggg==\n", "text": [ "" ] } ], "prompt_number": 26 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Much longer time steps can now be used without compromising stability. In analogy to ODE methods, the previous methods correspond to explicit and implicit Euler methods. We may develop a method which is second order in time and space by proper time-centering:\n", "$$\\frac{U^{n+1}_i-U^n_i}{\\Delta t} = \\frac{D}{2}\\left(\n", "\\frac{U^{n+1}_{i+1}-2U^{n+1}_i+U^{n+1}_{i-1}}{\\Delta x^2}+\n", "\\frac{U^{n}_{i+1}-2U^{n}_i+U^{n}_{i-1}}{\\Delta x^2}\\right)\n", "$$\n", "Again, substituting $\\alpha=\\frac{D\\Delta t}{\\Delta x^2}$ and collecting the advanced values on left side yields\n", "$$-\\alpha U^{n+1}_{i+1}+(2\\alpha + 2)U^{n+1}_i-\\alpha U^{n+1}_{i-1} = -\\alpha U^{n}_{i+1}+(2\\alpha + 2)U^{n}_i-\\alpha U^{n}_{i-1}$$\n", "Or expressing in vectorized form, where $K$ is the matrix of second differences:\n", "$$U^{n+1}_i-U^n_i = \\frac{\\alpha}{2} K (U^{n+1}+U^n)$$\n", "$$({1}-\\frac{\\alpha}{2} K) U^{n+1}=(1+\\frac{\\alpha}{2}K) U^n$$" ] }, { "cell_type": "code", "collapsed": false, "input": [ "def diffusion_solve_CN(U, dx, dt, nx, nt, D):\n", " alpha2 = D*dt/dx**2*0.5\n", " M1 = eye(nx)-d2matrix(nx)*alpha2\n", " M2 = eye(nx)+d2matrix(nx)*alpha2\n", " LU = factorized(M1.tocsc())\n", " for it in range(0,nt-1):\n", " U[it+1,1:-1] = LU(M2.dot(U[it,1:-1]))" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 27 }, { "cell_type": "code", "collapsed": false, "input": [ "diffusion_solve_CN(U, dx, dt, nx, len(t), D)\n", "pcolormesh(U, rasterized=True, vmin=-1, vmax=1.0)" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 28, "text": [ "" ] }, { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEACAYAAACj0I2EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEzJJREFUeJzt3XuMXOV9xvHn6ZpdbPAFk0AuNjaiEJE0SUE0IWkSpglt\nXYQgrSoVFFIS1Kitmobkj7RcRGIjETWJoiRqRS8EEE0ClQIIEQmpOJcxtDQkaQwBbC6puCNM4kts\nsL3rtX/9Y87MeXc9s7tzZsbjff39SKv5zTtnzrzv7PL4cObyc0QIADD//cawJwAA6A8CHQAyQaAD\nQCYIdADIBIEOAJkg0AEgEzMGuu2bbG+x/Ugy9mXbm20/bPtO20sHP00AwGxmO0K/WdKaaWP3Snpb\nRLxT0pOSrhzExAAA3Zkx0CPifknbp42tj4gDxdUHJa0Y0NwAAF3o9Rz6ZZLu6cdEAAC9qRzotq+W\nNBERt/ZxPgCAihZUuZPtj0k6T9KHZtiGL4kBgAoiwlXu1/URuu01kj4r6cKI2DvLpLL9+fznPz/0\nObA+1sb68vvpxWxvW7xN0gOS3mL7eduXSfpHScdKWm97o+3re5oBAKAvZjzlEhEXtxm+aUBzAQD0\ngE+KVlSr1YY9hYHKeX05r01ifUcy93rOpuOO7RjUvgEgV7YVh+pFUQDA4YlAB4BMEOgAkAkCHQAy\nQaADQCYIdADIBIEOAJkg0AEgEwQ6AGSCQAeATBDoAJAJAh0AMlGpY9FcfU5Xa0T7W9dHNJk8cDk+\nqolWPabxYtv2ty/S7mR8vBjb0/b2hUm9WLuS+tXG5f5k7NflY3hbsohfJ/UrOtjmNmOStDWpdxaX\nryZjryX1eJt6MhlL69SCNvVYMpbWxyT1sUm9pLg8vsNjnNrh8U4oLpeWQ7G8rHctHS3rkcWNy+SB\nd2lxq96jRa1695R6oSRpIlnIuMr9puMTyfikRiRJ+5MJ7y/GpteTSV3e3v4/i/1ttsXh41pdN+wp\nDB1H6ACQCQIdADJBoANAJgh0AMgEgQ4AmSDQASATBDoAZIJAB4BMEOgAkAkCHQAyQaADQCYIdADI\nBIEOAJmYMdBt32R7i+1HkrHlttfbftL2vbaXDX6aAIDZzHaEfrOkNdPGrpC0PiJOk/T94joAYMhm\nDPSIuF/S9mnDF0i6pahvkfThAcwLANClKufQT4yILUW9RdKJfZwPAKCinl4UjYiQFH2aCwCgB1Va\n0G2x/YaIeNn2G9W+MZskacPa++Qi71fXTtJv1t7UdrvZW4KVY2kLsmaburT92FjSz21qC7KD25Ht\nH+nQUqxTp7FmG7e07dwJSf1UUqdt5SbbjHVqQbe3NcmD7z9d+ttrzrlT67p026PbPHbaai9tO/dS\nUp/eZg4dnqv0ud3ftiXc7O3hmtt0ajvX6e+meb+5tJ3r1G6u3f2AQajX66rX633ZV5VAv1vSpZK+\nWFze1WnDc9Z+YEpv0M7JBABHplqtplqt1rq+bt26yvua7W2Lt0l6QNJbbD9v++OS/kHS79t+UtIH\ni+sAgCGb8Qg9Ii7ucNO5A5gLAKAHfFIUADJBoANAJgh0AMgEgQ4AmSDQASATBDoAZIJAB4BMEOgA\nkAkCHQAyQaADQCYIdADIBIEOAJkg0AEgEwQ6AGSCQAeATFTpWNSVblt4lS3Gyr5sIx1ahjVbz40k\nnZDGkxZlabekiSltzEYPmttE0pZtLG0Pl2q2aTspGftpUh/dZlupbDfXqe1cOr63ze1zaUHXXPb+\ndhtOe4x290vn/nRSn5XU6Zqa7fiSfaXP4ZTntni+2/0OpKm/s/G2v6f27erSbdu2GKTtHI4wHKED\nQCYIdADIBIEOAJkg0AEgEwQ6AGSCQAeATBDoAJAJAh0AMkGgA0AmCHQAyASBDgCZINABIBMEOgBk\nonKg277S9mO2H7F9q+2x2e8FABiUSoFue7WkT0g6MyLeLmlE0kX9mxYAoFtVvw99p6R9khbZ3i9p\nkaQX+zYrAEDXKh2hR8Q2SV+R9JyklyTtiIjv9XNiAIDuVDpCt32KpE9LWq1GH5vv2P5IRHw73W7D\n2vta9araKp1SW9G6PnvXmPL28Rm2k6Z2JlqQ1Gkno9GkK86oJor9lmMLR3a36jEdKHeeTnN5cflg\nMnZqUqfdi9JJN7sQpV2D2nU0Su+3NxmbS8eiZreguXQ6Sl/taG6/NRlLuxRtTup3d3js5sONlMcH\nU7sQjbUZa9+9aKLN/aZ2N5q5S1Faz6VLEd2JMGz1el31er0v+6p6yuUsSQ9ExFZJsn2npPdKmhLo\n56z9wLS7deqPBgBHplqtplqt1rq+bt26yvuq+i6XxyWdbXuhbUs6V9KmyrMAAPSs6jn0hyX9uxon\nGX5eDP9bvyYFAOhe1VMuiogvSfpSH+cCAOgBnxQFgEwQ6ACQCQIdADJBoANAJgh0AMgEgQ4AmSDQ\nASATBDoAZIJAB4BMEOgAkAkCHQAyQaADQCYIdADIBIEOAJmo/PW5VaTtvtK2ce3bg422qdq3o0tb\nzY1MaUFX1hNF27lGPTrlUpImxsrWZ8ccvaft/LWtuExbtN2Z1EuTOm3p9tq0y+l12o6uWOCepAXd\nvg6Nno5KuqctbD4xY203nfqbTrc5prg8IRn7UVL/SVJvS+pji8ujy6H0OZxo02JutlZzjfrg+3Vq\nO5fub7JNKznazuFIwxE6AGSCQAeATBDoAJAJAh0AMkGgA0AmCHQAyASBDgCZINABIBMEOgBkgkAH\ngEwQ6ACQCQIdADJBoANAJioHuu1ltm+3vdn2Jttn93NiAIDu9PL1uV+XdE9E/KntBSq/iBUAMASV\nAt32Uknvj4hLJSkiJjX1W70BAIdY1VMuJ0v6pe2bbf/M9g22F/VzYgCA7lQ95bJA0pmSPhkRP7H9\nNUlXSPpcutGGtfe16lW1VVpdW1V1ngCQpXq9rnq93pd9OSK6v5P9Bkn/ExEnF9ffJ+mKiDg/2Sau\niavmtL+0VVw5VraVW9Chrdxo0a9tLGkvt1C7W/Uila3kFmvXQfXx+lVrbJl2tOrX7SzPHi14JZlU\ns63cPcnYbyd1Op6egGq2bmu3L0l7km13Fq3p0iZ4+9TeUUm9sLhckrySsTBtiXd8Uqft5pYXl+m2\n5yX1Qx3Gi/1NJvv61ZJyJzu0rFVv1eskSbu0uDWW1rtbs5f2qPwfvfFW67qyRV3aPi5tO9eu3Ryt\n5o4s1+q6YU+hL2wrIlzlvpVOuUTEy5Ket31aMXSupMeq7AsA0B+9vMvlbyV92/aopP+T9PH+TAkA\nUEXlQI+IhyX9Th/nAgDoAZ8UBYBMEOgAkAkCHQAyQaADQCYIdADIBIEOAJkg0AEgEwQ6AGSCQAeA\nTBDoAJAJAh0AMkGgA0AmCHQAyASBDgCZ6OX70Pumm84ykxW70KQdkJr1aNI1J+16ND5WdjdaMHmg\n3Emzs9DFyY6vTeqTdPC2aZ12JtpW1tv2JuPFZTmDuXUsaq5kz2vl2PKkEdSS9Dc9ltTNjkVpR6Pv\nJHXaVPC5pC6aE42PlccEabehdh2J0s5Es3UpkspORelYKu1SRHcigCN0AMgGgQ4AmSDQASATBDoA\nZIJAB4BMEOgAkAkCHQAyQaADQCYIdADIBIEOAJkg0AEgEwQ6AGSCQAeATPQU6LZHbG+0/d1+TQgA\nUE2vR+iXS9okKfowFwBADyoHuu0Vks6T9A1J7tuMAACV9HKE/lVJn5V0YLYNAQCDV6ljke3zJb0S\nERtt1zptt2Htfa16VW2VVtdWVXk4AMhWvV5XvV7vy74c0f3pb9tfkPRRSZOSjpa0RNIdEfHnyTZx\nTVzVl0mmRrQ/qRut5BYkY6NJK7kxjbfqhdrTqpdpuyTpOO1ojR2vrW3r1+0u67GniuJfkgldktT/\nmtRPJXWxi32vlEMvJu3otiSbNjvT7UzGJtVe+q/xkuJyeTJ2YlK/eWlZH3VCckOz9dypydhfJvW3\nkvqvynK82P5Xi8redVvVvt6uZZKkHTquNbYnaUE3nvTEm0jazTXbDdJqDnNxra4b9hT6wrYiotJp\n7EqnXCLiqohYGREnS7pI0g/SMAcAHHr9eh8673IBgCGrdA49FREbJG3ow1wAAD3gk6IAkAkCHQAy\nQaADQCYIdADIBIEOAJkg0AEgEwQ6AGSCQAeATBDoAJAJAh0AMkGgA0AmCHQAyASBDgCZINABIBM9\nf33uoZZ2rNnf6mhTdizq1P2m3fhEMpZarF2teuy15Ibiobf/c9lt57h3lJ2Q9FvJtsck9XONi2eS\nLkXPJjenHYuanYqSvc7JtmmXkpJVSPuSxz41ndtJxWXaFumvy3L7z5O1PlrOqvm8TCwqn+O0S9GW\npF/SjqJj0W4tao11+t2k6E4EdIcjdADIBIEOAJkg0AEgEwQ6AGSCQAeATBDoAJAJAh0AMkGgA0Am\nCHQAyASBDgCZINABIBMEOgBkgkAHgExUCnTbK23/0PZjth+1/al+TwwA0J2qX5+7T9JnIuIh28dK\n+l/b6yNicx/nBgDoQqUj9Ih4OSIeKupXJW2W9KZ+TgwA0J2ez6HbXi3pDEkP9rovAEB1PQV6cbrl\ndkmXF0fqAIAhqdyCzvZRku6Q9K2IuKvdNhvW3teqV9VWaXVtVdWHm1GnVmVpy7O0zdluNdqqNVuj\nSTO0T3v9E636jC2NlwiOe0/Zim38R+Xjjf1FWe/877L+8d7G5QvJ3NJWcWn3t31txuaiOaO07Vza\n2i597GdfKut3FRNZ8sflWLqmdK26oSw3vv50SdITektrLH0Od2lxq55stQos/9xoLwc01Ot11ev1\nvuyrUqDbtqQbJW2KiK912u6ctR+oOi8AOCLUajXVarXW9XXr1lXeV9VTLr8r6RJJv2d7Y/GzpvIs\nAAA9q3SEHhH/JT6UBACHFUIZADJBoANAJgh0AMgEgQ4AmSDQASATBDoAZIJAB4BMEOgAkAkCHQAy\nQaADQCYIdADIBIEOAJkg0AEgEwQ6AGSCQAeATDgiBrNjO66Jqway78PJiL8w7CkM1P4j4HeIPFyr\n64Y9hb6wrYhwlftyhA4AmSDQASATBDoAZIJAB4BMEOgAkAkCHQAyQaADQCYIdADIBIEOAJkg0AEg\nEwQ6AGSCQAeATFQOdNtrbD9u+ynbf9/PSQEAulcp0G2PSPonSWskvVXSxbZP7+fEDnfP1J8d9hQG\n6ulhT2CAcv/dsb4jV9Uj9HdJ+kVEPBMR+yT9h6QL+zetw9+zmf9RPTPsCQxQ7r871nfkqhrob5b0\nfHL9hWIMADAkVQN9MF0xAACVVepYZPtsSWsjYk1x/UpJByLii8k2hD4AVFC1Y1HVQF8g6QlJH5L0\nkqQfS7o4IjZXmQQAoHcLqtwpIiZtf1LSf0oakXQjYQ4AwzWwJtEAgENrIJ8UzelDR7ZX2v6h7cds\nP2r7U8X4ctvrbT9p+17by4Y9117YHrG90fZ3i+vZrM/2Mtu3295se5Ptd2e2viuLv89HbN9qe2y+\nrs/2Tba32H4kGeu4lmLtTxV58wfDmfXcdVjfl4u/zYdt32l7aXJbV+vre6Bn+KGjfZI+ExFvk3S2\npL8p1nOFpPURcZqk7xfX57PLJW1S+Q6mnNb3dUn3RMTpkt4h6XFlsj7bqyV9QtKZEfF2NU6BXqT5\nu76b1ciOVNu12H6rpD9TI2fWSLre9uH+dSbt1nevpLdFxDslPSnpSqna+gax+Kw+dBQRL0fEQ0X9\nqqTNarzn/gJJtxSb3SLpw8OZYe9sr5B0nqRvSGq+up7F+oqjnfdHxE1S4/WfiPi1MlmfpJ1qHHQs\nKt6ssEiNNyrMy/VFxP2Stk8b7rSWCyXdFhH7IuIZSb9QI38OW+3WFxHrI+JAcfVBSSuKuuv1DSLQ\ns/3QUXE0dIYaT/qJEbGluGmLpBOHNK1++Kqkz0o6kIzlsr6TJf3S9s22f2b7BtvHKJP1RcQ2SV+R\n9JwaQb4jItYrk/UVOq3lTWrkS1MOWXOZpHuKuuv1DSLQs3yV1faxku6QdHlE7Epvi8Yry/Ny3bbP\nl/RKRGxUeXQ+xXxenxrv5DpT0vURcaak1zTt9MN8Xp/tUyR9WtJqNQLgWNuXpNvM5/VNN4e1zNt1\n2r5a0kRE3DrDZjOubxCB/qKklcn1lZr6r8y8Y/soNcL8mxFxVzG8xfYbitvfKOmVYc2vR++VdIHt\npyXdJumDtr+pfNb3gqQXIuInxfXb1Qj4lzNZ31mSHoiIrRExKelOSe9RPuuTOv8tTs+aFcXYvGP7\nY2qc9vxIMtz1+gYR6D+VdKrt1bZH1Tipf/cAHueQsG1JN0raFBFfS266W9KlRX2ppLum33c+iIir\nImJlRJysxotpP4iIjyqf9b0s6XnbpxVD50p6TNJ3lcH61HiB92zbC4u/1XPVeHE7l/VJnf8W75Z0\nke1R2ydLOlWNDznOK7bXqHHK88KI2Jvc1P36IqLvP5L+SI1Pkv5C0pWDeIxD9SPpfWqcW35I0sbi\nZ42k5ZK+p8ar0vdKWjbsufZhredIuruos1mfpHdK+omkh9U4gl2a2fr+To1/pB5R40XDo+br+tT4\nv8SXJE2o8Vrcx2dai6Sripx5XNIfDnv+FdZ3maSnJD2b5Mv1VdfHB4sAIBOH+3s2AQBzRKADQCYI\ndADIBIEOAJkg0AEgEwQ6AGSCQAeATBDoAJCJ/weq3nFCqunzOAAAAABJRU5ErkJggg==\n", "text": [ "" ] } ], "prompt_number": 28 }, { "cell_type": "markdown", "metadata": {}, "source": [ "What about other boundary conditions than 0 Dirichlet? In explicit schemes, the implementation is trivial. Otherwise, we need to put the BC to the RHS. For example nonzero Dirichlet BC conditions result:" ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 28 }, { "cell_type": "code", "collapsed": false, "input": [ "def diffusion_solve_implicit(U, dx, dt, nx, nt, D):\n", " alpha = D*dt/dx**2\n", " d2x = eye(nx)-d2matrix(nx)*alpha\n", " LU = factorized(d2x.tocsc())\n", " for it in range(0,nt-1):\n", " b = U[it,1:-1]\n", " b[[0,-1]] += alpha*U[it+1,[0,-1]]\n", " #U[it+1,1:-1] = LU(U[it,1:-1])\n", " U[it+1,1:-1] = LU(b)\n", "\n", "def diffusion_solve_CN(U, dx, dt, nx, nt, D):\n", " alpha2 = D*dt/dx**2*0.5\n", " M1 = eye(nx)-d2matrix(nx)*alpha2\n", " M2 = eye(nx)+d2matrix(nx)*alpha2\n", " LU = factorized(M1.tocsc())\n", " print(U[0,0]+U[0+1,0])\n", " for it in range(0,nt-1):\n", " b = M2.dot(U[it,1:-1])\n", " b[[0,-1]] += alpha2*(U[it,[0,-1]]+U[it+1,[0,-1]])\n", " U[it+1,1:-1] = LU(b)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 31 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Model of heat conduction across a rod. Note the \"huge\" time step" ] }, { "cell_type": "code", "collapsed": false, "input": [ "U, dx, dt, x, t = diffusion_init(maxx, maxt*100, D, nx, SC*40)\n", "U[:,0] = 1." ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 32 }, { "cell_type": "code", "collapsed": false, "input": [ "diffusion_solve_CN(U, dx, dt, nx, len(t), D)\n", "pcolormesh(U, rasterized=True, vmin=-0, vmax=1)" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "2.0\n" ] }, { "metadata": {}, "output_type": "pyout", "prompt_number": 33, "text": [ "" ] }, { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAEACAYAAABS29YJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXuwZFd13r/FndFo9EIWmNEzlsBSiqGwBcaCRE64SmQi\nOykJ/8OjgiMCSTmFA4RyHCSlymKcimKcQiFxCpdjHpZxUKwCrBrKUNaI6KYgKRDEEggGGQ3RJBqh\nucOMGGlGd+a+ZuWPPt29Tt+1ez96n1ef9avquvvus88+53Tf+53d3157HWJmGIZhGPPNi5o+AcMw\nDKN6TOwNwzB6gIm9YRhGDzCxNwzD6AEm9oZhGD3AxN4wDKMHTBV7IjqbiL5ORI8S0X4i+ndF/UVE\ntI+Ivk9EDxDRhWKf24noCSJ6nIjeVPUFGIZhGH7IF2dPROcw8woRbQPwVQD/EsDNAI4y8+8S0QcB\n/AQz30ZEuwF8BsDPA7gMwIMArmHmM5VehWEYhjEVr43DzCtF8SwACwB+jIHY31PU3wPgzUX5FgD3\nMvM6Mx8EcADAdTlP2DAMw4jHK/ZE9CIiehTAMoCHmPm7AHYx83LRZBnArqJ8KYBDYvdDGIzwDcMw\njAbZ5mtQWDDXEtGLAfwFEd0wsZ2JaJoXZPkYDMMwGsYr9kOY+Tki+nMAPwdgmYguZubDRHQJgCNF\ns6cBXCF2u7yoK+G5ORiGYRgOmJlS9vNF47x0GGlDRDsB/CKARwDsBXBr0exWAPcX5b0A3kZEZxHR\nVQCuBvCw44Tn9nXnnXc2fg52bXZ9dn3z95oF38j+EgD3ENGLMLgxfJqZv0xEjwC4j4jeDeAggLcU\nAr6fiO4DsB/ABoD38KxnaBiGYczMVLFn5scAvFapfxbAjY597gJwV5azMwzDMLJgK2grYHFxselT\nqIx5vjbArq/rzPv1zYJ3UVUlByUyd8cwDCMSIgJXMUFrGIZhzAcm9oZhGD3AxN4wDKMHmNgbhmH0\nABN7wzCMHmBibxiG0QNM7A3DMHqAib1hGEYPMLE3DMPoASb2hmEYPcDE3jAMoweY2BuGYfQAE3vD\nMIweYGJvGIbRA0zsDcMweoCJvWEYRg8wsTcMw+gBJvaGYRg9wMTeMAyjB5jYG4Zh9AATe8MwjB5g\nYm8YhtEDTOwNwzB6gIm9YRhGDzCxNwzD6AEm9oZhGD1g27SNRHQFgD8G8DIADOC/MPN/IqIPAfgn\nAH5UNL2Dmb9U7HM7gHcB2ATwPmZ+QOt7DxF2it/PF2VZf8GUusn687XtZ4v9zhUbdojysP5spW6y\nHLOf1hYAzvNs3zFjW9lGbGex35o455Vzxu/iKs4CAJzCOVvqJutXRHmtaLMSsN9qcVKnsHNL3aCP\ncf2aUl+uG/fr68+135o4T981rZ0RbU+Kazo1qD/zwrgOp2lcPjEuYlWUTw7bKnUh9a62qcc7Jcov\nBG6f7I+HhXXHCa046oedP+84oKv+hFInywDzneg7U8Ueg0/rA8z8KBGdB+B/E9E+DD7Ou5n5btmY\niHYDeCuA3QAuA/AgEV3DzGcqOHfDMAwjkKliz8yHARwuyieJ6HsYiDgAkLLLLQDuZeZ1AAeJ6ACA\n6wB8Ld8pl/HdreaeHr4Bm1ho+hS6xYZSt67UAYPv48ZcEuzZE9GVAF6DsXC/l4i+RUSfIKILi7pL\nARwSux3C+ObQfhbEyxixiW2jl9FiNsTLMCYIEvvCwvksgPcz80kAvw/gKgDXAngGwEem7M5TthmG\nkYKJuhGJd6hGRNsBfA7AnzDz/QDAzEfE9o8D+ELx69MArhC7X17UbeEhANuL8ssxuGvMHTYQnpne\nWDYuW8UEvdcsLS1haWkpS1++aBwC8AkA+5n5o6L+EmZ+pvj1VwA8VpT3AvgMEd2NgX1zNYCHtb5v\nQDmSxmgHTVk1vRF1wATcCGZxcRGLi4uj3/fs2ZPcl+8/+3oA7wDwbSJ6pKi7A8DbiehaDCyaJwH8\nGgAw834iug/Afgz+pN/DzGbj5MC+JXQPE3WjRfiicb4K3df/0pR97gJw14znZcwxXR7Fb2zo535m\n0+7GRruxFbSGYRg9wIYjPWVjwe7zPjY6/A3EMCYxsc9Bbk2Yk0+lK3H5XbaVDCOUbvw3GkYGTNSN\nPtNfse/Lldd0nW2zPEzYDaOMGbeGYRg9oC/jW2NOsRG8YYRhYt8kud/9GnSvSbumKxO+htFGzMYx\nDMPoATZUiqUn71iT9oiN4A0jP/ZfZTRG2yJ4DGOeMbHvIvaptZcN7QFucD8BSkuW5kqg5kqDbBgB\nmGzURYvf6Tosmy6P4l3vz6YjKVqt2OMFjUBaLEGGUR1OAe/wTSkIS7vcW0zse0QdWXjbKJY24evB\nbgC9oDP/Bdv9TfT92qc9cXT9/A3DaAWdEfvaMZH10tQo3uX/d3leoJXYhPBcYWKfm7a8o6Z7hmEI\n2iJNRotpdoGV3bVahUX5dBYT+7bQ80/CRN0wqqXnEjMHZP4E6xBdE/Y5hJs+AcOHJUIzDMPoATay\nn4UOv3ub26o/+baN4GOjeLT4/M0Nx/umraZ1xa9r9THpFAwjgQ7LVYdplwZmoW0Ll9p2o8lOJ24C\nnTjJ3tCu/9Cm6MK7MOfalcrci7phZKILMmfMQvYJ3Ob+ZEzYPVi2TGMKJvZG6zBRN4z8mNi3kRZ8\nKnWkHmibz98rbHFU75ir/7bUZGmVUfe72+JP04TdGGMTt00wNc6eiK4gooeI6LtE9B0iel9RfxER\n7SOi7xPRA0R0odjndiJ6gogeJ6I3VX0BtbBNvOaQDSyMXn1hEwujl69+88zC6CXZWF8YvcaV4tVF\nun7+hhPfoqp1AB9g5lcBeAOAXyeiVwK4DcA+Zr4GwJeL30FEuwG8FcBuADcB+BgR2cKtHtHGG4cm\n6rXjElETVqMmpgoxMx9m5keL8kkA3wNwGYCbAdxTNLsHwJuL8i0A7mXmdWY+COAAgOsqOO/5YEG8\nathvc9v4lQPX6Dgn2s1jE9tGr/L56PWdo22j67adj5FE8KibiK4E8BoAXwewi5mXi03LAHYV5UsB\nHBK7HcLg5jDfzLnNY7SETfEyjEiC5ImIzgPwOQDvZ+YTRDTaxsxMRNPSIKnbHsJ4QvXlAK4NOt0Z\naY+zEE7qDSTjjSf3yF3rr022T6eYK+G3BQGTLC0tYWlpKUtfXkkgou0YCP2nmfn+onqZiC5m5sNE\ndAmAI0X90wCuELtfXtRt4QYAO5NP26iKxr1tIw91WC6mzZWzuLiIxcXF0e979uxJ7ssXjUMAPgFg\nPzN/VGzaC+DWonwrgPtF/duI6CwiugrA1QAeTj67uqnSjumh1VOHpx9DrNevRutsLIxehtElfNJz\nPYB3APg2ET1S1N0O4HcA3EdE7wZwEMBbAICZ9xPRfQD2YzC2eA8zW6brKkm8eWwuVG/N5EY7huu4\ntdtCmvi7Rr5NTnTaJGtvmSoVzPxVuEf/Nzr2uQvAXTOeV3upamTekYGiPdykIkyEjYrpkalgGB3E\nvHcjEyb2ualqUJr4SXHi+eSOVY+xYIwA7JuAEYmJfRPM8bvexgVNs3r9GzYZa8wB7fvPNNJJnazt\n2J+BS6jr/qZwRluK3OSIO7cdM1cx/EYr/8t92St9J7095KpaeeUFbT43Qc6IF9cNp2s3IidNPlik\n06JtflUu5uQ/qQLa/M4kntvGQlpOutxhjLZato/YLHDTtFnS2o+9e0YVtHkw2+ZzM6ZiclUlvgFs\nle9+sn+fNupuKt+N6xia/RPr9asTu67J2g3S63OSU2hNtHuHiX3b6eEn1JpVsbNSh1duom0E0kMp\nMYC8I/i6z6HzmEAbDWBin4Mm30WPXqY+qKSOtMaGIPcNILU/uxHNLSb2IeTQKd87PSefRJOiHrVQ\nKsbr32jww5n7IJa5v8DWMCcSMyNdeBcqPMc2j7pbnWbBVtYaHaILMjcT2+b+CgUVaU/uBU8xkTKd\nw2WD5LZHUgfEKfvZ4Hsu6JMU9p7Nmu98c7P6tUo6vbrVh90l2oT9N+YmR6qGBuPzuxbe6P7WEW7/\nOOvPdOu9SL5xzPUNxxjSebH35dHJTlf+/yv6ZHPGwLv2iZk8rZuNddeiqsC6oINE1qf2Nyt2k+gU\nzf/3tIkevRt1j+DnxpM3jI7SI3lrISHvfkOWThsWXcUSlzoh3P7JTldi2btynkYQ/RL73FfbZO6b\nQEIWVbVB2KPy09T9rcSZDyfjedRt0ZiQ944WyJHRN+ba0mnLA0RaHQjjO7mQO5HdrWIxsXdR9zsz\na4ROQB+bC/WO4OsQ9RzfCGqZv+iVNrX6TtNbWi/2tUfbpJDDe28JORdKpe7TqJ8u0J49qz6KEKjf\nTqkqEibHoNtoJa0X+xiibgxtFt8MeXS4outLHQV3LX6/UaoUVAuX7C1zJfaVUffEbg34RuJtFvXY\nOPzGJnnrHgXbqNuYgol9bnK/ozOO8kOeO1u3J986P91zHtmfTpVqaVcl5r5+7dvAXGBi3yQtGOGH\n0OYJ2taRKoypQh5y47ARv4E5FfvofF8xmhTTd853N0POnRDxrVPYszwTtqIc9oAjj31qbH3v4+Et\nQqdpvN/xieiTRLRMRI+Jug8R0SEieqR4/ZLYdjsRPUFEjxPRm6o68bljm3i1mE1sG73q2K8zbCBN\nmFP387EpXhrr4uWiqnNrzwF7Rch/3qcA/B6APxZ1DOBuZr5bNiSi3QDeCmA3gMsAPEhE1zDzmUzn\n2y6a0q2AwaVv5WxVFktuMc+R6mAuJ2jboInc9AkYMXj/M5n5K0R0pbJJm526BcC9zLwO4CARHQBw\nHYCvzXKSQE26WuX/f0P2T2oO+xDRzhmx05bQTC2tcVSmy6CDeLZXOYFb+0jdaAuzyMp7iegfAfgm\ngN9g5uMALkVZ2A9hMMJvjO113CVidSrnOWXQSJ/QtmWCVuuvdcnPqvTm67C9LfJmbkmVnd8H8NtF\n+d8A+AiAdzvaql/2HsJ4EdTLAVybeCKV0RZruS3n4aGOCdq6cYZcpjAvg1ybZ62VpaUlLC0tZekr\nSUqY+ciwTEQfB/CF4tenAVwhml5e1G3hBgA7Uw4+C20b5ee2dnwx9xmicVK31xGTH1OfbBu5bgBV\nWS85YuBbnQIh5uDzcscMZ3FxEYuLi6Pf9+zZk9yXf8WNAhFdIn79FQDDSJ29AN5GRGcR0VUArgbw\ncPLZtYU6ImUWxCsDm9umT9JuYmH0Stmem7ZE62xsLIxeQ85sbhu94jpDTwNMennRrcf710tE9wJ4\nI4CXEtFTAO4EsEhE12Jg0TwJ4NcAgJn3E9F9APZj8Gm/h5mzz9l3IjkaUN3NoSXWTp1J01pJk4/7\nm/UbQQym23NBSDTO25XqT05pfxeAu2Y5qVpI1Zs69stoBYWkNc6RJyclwiY2RYJWX+kDx1M8+9SI\nmJD9ck6e2kRs75iTIVZDtE34PeQIp2xDZE6Om0QydefDafIYlZ2zfVVoAhN7YG7eharSGseQ+2ZQ\n9+IvNUVCWIdp1GnHVNGf0Rk6L3NZLyDHSL2qxVORx/Blu5Qi6rM8UuPT8++nxc7PvmpWW0hVwpn1\ncvputeTDyW3/tO5mYLGeuei82Es6fTEVnXyd0TSxtGXVrAvnytmpO2Vq4yOnBtaip627i/SOTutj\npeTQodQ+Kvp2EJb1sl0PNZl11ey0+iTaliLBRZORQkYrmUux3557wLjgKKf2EUOMmEe0zTER67d/\nqk9Mlur1R0XghETlpEbY5LRsqrJrkm8cZsG0ibkU+xJtdApSzilyn5gEaLOumk0lR/bKGJzROqlp\nEdoQYRMrxJ3Q306cZOeYf7GvklnfvdT9I7SpLl88ddK1zUSvmAXyjObbYm8nWzYm1m2ke/+Bs5Aa\nNVNHvzVE9MRYMyG+eGqq4hgfftaFVK76OBtHxNa7hLiOHDcx7aXemmVjoG9i76LNaQ0q8uzrGGnX\n4cNnuUlUFVuf+/mwcz/pajeMKjGxl1T5btRs+YSkSRhS1WStb3tbUhk7yZniOJVOx8gbbcLE3kVV\nlo/sL0OEjm/VrBzN5hbtmAeESOKsmdkWUrnqXQupvLH1PhsnxObJ6dnHDoaHN4/UG0P2tIY2mq8L\nE/tY5DuWc+CXIbzTt2pWUseEauo3hpjcN1luEr4RfB0j5thj1Gn/ZNdjE/gmMLHPQZXf9pPt5PCT\nyrGQqs7wzdry6Mx6E8jxYJHUtjl8+Kw3ORP4pmlU7OfyTuOyf2a1hRL3lwIYI9ohgpqeUyfNmknd\nT40acgi5Gm4pRc+lWTltnJhIGinqOXLxZNVkE/g20ajeyr8x1wNJtinb5Ul35kEmGhne/ZRQcCD/\nKDglnr+WlMSCrHnrc1OljaNRqQ7bTHEbaeXIPvWkthc7bnONgl2+eAv+13NMCMtVszGCWYd/nzMM\nM/YmkXz+vtz1OSZaY+yW1lkztXRsZGIunZRkYi2YbUpdVXl5ZL+JNk5MpEzIJOmsuXFckUJVPX1K\nRuA4Y+vVRVWljh37FT99Ns9kWSPVmnHZP7OO4lsTh2/Mgol9ldTw7SHVxsk58h+0qS9bZu7sls5w\ny1lH7nVEzOQeUOeYAzBaiYl9E2jfCHL49wtxo/hx2+Yjd6pcYKUeL9Wnz532IFVQc3ruJuq9wMS+\njXisopjHD7psFd/o2GfXTOs79HixEThapJCrnByBo7XJkfYgxlZxtU39pqCN1jtpzax7yhb9Mw0T\n+44jF1LV8TDwHJO5qQnUUkke2ddp4+TQqdYJuGuiwleuLBa01/RL7H0TrSHROr4J2pCyNunqO0Yk\nWvbKmHQJvuyWIX3HjPxD9tMmj13noH6TCJmU1SJwQnzsdc9230jbhavtrGkPJL6JXWeKBN9I22gT\n8y/283KF4jo2S+W0C4xdbKXtl7J90CYtp04qWh6coOfL1jmyz7FfSH++m5Ixt8yLFJbYPpdX5Sd2\nFK9R52RtTA4c1/GiwjTbZt3I9rlj6G1wbUzQeVksraytIy+9Zr3kiM+POEZq3vooy8MzSTqJLz4/\np6UTMykL6BOz3klZIN2a0bbnmGjV+jBRNwLpvNhnpcq0xhXhCrdMHdnXufI2+4rXUh8ZR/ZVpSd2\n7VelgGtev29ewIl5QV2iI5JWJns+nE6+C2HErqDV96t+UVWMpRN7k/A+icqVFqGqB4dUtbo1dRI4\nC9pF2deONuH97yeiTwL4+wCOMPOri7qLAPwpgJ8CcBDAW5j5eLHtdgDvwuDP7X3M/EA1p55AVSP3\nGqJ4ZGy9K9wyZtLVF38fG4Hjj4HX4uX9frtvP+fxfJOypQgc6GXfdk3LXG19N47YkXbnBtUWrdM0\nIU+7+BSAmybqbgOwj5mvAfDl4ncQ0W4AbwWwu9jnY0QU/kSNKWwTr3p2DOgvZ78VsYkF9bUhXnF9\nbBu9NEL6HfZV234bC1stnI2F8cvFBqaL6qZ4+dpq/W6I/XvLungZVeKVKmb+ChFdOVF9M4A3FuV7\nACxhIPi3ALiXmdcBHCSiAwCuA/C1TOdbHx3071OJtWyGzPrs2iofQuKdlHURs7o1pq8cI3FfCGVI\nIjRtYtfVNvsjCH2Y4FdJqoztYublorwMYFdRvhRlYT8E4LLEYwBI9OdjtUKzVVL2n1bW7J2QhVRF\n2RVb74uwCYnA8SUbS7VmXMfWIndy7FeqFz79aFTvWjzls0pCBDXm4SWp4tsK6lhI1bqLngtmHrMy\nMxPRtDGAuu0hADuK8isAvHbWEwGwPUbkZx881jriL4n9Qo6TdxzHE0Lp2y/39qgwzdwPEQ/d7mqb\n257xraBtjdffuQmF1rK0tISlpaUsfaXK1TIRXczMh4noEgBHivqnAVwh2l1e1G3hBgDnJx68dnJb\nOto3iZCJXQXfCDxkhO4jZrFWTLI1dx/h+yU/fSrH4wNj2vr6SA5/TKSW4+X2rvrH4uIiFhcXR7/v\n2bMnua9U6doL4FYAHy5+3i/qP0NEd2Ng31wN4OHkswvAewEuEZ21reskQsoR8MiBiIvAiRmhaxZK\nDoulfIyt+/kie2L3kz69FPuRV++KwIl5rmzMw0tSI3Nye+9ZSb3D5bgzGrMQEnp5LwaTsS8loqcA\n/BaA3wFwHxG9G0XoJQAw834iug/Afgw+0fcwc/Q0T2eeKzurx5/5xuAjTFxTY/Gnj7p98wIxjxqs\nNJ6+qtG8T5Sr1L9ajleHz243iVkIicZ5u2PTjY72dwG4a5aT8v3ZeE86Viyr8vo7HNFTVYROSH9R\n/n2IT+971KCLmJQEPjumzWSPwIkR/q68Sd2nlRJUy8jed+UuP10j1v4ZtgkY2Q8diJgIHNkm1o7x\nWyypcwTh5xazn9e6AcZiHxNVI8suAY9ZgBWTqjhmojV2UjbJk0+Nuokd7ZvwV0krxT4VbyK0WbdP\ntqkx+ickB44mkiF2TWo+nJz57LNE4zQ5Kdu1aMEsk7JVCr+Rm86IfaWj/ZzvQo6bgTKyrwtfvHzo\n/tPbpOXRqeyBJC5S/O0qFzmlLqqqBRuVt53WiH3qiXhz17u2x+zn6yPkGFrZoYtaHpyYCBxZjp2U\n1XPcjLev4Sy1XrNYVkcrKfRju7d7rCLh08vRvDf3jc+umSxrFktM27bYMTOnRo6dlfa18R08x0N/\njUlaI/ZZP7ocK2FT28QcO/MErh56Ge7TxxITxRNj+bj2G9WFPDi83MkAX4jlZNnXttGRtHLc1Mnh\nytIimCC3idaIfS0nUtVNILbfiPbDidmYB4/Mgj55Oj2EUhJzEwmLyVe+dWipECZJfZZsqY+ItkOq\ntGNSbe/kCKEcI/BZj2HkojVi70I7QZd/r1o6qaP1WGvGt10reyJwJOm5amKjccJvLqslS2e6bbQm\nLJsUywcY2zdO6+b0+HzU0XqIHXM6Yr+cdkysRZPq9cccQ+3MJc45MshpfdjNIBetEfvav/BpV54a\nMRMTbhlAah4cf9KwtBDKmGOH7KdFCkmSHy9Y7mRMTN75mL7qWBQac3MJoRbtNIFuI60R+9aciIZL\nb2Pi733bPW9ATMhjbmJvGENi5giiQi9jrJtBJ8MT0skh4FWlHI4h5ttBkE8fcxfUCJkkMeqi1Rob\nghpb77qqGNGOaesKt3TZOAtb63xPospvx0y/ecRE/Mg2vmgduV+I5SMjb9ZOD+qDrJtT2Frvsmt8\ni6Zk3aqnj9zROKlevyTKIvJ9Dcot4L6VaRK7ecxCK8Ve8+STTzRVtFPrM0T0uFbLqm0j7JgqV97q\n5zZ9IjkoN06qdaPVx0YNJoUpZiY19DKL/qVO0GoHN2unaVop9sl/FjkdjToEPmLxVGqoZIxPL9vE\nHs/3TUI7p6BonJiHkJQ72Upq2GTMqDt3NI6P2HnPpDBLE+p5oJVi78O7mtalTTkWR/nqY6JxBK5J\n2VQ7Zhgpk3sxlrRYZDTO2uh4+mhd2jvDyJsQ60aNvHGlKtasG2Bst4RYLNKaOaVsrzsaJ+bmkrwA\nK2bSIsZiibVj7KZSJa0U+9YlQvPVp84ROPbzjapzxNz7rJkc6RJSJ3D9eekjHikoSZ2IbTJePtW6\nSdbNmPDHWDFPwfz4XDQq9jlE3fsowtxpDyQLEz8ny56JW9ekrCaCpZFvgB0zbC/j19eUuHhXfyHp\nEsqj/K2jddc3Ce1bgDYRCzgmY0NG89poPSR7pexvVdmeO3vl8Nxi7R/fNwmJ17qJEe1Un76OG4Mx\njUbFPuRjjrohpKYq0EgNt4y4YbgeIl5q4zlgahimT5Rj7R/fdu/DTWJSIPgmYl1tQhZVae3reHJU\nrAXjS9KWFHXj6jC3UKeGbOY4dn/pzMi+MtEH0uPlM1pBrhTG4zq/qMeEUOZGs3+ibiiuFAjaZGzI\n6Dl1UVVVz3x1HTsms6bvPLNaN7N0aNZLG2mlZ+8j6qRDbJwhLgvG11/spGxRLo3sHcI4tEeCYtI9\nIZSa7TKo32rN+LZPthlO1sq2vnNeXRO2khB1aeOUxF6zPLT0BpPlVU9bV39DS2dVqXPVh7T12T8h\n56mFk8obQGXWTeoErWWybJrW2Dip/v22GDGXhAi7tt03QRth/0ifXqKvMJ0+2TlZ70tJkHtRVUwc\nvXdVrCu00udTxwhj6sRn3SmHs4/mY6JmYmauU7FHGNZFa2wc+ZHPfFKpHcTeOBJ9f9aEP0CUtbY+\nMY/19GMWVfnO1zkJPEpoJvpdd/QbYyH7hDFm4tPVh48cN4bUBatRMfSp1k0OP9389ibopI0TRUik\nTEw0jtafr9+J8tC+WduhWzBaNklXGgLNSpH7udq6rJmYxxlqcfau8ynZO8XIfU2kOjizKqyb0+LN\nkpbG6Sl1gNuO0eLsXWkPYo7naxtyc9Hsn5iY+6DBru9rTkw8aczKtNSEP1liSI0JWmnjyI/ap8OV\nEWLNaG0Dbhh6CmP/pKuvrS+E0tVvzKKq1JDN0nkWI/qg0XyqHaO1ibWNZ42dj1kDIMmS48a1Y0oa\n4VjBTbVbTNirpJM2jsxbv10T2siRdhSp3wg8Au8WUd9Ie7qYhxwjxxOuvMdTVsWWwipdq2LLF7WV\nEDvG5/W7NMaXLVMjZuAbcozkQW4OAU9dbRaD+fB1Mf82TgjaTSLmhuFqK+qlT7+6Y2BfuKNVttox\nq8rDPwb1ulUyFOvYRVVaVE1qFE8pqkhE2KyeKvaTGStPi4lYVxSLZpW47BgtukWLtJnsQ1uktR7Q\ndnVKnet8ZNl1jCzRNqkLAlIT9/j2sxF8E7TGxtmZs+OQ0bwm7Dm+BTiE3/f0qZgJWqc9kjGqJkdq\nZFeOm+S0Bz6d8pVjJ0993whiHibiG+WHfOvwCnyI/5MaWuk7no3Q205rbBw5ILogphPXaLwqtBtG\n5E1iuIAqdXFUavrh3G19NytXjpuRZRMj1LKcupAqd0bKVM9eOyefzRNEyAg8amZX6atJUY/5oIxJ\nGhV718G1mHtZ58gs4D9ITDSO71tAyH6ivLpjHFOvZaRcK9kjWy2b1VLWSN260eyWtaC2W9vEtJVt\nVtbG39Hk4qjSQqlhtE1MlIssu7b77J+YxVGyPmZxlOvctAgc2Z/L5okazefO65ASueM6Ro5zM2ah\nUbGXf/Nor92oAAAVeUlEQVQXiXLyYqtZQ3dibgyu4zl8eo2QSdfxBK3futHqNzy2y2TfG8p+URO7\nQWkPip857JjUMMWqniKVGr+fvDjKJYYxHnmO9JxG25lJ7InoIIDnMfjzXWfm64joIgB/CuCnABwE\n8BZmPq7tL316KfwviTmJ1FF+jsicCGSiM80X93n2sdaNnqtGX4Xru7nELZQS17TuEfscK2FjxDWH\njSPJ4dlr+yUvjooJ6bERc9+YVeYYwCIzPyvqbgOwj5l/l4g+WPx+24zHKeFNaxzi42uj8ZDR/A6l\nLOvOHhfXZHlBs0r8i6pWR4uVwq0b2d5l47gWP60Ut+CQtqU2hU0zirQBcOaFc0ZlNdomxI55QWkT\nEo0jyyeUtqk2jpae2HVuIRE2w/ZBAi9PxDehoLWV7auyfGR96n5GFeQY004+F+5mAG8syvcAWIJD\n7OWf4y5R1j7yIDtHu5qYEXyOkX8pdFzPUa/hj7MPGWmnRdXExfVvHc0Drth58afhezh3zAPAY1ep\n1v0wcN8xskbYpIY01jFDnYp966iCHCP7B4loE8AfMPMfAtjFzMvF9mWUdbyEtHGeF+WoaByfwLtw\njeJj2ibOEfgjbKZHyoT58OGLo7T4+6DQSy3axrU4Kocdk+rZzzpAjTleyDGyR9j42tY9ek6Nz5/W\nlzErs4r99cz8DBH9JIB9RPS43MjMTETOsYtrtO4bVDujcbTBc2o0zg5PW9lGZuQVZZn7RrNbXPaI\nVj4lbo2uyJ0VnBO836pnoZSrrYy2OXVyfLy1YfmkeINcdsvJ4ucJz/bJPk4qdZpdM9lfSqpiWR9j\n47iibryefOydaCWibR02TkwETepNwEb8szCT2DPzM8XPHxHRnwG4DsAyEV3MzIeJ6BIAR7R9H8JY\nI18B4EaxTftInScaM7L31cf4+47tITnqfZOuWjnMjgm3cVyrabVQz1JbV6ZKLXbel2s+dYLWZ9e4\n6uuwcVx9ZZl0bWrxQKrI1pVTZ35ZWlrC0tJSlr6SxZ6IzgGwwMwniOhcAG8CsAfAXgC3Avhw8fN+\nbf8bAFwmfpc2zkVIRLuaEAHX6kJuGEob37NkZTnGxonNVePbz/dgEZ83Dzj8eZdoa8IfcmPQJjal\niLqeQZtznjHVxgkSeN9BUkfP2jFykGOOIOYYs7TpPouLi1hcXBz9vmfPnuS+ZhnZ7wLwZ0Q07Oe/\nMvMDRPRNAPcR0btRhF7OcAx9gVXESDtLugSXpVPUs9guQyx9tonPuhn0sTWq5pSwa8rWjbR0tkbV\nyP189o/cb0XYNdK6wUlx4ScnfgJuu8XX1lcfYsH4onhiFjzF2DhRAi87DxF4X+KeHCN7ScwkcMw3\nAls01QTJYs/MTwK4Vql/FmVXxomM13ylKGsf//aQM42xY7Syb6Wso1yybhbSRt2usn+R03Qbx5VA\nzZcUzfXIQO8ErGu0njqy1+pdtooU/piRvUvMteO5Qiujomp8s8d1+N85+otpm3pudjPIRaMraF8m\nysdE+Wrfjr44e0nMTcJV57kJxPj0spy6gtZlx2gC7krJ4PPyvd78oJMxmi/us3RC2mp6EjJCT01o\npol5clRNzMRm7pj0VE9LI3YkXpfHb8TQqNiH4D1BTZRdQu2yY4Ztzg5oe/bW+mHKYsAd5aJF3qw4\n7JgVEUEzLMuoGmnHnHLaOOds6TfEClpZKfY7MT5eaXHUSRE7L+0WzZpxlbVFTjGROT57CJjdjpHl\nKK3z2S6TbdYj2qZOtGrnF3LjSLF8tP2ntdXahNg/Ie0NSaNiL8N0Xi3K2p+s88HiqXHvMfs5bh7D\n3Dcx1g0gJ0HDFzmFTa5u9fpLOeVdls6ZrZZNafLVlWs+JkmZNpJ2LaTy6V6Okb3LjtGOFxQXn3Oi\nNXaUH3pucj+fwIccI0aUTZCbplGxl9E4y6IsbZyZTzBkglbz7ANuLkM9DAlpTLVxtLh3lx2jlV3n\nU+pDWDPqg0VCQhY1H97li2t6E7KCdjhCzz0n6ctk6fTjfSIaIvA5wxtTbwwxN5+QPnzHNbumCVpp\n48iTGkbjbHcJcYxo++pdETiO3Dcr52oP2XYsRlLsFreNszU6RstfA7gtnVOKjVMqOxZHjSwbn10D\nAD8W5VVlu7RmZHSMtjjKZ8G49suRGjmrTePanjPu3XU8KNsn2/hEuSpbJUTgU/czQmhU7A+K8utF\nOTm+QLNjUu0fR1mmLdYeQhIyytcsFrcds6P0c9B2fEK+eQFnv65c88N4eZ9dM+h8a71r1aymiyGx\n9XVEG2b14WMEXtbnmKDNKfCTbVLahmD2Tl00KvY/LcqHRFmGYaopFbz5FEQ58zcCLfImxLqRI3Pf\nA0m0PtwPIQl/Bq0znLL0LNhRB1vrgDjP3mX/aKHlMfulpj1w9eskxYev8k40q8BPtg/dntpX7v1i\n+jAmaVTsXY8l3K6UnflwNFF2ibbDjhnVu6JxRH35iVNbc9y4bBM5Mh9aLC4LRu53AudP7ddVHu53\nYu288TFkLpvj548vSlo2wycPxNgxshwSjXNaqfNZN3K/VBsnOY1wyGj9lGd7zkVOsk2soMYseErp\nN3a/mD6MWWhU7B8TZbkKyzs+8S1+Ckl7EDGylytky8nNpueMd43sh+UNz3bZR9hjCZVoHBEvX340\noCfCJiRnvG8/l+Ohxdn7rBtAj8YJ0VavyLvE3CfgMfnl2+C9T7bxtZ22zyz7xfQR25/hojUTtNqk\nrCxvD/DTR4PnkNG8Noo/T9SdOy7Kh5CUJ0o9k6CeydOTOF/dLuuHI3TXt4ATSltgHC9/Uo7gj4s3\nQz47TBuB+7a7yiEpEIYaqaU0CNnPdWOIyizpuqOsKPUxnn2q2LvO0yfaqSP7kLa+c4jpI7a/lH6N\naTQq9teJ8vdE+WdEeafPT9eICbeUZcd2uWiqPFE6fYLWnftm0Ed5ZO/Lo+NKh6x78mtDH17z40PK\nMT69bBPzMPCQ89F8dpf37g2RjBH4kP1yTLRq22MWR0lyx73nGLmH9mVUTaNi/7Ao3+Ros+77e4oQ\nbe/KWod141o0pYVFrjri4ctlLXZ+uh3jehShM8Jm6M+fdjwtKsaOcY2ktTauY2h95AihDIqB12aE\nXQI/qzWTOtFah8BPa+/bL3T/mL5yHMMIpTU2jrRu1JMKiYH3TbQKa0a1dMT2lXPHE7HaxCegT7S6\nbBWtHNP2OC7Utz8vzqdk2RRvmMuOOe4oDydjXROxvsnaGBsn5IElXmGPsWZ8oj7ZJqfYt8GOmXac\nmD5C+wrBhL0uGhV7OSn7l6Is7Z0LhgKcmvUyJG2xMrIvpyqePum6qlg7k/utKZaOXBzltGaUY0i7\nRj7gG6fFBWijZynUvhG41EJX3ppZbZwQeyir9x4bA+9LPyzxraCt246pyhc3ge8qjYr9g6LsSnp/\nqhAG+bzaqMVRESkQSlE3C9P9dkAueNJXt7of/bd1oVQ51/zWhGWlOm3FK1AeHb+g1IWMpLWRvWs/\n7SYQY+O4UhJHee++SBrA76f7bgJV2jF1xLXH9BHbXygm8E3TGhvHFY2zTRPwGBvn3PDyiRePBVWz\nawb1521p47RYRFvZZljW6ib7OH5mUH9CWDRrR8Uj2WXKguewtV7WuWwcrT4ktl57WEjMM2Gd//8u\n0dZCIV0j+6rCImOEet4tmNRjG03QqNjL0fz/EOXrRXlk44Q8ANxn43jy3ciJWC3EEph8+tRWiyVk\ntK5aM66FWcUofk0+IUoKqkuINfH1pRyW7UNi67XQSZ/AAw5d8Am8LMeEQsr6EBtHO6fUFa1aXyHt\nbYRu5KdRsb9PlG8VZTmyH0bjqGkTAF3sXfH0jrz0G4XOrgZEvGix8yE2jrbfKSVvPVBOUjby5E+K\nk09dpery7H0Tra5RvjaKjxL4EAvGt6I1tx0Ts1hp3i0YE/Z5ojU2jhTzneKstvtsHM2acdk1L9bL\nJy4YiKvLSvmxw24Zlo/iJaO6Y3ip2la2Gds4PzGuWxm3PXl0XMZRJarmR6KsWTeyPsS6OaG0ifHp\nZRun364J+Ipje4wdEzPpGrtKdVYxNwE32kOjYv9ecfTPi7/B60X5Ik3AXaP1Ydnn6WM8mgfGo+7y\nQ7jHo2tX/YoSeqltn2wzGtmfEW3Fk6HUUbwUap91I+tDrJsXlLLLe3eN7EcinzpCr9J79+0nyW3T\n+PZL6SsEE3WjTKNi/3vi71EK/04p0BsTPwH3owZ9cfaiLBOarY1CIf2JyTRhd20vpS/QkpSJSdcz\nz4m7mRYDH5OyQJa1GPrJeu0m4fPmgYlRvE/Ac3jvvlDI3JOnKdZLrMjOKuwm6kYYrbFxZFbLnT5r\nxhVtM7RmRLCKcE+wIcrHdoztlqHFckw0PirsmCPi0ejL2DXuo2jvsnFK/a2Ny88dLeybw+LucxR6\n+fiUumnlE0qdL2Ol3M+pIS4BX/Fsj3nQR9tG6LZq1Og+jYr9neL5g3/0xLh8vRg9Xj3USIffHpOq\nePhkKUB/2pMrbLI8Qt+5pb4UKukKw/yxWN16dMew8Rjpt2sCHiLwOcQ+yo7Ryqlx77ELnkK3T7ZJ\naRuyXwwm7Eb9NCr2e4TA3yEEfLsYgY/E2vX/oU3WytG+yGS5sqA/zk/LGe9KX6CJuTNG/tlx/Zlj\nik0jJ1qPYet2YDyizy32UltLaKLtmkjVyqnRMalx71C2T1JVXHvMORhGs7TGximlML5IlIcaKUfz\nFzjaDm8Sl46rjv7kWO2fERt+qJTLdZeMykeEdSNtnKG9s7wyrjt5WNypDoskZIexteyzbmS9S9R9\n0TYlX91nwQDju0RsWGSKgMdEwUy28bX17ReCCbcxPzRr4/zquPylT4/LrxSj3CtfVxTkaP0ljnJx\nQ9gQbUNG6L4VrWVPfqsPHyXwgw631rlCK302jne0/ryjHGPHhIzWZ00aBmX7JDktFhNyo180a+MI\ngb9TivbrRHk4ipdi71pBW7Q5foGepuBoaQJ2q4Afc060jsvLa+PJ2ucOFSN6KfBPifNZFmU5cn9m\n1PEYbTQv612jeefI/fmJn5PbfSP0kJF9qsVSxwjdxNwwJK2xcWTe+u0vExsunfgJAFeI8lXj4o9/\neuDJP4krR3U/EI81P4BXjMpPiU4OFp1I6+aptfH25w5ePD6IfDL6IaVOjtafEWUp7EcnfgJlMS+N\n0IcCJ+MmY0brOZKGQdnualN3FkbDMEKoROyJ6CYAH8UgIv7jzPxhrd2dIiHOqV8elzc+Py7v/IWi\nIL15cTN44apxvPzBQuQPijuAFHXp2cv6ociXBP6AEPiD0MtDkXfZNZp1A+jJxkpi+KwoD8VcNg4Z\noWs53FMnQetO0mUYRm6yiz0RLQD4zxikq38awDeIaC8zf2+y7Z43jsv/TNTv+g3xy7XFTzGCx8vH\nxeM7RMqBwrKRdoycdHWN+A8+O6hfe0LM/P5AHO+gKGujeGndHANwagnYuVgW+JKwnxKNtQbaaD1E\n4H0LkCSpaQH+D8YfwDwK+JMo/7HNG3Z9faWKkf11AA4w80EAIKL/BuAWlB8zu4XHRHmXtHH+7uDH\nty+5ZlT1V/jro/IPhDUzrJfbv39mvN+xRy8b9/sd5eAugZdlGQ8/GoEfm6j7Ewwy8McIeN2Tman8\nAMBfq6jvNnAQ8y0WB2HX10+qEPvLUB7rHgLweq3hfv6jUfnOZ9453vA/RaPiFvGKiw6Mqp7aIS2Y\n8cj9m/g5AMB3vv3z4/2XRF/fFOVHRXnY9SkpnE+Lspxp1SyWyZH2Mxg8e6uObImGYRh+qhB7Z97D\nSXbTO0flPaL+zjvG5T//t4Oh/W/jt0Z1D3/xb48bfFbs+JXi5wFRVxqWS9E+Isou0R4Su6LzhYn+\nDcMwmoWYg7U5rEOiNwD4EDPfVPx+O4AzcpKWiPIe1DAMoycwM/lbbaUKsd8G4K8wcNt/COBhAG/X\nJmgNwzCMeshu4zDzBhH9cwB/gUHo5SdM6A3DMJol+8jeMAzDaB8v8jfJBxHdRESPE9ETRPTBOo9d\nBUR0BRE9RETfJaLvENH7ivqLiGgfEX2fiB4gogt9fbUZIlogokeI6AvF73NzfUR0IRF9loi+R0T7\niej183J9RHR78bf5GBF9hoh2dPnaiOiTRLRMRI+JOuf1FNf/RKE5b2rmrMNxXN+/L/42v0VEnyei\nF4ttUddXm9iLxVY3AdgN4O1E9Mq6jl8R6wA+wMyvAvAGAL9eXNNtAPYx8zUAvlz83mXeD2A/xpFW\n83R9/xHAF5n5lQB+BsDjmIPrI6IrAfxTAK9l5ldjYKm+Dd2+tk9hoB8S9XqIaDeAt2KgNTcB+BgR\n1Tq4TUC7vgcAvIqZfxbA9wHcDqRdX50XP1psxczrAIaLrToLMx9m5keL8kkMVgVcBuBmAPcUze4B\n8OZmznB2iOhyAL8M4OMAhlEAc3F9xSjpbzHzJ4HBfBMzP4f5uL7nMRiMnFMETZyDQcBEZ6+Nmb+C\niWWNcF/PLQDuZeb1YoHnAQw0qLVo18fM+5j5TPHr1wFcXpSjr69OsdcWW13maNs5ipHUazD4QHYx\n8zCofxkQSfC7x38A8JsAzoi6ebm+qwD8iIg+RUR/SUR/SETnYg6uj5mfBfARAP8PA5E/zsz7MAfX\nNoHrei7F1rSFXdebdwH4YlGOvr46xX5uZ4KJ6DwAnwPwfmYuZcHhwQx4J6+diP4BgCPM/AjGo/oS\nXb4+DKLRXgvgY8z8WgxWw5Vsja5eHxG9AsC/AHAlBsJwHhG9Q7bp6rW5CLiezl4rEf1rAGvM/Jkp\nzaZeX51i/zTKyYmvQPnO1EmIaDsGQv9pZr6/qF4moouL7Zegu8tp/yaAm4noSQD3Avg7RPRpzM/1\nHQJwiJm/Ufz+WQzE//AcXN/rAPwvZj7GzBsAPg/gb2A+rk3i+luc1JvLUc6B0hmI6J0YWKn/UFRH\nX1+dYv9NAFcT0ZVEdBYGkwt7azx+doiIAHwCwH5m/qjYtBfArUX5VgD3T+7bBZj5Dma+gpmvwmBy\n778z869ifq7vMICniGiYLe9GAN8F8AV0//oeB/AGItpZ/J3eiMEk+zxcm8T1t7gXwNuI6CwiugrA\n1Rgs8OwURbr43wRwCzOfFpvir4+Za3sB+CUMVtceAHB7nceu6Hp+AQMv+1EAjxSvmzDIvv8gBrPn\nDwC4sOlzzXCtbwSwtyjPzfUB+FkA3wDwLQxGvy+el+sD8K8wuHk9hsHk5fYuXxsG3y5/CGANg/m/\nfzztegDcUWjN4wD+XtPnn3B97wLwBID/K/TlY6nXZ4uqDMMwekDb404NwzCMDJjYG4Zh9AATe8Mw\njB5gYm8YhtEDTOwNwzB6gIm9YRhGDzCxNwzD6AEm9oZhGD3g/wMFhep7EHtAggAAAABJRU5ErkJg\ngg==\n", "text": [ "" ] } ], "prompt_number": 33 }, { "cell_type": "code", "collapsed": false, "input": [ "diffusion_solve_implicit(U, dx, dt, nx, len(t), D)\n", "pcolormesh(U, rasterized=True, vmin=-0, vmax=1)" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 34, "text": [ "" ] }, { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAEACAYAAABS29YJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnX+sZVd137+LN+PxL7BxKP6tjBtsFaM0hrQOLWl5NK7r\ntJVN/uGHSuUWWiGRBorUFJtKhWklF1Lh0qYiSoNBDi0uFiSWURPFhvpJoBYMxCaEwQFHTMMY/Iwx\nY2Y843nzZlb/uOe+u+59e9+11z77nHt+rI909fbbZ599zrlv5nv2+e611yFmhuM4jjNsXrDqE3Ac\nx3Gax8XecRxnBLjYO47jjAAXe8dxnBHgYu84jjMCXOwdx3FGwFKxJ6KziejLRPQoER0kov9Q1V9E\nRA8S0beJ6AEiulDsczsRfYeIHiOiG5u+AMdxHEeHtDh7IjqXmY8T0R4AXwTwrwDcDOBpZv4NInoP\ngBcz821EdC2ATwL46wAuB/A5ANcw85lGr8JxHMdZimrjMPPxqngWgDUAP8ZE7O+u6u8G8PqqfAuA\ne5j5FDMfAvA4gOtLnrDjOI5jRxV7InoBET0KYBPAQ8z8TQAXM/Nm1WQTwMVV+TIAh8XuhzEZ4TuO\n4zgrZI/WoLJgriOiCwD8ERG9bmE7E9EyL8jzMTiO46wYVeynMPOzRPS/APw8gE0iuoSZnySiSwE8\nVTV7AsCVYrcrqro5lJuD4ziOE4GZKWc/LRrnJdNIGyI6B8DfBfAIgPsB3Fo1uxXAfVX5fgBvIqKz\niOgqAFcDeDhywoP9vO9971v5Ofi1+fX59Q3vUwdtZH8pgLuJ6AWY3Bg+wcyfJ6JHANxLRG8DcAjA\nGyoBP0hE9wI4CGAbwDu47hk6juM4tVkq9sz8DQCvCtQ/A+CGyD53ALijyNk5juM4RfAVtA2wvr6+\n6lNojCFfG+DX13eGfn11UBdVNXJQInd3HMdxjBARuIkJWsdxHGcYuNg7juOMABd7x3GcEeBi7ziO\nMwJc7B3HcUaAi73jOM4IcLF3HMcZAS72juM4I8DF3nEcZwS42DuO44wAF3vHcZwR4GLvOI4zAlzs\nHcdxRoCLveM4zghwsXccxxkBLvaO4zgjwMXecRxnBLjYO47jjAAXe8dxnBHgYu84jjMCXOwdx3FG\ngIu94zjOCHCxdxzHGQEu9o7jOCPAxd5xHGcEuNg7juOMgD3LNhLRlQB+F8BLATCA/8bM/4WI3g/g\nnwH4YdX0vcz8h9U+twN4K4DTAN7JzA+E+j5AhHPE7y8UZVn/oiV1sf3k9nPOFvudJzbsE+XzFn4C\nwNmB7bH9LG0B4Hxl+z5DW+2cxX4s2m6JtsfPnX27J3EWAOAEzp1tF+WtanusXtadFG1PzNXvq+rO\n2VU36WNWvxWoD/WV0t98X3nXtHVGtD0mzuPEpP7Mc7M6PE+z8tFZESdF+di0baAupf5YZPsJUX4u\nsJ92PrFjaP0CE6UAAJwSlfKAxyP1085/EjlgrP5ooE6WAeb3YewsFXtM/lrvZuZHieh8AF8jogcx\n+XPeycx3ysZEdC2ANwK4FsDlAD5HRNcw85kGzt1xHMdJZKnYM/OTAJ6syseI6FuYiDgAUGCXWwDc\nw8ynABwioscBXA/gS6H+5b33haEGfUe7lToqp7G26lNYLdvGeseJkOzZE9F+AK/ETLh/jYi+TkR3\nEdGFVd1lAA6L3Q5jdnNwHCfGKfFxnAZIEvvKwvk0gHcx8zEAvwXgKgDXAfgBgA8t2Z2XbHOc8bEt\nPo7TEqrRQER7AXwGwH9n5vsAgJmfEts/CuCz1a9PALhS7H5FVbeLh0R5PyYzwM4SWrKETq/Iexqk\nXeNi7tRkY2MDGxsbRfrSonEIwF0ADjLzh0X9pcz8g+rXXwHwjap8P4BPEtGdmNg3VwN4ONT362qe\nuON0iTOnfYLGKc/6+jrW19d3fj9w4EB2X9q/0NcAeAuAPyGiR6q69wJ4MxFdh4lF810AbwcAZj5I\nRPcCOIjJuOYdzOw2jjNHn0fx29v9PXdn3GjROF9E2Nf/wyX73AHgjprn5TTM9pqvp9PY7vFNyXEW\n8WfPErgmBFmV/2+lz08ajpNKP/43jg3/qziOU5jxyspYrryl6+yD5eEjeGfMjEXynJHhwu4487jY\nO73GRd1x0nCxXyWlv/0WdG+Vdk1fJnwdp4t4/J3jOM4I8KGSlZF8Y6u0R3wE7zjl8ZG94zjOCPAh\nVB8ZyF+tD+GajjMUBiIbPaDD33Qblk2fhT32/ZwO5cnZDr3TB5OXdIYIZcaMZcsM5bqP5b+PHc8Z\nLR2WIMdpn8GHcnra5dHiYu+MksGLuuMs0Bux35u7X9//Txc8/zZSrq9SRGNWkUf3KPhofxT4/4IY\nfb9JtICPjgeOvw93ULjYO72jz5O9jrMqXOxL05VvtKT9s9IFVi7sncKjfHpLV6TJGflfwkXdcZrF\nV9A6juOMgJGPJwdA4b9gGyPsQY7ix/4icl71CTgaLvZOY/RF1C0hm6e3Df9lYiGNoXr3wp2GcbGv\nQ+631wENPL1nfH/6XkbxeAy8U4jx/Y8P0YdvoeM6taqFS9G8NV3/wupiyamzMjp3QqOmDzLnOABG\nIOCO0yAu9kOn+ASu/5PpHb4S1oGLfTfpwF+lDX87duPwEXwL+ITw6OiArJQjN1laY7T97Q7qr+kM\nF/fyV4HLg9MKbv84zmrx/4EpDPxbKmnZdDG8MWQLWaJ4Tp8Jt90+FXpTVeQkrPWromvn4xRjaboE\nIrqSiB4iom8S0Z8S0Tur+ouI6EEi+jYRPUBEF4p9biei7xDRY0R0Y9MX0GvWxKeN/QbCaezZ+aTU\n945t8XGcQmi5cU4BeDczvwLAqwH8KhG9HMBtAB5k5msAfL76HUR0LYA3ArgWwE0APkJEnn+nI5ze\nM/sU6Q9rOx9nwPjNZxAsFWJmfpKZH63KxwB8C8DlAG4GcHfV7G4Ar6/KtwC4h5lPMfMhAI8DuL6B\n886jqRHxHvEZIG2I+jbWdj4p9aPktPgsq3OcAMmjbiLaD+CVAL4M4GJm3qw2bQK4uCpfBuCw2O0w\nJjcHJ5fcm0iHbz7+NOCEOQVfFNAcSXJAROcD+AyAdzHzUSLa2cbMTETLct4Ftz0kyvsBvDTlRJqm\nSXHsqPAuUlKE+yLolnj/033JbtmG5eK63DgbGxvY2Ngo0pcqQUS0FxOh/wQz31dVbxLRJcz8JBFd\nCuCpqv4JAFeK3a+o6nbxuvxzdhyVkO3T6MRt6CbQRTF0371XrK+vY319fef3AwcOZPelReMQgLsA\nHGTmD4tN9wO4tSrfCuA+Uf8mIjqLiK4CcDWAh7PPztHJnC84vba28ylBG9ZM6Bi9nCT2CU9nBWgS\n8RoAbwHwJ0T0SFV3O4APALiXiN4G4BCANwAAMx8konsBHMTkn/I7mHlYrzVoanDYE63yl5s0xCqF\nv4tPIE5xlkoXM38R8dH/DZF97gBwR83zcgrBHdbNUYq6FR/9O4XoybRhj2hKv1r+S5X2t13YC+M3\nAceIi/0qaOpbz+y3pLCvcvWqKQVCpG1oYnc7EoFzJrQ6bUh2jMfuD4pOir2WvbKTJz1CVpUGuZdP\nCT4Sd1bMIHVzb8pVdfnKu3xuDdH7fDa5+OSogt8lSzHS/2EJdPmGkXnc7bW8NEWlR/Ce+kChtL6V\ntGOyz83vaqvGxb5JRq5pLuqO0x1c7FdJB9Mz9M0Pj91QLF6/aWJ3lekSSo743R0ZHS72Xaehv1Cu\nqFsFs+QxOv2ksE2769qIZnHRdhJxse87HdY/jb49RRSjywLd5XNzauFi7wQpLcSjFfZUSousi7az\ngIs9oH8Lmk6ViNypuz1CqbdSJR+vhZuE1UoKLpSyeP3bkS+xDf9+8EEsg7/AzuBi79Tw70uuvB3I\nyL/rI+ocbXU9HgSDF/s9g79CQU/0cpQLqErfBFyAHSMj/F/XMCUsnYZE+3TmnS9XnNsQddNbpqyh\nl2dauHuWjNjJ7auxqCG/I3UJF3srXRg9N/nCpcwLXFWenLbZPtXCP4CuW0FOL1n9/56aaEnTRksP\n/rKxG0Sn4+ljjFGgPStmr+iBJLRI299GyvEa0j1NUEsvnlrlBOzKnghiN4DcG0Nf9nM6iYt9k7Tx\n7SrHGFPopeWJwOT1rzJFgoYLspPIuMS+9NV2WAPaoI0UCZ2m5E2gbdH2m8ToGJfYW+iapZOiKw0l\nP+tyHp3OEQtAyRXXXF98pYEw2sG17Slflt+trLjYt0UH9O30WgdOogZ1wymt9k/sdYSjxjW2t7jY\nl2BA32IXYuq7EGIZI/je2Tq0PeLX8ND4wdLd/1UVltBKUxhmlwdtBfLksHJ9uYJaMg6/LVFvxTYK\niXYXo2c8XHK0dF7sO8EIJ3a7sriqbhKzGCudT8gdPbuF4tTAxb40HZvYzX3vbApjicaJhl6GXlgy\nRPxpYBC42K+SBiNsLPQlGqfupGuJJwKVXGHMHbWnPCVofbfyxOCTAavGxR6w2SqWb6zkt1sgwVqT\no+chjvKjeexz8FWszooZpNj3Mq1xR85ZH+WXi9Zp621YRY8TtXRy+2toP+0Jo8QTQXH8ztYkqqFL\nRB8jok0i+oaoez8RHSaiR6rPL4tttxPRd4joMSK6sakTb4Q18UlhD9JFeg9s7ZeRcJ6n9yxPlXAa\nazuf8PY9Ox8LuftZ+itxjNj1nz6ztvOpzbb4WDglPl2GxcfpPCn/Wz4O4DcB/K6oYwB3MvOdsiER\nXQvgjQCuBXA5gM8R0TXMfKbQ+Y6DjozyNUpG7HQ502XxtMY+4emsAFVWmPkLRLQ/sCkUinALgHuY\n+RSAQ0T0OIDrAXypzkkC+fq3N3dHy35WLSjp+2vx9AmeVm4GTI02JnZN75I1JD9Lou3Y+jayDBS9\nEbkt0yXqxOX9GhF9nYjuIqILq7rLABwWbQ5jMsJ3OsA21nY+MZZZO8vI2U+zklbN6e213WGX22uz\nj4VcS6dr9MVicnaRO+79LQD/rir/ewAfAvC2SNugo/eQKO8H8NLME+kcq4zs6UCYZun9Shyj0wul\nNEqMxAcjzH2/U9rZ2NjAxsZGkb6y5IGZn5qWieijAD5b/foEgCtF0yuqul28LufAdWnD0sndr/BN\nQkvhUiK2PjfCpq7FUkK8LcnPsvPhNJm8saTuFY/Dt+wwPgG3sL6+jvX19Z3fDxw4kN1Xlo1DRJeK\nX38FwDRS534AbyKis4joKgBXA3g4++wi7BWf1ikZVdMRSkfQTEmxjXL6k8TOPWQRFbeNcq0Zbb/T\n4pPLUGwjpxjq/24iugfAawG8hIi+B+B9ANaJ6DpMLJrvAng7ADDzQSK6F8BBTP6ZvYOZuxmY1bG0\nBnMUHOWXSGucItY5gl7ijVMliNo/OSmO+yKuHhE0OlKicd4cqP7YkvZ3ALijzkm1zkDsnRKZLutG\n5nQtRUKsPvspJpYPp+4iphIZK0svlGrM6+/LHXFY9N6M6MQFdOBds1aayoeTe+yuROQUTZGQQl2/\n3HXTSaQTWtkZumztGGky2+WUplIrlDiHTqRI6IoQWyybrpyzU5xBiX3Ri8nVCnkSHUiwlhYdU1+0\ntTYlslDWvXlEbwyR1AhZK2dLROC0tcAqta8iNJnW00lhUGI/GFYq/HbRth7DQt1Vs7H6oqtmU+hl\nYrIAPrHbWwYp9ntL278l+ivxpKBQ8vWo4/XsW74JlKRzdo2PyrvEIMV+jiaFP/fbm/bRYP4dLSeO\nHB2XiLDJEXBrrhrLMSxPBKGFVIBYTBW7AbSRq6akXZPSX4lj1MZvEk0wfLEfIgbhLx2zXjKffelz\naGeytub2Un20QbZl42LdRcYl9k1dbYl+W4jZLzFZq1H6ReVNLbAqupAqhdIC3pRl05UbjVOccYm9\nZE+k3IV+M/uQnr22ctYyEZu70lUbaVtH4iZrxrCQSo2tT1lIFRJJy6Rring3FWFTZNK1xGjenwia\nZLxiH6LJb6Nu303G5BcYPbe5SKu17JZdSJfgI3inEC72MZoa+Zfoz2TdWF8rWDddQllPv4TXH5zY\nzX3tYAmRLOnZtz0YLp7pykfzbeFibyV30VSItYSygsyHo62a1aJn5ts2s5BK1ltH6CZrxhBbL336\n4EKqmPhuK200myflGCXay7Y5lk1xPXaBXwUu9iUoEY4ZI7M/i6VRdyFVyvHaSIRWPItmU5O1Eovu\nte3DF7V6XOBXjYt9kzQ18Wsg51WBTW6ftCm4ojUTNbY+Rhux9W20jVFUk13gu0RnxF6+iCRkl8e2\ny/rsl4u3QeiiCkTrzEXgNLSQKiWuPSeNsDXc0mTNWJ4ItBF8ro2jbY+1LR3FE0Ieo/hkrc/+dpEu\ny2M2e0LCCuge+Sq/jcITwhbrwjKqrhtuWWK/IjeJlFTGsZDLkljEunPWTCsdO4UYpNgXISa+bX9j\naws/jVjtkTbeEtVGGGY22ii/7YRmJd5R626KAxf7eoTsmJQIm4Khl7kLqSwj9LTcOMstltxImpD1\nlLJf6PzNqYy3F37Gtsfa5No4ubH1JUQ9FEHkmS4HgYv9KmnwicESn57iz6duz43caXuCNurT5/r3\nOx0bT2RVad4tYaHOIHCx7yIFbwK5ItpGuGWJdAmWcMtYBI7KKt/tWnoi1YV9tLjYWykZVZN53NhC\nqtwIG9t+Wr4bmzWzrC7WX5GbhBD+aLhlSWvGksq4xAtS6i6kqsX0YmJfgFaOeVNa2ScnluFiD8y+\nhQIrWpMmdkOTroVdjKnAWSdc24jM0babRutdi8BpciWsZCrgbYzUoykSYuLrdJHhi31XomosKE8P\nlth6SWyErk2epsTnWzJZhkfr6TnqzSP7wMTs3KSs9kLx2MhYmxyNDWC1kXasbcnInNi5exTPYOmy\n5GXT6cVVDWJZNBWjzTQL1vQGtW8SKZ593Rw2uSP70jH0LtTOAoOSRcMgtzlSniRCXn9D524Jt5zf\nLz1CJ6XvkpaO9SahTszWjcCZnFQ6JSJwtgN1jrOELshjLfbqTdKx+O2yXPpbVI4Ri63Pn3TdfUMo\n8d5ZbdRtGcGnHE9tG5uUlT69ZnnkxtZbcuo0mcog5PVnT+x6aE+f6KXYFxV4oKffQhrWEbrWR7zN\nnqVtS1o65ptEaGI2xdJp48UhTcXLtx5HH7oof+zoEgOWOYW2r9xi4wTKnKBNFiy57SUlcu60Yemo\nk7IxSiYsKxH+uFIBL4kL/6pRJY+IPgbgHwB4ipl/tqq7CMCnAPw0gEMA3sDMR6pttwN4Kyb/NN/J\nzA+0cqIpO5YQ+FB/JUI2DcRi60MRNpYInEmbPbu2p8TnW6yZUEx+if1MWS1lnZYOIbY9JL7WFMea\nrTJIPGRzFSx/tdGEjwO4aaHuNgAPMvM1AD5f/Q4iuhbAGwFcW+3zESJKOUb32BP5rOrYDXIaa7s+\nKWxjbeej9WvZbj2P4Lltr+18ppw5vWfnE+WU+AQ7Fh8N2Zdlvxiyj1C/seOF6lPaMhp4DeEyln3x\nTl1UGWHmLxDR/oXqmwG8tirfDWADE8G/BcA9zHwKwCEiehzA9QC+VOh808jViNKiqk3srkXaBojF\n1msj8NwslpbJXtmmdO77InH2wZF9ZFJ2ro2yPbdtiN7EuvuovK/kytvFzLxZlTcBXFyVL8O8sB8G\ncHnmMQDYJmP3WnTNILRR2gyhNGS3jPahRODINtbFUVreeT1yJ3wz0K0pcQ7Cp5eTsjtevbZ4atLh\n8u1apEyKx54bbWOxf1r3+kuGEPlNpAlqSxMzMxEte9gLbntIlPej5h2haZry/bWRf8Lx1LS+xhF6\niPJvuFp+UaVfIh7E6q1r2+uubm07r00rx/PMbXXZ2NjAxsZGkb5ypWuTiC5h5ieJ6FIAT1X1TwC4\nUrS7oqrbxesyD7yI6QJyR/65dDg9Q8oo39JHTmilJQ1DbL+5eu29srEbQO5oXSM3XcIgGbdo57K+\nvo719fWd3w8cOJDdV64E3Q/gVgAfrH7eJ+o/SUR3YjJYvxrAw9lnV4ImBT7XxjGEXvKOA5EegSPr\nU+yRkPBbLB+J7dxs/r72tCKtG5NPHxvlWhZNaStaNU/emqsmdDxLH9naW/JxxmmTlNDLezCZjH0J\nEX0PwL8F8AEA9xLR21CFXgIAMx8konsBHMTkL/oOZjbP52s+ffFBclM3hA6P7DWscfiWCeGSo/yk\nt08FhV85Sdkm9wXgfdE0eX3Z0Te5PnvupIVjJSUa582RTTdE2t8B4I46J1Uby+g6hRXaP1MHwhKB\nI9vkjtBjlEjJYDnuIH16Sd2J1k7eXHyCtYv0bLy5HDURmmV7yjcTCqHU2sqy4XgpOXB0OyZl8nTS\nJsXyMVksSlTNFvZl7Sd9+njum6kXNquK2hyhcorNo+XRaVuUtRuGaQ4h5QUi2n4pdOZuNUh6I/bF\n8+FILN+C5YZRIBR02fqfUmiibbFmckfrsfPR9ou+kGRuhWzghSQp4lvX327be+/0gLrTJzcKeiP2\n2ZROkVDiGKE0C5loI/CUEboFSxSPZhU16tMHGyzfnNwm1LbNN0el0PpNwHLhJZ4IuvJF94fei33w\nRSW5nn2KraK1SYiw0baH3jFrnzBdPkLXFkdZ8uHINraFW7F+l+9nsm4A4PnqZ2yU/HykPiTgsm2o\nv9Lee64do60KNk3KWicqXIi7SO/F3kTdlbKl+jW0n07MWiZlrZRIoBY+t+WjdXPag9AoP7oqVnmX\nbBspELpix2TPF9QdgacczO2dthiX2Mdo4yZgaFvCp9dG6KVfYahH5qRbPrH9pvZNkk8/3/m0s+Xb\ngfy0xZYwTQsWr9+CaWFXykGaEm2/GZSiM2Kf4ppMiU3WmiydUH1kkrTIfqE2kbahPDgWKwXQR+ix\n/qZRMdYniWmbkwlRNZrlcxJnBfeb2jdyNL/1/Ox4QesGCI+eNetG1luieJq0Y7QRunVhVha5wu+j\n/FXTGbHvNLFBqzYgtmzXQi+VkMfFcl2sN5cQuattY2126tqwbmQ5d9LWGmGTQ65Fk7R4yrLowIJ7\n+qtgmGKfMrrW6nN1M2U/LWgkkhpBQ0txXCLNgnbs3Fj++e2R/SyvF9TiyEtaNyl9aJS+MRSN6y8R\ne+qsmt6IfdRV0a7AIvAp1o3BjlEjcyIROKHVslbvPTcaRzuetFi2RFmzjUL2Top1c3JL1FfCbrJu\nJie3e3uKgGtRPCXtGKt1Y8mpI8leQLWsbvGAoYOn7LesztqHs0hvxD6Jci6GzbMvsJ82KZtrpcz3\nkRdCaenbEjtvjsbRkpvNdyLaBLanCHiINhZHpWCJ62/Mp/eRe58Ylthr5Ap46X61CdpIaoRgW4Md\nk4s15t5ynsF+LS8hSRFqzXu3pEtog9zJ1ZTzNMXU18VH3F2ik2KvpUZQUyek6Ftu5I5ix6iWT6Tf\nuQic0EKiBO/dEo1TYuVtqI3McROPzNlTbV8edQMgvGgqJvZSDEOWjWx7MrA91l9su+zjxJI6oLlo\nHIvFFEWzR1IOYhnlW54U/IZRik6KvQX1VYQlVtPGsEzsri1vIydlQ6QsnrJN0Kavbs31+uPnpkTm\nxPLS7yQ0S8hLL9Fi0i0+fIlMvhb7J6W/ENkOSwlrxkW7i3Re7E0J0OouckoZlWsCn5IuoSprk7LA\nbBLUGpNuyUI5Pxo/a+64i9u1p4rQuU/63d3H3HEDE7FAZDI2NhH7fKQ81ZmTyvbF/rQUCHLkflJp\nmztBq1k6KcfITodgWd5rmZS1tim533jpvNhPKSr6i2hPB7HtNSdrY5OyWg53W24c3bqxYIn+iR1D\nfQqwxNG3vSrWMqIuscgpdyGVSlMjeMA2y23p16lDb8Q+G8sVlvb6lTahlbKLhCczdTvG0ras/ZNu\nMc1tD03EAmF/XptQ3X3w9LZ1V8Xmvsg7N4Qytt30xin30MdAL8V+zhGxXEHm5KkpqsYwWRtbPLU1\nZ82kL46at3f2JbcNHTtk7cTOTZ5fLLY+ZDdJ62br+Vl5Lnvl86KsxcvL+hOBNjGbR7N0ckMoY21z\n7R/N60+KGtJmnS0JeLTHipTHDstjk99o6tBLsY+iLXJK3T+1Tai9MhELzHv1U7RJV8tCqli/uZks\nY+j5bpbfXEyjeVm2jOZjbZpaKJXSNoR1QVSoTWOjecv21DZO2wxL7EMkiG+QlJG9pd+gjRPuWItW\nkZQR/nIrdpMWVe1kr4x58wmhlaHtWjnFuqm7UCpF5zT7R6LdBMxrAHImUkuEB1mOl9vWWUZnxF5O\nwLZyUrnx8qGybHu2cgwAW1Ub6dNL20QK5slANM58dMzyKJ5Y25htNG0fs25icfSh40Ujc6oIm7kY\n+pMi6ub5SGjl1PKw2jHT+pB9sti2bioDWZdrFaXE50/3SxrNh1IgWO98qdtjx7XcDNyuaYLOiH2M\nkA7LVMZ7Q4KaEgpZ94RS+lPaxiczl0fjxOpzF1WF+rPaRrZVsdVTQIp1o4VCIrA91iY3vLFEW4lm\n/8TaSkwjem00b1lUFaNEKk8fxTdJ58W+MbSbRFOWDxBcQGWJX4+3rWfdxM9HF37tfEKrYudfI5iw\nUCok9ineu+b1l1zdmmsVSVIsHZXcsEjLQUqIuo/i22L4Ym+xZgwLoubKsZuEcCZCC6i0aBVALqra\nHWkDxK2SrZ1onPACrLg1E1rElRLFM7V/hLVzRkbbiHOejuhFBE7SQqmpLsTsGM3eiQmutl9su+bZ\nW1InxGye2H5B+yZlcuGU0rbE44plP6ct+i/2mmiXQLthJIzsw2+fyo2tXx6zLtuUj7MPWzfBRVxq\njhvDaF6Wte2AbSWsRNMsicXf144R2y7JjrZpajTvAt4nein2ptj6VSIFXpY1fzsUppjg2W8Hyinz\nAtpK2OwoHi20srQdI7GshM0NDql7Y0Ckzhw7v6zjRSw5lS2efQlywp9kvVtCy+i8bO5d+LkUS1SN\nFo2zL6F89sLPJceWPv1UGLeii5VkNM5ueyTFjtGjasJ5a7TcONqxowul5qJtqi9GWhcWOya2/YRS\nH7NVtHpL21AU0GI5ZPWY7BqJ3DEmhscD9RbLxyK4sWOUuNM6dei82KuUWEBVdzGWRNxQYonOQlaJ\nRM8QmR7Tp2QmAAAUQklEQVQpY42qqZ0uIWWhlOUNUJod09REq6wvkXI4NlrPXhxVNyzSMmPsDIFa\n8kZEhwD8BJN/vqeY+XoiugjApwD8NIBDAN7AzEcs/ar57C2RMrHRvIYlBULkGLHcNzkrXS057OP7\n5d0kkiZrAwul1GibFO89Nxqn7uIoWZ8bSWPx7LMXR1nTENQdSVu+OB+1d4m6Y1kGsM7Mz4i62wA8\nyMy/QUTvqX6/reZxbDZO7GagifY+Q1tZH2mrvZAkZLsA4WibeH6a3RExsg/r4qiQjXMc5wbLJ86c\nM6s/Nqk/eUJG2ESibTQ75phSH7NVjkbqm7JxQpZNykKqkGWTNJoPWTaaXSPbyvrSy4I13HtfNSWM\ni8W8szcDeG1VvhvABgqIfRTLwiaJZcRv8foVnx5IGa3vHj3b7RjLCH338bQ8+UA42mbOm49F22jh\nhlrZ+nLupmyc3JxhElXkSyx4yn1cyfXsfUTfRUqM7D9HRKcB/DYz/w6Ai5l5s9q+CeDiOieSdYIp\nAq/tl2IVBZ4kNJ8emIm1TbRz7Zjdk71p+yUcQ3s/rBaTbo3G0VILtG3j5N4YTK8J1CyYLk98ljif\nLlzHMKgr9q9h5h8Q0V8C8CARPSY3MjMTkSk6OIU9msVSYnGU3C6jbaTVEzietG621pZbKHErRYuO\nmZVPCFsl1EdsP3nseZvmnKVt1WgbmddGCvFzohyyY+R2eZMIWTOaXROrb9vGMc17poi2ZuM0NYFh\nHc1rF5t7E3B7pw61xJ6Zf1D9/CER/T6A6wFsEtElzPwkEV0K4KnQvg+J8n4Af6WpE7X49ylpD5T9\n5uYksyddl0fHWJ4ILKP52PHm5gXESti5VwZOwylTrBJtFCwFs250DGAT7ZI2TnaMfIpwtjHrXNdD\nt47KXcwX2djYwMbGRpG+sjWUiM4FsMbMR4noPAA3AjgA4H4AtwL4YPXzvtD+r6txUntThNiyPXRD\nSLF0AmIfe3G4tuApVg6lS4hZLCGfXXu5SezYKe+HDSYvi8W6hzTE8h5Y2XdK21DfKTeXujaOpHgI\npSaebUy65nrzFiFP6Xcc9s76+jrW19d3fj9w4EB2X3VG9hcD+H0imvbzP5j5ASL6KoB7iehtqEIv\naxxDj8IxRMeoUTUhiwaYt3HO3t1mW+y3tS8cYRNe8JRu40ztFQA4IcqxY0ytmZNzUTVyv+XRNrG8\nNsHRPDATzNjiKGm9PKds1yJsUuwfS/SPFkFjsXGKRNVo1o1sX3pkLwkdIzc3DgLbl7VxSpMt9sz8\nXQDXBeqfAXBDnZMKsTdlpB3abrFjYtuVvrUQy1i5RDSOFjUTy2ev7ScjbYqEU2rJv2J54EP7ScFt\nMs5+qrMpNo4pqsay2lTrz2r/aFj8dEvb3JG43wxKMawVtHVvAplpFrQQy1h9/o0hbOnoNo4e0TMd\n0csRfFIqYi3CJjRijm3XrJnSIZSWkM0iUTWhNtY7kQXLnS91/2W05fE7Fjoj9nsj5R1nJsWaCdVp\nMfLAzL4JWDS76mXa4qocs26kPRJe8BSLiNldlhZMiqUz2293pM3uPkRETyXyJ47N6nBMXLRc8BQq\nh+yaxbYnAtst1kyoL6BMjpvQE4Y8B9MIPiUu9JShbe5Ea+j8Ukblde0Yi+UT227tz4nRGbEvQl3P\nPtRXbDtm9o1lhC7LsZWw+uKovAnaqKUTCqeMvTlK87pTUiCEkobFrBktL32TI/vpeZgmWrVFAotl\nTSRLjPItE7S5Aq61cUFeNcMS+7rEbgCR+p207AbrBpiJsiWEMn5jiFk6u9MlhKJ8gMjiqNDk6+Sk\nw+VQhI1mzcT61Tx57WYQK1t8eiDTpomJehvhjbk3hly/3QW8T3Re7HdSHKcI8dRuCaUhTqlPaMvn\nzcon9+1+yXbMugnVn1CsG1met3zCdkzI3oke47goH53td+bZ6gI1uwYIWzYp+03FNxaBE8uNkxtV\nk/PGqShadIwlP42st7SV9VarJMc2KWGrpDwR5O7npNBJsVdPKmXx0xTLaD2hbWiFrBYjv1ivvfov\ntOpVy1u/WJ7eHEKTr8B8vPx8PpvqZ4qIhsIsU6znkOVj2S9X7LMdkdwRepOhQpb4/BClffhc/Omg\nLTop9iGS3k4VSpcQ7TBQTvD6QytkpVUSW+l6MuCXx2LrQ33IY8TK6uIoEWEzH04ZsGxido3Fs4/d\nGEKWT0z4Q/vl6qnJookdUBsxl45718TQMtEa66+EgJfsN/cpwFlGZ8Q+tngqWK9F2CRE0gSjcbQ3\nUmFm3QDhHDfSNjkxZ93IfDYT2yQWVSP3O4YXLu1Xlo9WbWX56PFZXdCumRxkxpFAnaUci8YJLaCy\nWDeynDuyj060agZ+ysg+1DZX7FPOLSSIlmgcSVNtU/az9OHUoTNiLwkJfDRFgraC1mLjREbzcytk\nA8nNYiPt2Mh+Wo6lLwjns1++wnaxPLVsonaNZdWoRVBTompCoZCWp4eU8zEteModrZdMTIaE7V0T\neKdPdFLsJTsnGButhwQ8tl2O1s8P1Mv95iZiZ4umQrHxKSPtE4GJ1KPiJOQoX9ZP+9C2A8CxM2Jk\nf2RS3pLx8kfElyHfHRaadH02sl2boNVSJMg2sReW5GavzJ5ctaQk0CYXSvvtuZOguRE4dSN3SvWR\n06+zjM6LvYolVbElN46ciJ1bNLU8ln0+nDI8cg/ZP1o5aVI2lMNG+vEpo/WTynZtglbz6WNtcxdH\nJXnvdS2W0knzEWhjFfg2LBaLKDclxC7wpeiM2Mc0ea/FmrEsqgqN+CPWTSxqJrQSdl6IY7nkd9s/\nW8oEbGz7fMKyQA4bTZxj5ZQbQ6hNioCH2uYujkry4UOj7pQQSUuazdyJ1pBVJBmKwLsfv2o6I/Yx\n9mgCvhYoy+3CjlEnYEXb4+eFs0KGrJnQxGhK+VhC22mbH+PCWV3ArgGALVHGkSqHjbRrYtZNqD43\nzj4lk+WxQNuUclDYUyyW44E6zbqR9SXEHpF6SxRL3bYp+1n6sPZX8niOlU6KfTACxzJat6Yqrtqz\nHM2vhePeLW+DiodLLs+NE1qMFbVr5nLYiCRllslMywg918bRjtGYwMt6i8DL+ty2PlovezynDp0R\n+2joZSiTZa6NE3oKEG1kDH1InIGYHROOlIn77MutGTkZu5PPfktM8IoQyqgnHwqFjE2IhkbdloRm\nsu+UiVZN7FVrRnacIuAWz16b8dX6aEvg647cc0XWBb6vdEbsJXNarUyeztWft6RusSyjcS6Y/Dj6\nopmIplgsO7HsorMjwm7Ryj+Obn/x7Bhbk76ffXpWhyPiAmM2zZFAnRRwLdpG7he7SYQsndhoPhSN\nk+S3h4Q9RZybel9r7Dwt6QtifWhttf2sfeT0m4KLeRfpjNibTkSzcWKhl5EFVtPJ2JTR/PwIfLfF\norWV5ehCqTNigdXUppFphqVoxzx0ywg91EeKwGux87FFVUGRTwmL1EboTb3VqY3Reh8tmNxjO6ug\nM2Ivmcttn2PNyLqYfy/qp3H0KaIdahNLYnYimks+lGt+d055QHjyUixjoh2qj61SlbpoSWimxcZn\nC7wlBt6a+127SUi0KJ5Y26FPnuYc1+kSnRH72MtL1FQG5wfKMevmglmRRfnovokdE7ddZhbKj/BT\nO+Wn8ZJddbIsbZofVW1l33PHeGZW3nr6RbOT+/FOBzM060bu91xkuzayT/HpTdZMSNiPR7aXDIWU\n9bHJU03Mu5JOIKdfKy7mQ6UzYi8JnlTMpw+N1mNRN6J8/LzZqtjZC8BjtsvuUTkgc9xkjuxFmuH5\nqBpxzlNv3WLdADMhtgi83C9lZC/1eQfNgpFtYttLi71lxaqP0J1h0hmxj56Itio25MMnTNDKVbHH\nldzvWt75E5GbwXxMvqivJl2PzcXFK6kM5ITqj0VZ1lvi5TUv3yTwgC7gbXvvJSdP254w9UlSpzyd\nEXtp3ZwTylFzfqBusTy1Zn4qUAfguYtmo/l562VisTwdsWM2cfFO+Sm8dKc8bf+U2P50oF8AOPKT\nmU1z4skXTxvPiJVDUTWWskXggfDgeY6YgB9XtueGQmoTpk2O0F3AneHQGbGXBN9KFVscdV6gPjKy\nP7EvPEKfhVDGQizDicemXn50P5Fe+MTTM7HfGZlbRDt0AwD0lbCx0Et1IjXFjgmVS8e9547QV+Wn\nu6g73aQzYh+MrQfCmSxjk7XnLfwEcFJOxGamL5iLew+0mZto3RJpDaTAHxGrW6eTrdpoXtan3Bgs\n+eWjdkxItGMTqaGRe0pUjWXC1OK3x9pobbX9rH3UPYbjlKczYi/WhGJvyJoRASpzNs3Folw5LHzZ\nrGrz3Jnt8n1cFix/D1dWdZfu1ElrRrbdFDbOtM3Tm7MTOrMpTv5JhMtPL/wEgB+KsrY4SpZj0TZB\nMdcsGFmfMprXcsPEko0t7hPra5GSk6cpfWi4aDv9oTNiPxduGZqANXj2Ry+YRdXEwimlJz8LhZyN\n4J8ObJ/sN/PhpyKfJPA/EOWpyMfCKUOj/JhdI0fuc9oTmjD9iShrYm5NGhZKNqYJca4dE2sbo+5o\n3EXd6T/dFPuQTSNF/UWRciX2R9bSUxYAM+EPTdoC8xO0P3omMIo/LM5BG83Lcsy60Ub2UQsmJOYW\nv12WNVFfrF/VAqSUPjRczJ3h0xmxf6H05KVNM3VQrhR1f3lW5JfPyn9+0RUAgMfxslkdfman/F3s\n3yl/T3R4CFcBWLB5NmeWzplD4k7zhDiPqchbBB6YWTYxgZ/TnpBox0bosj7kocdG6CVSC4S2S5ry\nwl2oHSeFRsSeiG4C8GFMplc/yswf1PZ5kSWccmabY/Oi2QzsVKylkMfKU4GX9T/4C3FHOSS+mkPi\n2KFRvLRoNkVZG7nPjdA1MbdYMMBMBFNe0pE7QpeUFHMXcMcpTXGxJ6I1AP8VwA2YjIO/QkT3M/O3\nlu4oBf4iUZ4Ku9Dh566axctLAZ+O3A+JEXx0lL81a/Ps45dMG884JMrfE2U5sp+K/eII/uQGsG99\nftI1OFqXAi7LRwNtc/PIWEbli21Cbf8c2Pkehzga/y4gBgLDw69vrDQxsr8ewOPMfAgAiOh/ArgF\nwHKxF3YMflGU//7kxzdfNhPqL+MXdspfw8/vlB/BKwEAB7dmnT37yCWzvv5M9PuYKB9a+AnErZm5\n0fhUoJ9ZqPsUJneskGjLcmzUrYk2AttjbayCrAnxnwHC7hoehzBssTgEv75x0oTYX475sfBhQKhz\njOtE+cZZ8Usvm2z4HH5pp+4L+Ns75a+dmYn9j750+aTwVdHXY5FyaLQ+J+RSwH8UqZ+K+WIY4w8x\nubdZwhTbjhd3HGdMNCH20byHy/jxB2aR9r+Nt++UP4U3AgAe/eNXzxpviB3/ryh/s/p5SNSdkMIo\nPZjF0TgwPxK3LCRaHF0/C+D/YZ6msh46juPoEHOWNsc7JHo1gPcz803V77cDOCMnaYmo7EEdx3FG\nAjOT3mo3TYj9HkyM3V8C8H0ADwN4szpB6ziO4zRGcRuHmbeJ6F8A+CNMQi/vcqF3HMdZLcVH9o7j\nOE73eIHepBxEdBMRPUZE3yGi97R57CYgoiuJ6CEi+iYR/SkRvbOqv4iIHiSibxPRA0R0odZXlyGi\nNSJ6hIg+W/0+mOsjoguJ6NNE9C0iOkhEvzCU6yOi26t/m98gok8S0b4+XxsRfYyINonoG6Iuej3V\n9X+n0pwbw712h8j1/cfq3+bXiej3iOgCsc10fa2JvVhsdROAawG8mYhevnyvznMKwLuZ+RUAXg3g\nV6trug3Ag8x8DYDPV7/3mXcBOIhZpNWQru8/A/gDZn45gL+KSYBu76+PiPYD+OcAXsXMP4uJpfom\n9PvaPo6JfkiC10NE1wJ4IyZacxOAjxBRq4PbDELX9wCAVzDzzwH4NoDbgbzra/PidxZbMfMpANPF\nVr2FmZ9k5ker8jFMgusvB3AzgLurZncDeP1qzrA+RHQFJkvbPgpgGgUwiOurRkl/i5k/Bkzmm5j5\nWQzj+n6CyWDk3Cpo4lxMAiZ6e23M/AXMv5QTiF/PLQDuYeZT1QLPxzHRoM4Suj5mfpCZz1S/fhnA\nFVXZfH1tin1osdXlLR6/UaqR1Csx+YNczMzTLDmbmM+63zf+E4BfB3BG1A3l+q4C8EMi+jgR/TER\n/Q4RnYcBXB8zPwPgQwD+AhORP8LMD2IA17ZA7Houw3wmqyHozVsB/EFVNl9fm2I/2JlgIjofwGcA\nvIuZ5cos8GQGvJfXTkT/EMBTzPwIZqP6Ofp8fZhEo70KwEeY+VWYvApmztbo6/UR0c8A+JcA9mMi\nDOcT0Vtkm75eW4yE6+nttRLRvwGwxcyfXNJs6fW1KfZPYD5R8ZWYvzP1EiLai4nQf4KZ76uqN4no\nkmr7pQCeWtX51eRvAriZiL4L4B4Af4eIPoHhXN9hAIeZ+SvV75/GRPyfHMD1/TUA/4eZf8TM2wB+\nD8DfwDCuTRL7t7ioN1dgfgl9byCif4KJlfqPRLX5+toU+68CuJqI9hPRWZhMLtzf4vGLQ0QE4C4A\nB5n5w2LT/QBurcq3Arhvcd8+wMzvZeYrmfkqTCb3/jcz/2MM5/qeBPA9IrqmqroBk6Qbn0X/r+8x\nAK8monOqf6c3YDLJPoRrk8T+Ld4P4E1EdBYRXQXgakwWePaKKl38rwO4hZmfF5vs18fMrX0A/DIm\nq2sfB3B7m8du6Hp+ERMv+1EAj1SfmzBJefk5TGbPHwBw4arPtcC1vhbA/VV5MNcH4OcAfAXA1zEZ\n/V4wlOsD8K8xuXl9A5PJy719vjZMni6/D2ALk/m/f7rsegC8t9KaxwD8vVWff8b1vRXAdzBJtDXV\nl4/kXp8vqnIcxxkBXY87dRzHcQrgYu84jjMCXOwdx3FGgIu94zjOCHCxdxzHGQEu9o7jOCPAxd5x\nHGcEuNg7juOMgP8PvaTgGeK/oQgAAAAASUVORK5CYII=\n", "text": [ "" ] } ], "prompt_number": 34 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Neumann BC" ] }, { "cell_type": "code", "collapsed": false, "input": [ "def diffusion_solve_implicit_neumann(U, dx, dt, nx, nt, D, NBC = [0., 0.]):\n", " alpha = D*dt/dx**2\n", " n_unknowns = nx + 2\n", " d2x = eye(n_unknowns)-d2matrix(n_unknowns)*alpha\n", " d2x[0, 1] -= alpha\n", " d2x[-1, -2] -= alpha\n", " LU = factorized(d2x.tocsc())\n", " for it in range(0,nt-1):\n", " b = U[it,:]\n", " b[[0,-1]] -= alpha*2*dx*array(NBC)\n", " U[it+1,:] = LU(b)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 35 }, { "cell_type": "code", "collapsed": false, "input": [ "U, dx, dt, x, t = diffusion_init(maxx, maxt*50, D, nx, SC*20)\n", "U[0,:] = 0.\n", "U[0, (xmaxx*0.2)] = 1." ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "diffusion_solve_implicit_neumann(U, dx, dt, nx, len(t), D)\n", "pcolormesh(U, rasterized=True, vmin=-0, vmax=.5)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] } ], "metadata": {} } ] }