{ "cells": [ { "cell_type": "markdown", "id": "b2d834fe", "metadata": {}, "source": [ "## Floating Point Numbers Review\n", "\n", "Reminder that __normalised__ floating point numbers are stored as\n", "\n", "$$ (-1)^s \\times 1.d_1 d_2 ... d_n \\times 2^\\text{exp} $$ \n", "\n", "or equivalently\n", "\n", "$$ (-1)^s \\times \\bigg( 1 + d_1 \\cdot 2^{-1} + d_2 \\cdot 2^{-2} + ... + d_n \\cdot 2^{-n} \\bigg) \\times 2^\\text{exp} $$ \n", "\n", "where $s=0,1$ is the sign bit, $d_i = 0,1$ are the mantissa bits and the exponent $\\text{exp}$ ranges from $\\text{exp}_{\\text{min}}$ to $\\text{exp}_{\\text{max}}$. Normalised means that $d_0 = 1$." ] }, { "cell_type": "markdown", "id": "01271621", "metadata": {}, "source": [ "For example, the Float32 data structure is normalised and has 1 sign bit, 8 exponent bits, and 23 mantissa bits (also known as precision bits). The exponent is given as the number represented in the exponent bits - 127." ] }, { "cell_type": "markdown", "id": "fc35b883", "metadata": {}, "source": [ "### Question 1" ] }, { "cell_type": "markdown", "id": "14051675", "metadata": {}, "source": [ "What is $\\frac{1}{2}$ in Float32?" ] }, { "cell_type": "markdown", "id": "ee1cf02c", "metadata": {}, "source": [ "Sign bit $s = 0$\n", "\n", "$0.5 = 1.\\underbrace{0000000}_{\\text{mantissa}} \\times 2^{-1}$ in binary\n", "\n", "So $ \\text{exp} = 127-1 = 126 = 2^6+2^5+2^4+2^3+2^2+2^1 = 01111110 $ in binary\n", "\n", "So representation in Float32 $= \\underbrace{0}_{\\text{sign}} ~ \\underbrace{01111110}_{\\text{exponent}} ~ \\underbrace{0000...000}_{\\text{mantissa}}$" ] }, { "cell_type": "markdown", "id": "7ad78296", "metadata": {}, "source": [ "### Question 2\n", "\n", "What is $0.75$ in Float32?" ] }, { "cell_type": "markdown", "id": "2652fdf8", "metadata": {}, "source": [ "Sign bit $s = 0$\n", "\n", "$0.75 = 0.5 + 0.25 = 1.\\underbrace{100000}_{\\text{mantissa}} \\times 2^{-1}$ in binary\n", "\n", "So $ \\text{exp} = 127-1 = 126 = 2^6+2^5+2^4+2^3+2^2+2^1 = 01111110 $ in binary\n", "\n", "So representation in Float32 $= \\underbrace{0}_{\\text{sign}} ~ \\underbrace{01111110}_{\\text{exponent}} ~ \\underbrace{1000...000}_{\\text{mantissa}}$" ] }, { "cell_type": "markdown", "id": "195767b7", "metadata": {}, "source": [ "### Question 3" ] }, { "cell_type": "markdown", "id": "6b5758f5", "metadata": {}, "source": [ "What is the smallest (positive) integer not represented exactly using Float32?" ] }, { "cell_type": "markdown", "id": "b162a795", "metadata": {}, "source": [ "Any positive power of 2 less than $2^{255-127=128}$ is represented exactly. This is because you can just set the mantissa to zero and the exponent power to be what you want.\n", "\n", "This means that $2^{23} + 1$ is the largest integer not represented exactly. Any integer less than $2^{23}$ can be represented exactly using the 23 mantissa bits, while $2^{23}$ cannot be it is a power of 2 so can be done also by setting the mantissa to zero and setting the correct exponent." ] }, { "cell_type": "markdown", "id": "bf30f8b8", "metadata": {}, "source": [ "### Question 4\n", "\n", "What is the smallest Float32 number greater than 1?" ] }, { "cell_type": "markdown", "id": "327450c2", "metadata": {}, "source": [ "1 in Float32 $ = 1.\\underbrace{000000}_{\\text{mantissa}} \\times 2^{0}$ in binary\n", "\n", "So next biggest number is just increasing the mantissa by smallest thing possible\n", "$ = 1.\\underbrace{000000...01}_{\\text{mantissa}} \\times 2^{0}$ in binary $= 1 + 2^{-23} $" ] }, { "cell_type": "markdown", "id": "0b870195", "metadata": {}, "source": [ "## Differential Equations Review" ] }, { "cell_type": "markdown", "id": "a1125c4f", "metadata": {}, "source": [ "The most basic case is solving a scalar ODE with initial condition $y(x_0) = y_0$\n", "\n", "$$ y'(x) = f\\big(x,y(x)\\big) $$ \n", "\n", "The first step is to apply the fundamental theorem of calculus and write this as an integral equation\n", "\n", "$$ y(x) = y_0 + \\int_{x_0}^x y'(s) ds $$\n", "\n", "Plug in the formula\n", "\n", "$$ y(x) = y_0 + \\int_{x_0}^x f\\big(s,y(s)\\big) ds $$\n", "\n", "Then the goal is to approximate the integral term using whatever method you like." ] }, { "cell_type": "markdown", "id": "0ea4bb6f", "metadata": {}, "source": [ "### Question 1\n", "\n", "Approximate the integral using a left Riemann sum. What is the resulting scheme?" ] }, { "cell_type": "markdown", "id": "8fef125e", "metadata": {}, "source": [ "This is Forward Euler. Approximate integral is\n", "\n", "$$ \\int_{x_0}^x f\\big(s,y(s)\\big) ds \\approx (x-x_0) f\\big(x_0,y_0\\big) $$\n", "\n", "Which gives the scheme\n", "\n", "$$ y(x) = y_0 + (x-x_0) f\\big(x_0,y_0\\big) $$" ] }, { "cell_type": "markdown", "id": "60a37355", "metadata": {}, "source": [ "### Question 2\n", "\n", "Approximate the integral using the trapezoidal rule. What is the resulting scheme?" ] }, { "cell_type": "markdown", "id": "0029d248", "metadata": {}, "source": [ "Approximate integral is\n", "\n", "$$ \\int_{x_0}^x f\\big(s,y(s)\\big) ds \\approx (x-x_0) \\bigg( f\\big(x_0,y_0\\big) + f\\big(x,y(x)\\big) \\bigg) * 0.5 $$\n", "\n", "Which gives the scheme\n", "\n", "$$ y(x) = y_0 + (x-x_0) \\bigg( f\\big(x_0,y_0\\big) + f\\big(x,y(x)\\big) \\bigg) * 0.5 $$" ] }, { "cell_type": "markdown", "id": "52c0d3ef", "metadata": {}, "source": [ "### Question 3\n", "\n", "Consider the case now where $f'(x) = 2f(x)$. Write code to find the solution at $x=1$ using the trapezoidal rule. Use $\\Delta x = 10^{-2}$ and initial condition $f(0) = 1$.\n", "\n", "Plot the computed solution and the exact solution (you need to solve for the exact solution)." ] }, { "cell_type": "markdown", "id": "2f8b9dff", "metadata": {}, "source": [ "Plugging in to the equation above gives the scheme:\n", "\n", "$$ y(x) = y_0 + (x-x_0)( 2y_0 + 2y ) * 0.5 $$\n", "$$ y(x) = y_0 + (x-x_0)y + (x-x_0)y_0 $$\n", "$$ (1-x+x_0)y(x) = y_0 + (x-x_0)y_0 $$\n", "$$ y(x) = \\frac{(1+x-x_0)y_0}{1-x+x_0}$$" ] }, { "cell_type": "code", "execution_count": 1, "id": "3b280081", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "trapezoidal (generic function with 1 method)" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using PyPlot\n", "\n", "function trapezoidal()\n", " \n", " dx = 1e-2;\n", " N = round( 1/dx );\n", " xvals = [0.0];\n", " fvals = [1.0];\n", " x = 0.0;\n", " \n", " for ii = 1:N\n", " x += dx;\n", " push!( xvals, x );\n", " push!( fvals, ( fvals[end] * (1+dx) ) / (1-dx) );\n", " end\n", " \n", " plot( xvals, fvals )\n", " plot( xvals, [exp(2*x) for x in xvals] )\n", " \n", "end" ] }, { "cell_type": "code", "execution_count": 2, "id": "1704bb3b", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhYAAAGdCAYAAABO2DpVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA58ElEQVR4nO3dZ3hUZcLG8fukF0iogUBC6CWBAOkBFQvIqqhYEBUV7KygKPbVd5Vdd7HuWlZRUFFUiopgowgqRaUlhBpqCBBK6KSQMClz3g9AVlbahJk5mcn/d13zYc6c4dw8hpzb0x7DNE1TAAAATuBjdQAAAOA9KBYAAMBpKBYAAMBpKBYAAMBpKBYAAMBpKBYAAMBpKBYAAMBpKBYAAMBp/Ny9Qbvdrl27dqlu3boyDMPdmwcAANVgmqaKiorUrFkz+fic/riE24vFrl27FB0d7e7NAgAAJ8jLy1NUVNRpP3d7sahbt66kY8HCwsLcvXkAAFANhYWFio6OrtqPn47bi8WJ0x9hYWEUCwAAPMzZLmPg4k0AAOA0FAsAAOA0FAsAAOA0FAsAAOA0FAsAAOA0FAsAAOA0FAsAAOA0FAsAAOA0FAsAAOA0FAsAAOA0FAsAAOA0FAsAAOA0FAsAALzEr2NH6JdPRslmO2pZBrfPbgoAAJxv85olSt05QX6GXeuyUtQp7QpLcnDEAgAAD2fa7Tr6zePyM+zKqnORZaVColgAAODxMmZ9rM5lK3XU9FezAa9amoViAQCABys5UqiopS9IklbGDFGTmA6W5qFYAADgwVZMGqVI7Ve+0Vhdb37O6jgUCwAAPNWO3PVKyPtYkpSf9qyCQupanIhiAQCAx9rz5WMKMsqVHdhVXfvcYXUcSRQLAAA80soFXyvxyEJVmD6q0/9VGT41Y5deM1IAAIBzZrOVKvznv0iSMiNuUItOKRYn+i+KBQAAHiZjymi1NHfooMLUadCLVsc5CcUCAAAPsnvHFnXLeVeSlNvtCYXVa2RxopNRLAAA8CA7Jo9UqGHTRv+OSrjmAavj/AHFAgAAD7F64bdKLv5ZdtOQ39X/kuHja3WkP6BYAADgAcpsNtX9+WlJUmbj/mod39PiRKdGsQAAwANkfP6iWtrzdFh11eHWl62Oc1oUCwAAarg9O3MVv3mMJGlz/KMKaxBhcaLTo1gAAFDD5U16RHWMUm3066CEax+yOs4ZOVQsWrZsKcMw/vAaNmyYq/IBAFCrrZz3lZKKf1alacjv2tfl41vzLtj8PT9HVl62bJkqKyur3q9Zs0Z9+vTRgAEDnB4MAIDa7mhpiRrMf0aSlNFkgFK79LA40dk5VCwaN2580vsXX3xRbdq0Ua9evZwaCgAASJmTRqmnuUv7VU9xg16yOs45cahY/F5ZWZk+/fRTjRw5UoZhnHY9m80mm81W9b6wsLC6mwQAoNbIy8lW4rYPJEPKS35W3cMbWB3pnFT74s3p06fr8OHDGjJkyBnXGz16tMLDw6te0dHR1d0kAAC1gmm368AXIxRklGttYDd1u+JuqyOdM8M0TbM6X+zbt68CAgL07bffnnG9Ux2xiI6OVkFBgcLCwqqzaQAAvFrmrAlKXPygykxf7b3tJ0W162Z1JBUWFio8PPys++9qnQrZtm2b5s6dq6+++uqs6wYGBiowMLA6mwEAoNYpKjioqMXPS5Kyom5Xag0oFY6o1qmQ8ePHKyIiQldddZWz8wAAUKut+ewpNdEB7TSaqOugf1gdx2EOFwu73a7x48dr8ODB8vOr9rWfAADgf2zMWqiUPZ9Lkg5d/KKCQupYnMhxDheLuXPnavv27brrrrtckQcAgFqporxc+u4R+Rqmlte9RJ17XW91pGpx+JDD5Zdfrmpe7wkAAE5j2RevKL1yk4oUrJhb37Q6TrUxVwgAABbL37FFXTYcKxPrYh9Vw8gWFieqPooFAAAW2zlpxLFJxvw7KOmGkVbHOS8UCwAALJT5w0QlHlmgCtNHAde+VeMnGTsbigUAABYpLDioqN+elSRlNrtVLTunWpzo/FEsAACwyNpPHv/vMytuf9HqOE5BsQAAwALrMn5U6r6pkqTDl76soJC6FidyDooFAABuZrMdVeCMR+RjmMoM76u4C/tbHclpKBYAALhZ5qS/qbV9mw6rrtre/obVcZyKYgEAgBtt37RKibljJUk5Cc8ovFGkxYmci2IBAICb2CsrVfT5Awo0yrUmKEEJ/e63OpLTUSwAAHCTpVP/rbjy1SoxA9Xw5ndk+Hjfbtj7/kYAANRAu/NyFLf2VUnSmo4PKbJlJ4sTuQbFAgAAFzPtduVPfEB1Tzy2e8BTVkdyGYoFAAAuljHjQ3UvXawy01dBN4yRj5/Dk4t7DIoFAAAudHDvLrXJGCVJymp5j1p0TLQ4kWtRLAAAcKGcTx9SAxUq1ydGCYP+ZnUcl6NYAADgIsvnTlZy4RxVmoYq+r0p/4AgqyO5HMUCAAAXKDi4T1G/PC1Jymh2q9olXGxtIDehWAAA4ALrJzyoCB1UntFMXW9/2eo4bkOxAADAyVb8/IVSD8+U3TR05Io3FRRSx+pIbkOxAADAiQoLDipy/pOSpGVNBqhjSh+LE7kXxQIAACda99EINdEB7TSaKP6O16yO43YUCwAAnGTlgq+VeugbSVLh5a8ruE6YxYncj2IBAIATFBYcVMRPIyVJSxtdr07pV1qcyBoUCwAAnGDd+AcVqf3aZTRR5yH/tjqOZSgWAACcp6yfvlDq4e8kSYV9X1dInXrWBrIQxQIAgPNQcHCfmi14QpK0JOImdUyrnadATqBYAABwHjZ8NExNTjwIa8i/rI5jOYoFAADVlPnDZ0opnK1K01DplW8pKKSu1ZEsR7EAAKAaDu7dpZjf/iJJymg2SO2Te1ucqGagWAAA4CDTblfux/erkQ5rm0+0ug2uPXOBnA3FAgAABy37bqwSjyxQuemrimvHKDAo1OpINQbFAgAAB+Tn5ahj5t8kSctb3as2XS+0OFHNQrEAAOAc2Svt2vfpPQozjmijX3slDvq71ZFqHIoFAADnaMkXL6uLbbmOmv4KGfi+/PwDrI5U41AsAAA4B9s2rlK3dcdmK13VaaSi2nW1OFHNRLEAAOAsystsOvr53Qo2yrQmsLuSb3rS6kg1FsUCAICzWDbhL+pQsVGFClHE7R/I8PG1OlKNRbEAAOAM1i39USl5H0qSNif/XRFRbSxOVLNRLAAAOI2iwkOqO/MB+Rl2ZYb1VsJV91gdqcajWAAAcBrZHw5XlJmvfDVSu7veszqOR3C4WOzcuVO33XabGjZsqJCQEHXr1k2ZmZmuyAYAgGUyZ3+q1MPfyW4aOtT3LYXVa2R1JI/g58jKhw4dUs+ePXXJJZdo5syZioiIUE5OjurVq+eieAAAuN/endvUatFTkqRlzQYpNf1KixN5DoeKxUsvvaTo6GiNHz++alnLli2dnQkAAMvYKyuVP+FOxatIW3xbqfvgV62O5FEcOhXyzTffKCkpSQMGDFBERIS6d++ucePGuSobAABut3jSC4q3Zeqo6S+/mz5UQFCw1ZE8ikPFYsuWLRozZozatWun2bNna+jQoXrooYc0YcKE037HZrOpsLDwpBcAADXR5lW/KWnTm5Kk1Z2fVIsOCRYn8jyGaZrmua4cEBCgpKQk/fbbb1XLHnroIS1btkyLFi065Xeef/55jRo16g/LCwoKFBYWVo3IAAA4X8mRQu17NV0x5g6tCOmhro99L8OHmydPKCwsVHh4+Fn33w6NWGRkpGJjY09a1qlTJ23fvv2033n66adVUFBQ9crLy3NkkwAAuMWqD4Yrxtyhfaqvlnd+SKmoJocu3uzZs6c2bNhw0rKNGzcqJibmtN8JDAxUYGBg9dIBAOAGy3/4VGkHv5Yk7b3sdcU1jrQ4kedyqI498sgjWrx4sf75z39q8+bNmjhxosaOHathw4a5Kh8AAC61Oy9HrX87NqnYkqa3Ku7C/tYG8nAOFYvk5GRNmzZNkyZNUufOnfX3v/9dr7/+ugYNGuSqfAAAuExFebkOTrhD9VSszX5tlXDXv62O5PEcOhUiSf369VO/fv1ckQUAALdaOuEv6lG+RkfMIIXcMkH+AUFWR/J4XJkCAKiV1vw2U6nbjz2LaX3S82rWJs7iRN6BYgEAqHUO7d+jRj8Ml69hKiO8rxKv/rPVkbwGxQIAUKuYdrtyP7xTTbVfeUYzxd471upIXoViAQCoVRZPflEJJb+qzPRV2XXjFFKnntWRvArFAgBQa2zMWqjEDa9JklZ0elRt4i+wOJH3oVgAAGqFgsMHFPrN3QowKrQitKeSb3ra6kheiWIBAPB6pt2uTe/fpebmHu02GqvVPR/zyG4XYVQBAF5vyZevKql4nspNXxX3G6fw+o2tjuS1KBYAAK+2edUidV/7siRpefsRapd4icWJvBvFAgDgtQoLDipw2l0KNMq1MjhNKbf8n9WRvB7FAgDglUy7XZvGDla0uUv5aqSWXFfhFowwAMArLZ78ohKPLFCZ6auiq99XeMOmVkeqFSgWAACvsy7jJyVueFWStKLjo1xX4UYUCwCAVzm0f4/qfXevAoxKZdW5SMkDeV6FO1EsAABew15ZqW0f3K5I7dcOI1Jt7/mI6yrcjNEGAHiNxROeVbfSJbKZ/iq7frzq1mtodaRah2IBAPAKKxdMV+rWMZKkVV2fVesu6RYnqp0oFgAAj7d7+yZF/zRcvoapZfX7Kfn6h62OVGtRLAAAHs12tESFEwapgYqU49tGXe59z+pItRrFAgDg0bLGDVOHig0qVKhCbp+ooJA6Vkeq1SgWAACPtfTrd5V24CtJ0taL/q3Ilh0tTgSKBQDAI21e9Zu6LD8298eS6LsUf+lAixNBolgAADzQ4f17FDxtiIKNMq0OTlby4FesjoTjKBYAAI9SWVGh7e/fqubmHu00mijm3ony8fOzOhaOo1gAADzKkvGPKf5ohkrNAJXf8InCGkRYHQm/Q7EAAHiMzFmfqMfO8ZKkdcn/UMvOqRYnwv+iWAAAPMLW9cvVYdHjkqQlTQYqod99FifCqVAsAAA1XsHB/fKdMkh1jFJlB8Qr8e63rI6E06BYAABqtMqKCm0de7OizV3KNxor8t4p8gsItDoWToNiAQCo0ZZ8OFJdjy5TqRmgkusmqH7jZlZHwhlQLAAANday799Xj10fS5Kyk/+p1vE9LE6Es6FYAABqpM2rFilu6V8kSUsib1Niv3stToRzQbEAANQ4+/fsUOhXtyvEsGl1UKKS7n7D6kg4RxQLAECNYrOVau/7AxWpfdphRCrmvsny5cmaHoNiAQCoMUy7XSvG3KPY8jUqUrDMWybzZE0PQ7EAANQYiyf/U6mHv1OlaWjrxf9RdPtuVkeCgygWAIAaYdX8aUrZ8KokKaP9I+py8Y0WJ0J1UCwAAJbbtnGlWv48TL6GqYx6f1LKLf9ndSRUE8UCAGCpw/v3yGfSzQrTEW3w76QuQz+U4cPuyVPxXw4AYJkym015Y2+qelx343u+UGBQqNWxcB4oFgAAS5h2u7LevUddylboiBmkowMmqkGTaKtj4TxRLAAAllg86R9KPfSN7KahnF5vqGVsitWR4AQOFYvnn39ehmGc9GratKmrsgEAvFTWj1OUsvE1SdKy9iMVf+nNFieCszj8KLO4uDjNnTu36r2vr69TAwEAvNvmVYvUfsFD8jVMLWtwtVJuedbqSHAih4uFn58fRykAANWyZ2eu6n41SKHGUa0N7KZuQ9/nDhAv4/B/zU2bNqlZs2Zq1aqVbr75Zm3ZsuWM69tsNhUWFp70AgDUPsVFh1X04Q1qogPa5hOt6D9PlX9AkNWx4GQOFYvU1FRNmDBBs2fP1rhx45Sfn68ePXrowIEDp/3O6NGjFR4eXvWKjuaKXwCobSrKy7XpnZvVtjJHBxUm/9u/VFi9RlbHggsYpmma1f3ykSNH1KZNGz3xxBMaOXLkKdex2Wyy2WxV7wsLCxUdHa2CggKFhYVVd9MAAA+y6J37lL53io6a/tp29RR1SLrM6khwUGFhocLDw8+6/z6veWhDQ0PVpUsXbdq06bTrBAYGKjAw8Hw2AwDwYIsm/kPpe6dIktalvaLulAqvdl5XzNhsNq1bt06RkZHOygMA8CKZsz5R6oZXJElLWj+o7lfcaXEiuJpDxeKxxx7T/PnzlZubqyVLlujGG29UYWGhBg8e7Kp8AAAPtX7Zj4pb9Ih8DFNLG16rlNv+ZnUkuIFDp0J27NihW265Rfv371fjxo2VlpamxYsXKyYmxlX5AAAeKC9nrZp8P1hBRrlWBacogdtKaw2HisXkyZNdlQMA4CUO7t0lfXqj6qtIm33bqO2wL+TnH2B1LLgJ9REA4DSlR4q0d+z1ijZ3abfRWPXuna6QOvWsjgU3olgAAJyiorxM6/8zQB0r1qlQoSob+LkaNW1hdSy4GcUCAHDeTLtdmWPuUffSRbKZ/tp1xYeK6ZhgdSxYgGIBADhviz7+i1IPfi27aSi7x7/UMfVPVkeCRSgWAIDzsuSrN9Vj2xhJUkbsU+re9w6LE8FKFAsAQLVl/ThFiSufkyQtaTZYKQOfsjgRrEaxAABUS/bSueq4YLj8DLsywvsq5Z7XrY6EGoBiAQBw2JbsZWo+4w4FG2VaFZyibsM+4QFYkESxAAA4aPe2jarz+U0K1xFt8O+odsOmyi+AySZxDMUCAHDODu7dpbKPrlOEDmqbT7SaDv1awXVOP4U2ah+KBQDgnBQXHda+sdcqxtyhPWqooLumK7xhU6tjoYahWAAAzupoaYly37pWHSo26rDq6OjNX6pJVFurY6EGolgAAM6oorxM2W/dpC5lK3TEDNK+az7jqZo4LYoFAOC07JV2Zb49WAklC1Vm+mlrn7Fql3Cx1bFQg1EsAACnZNrtWjJ2mFIPz1ClaSi75+uKu+Baq2OhhqNYAABOadHHf1H6nomSpKxuf1e3y2+3OBE8AcUCAPAHiya+UDX/x5J2jyrpugctTgRPQbEAAJxkydQ3lL7xFUnS4hb3K3XQXy1OBE9CsQAAVMn4/gMlrTo+qVjTW5U65EWLE8HTUCwAAJKkrLmT1XXp4/I1TC1tcI1S7nub+T/gMH5iAABauWC6YhcOl79Rqcyw3kp6YDylAtXCTw0A1HJrfpup9j/eq0CjXCtCe6rr8Iny8fOzOhY8FMUCAGqx9ct+VKvZQ6qmP499kJlKcX4oFgBQS21asVDNvrtNocZRrQ3spvYPTlNAULDVseDhKBYAUAttXr1YjaffrDCjROv849T6wW8VFFLH6ljwAhQLAKhltqxdqoZTb1Q9FWujX3tFP/idguuEWR0LXoJiAQC1SG52pup9cYPqq0ib/dqq6fCZqhPWwOpY8CIUCwCoJbauz1Ldz69XAxUqx7eNIh6YqbB6jayOBS9DsQCAWmD7ppWqM7m/Gumwtvi2UqMHZiisQYTVseCFKBYA4OW2bVih4M+uVSMdVq5vSzUYOlPhDZtaHQteimIBAF5s24YVCpnUX411SLk+Map3/wzVaxxpdSx4MR6tBgBeatuGLIVOOnb6I9enper9eabqN25mdSx4OY5YAIAX2rp+eVWp2OLbUvUpFXATigUAeJnc7AzVmXydGumwcnxbqcHQWapHqYCbcCoEALzI5tWL1HDqANVXkbb4tlLDP89SvUZcqAn3oVgAgJfYtGKhGk+/WfVUrM1+bRXx55kKa8gtpXAvToUAgBdYn/mzmkwbqHoq1ga/jooYPptSAUtwxAIAPNy6pXMU/f3tqmOUap1/nKIf/I7HdMMyHLEAAA+2auE3avH9INUxSrU2oKtiRsygVMBSHLEAAA+VNXeyYhcOV6BRrtVBiWr30NcKCqlrdSzUchyxAAAPlDHjQ3Ve+IACjXKtCOmp9g9/R6lAjUCxAAAPs2T6f9R9yUj5G5XKDLtMnR+epsCgEKtjAZLOs1iMHj1ahmHo4YcfdlIcAMCZLJo8WqkrnpGvYWpZ/X7q/tDn8gsItDoWUKXa11gsW7ZMY8eOVXx8vDPzAABOwbTbtejjZ9Rj2zuSpCURNyll6LsyfHwtTgacrFpHLIqLizVo0CCNGzdO9evXd3YmAMDvmHa7Fr03/L+lIvpepQx9j1KBGqlaxWLYsGG66qqr1Lt377Oua7PZVFhYeNILAHBuKisqtPStO9Rjz2eSpCXtH1Xq3a/K8OESOdRMDp8KmTx5spYvX65ly5ad0/qjR4/WqFGjHA4GALWd7WiJVv/nFqUWz1OlaSir29+Uet1DVscCzsihypuXl6cRI0bo008/VVBQ0Dl95+mnn1ZBQUHVKy8vr1pBAaA2KS48pI3/vlJJxfNUZvpqVfrrSqJUwAMYpmma57ry9OnTdd1118nX97/n9SorK2UYhnx8fGSz2U767FQKCwsVHh6ugoIChYWFVT85AHipg3t3af9716p95UaVmIHK7T1OcRdea3Us1HLnuv926FTIZZddptWrV5+07M4771THjh315JNPnrVUAADObPe2jSr/qL/amzt1SHW1v/9niuvey+pYwDlzqFjUrVtXnTt3PmlZaGioGjZs+IflAADH5GZnKuTzAWqhA8pXI5XdOlXt2nezOhbgEOYKAYAaIHvJD2o+c4jCdUTbfKIUdNfXahHV1upYgMPOu1jMmzfPCTEAoPZa/sOniv31YQUZ5drg31FNh36t8IZNrY4FVAtHLADAQou/eE3Ja/4uX8PUyuA0dXjwSyYTg0ejWACABUy7XYvGP6keeWMlQ1pW/yp1f+Aj+fkHWB0NOC8UCwBws/Iym7LeGaIeh2dIkhZH3aXUu17jaZrwChQLAHCj4qLDynn7BqUczVClaSgz7hml3fS41bEAp6FYAICb7N+9XYfe76+ulTkqNQO06aI3lHLZrVbHApyKYgEAbrB1/XIFTBmoduZeHVSYDlz7ieITLrY6FuB0FAsAcLE1v3yrFnPvV5iOaIcRKZ/bpqpdmzirYwEuQbEAABdaOv0/6pb1VwUYlVrvH6um932leo0jrY4FuAzFAgBcwLTbtejDx9Vjx/uSIWXWvURxD3ymoOBQq6MBLkWxAAAnO1paotVj7lCPwjmSpMXN7lDK3a/Lh4kaUQtQLADAiQ7u3ak9425Ucnm2KkwfZcX/n9JuGGl1LMBtKBYA4CRb1y+X/5Sb1cnco0KFaHvvMUq+sL/VsQC3olgAgBOsWjBdLX/8s8KMEu0ymqh84BR17tjd6liA21EsAOA8LZ7yspKyR8vPsGu9f6ya3DdVzRo3szoWYAmKBQBUU3mZTZlj71fa/mmSIWWE9VGXByYoMCjE6miAZSgWAFANBQf3Ku+9AUqzrZDdNLS09XCl3v43JhJDrUexAAAHbduwQj6Tb1Fnc5dKzEBtvODfSuszyOpYQI1AsQAAB6z46XO1nj9CYUaJ8o3GKr3xU3XrnGZ1LKDGoFgAwDkw7XYt+vQ5peW8JR/D1Hr/WDW+53M1bRJtdTSgRqFYAMBZlB4p1pp371CPoh8lQ1ra4Gp1u/99BQQGWR0NqHEoFgBwBru3bVTxhFuUXLlZ5aavsuKeVMqAJyTDsDoaUCNRLADgNNb8+p2azfmzIlWoQ6qr3X3fU0qPq6yOBdRoFAsA+B+m3a7Fk/6p5I2vyc+wK8e3tULumKzYmA5WRwNqPIoFAPxO6ZFirXnvTqUX/nBsuvOwPoobOl5BIXWtjgZ4BIoFABy3c8t6lX52q5Irc1Rh+iizw0il3PwMD70CHECxAABJK3/+Ui3nP6TmOqJDCtPuvmOU2qOf1bEAj0OxAFCr2Ssrtfjjp5W2bax8DFMb/dorfMgkxUa1tToa4JEoFgBqrYKDe5U77nb1KF187PkUDa9V13vfZRIx4DxQLADUSptWLFTo13erm7lHNtNfq7s9p5TrHrQ6FuDxKBYAahXTbteSqf9W9zWjFWiUa5fRRKXXf6ik+AusjgZ4BYoFgFqjpLhAa8feo7Tjt5KuCOmhVvdOULP6ja2OBngNigWAWmHr+uUyPx+sZPv2Y7eSth2ulEGjuJUUcDKKBQCvt2z6O4rLel4hhk37VU97+45Rao8rrY4FeCWKBQCvVXqkSKvH3aeUwzMkQ1od2F2Rd36i2KZMdQ64CsUCgFc6dupjiFLs22Q3DS2NuU/Jd/xTvn782gNciX9hALyKabdr2fS31HnlPxRi2HRA9bS7z1tKu+Aaq6MBtQLFAoDXKCo4qPXv36uUormSIa0J7K6md05Q56YtrI4G1BoUCwBeYdOKhQr++l4lm7tVYfooo/Uwpdw2Sj6+vlZHA2oVigUAj2avrNSSiX9X4uY3FWBUKl+Ndfiqd5WW0tvqaECtRLEA4LH2796uXR8PUfrRTMmQlodeqLZ3faimDSOsjgbUWhQLAB5p5U9TFL3gMcWrUKVmgFZ3eUrJ1z/CA68Aizn0L3DMmDGKj49XWFiYwsLClJ6erpkzZ7oqGwD8wdGSYi1++251XXCfGqhQW3xbau8ts5Ry46OUCqAGcOiIRVRUlF588UW1bdtWkvTxxx/r2muvVVZWluLi4lwSEABO2Lx6kXyn3ac0+3ZJ0pKIAep65xsKCg61OBmAEwzTNM3z+QMaNGigV155RXffffc5rV9YWKjw8HAVFBQoLCzsfDYNoJawV1ZqyaS/K3HTWwowKnRA9bTz4tcUf/GNVkcDao1z3X9X+xqLyspKffHFFzpy5IjS09NPu57NZpPNZjspGACcq/y8zdr3yd1KL1txfEbSdLUY8oHiI5pbHQ3AKThcLFavXq309HQdPXpUderU0bRp0xQbG3va9UePHq1Ro0adV0gAtY9ptyvj+7HqkDFKXYyS4xdoPqnk60dyLQVQgzl8KqSsrEzbt2/X4cOHNXXqVL3//vuaP3/+acvFqY5YREdHcyoEwGkd3r9HOR/dp8TieZKkjX7tFTLwfUW162ptMKAWO9dTIed9jUXv3r3Vpk0bvffee04NBqB2Wvnzl4qc/7gidFAVpo8yY+5R4m0vyC8g0OpoQK3m8mssTjBN86QjEgBQHcWFh7T2o4eUevAbSdJ2n+Yqu/pdpXa/yOJkABzhULH4y1/+oiuuuELR0dEqKirS5MmTNW/ePM2aNctV+QDUAmt/m6H6c0Yo1dwrSVoScZO6DvmXgkLqWpwMgKMcKhZ79uzR7bffrt27dys8PFzx8fGaNWuW+vTp46p8ALxYyZFCrfz4MaXu+Vw+hqndRmMd6vO6Unv0szoagGpyqFh88MEHrsoBoJbJXjxbdWc/rHRzl2RISxtcrdghbykyrL7V0QCcB+YKAeBWpUeKtPLjR5Vy/CjFXjVQfq+XlXLJAKujAXACigUAt8leNFN1fxiptONHKZbVu0odhryp+HqNrI4GwEkoFgBcrrjwkNZMGKm0/V9JUtVRimSOUgBeh2IBwKVWzpuqJvOeUJr2S5KWNuinjne8wVEKwEtRLAC4xOH9e7Txk4eUUnDsdvRdRoQOXfqqUi681uJkAFyJYgHAqUy7XRkzP1TrZX9TigpkNw0tazJAXe54Rc3q1LM6HgAXo1gAcJr8vM3aPXG4kksXSZK2+USr9E+vKzWlt8XJALgLxQLAeausqNCyL15R5/VvqLtRqjLTV8tj7lbCrX9TQFCw1fEAuBHFAsB52bzqN1V+M0JpFRslQ1rv30nB17+ttE6JVkcDYAGKBYBqKSku0MpPn1Ly7snyM+wqUrCyY0cq+YZH5ePra3U8ABahWABwWNbcyWr6y7NK1z7JkJbX6aUWt76p1GYtrY4GwGIUCwDnLH/7Ju2ePELdS3499l6NtffCF5Rw2c0WJwNQU1AsAJxVeZlNGVP+qfjNY9TUsKnc9FVGs0HqOugFNa0TbnU8ADUIxQLAGa39bYZC5j6pdPt2yZDW+ccp+Lo3lB6bbHU0ADUQxQLAKe3fvV25k0YquXCOJOmQwrQ5/jEl9R8uw4eLMwGcGsUCwEnKy8uU+cXLitvwHyUbpceenNnoWnW89WUlN2xidTwANRzFAkCVNb98p9Cf/qI0+zbJkDb5tZOuek2p3XtZHQ2Ah6BYAFB+3mbtnPKoEovnSZIOqa42xT2spOselo8fvyYAnDt+YwC12NGSYmVNeUFdt36opoZNlaahjMbXqeMtLyqF0x4AqoFiAdRCpt2urB8mqOmSfyjd3Ft1t0fA1a8pNT7d6ngAPBjFAqhltqxZopJvHldC2UpJ0h41VF7y00q84m4ZPj4WpwPg6SgWQC2xf88O5Ux5WkkHvpWvYeqo6a8VLQar68DnlFQnzOp4ALwExQLwckdLS5T1xWh1yRmnVKO0am6PyAGvKC2mg9XxAHgZigXgpUy7XctnfaTIZS8q3dxz/PbRtqro8w8lpP7J6ngAvBTFAvBC65fNlWY/q8SKdZKkvWqg7d0fU0K/oUxpDsClKBaAF9m5Zb3yv3qy6nkUJWagVrW4Q/EDn1VSnXqWZgNQO1AsAC9waH++Nnz+VyXs+VLNjUrZTUOZ9a9Uq5v+qbRmLa2OB6AWoVgAHuxoSbGyvnxRcTkfKM0okQxpdWCCQq9+UcmdU62OB6AWolgAHqiivEzLv31XMateV7oOSIaU49tKJb3+qi4XXW91PAC1GMUC8CCm3a6sORPVcMmLSrHnSZLy1Vg7Eh5VwlX3cWEmAMtRLAAPsXbRTPn+OEoJx+/0OKw62tD2XnW94XE1DQ61OB0AHEOxAGq4TSsWqnTW84o/miFJKjUDtDLqVsUO+D+l1mtkcToAOBnFAqihtq1frgPfPqeEIwskSeWmr5Y3ulptbvib0prFWJwOAE6NYgHUMDu3ZGvn16OUeHi2YgxTdtPQ8vDeatb/b0ptHWt1PAA4I4oFUEPk523WtmmjlHDgezU3KiVDygrtqfpXjVJSbLLV8QDgnFAsAIvt27VVOdP+roS909XUqJAMaVVQkoIu/z91T7jY6ngA4BCKBWCR/bu2afO0F9Rt7zSlGeWSIWUHxMvnsmcVn9rX6ngAUC0UC8DN9u/ers3T/qFue6ZWFYr1/rGqvOgpxV1wjWQYVkcEgGqjWABusnfnVm2Z/o+TjlCs9++kiguPFQrDx8fqiABw3igWgIvt2ZGj3On/UPd931QVig1+HVV24ZPqfGF/CgUAr0KxAFxk55Z12vHdP9X9wPdqcvwuj3X+caq48El1vuBqCgUAr0SxAJxs24Ys7Z0xWt0Pz1Fzwy4Z0tqALjJ7Pam49KsoFAC8mkPFYvTo0frqq6+0fv16BQcHq0ePHnrppZfUoUMHV+UDPMamFb+ocM5L6l68UDGGefy20WT5X/K44rjLA0At4VCxmD9/voYNG6bk5GRVVFTomWee0eWXX67s7GyFhjIJEmof025X9pLZqpz/atVcHjKkFSE9FNrnKcV372VtQABwM8M0TbO6X963b58iIiI0f/58XXTRRef0ncLCQoWHh6ugoEBhYWHV3TRgKXtlpVb+NFnBS99Sx/Jjs41WmoaywnurYd8n1CouxeKEAOBc57r/Pq9rLAoKCiRJDRo0OO06NptNNpvtpGCApyqzHdWKGeMUsfo9dbfnHVtm+imrUT9F93tKSa06WZwQAKxV7WJhmqZGjhypCy64QJ07dz7teqNHj9aoUaOquxmgRigqOKg1376p1psnKEUHji0zg7W2+QC1vfoxpUYy2ygASOdxKmTYsGH6/vvv9csvvygqKuq0653qiEV0dDSnQuAR9u7MVc53r6rzrqmqa5RKkvapvnJa367Yax5WWL2GFicEAPdw6amQBx98UN98840WLFhwxlIhSYGBgQoMDKzOZgDL5KxerINz/62uh+co/fgzKLb5RGlP3L3qetV9SgsKsToiANRIDhUL0zT14IMPatq0aZo3b55atWrlqlyA29krK7V6/lfyWfK2utiy1EY6PjFYF5WlDlf8xQMU4+trdUwAqNEcKhbDhg3TxIkT9fXXX6tu3brKz8+XJIWHhys4ONglAQFXKz1SrFUzx6pJ9ofqevyCzArTRyvDeqnuJQ8rlqnLAeCcOXSNhXGaWRfHjx+vIUOGnNOfwe2mqCn27sxVzvevq+OuqaqvIklSsRmsNU37K+bKRxQZw4PfAOAEl1xjcR6PvABqjA0ZP6lo/tvqWvjzsesnJO0yIrS97W2Ku2q40rggEwCqjblCUCvYbKVaNftjha/6QB0qNh5bePz6iaOJ96vrZbeomR//HADgfPGbFF5t786typn1H7XL+0LJOizp2AOtVtbrrfqXPqjYrhdYGxAAvAzFAl7HtNu1bslslf76ruKLFlad7tirBsppOVDtrxiu5CZnvk0aAFA9FAt4jeLCQ1oz631FrP9UsfatxxYaUrZ/Z5V2u1PxfW5XegDPVAEAV6JYwONtyV6qvT+9o877Zint+NMxS80ArW7YV40uHa7YzmkWJwSA2oNiAY90tPSIVs/5RHVWT1Cn8rVqLUmGlGc00852t6rTn4YqpUFjq2MCQK1DsYBH2b5ppXbOHaOOe75T8vFnT1SYPlpdp6cC0u5VbM+rFe3jY3FKAKi9KBao8Y6WHtGauZ8qePWniitbpRbHl+erkbbG3Kg2lw9V9+Y8Xh4AagKKBWqsresylf/zWHXc+52SVCxJqjQNrQlJkZLuUudeN6opz54AgBqF38qoUYoLDyl7zkcKWzdZHSvWq+Xx5flqpNwW16tVn/vVNbqtlREBAGdAsYDlTLtd65fNUdGij9T50I9KMWySpHLTV2tC0+STfKc6X3gdRycAwAPwmxqW2bszVzlzxylq6zR1MncdW3jizo7WA9S2z73q3jTa2pAAAIdQLOBWR0uPaO1Pk+S/erLiSjMUYRyb2K7EDNSa+pepbvoQdUzuw50dAOChKBZwOdNu14bMn1SweII6HZijRJUc++D4UzGLYwcqrvcdSqlbz9KcAIDzR7GAy+zKXa9t88YrevvX6mjurlqer0bKjbpGURffpdi2XSxMCABwNooFnKrg4D5t+OkT1dkwVbHla9Ts+PISM1DZ4RcpKOUOxaZfpaa+vpbmBAC4BsUC5812tETZ87+QVn2uuOLFSjEqJEl201B2UDcdjR2gTpcOUhKnOgDA61EsUC2VFRVat3iGSjInq+OheequI8c+MKRcnxjtiblGrS67U52j2lgbFADgVhQLnDPTbtemFb/owJKJarNntjrrYNVne9VAWyKvUETPwWrdOVU8YBsAaieKBc5qS/ZS7fl1oqJ3zVR7M79qeYFCtaHBpQpNvEUdU/sqggdYAUCtx54Ap7R94wrt/GWSmu6Yodb27cemJdexizDXhfWQb/xN6nRhf6UEhViaEwBQs1AsUGXH5jXa8ctENd4+Q23suVWziJaZflobmiJ73A3q1GuAEuuEW5oTAFBzUSxque2bVmrnr5MUkTdLbSpzFXV8ebnpq+zgBJV17K/2vW5W9/qNLM0JAPAMFItaxrTbtXX9cuUvnqImO35Qa/vWqiMTFaaP1gV1U2n7a9Xh4lvUtWETS7MCADwPxaIWMO12bV71q/Yv/VLNd89RK3Nn1V0b5aav1gUnyNb+arW7aKC6NGpqaVYAgGejWHipyooKrV/6g4qypqnlvp/VTvvU7vhnZaafskOSVNa+nzpcdJPiOTIBAHASioUXOVpSrPW/fqOytd+q3eFfFKfCqs9KzECtr5Mqe6er1eHCG9UtvIGFSQEA3opi4eEO7dutTb9Mld+mmep4ZJm6GbaqzwoUqo3hF8q/87Xq2PMaJYTUsTApAKA2oFh4oLxNK7VzyVcK2zZXHcrWKsUwj31gSPlqrG2Neym067XqkNJXyQGB1oYFANQqFAsPUF5m08Zlc1W06ls13ztf0eYuRZ/40JByfFtrb7NL1TjperXpkq6mPj5WxgUA1GIUixrq0L7d2vzbdPlsnq12RUsVd2KSL0llpq82BHVVSeu+ikm/QW1atBNTfQEAagKKRQ1hr7Rry5pF2rf8W9XfOU/ty9cr+cQpDkmHFKbN4T3k2+kKtUu/Rl24+BIAUANRLCxUcGi/chZ/q4oNP6jV4UVqq0Nqe+JDQ8rxbaW9TXupftd+apdwiZKZ5AsAUMOxp3Ije2WlclYv0v4V36nezgVqV7ZOCYa96vMSM1AbQhNV1uoyxaT1V5votpziAAB4FIqFi+3Pz1Pukm+lnJ/UunCp2qmg6kFVMqRtPlHa3fgChcZdqfYpfdSd2UIBAB6MYuFkR0uPaNOyuSpe94Mi9v6mNpVb9Pvpu46YQdoYmqCyVpcpOvlqxbTsoBjL0gIA4FwUi/Nkr6xUbvYy7Vs5SyE7Fqh96Sp1McpPWmezbxvta9JTYXF/Uruky9Q9MMiitAAAuBbFohp2bd2gHZkz5bt1vloVZaiNCv97LYQh7VN9bQ1PkdHmUrVK7ae2TaL+e1EmAABejGJxDvbn52lr5izZN89T88PL1Nzco2a/+7zEDNSm4HiVRl+kpglXKqZDghrzkCoAQC1EsTiFggN7lJPxg8o3/awmB5eppX37SddJVJg+2hzQUYea9lB4XB+17X6xunJ6AwAAioUkFRzcpy0Zs2XbvECNDmSodcUWJfzu4VTSsWdK7GuUpuAOF6tNUl91DKtvUVoAAGouh4vFggUL9MorrygzM1O7d+/WtGnT1L9/fxdEc51D+3Yrd/lcleUsUOMDGWpVkavuvy8ShrTNJ1r5DZLl37aXWiX2VZvGkTxTAgCAs3C4WBw5ckRdu3bVnXfeqRtuuMEVmZxu785cbV8xV5W5v6nJwQy1tG/XSccbjj9PIr9+kvxaX6iYxMsV07QFt4ECAOAgh4vFFVdcoSuuuMIVWZzCtNuVt3mV8lfPk5G3WM0Klqu5uUcR/7PeVp9o7amfKN/WF6hlYl+KBAAATuDyayxsNptsNlvV+8LCQudv42iJsqa+qsBdSxRzZJVaqFAtfvd5pWko16+19jdMVGCbCxTTvbdaRjRXS6cnAQCgdnN5sRg9erRGjRrl0m0EBASpw6axqq8iSZLN9FdOQAcVNE5UaLuL1CrhUrUNb8CzJAAAcDGXF4unn35aI0eOrHpfWFio6Ohop27D8PHR+pa3yfDxV72OF6lVfE/FMucGAABu5/JiERgYqMDAQFdvRulDXnT5NgAAwJnxeEgAAOA0Dh+xKC4u1ubNm6ve5+bmasWKFWrQoIFatGhxhm8CAABv53CxyMjI0CWXXFL1/sT1E4MHD9ZHH33ktGAAAMDzOFwsLr74YpmmefYVAQBArcM1FgAAwGkoFgAAwGkoFgAAwGkoFgAAwGkoFgAAwGkoFgAAwGkoFgAAwGkoFgAAwGkoFgAAwGlcPrvp/zrx1M7CwkJ3bxoAAFTTif322Z6+7fZiUVRUJEmKjo5296YBAMB5KioqUnh4+Gk/N0w3T/xht9u1a9cu1a1bV4ZhOO3PLSwsVHR0tPLy8hQWFua0PxcnY5zdh7F2D8bZPRhn93DlOJumqaKiIjVr1kw+Pqe/ksLtRyx8fHwUFRXlsj8/LCyMH1o3YJzdh7F2D8bZPRhn93DVOJ/pSMUJXLwJAACchmIBAACcxmuKRWBgoJ577jkFBgZaHcWrMc7uw1i7B+PsHoyze9SEcXb7xZsAAMB7ec0RCwAAYD2KBQAAcBqKBQAAcBqKBQAAcBqPKhbvvPOOWrVqpaCgICUmJmrhwoVnXH/+/PlKTExUUFCQWrdurXfffddNST2bI+P81VdfqU+fPmrcuLHCwsKUnp6u2bNnuzGt53L05/mEX3/9VX5+furWrZtrA3oRR8faZrPpmWeeUUxMjAIDA9WmTRt9+OGHbkrruRwd588++0xdu3ZVSEiIIiMjdeedd+rAgQNuSuuZFixYoKuvvlrNmjWTYRiaPn36Wb/j9n2h6SEmT55s+vv7m+PGjTOzs7PNESNGmKGhoea2bdtOuf6WLVvMkJAQc8SIEWZ2drY5btw409/f3/zyyy/dnNyzODrOI0aMMF966SVz6dKl5saNG82nn37a9Pf3N5cvX+7m5J7F0XE+4fDhw2br1q3Nyy+/3Ozatat7wnq46oz1NddcY6ampppz5swxc3NzzSVLlpi//vqrG1N7HkfHeeHChaaPj4/5xhtvmFu2bDEXLlxoxsXFmf3793dzcs8yY8YM85lnnjGnTp1qSjKnTZt2xvWt2Bd6TLFISUkxhw4detKyjh07mk899dQp13/iiSfMjh07nrTs/vvvN9PS0lyW0Rs4Os6nEhsba44aNcrZ0bxKdcd54MCB5rPPPms+99xzFItz5OhYz5w50wwPDzcPHDjgjnhew9FxfuWVV8zWrVuftOzNN980o6KiXJbR25xLsbBiX+gRp0LKysqUmZmpyy+//KTll19+uX777bdTfmfRokV/WL9v377KyMhQeXm5y7J6suqM8/+y2+0qKipSgwYNXBHRK1R3nMePH6+cnBw999xzro7oNaoz1t98842SkpL08ssvq3nz5mrfvr0ee+wxlZaWuiOyR6rOOPfo0UM7duzQjBkzZJqm9uzZoy+//FJXXXWVOyLXGlbsC90+CVl17N+/X5WVlWrSpMlJy5s0aaL8/PxTfic/P/+U61dUVGj//v2KjIx0WV5PVZ1x/l+vvfaajhw5optuuskVEb1CdcZ506ZNeuqpp7Rw4UL5+XnEP9saoTpjvWXLFv3yyy8KCgrStGnTtH//fj3wwAM6ePAg11mcRnXGuUePHvrss880cOBAHT16VBUVFbrmmmv01ltvuSNyrWHFvtAjjlic8L/TrJumecap10+1/qmW42SOjvMJkyZN0vPPP68pU6YoIiLCVfG8xrmOc2VlpW699VaNGjVK7du3d1c8r+LIz7TdbpdhGPrss8+UkpKiK6+8Uv/617/00UcfcdTiLBwZ5+zsbD300EP661//qszMTM2aNUu5ubkaOnSoO6LWKu7eF3rE//o0atRIvr6+f2i+e/fu/UMTO6Fp06anXN/Pz08NGzZ0WVZPVp1xPmHKlCm6++679cUXX6h3796ujOnxHB3noqIiZWRkKCsrS8OHD5d0bOdnmqb8/Pz0ww8/6NJLL3VLdk9TnZ/pyMhINW/e/KTpoTt16iTTNLVjxw61a9fOpZk9UXXGefTo0erZs6cef/xxSVJ8fLxCQ0N14YUX6oUXXuCospNYsS/0iCMWAQEBSkxM1Jw5c05aPmfOHPXo0eOU30lPT//D+j/88IOSkpLk7+/vsqyerDrjLB07UjFkyBBNnDiR86PnwNFxDgsL0+rVq7VixYqq19ChQ9WhQwetWLFCqamp7orucarzM92zZ0/t2rVLxcXFVcs2btwoHx8fRUVFuTSvp6rOOJeUlMjH5+RdkK+vr6T//h81zp8l+0KXXRbqZCduZfrggw/M7Oxs8+GHHzZDQ0PNrVu3mqZpmk899ZR5++23V61/4habRx55xMzOzjY/+OADbjc9B46O88SJE00/Pz/z7bffNnfv3l31Onz4sFV/BY/g6Dj/L+4KOXeOjnVRUZEZFRVl3njjjebatWvN+fPnm+3atTPvueceq/4KHsHRcR4/frzp5+dnvvPOO2ZOTo75yy+/mElJSWZKSopVfwWPUFRUZGZlZZlZWVmmJPNf//qXmZWVVXVbb03YF3pMsTBN03z77bfNmJgYMyAgwExISDDnz59f9dngwYPNXr16nbT+vHnzzO7du5sBAQFmy5YtzTFjxrg5sWdyZJx79eplSvrDa/Dgwe4P7mEc/Xn+PYqFYxwd63Xr1pm9e/c2g4ODzaioKHPkyJFmSUmJm1N7HkfH+c033zRjY2PN4OBgMzIy0hw0aJC5Y8cON6f2LD///PMZf+fWhH0h06YDAACn8YhrLAAAgGegWAAAAKehWAAAAKehWAAAAKehWAAAAKehWAAAAKehWAAAAKehWAAAAKehWAAAAKehWAAAAKehWAAAAKehWAAAAKf5fygbmzssWE5TAAAAAElFTkSuQmCC", "text/plain": [ "Figure(PyObject
)" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "trapezoidal();" ] }, { "cell_type": "markdown", "id": "dbc721b7", "metadata": {}, "source": [ "So we see the solution matches up extremely well." ] }, { "cell_type": "markdown", "id": "45c1b6ec", "metadata": {}, "source": [ "### Question 4\n", "\n", "Consider now the problem of $f''(x) = -4f(x)$. Introducing a new variable $v(x) = f'(x)$ rewrite this as a system of first order equations." ] }, { "cell_type": "markdown", "id": "69d6583a", "metadata": {}, "source": [ "$$ v'(x) = -4f(x), f'(x) = v(x) $$\n", "\n", "$$ \\underbrace{\\begin{bmatrix}\n", " v(x) \\\\\n", " f(x)\n", "\\end{bmatrix}}_{u'}' = \\underbrace{\\begin{bmatrix}\n", " 0 &-4 \\\\\n", " 1 & 0\n", " \\end{bmatrix}}_{A} \\underbrace{\\begin{bmatrix}\n", " v(x) \\\\ \n", " f(x)\n", " \\end{bmatrix}}_{u} $$" ] }, { "cell_type": "markdown", "id": "83e7c985", "metadata": {}, "source": [ "### Question 5\n", "\n", "Code up the above using the trapezoidal rule. Set an initial condition of $u(0) = 0, u'(0) = 2$." ] }, { "cell_type": "markdown", "id": "edfd4d7d", "metadata": {}, "source": [ "Plugging in to the equation above gives the scheme:\n", "\n", "$$ u(x) = u_0 + \\frac{x-x_0}{2}\\big( Au(0) + Au(x) \\big) $$\n", "$$ u(x) = u_0 +\\frac{x-x_0}{2}Au(0) + \\frac{x-x_0}{2}Au(x) $$\n", "$$ (I-\\frac{x-x_0}{2}A)u(x) = (I+\\frac{x-x_0}{2}A)u_0 $$\n", "$$ u(x) = \\big(I-\\frac{x-x_0}{2}A\\big)^{-1}\\big(I+\\frac{x-x_0}{2}A\\big)u_0$$" ] }, { "cell_type": "code", "execution_count": 5, "id": "12f56490", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "trapezoidalsystem (generic function with 1 method)" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using LinearAlgebra\n", "\n", "function trapezoidalsystem()\n", " \n", " dx = 1e-2;\n", " N = round( 1/dx );\n", " xvals = [0.0];\n", " fvals = [ [2.0,0.0] ];\n", " x = 0.0;\n", " \n", " A = [ [0,1] [-4,0] ]\n", " \n", " for ii = 1:N\n", " x += dx;\n", " push!( xvals, x );\n", " push!( fvals, ( I-dx*A*0.5 ) \\ ( (I+dx*A*0.5)*fvals[end] ) );\n", " end\n", " \n", " plot( xvals, [ fval[2] for fval in fvals ] )\n", " plot( xvals, [ sin(2*x) for x in xvals] )\n", " \n", "end" ] }, { "cell_type": "code", "execution_count": 6, "id": "6722f5d4", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAA9hAAAPYQGoP6dpAABAw0lEQVR4nO3deVxVdeL/8dflsgkJKiiCIOKOu4LgktO00WhT2VRamVu2MNWYWk06zrfFacZfy1RjpZa5ZKlZtozNWGkz424uiDuuqICACsoiyHbv+f1hMUNqchE4917ez8fjPh7dD+dw3/eEnDfnnPs5FsMwDERERERM4mF2ABEREWnYVEZERETEVCojIiIiYiqVERERETGVyoiIiIiYSmVERERETKUyIiIiIqZSGRERERFTeZodoDrsdjuZmZk0btwYi8VidhwRERGpBsMwKCwsJCwsDA+Pyx//cIkykpmZSUREhNkxREREpAbS09MJDw+/7Nddoow0btwYuPBmAgICTE4jIiIi1VFQUEBERETlfvxyXKKM/HhqJiAgQGVERETExVzpEgtdwCoiIiKmUhkRERERU6mMiIiIiKlURkRERMRUKiMiIiJiKpURERERMZXKiIiIiJhKZURERERMpTIiIiIipnK4jKxdu5bbbruNsLAwLBYLX3755RXXWbNmDTExMfj6+tK2bVtmz55dk6wiIiLihhwuI0VFRfTs2ZO33367WssfPXqUIUOGMGjQIJKTk/nDH/7A+PHj+eyzzxwOKyIiIu7H4XvTDB48mMGDB1d7+dmzZ9O6dWvefPNNAKKjo9m2bRuvvfYad911l6MvLyIiIm6mzm+Ut2nTJhISEqqM3XLLLcydO5fy8nK8vLwuWqe0tJTS0tLK5wUFBXUdU0REakFFeRl5p7MozDtFSeEZSs+dpaIoD9v5fIzy8xgVJVBRiqWiFOzlF38DDy8MTx+wemPx9MHi7Ye1USCefk3w9m+CzzVNCQhqSZPgULy8fer/DUqdqPMykp2dTUhISJWxkJAQKioqyMnJITQ09KJ1pk+fzosvvljX0URExAF2m42TGYfJTT9A8enj2M6mYS3IwO98JteUnyHQnkegUUiwxSC4HvLkcQ35Hk0559mM4kah2Bq3wqNpaxo1j6RZeDQtW7fH6ukSN6dv8Orl/9JPbx1sGMYlx380ZcoUJk2aVPm8oKCAiIiIugsoIiKVystKOXFkNzlHtlOeuQef/KM0PX+cUFsmoZZyLv4T8n9YwGZYKLA0psjiz3nrNZRa/SnzbIzd6ovd6oNh9cHw9MXw8IT/3Q8YBhZ7BZaKEiz2Miy2UqwVJXhVFOJjK6KR7Rx+RhFNjAKsFoMmnKOJ/RyUpUPZTsgHMv777coMTzKsoZzxjaAksB2eod1oFtWb8I498PFpVEdbT2qizstIy5Ytyc7OrjJ26tQpPD09CQoKuuQ6Pj4++Pjo8JuISF0rKy3heMpWzhzchCUzmaDC/URUpNHGUkGbny5subCDz/YIIc+nJecbhWEPCMezaWt8m4Xh3yyUwOataBLUkqaenjSto8x2m40zudkU5GZxLjeT87kZ2M5mYCnMwLcokyZl2bS0ZeNjKSfSnk5kcToUb4QsYDuUG1aOWsPJadwZW8veNOkQT5uu8fg28q+jxHIldV5G+vfvz1dffVVlbOXKlcTGxl7yehEREak7Z3OyObr9O0qPrKdpbjJR5UfoYPnJtRsWKDJ8SfeOIr9xR4zgjjRq2YngNl0JiehAay8vWpsTHwAPq5VmLVrRrEWryy5jq6ggM+MIuWn7KM7cDzkHCcg/QKvyowRYiomyHycq/zjkfwsHoPwrK4c8o8ht1gfvtgNo3ftGglua+S4bFofLyLlz5zh8+HDl86NHj7Jjxw6aNWtG69atmTJlCidOnGDhwoUAJCYm8vbbbzNp0iQefvhhNm3axNy5c1myZEntvQsREbmk/LM5pG5ZQdnBfxFyNok29vSqRywskI8/x32jKQrqgU/rGEI6xBIa2YHOVqtZsa+a1dOTsDadCGvTqcq4YbdzMuMI2Qe3cf74NvxydhFxfj9NLQV0sB2mw+nDcPoT2AwZllAym8ZibX89bfsOoWnznz1BJVfBYvx4AUc1rV69muuvv/6i8dGjR7NgwQLGjBnDsWPHWL16deXX1qxZw8SJE9m7dy9hYWE8++yzJCYmVvs1CwoKCAwMJD8/n4CAAEfiiog0KLaKCg4lryZv5z9omr2R9uUHsVqq/po/5hHBySZ98IjsT2i3X9AqKhqLR8OdkPtCQTnMid1rqDi6keZnt9Om4jgeP9luh63tON2iPwE9fk3H2Bvx8vI2KbHrqO7+2+EyYgaVERGRyyvMy+XQpuXY9n9N+/xNNKXqdAjHPcLJbhaPV4cbaNP7Bpq1CDMpqevIP3uaY9v/xfmD/yHk9Cai7MerfL0Afw42jsdon0DHa39DYFDIZb5Tw6YyIiLixvJzT3Jw7VK8D35FdHES3hZb5dcu7CjjsLe9kdaxg2kZ0d7EpO4hJzuNY1v+CYf/RfuCTTThXOXXyg0r+317UtzuVtr9YjjBLfXpzx+pjIiIuJmCvFz2/2cRjfZ/QeeSnXj9TwFJs7TiRMh1BPS4jY59b9IphDpkq6jg0PbVnN35FS2z/lPlqInNsLDfpxvnOtxJp+sfoElwwz5iojIiIuIGykrOs3ftZxi7PqFr4UZ8/ueTL0c8ojgVfgth/YcRGR1jYsqGLePwbjI2LCUo/Rs6VByqHC8zrOz1i6Oi6910u/5eGvlfY2JKc6iMiIi4sNQ933N67ft0OvV1lVMCxzwiyIr4Na2uvZ/WHXqYmFAuJTvtIMfWfETzY8tpZztaOV5g+JESdDNNBo6lY+/rGswFwyojIiIupiAvl5SV7xN0YCntbUcqx0/RjCMht9B84CjadevXYHZkru54yjayNnxE6xP/JMw4VTl+1KM1J9veTadfPUrT4JYmJqx7KiMiIi7iyK6N5K6eSbfclfhZLtwktMywsrvxtXjGjKLboKG6x4oLs9tspHz/Nec3L6Bb/mp8fzjVVmp4savJDTS+9lE6xVzvliVTZURExImVlZaw69v5XLP7AzqXp1SOH/Vozcn2w+l08zhNsuWGCvJy2L9qAUH7F9Puf45+HbG2JafLaHoOeditpqVXGRERcUJ5p7PY/88ZtDu2hOacBS4cBdkVcB1+Ax4hOv4Wt/wLWX7CMDi0fTV562bT/ey/Ko+WnCGAA+H30OHWCQSHuv509CojIiJOJP3QLrK+eY0eOSsqdzynacrhyHtpP/gxmus+KA1Wfu5JUr6eRZvDH9GS08CFgrqzyU0E3fI0bbvEmZyw5lRGREScwKEd6yj87lV6Fa6tnF78sLUdZ3s8TM9fjcXbx9fkhOIsKsrL2P3dIhptf4/O5fsqx3c06ofPL58iOj7BxHQ1ozIiImISw25n36YVGGtfo1tpcuX4jkb98Br0JF36/UqnYuRnHdq+msJ/v16lxO7z6kpZ/wn0/OXdLvPzozIiIlLPDLudfRu+wmPty0SX7wWgwvAgOfAmgm55hrZdXfdwu5gj49AuMr9+hV65X+NtqQDgoLUjRf2fotcNw5y+lKiMiIjUE8NuZ++Gr/Bc+/8qD6+XGl4kN7+N1r+efNFt7EUclZN1nMPLX6ZH5rLKj38fsransN8ket94n9OWEpUREZF6sH/LSuzfTaNL2W7ghxLSYihRQ6cS0irK5HTibs6eOsGBL6bTI/OTylJywLMjZdf9kW7X3o7FYjE5YVUqIyIidSh1z/cU/vM5ep7fDPy3hLQd+kdatGpjbjhxe3mns9j/xV/ocWJpZSnZ490Lz5ufo3PfG01O918qIyIidSDz6H6yvvgDMQX/Ai5cE5LU7FYif/MiLSPamZxOGprck+kc/mwavU9+XnlNSbLfAJrdMZ3ITr3MDYfKiIhIrco/c5qUT56jT9Ynlb/0tzW+gRa3v6gb1onpstMOkv7FC/Q5swKrxaDC8GBb8B20H/YSwSHhpuVSGRERqQVlpSVs/+xVog/OIpAiAPb49MJ3yF9o33OgyelEqjq+P5mzy6fQq3gTAOeMRuxuM5Zew6bSyP+aes+jMiIicpV2/XspTde9QISRCVy4b0zBtc/R45d3Oe2nF0QA9m38J97/fo72FYcByKI5mfFT6XPL6Hr92VUZERGpofSDOzj7+dP0KNkKQA5NONJ9AjG3P46nl7fJ6USqx26zkfz1XMK3/T9CyAVgj3cPfH/9Cu179K+XDCojIiIOKio4y+7FU4nJ+hgvi40yw0pS6L10ve8lAgKbmR1PpEbOFxWy8+MX6ZW2AF9LOTbDwtbgO4ge8SqBzVrU6WurjIiIVJNht7Pj2w9otXkaLTgDwI5G8TT7zau07tDT5HQitSM77SCZn/6ePoX/AS7cIfhQz2fpe/tjeFjr5tSNyoiISDVkHN7N2WVP0r0k6cJzS0tyrp1GrxuHm5xMpG7s2/hP/L57ljb29AvPvbric/sbtOseX+uvpTIiIvIzykrOs33x8/Q+Pg8fSzmlhhfbI0bTe8Q0fBv5mx1PpE6Vl5Wwfemf6X74XfwspRfmy+nyLPHDJ9fq61R3/+1Zq68qIuIC9m9ZSaNvJtHPng4W2OUbS9O7/0b/9t3MjiZSL7y8fYkf+SdOZozmwMcT6VG4jqDO15qWR2VERBqMwrxcUj6cRFzulwDkEkhq7HPEDnlQH9WVBikkvD0hT39F2sEdtO/Yy7QcKiMi0iDs+vdSWq6dTNwPF6huaTKETiPfpG9QiMnJRMzX2sQiAiojIuLmCnJPcXDhE8TmfwtAhiWUvBtfJe7a20xOJiI/UhkREbe1Y9ViwjdMIZY87IaFzS3vpdeo1wg3YVpsEbk8lRERcTsFeTkcmP9b+uavBCDNoxXFg2fQv+9NJicTkUtRGRERt7J33d8J/tdE+pKLzbCwJWwEvUe+jK+fjoaIOCuVERFxCyXFhexcMJH4U58CF64NOXfr2zoaIuICVEZExOUd2bUBzy8fId6eAcDmoKF0GzOD8MaBJicTkepQGRERl2W32di65E/0PjQDb4uN0zQl87rXiL/+brOjiYgDVEZExCXlZB4j+4MxxJcmgwW2+11L2wfn0jO4pdnRRMRBKiMi4nJ2/msJrdc9QzcKOW94s6vbZOLumqhZVEVclMqIiLiMstISts97kn4nPwbgsLUdXsPmEt+pt8nJRORqqIyIiEs4kZpC8eKR9Ks4BMCmFsPp8+Cb+Pj6mZxMRK6WyoiIOL3tX8+nw/dTaGU5Tz7+pA58lf43jzA7lojUEpUREXFaZaUlJL//OPGnl4EF9nt1ocnIhfRu3cHsaCJSi1RGRMQpnUw7RN7CEcRXHABgU+goYse+hpe3j8nJRKS2qYyIiNPZvfozwlc/SScKKcCf1Gtfp/9N95odS0TqiMqIiDgNu83Glg8mE3d8Dh4Wg8PWdjQasZhebTubHU1E6pDKiIg4hYK8XI6+ez/9zn8PFtjc7A56PjwL30b+ZkcTkTqmMiIipjuesg2PT0bS08ik1PBiZ68XiL/zCbNjiUg9URkREVMlfz2fTt8/i5+llGyCKbxzAXG9BpkdS0TqkcqIiJjCbrOxZf5T9MuYDxbY492LsIcW06FFK7OjiUg9UxkRkXp3ruAMh2ePoF/xRgA2hdxH34dm4OnlbXIyETGDyoiI1KvM1H2UfTSMXvb0C9eH9PkT/e/4rdmxRMREKiMiUm/2rF9OxHeJhFHEKZpx5vZ5xMVcb3YsETGZyoiI1Iuty16j9+4/42mxc9CzI03GfkrnVm3MjiUiTkBlRETqlK2igm1zHiP+5FKwwNbGN9H9sYWaP0REKqmMiEidOVdwliOzhhN/fjMAmyIT6Td6OhYPD5OTiYgzURkRkTpxMu0QxQvupqf9GCWGF3vjX6H/kAfNjiUiTkhlRERq3ZFd6wn4/AGiOEsOTThzxwfE9Pml2bFExEmpjIhIrdr176W0X/M7/CylHPWIxHf0MjpGdjQ7log4sRqduJ05cyZRUVH4+voSExPDunXrfnb5RYsW0bNnT/z8/AgNDWXs2LHk5ubWKLCIOK8tn7xM1zWP4mcpZbdPb4LG/4dQFRERuQKHy8jSpUuZMGECU6dOJTk5mUGDBjF48GDS0tIuufz69esZNWoU48aNY+/evXz66ads3bqVhx566KrDi4hzsNtsfP/u48Tt+wtWi8GWJkPoNOlbApoEmR1NRFyAw2Xk9ddfZ9y4cTz00ENER0fz5ptvEhERwaxZsy65/Pfff0+bNm0YP348UVFRXHvttTz66KNs27btqsOLiPnKSkvY/rdh9Mv6CLjwiZm+4xfh7eNjcjIRcRUOlZGysjKSkpJISEioMp6QkMDGjRsvuc6AAQPIyMhgxYoVGIbByZMnWbZsGbfeemvNU4uIUyjMy+Xg678ituA7yg0rW3r9hf5jX9ZHd0XEIQ5dwJqTk4PNZiMkJKTKeEhICNnZ2ZdcZ8CAASxatIjhw4dTUlJCRUUFt99+O2+99dZlX6e0tJTS0tLK5wUFBY7EFJF6kJN5jIK5d9DNdowiw5fD188k7pd3mR1LRFxQjf58sVgsVZ4bhnHR2I/27dvH+PHjee6550hKSuKbb77h6NGjJCYmXvb7T58+ncDAwMpHRERETWKKSB1JP7iDivduoq3tGDk0IfPOz+ipIiIiNWQxDMOo7sJlZWX4+fnx6aefcuedd1aOP/nkk+zYsYM1a9ZctM7IkSMpKSnh008/rRxbv349gwYNIjMzk9DQ0IvWudSRkYiICPLz8wkICKj2mxOR2ndo+2qClz9AUwpJs4RhHfkFrdp2NjuWiDihgoICAgMDr7j/dujIiLe3NzExMaxatarK+KpVqxgwYMAl1ykuLsbjJ+ePrVYrcOGIyqX4+PgQEBBQ5SEi5tuz7kta/X0YTSnkoLUj/r/9TkVERK6aw6dpJk2axPvvv8+8efNISUlh4sSJpKWlVZ52mTJlCqNGjapc/rbbbuPzzz9n1qxZpKamsmHDBsaPH09cXBxhYWG1905EpE5tXzGPjt89+MMcIn0Im7CKoBatzI4lIm7A4RlYhw8fTm5uLtOmTSMrK4tu3bqxYsUKIiMjAcjKyqoy58iYMWMoLCzk7bff5qmnnqJJkybccMMNvPzyy7X3LkSkTm355BVi9/4FD4vBtmuup/sTi/Hx9TM7loi4CYeuGTFLdc85iUgtMwy+XziVfkffAeD74N/QN3EOVk/dSUJErqy6+2/9RhGRSzLsdjbPGU+/rA8B2Bj+EP0ffFVziIhIrVMZEZGL2Csq2DZzLP3OLAdgY7tJDBj5vMmpRMRdqYyISBXlZaXseute4gr/jd2wsLXH8wy4a6LZsUTEjamMiEil0pJiUmbcRUzxRsoNKzvjXiX+1nFmxxIRN6cyIiIAnC8q5PBbd9CrJIlSw4uUX7xN7I33mh1LRBoAlRERoajgLMfevp3uZbsoNnw4ctMceg26w+xYItJAqIyINHD5Z3PIfudWulbsp9BoRMaQhXSPT7jyiiIitURlRKQBy889yamZQ+hkO0w+/py842Oi+/zC7Fgi0sCojIg0UHmns8idPYQOtlTOEMDZu5bRsXu82bFEpAFSGRFpgM6czCD/3VtpZz9GDk0oHPY57brEmB1LRBoolRGRBiYnO51z7w0myp7OaZpSfN8XRHXqbXYsEWnAVEZEGpCcrOMUzRlCG3sGp2hGyYi/E9mhh9mxRKSBUxkRaSBystMomjOESHsG2QRTMXI5rdt1NTuWiIjKiEhDkJOdTtF7/y0itlFfEd62i9mxREQA0O03Rdxc7skL14hE2tM5SRAVI5fTSkVERJyIyoiIGztzMoNz7w6mjT2dUzSj7IHlhOvUjIg4GZURETeVdzqL/HeHEPlDESl94Csi2nczO5aIyEVURkTcUP7Z0+TOHkKU/TinaUrJiL+riIiI01IZEXEzhXm5nHxnCO1sqeQSSNG9X9BaH98VESemMiLiRooKznLinV/TseIgZ2lM/j3LaNNZE5qJiHNTGRFxE+eLCjn+9m10Lt9HAf7k3rmUtl3jzI4lInJFKiMibqC0pJjDb91Bl7LdFBqNyLp9Ce17DjQ7lohItaiMiLi4irJSUmbcRfeSJIoNH04MWUinPteZHUtEpNpURkRcmK2igh1v30+v4o2UGl6k3jSHzvEJZscSEXGIyoiIizLsdpJmjia24DvKDSspg96m26A7zI4lIuIwlRERF2TY7WyenUjcmX9gMyzsjHuVXjfda3YsEZEaURkRcUGbFzxLv1NLAUjq9Sdibx1nciIRkZpTGRFxMZuXvES/tPcA+L7Ts8Td+TuTE4mIXB2VEREXsu3Lt4g/8CoAm1o/Qr/7/mByIhGRq6cyIuIidqz8kN7J/wfA9y2G02/MyyYnEhGpHSojIi5gz/rldNkwAavFYGuTwcQ9OguLh/75ioh70G8zESd3eMdaolY9jLelgmT/a+n9+EI8rFazY4mI1BqVEREnln5wB0FfjsDfUsJe755EP/EJnl7eZscSEalVKiMiTur0iaN4Lb6bphRw2NqO1o9/iW8jf7NjiYjUOpURESdUkHuK4rm30ZLTpFnCaPrIchoHNjM7lohInVAZEXEy54sKyZx9O5H2dE7RDM/RXxIUEm52LBGROqMyIuJEKsrLOPj23XQuT6EAf84NW0pYm05mxxIRqVMqIyJOwrDbSZ45hp7nv6fE8OLE4AW07RJndiwRkTqnMiLiJDbPe4q+Z/+JzbCwb+CbRMcnmB1JRKReqIyIOIHNn7xMv4x5AGzr9n/0SXjA5EQiIvVHZUTEZMkrP6Tv3ukAbIp4hPh7njI5kYhI/VIZETHR/i2riN4wEQ+LweZmd9BvrO43IyINj8qIiEnSDu6g5Yox+FrK2eHXn5jfvq/7zYhIg6TffCImyMlOw3PJPTThHAc8O9HxcU3zLiINl8qISD0rKswjb85QwoxTZFhCCX7kC/z8A8yOJSJiGpURkXpUUV7GkZn30N52hDMEwAOfEdSildmxRERMpTIiUk8Mu53ts8bR4/wWzhvenL5tIeHtupodS0TEdCojIvVk80fPE3dmOXbDwv6Bb9Ap5nqzI4mIOAWVEZF6kPTPOfRLnQHAls6/p7cmNRMRqaQyIlLHUjZ/Q/ctkwHY1GI4/e77g8mJRESci8qISB1KP7ST0K/H4W2pYLvftcQ9MtPsSCIiTkdlRKSOnD2dhcfiYZVziUQ//jFWT0+zY4mIOB2VEZE6UFpSTPZ7d9HKyCbT0oKghz+nkX9js2OJiDgllRGRWmbY7eyZOZLo8r0UGH6UDfuY4JBws2OJiDgtlRGRWrZ5/u+JKfiOcsPK8Rtn0SY6xuxIIiJOTWVEpBZtWz6LfulzANje4//o/ouh5gYSEXEBKiMitWT/5m/pkfRHADaFjiT+rokmJxIRcQ01KiMzZ84kKioKX19fYmJiWLdu3c8uX1paytSpU4mMjMTHx4d27doxb968GgUWcUYnUlMI+fqhCx/h9R9E/EN/MzuSiIjLcPhzhkuXLmXChAnMnDmTgQMH8u677zJ48GD27dtH69atL7nOsGHDOHnyJHPnzqV9+/acOnWKioqKqw4v4gwK8nKp+OgemlLAIWt7Ov92ER5Wq9mxRERchsUwDMORFeLj4+nTpw+zZs2qHIuOjmbo0KFMnz79ouW/+eYb7r33XlJTU2nWrFmNQhYUFBAYGEh+fj4BAbrVujiPivIyUv76K7qXJHGKZvDwv2nRKsrsWCIiTqG6+2+HTtOUlZWRlJREQkJClfGEhAQ2btx4yXWWL19ObGwsr7zyCq1ataJjx448/fTTnD9//rKvU1paSkFBQZWHiDNKejeR7iVJFBs+5A/9UEVERKQGHDpNk5OTg81mIyQkpMp4SEgI2dnZl1wnNTWV9evX4+vryxdffEFOTg6PPfYYZ86cuex1I9OnT+fFF190JJpIvdv8ycvE53wGwP4Br9Gn17UmJxIRcU01uoDVYrFUeW4YxkVjP7Lb7VgsFhYtWkRcXBxDhgzh9ddfZ8GCBZc9OjJlyhTy8/MrH+np6TWJKVJn9qz7OzF7/x8Am6KeoM8to0xOJCLiuhw6MhIcHIzVar3oKMipU6cuOlryo9DQUFq1akVgYGDlWHR0NIZhkJGRQYcOHS5ax8fHBx8fH0eiidSb9MO7af2v3+JpsbM1MIF+I/9kdiQREZfm0JERb29vYmJiWLVqVZXxVatWMWDAgEuuM3DgQDIzMzl37lzl2MGDB/Hw8CA8XFNki2vJP5uDsXg4ARRxwLMz3RMXYPHQdD0iIlfD4d+ikyZN4v3332fevHmkpKQwceJE0tLSSExMBC6cYhk16r+HrO+//36CgoIYO3Ys+/btY+3atTzzzDM8+OCDNGrUqPbeiUgdqygvI+3dYbS2nyCbYIIe+hTfRv5mxxIRcXkOzzMyfPhwcnNzmTZtGllZWXTr1o0VK1YQGRkJQFZWFmlpaZXLX3PNNaxatYrf/e53xMbGEhQUxLBhw3jppZdq712I1IOkOY8T/8MnZ87d9SHtW156Xh0REXGMw/OMmEHzjIjZtn72Bn13vwDA9n4z6POr0eYGEhFxAXUyz4hIQ7R/87f03HXhItVNrRNVREREapnKiMjPyE47RPOvH8bbYiPJ/xfEj754lmEREbk6KiMil3G+qJBzHwwniHyOWKOI/u1HeFj1T0ZEpLbpN6vIJRh2O/tmPUB72xHOEIDfqKX4XRN45RVFRMRhKiMil7B54R+JObeaMsNK9i3vERrZyexIIiJuS2VE5Cd2/vtj4o7OBCC521S69B9sciIREfemMiLyP9IO7qDtmgl4WAw2Bw0l/p6nzI4kIuL2VEZEflCYlwsf309jy3lSvLrS+5F3zY4kItIgqIyIAHabjdR376e1/QQnCSJ43Md4+/iaHUtEpEFQGREBtsx/mp7nv6fE8CL/9gU011TvIiL1RmVEGrzklR/SL2MeALt7T6Njn1+YnEhEpGFRGZEG7fj+7XTc8DQAm1oMp+/Qx0xOJCLS8KiMSINVkJeLx9IR+FtK2Ovdg9iH3zY7kohIg6QyIg2S3WYj9b0RRBiZZBNMy4c+xsvL2+xYIiINksqINEhbPphMr+JNlBpeFNyxgKAWrcyOJCLSYKmMSIOz819L6Jf23oX/7v0CHXsPMjmRiEjDpjIiDUrG4d20XTcRgO+D7yJu6BMmJxIREZURaTCKz+VRvngEjTlPilcX+jw80+xIIiKCyog0EIbdTsq7Y4iyHyeHJgSP1QyrIiLOQmVEGoQtS14ipvA/lBtWTg2eQ/OwSLMjiYjID1RGxO3t2/hPYg6+AUBS52foEp9gciIREflfKiPi1k6fOErIyt/iabGzNeBm4oc/a3YkERH5CZURcVtlpSXkLriPIPJJ9WhDt0fnY/HQj7yIiLPRb2ZxW8nvP07n8hQK8MN7xGIa+Tc2O5KIiFyCyoi4pW3/eI/408sAOHLt64S362pyIhERuRyVEXE7x/ZtpcvWPwKwqdUYet90n8mJRETk56iMiFspzD+D9dNR+FlK2e3Tm7ixfzU7koiIXIHKiLgNw27n0HujK+/E22rcIqyenmbHEhGRK1AZEbexeclL9ClaS5lhJe/Xc2imO/GKiLgElRFxC/s3ryTm4JsAbI/+PZ1jbzA3kIiIVJvKiLi83JMZNPv6UbwsNpIa30D8sN+bHUlERBygMiIuzVZRQda8EbTgDMc9wun8iCY2ExFxNfqtLS5t64Kn6Va6gyLDB2PYh/g3bmJ2JBERcZDKiLisnf/5hH4Z8wFI6ftn2nTuY3IiERGpCZURcUkn0w4RuWYiAN8H30nsrx82OZGIiNSUyoi4nLLSEvIWjqAJ5zhk7UDvh2aaHUlERK6Cyoi4nOS5v6NTxQEK8Mf/gY/w8fUzO5KIiFwFlRFxKcnfLCD+1CcApA58jbCoziYnEhGRq6UyIi7jxJE9tN80GYBNLUfQ6+b7TU4kIiK1QWVEXELJ+SJKFo+kseU8KV5diX3wDbMjiYhILVEZEZew6/3HaGdL5SwBBI1ZhJe3j9mRRESklqiMiNNL+ucc4nK/xG5YSP/lm7RoFWV2JBERqUUqI+LU0g/tovOWPwKwOXwMPX55l8mJRESktqmMiNMqKT5H+ccj8beUsNe7O33HvGJ2JBERqQMqI+K0dr3/W9rajnGGAFqM+QhPL2+zI4mISB1QGRGnlPSPOcSdWY7dsJBx/Qyah7UxO5KIiNQRlRFxOumHd9N564/XiYylx3V3mpxIRETqksqIOJWS80WULRmFv6WEfV7d6TvmZbMjiYhIHVMZEaeyc+4TlfOJBI/5UNeJiIg0ACoj4jS2fz2f+JzPAUi/7g3NJyIi0kCojIhTyDyaQofvpwCwKXQUPa6/2+REIiJSX1RGxHRlpSUULbpw35n9XtHEjn3N7EgiIlKPVEbEdNvnTaBDxSHy8afJyA913xkRkQZGZURMtfPfH9Pv5BIAjgx8lZatO5icSERE6pvKiJjm1ImjtF77NADfN7+HPjePMDmRiIiYQWVETGGrqOD0ByNpSiGHre3oPW6G2ZFERMQkNSojM2fOJCoqCl9fX2JiYli3bl211tuwYQOenp706tWrJi8rbmTrB5PpWrabIsMXn/s+wMfXz+xIIiJiEofLyNKlS5kwYQJTp04lOTmZQYMGMXjwYNLS0n52vfz8fEaNGsWNN95Y47DiHvZuXEHftPcBSImdRkT77iYnEhERM1kMwzAcWSE+Pp4+ffowa9asyrHo6GiGDh3K9OnTL7vevffeS4cOHbBarXz55Zfs2LGj2q9ZUFBAYGAg+fn5BAQEOBJXnEze6SzK3hlAC86wpckQ4iYsMTuSiIjUkeruvx06MlJWVkZSUhIJCQlVxhMSEti4ceNl15s/fz5Hjhzh+eefr9brlJaWUlBQUOUhrs+w2zk2fywtOEOaRyu6PTTb7EgiIuIEHCojOTk52Gw2QkJCqoyHhISQnZ19yXUOHTrE5MmTWbRoEZ6entV6nenTpxMYGFj5iIiIcCSmOKktS6fTq3gTpYYX5UPfx++aQLMjiYiIE6jRBawWi6XKc8MwLhoDsNls3H///bz44ot07Nix2t9/ypQp5OfnVz7S09NrElOcyJFdG+i9/3UAkjs/RbseA0xOJCIizqJ6hyp+EBwcjNVqvegoyKlTpy46WgJQWFjItm3bSE5O5oknngDAbrdjGAaenp6sXLmSG2644aL1fHx88PHRLJzuoqgwD+8vHsLbUkGy3wDihz9rdiQREXEiDh0Z8fb2JiYmhlWrVlUZX7VqFQMGXPyXbkBAALt372bHjh2Vj8TERDp16sSOHTuIj4+/uvTiEvbNTSTCyOQkQUQ9OB+Lh6a3ERGR/3LoyAjApEmTGDlyJLGxsfTv35/33nuPtLQ0EhMTgQunWE6cOMHChQvx8PCgW7duVdZv0aIFvr6+F42Le0r6xxz65n2NzbCQc8vbdA1uaXYkERFxMg6XkeHDh5Obm8u0adPIysqiW7durFixgsjISACysrKuOOeINAyZR1PouPX/wAJbIsbRf8AQsyOJiIgTcnieETNonhHXU15WSuorg+hUcYAUry50+P0aPL28zY4lIiL1qE7mGRGprqQPnqFTxQEK8KfJAwtURERE5LJURqTW7Vm/nLiMhQAciv8zoZGdTE4kIiLOTGVEatXZ01m0+O5JPCwGW5rdRszgsWZHEhERJ6cyIrXGsNs5Pv9BWnCG4x7hdHvwHbMjiYiIC1AZkVqzZdlr9CreSJnhSfmdmu5dRESqR2VEasWxfVvpufcVALZ3mkj77v1NTiQiIq5CZUSuWknxOYxl4/C1lLPTty9xw/9gdiQREXEhKiNy1XbOf5Io+3FyCaTVmPl4WPVjJSIi1ae9hlyVnf/+mPjTywA4cd3rBLeMMDmRiIi4GpURqbGc7DRar30GgO9bDKfH9XebnEhERFyRyojUiN1mI3PBWJpSQKpHG3qNfdPsSCIi4qJURqRGtiydTo+SbZQYXljvmYdvIz+zI4mIiItSGRGHpe7ZTJ8DbwCws8szREbHmJxIRERcmcqIOKSk+Bwen4/D21LBDr/+xN3zjNmRRETExamMiEN2zvsdbezp5NCE1mPmYvHQj5CIiFwd7Umk2nb9eynxOZ8DkHn9GzRr0crkRCIi4g5URqRack+m0+p/P8Z73W9MTiQiIu5CZUSuyLDbyVjwIEHkc9SjDb3GvmF2JBERcSMqI3JFWz59hZ7nt1BqeGHcPQffRv5mRxIRETeiMiI/63hKEj33vQZAcqeJtO0SZ3IiERFxNyojclmlJcVU/HA33l2+scQNn2J2JBERcUMqI3JZyQuepp3tKGcJIGzUPN2NV0RE6oT2LnJJe9YvJy5rMQBHB0wnOCzS5EQiIuKuVEbkIvlnTtP8uwl4WAw2N7udPgkPmB1JRETcmMqIVGHY7Rye9xAh5JJuCaPb2LfNjiQiIm5OZUSqSPrHu8ScW02F4cH522bh3zjQ7EgiIuLmVEakUuaxA3RKehGArW0eoWOfX5obSEREGgSVEQHAVlFB3qIHaWw5z36vLvR94E9mRxIRkQZCZUQA2LLoebqU76HI8KXxffPw9PI2O5KIiDQQKiPC4Z3riU2dBcDeXn+kVdtokxOJiEhDojLSwJUUn8Pr74l4WWxs9x9E3zseNzuSiIg0MCojDdzO+U8SaU8nhyZEjZmDxUM/EiIiUr+052nAdq3+jPjTywDIvO6vNG0eanIiERFpiFRGGqi8nGxCVz8FwPfBd9Hj+rtNTiQiIg2VykgDZNjtpC54mOac5bhHOD3H/s3sSCIi0oCpjDRASV/Nos+5tZQbVspun00j/8ZmRxIRkQZMZaSByTp+gM7bL0xoti3qETr0GmRyIhERaehURhoQW0UFZxeN4xrLefZ7RdN3xDSzI4mIiKiMNCRbP/4TXcp2U2z40Pi+uZplVUREnILKSAORumczfQ69DcCeHlNo1baryYlEREQuUBlpAEpLirF8/gjelgqS/frT984nzY4kIiJSSWWkAUj+4Gmi7Mc4QwARo9/XLKsiIuJUtFdyc/s2fU1c5mIAjg/8fwSHhJucSEREpCqVETdWmH+Gpit/h4fFYEuTIfS+eYTZkURERC6iMuLGUuY/TqhxmkxLCNFj3zE7joiIyCWpjLip5JUfEZe3ArthIe+WGTQObGZ2JBERkUtSGXFDOdnpRG6cAsDmsBF06fcrkxOJiIhcnsqImzHsdtIXPkIzCkj1aEOf0a+aHUlERORnqYy4mW1fvkXv4o2UGZ7wm/fw8fUzO5KIiMjPUhlxI5lH99Nl518A2N7uMdp2izc5kYiIyJWpjLgJe0UFeUsewt9SQopXV/re/7zZkURERKpFZcRNbPn4pcqb4AXc9z5WT0+zI4mIiFSLyogbOLpvK30OvQXAnu6TadW2i8mJREREqk9lxMWVlZZg/+zCTfB2Noqn728mmB1JRETEISojLm77h5NpZ0vlLI1pNUo3wRMREdejPZcLO7DtX/RNXwBAavyfCA5tbW4gERGRGqhRGZk5cyZRUVH4+voSExPDunXrLrvs559/zs0330zz5s0JCAigf//+fPvttzUOLBecLyrE759PYLUYbAu4iZjBY82OJCIiUiMOl5GlS5cyYcIEpk6dSnJyMoMGDWLw4MGkpaVdcvm1a9dy8803s2LFCpKSkrj++uu57bbbSE5OvurwDdmuBROIMDI5RTM6jJltdhwREZEasxiGYTiyQnx8PH369GHWrFmVY9HR0QwdOpTp06dX63t07dqV4cOH89xzz1Vr+YKCAgIDA8nPzycgIMCRuG5pz7q/0+1fowDYdf18elz3G5MTiYiIXKy6+2+HjoyUlZWRlJREQkJClfGEhAQ2btxYre9ht9spLCykWbPL30W2tLSUgoKCKg+5IP9sDs3/NRGAzUFDVURERMTlOVRGcnJysNlshISEVBkPCQkhOzu7Wt/jr3/9K0VFRQwbNuyyy0yfPp3AwMDKR0REhCMx3drBBY8RQi4ZlpZ0HzvD7DgiIiJXrUYXsFoslirPDcO4aOxSlixZwgsvvMDSpUtp0aLFZZebMmUK+fn5lY/09PSaxHQ7ySs/om/+t9gMC+cGv43fNYFmRxIREblqDs0ZHhwcjNVqvegoyKlTpy46WvJTS5cuZdy4cXz66afcdNNNP7usj48PPj4+jkRze2dOZhC5cQoAW8JG0j/uZpMTiYiI1A6Hjox4e3sTExPDqlWrqoyvWrWKAQMGXHa9JUuWMGbMGBYvXsytt95as6QNmGG3c3zhozSjgKMebegz+mWzI4mIiNQah++mNmnSJEaOHElsbCz9+/fnvffeIy0tjcTERODCKZYTJ06wcOFC4EIRGTVqFH/729/o169f5VGVRo0aERio0wzVkfSP2cQWrafcsGIfOgsfXz+zI4mIiNQah8vI8OHDyc3NZdq0aWRlZdGtWzdWrFhBZGQkAFlZWVXmHHn33XepqKjg8ccf5/HHH68cHz16NAsWLLj6d+DmTmYcpuP2PwGwrc3D9O9x+SNQIiIirsjheUbM0FDnGTHsdva8fCPdS7dz0LMjbZ/dgKeXt9mxREREqqVO5hmR+rVl2Wt0L91OieGFzz1zVERERMQtqYw4qROpe+m+9zUAdnSaQGSnXuYGEhERqSMqI07IVlFBwZKH8bOUste7B3HDp5gdSUREpM6ojDihrR+/RHT5XooMX5re/z4eVqvZkUREROqMyoiTOZayjT6H3gJgb4/JhLXpZHIiERGRuqUy4kTKy0opX/Yo3pYKdvr2pe+dT5odSUREpM6pjDiRbR/9Hx1sh8nHn7BR72Px0P8eERFxf9rbOYnDO9cTe/x9AA7FPE/zsDbmBhIREaknKiNOoLSkGM+//xYvi43t1/yCmFsfNjuSiIhIvVEZcQLJC5+ljT2NXAJpM2q2Ts+IiEiDor2eyfZv/Y6+Jz4EIG3AX2jWopXJiUREROqXyoiJzhcV4r/id1gtBlsDE+id8IDZkUREROqdyoiJdn0wkQgjk1M0o+PomWbHERERMYXKiEn2bvgH8ac+BSDrl68R2Ky5yYlERETMoTJignMFZ2n63UQANje7nZ6/vMvkRCIiIuZRGTHBvgXjCTNOkWlpQdcxM8yOIyIiYiqVkXq2a/VnxJ1ZDsCZm97gmoCmJicSERExl8pIPco/m0PL1U8D8H3ze+g28NcmJxIRETGfykg9OrjgMVpwhnRLGD3HvGF2HBEREaegMlJPdqxaTN/8b7EZFooGz6CRf2OzI4mIiDgFlZF6kJeTTfiGKQBsCRtB57ibTU4kIiLiPFRG6sGRDxIJJo9jHhH0HvWK2XFEREScispIHdv+9XxiCv9DheFB+W0z8W3kb3YkERERp6IyUodyT2YQtfk5ALZGjKFD71+YnEhERMT5qIzUEcNuJ23hozSlgCPWKGJGTjc7koiIiFNSGakjSf94j95F6ykzrDB0Nt4+vmZHEhERcUoqI3XgdOYxOmyfBkBSm4dp172fyYlEREScl8pILTPsdjI/fIRAijhkbU/siGlmRxIREXFqKiO1bOvf36Hn+c2UGZ543f0uXt4+ZkcSERFxaiojtSg7/TDRO/4MQFK7x2gTHWtyIhEREeenMlJLDLudU4seobHlPAc8OxN3//NmRxIREXEJKiO1ZMtnb9CjJIkSw4tGw97F6ulpdiQRERGXoDJSCzKPHaDbngvTvO/oOJ7WHXuZG0hERMSFqIxcJbvNxpklj+BvKSHFqyt9h//B7EgiIiIuRWXkKm1d9irdSndQbPjQePh7Oj0jIiLiIJWRq3AidS/d970OwO7oiYS372ZyIhEREdejMlJDdpuN/CWP4GcpZa93D/re83uzI4mIiLgklZEa2rL0L3Qp30Ox4UPT++bgYbWaHUlERMQlqYzUQPqhnfQ68DcAdnd9hrCoziYnEhERcV0qIw6yVVRQ9Mmj+FrK2e3Tm7i7nzI7koiIiEtTGXHQ1o9fonN5CueMRjQfMQeLhzahiIjI1dCe1AHH9yfT+9DbAOzrMZmWrTuYnEhERMT1qYxUU0V5GaXLHsHHUs4u3770vXO82ZFERETcgspINW1bPI2OFQcpwI+QB97V6RkREZFaoj1qNRzbt5U+qbMAONBrKiHh7UxOJCIi4j5URq6gvKyU8s8S8bZUsKNRP2Jvf8zsSCIiIm5FZeQKti1+ng62w+TjT6uROj0jIiJS27Rn/RmpezYTc/Q9AA71+T+ah7UxN5CIiIgbUhm5jPKyUowvEvG22Ej2G0jMrx81O5KIiIhbUhm5jG0f/ZF2tlTyuIaIUbN1ekZERKSOaA97CYd3biD2+NwL/933BYJbtjY5kYiIiPtSGfmJstISrH9/DC+Lje3+vyBm8DizI4mIiLg1lZGfSPpwClH2Y5wlgMhRs3R6RkREpI5pT/s/DiWvpW/6AgBS46YRFBJubiAREZEGQGXkB6UlxXh99RieFjtJ1/ySmCFjzY4kIiLSIKiM/GD7wmdpY08nl0Dajp5tdhwREZEGo0ZlZObMmURFReHr60tMTAzr1q372eXXrFlDTEwMvr6+tG3bltmznWtnf3D7auJOfAjA8f5/pmnzUJMTiYiINBwOl5GlS5cyYcIEpk6dSnJyMoMGDWLw4MGkpaVdcvmjR48yZMgQBg0aRHJyMn/4wx8YP348n3322VWHrw0l54vw+cfjWC0G2wJuos8tI82OJCIi0qBYDMMwHFkhPj6ePn36MGvWrMqx6Ohohg4dyvTp0y9a/tlnn2X58uWkpKRUjiUmJrJz5042bdpUrdcsKCggMDCQ/Px8AgICHIl7RZvefZz+WR+RQxO8freFwKCQWv3+IiIiDVV1998OHRkpKysjKSmJhISEKuMJCQls3Ljxkuts2rTpouVvueUWtm3bRnl5+SXXKS0tpaCgoMqjLuzf+h1xmYsAyBg4XUVERETEBA6VkZycHGw2GyEhVXfaISEhZGdnX3Kd7OzsSy5fUVFBTk7OJdeZPn06gYGBlY+IiAhHYlaLYbfj/fVTWC0GWwNvodfN99f6a4iIiMiV1egCVovFUuW5YRgXjV1p+UuN/2jKlCnk5+dXPtLT02sS82dZPDyw3jOPZL8BdBz9Tq1/fxEREakeT0cWDg4Oxmq1XnQU5NSpUxcd/fhRy5YtL7m8p6cnQUFBl1zHx8cHHx8fR6LVSGR0DJHRX9f564iIiMjlOXRkxNvbm5iYGFatWlVlfNWqVQwYMOCS6/Tv3/+i5VeuXElsbCxeXl4OxhURERF34/BpmkmTJvH+++8zb948UlJSmDhxImlpaSQmJgIXTrGMGjWqcvnExESOHz/OpEmTSElJYd68ecydO5enn3669t6FiIiIuCyHTtMADB8+nNzcXKZNm0ZWVhbdunVjxYoVREZGApCVlVVlzpGoqChWrFjBxIkTeeeddwgLC2PGjBncddddtfcuRERExGU5PM+IGepynhERERGpG3Uyz4iIiIhIbVMZEREREVOpjIiIiIipVEZERETEVCojIiIiYiqVERERETGVyoiIiIiYSmVERERETKUyIiIiIqZyeDp4M/w4SWxBQYHJSURERKS6ftxvX2myd5coI4WFhQBERESYnEREREQcVVhYSGBg4GW/7hL3prHb7WRmZtK4cWMsFkutfd+CggIiIiJIT0/XPW/qmLZ1/dB2rh/azvVD27l+1OV2NgyDwsJCwsLC8PC4/JUhLnFkxMPDg/Dw8Dr7/gEBAfpBryfa1vVD27l+aDvXD23n+lFX2/nnjoj8SBewioiIiKlURkRERMRUDbqM+Pj48Pzzz+Pj42N2FLenbV0/tJ3rh7Zz/dB2rh/OsJ1d4gJWERERcV8N+siIiIiImE9lREREREylMiIiIiKmUhkRERERU7l9GZk5cyZRUVH4+voSExPDunXrfnb5NWvWEBMTg6+vL23btmX27Nn1lNS1ObKdP//8c26++WaaN29OQEAA/fv359tvv63HtK7N0Z/pH23YsAFPT0969epVtwHdhKPbubS0lKlTpxIZGYmPjw/t2rVj3rx59ZTWdTm6nRctWkTPnj3x8/MjNDSUsWPHkpubW09pXdPatWu57bbbCAsLw2Kx8OWXX15xnXrfFxpu7OOPPza8vLyMOXPmGPv27TOefPJJw9/f3zh+/Pgll09NTTX8/PyMJ5980ti3b58xZ84cw8vLy1i2bFk9J3ctjm7nJ5980nj55ZeNLVu2GAcPHjSmTJlieHl5Gdu3b6/n5K7H0W39o7y8PKNt27ZGQkKC0bNnz/oJ68Jqsp1vv/12Iz4+3li1apVx9OhRY/PmzcaGDRvqMbXrcXQ7r1u3zvDw8DD+9re/Gampqca6deuMrl27GkOHDq3n5K5lxYoVxtSpU43PPvvMAIwvvvjiZ5c3Y1/o1mUkLi7OSExMrDLWuXNnY/LkyZdc/ve//73RuXPnKmOPPvqo0a9fvzrL6A4c3c6X0qVLF+PFF1+s7Whup6bbevjw4cYf//hH4/nnn1cZqQZHt/PXX39tBAYGGrm5ufURz204up1fffVVo23btlXGZsyYYYSHh9dZRndTnTJixr7QbU/TlJWVkZSUREJCQpXxhIQENm7ceMl1Nm3adNHyt9xyC9u2baO8vLzOsrqymmznn7Lb7RQWFtKsWbO6iOg2arqt58+fz5EjR3j++efrOqJbqMl2Xr58ObGxsbzyyiu0atWKjh078vTTT3P+/Pn6iOySarKdBwwYQEZGBitWrMAwDE6ePMmyZcu49dZb6yNyg2HGvtAlbpRXEzk5OdhsNkJCQqqMh4SEkJ2dfcl1srOzL7l8RUUFOTk5hIaG1lleV1WT7fxTf/3rXykqKmLYsGF1EdFt1GRbHzp0iMmTJ7Nu3To8Pd32n3utqsl2Tk1NZf369fj6+vLFF1+Qk5PDY489xpkzZ3TdyGXUZDsPGDCARYsWMXz4cEpKSqioqOD222/nrbfeqo/IDYYZ+0K3PTLyI4vFUuW5YRgXjV1p+UuNS1WObucfLVmyhBdeeIGlS5fSokWLuornVqq7rW02G/fffz8vvvgiHTt2rK94bsORn2m73Y7FYmHRokXExcUxZMgQXn/9dRYsWKCjI1fgyHbet28f48eP57nnniMpKYlvvvmGo0ePkpiYWB9RG5T63he67Z9KwcHBWK3Wixr2qVOnLmp8P2rZsuUll/f09CQoKKjOsrqymmznHy1dupRx48bx6aefctNNN9VlTLfg6LYuLCxk27ZtJCcn88QTTwAXdpqGYeDp6cnKlSu54YYb6iW7K6nJz3RoaCitWrWqcqv06OhoDMMgIyODDh061GlmV1ST7Tx9+nQGDhzIM888A0CPHj3w9/dn0KBBvPTSSzp6XUvM2Be67ZERb29vYmJiWLVqVZXxVatWMWDAgEuu079//4uWX7lyJbGxsXh5edVZVldWk+0MF46IjBkzhsWLF+t8bzU5uq0DAgLYvXs3O3bsqHwkJibSqVMnduzYQXx8fH1Fdyk1+ZkeOHAgmZmZnDt3rnLs4MGDeHh4EB4eXqd5XVVNtnNxcTEeHlV3W1arFfjvX+5y9UzZF9bZpbFO4MePjc2dO9fYt2+fMWHCBMPf3984duyYYRiGMXnyZGPkyJGVy//4caaJEyca+/btM+bOnauP9laDo9t58eLFhqenp/HOO+8YWVlZlY+8vDyz3oLLcHRb/5Q+TVM9jm7nwsJCIzw83Lj77ruNvXv3GmvWrDE6dOhgPPTQQ2a9BZfg6HaeP3++4enpacycOdM4cuSIsX79eiM2NtaIi4sz6y24hMLCQiM5OdlITk42AOP11183kpOTKz9C7Qz7QrcuI4ZhGO+8844RGRlpeHt7G3369DHWrFlT+bXRo0cb1113XZXlV69ebfTu3dvw9vY22rRpY8yaNaueE7smR7bzddddZwAXPUaPHl3/wV2Qoz/T/0tlpPoc3c4pKSnGTTfdZDRq1MgIDw83Jk2aZBQXF9dzatfj6HaeMWOG0aVLF6NRo0ZGaGioMWLECCMjI6OeU7uW//znPz/7O9cZ9oUWw9CxLRERETGP214zIiIiIq5BZURERERMpTIiIiIiplIZEREREVOpjIiIiIipVEZERETEVCojIiIiYiqVERERETGVyoiIiIiYSmVERERETKUyIiIiIqZSGRERERFT/X8zedW1PRSBAgAAAABJRU5ErkJggg==", "text/plain": [ "Figure(PyObject
)" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "trapezoidalsystem();" ] }, { "cell_type": "markdown", "id": "696cecd8", "metadata": {}, "source": [ "So the solution also matches well here." ] }, { "cell_type": "code", "execution_count": null, "id": "8a50446a", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Julia 1.6.5", "language": "julia", "name": "julia-1.6" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.6.5" } }, "nbformat": 4, "nbformat_minor": 5 }