{ "cells": [ { "cell_type": "code", "execution_count": 217, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Introduction to finite differences\n", "\n", "The method of finite differences is used, as the name suggests, to transform infinitesimally small differences of variables in differential equations to small but finite differences, thus enabling solution of these equations by means of numerical calculations in a computer. This might be necessary if the equations at hand are too complicated to be solved analytically, or if the geometry of the problem setup is too complex.\n", "\n", "The starting point for the finite differences method is the definition of a derivative:\n", "\n", "$$\\frac{\\mathrm{d}f}{\\mathrm{d}x}\\Big\\rvert_{x=x_0} = \\lim_{\\Delta x \\rightarrow 0} \\frac{f(x_0 + \\Delta x) - f(x_0)}{\\Delta x}\n", "= \\lim_{\\Delta x \\rightarrow 0} \\frac{f(x_0) - f(x_0 - \\Delta x)}{\\Delta x}$$\n", "\n", "However, computers can only handle numbers of finite value, not infinitesimally\n", "small values such as $\\Delta x$ at the limit used in the definition. Thus,\n", "the derivative is approximated by dropping out the limit:\n", "\n", "$$\\frac{\\mathrm{d}f}{\\mathrm{d}x}\\Big\\rvert_{x=x_0} \\approx \\frac{f(x_0 + \\Delta x) - f(x_0)}{\\Delta x}\n", "\\approx \\frac{f(x_0) - f(x_0 - \\Delta x)}{\\Delta x}$$\n", "\n", "where $\\Delta x$ is *sufficiently small*. If we could choose arbitrary \n", "small values of $\\Delta x$, we would reach the exact value of the\n", "derivative, but in reality we are always limited by the numerical accuracy\n", "of the computers." ] }, { "cell_type": "code", "execution_count": 218, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "53\n" ] } ], "source": [ "n = 1\n", "while (1 + 2**(-n) != 1):\n", " n = n + 1\n", "print(n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The first approximation above is called the *forward difference*, second the *backward difference*. \n", "One can define also the *central difference* approximation:\n", "\n", "$$\\frac{\\mathrm{d}f}{\\mathrm{d}x}\\Big\\rvert_{x=x_0} \\approx \\frac{ f(x_0 + \\Delta x) - f(x_0 - \\Delta x)}{2\\Delta x}$$\n", "\n", "Other approximations exist, too.\n", "\n", "### Taylor series approach; accuracy\n", "\n", "Taylor's Theorem states that function $f$ can be presented with following series in the vicinity of $a$:\n", "\n", "$$ f(x) = f(a) + \\frac{df}{dx}\\frac{(x-a)}{1!} + \\frac{d^2f}{dx^2}\\frac{(x-a)^2}{2!} + \\frac{d^3f}{dx^3}\\frac{(x-a)^3}{3!} + \\ldots $$\n", "\n", "Finite differences approximation is based on truncating the series, typically after the second term, leading to \n", "\n", "$$ \\frac{f(x) - f(a)}{x-a} = \\frac{df}{dx} + \\mathcal{O}((x-a)^2) $$\n", "\n", "$\\mathcal{O}$ represents the *truncation error*, and shows that if $x-a$ (i.e. $\\Delta x$) is halved, the error is dropped to one-fourth. This applies to central difference. Forward and backward differences are *first-order accurate*.\n", "\n", "## Approximating derivatives\n", "\n", "Let's define a known function $f(x)=\\sin(2x)+2\\cos(x)$ and plot it on range $x=[-2,2]$:" ] }, { "cell_type": "code", "execution_count": 219, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xd4VGX+NvD7OzPpIb0Q0gOh94SEbgFXZLGLCIqgILKua1nX1V23+ttd19XVtaAuigVBRBHFhgqKjZoEQg0lkArpgSSkZ+Z5/0jcF5GSkJl5ptyf68qVgTmZc89h5ubJOc+cI0opEBGR6zDoDkBERNbFYicicjEsdiIiF8NiJyJyMSx2IiIXw2InInIxLHYiIhfDYicicjEsdiIiF2PSsdKwsDCVkJCgY9VERE4rKyurUikVfr7ltBR7QkICMjMzdayaiMhpiUhBZ5bjrhgiIhfDYicicjEsdiIiF8NiJyJyMSx2IiIXw2InInIxLHYiIhejZR47EZ1fRV0z8qvqcfR4I46eaERzqxkiAoMIQvw8EBfqh/gQX8SG+MJoEN1xyYGw2IkcRH1zG9bnlGHz4SpszatGXmV9p34uwNuEcX3CMCE5HJf2j0DPQG8bJyVHx2In0kgphc1HqrAqqxif7SlFQ4sZAd4mpCWGYFZaHJIj/RET7IvoIB/4eBqhlIJFtY/mC6rqUVDVgMyCanx3qBJr95TCIMCkAZG4ZXQ8JvQJg4EjebckSim7rzQ1NVXxlALkziwWhS/2leH5DYew52gteniZMG1YFK4dEYOU+OAu71pRSiG3/CTe234U72YWoaq+Bb3D/fDg5f1x+aBIiLDgXYGIZCmlUs+7HIudyL6+PViBv3+SgwNldUgI9cVdF/fBVcN7wdvDaJXHb24z47M9pXjmy0M4UlGPlPhg/H7qAKTEB1vl8UkfFjuRgymqbsDfPtmHz/eWISHUF/dN7otpQ6NgMtpmclqb2YJ3Movx9PqDqDzZjHnjEvGby/tZ7T8Qsr/OFjv3sRPZmMWi8PqmfPzr8/0QCB68vB/mT0iEl8m2BWsyGjArPQ5XD++Fx9bm4JXv8/D1wQr8e/owDIsNsum6SS/OYyeyodKaJtz66jY8+vE+jOsdhi8fuAi/vKSPzUv9VH5eJvztmiFYensaTja14YaXNmHFtkK7rZ/sjyN2IhtZv68MD7y7Ey1tFvzj2iGYmRar9SDmxL7h+Py+ifjV2zvwu9W7sb+kFn+YNhAeNtoVRPrwX5TIyiwWhf+sP4j5SzMRE+yDT+4Zj1npcQ4xMyXQ1wOvzR2FOyYk4o3NBZjz6jbUNbXqjkVWxmInsqK6plYseDML/1l/CNeNjMZ7vxiLpHB/3bF+xGgQPPLzgXhy+jBsy6vGLa9sxYmGFt2xyIpY7ERWUlbbhOkvbcaGA+X485UD8e/pwxx6BsoNKTF46ZYU5JTW4abFW1BR16w7ElkJi53ICnLL63DdC5tQVN2A1+aOwm3jEh1i18v5TB4YiVfnjEJBVQNm/HczyuuadEciK2CxE3VTVsFxXP/iZjS3WbDyzjGY2Pe8F5F3KOOTw7B0XhpKapow59UM1DRyn7uzY7ETdcOWI1WYvWQrgn09sPoXYzE4OlB3pAsyKiEEL81OQW55He54IxNNrWbdkagbWOxEF2hjbiXmvrYNvYJ88M6dYxAX6qs7Urdc1DccT904HBkF1bj7re1oM1t0R6ILxGInugDfHKzA7a9nICHUD28vGI2IANc4Ve6Vw3rh0asGYX1OOf7+aY7uOHSB+AEloi7acqQKC5ZmIincH8vnpyPEz1N3JKuaPSYBeZUNeHVjHvr37IEZo+J0R6Iu4oidqAuyi05g3usZiA3xxbJ5aS5X6j/4/dT+mJAchj98sAeZ+dW641AXsdiJOml/aS3mvLoNIf6eWDYvHaH+Xroj2YzJaMDzM0ciJtgXC5dloaSmUXck6gIWO1EnFB9vwK1LtsHbw4Dl80a7xeXnAn098PKtqWhsMeOeFTt4MNWJsNiJzuNEQwvmvLoNja1mLL093elnv3RFnwh//OO6IcjIP46n1x/UHYc6icVOdA5NrWbMfyMTRdWNePnWVPTr2UN3JLu7eng0bhoVixe+PoxvD1bojkOdwGInOguLReH+ldnIKjyOp2YMw+ikUN2RtPnzlYOQHOGP+1dmo7yWpx1wdCx2orN48osDWLunFI9MHYBpQ3vpjqOVj6cRi2aNxMnmNjy8ejd0XFKTOo/FTnQG72UV44WvD2NmWhzmjU/UHcchJEf2wENT+uOr/eV4N7NYdxw6h24Xu4jEisgGEckRkb0icq81ghHpsi2vGg+v3oWxvUPx6NWDnOIsjfYyd2wC0hND8OjH+1B8vEF3HDoLa4zY2wA8oJQaAGA0gF+KyEArPC6R3RUfb8DCZVmIDfbFizen8LJxpzEYBE9OHwalFH67ahcsFu6ScUTdftUqpUqUUts7btcByAEQ3d3HJbK3xhYz7nwzC61tFrwyJxWBvh66Izmk2BBf/GHaQGw6XIXlWwt0x6EzsOpwREQSAIwAsPUM9y0QkUwRyayo4JQpcixKKTy8ehf2ldTimZnDHe5ydo7mplGxGN8nDP/67ADKOEvG4Vit2EXEH8B7AO5TStWefr9SarFSKlUplRoe7lwXIiDX9/J3R7Am+xh+87N+uLR/pO44Dk9E8LdrBqPFbMFfPtyrOw6dxirFLiIeaC/15Uqp1dZ4TCJ72ZRbiX+u3Y8rBvfEXRf31h3HaSSE+eGeSclYu6cU6/eV6Y5Dp7DGrBgBsARAjlLqqe5HIrKfkppG/GrFDiSG+eGJ6cM4A6aL7piQhL6R/vjTmj2ob27THYc6WGPEPg7AbACXikh2x9dUKzwukU21tFlw1/LtaGo147+zU+DvxcsTdJWnyYDHrhuCYzVNeObLQ7rjUIduv5KVUt8D4DCHnM7fP9mHHYUn8MLNI9Enwv3OAWMtKfEhmJEai1e/z8ONqbHoE8EDz7pxki65pY92HsMbmwswf3wipg6J0h3H6T04pR98PI149ON9PN2AA2Cxk9s5UnESv1u9GynxwXjoiv6647iEMH8v3D+5L749WIH1OeW647g9Fju5laZWM+5avh0eRsFzM0fwk6VWNHtMPJIj/PF/H+9DU6tZdxy3xlc1uZW/frQX+0vr8NSM4egV5KM7jkvxMBrw5ysHobC6AUu+z9Mdx62x2MltrMk+ihXbivCLi3vjkn4RuuO4pPHJYbhsYCRe/PowKk82647jtljs5BbyK+vx+4796g9c1ld3HJf28BX90dhqxjPrOf1RFxY7ubzmNjPuXrEdJqMBz84cARP3q9tU73B/zEqLw1vbCnG44qTuOG6Jr3Byef9cux97jtbiyenDEM396nZx7+Rk+HgY8fja/bqjuCUWO7m0L3PK8NrGfMwdm4DLBvLkXvYS5u+FX1zcG1/sK8O2vGrdcdwOi51cVlltE37z7k4MjArA76Zyvrq93T4uET0DvPHY2hx+aMnOWOzkkswWhfvezkZTqwXPzRoBL5NRdyS34+NpxL2Tk7Gj8AS+2s8PLdkTi51c0kvfHMbmI1X461WD0JsXzdDmhpQYxIf64onPD/AyenbEYieXs73wOJ5adxDThkZhemqM7jhuzcNowK8v64v9pXX4ZHeJ7jhug8VOLqWuqRX3vr0DPQO88fdrh/D86g7gyqG90L9nDzy17iDazBbdcdwCi51cyh8/2IOjxxvx7MzhCPThxagdgcEg+PVlfZFXWY9VWcW647gFFju5jPd3FOOD7GO4d1JfpMSH6I5Dp7hsYCSGxQbhua9y0dLGUbutsdjJJRRWNeCPH+zFqIRg/PISXrfU0YgI7pucjKMnGvHedo7abY3FTk6v1WzBPW/vgAjw9IzhPGWAg7q4bziGxQZh0QaO2m2N7wByes9+eQjZRSfw2HVDEBPsqzsOnYWI4L5JySg+3ojVHLXbFIudnNq2vGos2pCL6SkxmDa0l+44dB4X9wvHsJhAPL8hF62cIWMzLHZyWjWNrbh/ZTbiQnzxl6sG6Y5DnSAiuHcyR+22xmInp6SUwu/f342y2iY8c9MI+HmZdEeiTrqkXwSGxgRi0YbDnNduIyx2ckrvZhXjk10l+PXP+mJYbJDuONQFIoK7L+mDwuoGfhrVRljs5HTyKuvxlw/3YnRSCO6cyKmNzmjygEj0jfTHCxsO8xwyNsBiJ6fS0mbBvW/vgIfRgKdnDIfRwFMGOCODQXDXxX1woKwOX/LMj1bHYien8u91B7CruAaPXz8EUYG8GpIzmzY0CrEhPnh+Qy7P125lLHZyGt8fqsR/vzmCmWmxmDI4Sncc6iaT0YCFF/XGzqIT2Hy4Snccl8JiJ6dQdbIZ97+Tjd7hfvjjtIG645CVXD8yBhE9vPD8hlzdUVwKi50cnlIKD67ahZqGVjw3cyR8PTm10VV4exgxb3wiNh2uwu7iGt1xXAaLnRzeG5vy8dX+cvxuan8M7BWgOw5Z2cz0OPTwMuG/3x7WHcVlsNjJoe09VoN/fLofl/aPwNyxCbrjkA0EeHtg1ug4fLq7BIVVDbrjuAQWOzmshpY2/GrFDgT5euCJG4byakgu7PZxiTAaBEu+P6I7iktgsZPD+vOavcirrMd/bhqOUH8v3XHIhiIDvHHN8GiszCxCdX2L7jhOj8VODmlN9lG8m1WMX17cB2N7h+mOQ3awYGISmlotWLo5X3cUp8diJ4eTV1mP36/ejZT4YNw3OVl3HLKT5MgemNQ/Aks3F6Cp1aw7jlNjsZNDaWo14+63tsPDZMBzM0fwakhu5o6JSaiub8H7O47qjuLU+K4hh/LYpznYe6wWT94wDL2CeMoAd5OeGILB0QFY8n0eTw7WDSx2chif7SnBG5sLMH98IiYPjNQdhzQQEcwfn4Tc8pP45lCF7jhOyyrFLiKviki5iOyxxuOR+ymoqseDq3ZhWEwgfjulv+44pNHUIVHoGeCNJd/l6Y7itKw1Yn8dwBQrPRa5maZWM+5avh0GETw/ayQ8TfxF0p15mgyYMzYB3+dWIqekVnccp2SVd5BS6lsA1dZ4LHI///fxPuw9VounbhyG2BBf3XHIAcxKi4OPhxFLvueo/UJwaERarck+iuVbC3HnRUmYNID71aldoK8HbkyNwZrsoyiva9Idx+nYrdhFZIGIZIpIZkUFD4oQcKC0Dg+/txujEoLx4M/66Y5DDmbO2AS0mhXe2lqoO4rTsVuxK6UWK6VSlVKp4eHh9lotOajaplYsXJYFf28TFs0ayfnq9BNJ4f64pF84lm0pRHMbP7DUFXw3kd1ZLAoPvLMTRdUNeOHmkYgI8NYdiRzUbeMSUXmyGZ/sKtEdxalYa7rjCgCbAfQTkWIRmWeNxyXX9OI3h7FuXxl+P3UARiWE6I5DDmxCchj6RPjjtY35vC5qF1hrVsxMpVSUUspDKRWjlFpijccl17PhQDme/OIArhzWC7eNS9AdhxyciGDu2ATsPlqDrILjuuM4De6KIbvJq6zHPSt2oH/PADx+/RCeX5065bqR0QjwNuG1jfm6ozgNFjvZRV1TK+5YmgmTQbB4dgqvW0qd5utpwk1pcfhsbymOnWjUHccpsNjJ5iwWhftX7kReZT0W3TySH0KiLps9Oh4WxamPncViJ5t74osDWJ9Thj/8fAAvmkEXJDbEF5P6R2LFtkKeq70TWOxkU6uyivHi14cxKz2OF6OmbpkzNh5V9S34dDenPp4Pi51sJiO/Gr9bvQvj+oTir1cN4sFS6pbxfcKQFO6HNzYX6I7i8FjsZBMFVfW4880sxAb74oVZKfDgJ0upm0QEc8YkYGfRCWQXndAdx6Hx3UZWV13fgrmvZUAphSVzRyHQ10N3JHIR142Mhp+nEUs35euO4tBY7GRVTa1m3LE0E0dPNOKVOalIDPPTHYlcSA9vD9yQEoOPd5Wg8mSz7jgOi8VOVmO2KNy/MhvbC4/jmRnDkRLP0wWQ9c0eE48WswUrM4p0R3FYLHayCqUU/rRmD9buKcUjUwfgiiFRuiORi+oT0QNje4fira2FMPOC12fEYiereHrdQSzfWoiFF/XG/AlJuuOQi5s9Oh5HTzTiq/3luqM4JBY7ddtrG/Pw7Fe5mJEai4em8IIZZHuXDYxEZIAX3tzCqY9nwmKnbnknswh//WgfLh8Uib9fO5hz1ckuTEYDZqXF49uDFcirrNcdx+Gw2OmCrck+iofe24UJyWF45qYRvAoS2dXMtFiYDILlHLX/BN+JdEE+3V2CX7+zE6MTQ7F4diq8PYy6I5GbiQjwxuWDe+LdrGI0tvD8MadisVOXrd1dgntW7MCI2CC8MicVPp4sddJj9uh41DS24qNdx3RHcSgsduqSNdlHcfeKHRgWG4TXbhsFPy+eV530SU8MQXKEP3fHnIbFTp32TmYR7luZjVEJwVh6exp6ePNUAaSXiODm9DjsLK7B7uIa3XEcBoudOuW1jXn47apdGN8nDK/NTeNInRzGdSkx8PEwYhlH7f/DYqdzUkrh8c/2/29K48u3cp86OZYAbw9cPbwX1uw8iprGVt1xHAKLnc6qzWzBb1ftwotfH8bMtDi8cHMKZ7+QQ7pldDyaWi14f3ux7igOgcVOZ1TT2IrbXs/Au1nFuHdSMv5x7WAYDfzwETmmwdGBGBYbhGVbC6EUzx/DYqefKKiqx3UvbMTmw1X41/VDcf9lffmJUnJ4t6THIbf8JLbmVeuOoh2LnX5k0+FKXLNoI6rqW7BsfjpuHBWrOxJRp1w5rBcCvE1YvrVQdxTtWOwEoP0g6YtfH8Ytr2xFiJ8nPrhrHEYnheqORdRp3h5G3JASi8/28CIcLHZCTWMrFryZhcc/248rhkRhzd3jkcArH5ETmpUeh1azwruZ7n0QlcXu5rblVWPqM99hw/5y/HHaQDw/cwT8OUednFSfCH+MTgrBW9sKYHHji3Cw2N1Uq9mCJz7fj5sWb4bJKHhn4RjMG5/Ig6Tk9G5Oj0dRdSO+y63UHUUbDs3c0O7iGjz03i7sK6nFjakx+NOVgzhKJ5dx+aCeCPXzxPItBbiob7juOFrw3exGGlvMeHr9Qbzy3RGE+Xvhv7NTcPmgnrpjEVmVp8mA6amxePm7IyipaURUoI/uSHbHXTFuQCmFj3Yew6R/f43F3x7BjFGxWPfri1jq5LJmpcXBbFFYmVGkO4oWHLG7uJ1FJ/C3T/YhI/84BkQF4OkZw5HOaYzk4uJCfTGxbzhWZhTh7kv6uN3VvVjsLiqnpBZPrTuIdfvKEOrniceuG4IbU2N5WgByG7PS4rBwWRY2HKjAZQMjdcexKxa7i9leeBwvf3sEa/eUooe3CfdP7ovbxyfw3OnkdiYNiEBEDy+8tbWAxU7Op6XNgvU5ZVjyfR6yCo4jwNuEX13aB/PHJyHQl4VO7snDaMCMUbF4fkMuio83ICbYV3cku2GxO7GDZXV4N7MIq7cfRVV9C2JDfPCXKwdiemosL4RBBOCmtDgs2pCLlRlFeOBn/XTHsRu++52IUgoHyuqwdncpPt1dgkPlJ2EyCCYPiMSMtFhMTA7nPnSiU0QH+eDifhF4O6MI90xKhoebHES1SrGLyBQAzwAwAnhFKfVPazwuAaU1TcjIr8Z3hyrw3aFKlNQ0QQQYlRCCv141CFOHRCG8h5fumEQOa1ZaHOYvzcSXOWWYMjhKdxy76Haxi4gRwCIAlwEoBpAhIh8qpfZ197HdiVIKpbVN2F9ah5ySWuw7VosdhSdw9EQjACDA24TxyWG4Jzkck/pHICLAW3NiIudwSf8I9Ar0xvKthSz2LkgDkKuUOgIAIvI2gKsBuH2xWywKDa1mNDS3oa65DTWNrahpaEV1fQtKa5tQXtuEYzVNKKxqQGF1Axpbzf/72eggHwyPC8K88YlIiQ/GoF4BbjcXl8gajAbBjFFxeHr9QRRWNSAu1PUPolqj2KMBnPrxrmIA6VZ4XKswWxQq6ppRUtOIstomVNW34Hh9C6rrW3GyuRX1zWacbG5DU6sZLWYLmlstMFsU2izt3y0KUFA49Wpbp195y6LU/5Y1WyxoabOg1azQYracM1uAtwmRAd6ID/XF+OQwJIT6ol/PAPTr2QOBPpzNQmQtM0bF4tmvDmFFRiEemtJfdxybs0axn+lo3U/OlykiCwAsAIC4uDgrrPbH2swW5FacxN6jtdhfWosjFfXIq6xHYXUD2s5w+k4/TyN6eHvAz8sIfy8TvDzav4f6GWAyGGA0CAwGgVEAEWl/kqc8046/gUj7X///5QWeJgM8jAZ4mQzw8zLC19MEPy8jgnw9EeTjgRA/T0T08IaPJy8MTWQPPQO9cWn/CLybWYT7J/eFp8m1f/u1RrEXAzj1+mkxAI6dvpBSajGAxQCQmpra7RMl1ze3ISO/GtvyqpGRX41dxTVobmsfIXuZDEgM80O/nj1w+eCeiA7yQVSgNyIDvBHm74VgPw94mViqRO5kVnoc1u0rwxf7SjFtaC/dcWzKGsWeASBZRBIBHAVwE4BZVnjcnyioqseXOeXYcKAcW49Uo8VsgckgGBQdiFtGx2NIdCAG9QpAUrg/p/0R0Y9MTA5HdJAP3tpayGI/H6VUm4jcDeBztE93fFUptbfbyc7gpW8OY8W2IvSJ8MecsfGY2DccKfHB8PXkdHwiOjejQTAzLRZPfnEQeZX1SHThyz+KOv1IoB2kpqaqzMzMLv9cYVUDALjFUW0isr7y2iaM/edXuG1cAh75+UDdcbpMRLKUUqnnW86pjiDEhfqy1InogkUEeOOygZFYlVWMplOmF7sapyp2IqLumpUeh+MNrfh8b6nuKDbDYicitzKudxjiQ32xfGuh7ig2w2InIrdiMAhmpsVhW141DpXV6Y5jEyx2InI701Ni4GEUvLXNNUftLHYicjuh/l6YMjgK77noQVQWOxG5pVlpcahtasNHO3/yQXmnx2InIrc0OikESeF+Lrk7hsVORG5JRHBzejx2FJ7A3mM1uuNYFYudiNzW9SOj4WUy4C0Xm/rIYicitxXk64lpQ3vhgx1HcbK5TXccq2GxE5Fbu3l0HOpbzFiTfVR3FKthsRORWxsRG4QBUQFYtqUQOk6KaAssdiJya+0HUeOQU1KL7KITuuNYBYudiNzeNSOi4edpxLItrnEQlcVORG7P38uEa0ZE4+Ndx3CioUV3nG5jsRMRAbhldDya2yxYlVWsO0q3sdiJiAAMiApASnwwlm8thMXi3AdRWexERB1uGR2HvMp6bDpcpTtKt7DYiYg6XDE4CsG+HnhzS77uKN3CYici6uDtYcSNqbFYn1OO0pom3XEuGIudiOgUs9LjYFEKK5z4rI8sdiKiU8SH+mFicjhWbCtEq9miO84FYbETEZ1m9uh4lNc144u9ZbqjXBAWOxHRaS7pH4HoIB8s3ZyvO8oFYbETEZ3GaBDcMjoeW/OqcbCsTnecLmOxExGdwYxRsfA0GfDm5gLdUbqMxU5EdAYhfp6YNjQKq7cXo66pVXecLmGxExGdxa1jElDfYsb7O5zrIhwsdiKisxgeG4ShMYFYurnAqS7CwWInIjqHW8ckILf8pFOdP4bFTkR0DtOGRiHEzxOvb8rXHaXTWOxEROfg7WHEzLRYfJlThqLqBt1xOoXFTkR0HreMjoeIYNkW55j6yGInIjqPqEAfXD4oEm9nFKGxxaw7znmx2ImIOmHOmATUNLZiTbbjT31ksRMRdUJaYgj69+yB1zflO/zURxY7EVEniAhuG5eA/aV12HKkWnecc+pWsYvIdBHZKyIWEUm1VigiIkd09fBoBPt64LWNebqjnFN3R+x7AFwH4FsrZCEicmjeHkbcnB6PdTllKKxy3KmP3Sp2pVSOUuqAtcIQETm62WPiYRTBG5vzdUc5K+5jJyLqgsgAb/x8aBRWZhQ57Fkfz1vsIrJeRPac4evqrqxIRBaISKaIZFZUVFx4YiIizW4bl4iTzW1YlVWsO8oZmc63gFJqsjVWpJRaDGAxAKSmpjr2XCEionMYHhuEEXFBeH1TPm4dkwCjQXRH+hHuiiEiugDzxieioKoB63Mc74LX3Z3ueK2IFAMYA+ATEfncOrGIiBzblEE9ER3kg1e+O6I7yk90d1bM+0qpGKWUl1IqUil1ubWCERE5MpPRgNvGJSAj/ziyi07ojvMj3BVDRHSBZoyKRQ8vk8ON2lnsREQXqIe3B2amx2HtnlIUH3ecDyyx2ImIumHO2AQAwOsb87XmOBWLnYioG6KDfPDzIVF4O6MItQ7ygSUWOxFRNy2YmISTzW1YvqVQdxQALHYiom4bHB2I8X3C8OrGPDS36b/CEoudiMgKFl7UGxV1zXh/u/4rLLHYiYisYFyfUAzqFYDF3x6B2aL3rCksdiIiKxARLLyoN45U1mPdPr2nGWCxExFZyRWDeyI2xAcvfXNY63VRWexERFZiMhqwYEISsotOYPORKm05WOxERFY0PTUWYf5eWLQhV1sGFjsRkRV5exixYGIiNuZWYUfhcS0ZWOxERFZ2c3o8gnw9tI3aWexERFbm52XCbWMTsT6nHPuO1dp9/Sx2IiIbmDs2Af5eJiz62v6jdhY7EZENBPp6YPaYeHy6uwS55Sftum4WOxGRjcwfnwhvkxHPfXXIrutlsRMR2UiovxduHRuPD3ceQ255nd3Wy2InIrKhOyf2ho+HEc98ab997Sx2IiIbCvHzxJyxCfh41zEcLLPPqJ3FTkRkYwsmJMHXw4hnv7TPvnYWOxGRjQX7eWLuuAR8srvELqN2k83XQEREuGNCEnYV16ClzWLzdbHYiYjsIMjXE2/OS7fLurgrhojIxbDYiYhcDIudiMjFsNiJiFwMi52IyMWw2ImIXAyLnYjIxbDYiYhcjCil7L9SkQoABRf442EAKq0Yx1qYq2uYq2uYq2scNRfQvWzxSqnw8y2kpdi7Q0QylVKpunOcjrm6hrm6hrm6xlFzAfbJxl0xREQuhsVORORinLHYF+sOcBbM1TXM1TXM1TWOmguwQzan28dORETn5owjdiIiOgeHL3YReUJE9ovILhF5X0SCzrLcFBE5ICK5IvKwHXJNF5GMxZMVAAAEkElEQVS9ImIRkbMe4RaRfBHZLSLZIpLpQLnsvb1CRGSdiBzq+B58luXMHdsqW0Q+tGGecz5/EfESkZUd928VkQRbZelirrkiUnHKNppvp1yviki5iOw5y/0iIs925N4lIiMdJNfFIlJzyvb6kx0yxYrIBhHJ6Xgv3nuGZWy7vZRSDv0F4GcATB23Hwfw+BmWMQI4DCAJgCeAnQAG2jjXAAD9AHwNIPUcy+UDCLPj9jpvLk3b618AHu64/fCZ/h077jtph2103ucP4C4AL3XcvgnASgfJNRfA8/Z6PZ2y3okARgLYc5b7pwJYC0AAjAaw1UFyXQzgYztvqygAIztu9wBw8Az/jjbdXg4/YldKfaGUauv44xYAMWdYLA1ArlLqiFKqBcDbAK62ca4cpdQBW67jQnQyl923V8fjv9Fx+w0A19h4fefSmed/at5VACaJiDhALi2UUt8CqD7HIlcDWKrabQEQJCJRDpDL7pRSJUqp7R236wDkAIg+bTGbbi+HL/bT3I72/+VOFw2g6JQ/F+OnG1IXBeALEckSkQW6w3TQsb0ilVIlQPsLH0DEWZbzFpFMEdkiIrYq/848//8t0zGwqAEQaqM8XckFANd3/Pq+SkRibZypsxz5PThGRHaKyFoRGWTPFXfswhsBYOtpd9l0eznENU9FZD2Anme46xGl1JqOZR4B0AZg+Zke4gx/1+3pPp3J1QnjlFLHRCQCwDoR2d8xytCZy+7bqwsPE9exvZIAfCUiu5VSh7ub7TSdef422Ubn0Zl1fgRghVKqWUQWov23ikttnKszdGyvztiO9o/hnxSRqQA+AJBsjxWLiD+A9wDcp5SqPf3uM/yI1baXQxS7Umryue4XkTkApgGYpDp2UJ2mGMCpI5cYAMdsnauTj3Gs43u5iLyP9l+3u1XsVshl9+0lImUiEqWUKun4lbP8LI/xw/Y6IiJfo320Y+1i78zz/2GZYhExAQiE7X/lP28upVTVKX98Ge3HnRyBTV5T3XVqoSqlPhWRF0QkTCll0/PIiIgH2kt9uVJq9RkWsen2cvhdMSIyBcBDAK5SSjWcZbEMAMkikiginmg/2GWzGRWdJSJ+ItLjh9toPxB8xqP3dqZje30IYE7H7TkAfvKbhYgEi4hXx+0wAOMA7LNBls48/1Pz3gDgq7MMKuya67T9sFehff+tI/gQwK0dsz1GA6j5YdebTiLS84djIyKShvbOqzr3T3V7nQJgCYAcpdRTZ1nMttvLnkeLL/AIcy7a90Vld3z9MFOhF4BPTzvKfBDto7tH7JDrWrT/r9sMoAzA56fnQvvshp0dX3sdJZem7RUK4EsAhzq+h3T8fSqAVzpujwWwu2N77QYwz4Z5fvL8ATyK9gEEAHgDeLfj9bcNQJKtt1Encz3W8VraCWADgP52yrUCQAmA1o7X1zwACwEs7LhfACzqyL0b55gpZudcd5+yvbYAGGuHTOPRvltl1ym9NdWe24ufPCUicjEOvyuGiIi6hsVORORiWOxERC6GxU5E5GJY7ERELobFTkTkYljsREQuhsVORORi/h8gvJjahMNOvAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "def f(x):\n", " return np.sin(2*x) + 2*np.cos(x)\n", "\n", "x = np.linspace(-2,2,100)\n", "\n", "plt.plot(x, f(x))\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "From maths we know that $\\frac{df}{dx} = 2\\cos(2x) - 2\\sin(x)$. Let's define this as a python function and plot it, together with $f$.\n", "\n", "### Do\n", "\n", " - Define function `dfdx` that returns the analytically calculated derivative of $f$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# TODO: Define function `dfdx` here" ] }, { "cell_type": "code", "execution_count": 220, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xd0VcXax/HvnPQeSOgQQu891NCUIqLSBAGRakMswLXra7tXrvXaG4qABUGlIyCgoPQSQie0QIBAgJCQhPQ27x87KCqQdnL2OSfPZ60sErKz53c24clk9uwZpbVGCCGE87CYHUAIIYR1SWEXQggnI4VdCCGcjBR2IYRwMlLYhRDCyUhhF0IIJyOFXQghnIwUdiGEcDJS2IUQwsm4mtFocHCwDg0NNaNpIYRwWDt37ryota5U2HGmFPbQ0FAiIiLMaFoIIRyWUupkUY6ToRghhHAyUtiFEMLJSGEXQggnY8oYuxBCFFVOTg6xsbFkZmaaHcVmPD09qVmzJm5ubiX6einsQgi7Fhsbi5+fH6GhoSilzI5T5rTWJCQkEBsbS506dUp0DhmKEULYtczMTIKCgspFUQdQShEUFFSq31CksAsh7F55KepXlPb1ylCMEHZGa01CWjbH49OIv5xFYno2l9Kyyc3LB6VQgJ+nK0G+7gT7elCzgjchFb1xsZSv4ieuTwq7ECa7lJZNxMlLRMQkEnnqEkfOp5KckVOsc3i6Wahf2ZcWNQLpUi+IzvWCCPb1KKPE5c8HH3zAp59+Stu2bZkzZ47ZcQolhV0IExy7kMqag+dZc/Acu04noTW4uSha1Ajg9pbVqFfJl7qVfKgW4EUFHzcqeLvj5mKMnGqtScnMJSE1i4up2cQkpHHk3GUOn7/MT3vOMnf7KQCa1/DnjpbVub1VdWoEepn5ch3eJ598wsqVK0t8M9PWpLALYSNJ6dks2X2WH3eeZv+ZFMAovpN7NaBLvWBa1gzA082l0PMopQjwciPAy426laBDnYp/fC43L599Z5LZHJ3A6oPneW3lIV5beYgOdSoyrksofZtWwdVFbq0Vx8SJEzl+/DgDBgxgwoQJTJ061exIhVJaa5s3GhYWpmWtGFFeRMWl8OXGEyzdfZbsvHyaVPNnaLua3Nq8KtXLuCd9KiGdZXvPMm/HKU4nZlA9wJOxXUIZ3bk23u6O0a+LioqiSZMmALyy7AAHz6ZY9fxNq/vz0h3NbnjMlfWtgoODrdr2jVz9uq9QSu3UWocV9rWO8S8rhAPaEp3Ax+uOsfHYRbzcXLirfU1GtA+heY0Am2UICfLm4ZvqM7FHPX6NOs+sTTG8tvIQX2w4wcM31ePujiF4uBb+W4JwLFLYhbCy7ScSeWfNYbYeT6SynwdP9WvE3R1CCPR2Ny2Ti0XRt1lV+jarys6Tiby16jCvLDvIjA0neOH2JtzSrKpDTCksrGctDFLYhbCSYxcuM215FOsOxxPs68GLtzfl7o4hRRo3t6V2tSsy9/5ObDqWwKvLDzLx20h6NKzEKwOaERrsY3Y8YQVS2IUopUtp2bz3yxG+3XYKbzcXnrm1MWM7h+Llbl8F/WpKKbo2COanR7vy1ZaTvLvmCH3fW89TtzRiQngdLDIn3qFJYReihLTWzN8Zy39XRJGckcPdHUOY2rshQQ40f9zVxcK9Xetwe8tqPL9oP68uj+KXqPO8PawVNSt4mx3PbsTExJgdoViksAvHlZsNWSng4gaunuDiDjYaJ46OT+W5hfvYdiKRdrUrMG1wcxpX9bdJ22Whir8nX4xpx48Rsbyy7AD93tvAm0Nb0r9FNbOjiRKQwi7sX1YqnN4GZ3dB3G6IPwypFyAz6a/HWVyhQigENYBKDSG0G9QOB3fr9Tzz8jVfbjzO26uP4Olq4bUhLRgeVssphi6UUtzVvhad6wXx6NxdTJoTyf3d6vBUv8Z/PBwlHIMUdmGf0hMhaikcWg7Hf4e8LOPvK9aDKs2gbk/wqQSegZCfA7mZkHUZEqIh4RhE/wqb3gcXD6jdGVoOh2aDwa3k88aPx6fyxI97iDyVRJ+mVZg2uDmV/Tyt8nLtSa2K3vzwYGemLT/IFxtOsOd0Mh+PakslP8cZYirvpLAL+6E1nNkJO2bA/oVGMQ+sDe3vgwZ9oHob8Aos2rlyMuDkZoheC4dXwuKH4OdnodVI6DwJAkOKEUszb8dp/r3sIO6uFt4b3pqBras7xPTAknJ3tfDKwOa0rV2BpxfsZfAnm5g5rj0Nq/iZHU0UgTx5KsynNRz/Ddb9F2K3g7svtBoBbcdC1RalHzfXGmI2ws5ZcHCpcb7290G3x8Hnxk8SXkrL5pmFe1l14Dxd6wfzv7taUcXf+XrpN7I3Nol7v4ogMzuPT+5pS7cGlWza/rWewCwPSvPkaakHzpRStZRS65RSUUqpA0qpyaU9pyhHTm2D2bfBN4Mg5Qzc+hY8fghu+x9Ua2mdm6FKQZ1uMHQmTN5tDMts+wzebwWbP4T8vGt+WURMIre+v4G1hy7wfP8mfD2hQ7kr6gAtaway+OFwalTwYtysHczfGWt2JFEIa9wRyQUe11o3AToBDyulmlrhvMKZpV6ARRNhZl9jXPzWt+CxXdDxAfAow1/3A2rCwI9g0jYI7Qqr/w9m3mLckC2Qn6/57Pdohn++FQ83C4smhXN/97pOcYO0pGoEevHjxM50rhvEEz/uYdamE2ZHMtXLL7/M22+/zaFDh2jdujVt2rQhOjoagAcffJBNmzb95fiYmBiaN29us3ylLuxa6zitdWTB+5eBKKBGac8rnFR+vjGG/mEY7JsPXf8Fj0UaBd3VhjfnKjWEkfNgyAzjZutn3WDrpySnZXP/1xG8vvIQ/ZpV5adHu9p0bRd75ufpxpfjwujXrCqvLDvIe78cwYyhXHuyePFiBg4cyK5du6hXrx4A27Zto1OnTqbmsuocJqVUKNAG2GbN8wonkRwL3w6G5Y9D9dYwaQv0fgncTXqMXSloOQwe3g71e8HPz7D9nSHsOHqaVwY046O72+DnWbJd4p2Vh6sLH93dhqHtavLeL0d5beWhclPcp02bRqNGjejduzeHDx8mPT2d9957jxkzZnDTTTcBxrh4w4YNcXFxYefOnbRq1YrOnTvz8ccf/3Ged955hwkTJgCwb98+mjdvTnp6ulWzWm1WjFLKF1gATNFa/2NdTaXUA8ADACEhRZ+RIJzEnu9hxZOQnwu3vwftxtnsYaJC+VZmSeM3iY7yZwrfs7XyBbwbz7OffHbG1cXCm3e2xMvNhc/XH8eiFE/3a2SbWUIrn4Fz+6x7zqot4NbXb3jIzp07mTdvHrt27SI3N5e2bdvSrl07Jk6ciK+vL0888YQRb+VK+vXrB8D48eP58MMP6dGjB08++eQf55oyZQo9e/Zk0aJFTJs2jenTp+Ptbd2nfK3SY1dKuWEU9Tla64XXOkZr/bnWOkxrHVapkm3vqgsTZafDkodh0QNQpSk8tBHCxttN0czL17y2MorJ3+9la/VxpNw5D+/MCzCjD5zdbXY8u2WxKP49sBmjOobw2e/RvLXqsFP33Dds2MDgwYPx9vbG39+fAQMGXPO4VatW0a9fP5KTk0lKSqJHjx4AjB49+o9jLBYLs2fPZvTo0fTo0YPw8HCr5y11j10ZP6a/BKK01u+UPpJwGhePwQ9j4MIB6P4k9HwWLPazMFZyRg6T5+3it8Px3NMphJfuaGY8YVntF/hmMMy+HUbMgbo9zI5ql5RS/Gdgc/I1fPJbNO6uFqb0bli2jRbSsy5Lhf1Gkp6eTlJSEtWrVycpKemGxx89ehRfX1/Onj1r7ZiAdXrs4cBo4Gal1O6Ct/5WOK9wZEdWw+c94XIcjFoAN/+fXRX1mItpDP5kExuPXmTa4Oa8OqjFn4/NB9eHe1cZM2jmDDXmvotrslgU0wY1/2PMfbaTzpbp3r07ixYtIiMjg8uXL7Ns2bJ/HLNu3bo/xtoDAwMJCAhg48aNAH/ZADs5OZnJkyezfv16EhISmD9/vtXzlrrHrrXeCNjH79XCfFobc8PXvGiMXY74DgJrmZ3qL7ZEJzDx251YFHx7X0c61Q3650H+1WH8CvhuOMwfD3d9A42lv3ItFovi9SEtSMnI4eVlBwn0dmdQG+eaGNe2bVuGDx9O69atqV27Nt26dfvHMStXrmTo0KF/fDxr1iwmTJiAt7c3t9xyyx9/P3XqVCZNmkTDhg358ssvuemmm+jevTuVK1e2Wl558lRYT242/DQFds+BpgNh0KfmzXi5ju93nOL5RfsJDfbhy7Fh1A4qJF9mCnw9EM7vN6ZH1u9lm6AOKDMnj3GzthMRc4kvxoRxU2PrFCpHefK0bdu2bNu2DTc368ykMvXJUyEAowB+N8wo6j2egaGz7aqo5+drXl95iKcX7KNzvSAWTupSeFEH8PSHexZAcCOYN8pYmkBck6ebC1+MCaNxNT8e/i6SfbHJZkeyqcjISKsV9dKSwi5KL+UszLrVKHqDPoWbngWL/XxrZebk8cjcSD77PZpRHUOYNa49/sWZn+5dEcYsNhYOmzsSzh8su7AOzs/TjZlj21PB250JX+3gTFKG2ZHKJfv53yccU/wRY2rgpRi4+wdofbfZif7iYmoWIz7fysr95/i/25rw6qDmuJZkbXGfYKPn7uYF390Fl89ZP6yTqOzvyazx7cnMyWP8rO0kZ+SU+pzOPJXyWkr7eqWwi5KL22P01POyYNxyuxt/jo5PZcgnmzl0LoVPR7Xjvm51S/cQTWAtuPt7SE8wbqpmp1kvrJNpWMWP6fe048TFNB6eE0luXn6Jz+Xp6UlCQkK5Ke5aaxISEvD0LPmCc3LzVJTMyS1Gz9UzAEYvNqYI2pEdMYnc/3UELkoxY2wYbUIqWO/kh1fCvLuhUX9jtowdDTvZmx92nOapBXsZHx7KS3c0K9E5cnJyiI2NJTMz08rp7Jenpyc1a9b8x5h9UW+eykYbovii18Lcu4153mMWG3/akeV745j6w25qBHoxe3z7ot0kLY5Gt0LfabDqWdj4P+PhK3FNd7WvxaFzl5m56QSNq/oxvH3xlxNxc3OjTp06ZZDOeUlXQxTPkdXw3QgIqg/jV9pdUZ+x4TiPzI2kRY0AFj5UxJkvJdHpIWgxDNZOg2O/lk0bTuK5/o3p3rAS/7d4PztiEs2OUy5IYRdFd2i5MQRRuQmMXQq+9rPmT36+5t/LDvLq8ihuaVqVOfd1pIKPe9k1qBTc8b5xLRbcC5dOll1bDs7VxcKHI9tQq4I3D30byfmU8jOkYhYp7KJoopYZ675UawVjlhhTAO3ElemMMzedYFyXUD4e1RZPNxssX+DuA8O/NXZg+mEM5GaVfZsOKsDLjemj25GenctD3+4kO7fkN1NF4aSwi8JF/QQ/jjM2kx69qOgbSttAUno2Y77czop953i+fxNeuqMpLrbc6SionjF3P243rP2P7dp1QA2q+PHW0FZEnkri1eXyLEBZksIubuzQcvhxLFRrbczj9vQ3O9EfziRlMPSzLew+ncQHI9twf/dSTmcsqSa3Q9gEY42c6HW2b9+B3NayGvd3q8PXW06yMFL2Ti0rUtjF9R3+GX4Yawy/jF5oTG20EwfOJjP4402cT8nkqwkdGNCqurmB+k4zlh1YNBHSLpqbxc493a8xHetU5LlF+zhy/rLZcZySFHZxbcd+gR9GQ9XmcI99FfUNR+MZPn0rLhbFgoe60LneNVZntDV3bxj6JWQkwpJHjFUuxTVduZnq6+HKpDmRpGfnmh3J6UhhF/90Yr2x4FWlRkZRt6Mx9QU7Yxk/awc1K3ixaFI4Dav4mR3pT1VbQO9X4MhKYzE0cV2V/T15f0QbouNTeWHxAbPjOB0p7OKvTm4xHpevUAdG28/sF601H/x6lMd/3EOHOhX5YWJnqgaU/JHrMtNxItQOh5+fheQzZqexa+H1g3ns5gYsiIzlh4jTZsdxKlLYxZ/O7IQ5w4xNJsYuBR87GOIAcvLyeXbhPt5Zc4QhbWswe3yH4q3OaEsWCwz8yNi0e9ljMiRTiMd6NaBLvSBeWnKAYxdkvN1apLALw7n98M2QgiVql4Kv9XZzKY3UrFzu+yqCeTtO8+jN9fnfsFa4u9r5t23FutD7ZeM+hQzJ3JCLRfHe8NZ4ubvw6NzdZObkmR3JKdj5/xBhE/FHjF2C3H2MnnqAfWxrFpecwbDPtrDx2EVeH9KCx/s2Mmc6Y0m0v//PIZmUOLPT2LXK/p68NbQlUXEpvPHzIbPjOAUp7OVd4gn4egAoi9FTrxBqdiIAouJSGPzxZk4npjNzXHtGdCj+4lGmslhgwIfG06g/P212GrvXq0kVxnUJZdamGNYeOm92HIcnhb08S441inpuprFMgJ0svbvu0AWGfroZgB8e7EyPhvazJk2xBNWDHk/CwSXGMwHihp65tTGNq/rxxI97ib8syzOUhlUKu1JqplLqglJqvzXOJ2wg9YIx/JKRZCwTUKWp2YkA+HpLDPd+tYPQYB8WPxxO0+r286RriXSZbDy4tOIJ2ZijEJ5uLnw4sg1pWbk8vWBvudlYoyxYq8c+G+hnpXOJspaWYBT1lLMw6kdjDRiT5ebl8/LSA7y45AA3N67MDw/a6XTG4nJ1hzveg+TT8NtrZqexew2q+PHMrY1Ze+gCc7adMjuOw7JKYddarwdkoWVHkJEE3w6GhGgYORdCOpmdiJTMHCZ8FcHszTHc27UO00eH4ePhRHvA1O4CbcfAlk/gvDyMU5ixnUPp1iCYacujOB6fanYchyRj7OVJ1mWYMxTOH4QRc6BuT7MTcTIhjSGfbGZzwcyXF2638eqMttL7FWMBtRVPydz2QlgsireHtcLDzcLU73eTU4r9UssrmxV2pdQDSqkIpVREfHy8rZoVV2SnGU+UnomEYbOgQR+zE7Hx6EUGfLSJi6lZfHNvR8eb+VIc3hXh5hfg5EY4sNDsNHavir8nrw1uwZ7YZD5ZF212HIdjs8Kutf5cax2mtQ6rVMlBZzk4qpwMmDsSTm2BIZ9DkztMjaO15suNJxgzcxtV/D1Y8nC4fSzkVdbajYOqLWH1C3IjtQhubVGNQa2r8+Hao+yLTTY7jkORoRhnl5NpLOh1Yr2xIUSLoabGyczJ44kf9/Kfnw7Sp2kVFk4KL7t9Se2NxQX6vwUpZ2DD/8xO4xBeGdCcIF93Hv9RnkotDmtNd5wLbAEaKaVilVL3WuO8opRys4wt26J/NR6WaTXC1DinE9O589PNLIiMZUrvBnw6qh2+znSTtChCOkHL4camHInHzU5j9wK83Xj9zpYcOZ/Ku2uOmB3HYVhrVsxIrXU1rbWb1rqm1vpLa5xXlEJuFnw/Go6ugtvfhbajTY3z+5F47vhoI6cS0/lybBhTejfE4ow3SYui9ytgcYM1L5mdxCHc1KgyIzuE8PmG40TEyOS7opChGGd0dVG/7R1j2zaT5OVr3llzhHGztlPV35Nlj3SlV5MqpuWxC/7VIHwyRC01lkkWhXr+tiZUD/Diyfl7ZUimCKSwO5ucTPj+nj+LenvzRsXiL2cxZuY2Pvj1KEPa1GThpC6EBpeT8fTCdHkE/KrDqucgX6bzFcbXw5U3h7bkxMU03l512Ow4dk8KuzPJToe5w+HoamP4xcSivvHoRfp/sIGImEu8ObQl/7urFd7u5Ww8/UbcfaDXC3A2EvYvMDuNQwivH8yojiF8uemEDMkUQgq7s7jy8NGV2S8mDb/k5OXz+spDjJ65jQAvN5Y8Es5dYbVMyWL3Wo4wNgr/5WVjSqoo1LP9ZUimKKSwO4P0RPhmMJzaCnfOgNZ3mxLjeHwqQz/dzGe/RzOyQwjLHulK46oOvohXWbJYoO80SImFbZ+ZncYhXD0kI7Nkrk8Ku6NLiYPZt0HcHrjra2h+p80jaK35eksM/T/YwMnEdD4d1Zb/Dm6Bl7uLzbM4nDrdoEFf2PguZFwyO41DCK8fzMgOtfhiw3F2n04yO45dksLuyBKiYWZfSDplrNLY5HabRziTlMGYmdt5cckBOtUNYtWU7tzaoprNczi0Xi9BZopR3EWRPNu/CZX9PHlq/h6ycmVI5u+ksDuqM5Ewsx9kpRrb2dXtadPm8/M132w9Sd93fmfnyUu8Oqg5s8a1p4q/Eyy1a2tVmxsPLW2bDslnzE7jEPw93fjvkOYcOZ/Kx7KWzD84VmHPyzHGk8u7I6uN4RdXT5jwM9RoZ9Pmj11IZeQXW3lh8X7ahFRg1ZTu3NOptuPsR2qPbnoOdL6s2V4MNzeuwpA2Nfhk3TEOnk0xO45dcazCvuo5+LwnXIgyO4l5ds6GuSMguAHctwYqNbJZ05k5eby96jC3vr+eqLgU3ryzJd/c24FaFb1tlsFpVagN7e+D3XMgXuZpF9ULtzcl0NuNZxbuJVeW9/2DYxX2liOM/Tln9IGja8xOY1t5ucaO98smQ/1eMG4F+FW1SdNaa1YfOEffd9fz0bpj3NGyOmuf6Mld7WtJL92auj0Bbj6wbprZSRxGBR93Xh7QjL2xyczcdMLsOHbDsQp7zXZw/zqoWAe+uwu2fFw+Ni3IuGTMUd/6CXSaBCPmgoevTZo+fO4y93y5jQe+2YmHq4Xv7u/IO8NbE+zrYZP2yxWfIOg8ydj8Om6P2Wkcxm0tqtG7SRXeWXOEkwmyHDI4WmEHCKhhjCs3vs0Ymllwn3ED0VmdPwBf9IKYjTDgI+j3GriU/ROccckZPLNgL7e+v579Z1J4ZUAzVk7uRpd6wWXedrnW+WHwDIS10msvKqUUrw5qjpvFwrML98km2DhiYQfjcexhX0OvF43daGb0gngnfFhh17fwxc3GpgzjfrLJCo2X0rJ5bUUUPd/6jQWRsYztEspvT/RkbJdQXF0c89vFoXgGGAuEHV0Fp7ebncZhVA3w5Nn+TdgcncAPEafNjmM6ZcZPt7CwMB0REWGdkx3/Debfa4y93/qm8dSlo4/7ZqXCiidhz3dQp4fxNKlv5TJt8mJqFjM2nOCbLTGk5+QxuHUNpvZpKDdGzZCdBu+3gspNYOwys9M4jPx8zYgvtnIoLoVfHu9BZT/nm3qrlNqptQ4r7DjH74LV7QkProfqbWDJJPhxrGNPiTy1FT4Lhz1zocfTMHpRmRb1UwnpvLz0AN3eWMf09dHc3KQKq6Z0553hraWom8XdB7r+y1j35/jvZqdxGBaL4rUhLcjMzeeVpQfNjmMqx++xX5GfZ+xKs/ZV8Ak2lqxt3N+6bZSlnEz4/XXY9D4E1ILB06F25zJpSmtNxMlLzNp0gp/3n8OiFANaV2dSz/rUr2ybm7KiEDmZ8EFrqFAHxq9w/N9Cbejjdcd4a9VhPh/djr7NbDNzzFaK2mN3nsJ+xdndsHgSXDgATQYYe0zaaFpgiUWvheWPG1ultRlt3CD18LN6MymZOSzedYY5W09x+Pxl/D1dubtjbcZ1CaVqgPP92urwtn0OK5+EMUts/mSxI8vJy+eODzdyKT2bNf/qgb+nm9mRrKb8FnYwnlDd/AH89ga4ekD3J6DDg+BmZ8Ur6TT88pKxHnfFenDb21DvZqs2kZOXz4aj8SyMPMOag+fJys2nRY0A7ukUwh2tqssa6fYsJxM+aGM8vDR+pfTai2H36SQGf7KJUR1DeHVQC7PjWE1RC7tz/q92cYNuj0PTQbDyaVjzImyfAb1fgmZDjOVSzZSeaOxSv/0L4+OezxkzIaz0gycrN4/NxxJYuT+ONQfPcyk9hwreboxoX4s729WkZc1Aq7QjypibJ3T7F6x4Ak78Lr32YmhdK5BxXUKZvTmGwW1q0K52RbMj2ZRz9tj/LnodrHkBzu2D4EbQdQq0GGb8ALCly+eMhZ52zIDsVGh1N/R8BgJLvxHF6cR01h+N57fD8Ww+dpG07Dz8PFzp1aQyt7WsTo+GlXB3dfx75eXOlV57YIjx/Ib02ossLSuXvu+ux9vdhZ8e64qHq+MvI23ToRilVD/gfcAFmKG1fv1Gx9u8sIOxr+SBhcbSqOf3Gzco242F1qPAv3rZtas1nNkJETNh7w+QnwtN7jAWfarcpESnzM3L5+iFVPacTmL7iUS2nUjkTJKxA0+NQC96NqpEryaVCa8f7BTfzOXe9i+MXvvoRVYfqnN26w5dYPzsHUzt3ZDJvRuYHafUbFbYlVIuwBGgDxAL7ABGaq2vO9/IlMJ+hdbGOjObP4CYDaAsUK8XNBtkbHhgjamFWhsLlUUtg73fQ2I0uHlDm3ug00NQsW4RT6O5mJpNTEIah89d5tC5FA7FXebA2RQyCrYFC/Jxp2PdinSsE0R4/SDqVfKV9VucTW7Wn712GWsvtke+i2T1gfOsmNzN4Wd92bKwdwZe1lrfUvDxswBa6+uuP2pqYb9a4nHY/R3snmtsTwZQvS3U7mLMi6/eBiqEgqWQXm/WZaOQnz9gPC0YvRZSzxmfC+1mrLXddIDxVCFGwb6clUtKRg5J6TkkpGVz8XIWF1OziEvO5ExSBmeTMjiZkE5qVu4fzfh5uNKoqh/NawTQulYgLWsGUCfYRwp5ebBtOqx8Csb+ZOy6JIos/nIWvd/5nUZV/Jj3QCcsFsf9/2LLm6c1gKuf4Y0FOlrhvKWWlZtHQmo2Sek5JKVnk5yRw+WsXNKzcknLziMrJ4/M3DvJqjuQoLSjNEjeRONLW6kZNx03nQNAHhZSXCqS5FqJDIs3WrmQhwX3/Ez88pLwy0/CP+/P7bkuW/w54NGGPQEjiXRvw7nUILLW55O1dhfp2bmkZ+WRlp1L/nV+nvp5uFI90ItqgZ6E1a5AaLAPocE+NKziR/UATyni5VXbMbD+bVj/phT2Yqrk58Hz/Zvw1IK9fB9xmpEdQsyOVOasUdivVWn+UbaUUg8ADwCEhFjnwmbm5HHiYhrH49OISUgj9lIGZ5IyiEvKID41i6T0nBsHV+Dp6oK7qwU3F3/cXPrjYrkND6886urTNMqPppq+QLBOoFLeRTxzU3EhH3edT7ZyJ0ZVJdmlEfFulYhxqUOMaygXXarg6mrB1aJwtVgIdLXg7mrBw9WCj7sr3h4u+Li7EuDlhr+XKwFe7gT7uhPs60GQrzt+TjTnVliRmxeEPwYK0nfnAAAWdElEQVSr/w9ObYMQu+g7OYxhYTVZuCuW/66IolfjylR28p2+HGooZufJRDYdSyAqLoWouBROJqb/ZdXeij7u1Aj0onqgJ5X9PKnk50GwrwcVfdwI8HIn0NsNXw9XfD1c8XJ3wcPVIj1g4Tiy0+C9FsYQ4T0LzE7jcI7Hp9Lv/Q30aVKFj0e1NTtOidhyKGYH0EApVQc4A4wA7rbCef9h2Z44Zm+OITTImybV/BnUpgb1KvlSt5IPoUE++Hg457R8IQBjDZnOD8Ov/zZmWtl4S0RHV7eSL4/dXJ+3Vx9h8MHz9G5axexIZcZa0x37A+9hTHecqbW+4WLSJe2xJ6Rm4eHmgq8UcFFeZaYYvfba4TDyO7PTOJzsXGO5gZTMHFZP7e5wQ582Xd1Ra71Ca91Qa12vsKJeGkG+HlLURfnm6Q8dJ8Lh5XC+fK9gWBLurhZeu7MF51IyeXuV8+4tK48iCuFoOj5o7I268V2zkziktiEVGNOpNl9vPcnOk5fMjlMmpLAL4Wi8K0LYeNg/33gWQxTbk/0aU9Xfk2cX7iU7N9/sOFYnhV0IR9TlUbC4Guv3i2Lz9XDlPwObc+R8KtN/jzY7jtVJYRfCEflVNZao2P0dpJw1O41D6t20Cre1rMaHa49x7EKq2XGsSgq7EI4qfHLBzmEfmZ3EYb18RzO83F14buE+8q/3OLgDksIuhKOqEAothsLO2Y69z6+JKvl58PxtTdgek8jcHafMjmM1UtiFcGRdp0JOmrFImCiRYe1q0qVeEK+vOMS55Eyz41iFFHYhHFnlJtCoP2yfDlnONU5sK0op/ju4Bdl5+bywZD9mbD5kbVLYhXB0XadCxiWI/MrsJA4rNNiHf/VpyJqD51mx75zZcUpNCrsQjq5WB6jd1biJmptldhqHdW/XOrSoEcBLS/dzKS3b7DilIoVdCGfQbSpcPmtsvyhKxNXFwht3tiQpPYf/LHfs5RqksAvhDOr1gqotYdN7xhRIUSJNq/vzUM96LIw8w2+HL5gdp8SksAvhDJSCrlMg4RgcWm52Gof2yM31qV/Zl+cW7uNy5o0367FXUtiFcBZNBkKFOkav3QlmdpjFw9WFN+5sSVxKJm/8fMjsOCUihV0IZ+Hiaqwhc2YnxGwwO41Da1e7AhPC6/Dt1lNsiU4wO06xSWEXwpm0HgU+lWHje2YncXhP9G1E7SBvnl6wl/TsXLPjFIsUdiGciZsndJoI0b9C3B6z0zg0L3djSOZUYjpvOdimHFLYhXA2YfeCu5/02q2gU90gxnSuzezNMWw/4Tjr8UhhF8LZeAUaG3EcXAyJJ8xO4/Ce7teYmhW8eGr+HocZkpHCLoQz6jTJ2IhjiyzpW1o+Hq68eWcrYhLSefNnxxiSkcIuhDPyrwYth8OubyE13uw0Dq9zvSDGdQll9uYYth63/1kypSrsSqlhSqkDSql8pVSYtUIJIawgfLKxdsx2WdLXGp7qZ8ySeXL+HlKz7HtIprQ99v3AEGC9FbIIIawpuAE0vg22fyFL+lqBt7srbw9rReylDKYtjzI7zg2VqrBrraO01o4x6CREeRQ+BTKTIPJrs5M4hfahFXmgW13mbj/FukP2u5aMjLEL4cxqtYfa4bDlY8hzzHVP7M3UPg1pVMWPpxfstdvlfQst7EqpX5RS+6/xNrA4DSmlHlBKRSilIuLj5WaOEDYTPgVSYmHffLOTOAVPNxfeGd6KS+nZ/J+d7rhUaGHXWvfWWje/xtuS4jSktf5cax2mtQ6rVKlSyRMLIYqnQR+o3Aw2vQ/5+WancQrNqgcwpXdDlu+NY/HuM2bH+QcZihHC2SllzJCJj4Kjq81O4zQm9qhHWO0KvLj4AKcT082O8xelne44WCkVC3QGliulVlknlhDCqpoPgYBaxpK+wipcLIp3h7dGA4//sIe8fPsZkintrJhFWuuaWmsPrXUVrfUt1gomhLAiFzfo/Aic2gKntpmdxmnUqujNvwc2Y3tMIp/9Hm12nD/IUIwQ5UXb0eBVQXrtVja4TQ1ub1mNd9ccYdepS2bHAaSwC1F+uPtAhwfh8Aq44Jg7A9kjpRTTBregir8nj83bRYodbKcnhV2I8qTDA+DqBZs/MDuJUwnwcuODkW04m5TJ84vMnwIphV2I8sQnCNqNhb3fQ3Ks2WmcSrvaFfhXn4Ys23OWHyPMvbZS2IUobzo/bGx2veVjs5M4nYk96tGlXhAvLT3A0fOXTcshhV2I8iYwBFoMg51fQbrj7ArkCFwsiveGt8bHw4WH5kSatjGHFHYhyqPwyZCTZqz8KKyqsr8n749oQ3R8qmnj7VLYhSiPqjSFhv1g22eQnWZ2GqcTXj+YKb0asmjXGebtOG3z9qWwC1FedZ0KGYkQ+Y3ZSZzSozfXp1uDYF5aeoD9Z5Jt2rYUdiHKq5BOENIFNn8Iufa5/KwjsxSMtwf7uPPgNztJtOESv1LYhSjPuv2rYEnfH81O4pSCfD349J52xKdm8ejcSHLzbLO6phR2Icqz+r2hagtjmQFZ0rdMtKoVyKsDm7PpWAJvrbLNhnNS2IUoz5QyxtovHoFDP5mdxmnd1b4WozqGMH39cdYcPF/m7bmWeQtCCPvWZCBUqAMb34UmdxjFXljdS3c0o6KPO53qVizztqTHLkR55+IKXafA2Ug4/pvZaZyWu6uFx/s2ws/TrczbksIuhIBWI8GvOmz4n9lJhBVIYRdCgKsHdHkUYjbIRhxOQAq7EMLQbix4B8GGt81OIkpJCrsQwuDuA50eMja8jttjdhpRClLYhRB/an8/ePjLWLuDk8IuhPiTVyB0uB8OLoV42zxMI6yvVIVdKfWWUuqQUmqvUmqRUirQWsGEECbpNAncvGC9jLU7qtL22NcAzbXWLYEjwLOljySEMJVPMLS/F/bPh4Ros9OIEihVYddar9ZaX9kiZCtQs/SRhBCm6/wouLjLWLuDsuYY+wRgpRXPJ4Qwi18VaDce9syDSzFmpxHFVGhhV0r9opTaf423gVcd8zyQC8y5wXkeUEpFKKUi4uPjrZNeCFF2wieDxRU2vGN2ElFMhS4CprXufaPPK6XGArcDvfQNNvfTWn8OfA4QFhZm+00AhRDF418N2o42Nr3u/oSxCbZwCKWdFdMPeBoYoLVOt04kIYTd6DrV+FN67Q6ltGPsHwF+wBql1G6l1GdWyCSEsBcBNaHtGNj1DVw6aXYaUUSlnRVTX2tdS2vduuBtorWCCSHsRLfHQVlkDRkHIk+eCiFuLKAGtBsHu7+TGTIOQgq7EKJwXf8FygXWv2V2ElEEUtiFEIXzrwZh42H3XHka1QFIYRdCFE3XqcbTqL+/YXYSUQgp7EKIovGraqz8uPcHuBBldhpxA1LYhRBF13UquPvC2lfNTiJuQAq7EKLovCsae6Me+gnO7DQ7jbgOKexCiOLpPMnYG1V67XZLCrsQong8/Izpj9Fr4cQGs9OIa5DCLoQovvb3gn8N+OUluP7af8IkUtiFEMXn5gU3PWeMsx9cbHYa8TdS2IUQJdNqJFRuCr/+G3KzzU4jriKFXQhRMhYX6P0yJB6HnbNNDiOuJoVdCFFyDfpC7a7G06iZKWanEQWksAshSk4p6PNvSL8Im943O40oIIVdCFE6NdtB86Gw5SNIOmV2GoEUdiGENfR5BVCw5iWzkwiksAshrCGgJoQ/BgcWwqmtZqcp96SwCyGsI3wy+FWHlU9Dfr7Zaco1KexCCOtw9zGmP8bthj1zzU5TrklhF0JYT4thULO9sdRARpLZacqtUhV2pdR/lFJ7lVK7lVKrlVLVrRVMCOGALBbo/zakJ8C6aWanKbdK22N/S2vdUmvdGvgJeNEKmYQQjqx6awi7F3bMgLO7zU5TLpWqsGutr37UzAeQZd6EEHDz/xlrti9/XG6kmqDUY+xKqWlKqdPAKKTHLoQA8AqEPv+BMxGw6xuz05Q7hRZ2pdQvSqn913gbCKC1fl5rXQuYAzxyg/M8oJSKUEpFxMfHW+8VCCHsU6sRUDsc1rwAl8+bnaZcUdpKi+QrpWoDy7XWzQs7NiwsTEdERFilXSGEHbt4FD4Nh4a3wHDpuZeWUmqn1jqssONKOyumwVUfDgAOleZ8QggnE9wAej4NUUvh4FKz05gv9YJNmintGPvrBcMye4G+wGQrZBJCOJMuj0HVFrDiCci4ZHYa85zeDu82h8Mry7yp0s6KuVNr3bxgyuMdWusz1gomhHASLm4w4CNIuwirnjc7jTmy02DRRPCtYtx3KGPy5KkQouxVbw1dp8DuORD1k9lpbO+XlyExGgZ9DJ7+Zd6cFHYhhG30eAaqtoRlj9lsrNkuHP8Ntn8OHR+COt1t0qQUdiGEbbi6w5AvjGGJpY+ClWbk2bXMZFj8MAQ1gN62W6teCrsQwnYqN4ber8CRn51/A2ytjR9gl+Ng8HRw87JZ01LYhRC21eEBqHsT/PwsnD9gdpqyE/ElHFwCvV40tg+0ISnsQgjbslhgyOfgGQDfj4bMlMK/xtHE7YWfn4P6fYzpnjYmhV0IYXu+lWHoTLgU43zj7VmX4cdx4F0RBn9m/CCzMSnsQghzhIZDrxfg4GLYNt3sNNaRnw8LH4RLJ+DOL8En2JQYUtiFEObpMhka3gqrnzemBTq6ddPg8HK45b/GDy6TSGEXQpjnynh7UAP4YQzEHzE7Ucntmw8b3oY2o6HjRFOjSGEXQpjL0x/u/h5c3OG7uyA90exExXdmJyx5GEI6w23vgFKmxpHCLoQwX4XaMOI7SDkL80ZBTobZiYou/gjMGWbcEL7rG+NBLJNJYRdC2IdaHWDwp3BqizENMjfL7ESFSzoN3wwCZYHRi8G3ktmJACnsQgh70vxOuON9OLYGFtwLeblmJ7q+tIvwzWDISoV7FkJQPbMT/UEKuxDCvrQbC/1eh6hlsPgh+yzuKXEw+zZIPm3cH6jW0uxEf+FqdgAhhPiHTg9BTjr8+m/ITjXmhLt7m53KcCkGvh5o9NhHzYfanc1O9A/SYxdC2Kduj0P/t40dh74ZbB+zZS5Ewcx+kJEEY5ZCnW5mJ7omKexCCPvV4X4YNhvORsKsW+HiMfOyHFwKM3qDzofxK2y+sFdxSGEXQti3ZoPgngXG5hyf94T9C23bfn6eMST0w2io1AjuXwdVmtk2QzFJYRdC2L863WHiBmM99/njYfkTkJ1e9u0mHjfG0zf8z3iidNwKCKhR9u2WkhR2IYRjCKhpFNZOD8OOL+DjjnBoRdm0lZcLmz6AT7pA3B4Y8KHx5uZZNu1ZmVUKu1LqCaWUVkqZs5SZEKJ8cHWHfv+F8SvBwxfmjYTvhhvF1xry84yhnundYM0LUO8meHgbtB1j+jIBxVHq6Y5KqVpAH+BU6eMIIUQR1O4CD66HbZ/Bb2/A9O7GcE3nR6F+L7C4FO98GZeMefObP4KLhyG4IQz7CpoOdKiCfoU15rG/CzwFLLHCuYQQomhc3KDLo8bYd+RXsPUz+G4YeFWEBn2gQV+o1goCQ8DV469fm51mTF08txcO/wzRayE/Byo1MTYAaTqo+D8c7EipCrtSagBwRmu9RzngTzUhhBPwCoTwydDxIWMt9MM/w9HVsPd74/PKAv41jNUj83MgLwcunwMKdm3yrwkdH4RmQ6BGW4fsof9doYVdKfULUPUan3oeeA7oW5SGlFIPAA8AhISEFCOiEEIUgas7NBtsvOXnQdxuuHjUmNlyKQbyc8HiBi6uEFgbKjeFKk2hQh2nKOZXU7qEew0qpVoAvwJX5hzVBM4CHbTW5270tWFhYToiIqJE7QohRHmllNqptQ4r7LgSD8VorfcBla9qMAYI01pfLOk5hRBClJ7MYxdCCCdjtdUdtdah1jqXEEKIkpMeuxBCOBkp7EII4WSksAshhJORwi6EEE5GCrsQQjiZEj+gVKpGlYoHTpbwy4MBe5wrL7mKR3IVj+QqHnvNBaXLVltrXamwg0wp7KWhlIooypNXtia5ikdyFY/kKh57zQW2ySZDMUII4WSksAshhJNxxML+udkBrkNyFY/kKh7JVTz2mgtskM3hxtiFEELcmCP22IUQQtyA3Rd2pdRbSqlDSqm9SqlFSqnA6xzXTyl1WCl1TCn1jA1yDVNKHVBK5SulrnuHWykVo5Tap5TarZQq80Xoi5HL1terolJqjVLqaMGfFa5zXF7BtdqtlFpahnlu+PqVUh5Kqe8LPr9NKRVaVlmKmWucUir+qmt0n41yzVRKXVBK7b/O55VS6oOC3HuVUm3tJFdPpVTyVdfrRRtkqqWUWqeUiir4vzj5GseU7fXSWtv1G8YOTa4F778BvHGNY1yAaKAu4A7sAZqWca4mQCPgN4x16K93XAwQbMPrVWguk67Xm8AzBe8/c61/x4LPpdrgGhX6+oFJwGcF748AvreTXOOAj2z1/XRVu92BtsD+63y+P7ASUEAnYJud5OoJ/GTja1UNaFvwvh9w5Br/jmV6vey+x661Xq21zi34cCvGTk1/1wE4prU+rrXOBuYBA8s4V5TW+nBZtlESRcxl8+tVcP6vCt7/ChhUxu3dSFFe/9V55wO9VNlv7GvGv0uRaK3XA4k3OGQg8LU2bAUClVLV7CCXzWmt47TWkQXvXwaigBp/O6xMr5fdF/a/mYDxU+7vagCnr/o4ln9eSLNoYLVSamfBvq/2wIzrVUVrHQfGNz5X7b71N55KqQil1FalVFkV/6K8/j+OKehYJANBZZSnOLkA7iz49X2+UqpWGWcqKnv+P9hZKbVHKbVSKdXMlg0XDOG1Abb97VNler2sttFGadxow2yt9ZKCY54HcoE51zrFNf6u1NN9ipKrCMK11meVUpWBNUqpQwW9DDNz2fx6FeM0IQXXqy6wVim1T2sdXdpsf1OU118m16gQRWlzGTBXa52llJqI8VvFzWWcqyjMuF5FEYnxGH6qUqo/sBhoYIuGlVK+wAJgitY65e+fvsaXWO162UVh11r3vtHnlVJjgduBXrpggOpvYoGrey5XNtYu01xFPMfZgj8vKKUWYfy6XarCboVcNr9eSqnzSqlqWuu4gl85L1znHFeu13Gl1G8YvR1rF/aivP4rx8QqpVyBAMr+V/5Cc2mtE6768AuM+072oEy+p0rr6oKqtV6hlPpEKRWsy3hvZqWUG0ZRn6O1XniNQ8r0etn9UIxSqh/wNDBAa51+ncN2AA2UUnWUUu4YN7vKbEZFUSmlfJRSflfex7gRfM279zZmxvVaCowteH8s8I/fLJRSFZRSHgXvBwPhwMEyyFKU13913qHA2ut0Kmya62/jsAMwxm/twVJgTMFsj05A8pWhNzMppapeuTeilOqAUfMSbvxVpW5TAV8CUVrrd65zWNleL1veLS7hHeZjGGNRuwversxUqA6s+Ntd5iMYvbvnbZBrMMZP3SzgPLDq77kwZjfsKXg7YC+5TLpeQcCvwNGCPysW/H0YMKPg/S7AvoLrtQ+4twzz/OP1A//G6EAAeAI/Fnz/bQfqlvU1KmKu1wq+l/YA64DGNso1F4gDcgq+v+4FJgITCz6vgI8Lcu/jBjPFbJzrkauu11agiw0ydcUYVtl7Vd3qb8vrJU+eCiGEk7H7oRghhBDFI4VdCCGcjBR2IYRwMlLYhRDCyUhhF0IIJyOFXQghnIwUdiGEcDJS2IUQwsn8PwK+409lisMNAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "x = np.linspace(-2,2,100)\n", "\n", "plt.plot(x, f(x), label='f')\n", "plt.plot(x, dfdx(x), label='df/dx')\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To calculate $\\frac{df}{dx}$ numerically, using finite differences, we can define a function that computes the derivative of a given function $g$ at given position $x_0$, using *step size* $dx$." ] }, { "cell_type": "code", "execution_count": 221, "metadata": {}, "outputs": [], "source": [ "def dfdx_fd(g, x0, dx):\n", " return (g(x0+dx)-g(x0))/dx" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's plot this on top of the analytical solution." ] }, { "cell_type": "code", "execution_count": 222, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3XlYlGX3wPHvPQPDLiqbIgruKzoq7lpmLpVbmi1mprmUlWWLZfb+esvebN/M0nIvs0xNLXM3913Ucd8RFUVAUAFZBmae3x+PkCbMoAyzeX+ui+vVt3lmDgiHe85z7nMLRVGQJEmS3IfG0QFIkiRJtiUTuyRJkpuRiV2SJMnNyMQuSZLkZmRilyRJcjMysUuSJLkZmdglSZLcjEzskiRJbkYmdkmSJDfj4YgXDQ4OVqKiohzx0pIkSS5r9+7dlxRFCbH2OIck9qioKGJjYx3x0pIkSS5LCHGmJI+TpRhJkiQ3IxO7JEmSm5GJXZIkyc04pMZelLy8PBISEsjJyXF0KE7J29ubiIgIPD09HR2KJElOzmkSe0JCAgEBAURFRSGEcHQ4TkVRFFJTU0lISKB69eqODkeSJCfnNKWYnJwcgoKCZFIvghCCoKAg+W5GkqQScZrEDsikboH82kiSVFJOU4qRJEmlKAqp14zEpVwjJSOXtCwjl68ZyTeZQQgEEODtQZC/jmB/LyIq+FKtoi9ajfzlL6lkYr/BN998w+TJk2nWrBlz5sxxdDjSXeLyNSOxZy4TG5/GnrOXOZ6UydXsvNt6Dm9PDbVC/YmuUp62NYNoUzOIYH+vMopYcnYysd9g0qRJLF++XN6glMrcyeRMVh9OYvXhi+w9dwVFAU+tILpKID0aV6ZmiD81QvyoHOhDBT9PKvjq8NSqlVNFUUjPySc1M5dLmUbiU69x/GIGx5Iy+GvfBX7deRaARlXK0bNxOD2ahFOlvI8jP13JzmRiv27EiBHExcXRq1cvhgwZwquvvurokCQ3cyXLyB+GC8zffY6D59MBNfmOur82bWsG0zgiEG9PrdXnEUIQ6ONJoI8nNUKgZfWKhf8t32TmwPmrbD2VyqrDSXy0/CgfLT9Ky+oVGdw2iq4NwvDQOtWtNakMCEVR7P6iMTExyr9nxRw5coT69esDMG7JIQ5fSLfpazYIL8e7PRtafEzBDJvg4GCbvrat3Pg1klzHkcR0pm8+zZ+GCxhNZupXLke/5hE82KgS4WW8kj6bmsWS/ReYu+ss59KyCQ/0ZlDbKAa2icRXJ9d1rkYIsVtRlBhrj5P/spJURradSuW7dSfZfPISPp5aHmsRwRMtqtGoSqDdYqgW5MuL99VixL01+ftIEjO3xPPR8qNM3XSaF++ryZOtquHlYf1dguRanDKxW1tZS5Iz23k6jS9XH2N7XBqhAV68+UBdnmxZjfK+OofFpNUIujasRNeGldh9Jo3PVh5j3JLDTNt0mnd61Kdbw0qypdaNOGVilyRXdDI5g/FLj7DuWArB/l78t0cDnmxVrUR1c3tqHlmRX4e3ZsvJVD5YepgRP+/h3johjOvVkKhgP0eHJ9mATOySVEqXrxn5es1xft5xFl9PLW89WI9BbaLw0TlXQr+REIL2tYP566X2/LjtDF+tPk7XrzfyZre6DGlXHY3siXdpMrHfID4+3tEhSC5EURQW7E7gw2VHuJqdx5OtqvFq5zoEuVD/uIdWw9D21enRuDL/WXSQD5YeYc2RJD5/tAkRFXwdHZ50h2Ril1xXvhFy00HrCR7eoNWBnerEp1IyeXvhAXacTqN5ZAXG92lEvUrl7PLaZSGsnDdTn27O/NgExi05xANfb+LTfo15KLqyo0OT7oBM7JJDGZINxCbF0iCoAWnZaexJ3kP9ivURQnDyykk6RnTEW1GIPfEHMUYz+ssXIOUYZCZDzpWbn0zjARWiIKg2hNSBqA4Q2Q50tlt5mswK0zfH8fmq43h7aPiobzSPx1R1i9KFEILHWlSlTc0gXvp1Ly/M2cPwDtV584F6hZujJNcgE7tkd+vOrmPp6aWk5aQRezEWheL3Usw5Mgeu77XQAo8ZNQSWD6Vd1Sboy9cB7/JgzoP8HMjNgNRTkHoSTv0NWyaA1gsi20Djx6FhH/C8877xuJRMRs/fx56zV+jSIIzxfRoRGuB9x8/nrKpW9GXec20Yv/QwUzedZt+5q3w3oBkhAa5TYrrbycQu2cXfZ/7m9xO/E3c1jvOZ5y0+VqDWrxECFArLKybgVy8FTElMuZpCv9BIQr0VWlVuhz5Uf/OT5GXDma1wai0cWw6Ln4cVY6FJf2jzApSvVuLYFUVh7q5zvL/kMDoPDV8/rqe3Ptyt2wN1HhrG9W5Es8gKjPl9P30mbWHG4BbUCQtwdGhSCTjlzlOpaK72NTIkG1h0YhEHUw9y/PLxW/67Bg0aoUFBQQtgNmFCQasAGi0mQKtR1x4mxQSAWTHfssL3EB68EfMGWaYsYsJibk3yigLxm2H3TDj8p/qLosUw6PA6+FneZXz5mpG3Fu5n5aEk2tcK5ovHmhBWzv1W6ZbsT7jC0B9jyTGamPRUMzrUDnF0SHetku48LXViF0JUBX4CKgFmYIqiKBMsXSMT+51xha+RIdnArou7uHjtIvOPz78lCd+YzD01nrxZ81GuHv2TmIvHwC+E2FrtiWkyGDx9iU2KJSZM/R6OTYolUBfIp7s+Jc+sTj78d5IXCHRaHdO6Trs1uRe4mgDrPwbDHPD0hY5vQesXQHNra2JsfBojf9lL6rVc3uxWj6Ht7942wPNXshk6axcnkjP55JHG9Gse4eiQ7kr2HCmQD7yuKMoeIUQAsFsIsVpRlMM2eG6Hee+99/D396dHjx488cQTCCFYsGABNWvWvOUxo0ePdmCkzsOQbGDIyiGFibfALcm8xZtczThPzOmd6FeOB/9K0Gk8NB+E3uOfOu6Nybngz7Ur1C4yyZsUEwoKuaZcFp5YCFD4i+GmJB8YAb2/hbYvw+p3YNX/weE/oPd3EFIXALNZYcqmOD5beYyICj4seqGdXccAOKMq5X2YP6INz/+8h9Hz95GRk8cz7eQUVGdV6sSuKEoikHj9zxlCiCNAFcClE3uBxYsX07t3b8aNG+foUJyWIdnAunPrWBm/8qak7iE8bk7mxqvEhDZHf3Y3rPkS8rKg/Wtwz2jQlWzHoz5UX2ySN5qMKCgsOrmIP079AQrotDqmdp166wo+pA70nwsHFsDyN+D7DtBlHFejh/La/H38fTSZ7tGV+fiRaAK85QHiAAHenkwfHMOoXw2MW3KYq9l5jLq/tlvfa3BVNr15KoSIApoCO2z5vPYyfvx4fvrpJ6pWrUpISAj169dn0qRJaLVaNm7cyLp16255TPPmzcnPz6dNmzZ89tlndOzYkbFjx6LRaBg/fryjP6UytztpN8NWDiNfyQdAK9SSxk3JvGDVfDUB/ngR4tZD9Xuh+xcQXPuOX7uoJN8kpAnf7/uenRd3AmA0GYlNii26NCMENH4UatwLS0bBirfYuXopu3KGMK5XM55uEymT1r94eWj59smmvLXwAF+vOUGW0cTYB+vJr5OTsVliF0L4A78DryiKcsvMXSHEs8CzANWqWelIWP4WXDxgq9BUlaLhwY+L/c+7d+9m7ty57N27l/z8fJo1a0bz5s0ZMWJEYbmluMd4eHgwa9Ys+vXrxzfffMOKFSvYscMlf7eVmCHZwJoza1gSt6QwqWvQ8EjtR6jsX/nWEsi+32DZG2DOhx5fQ/PBNt1MdGOS99R4MnTlUIxmI2bMHE49zNT9U2lRqUXRCd4/lD/qfcqpI+V4hd/YHpqMb725dtvs5Go8tBo+faQxPp5apmyMQyMEYx6oK5O7E7FJYhdCeKIm9TmKoiws6jGKokwBpoB689QWr2tLmzZtok+fPvj6qptZevXqdVuPadiwIQMHDqRnz55s27YNnc5xk/zKmiHZwDMrnyHfrCb0G0suPWv2vDl5GrPUUsfen6FaG3h4ElSsUabx6UP1TO82nc3nN7Ph3AZWn1nN6jOr8dJ63XJj1WRW+HTlUX7YEEfLqMGkt+tD+WUjYFoXeOp3CC/mJuxdTqMRvN+7IWZF4fsNp9AIeKObTO7OotSJXaj/ktOBI4qifFn6kLC4si5LJfmmtPSYAwcOUL58eZKSkmwZllPJMGYwbtu4wqSuQUPf2n2LXqVfOgnznobkQ3DPG9BxbJHdJ2WhYAXvpfXi2OVjhTdWFxxfUHhTtXpAQ0bN3cv6Yyk81boa7/ZsqO6wrLwGZveBWT3giTlqqUa6hRCC//VuhFmBSetPofPQ8ErnOo4OSwJssU+4HTAQ6CSEMFz/eMgGz2tX99xzD4sWLSI7O5uMjAyWLFlyW49ZuHAhqampbNy4kZdffpkrV67ccr0rMyQb+GjHR/T5ow9xV+LwEB5ohRadVkfPmj0ZFj3s5qR+fBVM6QgZiTDgd+j0f3ZL6jdqUakFXlovBOov5D9O/cHEPRMZunIYPab8zOYTlxjfpxEfPBz9z7b54FowdKXaQTOnn9r7LhVJoxGMf7gR/ZpH8PWaE8zactrRIUnYpitmM+Dy77+aNWvG448/jl6vJzIykg4dOpT4MZcuXeKtt97i77//pmrVqowcOZJRo0bx448/2vvTKBP/bmN8p9U71KlYp+h2QkWBrRNh9X/V+xpP/ALlqzoocnXlPrXrVGKTYjEkG9iQsAEzZowmI/nKUX4e1ovWNYJuvbBcODyzDH55HBY8A4/Nhnout16xC41G8HHfaNKz83hvyWHK++p4uGkVR4d1V5M7T12Ivb9GhmQDOxN3svrMao5ePgqARmh4qelLDIsedusF+Ub46xV180+D3vDw5BK3MdqDIdnA4BVDyDfnIYCn642ggp9X0btVC+Skw0+9Iemg2h5Z6367xuxKcvJMDJ65k9j4y0x9Oob76oU6OiS3U9INSnJkm1QkQ7KBYauGMdEwkaOXj6IRGrX0otEV7ga9SU46/PKomtTvfQv6zXKqpG42K6zY7U366WGEmO/DS+vF7GM/MHHPRIavGo4h2VD0hd7l1JuowXVh7gB1NIFUJG9PLVOfjqFe5QBe/GUPBxKuOjqku5ZM7FKR1p9bT64pF1C36ver3Y+RTUcWvdkn/QLMfFBNeg9PhvvGgsZ5vrVy8kyM/HUP3284Rf/GHVg96Gv61OmDgqKWZcxqr3uxfCvC04vVwWG/9ockt9h7VyYCvD2ZMagFFXx1DPlxF+evZDs6pLuS8/z0SU7BkGzgk52fMP/4fEDtevHSehV9gxQg5bjaGng5Hp6cB/on7R+0BZcyc3liynaWH7zI/3WvzwcPN8JDq6F79e7oNGpLqqIoRAdHW34iv2B15e7pA788BhkX7RC9awot583MZ1qQk2fimZk7uZqdZ/0iyaZkYpcKGZINDF05lJ+P/Ey6MZ1no5/lpWYvFb1KB0jcp67UTbkweKnT1Z9PpWTSd9JWjl5MZ/KA5gzrUKOwXbWg17179e4oKHxn+I4f9v1QfEkG1JvAT/4GWanqTVXjNTt9Jq6nTlgAPzzVnNOXrvHinD3km8yODumuIhO7VGjesXkYzUZAXan7ePoUvUoHOLNN7fP29IFnVjjdRp5d8Wk8Mnkr13Lz+XV4ax5oVOmWx+hD9Xx8z8cMajCIvcl7+dbwreV6O0B4U+g3Ay7uh4XPglkmrOK0rRXM+Iej2XzyEuOXHXF0OHcVmdglDMkGRm8YzV9xfyEQaIQGnbaYm6SgHl4xuw/4h8GQFWrftxNZuj+RAdN2UMFXx8IX2tK0WgWLjy/vXb7wz7mmXMv1doC6D0LX8XD0L9j8hS1CdluPtajKkHbVmbklnt92nXV0OHcNeYLSXe7G8QAaNLzZ4k2yTdnFtwAeXwW/PQXBdWDgIvB3rkMXpm2KY/yyIzSrVoFpT8dQwc/6aIeYsBi8tF7kmnJRUAjUlWBEb+vn4cIeWDsewps5XRnKmbz9UD1OpmTyf4sPUiPEnxZRFR0dktuTK/ZivPfee3z++eccPXoUvV5P06ZNOXXqVJGPuV2bNm2iYcOG6PV6srNv7hoYPHgwCxYsKFXst2PK/imF4wGEEGSbsosvvxxdCnOfhND6MOhPp0rqZrPC+0sO88HSI3RrUIk5w1qVKKmDWpKZ1nUazzV+jjDfML7e8zVf7f7KcklGCOg5Qf1a/D4ULp+x0Wfifjy0Gib2b0rVCr48//MektJzHB2S23PpxG5INjDtwDTLP4ClVDCPfe/evTcdslEac+bMYfTo0RgMBnx87vxw5dIwJBsYvmo4m85vKuxR99R4Fl9+ObJEnftSuQk8/YfaAugkCtoZZ2w5zeC2UXw3oBnenrc3vkAfqmdk05G82vxV0o3pzDg4g2Grhln+3tL5weM/g9mkfm3yc0v5mbivQB9PfhjYnCxjPs//vBtjvrw3UZacshTzyc5POJp21OJjMo2ZhcOdBIK6Ferir/Mv9vH1KtZjTMsxFp/THvPYp02bxrx581i5ciVr1qzh559/5qWXXmLt2rVUr14de+wEVndgDsakmNAIDW+3fJuMvIziyy9H/oL5g9Ubh08tVDftOIkrWUae/Wk3O+PT+M9D9RnWoXqpJgwmXktEIFBQLM9yLxBUU+3d/20ArP0fdP3gjl/b3dUOC+Czfk148Zc9fLD0MO/3buTokNyWUyb2ksjIyyg871JBISMvw2Jit8Ze89iHDRvG5s2b6dGjB/369WPhwoUcO3aMAwcOkJSURIMGDRgyZMgdfx7WKIrCN3u+KTwcWiDIyMsoekQAqOWX+YOgsl7t43aipH7+SjaDZuzkbGoW3/RvSq8m4aV+zoJ6e44pp3AUsVX1e0DMEHVGTs37oeZ9pY7DXXVvXBnDuepM3XQafdXy9G0mz04tC06Z2K2trOGfUkKeOQ9PjScfd/jY8srKCkfNY9+4cSP9+/dHq9USHh5Op06d7vhzsMaQbGDCngnEJsXedNJRseWXYytg3iC1/DJwIXg7z7mfhy5c5ZmZu8jOM/HjkJa0qVnEIK87UDA0bNuFbSyJW8IP+34gLSeN+6reZ/n7q+t4iN8Ci0bA81vUDU1SkcY8UI/9CVd5e9EBGlUJpE5YgKNDcjsuW2Mv+AEsdpv7HXDUPHZ7HE5QUH4pSOpvt3rb8tfu5BqYNxAqNbpefnGepL7pRAqP/7AdrUbw+/NtbZbUC+hD9Tyvf56RTUeSkZdRwnq7L/SbDtlp8MdIdcqlVKSCm6n+Xh68MGcPWcZ8R4fkdlw2sYP6A1hsB8dtctQ89nvuuYe5c+diMplITExk3bp1pf5c/k1RFCbsmVBYfgFIN6YX/7U7vVEdeBVSV03qPuVvfYyD/L47gWdm7iKigg+LXmhXpqu9C5kXCue4F9TbLaoUDZ3HwfHl6jA0qVih5byZ8ERTTqVk8s7iQ44Ox+04ZSnGERw1j71Pnz6sXbuW6Oho6tSpw7332va0HkOygYl7J5a8/HJmm7pdvkJ1GOg83S+KojBx7Um+XH2ctjWD+H5gc8p5l6D+XQr/rreH+5Wght9qhLpxacVYqHEfBMq55MVpVyuYlzvVZsLfJ2hVoyKPxThubr+7kfPYXcjtfo0MyQaeWfEM+Up+Yfkl3ZhefPfL+d3wY28ICINnloO/c8zTzjOZeWfxQebuOkffZlX4uG9jdB72ebNZcDjH3KNzCfEJoXuN7rSq3Mryu8S0OJjcDiLbwoAF8lBsC0xmhYHTd7D37BWWvNSOWqGy3m6JnMcuMckwiXzln/qlxfLLxYMwu+/1EbV/Ok1Sz8zNZ9iPsczddY6XOtXii0eb2C2pg1ruG9VsFEMaDeF0+mm+M3xnfZ5MxRrQ+T31PoUsyVik1Qi+flyPj07LS78ayMkzWb9Iskom9jLUp08f9Hr9TR8rV64s89c1JBt4fvXzbEvcVrLNRynH1VOCdH7qjlInKR8kXs3m0e+3sfnkJT7uG83rXeva5UZzUQpe98b+dotaDIfIdmpJJj3RDhG6rtBy3nzWrzFHEtP5ZIXl/StSyThVjV1RFIf94JaFRYsW2ey5Sloyu2n2S0k2H6Wdhp96gdCoK/UKUTaLuTSOJKbzzMxdZObmM2NwC+6t49jxBf+eJ9MstJnlCzQa6DURJrWBFWPgsZ/sE6iLur9+GIPbRjFzSzwdagfTqV6Yo0NyaU6zYvf29iY1NdUuOy9djaIopKam4u3tbfWx0w9M/2f2yw2bj4pM6lcT1KSen6OOCXCSKY3rjibTb/JWAOY918bhSR3+mSfTLaobCgqHU0twilJQTbj3DTj8h7onQLLorQfrUa9SAKPn7yclQ45nKA2brNiFEDOAHkCyoih3tE84IiKChIQEUlJSbBGS2/H29iYiovhdeoZkAzMPzWR9wno0QoNAWC6/ZCar5ZfsK2r5JaxBGUV+e37aFs97fx6ifuVyTB/UgkqB1n+Z2Ys+VE+TkCZk5WXx1Z6vSMpK4v5q91u+kdp2FOyfD8tGQ/UOTnUOrLPx9tQysX9TekzczJjf9zN9UIxbvYO3J5t0xQgh7gEygZ9KktiL6oqR7ty/yy9jW44lMy+z+PLLtVT4sYd6nN3ARVCttd1j/rd8k5kPlh5h1tZ4OtcPZcITTfHzcqpKYaG/z/zNK+tfAcBb6219g9yZrepJU21fkrNkSmDmltOMW3KYDx5uxFOtIx0djlOxa1eMoigbgTRbPJd0+2Yfnn1T+SUzL7P48kv2Ffi5D6Segv6/OkVST8/JY8iPsczaGs/Q9tX5YWCM0yZ1gNPppws3LpXoYI7IttDsadg2CZLkZhxrBrWJokPtYMYvPUJcSqajw3FJTlNjl26fIdnAGxveYNWZVWgoQfdLbgbM6QdJh+GJOVCjoz3DLdKZ1Gv0nbSVrdc7X97p0QCtxrnffseExaDTXj8IG4XoICsHYYO6I9W7HCx7U44bsEKjEXz+aBO8PDW8+puBPHle6m2zW2IXQjwrhIgVQsTKOnrpGZINDFk5hBXxKxAIxrQYY3n2i/GauqP0/B54dCbU7mL/oP9l84lL9Pp2C5cyc5k9tBVPtKzm6JBKpOBG6qN1HgVgy4Ut1i/yrQid3oEzm+HQwjKO0PWFlfPmoz7R7Eu4yqR1p6xfIN3Ebu93FUWZAkwBtcZur9d1V78c/YU8cx4AGqEhy5RV/OjdvGz4tT+c3QZ9p0L9nnaM9FaKojBjSzzjlx6mVqg/U5+OITLItW4q6kP16EP15Jvz+enQT+SZ8+gW1c1yrb35YNg9C1a9A3UekDdSrXgwujIP68OZuPYEneqFEh3hPIPonJ0sxbigFfErWHl6JQJhvfySl6MO9Dq9UT0QIrqffYP9l5w8E6Pn7+d/fx2mS4MwFr7QzuWS+o06R3bGhImfj/xsfUeqRgsPfQbp52GTPAS7JMb1akSQv47X58tdqbfDJoldCPErsA2oK4RIEEIMtcXzSrf6ft/3vLnhTWpVqMWULlMsl1/yc9Uj2079rW6WafKE/QO+wbm0LB6ZvJXf9yTwSufaTB7QHH8nvklaEscvH7+9G6nVWkPjx9VDOdLi7BChawv09eTjRxpzPCmTr1Yfd3Q4LsNWXTH9FUWprCiKp6IoEYqiTLfF80o3m2SYxHeG71BQOJN+Bm8P7+K7X/Jz4beBcGIl9PgKmg20f8A32HA8hZ7fbuZsWhbTB8XwSuc6aJz8JmlJ3HgjFaBJSBPrF3UeBxpPWP1uGUbmPu6rG0r/ltWYsimO2HjZfFcSshTjIhafXMzkfZML/55vzi9+dXhjUu/+pXpsm4OYzApfrj7O4Jk7qVTOmyUj23N/fffZLl5wI7Vvrb4oKCU7WL1cZWg3Co78qY5Jlqz6T/f6hAf68MaC/bIkUwIysbuAL3d/yTtb3qFGuRp4ab0s19XzcuC3p/5J6i0cVxVLycjl6Rk7+ObvE/RtGsHCF9oSFey69fTi6EP1jGs3jk5VO/H9/u/5evfX1hN825EQEA4r3wazbOezxt/Lg0/7Neb0pWt8vvKYo8NxejKxO7mPdnzEzIMzATh/7bzltkZjFvz6OJxYpZZfHJjUN5+4xEPfbCI2/jKf9mvMF481wVfn2vV0a7pX747RZGT6wenWb6Tq/OD+d+DCHjj4u/2CdGHtagUzoFU1pm85LUsyVsjE7qQMyQaGrhzKL0d/Kfz/8s35XDVeLbquXrD5qKD7xUHllzyTmY+XH2XgjB0E+njyx8h2d83JOGczz97eUXqNn1APCl/zntqSKlk19iFZkikJmdid0N7kvQxeMZidF3eiQYNOo7NcfslKg9l94Ox2eGQa6J+0f9BAXEom/SZv5fsNp+jfshpLRranXqVyDonFEf59I7XYFtQCGg10HQ/pCbDj+zKOzj3cWJKRXTLFc+/3xi7IaDLywfYPCg+eFkLwcK2HqexfueihXumJ8HNfSD2pzvyu38PuMSuKwuztZ/hw2RG8PbVMHtCMB6Mr2z0ORyu4kTrtwDQ2JGwgM68Ec06qd4DaXWHzV+oGJp8KZR6nq2tXK5j+LasydVMcD0ZXRl/VeQ5bdxZyxe5Etl3YRu/FvTl++TgewqNwld6zZs+iyy+pp2BGV7hyFgbMd0hSP38lm6dn7OS/fxyidY0gVr5yz12Z1AvoQ/V81fErqgZU5YvYLzCZS1AuuP9dyElXk7tUImMfqk9ogDdvLthHbr4syfybTOxOYt25dTy3+jkSMhPw1Hjydqu3LW8+Or8HZjwAuZnqPPUaHe0ar9msrtK7frmB3Wcu88HDjZg5uAVh5ZxnfrqjeGo9GdVsFCevnGT0htHWO2QqNVI3Le34Aa6et0+QLq6ctycf9m3E8aRMvpOzZG7hWondlKfWk93M4dTDjN00FgV1hI5ZMRd/kxTg+CqY1R08vGHICqjS3K7xnkzOpP/U7byz+CBNq1Vg5Sv38FTrSHkowg3CfMMQCNacXcOwVcOsJ/f73gbFDOs/sk+AbqBTvTD6Nq3CpHUnOXwh3dHhOBXXSuwr34YpHSH5iKMjsQlDsoG3Nr3FwOUD8dZ6W79JCuoQqV+fgODaMGw1hNS1W7w5eSY+X3mMByds5EhiOp8+0pjZQ1tStaKqQPFSAAAgAElEQVSv3WJwFTd2xJSoQ6ZCJLQYBoY5kCL7tEvqnR4NKO/ryVsL95Mvx/sWcq2bp42fUM+PnNbFaUbP3qm9yXsZsmII+Uo+AsG4e8cR6BVIbFJs0TdJTfmw+h3YPkm92dZvJnj52yVWRVFYfTiJD5Ye4WxaFn2bVuHt7vUJ9veyy+u7ooLDr3NMOSgoNAgqwdGDHUbDntmwbrw8/LqEKvjpeK9XQ0b+spcZW07z7D01HR2SU3CtFXtEcxi+DipWh18eg23fudyhBYZkA9/t/Y4xG8eQr6inHmmEhhNXTqAP1Rddfsm+rPaob58ErV+AJ361W1I/djGDp6bv4NnZu/Hy0PDL8FZ8+bheJnUr9KF6pnadSv96/QGIvViCoyD9gqDNC+riJXFfGUfoPrpHV6Zz/TC+XH2cM6nXHB2OU7DJmae3q9RnnhqvwaLn4MgSaNQPek6wW6IrjYJNR0azEQCt0ALgqfEs/iZp0iF17suVs3Yd5pV4NZsJa04wL/YcAd6evNalDgNaVcND61prAWfwxoY32JCwgWV9lxHsE2z5wTlX4evGULUVDJhnnwDdwMWrOXT5cgPREYHMGdbKbe/3lPTMU9cqxRTQ+cGjP8GWr2DtB5B0EB6bDSF1HB1ZkQzJBnZd3MWm85sKk7oGDY/UfqT4/nSAvT/D0tfBuzwM/ssu55Nevmbk+w2nmLU1HrOiMKhtFC93qk0FP531i6UijWw6klXxqxi1bhRvxLxh+TAO70B1QNjf4+DcTqja0n6BurBKgd6Mfag+by86wLzYczzewjVO4yorrrliv1HcelgwFPJz4MFP1V2XTvTb2pBsYNiqYeSacgH1sGkhBDqNrvhVem4mLHsD9v0C1e9Vd5P6h5ZpnJcyc5m26TSzt8WTlWeij74Kr3apI2+M2oAh2cDgFYMxKSa8tF5M6zrNcnI3XoMJTSC0PgxaYr9AXZzZrPDE1O0cTUxnzev3Ehrgfq237r1iv1GNjvDcRrU088cL1+ePf62eMelABav0nRd33pTU+9XpR7h/ePGr9LPb1c/l8hm4d4z6odGWWZxnU7OYseU0v+06R06+iR6Nw3mpUy3qhAWU2WvebWKTYilYQBUcxmExsev8oP1rsHIsxG2AGvfaKVLXptEIPuobzYMTNjHuz8N8N6CZo0NyGJdK7NsvbGdv8l7ahLe5+QcjsAo8/Yd6Ks3aD9S3sN2/hHoPOSTOIlfpCHRaHb1q9ir6hzovBzZ8DFsmQGBVeGY5RLYpk/gURSH2zGVmbjnNioMX0QhBL304L3SsRa1Q579X4WoKZsjkmnJRUKjqX4KhaDFDYOs3sP5jqH6PU70LdWY1Q/wZdX9tPlt5jN6HLtK1YSVHh+QQLlOKufHtrKfGk7davkW6Mf3Wle8FAyx+AZIPQf1e6hmTAfb5xzUkG9ieuJ2NCRs5cOkAoCb1vrX7EhEQUfwq/dRatZaeFgdNB8IDH4GX7VfM6Tl5LN57njnbz3IsKYNy3h482SqSwW2jqBTofm9bnYkh2cCGhA38dOgnukV148MOH1q/aMcUWP6Gumip0bGsQ3QbeSYzPSdu5nKWkdWv3Us5b09Hh2QzJS3FuExin3ZgGhP3TMTMP5sQNGjQaYuoVZvyrq92PgEPL7hnNLR8DjzLLnntTd7L0JVDyTPnATev0outpV85B2veVedxV6wJ3T+Hmp1sGleeycymEyks3HOe1YeTyM03E10lkKdaV6Nnk3C3n5HubL6I/YKfDv/Eot6LqBFYw/KD83Lgm6bq5qVnlstV+20wnLtCn0lbGNCqGh88HO3ocGzG7RK7IdnA8FXDyTPnoaBgVv5J8B0jOtIktMmtK+LUU7B8DJxcDYHVoPO70LCvOi7VRvYm72Xu0blsPr+ZdKO6rVmDhn51+hXf8ZKVpp5Sv3Oq+vcOr6udEDb6xZObb2LryVSWH0xk9eEkLmflUcHXk15NwnmkeQSNI+Q0PEdJy0mjy4IuRAZE8t82/7Vcawf1e2TZaLlqvwPjlhxi1tZ4FoxoQ/NIx95zsxW3S+ygJvfYpFgCdYF8susTjCZj4XwVUPvB32zxJpl5mTcn1FPr1F2bFw9AcF1o/wpEPwraO3uLZkg2sOX8FuKvxrPyzMrCGKz2pWdcVAc97ZoGxkxo8iR0fAvKl/4ginNpWWw8kcL6YylsPXmJa0YTAV4e3F8/lO6Nw7m3Tgg6D9mD7mg3lhR1Gh3Tu023nNwLVu3lq6lzgeSqvcSu5ebT9auN+Oq0/PVye7w8yq4JwV7s2hUjhHgAmABogWmKonxsi+f9N32ovvCHoHaF2sQmxXIu/RyLTi5CQSHPnMf4HeMB0Gl0jGk55p86/LMb4dBCdTTq4udh3YfQfBDoB0C5cKuvbUg2sCNxB5dzLjP32NzCeekFtEJbdF+6osD53RA7A/bPA3M+1O+pDn0KrX9HX4d8k5kTyZnsO3eFnafT2HE6jfNX1BN4qpT34eGmVbi/fijtagW7xTezO7mxQ8ZoNlrvkPH0hg6vqav2uHU2L9W5Mz8vDz54uBHPzNrF9+vjGNW5tqNDsptSr9iFEFrgONAFSAB2Af0VRTlc3DW27GO3VKIp4Knx5NN7PiXYJ5jYi7uIyVPgwHxiLx8lJtcIVWKIDY0ipl4/8KlQOK8FYNWZVVzOucyy08tueW6BQCu0KCg3r9IVRR1UdmQJ7P8N0k6Bpy80fQpaPw8VrdRWr1MUhUuZRuJTr3HsYgZHL6ZzNDGDQxfSyb5+LFiQn45WNSrSqnoQ7WoFUTPE32133bmDgu/Xgg6Zj9p/RI+aVubo5+f+s2qXtfbbNvKXPaw6lMSyUR1cvuvLbqUYIUQb4D1FUbpd//tYAEVRip0/atMNStxcovl016eFuzuLSvKgJmSN0GBWzGhQQFEwow7OURCYLfzc/DuZv9niTa5mJRPjWQF9Tq7aanlqLWReVC+I6qDO2m7QS91ViJqwM3LzSc/O40pWHqnXjFzKyOVSZi6JV3M4fyWbC1eyOZOaRWZufuFrB3h5ULdSAI2qBKKvWp7GEYFUD/aTidzFGJINbErYxI+Hf+S+qvfx2b2fWb9oxw+w/E0Y9Jd66pJUYikZuXT+cgN1wwKY+2xrNBrX/XmxZ2LvBzygKMqw638fCLRSFGVkcdfYOrHf6N91+DxTHlqNlmp+dTmVcbDkT1TwdRECoShoECiotauBGYFkYyQmJ4dWOVcpZ7pSeFmGphyHvJqyT9eMPbqmXFSCyM03k5tvJsuYT1auiWvGfMzFfNkDvDwIL+9D5fLeRFb0JSrYj6hgP+qEBRAe6C2TuBuZsGcC0w9MZ1HvRdQsb2UqYV62OkMmtJ7cjXoH5u06x5u/7+ejvtH0b+m64wbsWWMvKtPckraEEM8CzwJUq2abL2xOnonTl64Rl3KN+NRrJFzO5vwVI4lXapOSmUu6MgQP3zjys2pwBfCtdgxEPijXbyIKMyia6+9s1T+rfzEh0CJQUBQzHoqGp674kK0x0ihH0CA3C6PQcUVUZrO2HimeIcRrqxPvEcUlbRgeHho8NAIPjYbyHhp0Hhq8PDT46Tzw9dLip/Mg0MeTcj4eBProCPbXEezvRZC/jgA36rmVLHu6wdPMOTKHj3d8TKvwVsXvcwDw9IF2L8Oq/4OzO6BaK/sG6+IejYlg4d4EPlx2hPvrhRLq5id9uVQpZveZNLacTOVIYjpHEtM5k5Z109Tein46qpT3Iby8N6EB3oQEeBHs70VFP08CfXQkG49xIt1Aq8ot8fLQsO/SblpUagFwU129qD9bbUuTpDswdtNY/or7q/g9GTcyXoOvoyG8KTz1u30DdQNxKZk8MGETXeqHuey4AXuu2HcBtYUQ1YHzwBPAkzZ43lss2ZfIrK3xRAX5Ur9yOR5uWoWaIf7UCPEjKsgPPy9rn07b6x+qluH//OPe+MNU3J8lydYq+6kHf5sxk2fOs9wlo/ODNi/C3++rnVZ2PhLR1dUI8eflTrX4fNVx+hxOonODMEeHVGZs0scuhHgI+Bq13XGGoijjLT3+TlfsqZm5eHlq8beawCXJNdzU167VMb2rlb72nHR11R7ZDvr/Yr9A3YQxXx03kJ6Tx6pX73G50mdJV+w22bGiKMoyRVHqKIpS01pSL40gfy+Z1CW3og/V8/V9X6MVWtpUbmP9HaJ3OWg1Ao4thaRiO4qlYug8NHz0SDQX03P4fKX7ni0rtyJKkoN1rNqRx+o+xpbzW0jMTLR+QavnwNNP3Wwn3bZm1SrwdOtIftp+ht1nLjs6nDIhE7skOYEhjYaAgE92fcK0A9MwJBuKf7BvRYh5Bg4uUCeCSrftjQfqUamcN2MX7seYX/R+F1cmE7skOYFKfpVoH96ev8/+zcQ9Exm+arjl5N72JdB4qPP7pdvm7+XB/3o34nhSJj9sOOXocGxOJnZJchKRgZHAzR0yxQqopI6oMPwC6RfsFKF76dwgjO6NKzNx7UlOJmc6OhybkoldkpxE52qd0Qj1R9JD41G4l6JY7UaB2QRbv7VDdO7pvZ4N8dFpeXvhAczFbQd3QTKxS5KT0Ifq+bC9erLSA1EPWO+QqRAF0f1g9yx1xr9020ICvPhP9/rsjE/j111nHR2OzcjELklOpHuN7txf7X7WnltLprEE5YH2r0LeNXVImHRHHm0eQduaQXy87CgXr+Y4OhybkIldkpzM8OjhZBgz+O3Yb9YfHFof6j4EO3+AXPeqE9uLEIIP+0RjNJl554+DOOLwIVuTiV2SnEzD4Ia0DW/L9IPTmWyYbLk7BtRVe/Zl2POjfQJ0Q1HBfrzWpQ6rDyex7MBFR4dTajKxS5ITuq/qfWQYM5i8b7L11seqLSGyvXoTNT/XfkG6maHtqxNdJZB3/zzI5WtGR4dTKjKxS5ITyjBmABQe+Wix9RGgw6uQcUE9flG6Ix5aDZ880pgrWXn8b6lrj2uQiV2SnFCLSi3w1KgDqjRCY731seb9UKkxbPlabYGU7kiD8HI837EmC/ecZ/2xZEeHc8dkYpckJ6QP1TOj2wyCvYMJ8QmhcUhjyxcIAe1fgdSTcHSpfYJ0UyM71aJWqD9vLzxARk6eo8O5IzKxS5KT0ofqeaPFG1y4doF159ZZv6B+b6hQXV21u0Fnh6N4eWj55JHGJKbn8MmKo44O547IxC5JTqxrVFci/CP4Zvc3TN0/1fJNVK2HOkPm/G6I32S/IN1Q88gKDGlXnZ+3n2XbqVRHh3PbZGKXJCfmofGgS2QX4tLjmLi3BMPB9APALxQ2f22/IN3U6K51iQzyZczv+8ky5js6nNsiE7skOTlfT1+ghB0ynt7QegSc+hsS99kpQvfko1NLMmfTsvjMxQ7lkIldkpxc68qt8RDqyWFaobXeIRMzFHQBctVuA61rBPF0m0hmbY1n52nXmccjE7skOTl9qJ5JnSeh0+poEtrE+nAwn/LqQRyHF0PaafsE6cbGPFCPiAo+vLlgn8uUZGRilyQX0Ca8DYMaDCL2Yiynr5YgWbd+QT2IY5sc6Vtafl4efPpIE+JTs/h0hWuUZGRilyQXMaD+AHRaHTMPzrT+4HKVofHjsPdnyEwp++DcXJuaQQxuG8WsrfFsj3P+LplSJXYhxKNCiENCCLMQwkrhT5Kk0gjyCaJv7b78eepPvtr9lfXhYO1GqbNjdsqRvrbw5gNql8wbC/aRmevcJZnSrtgPAn2BjTaIRZIkK1pUaoFJMTHj4AzrrY/BtaFed9g5VY70tQFfnQefP9qEhMvZjF96xNHhWFSqxK4oyhFFUVyj6CRJbuBM+pnCPxvNRuvDwdq9AjlXYM9PZRzZ3aFFVEWe7VCDX3eeZd1R550lI2vskuRCYsJi0Gl0AGgowXCwqi0gsh1s+w5Mrjn3xNm82qUOdcMCGPP7fqcd72s1sQsh1gghDhbx0ft2XkgI8awQIlYIEZuSIm/mSNKd0Ifqmd5tOtXLVcfbw5s6FepYv6jdK5CeAAcWlH2AdwFvTy1fPt6Ey1lG/s9JT1yymtgVRemsKEqjIj7+uJ0XUhRliqIoMYqixISEhNx5xJJ0l9OH6nm/3ftk5mWy6OQi6xfU7gKhDWHLBDCbyz7Au0DD8EBe6VyHpfsTWWw47+hwbiFLMZLkgvShepqFNmPK/in8sO8HyzdRhVA7ZFKOwIlV9gvSzY24tyYxkRX47+JDnEvLcnQ4Nyltu2MfIUQC0AZYKoRYaZuwJEmyplO1TqTlpPGd4TvrHTKN+kJgVXWkr2QTWo3gq8f1KMDr8/ZhMjtPSaa0XTGLFEWJUBTFS1GUMEVRutkqMEmSLDOa1Bt3JRoOpvWENiPh7DY4u8NOEbq/qhV9eb93Q3bGp/H9hlOODqeQLMVIkou67ePzmg0Enwpy1W5jfZpWoUfjyny1+jh7z152dDiATOyS5LL0oXqmdp1KOV05qgZUpUlIE8sX6Pyg5XNwbBkku+bJQM5ICMH4PtGElfPm5bl7SXeC4/RkYpckF9Y8rDmjmo0i7mocOy/utH5By2fBwwe2flP2wd1FAn08+aZ/Uy5cyeE/ixzfAikTuyS5uN61ehPsE8y0A9OsP9gvCJoPgv2/wdWEsg/uLtI8sgKvdanDkn0XmB/r2K+tTOyS5OK8tF4MbDCQ7YnbGb99vPXhYG1eVA+73vadfQK8i4y4tyZtawbx7p+HOJGU4bA4ZGKXJDdQv2J9AOYem2u99bF8NYh+FHb/CFmucyqQK9BqBF8/rsfPS8vzc/Y47GAOmdglyQ0cSj2EQABqG6T14WCjIO+aOvlRsqnQct5MeKIpp1IyHVZvl4ldktxATFgMOq06HAyB9dbHsAZQ5wHY8T0Yr5V9gHeZdrWCeeX+Oizae565u87Z/fVlYpckN6AP1TOt6zT0IXpQIMS3BPOY2r8K2WmwZ3bZB3gXeqlTLTrUDubdPw9x8PxVu762TOyS5Cb0oXo+u/czNBoNsw7Osn5BtdZQrS1snQj5zjl+1pVprtfbg/10PDd7N2l2HPErE7skuZFKfpXoXbM3C44vYMKeCdY7ZDq8dn2k73z7BHiXCfL3YvJTzUnJzOWlX/eQb7LPdE2Z2CXJzbSu3Jp8JZ/pB6Zb75Cp1RkqRatjBuRI3zLRpGp5PujdiC0nU/lspX0OnJOJXZLcTEKmujlGQbF+fJ4Qaq390nE4+pedIrz7PNaiKgNaVeOHjXGsPpxU5q/nUeavIEmSXRUcn2c0G0t2fF793lChOmz+Cur3VJO9ZHPv9mxIRT8drWtULPPXkit2SXIzBcfn1QisgU6ro3aF2pYv0HpA+1fgwh6IW2+XGO9GOg8Nr3etS4C3Z5m/lkzskuSG9KF6xrcfT1Z+Fr8d+836BU36Q0A4bPqi7IOTypxM7JLkphoFN6JteFt+PPQjOfk5lh/s4QVtX4L4TfIgDjcgE7skubHh0cNJy0lj9IbR1lsfmw8C3yDY9Ll9gpPKjEzskuTGPDQeCAQbEjYwbNUwy8ld5wetn1cPvE7cZ78gJZuTiV2S3NiNrY4lGg7WYjh4lZO1dhcnE7skubGYsBi8tF6A2teuD9FbvsCnPLQcDof/hBT7bKaRbK9UiV0I8ZkQ4qgQYr8QYpEQorytApMkqfQKzkXtXbM3APHp8dYvav0CePrARllrd1WlXbGvBhopitIYOA6MLX1IkiTZkj5Uz//a/Y9GQY2YdmAaeWYrhy37BUOLoXBwAaSesk+Qkk2VKrErirJKUZSCI0K2AxGlD0mSJFsTQjCiyQjOZ57nzQ1vluD4vJdAq5O1dhdlyxr7EGC5DZ9PkiQbCvQKRCBYc3aN9Q6ZgDBo/gzsmwuX4+0Wo2QbVhO7EGKNEOJgER+9b3jMf4B8YI6F53lWCBErhIhNSUmxTfSSJJXYbXfItBsFGg/Y9GUZRybZmtUhYIqidLb034UQg4AewP2KhcP9FEWZAkwBiImJsf8hgJJ0lyvokMkxqbtQ9aFWOmTKVYZmA9VDr+8ZrR6CLbmE0nbFPACMAXopipJlm5AkSSoLBR0yvWr2QkEh/mq89Yvav6r+r1y1u5TS1ti/BQKA1UIIgxDiexvEJElSGdGH6vmg3QdEB0czZf8U8kxWOmQCI6DZ07B3Nlw+Y58gpVIrbVdMLUVRqiqKor/+McJWgUmSVDaEELyof5HEa4m8vuH1Ehyf9zoIjZwh40LkzlNJugv5efohEKw7t856h0xgFWg+GAy/yA4ZFyETuyTdhW67Q6b9ayC0sPGzMo5MsgWZ2CXpLvTvGTLRQdGWLyhXGWKeAcOvcjeqC5CJXZLuQgUdMo/VeQyAfZdKMKa3/avqbtQNn5RxdFJpycQuSXcpfaied9q8Q8eqHZl1cBZXc69aviCgkjr5cf88SD5inyClOyITuyTd5UbqR5KRl8Er616x3iHT/lXQ+cPaD+wTnHRHZGKXpLtcdn42GqEhNinWeoeMb0X1bNSjf8H53fYLUrotMrFL0l0uNikWrg/5yDXlWu+QafOCejaqXLXftrVn1zLtwDTr74xKyeqsGEmS3FtMWAw6rY5cUy4KCuH+4ZYv8ApQ2x9X/QdOb4LqHewTqIv7/fjvvLftPQQCL60XU7tOtT6v5w7JFbsk3eUKOmSGRQ/DW+vNytMrrV/UYiiUqwJr3oXiZ/9J1ymKwpT9U9Q/o5BnzrP+zqgUZGKXJAl9qJ6Xm73MsOhhrD23lve3vW+5XODpA/e9rdbZDy+2X6Auau25tVy4dgEPjQdaocVT40lMWEyZvZ5M7JIkFSo47Hr+8fkMXzXccnJv0h9CG8Df70O+0U4Rup48cx5f7f6KGoE1mNZ1GiObjizTMgzIxC5J0g0OpB5AIIAS3EjVaKHze5AWB7tn2SM8l2NINvDautc4k36G15q/RvOw5gyLHlamSR1kYpck6QYFN1ILNAlpYvmC2l0hsr26GzUnvYyjcy2GZAPDVg1jfcJ6BIJAr0C7vbZM7JIkFdKH6pnWdRp9avVBQWF/yn7LFwgBXd6HrEuwZYJ9gnQRsUmxGE1qiUogyvRm6b/JxC5J0k30oXreb/c+HSM6MvXAVC5lX7J8QURzaNQPtn0LV87aJ0gXEOEfgXJ9g4BOqyvTm6X/JhO7JElFej3mdXLzcxm3dZz1TTVdxgECVr9rt/ic3V9xf+Gt9WZY9LAyv1n6bzKxS5JUpKjAKDpHdmZ9wnom7plouUsmMALavQyHFsLZ7fYN1AltPr+ZDQkbeEH/AqOajbJrUgeZ2CVJsiCyXCQAZszWN9W0GwUB4bB8DJjNdorQ+cQmxfL2prcJ8w1jQP0BDolBJnZJkorVvkp7PDTq5BEhhOU6sc5PbX9MNMC+X+0Sn7MxJBsYvnI4l3Mvk5aTxuHUww6JQyZ2SZKKpQ/VM6PbDCr7VcZH60P1wOqWL4h+FCJaqKMGsq/YJ0gnsu7cOvKVfADMitmunTA3KlViF0L8TwixXwhhEEKsEkJYmR4kSZKraRralG86fcO1/Gu8u/VdyzdSNRp46HPISoV14+0bqBMoaA/VoCnzsQGWlHa642eKorwDIIR4GfgvMKLUUUmS5FTqVaxH52qdWXVmFevOrkOn1RXf6RGuh5ihsGsa6Aeof78LbEzYSGxSLI/VeYzK/pWJCYux+03TAqVasSuKcuNWMz8KpzpLkuRuagTWAEp4I7XT/6kz25e+flfcSN2RuIOxm8YS7hfOWy3fssvYAEtKXWMXQowXQpwDBqCu2CVJckPtqrTDU+MJqDspLZYZfMpDl//B+VjYO9tOETqGIdnAc6ufI92YTkp2CodSDzk6JOuJXQixRghxsIiP3gCKovxHUZSqwBxgpIXneVYIESuEiE1JSbHdZyBJkl0U3EitGVgTIQQVvCtYvqDJExDZDla/AxlJ9gnSAZbGLcWkmADH3jC9kVBsNCRfCBEJLFUUpZG1x8bExCixsY7/5CVJun3JWck8vPhhIgIi6BLZhRaVWhRfdrh0Aia3gzrd4HH3W7nnmfLoubgn5zPPoxEadBoL9x5sQAixW1EUq3dkS3XzVAhRW1GUE9f/2gs4WprnkyTJ+YX6hvJo3UeZcXAGR9OOWj7mLbg2dByjzmw//Cc06GX/gMuIIdnAt3u/5XzmeV5r/homxWT9hmlmMviHlnlspa2xf3y9LLMf6AqMskFMkiQ5OX9Pf0A95s1oMlouP7R9GSpFw7LRkH3ZThGWLUOygaErh7Lj4g40QkPT0KbWb5ie2wlfNYJjy8s8vtJ2xTyiKEojRVEaK4rSU1GU87YKTJIk59WiUgu8tF6AmtybhTYr/sFaT+j1LVy7BCv/Y6cIy9b2xO0Yzf+cGmW1rm68BotGgH+Yet+hjMmdp5Ik3baCue1dI7uioLA0bqnljUvhemj/ChjmwJG/7BtsGTh15RSgbkTSaUowknfNe5B2Ch7+DrzLlXl8pd2gJEnSXUofqqdJSBOGrhrKvOPz1CRnaePSvW/BidWw5GWo2tIuteaysPXCVlbEr6Bztc40DG5ova4etx52ToFWz0P1e+wSo1yxS5J0x4QQNA1pCpRg45KHDvpOVcsSf74ENurIs6dNCZt4bd1rhPuF81GHj6zX1XOuwuIXIag2dLbfrHqZ2CVJKpUOER0KNy4pKJbLEqH1oPM4OL7C5Q7A3pu8l5FrR3It/xqXsi9xNM1KE6CiqL/AMhKhzw/g6WOfQJGJXZKkUirYuNS6cmvMipn159Zbrre3fBZq3AcrxkKS43dpltRkw2TMijoewaSYrN8wjZ0Oh/+A+/+rHh9oRzKxS5JUavpQPd93/p6GQQ2ZfnC65ROXNBroOwW8A+G3gZCTfutjnMy2C9vYnrgdjdCgFVrrkxsT98OKt6FWF7Xd085kYpckySa0Gi1tw9sCar3daLbQ3+4fCv1mwOV4p6+3rz27lpfXvky4XzhTOk9hZNORlneX5mbA/GmtWAcAAAvdSURBVMHgWxH6fK/+IrMzmdglSbKZeyLuQafRAaAoCvoQCzcWo9rB/e/A4cWw4wc7RXh7dlzYwSvrXiHHlENKdgpeHl6Wb5iazbDwObh8Gh6ZDn7B9g34OpnYJUmyGX2onundptMtqhsKCj8e+pGp+6cWX29vOwrqPAir/qO2BToRk9nEhzs/RLk+jbxEdfV14+HYUuj2ofqLy0FkH7skSTalD9WjD9XjpfHiz7g/2ZCwofh5MgX19uldYd7TMHQNhNRxTOA3MCQb+CL2C+KuxuGh8UBRFOt19QMLYNPn0HQgtHLseUMysUuSVCaiAqMAtQUy15RLbFJs0SUM73Lw5G8w7X745TEYvlatTzuIIdnAMyufId+cj1Zoebvl21w1XrW8Een8bvjjRajWBrp/CULYN+h/kaUYSZLKxL/nySRmJhbfBlkhEp74BdIvwNwBkJdt52j/8eOhH8k35xf+/arxquW6espxmPOoekP4sdnqRiwHkyt2SZLKRME8ma0XtvLnyT+Zd3weAlF8WaZqS+gzGRYMVdsgn5gDHl52i9eQbGD24dmsObsGDRqEENbLL1fOweyHQWhg4GLwD7FbvJbIxC5JUpkpqLebFTM/7P/hpjG/Ra6AGz0CuZnqPJnfh0K/WaAt+zRlSDYwZOUQ8sx5aNAwpuUYsvKzLJdfrl2C2X3UeAf/BUE1yzzOkpKlGEmSylz7Ku0LyzJmzJxNP1t8Wab5IHjgYziyBBY/D6b8Wx9jY9MPTCfPnAeo82+y8rMsl1/SE2FWd7h6Tr0/ULlxmcd4O+SKXZKkMldQltlyfgvLTy9n0clFlssyrZ+HvCz15CVjptoTrvO1eVyGZAMT9kwgNikWjdAgKEH55XI8/NRbXbEPWACRbWweV2nJxC5Jkl0UlGU0QsOkfZNQUMgx5bA9cXvRK+MOr4NXOVj2hlry6P+rTbtldlzYwbNrnsWsmNEKLWNbjiUjL8Ny+SX5iBpLXjY8/afdZ8CUlEzskiTZVZvwNkw/OB2jyYiCwqLji8jOz+a+qvfdmlBbDge/EFg4HGY+CI/PgeBapXp9Q7KBFfH/3965xkhVnnH8999lV1IvgMLiiqKQIAg1pctKFZtWK1WgLVihST8Vi2BJQ9ombS2RpB9sbAPUfvCWpsUamxiFQtVtQV0rGGgMyIoLK/d7ilxWpFmgtLK78/TDObMdhp3ZmZ09Z4bp80sm886cd973f54588xz3usbNOxr6FrUC+BM+xnm3jo38wd3NARNQ9WXw3fXwNBxBemIElkR1mior6+3pqYeZnA5jlO2NLc203SiiWNnj7FizwoA+lX04/n7nu8+Wj64HlbMhs52mP4kfPaBXtX7QesHzHlzTtdwxtTJRxnXf0l0BjNKNzwBwyYEQxoHDOtV/YUi6X0z62G7JnfsjuMUkWUty3hqy1MkCCLnms/UMHXEVCYPn3yxk207EiyudWQz3DYPvvpYzu3uza3NrD64msZDjZz6zykAKlTBrFGzqL2iNnPzy6kD0PADOLQhmFE67ddQ1b+QUy4Id+yO45Q8za3NzGuc1zUipdM6AahUJQvGLwBxodPtOB/sH7rxGRgwHKYuhjHTstbReLCRRzY8ckHZQPYovbMDNj4L634ZbMZ93+OBYy/yjNJYHbuknwBLgSFmdrKn/O7YHcdJkmyWOXr2KKv2rOqK3pNUVVSx8LaFnG4//T8nf/hdWP1jaN0BN0+Bux+F2s91lVVXU8euU7tYvns5B9oOdJXVY5Se6Aw2x1i/NCh79DT42hNw1XVxmKJHYnPskm4AlgFjgAnu2B3H6Q2ZovdU+lX0Y/bY2VRXVjNp6ETY8Sqbt7/M6HNnaK0dw+NqoyPtcxVUUFlRScISmaP0f/8zGDf/7tNwcjcMvhnuXgRjZxQ9Sk8lTse+EvgF8BpQ747dcZzekoy4B1QPYMnmJbQn2jHsgtEr+VKpSmaOmnlhlH7+X8HQxePbYPcbsH8tJNphyC3w5Z/C2PuhorIPz6xvyNWxFzTcUdJ04CMz26oS+ldzHOfSJDnWHWDUoFEXOXmAhCW61khPIsSk2jvYfPw9Oq2DSjMw6BRUJdr5RtNyxlt14Lw72+HMcUiWcdX18IXvwbgHYFhdSUXovaXHiF3S34Bruzm0CHgUuNfM2iQdIkvELulh4GGA4cOHTzh8+HAhuh3H+T+iu0g+2QnaaZ1dTSwATSeaqK+pg5N7afro79Qnqhh/7iwkOqCiKlh7ZuCNUDMWho6FQSMuGWceeVOMpFuBt4Fz4VvXA0eBiWZ2PNtnvSnGcZzeknTyyWn/yXTG2aJlRORNMWbWAtSkVHiIHNvYHcdxektqc03ytXMhvrqj4zhOmdFna8WY2U19VZbjOI7TezxidxzHKTPcsTuO45QZ7tgdx3HKDHfsjuM4ZYY7dsdxnDKjKMv2SvoY6O3U08FAKY6Vd1354bryw3XlR6nqgsK03WhmQ3rKVBTHXgiSmnKZeRU3ris/XFd+uK78KFVdEI82b4pxHMcpM9yxO47jlBmXomP/XbEFZMB15Yfryg/XlR+lqgti0HbJtbE7juM42bkUI3bHcRwnCyXv2CUtlbRL0jZJr0gamCHfFEm7Je2TtDAGXd+StF1SQlLGHm5JhyS1SGqWFPki9HnoitteV0t6S9Le8HlQhnydoa2aJTVEqCfr+Uu6TNLy8PgmSTdFpSVPXQ9K+jjFRnNj0vUHSa2SPsxwXJKeDHVvk1RXIrruktSWYq+fx6DpBknrJO0Mf4s/7CZPtPYys5J+APcC/cL0YmBxN3kqgf3ASKAa2AqMjVjXLcBo4B2Cdegz5TsEDI7RXj3qKpK9lgALw/TC7r7H8NjZGGzU4/kD3wd+G6a/DSwvEV0PAk/HdT2l1PsloA74MMPxacDrgIDbgU0lousu4K8x26oWqAvTVwJ7uvkeI7VXyUfsZtZoZh3hy40EOzWlMxHYZ2YHzOw88DIwI2JdO81sd5R19IYcdcVur7D8F8L0C8D9EdeXjVzOP1XvSuAeRb+xbzG+l5wws/XAqSxZZgB/tICNwEBJtSWgK3bM7JiZbQnTZ4CdwLC0bJHaq+QdexpzCP7l0hkG/CPl9REuNmSxMKBR0vvhvq+lQDHsNdTMjkFw4ZOy+1Ya/SU1SdooKSrnn8v5d+UJA4s24JqI9OSjC2BmePu+UtINEWvKlVL+Dd4haauk1yWNi7PisAnv88CmtEOR2qvPNtoohGwbZpvZa2GeRUAH8GJ3RXTzXsHDfXLRlQN3mtlRSTXAW5J2hVFGMXXFbq88ihke2msksFZSi5ntL1RbGrmcfyQ26oFc6vwL8JKZfSppPsFdxVci1pULxbBXLmwhmIZ/VtI04FVgVBwVS7oCWAX8yMxOpx/u5iN9Zq+ScOxmNjnbcUmzga8D91jYQJXGESA1cklurB2prhzLOBo+t0p6heB2uyDH3ge6YreXpBOSas3sWHjL2ZqhjKS9Dkh6hyDa6WvHnsv5J/MckdQPGED0t/w96jKzT1Je/p6g36kUiOSaKpRUh2pmayQ9K2mwRbw3s6QqAqf+opn9uZsskdqr5JtiJE0BfgZMN7NzGbJtBkZJGiGpmqCzK7IRFbki6XJJVybTBB3B3fbex0wx7NUAzA7Ts4GL7iwkDZJ0WZgeDNwJ7IhASy7nn6p3FrA2Q1ARq660dtjpBO23pUAD8J1wtMftQFuy6a2YSLo22TciaSKBz/sk+6cKrlPAc8BOM/tNhmzR2ivO3uJe9jDvI2iLag4fyZEK1wFr0nqZ9xBEd4ti0PVNgn/dT4ETwJvpughGN2wNH9tLRVeR7HUN8DawN3y+Ony/HlgWpicBLaG9WoCHItRz0fkDjxEEEAD9gT+F1997wMiobZSjrl+F19JWYB0wJiZdLwHHgPbw+noImA/MD48LeCbU3UKWkWIx61qQYq+NwKQYNH2RoFllW4rfmhanvXzmqeM4TplR8k0xjuM4Tn64Y3ccxykz3LE7juOUGe7YHcdxygx37I7jOGWGO3bHcZwywx274zhOmeGO3XEcp8z4LyTywDbM4KtpAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.plot(x, f(x), label='f')\n", "plt.plot(x, dfdx(x), label='dfdx')\n", "plt.plot(x, dfdx_fd(f, x, 0.3), '.-', label='dfdx_fd')\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Do\n", "\n", " - Vary the third parameter `dx` of the last plot. Change it from 0.3 to larger and smaller values. What happens? \n", " - Modify the definition of `dfdx_fd` to a) backward difference, b) central difference form. How does the derivative approximation change?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Calculating the error" ] }, { "cell_type": "code", "execution_count": 251, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xd8VFX6x/HPSa+EhCQQCBAgoYQOoYMICCgqoKBgoSgqKNbVtayruz9X1+4qIigogoqIoICKoAjSa0IvobfQAoGQXuf8/rgBKQmkzMydmTzv1ysvQ2bm3m+uyZMz9577HKW1RgghhOtwMzuAEEII65LCLoQQLkYKuxBCuBgp7EII4WKksAshhIuRwi6EEC5GCrsQQrgYKexCCOFipLALIYSL8TBjp6GhoToqKsqMXQshhNNKSEg4o7UOu97zTCnsUVFRxMfHm7FrIYRwWkqpw6V5npyKEUIIFyOFXQghXIwUdiGEcDGmnGMvTn5+PklJSeTk5Jgdxen4+PgQGRmJp6en2VGEEA7AYQp7UlISgYGBREVFoZQyO47T0FqTkpJCUlIS9erVMzuOEMIBOMypmJycHKpVqyZFvYyUUlSrVk3e6QghLnKYwg5IUS8nOW5CiEs5zKkYIYSVWCxwehec3g05qZCdany9ah0IjoJqDcA32NSIwraksF9i3LhxTJw4kTZt2jB9+nSz4whReoX5sGMu7J4PB1dA1pmSn6vcoG4XiB0AjW+DKhH2yynsQgr7JSZMmMCCBQvkIqRwHtmpkDAV1n0G6cchMAKib4L63SGipTEy96kK2gKpR4yPY/Gw8yf49TlY+CK0HgY3/B2Capn93QgrkcJeZMyYMRw4cID+/fvz4IMP8swzz5gdSYiSaQ1bvoPf/gHZZ6HeDXD7hxDdG9xKuHRWPdb4aHQz9Pyncapm/WTjD8Pmb6H9w3DjS+AdYNdvRVif0lrbfadxcXH6yl4xu3btokmTJgD838872Hk8zar7jK1ZhX/d3vSaz7nQwyY0NNSq+7aHS4+fcHFnD8Avz8CBpRDZHvq9AzVbl397qUdg2duwaTqE1IdBk6FWW6vFFdajlErQWsdd73kONStGCHEdifPh025wbCPc+j48+FvFijoYF1UHfAIj50NBLnzRB1a8b1yEFU7JIU/FXG9kLYSrO5eZx5akVI6ey+ZEajanzmfRI3kat6VM5bBPY36IfpPA7Cgid5yiUY1A6oX6V3zaa1QXeHSl8W5g8WuQvAsGTAAPL+t8U8JuHLKwC1HZ5OQXsnLvGX7feZL4Q+c4cCbz4mN+bgWM95lIT8safvPowXtuj3JiWz4ZubsuPic0wJsO9ULo3jCMm5vXoIpPOdtL+AbD4C+hRgtY/H+QeQaGfA3egRX9FoUdSWEXwiRaaxIOn2P6uiMs2nmKjNwCqvh40L5eNQbHRdK6djD1g90J//Uh1L410Od1+nZ6nL5KobUmLbuAo+ey2Jp0nvUHU1h38Czzt53glXnbuSm2OkPb1aZrdGjZR/JKQbe/QUB1+OkJmHobDJsDfiG2ORDC6qSwX+LQoUNmRxCVQEGhhZ+3HufLVYfYmnSeQB8P+jWvQb/mEXRuEIqXR9Glr/xs+O5e2P8n3P4RtB15cRtKKYL8PAnyC6JZrSDu7VAHrTWbj6YyZ9Mxft5ynPlbT9AyMojHe8ZwU5Pwshf41veBfyjMHAbT74Lh82TGjJNwyFkxouzk+Dk+rTULtp/kvd92c+BMJg3C/BnZpR6D2tTCz+uKMVZBLnx7NxxYZlzYbH1fmfaVW1DIjxuPMWHpPo6ezaZ5rSBeG9CU1nXKccdp4nyYeT/UvxHumSnn3E0ks2KEcCAbj5xj4CereGz6RtzdFJ8Na8uiZ7ozrGPdq4u6xQJzHzWmMw6cWOaiDuDt4c497euw5NkbeXdwC5LTc7hz4mpe+nEbqVl5ZdtY41vh9nGwfwnMHSOzZZyAnIoRwobScvJ5d+Fuvll3mOqBPrw7uAV3tonE3e0ap0WWvAbbf4Cb/g2t7qnQ/j3d3bgrrjY3N6vBh3/sZerqQyzaeZL37mrJjY3CS7+hNsMgKwX++JfRb6bXqxXKJWxLCrsQNrIk8RQv/rCNMxm5jOwcxbN9GhHgfZ1fufgvYeX/oO0D0OVpq2UJ9PHkldtiGdQmkmdmbmbklxt45Ib6PNen0V/n9K+n69Nwdr8xx71WW2MkLxySnIoRwspy8gt5dd52HpwaT4i/F3PHduFftze9flE/tBLmPwsxfaDfe8bsFCuLrVmFeY934b4OdZi0/AB3fbaGU2ll6OV/y7vGDVFzxsCZfVbPJ6xDCrsQVrTnVDr9x6/kqzWHGdW1HvMe70KLyKrXf2H6SZj1gHFL/+Ap4G67N9M+nu68cUdzJtzXhn2n0hn4yarSt/Dw9IG7vwY3D5h5H+Rm2CynKD8p7FamtWbp0qUsXbqU4mYcJSYm0qlTJ7y9vXnvvfcue2zhwoU0atSI6Oho3nrrrYtfP3jwIB06dCAmJoYhQ4aQl1fGi1/CLhZsO8HAT1ZxNjOPqQ+045XbYvH2cL/+CwsLYPaDkJdh15uB+jWPYNaYzgDc9elq/kxMLt0Lq9Y2/vic2QMLX7BhQlFeUtjLqLCwsMR/Z2dnM3LkSLZv38727dsZOXIk2dnZlz0/JCSEcePG8dxzz121nbFjx7JgwQJ27tzJjBkz2LlzJwAvvPACzzzzDHv37iU4OJgvvvjCRt+dKI9Ci+adhYk8On0jjWoE8ssT3cp2YXLJf+DwKrjtfxBu3ymrsTWrMHdsF+qF+TNq2gbmbEoq3Qsb9DCuAWz6BhJ/tW1IUWZS2C/xzTff0L59e1q1asXo0aMvFu2AgABeffVVOnTowJo1a4iKiuK1116ja9euzJo16+LrfX19mThxIl9++SVffvklEydOxNfX97J9hIeH065dOzw9L7/le/369URHR1O/fn28vLwYOnQo8+bNQ2vNkiVLGDx4MAAjRoxg7ty5Nj4SorSy8goY/XU8E5bu5572tfnukY7UCPIp/Qb2/QGrPjQulrYcarug11C9ig/fj+5Ex/rV+Nv3W/h+w9HSvfDGl6B6c/j5SaP1gHAYjjkrZsGLcHKbdbdZoznc8laJD+/atYuZM2eyatUqPD09eeyxx5g+fTrDhw8nMzOTZs2a8dprr118vo+PDytXrrxsG9nZ2YwdO5YHHngAgLFjxzJhwoSrintxjh07Ru3atS/+OzIyknXr1pGSkkLVqlXx8PC4+PVjx46V6VsXtnE6PZdR0zaw/dh5XhvQlOGdosq2gayzMHcshDWGm0v+2bQHPy8Ppoxsx+ivE3j+h63kFloY1rHutV/k4QV3ToJJ3eHnp2DINza54CvKzjELuwkWL15MQkIC7dq1A4wiHR5uvJ12d3dn0KBBlz1/yJAhV23D19eXKVOmsGzZMsAo7KW9jbu48/GqqCdIcV8X5tp/OoMRU9aTkpHHpGFx3BRbvewbmf+ssYTdfd8bFyVN5uPpzqThbRk7fSOvzN2Ol7tiSLs6135R9Vjo+QosegW2zjTtXYe4nGMW9muMrG1Fa82IESN48803r3rMx8cHd/fLL4L5+/sXux2lFDfeeGOZ9x8ZGcnRo3+9BU5KSqJmzZqEhoaSmppKQUEBHh4eF78uzLP92HmGT1mPm4LvHulIy9qlmPVypW2zYcePxkpGES2tH7KcvD3cmXBfWx76Kp6XftxGVT8v+jatce0XdRoLu342VnOK6SPNwhyAnGMv0qtXL2bPnk1ysjEz4OzZsxw+fNhu+2/Xrh179+7l4MGD5OXl8d1339G/f3+UUvTo0YPZs2cDMG3aNAYMGGC3XOJyGw6d5Z5Ja/H1dGf2mM7lK+ppJ4zRemQ76OJ4SzB6ebjx6f1taBFZlSdmbGLN/pRrv8DN3bjwm51q3JkqTCeFvUhsbCyvv/46ffr0oUWLFvTu3ZsTJ05YfT8nT54kMjKSDz74gNdff53IyEjS0tLw8PBg/Pjx9O3blyZNmnD33XfTtKmx4Mjbb7/NBx98QHR0NCkpKYwaNcrqucT1Ld9zmmFfrCMs0JtZYzoRFVr8u7brWvB3KMiBgZ/adL56Rfh5efDlyHbUCfHj4a/i2XXiOvPcazQzRu4bv4Ija+0TUpRIuju6CDl+trV8z2ke+iqeBmEBfD2qPaEB3uXbUOKv8N09Rq+Vbs9aN6QNnDifzcBPVuHh5sa8x7tc+/vOy4RPOhjz8EcvB/dyLvYhSiTdHYWwkhV7T/PwV/FEhwXw7UMdyl/Uc9Ph1+cgPBY6P2ndkDYSEeTL5OFxpGTm8ug3CeQWFJb8ZC9/uOUdSN4JayfYL6S4ihR2Ia5h5d4zPDQtnvphAUx/qAPB/hXoRb7kdUg7biya4USj2RaRVXnvrpZsOHSOf87ZXuxMrYsa94OYvrD8PZnbbiKHKuxmnBZyBXLcbCP+0Fke+moD9UL9K17Uj22EdZ9Bu1FQu731QtrJbS1q8mSvGGYlJPHVmutMKujzunFaZunVM8yEfThMYffx8SElJUWKVBlprUlJScHHx/x50K5k+7HzPPDlBmoG+fLNQx0IqUhRt1hgwfPgH+bUfcyf7mUssff6/J1sPppa8hPDGhp/wOK/hORE+wUUF1X44qlSqjbwFVADsACTtNYfXes1xV08zc/PJykpiZycMrQQFYDxRzEyMvKqNgWifPYlZ3D3Z2vw9XRn1phO1Kx6/TuHr2nzDGPloYETodW91glpktSsPG4dZ9xxPf/JrlT1K+EPXmYKjGsNdTrAfbOKf44os9JePLVGYY8AIrTWG5VSgUACMFBrvbOk1xRX2IVwBMdTsxk0cTX5hZpZYzpRr7xTGi/ITYeP20JQbRi1CNwc5k1yuW0+mspdn67mhpgwJg+Pw62k1aBWjTPuSB02Bxr0tG9IF2W3WTFa6xNa641Fn6cDu4BaFd2uEPZ2PiufEVPWk5FTwFcPtq94UQdY/i5knDJmi7hAUQdoVbsq/7w1lsWJyXyx8mDJT+wwGqrWhUX/knVS7cyqP2lKqSigNbCumMceUUrFK6XiT58+bc3dClFhOfmFPPTVBg6nZPHZ8LbE1qxS8Y2e2QdrJkCr+yGybcW350CGd6pLn9jqvPvb7pIX6fDwhh7/gJNbYddP9g1YyVmtsCulAoAfgKe11lf9n9ZaT9Jax2mt48LCwqy1WyEqrNCieXLGJuIPn+N/Q1rRuUGodTb8x7/Aw8epL5iWRCnFW4NaEOTnydMzN5GTX8L89uZ3QWgj+PO/YLnGHHhhVVYp7EopT4yiPl1r/aM1timEvfznl538vvMUr94Wy60tIqyz0cNrIPEX6PoUBJaj86MTCPH34t3BLdhzKoN3f9td/JPc3I1R+5ndsE0uotpLhQu7MnrIfgHs0lp/UPFIQtjPlJUHmbr6EKO61uOBLvWss1GtjYuGgRHQcax1tumgbmwUzohOdfli5UFW7i3hhqQm/aFGC2Nee2G+fQNWUtYYsXcBhgE9lVKbiz76WWG7QtjUbztO8p/5O7m5aQ1e7mfFPjs750HSBujxMnj5WW+7Duqlfk2oH+bPCz9sJSO34OonuLkZ7YnPHYJNX9s9X2VkjVkxK7XWSmvdQmvdquhDFkEUDm1rUipPfbeJlpFV+d+QViVP2SurgjxY/H8Q1sTp56yXlo+nO+8ObsHx89m8s7CEG5Ji+kBke1jxgYza7cA15l8JUQYnzmfz0LR4qvl78/mIOHy93K//otJKmApnD0Dv14zzy5VE27ohjOwcxVdrDrPuQDH925WC7s/D+aPGSkvCpqSwi0olK6+Ah6bFk5VXyJSR7crfqbE4eZnGvPW6XSCmt/W26yT+3rcRtUN8eeGHrWTnFTMDJvomY7WoFR/IDBkbk8IuKg2LRfP0d5vZdSKNj+9tTaMagdbdwfpJkJlsTG+shOvS+nl58PadLTiUksWHf+y5+glKGT3oz+6HHXPsH7ASkcIuKo33F+3m952neOW2WHo0CrfuxnPOw8oPIbo31Olo3W07kc7RoQyJq83nKw+SeLKYG5ca327Ma1/xvtyNakNS2EWlMG/zMT75cz/3tK/NyM5R1t/Bmk8gJ9WY/VHJvXhLY6r4ePDynO1YLFf0onJzM0btyTthzwJzAlYCUtiFy9t8NJW/z95K+3oh/F//ZihrnybJTDEKe+wAqNnKutt2QsH+XvyjXxMSDp9jVsLRq5/QbBAERxmjdmnTbRNS2IVLO5WWwyNfxRMe6M3E+9rg5WGDH/lVH0J+ljFvXQAwuG0k7euF8OaCRFIyci9/0N0DOj8BxxLgyBpzAro4KezCZeXkFzL66wQycguYPDyOatacAXNBxmnY8LnREyWskfW376SUUrwxsBkZOQW8taCYue0t7wXfEFj9sf3DVQJS2IVL0lrzz7nb2Xw0lffvakmTCCt0ayzOmo+hIAdu+Ltttu/EYqoHMqprPWYlJF294pKXH7R/GHb/CqeLmUEjKkQKu3BJU1cfYnZCEk/2iuGW5lZq7HWlzDOwfrJxzjg0xjb7cHKP94wmLNCbf/+04+oLqe0eNrpfrhlvTjgXJoVduJzV+87w+vxd9ImtztO9bFhwV38M+dlww/O224eTC/Tx5IWbG7P5aCpzNh27/MGAMGh5D2z5DjKSzQnooqSwC5dy9GwWY7/dSP1Qfz6wZg+YK2Wm/DVaD2tom324iDtb16Jl7aq8tTDx6iZhncZCYZ5xLIXVSGEXLiM7r5BHvk6gwKKZNDyOAG8P2+1szXhjJkx3Ga1fj5ub4t+3x3I6PZePl+y9/MHQGGjUz7gAnZ9tTkAXJIVduAStNc//sJXEk2mMu6e1ddYrLUn2OWOE2XSgzIQppdZ1grmzTS2+XHmIo2ezLn+w46OQfRa2zTYnnAuSwi5cwqTlB/h5y3Ge69PI+u0CrrR+MuSlG3dQilL7e99GuLlx9WpLUV0hvCms+0xuWLISKezC6a3Ye5q3FybSr3kNHruxgW13lpsBaydAw5uhRnPb7svFRAT58nC3+vy05fjl0x+Vgg6j4dQ2OLzKvIAuRAq7cGpHUrJ4/NtNxIQH8u7gltZvF3ClhKnGqZhuz9l2Py5qdPcGhAZ488b8nehLR+ct7gbfYFj3qXnhXIgUduG0svIKeOTreLTWTBreFn9bXiwFyM8xpjjWuwFqt7PtvlxUgLcHf+vdkA2HzvHbjlN/PeDpC21HQuJ8SD1iWj5XIYVdOCWtNc/P3sruU+mMu6c1davZ8GLpBZunQ8ZJObdeQXfHRdKwegBvL0wkv/CS1r1xowBlzJARFSKFXTilz5Yf4JetJ3i+b2NutPXFUoDCAlj1EdSKg3rdbb8/F+bh7sbzfRtz8Ewms+KT/nqgam1ochskTJOpjxUkhV04neV7TvPOwkRubRHBmO717bPTnXMh9TB0+1ulXB3J2no1CSeubjAf/rHn8mX02j1s9LXf/qN54VyAFHbhVA6nZPLEjE00rB7Iu4Nb2P5iKRhT8FZ+aKz80/AW2++vElBK8cItjUlOz2Xq6kN/PRDV1TjO8V+Yls0VSGEXTiMzt4BHvkoAYNKwOPy8bHyx9IL9i42peF2eMlYAElbRLiqEno3Dmbh0H+ez8o0vKgXtRhm92o9vMjegE5OfUuEUtNY8N2sLe5PT+eTeNtSp5me/na/8EKrUMnquC6v6e99GpOcWMGHZvr++2HIoePrBBhm1l5cUduEUPvlzHwu2n+Qf/ZrQNSbUfjtOiodDK4xmVR5e9ttvJdEkogoDWtZk2upDJKfnGF/0CYLmg40WA9mp196AKJYUduHwFu86xfuL9jCwVU1Gda1n352v/B/4VIU2I+y730rk6Zsakl+omfDn/r++GDcKCrJhywzzgjkxKezCoe09lc5T322mac0qvDXIThdLLzizz7hhpv3D4B1gv/1WMlGh/gxuE8m3645wPLVommPNVsbU0g1fSP+YcpDCLhzW+ax8Hv4qHh9PdyYNi8PH092+AdZ8DO5e0P4R++63EnqiVzQazfg/LznX3m4UpOyV/jHlIIVdOKSCQguPz9jIsdRsPr2/DTWr+to3QEYybJ4Bre6BADvcAFXJRQb7MbRdHb7fcPSvtr6xA8E7yLhhSZSJFHbhkP77ayIr9p7h9YHNiIsKsX+A9ZOMlX06PWH/fVdSY3tE4+am+Ghx0WIcXn5Gc7Cd8yDrrLnhnIwUduFwvlt/hCmrDjKycxRD2tWxf4C8TKNfSeNbITTa/vuvpGoE+XB/h7rM2XSMQ2cyjS+2HQGFubB1prnhnIxzFfbCAkg7bnYKYUPrDqTwyrztdIsJ5Z+3NjEnxKZvjNa8nZ80Z/+V2Jgb6+Phpv46116juXERNWGqXEQtA+cq7L8+B1/0gXOHzU4ibODo2Swenb6R2iF+jL+3DR7uJvx4WgphzScQ2R7qdLD//iu58EAf7isatR9OuWTUfjoRjq43N5wTscpvjlJqilIqWSm13RrbK1HcA5CbBtNug9SjNt2VsK+0nHxGTdtAoUXzxYh2BPl6mhNk189Gs68uMlo3y5juRaP2JUWj9qZ3glegMWoXpWKtIdFU4GYrbatkES1h2FzIPg/TbpfTMi6ioNDC2OkbOXA6k4n3tbHtQtTXorWxkEZwPWjUz5wMgvAqPtzTvg4/bjrGkZQs4x6C5oNhxxzIOW92PKdglcKutV4O2Oeyda02MOxHyDxjFPfMFLvsVtiG1pp//7yDFXvP8MYdzegcbcd2AVc6ug6OxRvtA9zsPGdeXObRGxvg7qb45MK59rYjjDtRt802N5iTcK5z7BdExsF9s4zTMTOGQF6W2YlEOU1ZdYhv1h5hdPf65syAudTqj411N1vda24OQfUqPtzbvg4/bEwy5rVHtILqzYwL2+K67FbYlVKPKKXilVLxp0+frvgG63aCQZONJk0/Pmxc9BJOZeH2k7w+fyd9m1bnhb6NzQ2Tst9oHxA3CrxMOhUkLjO6e32Ugs+W7zfa+ba+H45vhFM7zI7m8OxW2LXWk7TWcVrruLCwMOtsNHYA3PwWJP4CC1+yzjaFXWw6co6nvttEy8iqfDikNW5uJq9KtHYCuHtK+wAHEhHky+C2tfl+QxKn0nKgxRCjxYOM2q/LOU/FXKrjGOj0OKz/DOKnmJ1GlMLhlEwemhZP9So+fD4iDl8vk89nZ52FTdONuxwDq5ubRVzm0e4NKNSaScsPgF+IcVF7y3dQkGd2NIdmremOM4A1QCOlVJJSapQ1tltqvV+D6N7w6/NweI1ddy3KJiUjl5FfbqBQa6Y+0I7QAG+zIxkDgoJsY4AgHEqdan4MaFWT6esOk5KRC62HQfZZ2P2r2dEcmrVmxdyjtY7QWntqrSO11vZd+sTNHQZ9DlXrwPfD4fwxu+5elE5mbgEPTt3A8dRsPh8eR/0wB2iFW5Br9IVp0AvCTbrTVVzTYzdGk1tgYcqqg9Cgh7GalZyOuSbnPxVzgW9VGPot5GfBzPshP8fsROIS+YUWHp2+kW3HzjP+3jbmNPYqzvYfIOMUdJbRuqOKDg+gX7MIpq0+zPlcizFraf9iGcBdg+sUdoDwxnDHZ8aV89//aXYaUcRi0bwweyvL95zmv3c0p3esg5zH1hpWj4fwplC/h9lpxDU81qMBGbkFfL3mkFHYtQW2fmd2LIflWoUdoMltxrnSDZNhx1yz01R6Wmte+2UnP246xrO9GzK0vclz1S91YCkk7zBuSLLnykyizJrWDOLGRmFMWXWI7IC6ULcLbP5WGoOVwPUKO0Cvf0GttvDTE3D2gNlpKrX//bGXqasP8VDXejze08Fa4K4ZDwHVjdvVhcN77MZozmbmMXPDEWPUnrIPkjaYHcshuWZh9/CCwV8ao7BZDxgXyITdfb7iAOMW72VIXG1evrWJfdcrvZ7kRNj3B7R7GDwcYGaOuK729UKIqxvM5BUHyW90O3j6webpZsdySK5Z2AGC68KACXBiM/z5htlpKp2v1x7m9fm76Ne8Bv+9s7ljFXUwRusevhD3oNlJRBk81qMBx1Kzmbcr3bhBcfuPkJ9tdiyH47qFHYzz7W0fgFXj4OAKs9NUGjM3HOGVudu5qUk4Hw5pjbvZd5VeKSMZtn5vrGfqX83sNKIMejQKp3GNQD5dth9Li3uMNt6J882O5XBcu7AD9H0DQurDnDGQnWp2Gpf3Q0ISL/64je4Nw/jkvjZ4eTjgj9iGL4zl1jqONTuJKCOlFI/1iGZfcgaLsmMgqI6cjimGA/7WWZmXv9EsLP2EsQKTsJlZ8Ud5bvYWujQI5bNhbfH2cMDWt/nZxoyphrfIeqZOql+zGtQO8eXT5QfRLYfC/j/hfJLZsRyK6xd2MGbI3PgibJtl3JAirO7bdUf4++ytdI0OZfLwOHw8HbCog7EoclaK3JDkxDzc3XikW302HUlla+itgDZOrYmLKkdhB+j6N6PAz3/OOMcqrGba6kP8Y842ejYOZ/JwB2jqVRKLxVjPNKKlMQ9aOK274mpTzd+LjzbmQ51ORmMwmdN+UeUp7O4exiyZvAz45Rn5IbACrTXjl+zlXz/toE9sdT69v63jjtQB9i2CM3uMG9gcbZaOKBMfT3dGdI5iSWIyJ6MGwpndcHyT2bEcRuUp7GC0HOjxstG/XU7JVIjFovnPL7t47/c93NG6luNeKL3U6o+NBlJN7zA7ibCC4Z3q4uflzriTTcHd2xi1C6CyFXaAzk9ArTjjQmr6KbPTOKX8QgvPzdrClFUHeaBLFO/f1RJPdwf/UTqxBQ6tgA6jjQU1hNOr6ufF0HZ1+H57Oln1+xrX0KRPO1AZC7ubOwycaKyTuuB5s9M4nbScfB6cuuFi75dXb4s1f/Wj0lg9HrwCoM0Is5MIKxrVrR4amKe7G33a9y0yO5JDqHyFHSCsIXR/HnbOlZsbyuB4ajZ3f7qGNftTeGdwC57oFeN4d5QW5/wx2PEjtBlutHcWLqNWVV9ubxHBf/dEYPELgy0zzI7kECpnYQfo8pSx6vn8ZyHnvNlpHN7WpFTumLCKY+eymfpAe+6Oq212pNJb96nR5rXDGLOTCBt45IYGpOfB1uA+sHuhsdRhJVd5C7v22xjIAAAde0lEQVS7J/T/2FhkYdG/zE7j0OZsSuKuT9fg4ebGrEc70TUm1OxIpZebDglTjb4iwXXNTiNsILZmFbrFhPLuqdZgyTfenVVylbewA9RqAx0fg4Qv4dBKs9M4nIJCC//9dRfPzNxCq9pV+enxLjSuUcXsWGWTMM3oJ9LpCbOTCBsafUMDVmVEkBoYDVtmmh3HdJW7sAP0+AdUrWvMbZf2vhclp+Vw3+frmLT8ACM61eWbhzpQzREWni6LwnxYO9G4GSmyrdlphA11ia5GbEQQ3+d2gaT1kLLf7EimksLu5Q+3vm/cuLLqI7PTOISVe8/Qb9wKtiad5/27WvJ/A5o5/nTG4uyYA2lJ0PlJs5MIG1NKMbp7faakxaFRlb7FgBP+ttpATG/jppXl71Xqv/R5BRbeWZjIsCnrCPbz4qfHuzCobaTZscpHa1g9DkIbQkwfs9MIO7i1eQTuVSPZ5tXS6AlUie8ul8J+wc1vGSvpVNJ2A7tPpjPwk1VMWLqfu9vWZt7jXYipHmh2rPI7sBRObjNuSHOTH/PKwMPdjQe71mNaRkc4dxCOrjc7kmnkJ/6CwBpw07/g4LJK9TYuv9DChKX7uP3jlSSn5zB5eBxvD26Bn5eH2dEqZvXH4B8Oze82O4mwoyHtarPSqxN5yhu2Vt4WA1LYL9X2QaPdwO8vV4pFObYcTaX/+FW8s3A3PRuHs/DpG+gdW93sWBV3cjvsX2y0D/D0MTuNsKMAbw/u6NCYBQVtKdz2Y6WdECGF/VJubsaF1KwUWPK62Wls5lxmHq/M3c4dE1ZxLjOPz4a15dNhbQl1tlkvJVk9Djz9ZT3TSmpk5yjm6W6456bC3t/NjmMKKexXqtnKWLl+w+dwbKPZaawqv9DClJUH6f7un3y7/gjDOtZl0d9uoG/TGmZHs57UI7BtNrQdCX4hZqcRJqgR5EO15n05rYPI21g5WwxIYS9Oz5chIBzm/w0shWanqbBCi2bupmP0/mAZr/2yk5a1q7LgqW7834BmBPq4WKfDNROMXuudHjM7iTDRqO4x/FTYGfd9v1fKFgNS2IvjEwR9/2s07o+fYnaacisotPDL1uPc8tFynp65GR9Pdz4fHsdXD7anoTPPeClJ1lnYOA2a3wVBTjpNU1hF4xpVOBR5O+46n/ztc8yOY3dS2EvSbBDUuwGW/AcyTpudpkyy8wr5as0her6/jMe/3USBRTP+3tb8+mQ3boqt7hwdGctj/WTIz5IbkgQAfXr2Zq+lFqlrvzE7it05+Zw2G1IK+r0HE7vAH/+GgZ+Ynei69p5KZ/q6I/y4MYm0nAJa16nKP/o1pndsDdydoWd6ReRlwfrPIKYvVI81O41wAF1jwpjq15MHzn6NPnsQFVLP7Eh2I4X9WsIaQaexsOpDaDMM6nQ0O9FVTqfn8uu2E/y05TgJh8/h6a64uVkEwzvVJa5usOuOzq+06RtjNlOXp8xOIhyEUooaXYbD4q85uHQa9e/8t9mR7EZpE+6yjIuL0/Hx8Xbfb7nkZsAn7cE3BB5ZaiyKbbLDKZksSUxm8a5kVu8/g0VDo+qB3NmmFoPbRjpfs66KKsyHca2N9UwfXCgLVYuL8gosbH+jCzXc06j58nan/9lQSiVoreOu9zyrVCml1M3AR4A78LnW+i1rbNcheAcYF1JnjTCmQHa0/2INyWk5rDt4lrUHUlhzIIUDpzMBaBDmz5juDRjQqhaNarjgxdDS2jYLzh817kFw8l9cYV1eHm6kxQyizZ7/sH/rChq0vMHsSHZR4cKulHIHPgF6A0nABqXUT1rrnRXdtsOIHQANesKfb0CzO42pkDaQX2jhcEoW+5Iz2HMqne3HzrPt2HlOnM8BjLvq4qKCub9DXXo2Dicq1N8mOZyKxQIr/2eshiXNvkQxWt8ykrzdb3Js6VQp7GXQHtintT4AoJT6DhgAuE5hVwpueRcmdIRFr8Idn5bqZYUWTU5+IVl5hWTlFZCWXcD57HxSs/M4k57L6YxcktNyOZaaTdK5bI6nZlNg+evUWP1Qf9rXC6F5rSDaRYXQtGYVPJyxfa4tJf5itFwePEVG66JYQcGh7AzuSuzZRZw4l05EsOu/u7VGYa8FHL3k30lABytst8IsFs2ZjFxOpeVyKi2H0xm5nMvK43xWPqlZ+WTkFZCdV0hmbgG5BRbyCy3kFVgotGgKLJpCi0ZrzYVSO0bdzogtMxi9oxlb3Jpw4RGLNvZl0cbrCgo1+YWWy4p0cdzdFNX8vYgM9qVV7arc1iKCBmEBxFQPoEFYAP7e5p/Pd2haw4r3IaQ+xA40O41wYOFdhhM6fykzF85iyD2u32rCGpWjuGHSVRVNKfUI8AhAnTp1rLDbv+QWFLLnZAa7T6Wz91Q6e5MzOHI2i6Nns8gtsFz1fB9PN4J8PfH39sDfywNfL3cCfTzw9nDD090NdzeFh5vCzU3hphQKYzC4xzKac3tX8i/3KXzc4HMsygOljMfclMK96PleHm54uis83d3w83LH18sDP093qvh6EuTrSRVfD8ICvAn288LN1ach2tL+JXBiM9w+DtzczU4jHFho69vJWFCFgN0/kJE7nAAXHzRZ47tLAi5dsj4SOH7lk7TWk4BJYMyKqcgOz2bmsWZ/ChsOnWXT0VR2HU8jr9Ao4F4ebtQP9Sc6LIBejcOJDPalRpAv4YHehFcxiqmPZwWKwM73CP5+OG/WXm/KhVRRRGtY9o4xE6blULPTCEfn4UV2w9vpuWsWs9YkMvzGZmYnsilrFPYNQIxSqh5wDBgK3GuF7V7lm7WHmb7uCLtOpAHg5+VO81pBPNA1iha1qtI4IpC6IX62PQ/dpP9fF1Kb3gGBLtDm1hkdWgFH1xo3kXlUsumdolzCOg+HxOkcXvU9Bd1iXfp6VYULu9a6QCn1OPAbxnTHKVrrHRVOVoyM3AKq+nryXJ+GdGoQSovIIPuvxXnhQurETsaF1Ds/s+/+hWHZOxBQA1oPMzuJcBa1O5DlV4vu6UtYuGMMt7WoaXYim7FKVdRa/6q1bqi1bqC1fsMa2yzOmO4NmPFIRx7vGUPbusHmLbAcGm30I9n6HRxaZU6GyuzwamPE3uUpWUhDlJ5S+LS9hy7uO5i9NAEzbs60F9d9L2Jr3Z6FoNrw63PGnY/Cfpa9A/5hRs91IcrAreVQ3LEQfWoB8YfPmR3HZqSwl5eXn7EAdvJOWCenY+zm6AY48KfxjsnLz+w0wtmExmCJaMMgz9VMWn7A7DQ2I4W9IhrfCtG9YembkHbVRCBhC0v/C37VZNk7UW5uLYfQhIMcTkzgwOkMs+PYhBT2ilAK+r1jnIr57R9mp3F9h9cYc9e7PG308BGiPJoNQit3Bnus5IuVB81OYxNS2CsqpL5xvn3HHNi32Ow0ru3PNyCgOrR7yOwkwpkFhKGie3G391p+SDhCSkau2YmsTgq7NXR5CkIaGBdS83PMTuOaDi43ZsJ0/ZucWxcV12IIVfOTaW3ZwddrD5udxuqksFuDpw/0exfOHoBVH5mdxvVoDUvegMCaMhNGWEejfuAVyGMhCXy95jA5+c6/aP2lpLBbS3Qv407UFe9Dyn6z07iW/YuNu0xveFbmrQvr8PKD2P50yl1JRmYGP248ZnYiq5LCbk193zRub5//rDHKFBVnscDi1yCoDrQebnYa4UpaDMEjP4MHQ3fx+YoDWK7TjdWZSGG3pioR0PMVY571ttlmp3ENO+fAiS3Q82Xw8DI7jXAlUV0hsCYj/ddy4Ewmf+w6ZXYiq5HCbm3tRkHNNvDbS5Dtune22UVBHiz+D4Q3heZ3mZ1GuBo3d2hxN+HJK2kalMvkFa5zw5IUdmtzc4fbP4SsFPjj32ancW4bp8G5g3DTv6XfurCNlkNRupBX6u5iw6FzbDziGoMxKey2ENESOj4GCVONm2pE2eVmGD1h6naBmN5mpxGuKrwJRLSkXdpvVPHxYLKLtBmQwm4rPf4BVevAz0/K3PbyWDsBMpPhpv+TtUyFbbUYivvJLTzdspCFO05y6Eym2YkqTAq7rXj5w23/MxZaXvG+2WmcS9oJWPmhsahJ7XZmpxGurvlgUO4M8V6Np5ubS5xrl8JuS9E3QYuhsPIDOGWTtUdc05L/gCUfer9mdhJRGQSEQ/RN+Cf+yOA2EcxOSOKMk7cZkMJuaze/CT5V4acnwOJad7fZxPHNsPlb6DAGQuqZnUZUFi2HQNoxHq9/krxCC1+tPmR2ogqRwm5rfiFwy9twLAHWfGJ2GsemtdEl068a3PCc2WlEZdKoH3hXoeahefRuUp2v1h4mK6/A7FTlJoXdHpoNgsa3wZLX4fQes9M4rl0/weFVxs1IPkFmpxGViacvNB0IO+fxaOcIUrPy+X7DUbNTlZsUdntQCm79wOhPMe8xOSVTnLws+O2fEB4rrQOEOVreA/mZtM5cQbuoYCavOEh+ocXsVOUihd1eAqvDLe9C0gY5JVOcFe/B+SNw6/vg7mF2GlEZ1ekEwVGw+VtG39CAY6nZzN96wuxU5SKF3Z6aD4ZGtxqnZJITzU7jOM7shVXjjBFT3c5mpxGVlVLGz+DB5fSMyCMmPIBPl+1HO2FDPyns9qSU0W7AOwB+fNjohVLZaW10w/T0k+mNwnwthwIat20zGdO9AYkn01m6+7TZqcpMCru9BYRD/4/h5FZjYebKbsePcHAZ9HrFODZCmCk4ymhjsWUG/VtGUDPIh4nLnG99BSnsZmh8K7QeZtxdeXi12WnMk3UWFr5k9NaJe9DsNEIYWg6FlH14ntjIqG71WX/wLAmHnas5mBR2s9z8JgTXhR9HQ855s9OY47eXjS6Y/cdL90bhOGIHgocvbPmWoe1qU9XPk4lLnWvULoXdLN6BcOdkSDsGPz9V+VZc2rsItnwLXZ6GiBZmpxHiLz5VoMntsP0H/N0KGNEpij92nWLPqXSzk5WaFHYz1W4PPf8JO+ZAwpdmp7GfnDTjj1loI+j+vNlphLhaq3uNd9KJvzCycxR+Xu5ONWqXwm62Lk9Dg16w4EU4uc3sNPax6FVIPwEDJxhrxArhaOp1h6DasHk6wf5e3Nu+Dj9tOc7Rs1lmJysVKexmc3ODOz4D32CYNRJyneftXrnsXmC8O+k0FiLjzE4jRPHc3Iw57fv/hPNJPNStPu5K8dly5xi1S2F3BAFhMOhzOHsA5j7muufb004Y31+N5sai30I4slb3Ahq2zKBGkA+D2tbi+/gkktMcf+EcKeyOol434wadXT+55sIcFgvMHQP52TBoipyCEY4vpB7U7Wq0kdaa0Tc0oKDQwhcrD5qd7LqksDuSTo9D87uMlgN7F5mdxrrWfgIHlhrTPMMamp1GiNJpfZ/xTvrIWqJC/bmtRU2+WXuYc5mOfdd4hQq7UuoupdQOpZRFKSUnTCtKKbh9HNRoBj+MgjP7zE5kHYdXwx//NloXtx1pdhohSi92AHgFwOZvAHisRwMy8wr50sEX4qjoiH07cCew3ApZBBitfYdMBzcPmD4IMs+Ynahizh+D74cbt2oPnCALUwvn4uUPTe+A7XMgN53GNarQt2l1vlx1kLScfLPTlahChV1rvUtrvdtaYUSR4Lpwz0xIPwnfDjF6lTujglz4fphxXn3IdFk8QzinNsMhP9O43wR4omcM6TkFfL3msMnBSibn2B1V7XbGTJljCUYnSGdbnONC18ZjCTBwIoQ3NjuREOUT2c64mW7j1wA0qxVEj0ZhfL7iAJm5jrl83nULu1LqD6XU9mI+BpRlR0qpR5RS8Uqp+NOnna8Npima3A43vwWJv8AvzzjXNMgV78Gmr6HbcxDb3+w0QpSfUsaoPWn9xXUUnugVw7msfL5dd8TkcMW7bmHXWt+ktW5WzMe8suxIaz1Jax2ntY4LCwsrf+LKpuMYozhunAYLnneO4r7xK2NmT4uh0ONls9MIUXEth4KbpzFYAdrUCaZrdCifLT9Adp7jvZuWUzHOoOc/jamQ6yfB7/907OK+e4HRByb6Jhgw3riDTwhn5x8KjfvBlhnGtSPgyV4xnMnI5dv1jjdqr+h0xzuUUklAJ2C+Uuo368QSl1EK+rwO7UfDmvFGu1uLAy6yu/cPoy1CRCu4axq4e5qdSAjraT3caDO9+1cA2tcLoXODany6bD85+Y41aq/orJg5WutIrbW31rq61rqvtYKJKygFt7wNHcYYN/vMGe1YS+vt+hlmDIXQhnDfbGP5PyFcSYMeUCXSONVY5KleMZxOz+WbtY41Q0beJzsTpYyLqb1ehW3fw4whkJthdirYOgu+HwE1W8GIn8G/mtmJhLA+N3dofb/RGOycUcg71K9WNGp3rHPtUtidjVLQ7Vlj1aEDy+Dzm+DMXnOyWCyw7F1jOmbdzjBsDvhWNSeLEPbQZpjxO7hx2sUvPX1TQ85k5DJ9neOM2qWwO6s2w+D+2ZCZDJN6wI659t1/brpx89Gfrxv9be6bZawKJYQrC4qEmD7GnPZC487T9vVC6BLtWKN2KezOrEFPGL3cuPln1gj4+Wn7rJ96fDNM7mXMgOn7X7hzEnj62n6/QjiCtg8YA6rE+Re/9EzRqH3amkOmxbqUFHZnFxQJI381pkNunAbj28POebaZEpmXBb+/ApN7Qk6qceql01jp/yIql5jexkXUS5azjIsKoXvDMD5dtp90B+ghI4XdFXh4Qd834KHFEBBuNN36eiAcWmWd7RcWwNbvYWInWD3OuIA0dj3U726d7QvhTNzcjTtRDyyFlL9WVHquTyNSs/Idol+7FHZXUqsNPPynMXPm1A6Y2g+m3AK7foH8cqz6kpMG8VPg4zbGBVJPPxjxC/QfJxdJReXWZhgo98suojaPDKJv0+p8vuKg6f3alTbhLsa4uDgdHx9v9/1WKvnZxnzbVR9B2jGjp3TDm40LP9WbQmjM1asY5WbA2f1G467EX+HgMijMg1ptjbYGDW+WO0mFuGDGvXB0Hfxt58Xfpd0n07n5o+WMvqEBL95i/cZ3SqkErfV1177wsPqehWPw9IUOoyHuQTi4HHbONUbu22cbjyt3CIww3la6eRh/CNKP//X64Cho/wg06Q+128t5dCGu1O5B2D3fuKbV4m4AGtUIpH/LmkxdfZAHu0QRXsXHlGgyYq9MCgvgzB5I3gnJuyDtOOhCoyWwuxdUqw/VYiA81hjRSzEXomQWC4yPA79q8NBfS1keOpPJTR8sY2j72rw+sLlVdykjdnE1dw+oHmt8CCEqxs0N2j0Ev71kTAGu2QqAqFB/7mlfh2/XH+HBLvWoH2b/9hpywlQIIcqr1b3GpIINky/78hO9ovH2cOP93/eYEksKuxBClJdvVeP8+rbZkHX24pfDA314qFt95m87wZajqXaPJYVdCCEqot3DUJADm6df9uWHu9Wjmr8Xby1IxN7XMqWwCyFERdRoBnU6wYbPL1snIdDHkyd6RrPmQApL99h3OVAp7EIIUVHtH4Fzh2Dv5WsN3duhLlHV/Pjv/F0UFNpvcRwp7EIIUVFN+hv9Y9ZOuOzLXh5uvNSvCXuTM/huw1G7xZHCLoQQFeXuAe0fNm4GPLntsof6xFanQ70Q/rdoD2l2ahAmhV0IIayh7Qhj6uPaTy/7slKKV26L5WxWHhP+3F/Ci61LCrsQQliDb7Axr33b95CRfNlDzWoFcWfrSKasPMjRs1k2jyKFXQghrKXDGKNxXvyUqx76e99G1K3mx8m0cnRaLSMp7EIIYS2hMRDT15j6eEWr7BpBPvz+zA20iwqxeQwp7EIIYU2dxkLmadj63VUPKTs11pPCLoQQ1lTvBohoBavGGZ1TTSCFXQghrEkp6Pq0sWhN4i+mRJDCLoQQ1takP4TUh5Uf2mZh+euQwi6EENbm5g6dn4DjG+HQCvvv3u57FEKIyqDlPeAfZoza7UwKuxBC2IKnrzGvff9iY4UlO5LCLoQQttL+YfAOguXv2nW3UtiFEMJWfIKg46PG7JgrmoPZkhR2IYSwpY5jwLsKLHvbbruUwi6EELbkG2yca9/1M5zaYZddSmEXQghb6/goeAXCsnfssrsKFXal1LtKqUSl1Fal1BylVFVrBRNCCJfhFwIdRsPOeZC8y+a7q+iIfRHQTGvdAtgDvFTxSEII4YI6jYUGPYy2vjbmUZEXa61/v+Sfa4HBFYsjhBAuyi8Ehs2xy66seY79QWBBSQ8qpR5RSsUrpeJPnz5txd0KIYS41HVH7EqpP4AaxTz0stZ6XtFzXgYKgOklbUdrPQmYBBAXF2f/rjhCCFFJXLewa61vutbjSqkRwG1AL61NaGMmhBDiMhU6x66Uuhl4Aeiutbb9Cq1CCCGuq6Ln2McDgcAipdRmpdSnVsgkhBCiAio6KybaWkGEEEJYh9x5KoQQLkYKuxBCuBhlxkQWpdRp4HA5Xx4KnLFiHGuRXGUjucpGcpWNo+aCimWrq7UOu96TTCnsFaGUitdax5md40qSq2wkV9lIrrJx1Fxgn2xyKkYIIVyMFHYhhHAxzljYJ5kdoASSq2wkV9lIrrJx1Fxgh2xOd45dCCHEtTnjiF0IIcQ1OHxhL+0qTUqpm5VSu5VS+5RSL9oh111KqR1KKYtSqsQr3EqpQ0qpbUUtF+IdKJe9j1eIUmqRUmpv0X+DS3heYdGx2qyU+smGea75/SulvJVSM4seX6eUirJVljLmGqmUOn3JMXrITrmmKKWSlVLbS3hcKaXGFeXeqpRq4yC5blRKnb/keL1qh0y1lVJ/KqV2Ff0uPlXMc2x7vLTWDv0B9AE8ij5/G3i7mOe4A/uB+oAXsAWItXGuJkAjYCkQd43nHQJC7Xi8rpvLpOP1DvBi0ecvFvf/seixDDsco+t+/8BjwKdFnw8FZjpIrpHAeHv9PF2y3xuANsD2Eh7vh7EegwI6AuscJNeNwC92PlYRQJuizwMxVpe78v+jTY+Xw4/Ytda/a60Liv65Fogs5mntgX1a6wNa6zzgO2CAjXPt0lrvtuU+yqOUuex+vIq2P63o82nAQBvv71pK8/1fmnc20EsppRwglym01suBs9d4ygDgK21YC1RVSkU4QC6701qf0FpvLPo8HdgF1LriaTY9Xg5f2K9Q0ipNtYCjl/w7iasPpFk08LtSKkEp9YjZYYqYcbyqa61PgPGDD4SX8DyfopW21iqlbFX8S/P9X3xO0cDiPFDNRnnKkgtgUNHb99lKqdo2zlRajvw72EkptUUptUAp1dSeOy46hdcaWHfFQzY9XhXq7mgtVlilqbiRVIWn+5QmVyl00VofV0qFY7Q3TiwaZZiZy+7HqwybqVN0vOoDS5RS27TW+yua7Qql+f5tcoyuozT7/BmYobXOVUqNwXhX0dPGuUrDjONVGhsxbsPPUEr1A+YCMfbYsVIqAPgBeFprnXblw8W8xGrHyyEKu674Kk1JwKUjl0jguK1zlXIbx4v+m6yUmoPxdrtChd0Kuex+vJRSp5RSEVrrE0VvOZNL2MaF43VAKbUUY7Rj7cJemu//wnOSlFIeQBC2f8t/3Vxa65RL/jkZ47qTI7DJz1RFXVpQtda/KqUmKKVCtdY27SOjlPLEKOrTtdY/FvMUmx4vhz8Vo/5apam/LnmVpg1AjFKqnlLKC+Nil81mVJSWUspfKRV44XOMC8HFXr23MzOO10/AiKLPRwBXvbNQSgUrpbyLPg8FugA7bZClNN//pXkHA0tKGFTYNdcV52H7Y5y/dQQ/AcOLZnt0BM5fOPVmJqVUjQvXRpRS7TFqXsq1X1XhfSrgC2CX1vqDEp5m2+Nlz6vF5bzCvA/jXNTmoo8LMxVqAr9ecZV5D8bo7mU75LoD469uLnAK+O3KXBizG7YUfexwlFwmHa9qwGJgb9F/Q4q+Hgd8XvR5Z2Bb0fHaBoyyYZ6rvn/gNYwBBIAPMKvo5289UN/Wx6iUud4s+lnaAvwJNLZTrhnACSC/6OdrFDAGGFP0uAI+Kcq9jWvMFLNzrscvOV5rgc52yNQV47TK1kvqVj97Hi+581QIIVyMw5+KEUIIUTZS2IUQwsVIYRdCCBcjhV0IIVyMFHYhhHAxUtiFEMLFSGEXQggXI4VdCCFczP8Dokc1AWZBL0EAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "def f(x):\n", " return np.sin(2*x) + 2*np.cos(x)\n", "\n", "def dfdx(x):\n", " return 2*np.cos(2*x) - 2*np.sin(x)\n", "\n", "def dfdx_fd(g, x0, dx):\n", " return (g(x0+dx)-g(x0))/dx\n", "\n", "def err(x):\n", " return dfdx(x) - dfdx_fd(f, x, 0.01)\n", "\n", "plt.plot(x, f(x), label='f')\n", "plt.plot(x, 100*err(x), label='err * 100')\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Do\n", "\n", " - Compare error to the function itself. When is the error (its absolute value) largest?\n", " - Choose another function $f$ and modify definitions of `f` and its analytically solved derivative `dfdx` above accordingly. What kind of function gives the smallest error?\n", " - Modify `dfdx_fd` and implement backward or central difference approximation. How does the error change?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Applying finite differences to a physical problem" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If we know the analytical form of $f$, there is often little interest in approximating its derivative numerically. Instead, derivative approximations are used indirectly to approximate the values of $f$ itself. Let's see how this works in a simple example of a falling sphere.\n", "\n", "Consider a solid sphere in a viscous fluid, with $\\rho_\\mathrm{sphere} > \\rho_\\mathrm{fluid}$. The gravity pulls the sphere downwards, causing a buoyancy force \n", "$F_\\mathrm{b} = V\\Delta\\rho g = V (\\rho_\\mathrm{sphere} - \\rho_\\mathrm{fluid}) g$. The viscous fluids resists the movement of the sphere downwards. The drag caused is given by the Stokes's law, $F_\\mathrm{d} = -6\\pi\\eta R v$. The sum of forces acting on the ball is then \n", "$\\sum F = m_\\mathrm{sphere} a = V\\rho_\\mathrm{sphere} a = F_\\mathrm{d} + F_\\mathrm{b} = -6\\pi\\eta R v + V (\\rho_\\mathrm{sphere} - \\rho_\\mathrm{fluid}) g$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Simplifying, and noting that $V = \\frac{4}{3}\\pi R^3$ and acceleration $a=\\frac{dv}{dt}$, gives \n", "$$\\frac{dv}{dt} = -\\frac{9 \\eta v}{\\rho_\\mathrm{sphere} 2 R^2} + \\frac{(\\rho_\\mathrm{sphere} - \\rho_\\mathrm{fluid}) g}{\\rho_\\mathrm{sphere}}$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " When the sphere has reached its terminal velocity, acceleration is zero, so $dv/dt = 0 \\Rightarrow v = \\frac{2}{9}\\frac{\\Delta \\rho R^2 g}{\\eta}$" ] }, { "cell_type": "code", "execution_count": 225, "metadata": {}, "outputs": [], "source": [ "def vel_term_stokes(rhosphere, rhofluid, radius, viscosity):\n", " return (2.0/9.0)*(rhosphere-rhofluid)*radius*radius*9.81 / viscosity" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A marble with radius of 2 cm in honey:" ] }, { "cell_type": "code", "execution_count": 226, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.0052320000000000005" ] }, "execution_count": 226, "metadata": {}, "output_type": "execute_result" } ], "source": [ "vel_term_stokes(2800, 1300, 0.02, 250)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We could solve the problem by integrating both sides of the expression for acceleration. However, as an example, we will solve this numerically. To apply finite differences, we replace all derivatives with their finite approximations at time $t_0$. In this case,\n", "$$\\frac{dv}{dt}\\Big\\rvert_{t=t_0} \\approx \\frac{v(t_0+\\Delta t) - v(t_0)}{\\Delta t}$$\n", "\n", "Inserting this into the previous formula gives $\\frac{v(t_0+\\Delta t) - v(t_0)}{\\Delta t} = -\\frac{9 \\eta v(t_0)}{\\rho_\\mathrm{sphere} 2 R^2} + \\frac{(\\rho_\\mathrm{sphere} - \\rho_\\mathrm{fluid})g}{\\rho_\\mathrm{sphere}}$. Note that on the right hand side $v$ has been replaced by $v(t_0)$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In other words, if we know velocity at some time $t_0$ (i.e. $v(t_0)$), we can calculate the velocity after time $\\Delta t$ has elapsed: \n", "\n", "$$v(t_0 + \\Delta t) = \\left(-\\frac{9 \\eta v(t_0)}{\\rho_\\mathrm{sphere} 2 \\pi R^2} + \\frac{(\\rho_\\mathrm{sphere} - \\rho_\\mathrm{fluid})g}{\\rho_\\mathrm{sphere}}\\right) \\Delta t + v(t_0)$$\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's write this as a python function:" ] }, { "cell_type": "code", "execution_count": 227, "metadata": {}, "outputs": [], "source": [ "def vel_next(vel_prev, dt, rhosphere, rhofluid, radius, viscosity):\n", " return (-9*viscosity*vel_prev / (rhosphere*2*radius*radius) + (rhosphere-rhofluid)*9.81/rhosphere)*dt + vel_prev" ] }, { "cell_type": "code", "execution_count": 242, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(0.0005255357142857143,\n", " 0.000998283242984694,\n", " 0.0014235449708098922,\n", " 0.0018060906768669342,\n", " 0.002150211032985211,\n", " 0.0024597657283326785,\n", " 0.0027382267600849766)" ] }, "execution_count": 242, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dt = 0.0001\n", "v0 = 0.0\n", "v1 = vel_next(v0, dt, 2800, 1300, 0.02, 250)\n", "v2 = vel_next(v1, dt, 2800, 1300, 0.02, 250)\n", "v3 = vel_next(v2, dt, 2800, 1300, 0.02, 250)\n", "v4 = vel_next(v3, dt, 2800, 1300, 0.02, 250)\n", "v5 = vel_next(v4, dt, 2800, 1300, 0.02, 250)\n", "v6 = vel_next(v5, dt, 2800, 1300, 0.02, 250)\n", "v7 = vel_next(v6, dt, 2800, 1300, 0.02, 250)\n", "v1, v2, v3, v4, v5, v6, v7" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can automate the calculation of the steps:" ] }, { "cell_type": "code", "execution_count": 245, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAD8CAYAAACYebj1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8VfWd//HXJxthDVtkN4CACIoLEVFb61IVl5HOjFOxm53q2AV/9ffQTn/6a2ep/c382v5musyMnRn3vUhtrdRWRatWW9YEWcMWQmISCIQkhCVkufd+fn/c22sIgVzg3pwk9/18PPLg3nO+59zPl8B533O+ZzF3R0REBCAj6AJERKTnUCiIiEicQkFEROIUCiIiEqdQEBGROIWCiIjEKRRERCROoSAiInEKBRERicsKuoCTMXLkSJ84cWLQZYiI9BrFxcX73D0/0fa9KhQmTpxIUVFR0GWIiPQaZlZxMu11+EhEROIUCiIiEqdQEBGROIWCiIjEKRRERCROoSAiInEKBRGRHqC4ooGH3ymluKIhoemp0quuUxARSYbiigZWlNUxd/IIZhcMS9r0cMRpDUVoDUVYXVHP6vJ6zhuXx7RRg6PTwxE2VTeyobqRqWcMpmDEANrCzrY9B/npu6WEwk5mhvH5SwsYPSSXivomFq+uJBxx+mVn8Pxdc4/6/FRQKIhIr3DSG+zyev5Quo/zJwxl6qjBHGkNc6Q1zNrKBr77agltYScr07jzY5M4Y3AupbWH4hvgjAzjE9NGMjg3m5rGI6wubyDiYAYFwweQmWEcbG5j78HW+OdlGESS8Mj7UMR58o/lx0xvC0VYUVanUBCRvqmzjbm7s3xHHX8o3cf00YOZMHwAB5tDrKvcz7+9vT3+TfrG80YzsF82lfVNLNuxL77BHpOXSyQCB5rbaGoNd1lDW9j5r9+XHTM9HHFW7WxgxKAcDreE4ht7d8jMMKaPHkL5vsPxUDBgzqThXDp5JDlZGawqr+fdLXtxomFx06wx/NmssbxRUsPLa6qJeHT65+YW8NlLCti25wDf+Pl6QuEIWZkZPPr5QmZPHMbG6kbueHIVbaEI2VkZzJ08Ill//celUBCR09ZxAx+JOAea23h/ey0ryuo5c/gARgzqx/6mVvY3tbF970HeLNnz0cZ8SC4toQj7j7QSjpz4s0IR57cbahg6IIdQJHLUBntIbjbnjx9Kae1B1lTsx4lusK+bOYobzxvDgJwsquqb+L+vbyEUjpCdmcG/334hcyYNZ/PuA/z1U6vjG+CnvzSH2QXDKK5o4LOPrYhP/8Gt53c6/W+vnx4PtzkVw1m+Y1983hcvm8TsgmGMGNSP36zfHZ8+/4JxnD16MGePHszYoQOOCclLJo/g+bvmdronlCrmnoT9nW5SWFjouveRSOp19i3+cEuId7buZVlpHeOH5TKkfw77DrWwZfdBlpbURDfwwODcLA61+3bdUYZBTlYGzW0fbf2njx7MRQXD2LH3EKt21se/Yf/5heO4fc6ZVNYf4YFfRr9JZ2d9dGy944a5q+kn6l8yp5/qMqlgZsXuXphw+0RCwczmAT8BMoHH3P17Heb3A54BZgN1wG3uXh6b9yBwJxAGvu7ub8SmlwMHY9NDiRStUBBJruKKBpbt2Me0UYMZOSiH3Y3NFJXX8+yKD6PH1g1GD8nlQHOIQy2hTteRm330Bv788XlcMS2fTbsO8E67Qyh3fmwS91w9lcH9svigcv9Jb8yTucFOJ0kPBTPLBLYB1wJVwGrgdncvadfma8Asd/+KmS0A/tzdbzOzGcDPgDnAWOAtYJq7h2OhUOju+xItVqEgcnzH2wj+sXQfb5bsYUxeLv1zMqluOELV/iNsqznI9r2Hulzv2aMGc+lZIyirPcT72/fFN/Jf/sRZ3HftNNZXNXbLt3U5NScbComMKcwBSt29LPYBi4D5QEm7NvOBf4y9fgn4DzOz2PRF7t4C7DSz0tj6lidaoIgcreNG0915d+tevvzsGtrCETIzjEvPGsGhlhBltYdoPHL0N/yczAzGDs2l/fdBA/5i9jj+5uOT2dvYwt3PFcU35v/8F+fFN/Kryuvj0z95ziiyMzOYXTCs0+Pex5v+J7MLhnW60T/edOkeiYTCOKCy3fsq4JLjtXH3kJk1AiNi01d0WHZc7LUDS83Mgf9290dOvnyRvqvjxr81FOG1Dbv5xkvrCIWdDDMm5Q9k74FmDjR/tOEPRZx1lfs5d1weBSMGsqGqMf7t/iufOItvXHc2GRl2zDf5z8wpYProIUwfzUlv5LWB7zsSCQXrZFrHY07Ha3OiZS93911mdgbwppltcff3jvlws7uBuwHOPPPMBMoV6d3CEec3G3Zx/+Loxt8Mxg7tT01jM6F2o7dhd9pCEeZfMI6sTOP5FRWEIk5OVgZP/nXnZ85cc84oMjKi/y21kZfOJBIKVcCEdu/HA7uO06bKzLKAPKD+RMu6+5/+3GtmLxM9rHRMKMT2IB6B6JhCAvWK9BpNrSF+9UE1b2/Zi2HUHmpha81BjrR9dI69O/TLyuTuKyaTnZnBf/5+B+HYWTg/vO2C+Eb65lljk3YIR9JXIqGwGphqZpOAamAB8JkObZYAdxAdK7gVeNvd3cyWAC+Y2Q+JDjRPBVaZ2UAgw90Pxl5fBzyUlB6J9EDFFQ38sbSWUUNyaQ076yv3s6G6ka01B4/a7Z4xZjAL5kxgQE4mj76/M77x/8Gts+Ib7yum5evbvaRMl6EQGyO4B3iD6CmpT7j7JjN7CChy9yXA48CzsYHkeqLBQazdYqKD0iFgYezMo1HAy9GxaLKAF9z99RT0T6RbtR8HmDCsP8UVDfx2w25e3bD7qIHd4QNzmDU+j2EDsllRFj0vP9PgplljWXjVFACunj7qpDb+Ismgi9dEksDdeXX9Lu5bvI62sGN8NHiWmWGEY2MBGQZ/8/HJPHDDdMyOHeztjhueSXpJxSmpItJBcXk9r2+qITvT2LW/meVldew50BKf78BVZ+fz9Wum0hqKHHX/mutmjia2l9zlMX+R7qZQEEnQweY2lu2o46XiKt4q2RPfE8jrn83Hp45k3ND+PLWsPH4rhnuunsqFZ0Y38hrsld5CoSByHO7OL9dU88raavYdamHbnkOEIk52psUDIcPg7ismsfCqqQBcN3O0xgGkV1MoSNprPzh8/vg8Vpc3sLSkhl+v28W+Qx/dGvlTF47jtosnYNDhdsYj4+vSxl96O4WCpLXiigY+++gKWkIRMswYkJPBwZYwOVkZTBjWn7pDrfGrgaecMSh+P3uNA0hfpVCQtBQKR/jjjjp+8PoWmkPRO3yG3ZkwfCBfv2YKH5+az5aag0edGdT+ASfaI5C+SqEgaaOovJ5frqmm/nArq8vrqTvcyoCcTDLNcKK3h/jup85N+Gpgkb5IoSB9Xk1jM//2u238bFVlfID4sskjuOPyiXwidt9/nRkkEqVQkD4nekuJfWRlwOryBn6/rfaop4BlGlw+dSTXzxwNaMMv0p5CQfqUNzbV8LXn18SvIB42IJuvXnkWM8YM4f6fr+vWB6CL9EYKBen13J1VO+t5alk5r2+sOeoagi99bBL/4+roNQSj8/prfECkCwoF6bWW79jH08sq2FxzgIq6JvL6ZzP/gnG8tnF3/Kriy87SNQQiJ0OhIL3OgeY2vvfaFl5Y+SEQvbDsq5+YzNevmUb/nEw+X1GgPQKRU6RQkF6j7lALT/6xnKeXl3Ow3eMnMwwG5WbTPycT0B6ByOlQKEiPt3RTDf/93g42VDXSFnHmzRzN1dPP4O9e2aiBY5EkUyhIj7XvUAv/8MpGfrOhBoBMM35y2wXccsE4ACbnD9JhIpEkUyhIj3OwuY3H3t/JY++X0dQabjfHqWw4En+nw0QiyadQkB7hT88wbmhq45W1u6g/3MqN541m3szRfPMX63WYSKSbKBQkcMXl9Sx4dAVt4egVBrPG5fHkFy/m/AlDARg3bIAOE4l0E4WCBGrz7gPc++LaeCBkGFx/7uh4IIAOE4l0J4WCBKKxqY0fvbWNZ5aXMzAnk+xMIxJxHSISCZhCQbpVUXk9j75fxrIddRxuCfHZSwq4/7pp7Kg9rENEIj2AQkG6zZK11dz74lrcwQz+362zuHX2BABmF+QoDER6AIWCpFwoHOHR93fyL0u34rG71WUAew60BFqXiBxLoSAptaXmAH/78/VsqG7kkknDWVu5P36zOo0diPQ8CgVJiZVldfzkd9tZubOOof1zePgzF3HTrDEUVzRo7ECkB1MoSNItWVfNvT9bixM9xfRfPn0+V519BqDTS0V6uoygC5C+w935eVEl9y9eF3/QjQEluw4EWZaInATtKUhSHGxu49u/2sgra3cxc+xgSvce1tiBSC+UUCiY2TzgJ0Am8Ji7f6/D/H7AM8BsoA64zd3LY/MeBO4EwsDX3f2NdstlAkVAtbvffNq9kW5XXNHArz6o5s2SGmoPtXL/tdP42lVTWFu5X2MHIr1Ql6EQ23A/DFwLVAGrzWyJu5e0a3Yn0ODuU8xsAfB94DYzmwEsAGYCY4G3zGyau//p1pf3ApuBIUnrkXSb4vJ6bntkBaGIY8D/+dS5fHZuAaCxA5HeKpExhTlAqbuXuXsrsAiY36HNfODp2OuXgGvMzGLTF7l7i7vvBEpj68PMxgM3AY+dfjekux1pDfOtX20kFPnonkX7j7QFXJWInK5EQmEcUNnufVVsWqdt3D0ENAIjulj2x8A3gchJVy2Bqmpo4i//cxlbag6SlWFkGho7EOkjEhlTsE6meYJtOp1uZjcDe9292MyuPOGHm90N3A1w5plndl2tpNSKsjq+9vwa2sIRnvzixQzpn62xA5E+JJFQqAImtHs/Hth1nDZVZpYF5AH1J1j2FuAWM7sRyAWGmNlz7v65jh/u7o8AjwAUFhZ2DCPpJsXl9fzne2W8s3kPBSMH8tgXCpmcPwhAYSDShyRy+Gg1MNXMJplZDtGB4yUd2iwB7oi9vhV42909Nn2BmfUzs0nAVGCVuz/o7uPdfWJsfW93FgjSM6zaWcenH1nBWyV7cOA7t8yMB4KI9C1dhkJsjOAe4A2iZwotdvdNZvaQmd0Sa/Y4MMLMSoH7gAdiy24CFgMlwOvAwnZnHkkvcLglxDd/sYFwbEDZgPVVjcEWJSIpY+6954hMYWGhFxUVBV1G2qg71MKXnlrN+qpGsto9BOf5u+bqkJFIL2Fmxe5emGh7XdEsnaqsb+ILT6xi1/4jPPqFQoYNzNGAskgaUCjIMTbtauSLT66mNRThhb+5hNkFwwENKIukA4WCxBVXNLC46EOWrN3N0AHZvPCVS5k6anDQZYlIN1IoCBANhAWPLKctHL1lxXduuUCBIJKGdOtsAeDZ5eW0hT+6ZcX2vYeCLUhEAqE9BWHpphp+vW4XZtFvCbplhUj6UiikuaWbalj4whrOGz+U+6+bxvqqRp1hJJLGFApp7M2SPSx8YQ0zxubxzJ1zGJKbzcen5gddlogESKGQhoorGnh+RQWvrKvm3HFDeTYWCCIiCoU0c9RZRgb3XztNgSAicTr7KM38ck3VR2cZARuqdR8jEfmI9hTSSOneg7yythojetqpzjISkY4UCmliz4Fm7nhiNbnZWfzgL89nZ91hnWUkIsdQKKSBxiNt3PHEKvY3tfLily/l3HF5QZckIj2UQqGPawmF+fKzRZTuPcSTf32xAkFETkih0IcVldfzrZc3snXPQX582wW6BkFEuqRQ6KOKKxq47ZEVhCNOVoYxYfiAoEsSkV5Ap6T2UY++VxZ/hKa7s6KsLuCKRKQ30J5CH7Tmwwbe2ryHDIs+U1mnnopIohQKfczeA8185dlixgzN5bvzz2XTrgM69VREEqZQ6ENaQmG+/FwxB5tDPHPnZUwfPYQrzz4j6LJEpBdRKPQR7s7f/2oTH3y4n59+9iKmjx4SdEki0gtpoLmPeG7lh7xYVMk9V03hxvPGBF2OiPRS2lPoA55dXs4/LNnE7IJh3HfttKDLEZFeTHsKvdxbJXv4u1c2EXHYVN3IB5X7gy5JRHoxhUIvFok43321JP6+LRzR9QgicloUCr3Yf79XRkV9E9mZRqZuhS0iSaAxhV5qbeV+/nXpVm46bwxfunwiK3bW63oEETltCoVe6GBzG1//2QeMGpLLP//FeeT1z2b2xOFBlyUifUBCh4/MbJ6ZbTWzUjN7oJP5/czsxdj8lWY2sd28B2PTt5rZ9bFpuWa2yszWmdkmM/tOsjrU17k73/7VRqr3H+Hfbr+AvP56vrKIJE+XoWBmmcDDwA3ADOB2M5vRodmdQIO7TwF+BHw/tuwMYAEwE5gH/DS2vhbganc/H7gAmGdmc5PTpb7tl2uqeWXtLu69ZiqzC7R3ICLJlcjhozlAqbuXAZjZImA+UNKuzXzgH2OvXwL+w8wsNn2Ru7cAO82sFJjj7suBQ7H22bEfP82+9GnFFQ28tnE3zy6vYM6k4Sy8akrQJYlIH5RIKIwDKtu9rwIuOV4bdw+ZWSMwIjZ9RYdlx0F8D6QYmAI87O4rT6UD6aC4ooHPPraC5rYIAHd9bBKZGRZwVSLSFyUyptDZ1qfjt/rjtTnusu4edvcLgPHAHDM7t9MPN7vbzIrMrKi2tjaBcvueFWV1tMQCIcNg+95DXSwhInJqEgmFKmBCu/fjgV3Ha2NmWUAeUJ/Isu6+H3iX6JjDMdz9EXcvdPfC/Pz0fJzkmLzceArn6FoEEUmhREJhNTDVzCaZWQ7RgeMlHdosAe6Ivb4VeNvdPTZ9QezspEnAVGCVmeWb2VAAM+sPfBLYcvrd6XtC4QhP/rGcvP5Z3HPVFJ6/a66uRRCRlOlyTCE2RnAP8AaQCTzh7pvM7CGgyN2XAI8Dz8YGkuuJBgexdouJDkqHgIXuHjazMcDTsXGFDGCxu7+aig72do/9YScbqht5+DMXcdMs3f1URFLLol/oe4fCwkIvKioKuoxuU1Z7iBt+8j5Xnp3Pf31uNtETukREEmdmxe5emGh73fuoh4pEnP/1i/X0y8rgu/PPVSCISLdQKPRQz62sYHV5A3938wzOGJIbdDkikiYUCj1QZX0T33ttC1dMy+fW2eODLkdE0ohCoYcpLq/n84+vJBJx/vnPddhIRLqXQqEHKa5oYMGjKyivayIUcfYcaAm6JBFJMwqFHuT32/bSFo6eDebueoqaiHQ7hUIPsrP2MBC9lYWeoiYiQdBDdnqI0r0HeW1jDZ885wwuPHOYnqImIoFQKPQA7s53fl1C/5xMvv+XsxgxqF/QJYlImtLhox5gacke3t++j/uunaZAEJFAKRQC1twW5ruvljBt1CA+N7cg6HJEJM3p8FHAHn2vjKqGI7xw1yVkZyqjRSRY2goFqHr/ER5+t5QbzxvNZVNGBl2OiIhCISjFFQ186anVhCPO/77xnKDLEREBdPgoEMUVDdz+yApawxGyMow9B1oYP2xA0GWJiGhPIQjLd+yjNRx95rKuXBaRnkShEIBI7LlGunJZRHoaHT7qZm3hCC9/UM2Zwwbw6YvHc+lZI3Xlsoj0GAqFbvZScRU79x3m0S8Ucu2MUUGXIyJyFB0+6kbNbWF+8tZ2LjpzKJ8854ygyxEROYZCoRs9s7ycmgPN/O310/XwHBHpkRQK3eRAcxs/fXcHV0zL59KzNLAsIj2TQqGbPPZeGfub2vjm9WcHXYqIyHEpFLrBvkMtPPaHndx03hjOHZcXdDkiIselUOgGD79TSksown3XTQu6FBGRE1IopNjrG3fzzLIKrjo7n7PyBwVdjojICSkUUqi4ooGFz39A2J33t++juKIh6JJERE5IoZBCb2ysIezRe1qEwhHd40hEejyFQgqV1x0GIFP3OBKRXiKhUDCzeWa21cxKzeyBTub3M7MXY/NXmtnEdvMejE3fambXx6ZNMLN3zGyzmW0ys3uT1aGeYu+BZt7dVsu154zivuvO5vm75uoeRyLS43V57yMzywQeBq4FqoDVZrbE3UvaNbsTaHD3KWa2APg+cJuZzQAWADOBscBbZjYNCAH3u/saMxsMFJvZmx3W2as9/oedhMIRvn3zORSMGBh0OSIiCUlkT2EOUOruZe7eCiwC5ndoMx94Ovb6JeAai97HYT6wyN1b3H0nUArMcffd7r4GwN0PApuBcaffnZ6hsamN51ZUcPOssQoEEelVEgmFcUBlu/dVHLsBj7dx9xDQCIxIZNnYoaYLgZWJl92zPb28nMOtYb565VlBlyIiclISCYXO7tzmCbY54bJmNgj4BfA/3f1Apx9udreZFZlZUW1tbQLlBqupNcSTf9zJJ885g3PGDAm6HBGRk5JIKFQBE9q9Hw/sOl4bM8sC8oD6Ey1rZtlEA+F5d//l8T7c3R9x90J3L8zPz0+g3GD9bFUlDU1tfPXKKUGXIiJy0hIJhdXAVDObZGY5RAeOl3RoswS4I/b6VuBtd/fY9AWxs5MmAVOBVbHxhseBze7+w2R0pCdoCYV59L0y5k4erjONRKRX6jIUYmME9wBvEB0QXuzum8zsITO7JdbscWCEmZUC9wEPxJbdBCwGSoDXgYXuHgYuBz4PXG1ma2M/Nya5b93u5TXV1Bxo5mvaSxCRXsrcOw4P9FyFhYVeVFQUdBmdCkeca/71XQbnZrPknsv1EB0R6RHMrNjdCxNtryuak+Thd0opr2ti3sxRCgQR6bUUCklQXF7Pj97cBsC/v1OqG9+JSK+lUEiCF4sq4+fZtoV04zsR6b0UCklQVqsb34lI39DlvY/kxD6sa6L4wwZunT2OSSMHMXfyCJ2OKiK9lkLhND27opxMM75x3XRG5+UGXY6IyGnR4aPT0NQa4sXVlcw7d7QCQUT6BIXCaXj5g2oONIf44mUTgy5FRCQpFAqnyN15elk5M8cO0RiCiPQZCoVTtHxHHdv2HOKLl03UxWoi0mcoFE7RU8vKGT4whz87f2zQpYiIJI1C4RRU1jfx1uY93D5nArnZmUGXIyKSNAqFU/DcigrMjM/NLQi6FBGRpFIonKQjrWEWra5k3szRjMnrH3Q5IiJJpVA4Sb9aW03jkTbu0GmoItIHKRROQnF5PT9cupWJIwZw8USdhioifY9CIUHFFQ3c/uhKag+1Ur3/CGs+3B90SSIiSadQSNCKsjpawxEAIhHX7bFFpE9SKCRo1vg8AAzdHltE+i7dJTVBH9Y3AfC5uQV86sJxurWFiPRJCoUELV5dyfTRg3lo/kzd1kJE+iwdPkrAlpoDrKtq5NOFExQIItKnKRQSsHh1FdmZxqcuHBd0KSIiKaVQ6EJrKMLLH1Rx3YzRDB+YE3Q5IiIppVDowlub99DQ1MZfFY4PuhQRkZRTKHRhcVElY/Jy+fjU/KBLERFJOYXCCexuPMJ722q5dfZ4MjM0wCwifZ9C4QR+UVxFxOHW2Tp0JCLpQaFwHJGIs7ioirmTh1MwYmDQ5YiIdIuEQsHM5pnZVjMrNbMHOpnfz8xejM1faWYT2817MDZ9q5ld3276E2a218w2JqMjybZyZz0f1jdx28UTgi5FRKTbdBkKZpYJPAzcAMwAbjezGR2a3Qk0uPsU4EfA92PLzgAWADOBecBPY+sDeCo2rUf6eVElg/tlMW/mmKBLERHpNonsKcwBSt29zN1bgUXA/A5t5gNPx16/BFxj0Ut/5wOL3L3F3XcCpbH14e7vAfVJ6EPSvb+9liXrdnHpWSPon6NnMItI+kgkFMYBle3eV8WmddrG3UNAIzAiwWV7lOKKBr701GpCEefdbbUUVzQEXZKISLdJJBQ6OxfTE2yTyLIn/nCzu82syMyKamtrT2bRU7KirI62cLTEcDii5yaISFpJJBSqgPajreOBXcdrY2ZZQB7RQ0OJLHtC7v6Iuxe6e2F+fuovIJt6xiBAz00QkfSUSCisBqaa2SQzyyE6cLykQ5slwB2x17cCb7u7x6YviJ2dNAmYCqxKTumpUVEXfW7CnR+fxPN3zdVzE0QkrXT5PAV3D5nZPcAbQCbwhLtvMrOHgCJ3XwI8DjxrZqVE9xAWxJbdZGaLgRIgBCx09zCAmf0MuBIYaWZVwD+4++NJ7+FJ+vX6XZw3Lo9v39TxBCsRkb4voYfsuPtvgd92mPb37V43A391nGX/CfinTqbfflKVdoOd+w6zvqqRb914TtCliIgEQlc0t/Pquuhwx02zdG2CiKQnhUKMu7Nk3S7mTBzO2KH9gy5HRCQQCoWYLTUH2b73EH92vvYSRCR9KRRifr1uF5kZxg3nKRREJH0pFIgeOvr1+l1cdtYIRg7qF3Q5IiKBUSgAayv3U1l/hFvOHxt0KSIigVIoAEvW7SInM4Przx0ddCkiIoFK+1AIR5zfrN/NlWfnMyQ3O+hyREQClfahsHJnHXsPtnDLBTp0JCKS9qHw63W7GZCTyTXTRwVdiohI4NI6FFpDEV7buJtrZ4zSw3REREjzUHh6WTn7m9qYMWZI0KWIiPQIaRsKxRUNfO+1LQD86K1tesKaiAhpHArLduwj7NEnrLWF9IQ1ERFI41AYlBO9a3iG6QlrIiJ/ktDzFPqi8rrD5GQaC6+awsem5usJayIipGkouDtLS/Zw5dlncO8npwVdjohIj5GWh482Vh9gd2Mz187QtQkiIu2lZSi8WVJDhsE15ygURETaS8tQWFqyh4snDmf4wJygSxER6VHSLhQ+rGtiS81BHToSEelE2oXC0pIaAK6bodtki4h0lIahsIfpowdz5ogBQZciItLjpFUo1B9upai8nut06EhEpFNpFQq/27yHiMN1M3XoSESkM2kVCktL9jA2L5eZY3VXVBGRzqRNKBxpDfP+9lqunTEKMwu6HBGRHiltQuH97bU0t0V06EhE5ATSJhSWluxhSG4WcyYND7oUEZEeK6FQMLN5ZrbVzErN7IFO5vczsxdj81ea2cR28x6MTd9qZtcnus5kCoUj/G7zHq6efgbZmWmTgyIiJ63LLaSZZQIPAzcAM4DbzWxGh2Z3Ag3uPgX4EfD92LIzgAXATGAe8FMzy0xwnUlTXNFAQ1ObDh2JiHQhka/Nc4BSdy9z91ZgETC/Q5v5wNOx1y8B11h0NHfSqxlSAAAFPUlEQVQ+sMjdW9x9J1AaW18i60yaZ1dUkJlhDOmfnaqPEBHpExIJhXFAZbv3VbFpnbZx9xDQCIw4wbKJrDMpisvr+c363YQjzl1Pr9azmEVETiCRUOjs/E1PsM3JTj/2w83uNrMiMyuqra09YaGd+UPpvviK9SxmEZETSyQUqoAJ7d6PB3Ydr42ZZQF5QP0Jlk1knQC4+yPuXujuhfn5+QmUe7SPTc0nNzuDTD2LWUSkS4k8jnM1MNXMJgHVRAeOP9OhzRLgDmA5cCvwtru7mS0BXjCzHwJjganAKqJ7Cl2tMylmFwzj+bvmsqKsjrmTR+hZzCIiJ9BlKLh7yMzuAd4AMoEn3H2TmT0EFLn7EuBx4FkzKyW6h7AgtuwmM1sMlAAhYKG7hwE6W2fyuxc1u2CYwkBEJAHm3umh/B6psLDQi4qKgi5DRKTXMLNidy9MtL2u5BIRkTiFgoiIxCkUREQkTqEgIiJxCgUREYnrVWcfmVktUHGKi48E9iWxnN5C/U4v6nd6SaTfBe6e8JW/vSoUToeZFZ3MaVl9hfqdXtTv9JKKfuvwkYiIxCkUREQkLp1C4ZGgCwiI+p1e1O/0kvR+p82YgoiIdC2d9hRERKQLvTIUzGyemW01s1Ize6CT+f3M7MXY/JVmNrHdvAdj07ea2fWJrrMnSFG/nzCzvWa2sXt6cfKS3W8zm2Bm75jZZjPbZGb3dl9vTk4K+p5rZqvMbF2s79/pvt4kLhX/1mPzMs3sAzN7NfW9OHkp+j9ebmYbzGytmXV9R1F371U/RG+1vQOYDOQA64AZHdp8Dfiv2OsFwIux1zNi7fsBk2LryUxknUH/pKLfsXlXABcBG4PuYzf+vscAF8XaDAa29bTfdwr7bsCgWJtsYCUwN+i+dse/9dj8+4AXgFeD7md39RsoB0YmWkdv3FOYA5S6e5m7twKLgPkd2swHno69fgm4xswsNn2Ru7e4+06gNLa+RNYZtFT0G3d/j+gzMHqqpPfb3Xe7+xoAdz8IbCZFzwg/Tanou7v7oVj77NhPTxtYTMm/dTMbD9wEPNYNfTgVKen3yeqNoTAOqGz3vopj/0PH27h7CGgERpxg2UTWGbRU9Ls3SGm/Y7vfFxL9xtzTpKTvsUMoa4G9wJvu3tP6nqrf+Y+BbwKR5JecFKnqtwNLzazYzO7uqojeGArWybSO33SO1+Zkp/ckqeh3b5CyfpvZIOAXwP909wOnXGHqpKTv7h529wuIPht9jpmde1pVJl/S+21mNwN73b34dItLoVT9W7/c3S8CbgAWmtkVJyqiN4ZCFTCh3fvxwK7jtTGzLCCP6CGS4y2byDqDlop+9wYp6beZZRMNhOfd/Zcpqfz0pfR37u77gXeBecksOglS0e/LgVvMrJzoYZmrzey5VBR/GlLy+3b3P/25F3iZrg4rBT24cgqDMVlAGdHBlD8Nxszs0GYhRw/GLI69nsnRgzFlRAd3ulxn0D+p6He75SbScweaU/H7NuAZ4MdB9y+AvucDQ2Nt+gPvAzcH3ddU97vDslfSMweaU/H7HggMjrUZCCwD5p2wjqD/Ik7xL+9GomeM7AC+FZv2EHBL7HUu8HOigy2rgMntlv1WbLmtwA0nWmdP+0lRv38G7AbaiH7buDPofqa638DHiO5arwfWxn5uDLqf3dT3WcAHsb5vBP4+6D5217/1dvOvpAeGQop+35OJhsU6YFMi2zZd0SwiInG9cUxBRERSRKEgIiJxCgUREYlTKIiISJxCQURE4hQKIiISp1AQEZE4hYKIiMT9f8YuP9Sz9O8HAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "Velocity at time 0.0049 s is 0.005202758993971469\n" ] } ], "source": [ "dt = 0.0001 # time step size\n", "v0 = 0.0 # initial condition\n", "nt = 50 # num of tsteps to calculate\n", "\n", "v = np.zeros(nt) # stores calculated velocities\n", "t = np.zeros(nt) # stores total time in seconds\n", "\n", "v[0] = v0\n", "t[0] = 0.0\n", "\n", "for it in range(1,nt):\n", " v[it] = vel_next(v[it-1], dt, 2800, 1300, 0.02, 250)\n", " t[it] = t[it-1] + dt\n", " \n", "plt.plot(t, v, '.-')\n", "plt.show()\n", "\n", "print(\"Velocity at time \", t[-1], \" s is \", v[-1])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Do\n", "\n", " - Use the code above to investigate what happens if the ball is being dropped from an elevation just slightly above the surface of the fluid. How do you need to modify the code?\n", " - Take the finite differences approximation of the equation and replace $v(t_0)$ on the right hand side by $v(t_0 + \\Delta t)$. Rearrange the equation so that you have an expression for $v(t_0 + \\Delta t) = \\ldots$. \n", " - Implement a function `vel_next_implicit`, similar to `vel_next` that calculates the velocity at next time step, but uses this modified expression.\n", " - Copy and modify the code above so that you create a storage for both velocity values calculated using `vel_next` and `vel_next_implicit`, and plot the outcome on top of each other. Vary the value of `dt` and see what happens. " ] }, { "cell_type": "code", "execution_count": 246, "metadata": {}, "outputs": [], "source": [ "# Implement `vel_next_implicit` here\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The first approximation (`vel_next`) is the forward difference approximation. `vel_next_implicit` is the backward difference approximation. When applied to functions of time, these are also called *explicit* and *implicit* approximations." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.4" } }, "nbformat": 4, "nbformat_minor": 2 }