{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\"Open" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# One Variable Equations" ] }, { "cell_type": "markdown", "metadata": { "cell_id": "00000-2303c462-3a38-4984-a18b-ad0743709cb9", "colab_type": "text", "deepnote_cell_type": "markdown", "id": "cx_2LAjUw7VH", "slideshow": { "slide_type": "slide" }, "tags": [] }, "source": [ "Throughout this section and the next ones we shall cover the topic of solutions to one variable equations. Many different problems in physics and astronomy require the use of complex expressions, even with implicit dependence of variables. When it is necessary to solve for one of those variable, an analytical approach is not usually the best solution, because of its complexity or even because it does not exist at all. Different approaches for dealing with this comprehend series expansions and numerical solutions. Among the most widely used numerical approaches are the Bisection or Binary-search method, fixed-point iteration, Newton's methods.\n", "\n", "For further details see for example Chap. 5 of Ref. ([1](../index.ipynb#Bibliography)) or Chap. 4 of Ref. ([3](../index.ipynb#Bibliography))\n" ] }, { "cell_type": "markdown", "metadata": { "cell_id": "00001-9b5c5089-cfa0-452e-9212-f426a9b2eab6", "colab_type": "text", "deepnote_cell_type": "markdown", "id": "sExxbMPEw7VM", "tags": [] }, "source": [ "- - -\n", "```{contents}\n", ":depth: 2\n", "```\n", "- - -" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Bibliography" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[1] [Kiusalaas, Numerical Methods in Engineering with Python](https://drive.google.com/file/d/0BxoOXsn2EUNIQUdFVkctR2xWRUk/view?usp=sharing)
\n", "[2] [Jensen, Computational_Physics](https://drive.google.com/file/d/0BxoOXsn2EUNIekRUMFVPYXVsSTg/view?usp=sharing). Companion repos: https://github.com/mhjensen [Web page](http://compphysics.github.io/ComputationalPhysics/doc/web/course)
\n", "[3] E. Ayres, [Computational Physics With Python](http://bit.ly/CPwPython)
" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "cell_id": "00002-49209fff-523b-4f09-a260-8afe344a3a8f", "colab": {}, "colab_type": "code", "deepnote_cell_type": "code", "id": "GD4z2RTZw7VN", "outputId": "2780c361-b150-411f-b8ac-88567710c34f", "slideshow": { "slide_type": "slide" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "In /home/usuario/.local/lib/python3.6/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle: \n", "The text.latex.preview rcparam was deprecated in Matplotlib 3.3 and will be removed two minor releases later.\n", "In /home/usuario/.local/lib/python3.6/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle: \n", "The mathtext.fallback_to_cm rcparam was deprecated in Matplotlib 3.3 and will be removed two minor releases later.\n", "In /home/usuario/.local/lib/python3.6/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle: Support for setting the 'mathtext.fallback_to_cm' rcParam is deprecated since 3.3 and will be removed two minor releases later; use 'mathtext.fallback : 'cm' instead.\n", "In /home/usuario/.local/lib/python3.6/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle: \n", "The validate_bool_maybe_none function was deprecated in Matplotlib 3.3 and will be removed two minor releases later.\n", "In /home/usuario/.local/lib/python3.6/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle: \n", "The savefig.jpeg_quality rcparam was deprecated in Matplotlib 3.3 and will be removed two minor releases later.\n", "In /home/usuario/.local/lib/python3.6/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle: \n", "The keymap.all_axes rcparam was deprecated in Matplotlib 3.3 and will be removed two minor releases later.\n", "In /home/usuario/.local/lib/python3.6/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle: \n", "The animation.avconv_path rcparam was deprecated in Matplotlib 3.3 and will be removed two minor releases later.\n", "In /home/usuario/.local/lib/python3.6/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle: \n", "The animation.avconv_args rcparam was deprecated in Matplotlib 3.3 and will be removed two minor releases later.\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Populating the interactive namespace from numpy and matplotlib\n" ] } ], "source": [ "%pylab inline\n", "rc('animation', html='jshtml')" ] }, { "cell_type": "markdown", "metadata": { "cell_id": "00003-8734dc8e-26e8-48ae-b4bf-33ecdf4eff54", "colab_type": "text", "deepnote_cell_type": "markdown", "id": "FY0-QkUnw7VS" }, "source": [ "JSAnimation import available at https://github.com/jakevdp/JSAnimation" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "cell_id": "00004-e90f9125-dcd9-4194-bce3-7cce10378e9a", "colab": {}, "colab_type": "code", "deepnote_cell_type": "code", "id": "5rHmuj6mw7VT" }, "outputs": [], "source": [ "#from JSAnimation import IPython_display\n", "from matplotlib import animation" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "cell_id": "00005-0a61ac2b-4f46-44ea-8318-57348dcc3944", "colab": {}, "colab_type": "code", "deepnote_cell_type": "code", "id": "zM--Q1_nw7VX" }, "outputs": [], "source": [ "import numpy as np\n", "from scipy import integrate\n", "from scipy import optimize" ] }, { "cell_type": "markdown", "metadata": { "cell_id": "00006-3bad2fa9-b3cf-46d5-9042-1f7768cc9854", "colab_type": "text", "deepnote_cell_type": "markdown", "id": "RNEo8n--w7Vc", "slideshow": { "slide_type": "slide" }, "tags": [] }, "source": [ "## Bisection Method" ] }, { "cell_type": "markdown", "metadata": { "cell_id": "00007-840214c8-76fc-4a9e-ae73-d63aa27a597e", "colab_type": "text", "deepnote_cell_type": "markdown", "id": "2ZgAUxxyw7Vf" }, "source": [ "The Bisection method exploits the [intermediate value theorem](http://en.wikipedia.org/wiki/Intermediate_value_theorem), where a continuous and differentiable function $f$ must have a zero between an interval $[a,b]$ such that $f(a)f(b)<0$, or equivalently, there must be a value $p\\in[a,b]$ such that $f(p)=0$. Below the algorithmm is stated explicitly." ] }, { "cell_type": "markdown", "metadata": { "cell_id": "00008-c9bd9d87-c169-4d2d-95f0-9cea2e1f4b1d", "colab_type": "text", "deepnote_cell_type": "markdown", "id": "tCrKSgUrw7Vh", "slideshow": { "slide_type": "slide" } }, "source": [ "### Steps BM" ] }, { "cell_type": "markdown", "metadata": { "cell_id": "00009-4c1446e1-f7d5-4e2d-918a-a6e3520086d1", "colab_type": "text", "deepnote_cell_type": "markdown", "id": "8uRIVSQOw7Vj" }, "source": [ "
\n", " \n", "
\n", "\n", "\n", "1. There must be selected two values $a$ and $b$ such that $f(a)f(b)<0$ and $p\\in[a,b]$ where $f(p)=0$. In other words, though we do not know the value of the root, we must know that there is at least one within the selected interval.\n", "\n", "2. To begin, it must be set $a_1=a$ and $b_1=b$.\n", "\n", "3. Calculate the mid-point $p_1$ as\n", "\n", " $$p_1 = a_1 + \\frac{b_1-a_1}{2} = \\frac{a_1+b_1}{2}$$\n", "\n", "4. Evaluate the function in $p_1$, if the [stop condition](#Stop-Condition) is true, go to step 6.\n", "\n", "5. If the [stop condition](#Stop-Condition) is not satisfied, then:\n", "\n", " 1. If $f(p_1)f(a_1) > 0$, $p\\in(p_1,b_1)$. Then set $a_2=p_1$ and $b_2=b_1$\n", "\n", " 2. If $f(p_1)f(a_1) < 0$, $p\\in(a_1,p_1)$. Then set $a_2=a_1$ and $b_2=p_1$\n", "\n", " 3. Go to step 3 using $p_2$, $a_2$ and $b_2$ instead of $p_1$, $a_1$ and $b_1$. For next iterations the index increases until the [stop condition](#Stop-Condition) is reached.\n", "\n", "6. The End!" ] }, { "cell_type": "markdown", "metadata": { "cell_id": "00010-843ff3e7-0e10-4488-9518-08a92129fcd5", "colab_type": "text", "deepnote_cell_type": "markdown", "id": "7uoUfZn9w7Vk", "slideshow": { "slide_type": "slide" }, "tags": [] }, "source": [ "### Stop condition BM" ] }, { "cell_type": "markdown", "metadata": { "cell_id": "00011-b90552b0-1873-490e-ba4a-8972b4b03191", "colab_type": "text", "deepnote_cell_type": "markdown", "id": "qtHv5IGIw7Vl", "slideshow": { "slide_type": "subslide" } }, "source": [ "There are several different stop conditions for this algorithm. The most used are stated below:\n", "\n", "* A fixed distance between the last two steps (absolute convergence): → `xtol`\n", "\n", " $$|p_i - p_{i-1}|<\\epsilon$$\n", "\n", "* A fixed relative distance between the last two steps (relative convergence): → `rtol`\n", "\n", " $$\\frac{|p_i - p_{i-1}|}{|p_i|}<\\epsilon\\ \\ \\ \\ \\ p_i \\neq 0$$\n", "\n", "* Function tolerance:\n", "\n", " $$f(p_i)< \\epsilon$$\n", "\n", "All these conditions should lead to a desired convergence expressed by the $\\epsilon$ value. However, the first and the third conditions present some problems when the function has a derivative very large or close to $0$ as evaluated in the root value. When the function is very inclined, the first condition fails as a convergence in the $x$ axis does not guarantee a convergence in the $y$ axis, so the found root $p$ may be far from the real value. When the function is very flat ($dF/dx\\rightarrow 0$), the third condition fails due to an analogous reason.\n", "\n", "A final stop condition which does not have mathematical motivation but a computational one, is a maximum number of allowed iterations. This condition should be used not only for this algorithm but for all iteration-based numerical methods. This condition guarantees a finite computing time and prevents undesired infinite bucles.\n", "\n", "* If $N>N_{max}$, stop!" ] }, { "cell_type": "markdown", "metadata": { "cell_id": "00012-ff797fa2-7ac6-4758-babd-e6465b4c2cf0", "colab_type": "text", "deepnote_cell_type": "markdown", "id": "nBQBHO23w7Vm", "tags": [] }, "source": [ "### Error analysis BM" ] }, { "cell_type": "markdown", "metadata": { "cell_id": "00013-bb303a1e-05d9-45b4-b2ac-16526adafb23", "colab_type": "text", "deepnote_cell_type": "markdown", "id": "2u-a2qPOw7Vn" }, "source": [ "If we suppose $f\\in C[a,b]$ and $f(a)f(b)<0$, the Bisection method generates a sequence of numbers $\\left\\{p_i\\right\\}_{i=1}^\\infty$ approximating a root $p$ of $f$ as:\n", "\n", "$$|p_i-p|\\leq \\frac{b-a}{2^n},\\ \\ \\ \\ \\ i\\geq 1$$\n", "\n", "From this, we can conclude the convergence rate of the method is\n", "\n", "$$p_i = p + \\mathcal{O}\\left( \\frac{1}{2^i} \\right)$$\n", "\n", "This expression allows us to estimate the maximum number of required iterations for achieving a desired precision. The next figure sketches the number of iterations required for some precision." ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "cell_id": "00014-56003021-94d2-44cf-8c2e-9c6836b21333", "colab": {}, "colab_type": "code", "deepnote_cell_type": "code", "id": "OsQu8x8Jw7Vo" }, "outputs": [], "source": [ "import time" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "cell_id": "00015-07db623f-b3fd-4fd8-98ec-8aac5df309a9", "colab": {}, "colab_type": "code", "deepnote_cell_type": "code", "id": "NRO_1bQGw7Vt", "outputId": "1816cac5-f362-471a-ef17-40cda939e068" }, "outputs": [ { "data": { "text/plain": [ "array([1, 2, 3])" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "i=np.arange( 1, 101)\n", "i[:3]" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "100" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "i[-1]" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "cell_id": "00017-6ca8417b-1bcd-4a50-8698-2c4ba327f63d", "colab": {}, "colab_type": "code", "deepnote_cell_type": "code", "id": "j02SkcI1w7V4", "outputId": "6845b6c1-d366-4b39-c200-612315cd0cee" }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAD4CAYAAAAKA1qZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAo/klEQVR4nO3dd3xUdfb/8ddJh0Q6RKqhBDB0EjokRClBiqgoIIuoKIIiSHbXRVd3cb+u6K4buiAiYEEiolJipCkk9BKlBEIXJHQEwYCUkM/vj4y/b75ZgikzuZO55/l4zEPuZ+bOnLPXnTd35npGjDEopZSyJy+rC1BKKWUdDQGllLIxDQGllLIxDQGllLIxDQGllLIxH6sLKIhKlSqZkJCQAu1z+fJlAgMDXVOQm9Ke7cOOfduxZyha3ykpKeeMMZVvdV+JCoGQkBC2bdtWoH3WrFlD586dXVOQm9Ke7cOOfduxZyha3yJyNK/79OMgpZSyMQ0BpZSyMQ0BpZSyMUu/ExCRQOAd4Dqwxhgzz8p6lFLKbpx+JiAis0XkjIik5lqPEZF9InJQRMY6lh8EFhpjngb6OLsWpZRSt+eKj4PmAjE5F0TEG5gG9ADCgIEiEgbUAI45HnbTBbUopZS6DXHFFFERCQESjDGNHdvtgHHGmO6O7ZccD00HLhhjEkQk3hgz4BbPNQwYBhAcHBweHx9foFoyMjIICgoqdC8lkfZsH3bs2449Q9H6jo6OTjHGRNzqvuL6TqA6//s3fsh+828DTAamikhPYOmtdjTGzARmAkRERJiCXie7Zs0aOnSK5I3ENJ7uVIdq5UoVovySxY7XUduxZ7Bn33bsGVzXt6VfDBtjLgNPuPp1dp+4RPyWY3y2LZ2xPRryaOtaeHmJq19WKaXcXnFdInocqJlju4ZjrVg0r1mO5S9E0qxmWV5ZlMqA9zbxw7nLxfXySinltoorBLYCoSJSW0T8gAHAkmJ6bQBqVSzNx0Pb8K+HmpJ28hIxE5N5N+kQmTezirMMpZRyK664RHQ+sBFoICLpIjLUGJMJjASWA2nAAmPMbme/dj5q45FWNVkVG0VU/cqM/3ovD7yzgbSTl4q7FKWUcgtO/07AGDMwj/VEINHZr1cYwWUCeHdwOIm7TvH3Jan0nrKOEZ3rMvKeevj7eFtdnlJKFRvbjo0QEXo2rcrKMVH0aV6NKd8epOfkdaQcvWB1aUopVWxsGwK/KR/oR9wjzZn7RCuuXMuk34wNvLZ0N1euZ1pdmlJKuZztQ+A3nRtUYUVsFIPb3sWc9UfoNiGZdQfOWV2WUkq5lIZADkH+Pvzj/sYseKYdft5e/OH9zby4cAcXr9ywujSllHIJDYFbaF27AomjOzGic10+/+44XSYksSz1lNVlKaWU02kI5CHA15u/xDRk8XMdqBzkz/CPU3hu3nec/eWa1aUppZTTaAj8jsbVy7J4ZAf+3L0BK9NO0yUuic9T0nHF4D2llCpuGgL54OvtxXPR9Ugc1YnQKkH88bMdDJmzlfQLV6wuTSmlikRDoADqVQliwTPteK1PI7YdOU/3Ccl8uPEIWVl6VqCUKpk0BArIy0sY0j6EFWMiaXlXef62eDf9Z27k0NkMq0tTSqkC0xAopBrlS/Phk615++Fm7D+dQY9Ja3lnzUEdSKeUKlE0BIpAROgXXoOVsZHc06AK/1q2j77vrGf3iYtWl6aUUvmiIeAEVe4IYMbgcKYPasmpi9foM3U9/16+l6s39GeTlVLuTUPAiXo0qcqq2EgeaFGdaasPcd/ktWw7ct7qspRSKk8aAk5WrrQfbz/cjA+fbM21G1k8/O5Gxi3ZzeVrOpBOKeV+NARcJLJ+ZVaMiWRIuxA+2Jg9kC5p/1mry1JKqf9DQ8CFAv19GNenEZ89044AXy+GzN7CHxfs4Ocr160uTSmlADcIAREJFJFtItLL6lpcJSKkAl+N6sRz0XVZtP04XeKS+XrXSavLUkqpwoeAiMwWkTMikpprPUZE9onIQREZm4+n+guwoLB1lBQBvt78uXtDlozsQHAZf0bM+47hH6Vw5tJVq0tTStlYUc4E5gIxORdExBuYBvQAwoCBIhImIk1EJCHXrYqIdAX2AGeKUEeJ0qhaWRY/14G/xDTk231n6BKXxGfbjulAOqWUJaQobz4iEgIkGGMaO7bbAeOMMd0d2y8BGGPG57H/P4FAsgPjV+ABY0xWrscMA4YBBAcHh8fHxxeoxoyMDIKCggq0T3E5dTmL2anX2H8hi0YVvXi8kT+VSxf9Ezp37tlV7Ngz2LNvO/YMRes7Ojo6xRgTccs7jTGFvgEhQGqO7X7ArBzbg4Gp+Xiex4Fev/e48PBwU1CrV68u8D7F6ebNLPPhxiMm7NWvzd2vfm1mrztsbt7MKtJzunvPrmDHno2xZ9927NmYovUNbDN5vK9a/sUwgDFmrjEmweo6rODlJQxuexcrYqNoFVKB15bu4eF3N3LwjA6kU0q5nrND4DhQM8d2Dcea+h3Vy5Vi7hOtiHukGYfOZnDfpLVMW32QGzqQTinlQs4Oga1AqIjUFhE/YACwxMmv4bFEhAdb1mDlmCi6Ngrm38v30WfqelKP60A6pZRrFOUS0fnARqCBiKSLyFBjTCYwElgOpAELjDG7nVOqfVS+w59pj7bk3cHh/JRxjfunreetZTqQTinlfD6F3dEYMzCP9UQgsdAVqf+ve6M7aVu7Im8kpjF9zSGWp57izYea0rp2BatLU0p5CLf4YljlrWxpX97q15SPh7bhRlYWj7y7kVcXpZKhA+mUUk6gIVBCdAytxPIXInmyQ20+3nyUbnFJrN5nm//GTinlIhoCJUhpPx/+1juMz0e0J9DfhyfmbCX20+1cuKwD6ZRShaMhUAK1rFWehFEdGXVPPZbsOEGXuCQSdp7Q0RNKqQLTECih/H28ie3WgKXPd6R6+VKM/OR7nvkohdM6kE4pVQAaAiXc3VXL8MWI9rx8X0OS9p+lS1wSSek39KxAKZUvGgIewMfbi2GRdVn+QiRhVcswJ/U6g2Zt5sefrlhdmlLKzWkIeJCQSoHMf7otQ8L82Jl+ke4Tk5m19jA3s/SsQCl1axoCHsbLS4iu5cvK2Eja1a3I61+l8dD0Dew//YvVpSml3JCGgIeqWrYU7w+JYNKA5hz96TI9J69l0qoDXM/UgXRKqf+lIeDBRIT7m1dnVWwUMY2rMmHVfvpMXceOYz9bXZpSyk1oCNhAxSB/pgxswXuPRXDhynUeeGc9bySm8et1HUinlN1pCNhI17BgVsZG0b9VTWYmH6bHpGQ2HvrJ6rKUUhbSELCZMgG+jH+wKZ881YYsAwPf28TLX+7i0tUbVpemlLKAhoBNta+XPZDuqY61id/yI93ikvkm7bTVZSmlipmGgI2V8vPmlV5hfPFsB8qW8mXoB9sYNf97fsq4ZnVpSqliYmkIiEgtEVkkIrNFZKyVtdhZ85rlWPp8R0bfG8rXqSfpOiGZxduP6+gJpWygKD8vOVtEzohIaq71GBHZJyIH8/HG3gRYaIx5EmhR2FpU0fn5eDGma30Snu9EzQqlGR2/nac+2MbJi79aXZpSyoWKciYwF4jJuSAi3sA0oAcQBgwUkTARaSIiCbluVYBNwFAR+RZYVoRalJM0uPMOvhjRnld63s36Q+foFpfMJ5t/JEtHTyjlkaQop/wiEgIkGGMaO7bbAeOMMd0d2y8BGGPG57H/n4AtxphkEVlojOl3i8cMA4YBBAcHh8fHxxeoxoyMDIKCggq0T0nnrJ7PXMliTuo10s5n0bCCF0808ic40D2/RrLjcQZ79m3HnqFofUdHR6cYYyJueacxptA3IARIzbHdD5iVY3swMPU2+zcGFgIzgLd/7/XCw8NNQa1evbrA+5R0zuw5KyvLxG85ahr/fZlp8EqimZl0yGTezHLa8zuLHY+zMfbs2449G1O0voFtJo/3VZ9CxYqTGGNSyQ4O5aZEhP6tahFVvwqvLErln4lpJOw8wb/6NaPBnXdYXZ5SqoicfW5/HKiZY7uGY02VcHeWDeC9x8KZMrAFxy78Sq8pa5mwcr8OpFOqhHN2CGwFQkWktoj4AQOAJU5+DWUREaF3s2qsio2iZ5OqTPrmAL2mrOX7Hy9YXZpSqpCKconofGAj0EBE0kVkqDEmExgJLAfSgAXGmN3OKVW5iwqBfkwc0ILZj0fwy9VMHpq+gdcT9nDleqbVpSmlCqjQ3wkYYwbmsZ4IJBa6IlVi3NMwmBVjKvDWsr3MWvcDK/ac5s0Hm9C+XiWrS1NK5ZN7Xu+nSow7Anx5vW8T4oe1xUvg0VmbGfv5Ti7+qgPplCoJNASUU7StU5FlL0TyTFQdFmw7RrcJSazcowPplHJ3GgLKaQJ8vXmpx90seq4D5Uv78fSH2xj5yXec04F0SrktDQHldE1rZA+k+2PX+qzYfZoucUl8+X26DqRTyg1pCCiX8PX24vl7Q/lqVEdqVwpkzKc7eHLuVk78rAPplHInGgLKpUKD72Dh8Pb8vXcYmw6fp9uEZD7adFQH0inlJjQElMt5ewlPdKjNijGRNK9ZjlcXpTJg5iYOn82wujSlbE9DQBWbmhVK89HQ1vyrX1P2nrpEj0lrmZF0iMybOnpCKatoCKhiJSI8ElGTVbFRdG5QmTe/3ssD72xgz4lLVpemlC1pCChLVCkTwIw/hPPOoJacvPgrfaau4z8r9nEt86bVpSllKxoCyjIiwn1NqrJyTBR9mlVjyrcH6Tl5HSlHdSCdUsVFQ0BZrnygH3H9mzP3iVb8ev0m/WZs4LWlu7l8TQfSKeVqGgLKbXRuUIXlYyIZ3PYu5qw/QveJyaw9cNbqspTyaBoCyq0E+fvwj/sbs+CZdvh5ezH4/S28uHAHF6/oQDqlXEFDQLml1rUrkDi6EyM61+Xz747TZUISy1JPWV2WUh5HQ0C5rQBfb/4S05DFz3WgcpA/wz9O4dl5KZz55arVpSnlMTQElNtrXL0si0d24M/dG7Aq7Qxd45L5PEUH0inlDMUaAiJSR0TeF5GFOdb6ish7IvKpiHQrznpUyeHr7cVz0fVIHNWJelWC+ONnOxgyZyvpF65YXZpSJVq+Q0BEZovIGRFJzbUeIyL7ROSgiIy93XMYYw4bY4bmWltkjHkaGA70L0jxyn7qVQnis2fa8VqfRmw7cp7uE5L5cOMRsvSsQKlCKchvDM8FpgIf/rYgIt7ANKArkA5sFZElgDcwPtf+Txpjztzm+V9xPJdSt+XlJQxpH8K9d1fh5S9T+dvi3YSW8+KuxhnUrRxkdXlKlShSkM9VRSQESDDGNHZstwPGGWO6O7ZfAjDG5A6A3M+z0BjTz/FnAd4EVhpjVt3iscOAYQDBwcHh8fHx+a4XICMjg6Age70x2KlnYwwbTmQyL+0a17OEvnV9ianti4+XWF1asbDTsf6NHXuGovUdHR2dYoyJuNV9BTkTuJXqwLEc2+lAm7weLCIVgX8CLUTkJUdYPA90AcqKSD1jzIyc+xhjZgIzASIiIkznzp0LVOCaNWso6D4lnd16jgYaLf+W5WfLsjD1FGmXS/HWQ01pXL2s1aW5nN2ONdizZ3Bd30UNgQIxxvxE9mf/OdcmA5OLsw7lecr5ezH9D+F8veskry7ezf3T1vNMZB1G3RtKgK+31eUp5baKenXQcaBmju0ajjWlLNGjSVVWxUbyQIvqvLPmEPdNXsu2I+etLkspt1XUENgKhIpIbRHxAwYAS4pellKFV660H28/3IwPn2zNtRtZPPzuRsYt0YF0St1KQS4RnQ9sBBqISLqIDDXGZAIjgeVAGrDAGLPbNaUqVTCR9SuzYkwkQ9qF8MHGI3SbkEzSfh1Ip1RO+f5OwBgzMI/1RCDRaRUp5USB/j6M69OI3s2q8uLCnQyZvYWHWtbg1V53U660n9XlKWU5HRuhbCH8rgp8NaoTI6PrsWj7cbrEJfP1rpNWl6WU5TQElG0E+Hrzp+4NWDKyA3eW9WfEvO8Y/lEKZy7pQDplXxoCynYaVSvLomc78JeYhny77wxd4pJYsO2YDqRTtqQhoGzJx9uLEZ3rsmx0JxreWYYXF+7ksdlbOHZeB9Ipe9EQULZWp3IQ8cPa8j99G/Pd0Qt0n5jMnPU/cDNLzwqUPWgIKNvz8hIGt72LFbFRtK5dgdeW7uHhGRs4cPoXq0tTyuU0BJRyqF6uFHMeb8WE/s04fO4yPSevY8o3B7hxM8vq0pRyGQ0BpXIQER5oUYNVsVF0bRTMf1bup/eUdexKv2h1aUq5hIaAUrdQKcifaY+25N3B4Zy/fJ2+76xn/NdpXL1x0+rSlHIqDQGlbqN7oztZGRtFv5Y1eDfpMD0mrWXz4Z+sLkspp9EQUOp3lC3ly1v9mvLx0DZkZmXRf+YmXlm0i1+u3rC6NKWKTENAqXzqGFqJ5S9E8mSH2szb/CPdJySzeu/tfjFVKfenIaBUAZT28+FvvcP4fER7Av19eGLuVsZ8up3zl69bXZpShaIhoFQhtKxVnoRRHRl1byhLd5yga1wSCTtP6OgJVeJoCChVSP4+3sR2rc/S5ztSvXwpRn7yPcM+SuG0DqRTJYiGgFJFdHfVMnwxoj0v39eQ5P1n6RKXRPyWH/WsQJUIxRYCIlJHRN4XkYU51rxE5J8iMkVEhhRXLUo5m4+3F8Mi67L8hUjCqpZh7Be7GDRrM0d/umx1aUrdVr5CQERmi8gZEUnNtR4jIvtE5KCIjL3dcxhjDhtjhuZavp/sH6e/AaQXpHCl3FFIpUDmP92Wfz7QmF3pF+k+MZlZaw/rQDrltvJ7JjAXiMm5ICLewDSgBxAGDBSRMBFpIiIJuW5V8njeBsAGY0wsMKJwLSjlXry8hEFt7mJFbCQd6lbi9a/SeGj6Bvad0oF0yv1Ifj+3FJEQIMEY09ix3Q4YZ4zp7th+CcAYM/53nmehMaaf489/AK4bYxaIyKfGmP63ePwwYBhAcHBweHx8fH57AyAjI4OgoKAC7VPSac/uwxjD5lM3mbfnGlcyoXddX3rV8cXHS5zy/O7atyvZsWcoWt/R0dEpxpiIW95pjMnXDQgBUnNs9wNm5dgeDEy9zf4VgRnAIeAlx1pp4H1gCvDc79UQHh5uCmr16tUF3qek057dz7lfrprnP/nO3PWXBNMtLsls//GCU57X3ft2BTv2bEzR+ga2mTzeV30KFSuFYIz5CRiea+0KkPt7AqU8TsUgfyYPbEGfZtX466JdPPDOeoZ2rE1s1waU8vO2ujxlY0W5Oug4UDPHdg3HmlIqD13CglkZG0X/VrV4b+0PxExKZuMhHUinrFOUENgKhIpIbRHxAwYAS5xTllKeq0yAL+MfbMInT7cBYOB7m3j5y11c0oF0ygL5vUR0PrARaCAi6SIy1BiTCYwElgNpwAJjzG7XlaqUZ2lftxLLRkcyLLIO8Vt+pFtcMt+knba6LGUz+fpOwBgzMI/1RCDRqRUpZSOl/Lx5+b676dmkKi8u3MnQD7bRp1k1/t47jIpB/laXp2xAx0Yo5Qaa1SzH0uc7MqZLfb5OPUnXCcks3n5cR08ol9MQUMpN+Pl4MbpLKF+N6kStCqUZHb+dpz7YxsmLv1pdmvJgGgJKuZn6wXfw+Yj2vNLzbtYfOkfXuGTmbT5Klo6eUC6gIaCUG/L2Ep7qVIcVL0TRtEZZ/vplKo/O2sSRczqQTjmXhoBSbqxWxdLMe6oNbz3UhN0nLtF9YjIzkw+ReTPL6tKUh9AQUMrNiQj9W9ViVWwUkfUr80biXh6avoG9py5ZXZryABoCSpUQwWUCmDk4nKmPtiD9wq/0mryOLw9c51rmTatLUyWYhoBSJYiI0KtpNVbFRtG7WTUWH7pBr8nr+P7HC1aXpkooDQGlSqDygX5M6N+cMeH+ZFzL5MHpG/ifhD1cuZ5pdWmqhNEQUKoEa1bZhxVjIhnUphbvr/uBmIlr2XDwnNVlqRJEQ0CpEu6OAF9e79uET4e1xdtLeHTWZsZ+vpOLv+pAOvX7NASU8hBt6lTk69GdeCaqDgu2HaNrXBIrdp+yuizl5jQElPIgAb7evNTjbhY914EKgX4M+yiFkZ98x7mMa1aXptyUhoBSHqhpjeyBdH/sWp8Vu0/TJS6JL79P14F06r9oCCjloXy9vXj+3lC+GtWROpUCGfPpDp6cu5UTP+tAOvW/NASU8nChwXfw2fD2/L13GJsOn6fbhGQ+2qQD6VQ2DQGlbMDbS3iiQ21WjImkec1yvLoolQEzN3H4bIbVpSmLFVsIiEhfEXlPRD4VkW6OtUAR+cCxPqi4alHKrmpWKM1HQ1vzr4easvfUJXpMWsuMJB1IZ2f5/Y3h2SJyRkRSc63HiMg+ETkoImNv9xzGmEXGmKeB4UB/x/KDwELHep9C1K+UKiAR4ZFWNVkVG0XnBpV58+u99H1nPXtO6EA6O8rvmcBcICbngoh4A9OAHkAYMFBEwkSkiYgk5LpVybHrK479AGoAxxx/1ilYShWjKmUCeHdwBNMHteTUxav0mbqOt5fv4+oN/b+inUh+LxkTkRAgwRjT2LHdDhhnjOnu2H4JwBgzPo/9BXgTWGmMWeVYGwxcMMYkiEi8MWbALfYbBgwDCA4ODo+Pjy9QgxkZGQQFBRVon5JOe7YPZ/Wdcd0wf+911p/IpGqg8GRjf0LLezuhQufTY11w0dHRKcaYiFveaYzJ1w0IAVJzbPcDZuXYHgxMvc3+o4AUYAYw3LEWCMwBpgODfq+G8PBwU1CrV68u8D4lnfZsH87ue82+M6b9+G9MyNgE8/fFqSbj6g2nPr8z6LEuOGCbyeN91adQsVIIxpjJwORca5eBJ4qrBqXU7UXVr8zyMZH8e9le5m44wqq004x/sAmdQitbXZpykaJcHXQcqJlju4ZjTSlVggX5+/Da/Y35bHg7/Ly9GPz+Fv782Q4uXtGBdJ6oKCGwFQgVkdoi4gcMAJY4pyyllNVahVQgcXQnnu1cly++P06XCUksS9WBdJ4mv5eIzgc2Ag1EJF1EhhpjMoGRwHIgDVhgjNntulKVUsUtwNebF2Masvi5DlQO8mf4xyk8Oy+FM79ctbo05ST5+k7AGDMwj/VEINGpFSml3E7j6mVZPLIDM5MPM+mbA6w/+BOv9grjoZbVyb7wT5VUOjZCKZUvvt5ePBddj8RRnQitEsSfPtvBkDlbSb9wxerSVBFoCCilCqRelSAWPNOO1/o0YtuR7IF0c9f/oAPpSigNAaVUgXl5CUPah7BiTCQRIRUYt3QPj7y7kYNndCBdSaMhoJQqtBrlS/PBE614++FmHDiTwX2T1jJt9UFu6EC6EkNDQClVJCJCv/AarIyNpEtYFf69fB/3T11P6vGLVpem8kFDQCnlFFXuCOCdQeHM+ENLzmZc4/5p63lr2V4dSOfmNASUUk4V07gqq8ZE8WCL6kxfc4j7Jq1l65HzVpel8qAhoJRyurKlffn3w834aGhrrmVm8fCMjfxtcSoZ1zKtLk3loiGglHKZTqGVWTEmksfbh/DRpqN0n5BM0v6zVpelctAQUEq5VKC/D+P6NGLh8HYE+HoxZPYWYhds5+cr160uTaEhoJQqJuF3VeCrUZ0YGV2PJdtP0CUuicRdJ3/7vRFlEQ0BpVSxCfD15k/dG7B4ZAfuLBvAs/O+Y/jHKZy5pAPprKIhoJQqdo2qlWXRsx0Y26Mha/adpUtcEgu2HdOzAgtoCCilLOHj7cXwqLp8PboTDauW4cWFOxn8/haOndeBdMVJQ0ApZak6lYOIf7ot/9O3MduP/Uy3CcnMXvcDN3UgXbHQEFBKWc7LSxjc9i5WjImkTZ0K/CNhD/1mbODA6V+sLs3jFVsIiEhfEXlPRD4VkW451gNFZJuI9CquWpRS7qlauVLMebwVE/s358i5y/ScvI4p3xzQgXQulN+fl5wtImdEJDXXeoyI7BORgyIy9nbPYYxZZIx5GhgO9M9x11+ABQUtXCnlmUSEvi2qszI2iu6N7+Q/K/fTe8o6dqb/bHVpHim/ZwJzgZicCyLiDUwDegBhwEARCRORJiKSkOtWJceurzj2Q0S6AnuAM0XsQynlYSoF+TNlYAveeyyCC1eu03faesYnpnH9pn5X4EyS30uyRCQESDDGNHZstwPGGWO6O7ZfAjDGjM9jfwHeBFYaY1Y51v4JBJIdIr8CDxhjsnLtNwwYBhAcHBweHx9foAYzMjIICgoq0D4lnfZsH3bp+/INw4J910lKz6RygGFo01I0rOBtdVnFqijHOjo6OsUYE3Gr+/L1Q/N5qA4cy7GdDrS5zeOfB7oAZUWknjFmhjHmrwAi8jhwLncAABhjZgIzASIiIkznzp0LVOSaNWso6D4lnfZsH3bqu2dX2HDwHKPnbeHNLVcZ1KYWY3s05I4AX6tLKxauOtZFCYECMcZMBibncd/c4qpDKVVyta9Xidc7lGLL1WBmr/+Bb/ee4Y0HmhDdsMrv76xuqShXBx0HaubYruFYU0opl/H3EV7tFcbnI9oT5O/DE3O38kL895y/rAPpCqMoIbAVCBWR2iLiBwwAljinLKWUur2WtcqTMKojo+8NJWHnSbrGJbF0xwkdPVFA+b1EdD6wEWggIukiMtQYkwmMBJYDacACY8xu15WqlFL/l7+PN2O61idhVEdqlC/F8/O/5+kPUzitA+nyLV/fCRhjBuaxnggkOrUipZQqoIZ3luHzEe2Zs/4Ib6/YR5e4JP563930b1WT7AsTVV50bIRSyiP4eHvxdGQdlr8QSaNqZRj7xS4GzdrM0Z8uW12aW9MQUEp5lJBKgXzyVFveeKAJu9Iv0n1iMrPWHtaBdHnQEFBKeRwvL+HRNrVYERtJh7qVeP2rNB6cvoF9p3QgXW4aAkopj1W1bClmDYlg0oDmHDt/hV5T1jJx1X6uZ+pAut9oCCilPJqIcH/z6qwcE0mPxlWZuOoAvaesY/uxn60uzS1oCCilbKFikD+TB7Zg1mMRXPz1Bg++s57XE/bw6/WbVpdmKQ0BpZStdAkLZkVsJP1b1WLWuh/oPjGZDYfOWV2WZTQElFK2UybAl/EPNmH+020RgUff28xLX+zk0tUbVpdW7DQElFK21a5uRZaNjmRYZB0+3XqMrnFJrNpz2uqyipWGgFLK1kr5efPyfXfz5bMdKFfKj6c+3Mao+d/zU8Y1q0srFhoCSikFNKtZjqXPd+SFLqF8nXqSLnFJLN5+3OMH0mkIKKWUg5+PFy90qc9XozpxV8VARsdv56kPtnHy4q9Wl+YyGgJKKZVL/eA7+HxEe17tFcaGQz/RNS6ZeZuPkuWBoyc0BJRS6ha8vYShHWuz/IVImtUsy1+/TOXRWZs4cs6zBtJpCCil1G3Uqliaj4e24c0Hm7D7xCW6T0xmZvIhMm96xugJDQGllPodIsKA1rVYFRtFZP3KvJG4lwenbyDt5CWrSysyDQGllMqn4DIBzBwcztRHW3Di51/pPWUdcSv2cS2z5I6eKLYQEJG+IvKeiHwqIt0ca7VEZJGIzBaRscVVi1JKFZaI0KtpNVaOiaJ3s2pM/vYgvSav47sfL1hdWqHk9zeGZ4vIGRFJzbUeIyL7ROTg772JG2MWGWOeBoYD/R3LTYCFxpgngRaFqF8ppSxRPtCPCf2bM+fxVmRcy+Sh6Rv4x9I9XLmeaXVpBZLfM4G5QEzOBRHxBqYBPYAwYKCIhIlIExFJyHWrkmPXVxz7AWwChorIt8CyojSilFJWiG5YhRVjIhnUphaz12cPpFt/sOQMpJP8/tdwIhICJBhjGju22wHjjDHdHdsvARhjxuexvwBvAiuNMasca38CthhjkkVkoTGm3y32GwYMAwgODg6Pj48vUIMZGRkEBQUVaJ+STnu2Dzv27c497zt/k9mp1zh9xRBZw4f+DfwI9HXOD90Xpe/o6OgUY0zEre7zKUJN1YFjObbTgTa3efzzQBegrIjUM8bMIPtv/+NE5FHgyK12MsbMBGYCREREmM6dOxeoyDVr1lDQfUo67dk+7Ni3O/fcGRjS+yYTVx3gvbWH2XvxJq/3bUy3RncW+bld1XdRQqBAjDGTgcm51lKB//rbv1JKlVQBvt6M7dGQnk2q8uLnOxn2UQq9mlZlXJ9GVAryt7q8/1KUq4OOAzVzbNdwrCmllO01qVGWJSM78Kdu9Vmx+zRd4pL48vt0txtIV5QQ2AqEikhtEfEDBgBLnFOWUkqVfL7eXoy8J5TE0R2pUymQMZ/u4Im5Wzn+s/sMpMvvJaLzgY1AAxFJF5GhxphMYCSwHEgDFhhjdruuVKWUKpnqVbmDz4a3Z1zvMLb8cJ5ucUl8tPGIWwyky9d3AsaYgXmsJwKJTq1IKaU8kLeX8HiH2tx7dzAvf7mLVxfvZumOk7z5UBPqVLbuaicdG6GUUsWoZoXSfPhka/7dryl7T10iZtJapq+xbiCdhoBSShUzEeHhiJqsio3ingZVeGvZXvq+s549J4p/IJ2GgFJKWaRKmQBmDA5n+qCWnLp4jT5T1/H28n1cvVF8A+k0BJRSymI9mlRlVWwk9zevztTVB+k5eS0pR88Xy2trCCillBsoV9qP/zzSjA+ebM3VG1n0m7GRcUt2c/maawfSaQgopZQbiapfmeVjInms7V18sPEI3SYkk7z/rMteT0NAKaXcTJC/D6/d35jPnmmHv68Xj83ewoe7r7nktTQElFLKTUWEVCBxVCeei65LldKuebvWEFBKKTcW4OvNn7s3JKa2r0ueX0NAKaVsTENAKaVsTENAKaVsTENAKaVsTENAKaVsTENAKaVsTENAKaVsTENAKaVsTNztR49vR0TOAkcLuFsl4JwLynFn2rN92LFvO/YMRev7LmNM5VvdUaJCoDBEZJsxJsLqOoqT9mwfduzbjj2D6/rWj4OUUsrGNASUUsrG7BACM60uwALas33YsW879gwu6tvjvxNQSimVNzucCSillMqDhoBSStmYx4aAiMSIyD4ROSgiY62uxxVEpKaIrBaRPSKyW0RGO9YriMhKETng+Gd5q2t1BRHxFpHvRSTBsV1bRDY7jvmnIuJndY3OJCLlRGShiOwVkTQRaWeHYy0iYxz/fqeKyHwRCfDEYy0is0XkjIik5li75fGVbJMd/e8UkZaFfV2PDAER8QamAT2AMGCgiIRZW5VLZAJ/NMaEAW2B5xx9jgW+McaEAt84tj3RaCAtx/ZbwARjTD3gAjDUkqpcZxKwzBjTEGhGdu8efaxFpDowCogwxjQGvIEBeOaxngvE5FrL6/j2AEIdt2HA9MK+qEeGANAaOGiMOWyMuQ7EA/dbXJPTGWNOGmO+c/z5F7LfFKqT3esHjod9APS1pEAXEpEaQE9glmNbgHuAhY6HeFTfIlIWiATeBzDGXDfG/IwNjjXgA5QSER+gNHASDzzWxphk4Hyu5byO7/3AhybbJqCciFQtzOt6aghUB47l2E53rHksEQkBWgCbgWBjzEnHXaeAYKvqcqGJwItAlmO7IvCzMSbTse1px7w2cBaY4/gIbJaIBOLhx9oYcxx4G/iR7Df/i0AKnn2sc8rr+DrtPc5TQ8BWRCQI+Bx4wRhzKed9JvsaYI+6DlhEegFnjDEpVtdSjHyAlsB0Y0wL4DK5Pvrx0GNdnuy/9dYGqgGB/PdHJrbgquPrqSFwHKiZY7uGY83jiIgv2QEwzxjzhWP59G+nho5/nrGqPhfpAPQRkSNkf9R3D9mfl5dzfGQAnnfM04F0Y8xmx/ZCskPB0491F+AHY8xZY8wN4Auyj78nH+uc8jq+TnuP89QQ2AqEOq4g8CP7i6QlFtfkdI7Pwd8H0owxcTnuWgIMcfx5CLC4uGtzJWPMS8aYGsaYELKP7bfGmEHAaqCf42Ee1bcx5hRwTEQaOJbuBfbg4cea7I+B2opIace/77/17bHHOpe8ju8S4DHHVUJtgYs5PjYqGGOMR96A+4D9wCHgr1bX46IeO5J9ergT2O643Uf25+PfAAeAVUAFq2t14f8GnYEEx5/rAFuAg8BngL/V9Tm51+bANsfxXgSUt8OxBl4D9gKpwEeAvycea2A+2d973CD7zG9oXscXELKvgDwE7CL76qlCva6OjVBKKRvz1I+DlFJK5YOGgFJK2ZiGgFJK2ZiGgFJK2ZiGgFJK2ZiGgFJK2ZiGgFJK2dj/A7d8k67s3QkiAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# 1/2**i → exp2(-i)\n", "plt.semilogy(i,np.exp2(-i))\n", "plt.grid()" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "cell_id": "00018-e8d4f228-6cce-4727-9bf2-f0768560ba70", "colab": {}, "colab_type": "code", "deepnote_cell_type": "code", "id": "_Eye-yOJw7V7", "outputId": "f83011e3-ca2f-438b-c6e7-280f27334322" }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZkAAAF4CAYAAACGvYTpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAABCzElEQVR4nO3dd3hUZfr/8fedRi+CgAgiiMqKWGGtK0kAJSgCAtKkiqKuuuuuu65d1rqW1bWwuhYICBIB6SI9CXaxrYKKIqKCBRFRaQkk9++PGb6/mE0kQ2ZyJsnndV1zwXnmzDmfJ4crN6fM85i7IyIiEgsJQQcQEZGqS0VGRERiRkVGRERiRkVGRERiRkVGRERiRkVGRERiJinoAPHkwAMP9NatW5d5/e3bt1OnTp3YBYpj1bXv6nf1on6XzVtvvbXZ3ZuU9J6KTBGtW7fmzTffLPP6OTk5pKWlxS5QHKuufVe/qxf1u2zM7PPS3tPlMhERiRkVGRERiRkVGRERiRkVGRERiRkVGRERiRkVGRERiZkq/QizmdUB/g3kAznuPiXgSCIi1UqlO5Mxs/FmtsnMVhVrzzCzNWa21syuDTf3BWa4+8VArwoPKyJSzVW6IgNkAhlFG8wsERgH9ADaA4PNrD3QEvgyvFpBBWYUERHAKuPMmGbWGpjv7h3Cy6cCY929e3j5uvCqG4Af3H2+mWW5+6AStjUGGAPQrFmzjllZWWXOsW3bNurWrVuuvlRW1bXv6nf1on6XTXp6+lvu3qmk96rKPZkW/P8zFggVl5OBh4BHzOwcYF5JH3T3x4HHATp16uSRDKVQXYecgOrbd/W7elG/y68yXi4rM3ff7u6j3P2yWN7037JzC72m9mLtlrWx2oWISKVUVYrMRuCQIsstw20V4sblNzLv43l0ntCZjzZ/VFG7FRGJe1WlyKwEjjCzNmaWAgwC5lbUzu858x7SWqfx9bavSc1MZdWmVfv+kIhINVDpioyZTQVeBdqZ2QYzG+3ue4ArgEXAh8A0d19dUZnqptTl+SHP0+2wbmzavom0zDTe/ebditq9iEjcqnRFxt0Hu3tzd09295bu/lS4fYG7H+nubd39jorOVTu5NvMGz6PH4T34fuf3dJnYhTe/KvvcNCIiVVGlKzLxrGZSTWYNnEXvdr35YdcPdJ3UlVe/fDXoWCIigVGRibIaSTWYfv50+rfvz095P3HW5LN48fMXg44lIhIIFZkYSE5MZmq/qQw5Zgjb8reRMSWD5Z8tDzqWiEiFU5GJkaSEJCb1mcTI40eyY/cOznnmHBatXRR0LBGRCqUiE0OJCYk81espxpw4hl17dtErqxfzP54fdCwRkQqjIhNjCZbAYz0f44rfXkF+QT59n+3LrA9nBR1LRKRCqMhUADPjoR4P8edT/szuwt2cP/18pq2eFnQsEZGYU5GpIGbGfWfdx3W/u44CL2Dwc4OZ/N7koGOJiMSUikwFMjPu6HIHY1PHUuiFDJ81nAnvTAg6lohIzKjIVDAz45a0W7izy504zoVzL+Q/b/4n6FgiIjGhIhOQ6864jvvOvA+AS5+/lIdffzjgRCIi0aciE6CrT7uahzIeAuAPC//AP1/5Z8CJRESiS0UmYFeefCX/6Rm6XPaXJX/hzhfvDDiRiEj0qMjEgTEdxzC+13gM44blN3BL9i24e9CxRETKTUUmTow6YRRPn/c0CZbArStu5fpl16vQiEilpyITRy449gKm9ptKoiXyj5f/wdWLr1ahEZFKTUUmzgw4egAzBswgOSGZB157gCtfuJJCLww6lojIflGRiUN9ftOHWQNnUSOxBuNWjuPS+Zeq0IhIpaQiE6fOOfIc5g6eS82kmjzx9hNcOOdCCgoLgo4lIhIRFZk4dlbbs1gwZAG1k2sz8b8TGT57OHsK9wQdS0SkzFRk4lx6m3QWXrCQuil1eeb9Zxj83GB2F+wOOpaISJmoyFQCZxx6BkuGLaFBjQbM+GAG/af3J29PXtCxRET2SUWmkjil5SksHb6UA2oewNw1cznv2fPYtWdX0LFERH6Vikwl0ungTiwfsZwDax/IC2tf4Nyp57Jj946gY4mIlEpFppI5/qDjyR6RTbM6zVi6binnPHMO2/K3BR1LRKREKjKVUIemHcgZmUPzus3JWZ9DxuQMfsr7KehYIiL/Q0WmkvrNgb9hxagVHFL/EF7+8mXOfPpMtu7aGnQsEZFfUJGpxA5vdDi5I3Np3bA1b2x8g66TuvL9ju+DjiUi8n9UZCq5Nge0IXdkLm0PaMvbX79Nl0ld+G77d0HHEhEBVGSqhFYNWrFi1AraNW7He9++R9rENL7Z9k3QsUREqn6RMbM6ZvammfUMOkssHVzvYHJH5nJ0k6P54LsPSM1MZcNPG4KOJSLVXNwWGTMbb2abzGxVsfYMM1tjZmvN7NoybOpvwLTYpIwvzeo2I3tENsc1O46Pv/+Y1MxUPt/6edCxRKQai9siA2QCGUUbzCwRGAf0ANoDg82svZkdY2bzi72amtmZwAfApooOH5QmdZqwfMRyOjbvyLof1pGamcq6H9YFHUtEqqm4LTLuvgLYUqz5JGCtu69z93wgC+jt7u+7e89ir01AGnAKMAS42Mzitr/R1KhWI5YOX8opLU/h8x8/p/OEznz8/cdBxxKRasjieXpfM2sNzHf3DuHl/kCGu18UXh4GnOzuV+xjOyOBze4+v4T3xgBjAJo1a9YxKyurzPm2bdtG3bp1y7x+RduxZwfXvn8t7//0Po1SGnH/sfdzaJ1Do7LteO97rKjf1Yv6XTbp6elvuXunEt9097h9Aa2BVUWW+wNPFlkeBjwSrf117NjRI5GdnR3R+kHYlrfN0zPTnbF4k3ua+HvfvBeV7VaGvseC+l29qN9lA7zppfxerWyXjzYChxRZbhluk1LUSanD/CHzOavtWXy34zvSJ6bzztfvBB1LRKqJylZkVgJHmFkbM0sBBgFzA84U92on12bOoDmcc8Q5fL/ze7pM6sLKjSuDjiUi1UDcFhkzmwq8CrQzsw1mNtrd9wBXAIuAD4Fp7r46yJyVRc2kmswcOJM+v+nD1l1b6fZ0N1758pWgY4lIFRe3RcbdB7t7c3dPdveW7v5UuH2Bux/p7m3d/Y6gc1YmKYkpTOs/jQFHD+CnvJ/oPrk7Kz5fEXQsEanC4rbISGwkJyYzpe8Uhh47lG352+gxpQfL1i0LOpaIVFEqMtVQUkISmb0zGXX8KHbs3kHPqT1ZuHZh0LFEpApSkammEhMSebLXk1za8VJ27dlF76zezFszL+hYIlLFqMhUYwmWwL/P+Td/OOkP5Bfk03daX2Z+ODPoWCJShajIVHNmxr8y/sVfTv0Lewr3MGD6AJ5d9WzQsUSkilCREcyMe868hxvOuIECL2DIzCE8/d+ng44lIlWAiowAoUJze5fbuTXtVgq9kBGzRzD+nfFBxxKRSk5FRn7hptSbuKvrXTjO6LmjeXTlo0FHEpFKTEVG/se1v7uW+8+6H4DfL/g9D772YMCJRKSyUpGREv3p1D/xSI9HALhq0VXc+/K9AScSkcpIRUZKdflJl/N4z8cxjGuWXsMdKzSKj4hERkVGftXFHS9mfO/xGMaN2Tdyc/bNe+fyERHZp6SgA0j8G3n8SFISUxg+azi3rbiNvD15ZCRlBB1LRCoBnclImQw5ZghZ/bNISkjinlfuYdyn43RGIyL7pCIjZda/fX9mnD+D5IRkntv4HJcvuJxCLww6lojEMRUZiUjv3/Rm9qDZJFsyj775KJfMu0SFRkRKpSIjETv7iLO5s8Od1EqqxZPvPMmoOaMoKCwIOpaIxCEVGdkvnRp1YsEFC6iTXIdJ/53E0FlD2V2wO+hYIhJnVGRkv6W1TmPR0EXUS6lH1qosBj03iPyC/KBjiUgcUZGRcjm91eksGbaEBjUaMPPDmfSf1p+8PXlBxxKROKEiI+V2csuTWT5iOY1qNWLex/Po82wfdu7eGXQsEYkDKjISFSc2P5Hlw5dzYO0DWbh2IedOPZcdu3cEHUtEAqYiI1Fz3EHHkTMih2Z1mrHss2X0mNKDn/N+DjqWiARIRUai6uimR5M7MpeD6x3Mis9XkDElgx93/Rh0LBEJyD7HLjOzVuXZgbt/UZ7PS+XT7sB25I7MpcvELrzy5SucNfksFl6wkANqHRB0NBGpYGUZIHM9sD+DVFn4c4n78Vmp5A5vdDgrRq2gy8QuvLHxDbpO6sqSYUtoXLtx0NFEpAKVpci0iXkKqZJaN2xN7shcuk7qyjvfvEP6xHSWDl9K0zpNg44mIhVkn0XG3T+viCBSNR3S4JDQpbNJXXh/0/ukZaaxbPgymtdrHnQ0EakAEd34N7NkM7vKzBaa2RtmNtvMrjCz+rEKKJVf83rNyRmRQ4emHfhw84ekZqay4acNQccSkQoQ6dNlDwP3E7rPshIoAG4D1ptZvyhnkyqkWd1mZI/I5viDjueTLZ+QmpnK51t1kixS1UVaZM4Hbnb3M939cnfvBxwC/BOYYmY9o56wHMysVfhsa7yZXRt0nuruwNoHsmz4Mjod3Il1P6yjc2ZnPt3yadCxRCSGIi0yDrz0iwb3be5+B/AAcHu0goULwyYzW1WsPcPM1pjZ2jIUjmOAGe5+IXBCtLLJ/mtUqxFLhy3l1Jan8sWPX5CamcqazWuCjiUiMRJpkVkInFfKe4uAduWL8wuZwC8mkjezRGAc0ANoDww2s/ZmdoyZzS/2agq8Bow2s+Xh7BIHGtRswKKhi+h8aGc2/ryR1MxUPvjug6BjiUgMRFpk1gGjzOwOMyv+zbozgNXRiQXuvgLYUqz5JGCtu69z93wgC+jt7u+7e89ir03AKOAWd+8CnBOtbFJ+9WrUY8GQBXRp04Vvt39LWmYa7337XtCxRCTKzL3s37M0s++BvcVlO7AC+ApoG36d4+6rSvl45OHMWgPz3b1DeLk/kOHuF4WXhwEnu/sVpXy+AzAW2Axsc/e/lLDOGGAMQLNmzTpmZWWVOd+2bduoW7duJF2qMqLV97yCPG5afRMrf1hJ/aT63HvsvRxZ78goJIyN6nrM1e/qJdJ+p6env+XunUp8090jehG60d8TuAGYBqwB9gCFwA+ECs/DwEWRbruEfbUGVhVZ7g88WWR5GPBIefez99WxY0ePRHZ2dkTrVyXR7PvO3Tu95zM9nbF4w3809Nc3vB61bUdbdT3m6nf1Emm/gTe9lN+rEQ+Q6e5fuvt8d7/D3Qe4ezugPnAq8DfgfUI32e+LdNtlsJFQkdurZbhNKrGaSTV5bsBz9D2qL1t3baXbpG68/MXLQccSkSiIyijM7r7D3V9398c99Gjz79y9YTS2XcxK4Agza2NmKcAgYG4M9iMVLCUxhax+WQw8eiA/5/9M98ndyV2fG3QsESmnuB3q38ymAq8C7cxsg5mNdvc9wBWEnmT7EJjm7lF72ECClZyYzOS+kxl27DC2795Ojyk9WLpuadCxRKQcyjJAZiDcfXAp7QuABRUcRypIUkISE3pPIDkhmfHvjqfnMz2ZNXAWPY7oEXQ0EdkPcXsmI9VXYkIiT/R6gss6XUZeQR59nu3D3DW6KipSGanISFxKsATGnT2OP578R/IL8uk3rR8zPpgRdCwRidB+Fxkza2dmBdEMI1KUmfFA9wf462l/ZU/hHgbNGMTU96cGHUtEIlDeMxmLSgqRUpgZd3e7mxvPuJECL2DorKFMfHdi0LFEpIzKW2T2Z1pmkYiYGbd1uY1b026l0AsZNWcUT779ZNCxRKQMdE9GKo2bUm/i7m534zgXz7uYcW+MCzqSiOyDioxUKtecfg0PdH8AgCteuIIHXn0g4EQi8mtUZKTSueqUq/j32f8G4M+L/8zdL90dcCIRKY2KjFRKl/32Mp4890kM49pl13Jb7m1BRxKREqjISKU1+sTRZPbJJMESuDnnZm5cfuPe0blFJE7oEWap1IYfN5wpfaeQaInc8eId/G3p31RoROJIeYrM18DF0Qoisr8GdRjEs/2fJSkhiXtfuZc/LfqTCo1InNjvIuPuP7n7U9EMI7K/+rXvx8wBM0lJTOHB1x/k98//nkIvDDqWSLWnezJSZZzb7lzmDJpDzaSaPPbWY4yZN4aCQo18JBIkFRmpUjIOz2D+4PnUSqrFU+88xag5o9hTuCfoWCLVloqMVDldD+vKCxe8QJ3kOjz93tMMnTmU3QW7g44lUi2pyEiVlNo6lcXDFlMvpR7Prn6WgTMGkl+QH3QskWqn3EXGzK43s6/M7D0ze9rM/mxmXaIRTqQ8TjvkNJYOX0rDmg2Z9dEs+k3rR96evKBjiVQr0TiTGQMcD3QHngFqApdGYbsi5XZSi5NYNnwZjWo1Yv7H8+md1Zudu3cGHUuk2ohGkfnE3Te5+9fu/oK73+nuA6KwXZGoOLH5ieSMyKFJ7SYs+nQRPaf2ZHv+9qBjiVQL0Sgy75rZ/WZWKwrbEomJY5odQ+7IXA6qexDLP1tOjyk9+Dnv56BjiVR50SgyDYHfARvM7FUz+7eZjYnCdkWi6qgmR5E7MpcW9Vrw4hcv0n1yd37c9WPQsUSqtHIXGXe/2N1PApoAFwIrgDbl3a5ILBzZ+EhyR+bSqkErXt3wKt2e7sYPO38IOpZIlRVRkTGzZDO7yswWmtkbZjbbzK4ws3ruXujuH7p7lrtfF6vAIuXVtlFbVoxcwWEHHMabX71Jl0ld2Lxjc9CxRKqkSM9kHgbuBxKBlUABcBvwuZn1i3I2kZg5tOGh5I7M5cjGR/LuN++SPjGdb7d9G3QskSon0iJzPnCzu5/p7pe7ez/gEOCfwBQz6xn1hCIx0rJ+S3JG5NC+SXtWbVpF2sQ0vvr5q6BjiVQpkRYZB176RYP7Nne/A3gAuD1awUQqQvN6zckekc0xTY/ho80fkZqZypc/fhl0LJEqI9IisxA4r5T3FgHtyhdHpOI1rdOU7BHZnHDQCazdspbUzFTWb10fdCyRKiHSIrMOGGVmd5jZAcXeOwNYHZ1YIhWrce3GLBu+jJNanMRnWz+j84TOfLrl06BjiVR6kRaZy4G6wHXAF2b2vJk9YWbLgYuAkVHOJ1JhDqh1AEuGLeG0Q07jy5++pHNmZ9ZsXhN0LJFKLaIi4+6NgUOBXsA/gO1A5/DrEOBFM1thZg+b2UXRDisSa/Vr1GfhBQvpfGhnvvr5K1IzU1m9SSfoIvsr4i9juvuX7j7f3e9w9wHu3g6oD5wK/A14HzgBuC+6UffNzA4zs6fMbEaRtj7hs61nzeysis4klU+9GvVYMGQBXdt05dvt35I2MY3/fvPfoGOJVEpRmU/G3Xe4++vu/nj40ebfuXvDSLZhZuPNbJOZrSrWnmFma8xsrZldu48c69x9dLG22e5+MaGRoQdGkkmqrzopdZg3eB4Zh2ewecdm0iem89ZXbwUdS6TS2WeRMbNW5XlFkCUTyCi270RgHNADaA8MNrP2ZnaMmc0v9mq6j+3fGN6WSJnUSq7F7IGzOffIc/lh1w90ndSVD376IOhYIpWKufuvr2BWSOj7MRFvG3B3TyzzB8xaA/PdvUN4+VRgrLt3Dy9fR2ijd+1jOzPcvX/470bo/tESd19awrpjCM2JQ7NmzTpmZWWVNS7btm2jbt26ZV6/KqlOfd9duJvbP7ydFZtXUCuhFncfezfHNDgm6FgVqjod76LU77JJT09/y907lfReUhk+H+Rgly2Aot+M2wCcXNrKZtYYuAM4wcyuCxejK4FuQAMzO9zdHyv6GXd/HHgcoFOnTp6WllbmcDk5OUSyflVS3fqenpbOsFnDyFqVxbWrr+X5Ic+T1jot6FgVprod773U7/LbZ5Fx98+jsqcK4O7fU2xWTnd/CHgomERSVSQlJDH5vMls+W4Li79dTI8pPZgzaA5ntdWzJCK/Jio3/mNoI6FHo/dqGW4TqXCJCYn8rd3fuOiEi9i1Zxe9pvZiwScLgo4lEtfivcisBI4wszZmlgIMAuYGnEmqsQRL4D/n/offd/o9eQV59Mnqw5yP5gQdSyRuxU2RMbOpwKtAOzPbYGaj3X0PcAWhcdE+BKa5u74ZJ4FKsAQeOfsRrjr5KnYX7qb/9P5MXz096FgicaksN/4rhLsPLqV9AaBrEhJXzIz7u99PjaQa3P3y3Qx6bhC7C3cz5JghQUcTiStlPpMJz4p5upkdHMtAIpWFmXFX17u4qfNNFHohQ2cOJfPdzKBjicSVSC6XFQDLgd/EKItIpWNm3Jp+K7el34bjjJozisffejzoWCJxo8xFxt0LgU+Ag2IXR6RyurHzjdzT7R4ALpl/CY+88UjAiUTiQ6Q3/m8Abjaz6vV1Z5Ey+Ovpf+Vf3f8FwJUvXMn9r94fbCCROBDpjf8bgcbAu2a2EfiWYkPOuPtJUcomUun88ZQ/UiOpBpc9fxlXL76avD15XHfGdUHHEglMpEVmVfglIqW4tNOlpCSmcNHci7h++fXkF+Rzc+rNhIbRE6leIioy7j4qVkFEqpILT7iQlMQURswewdjcseQX5HN7l9tVaKTa2a/vyYQfYz4VaAR8D7zm7l9FM5hIZTf02KEkJyRzwcwLuPOlO8kryOPeM+9VoZFqJaIiE57f5WHgYqDoEP4FZvY4cGX4KTQRAQZ2GEhKYgoDZwzkn6/+k/yCfB7MeFCFRqqNSJ8u+ztwIXA90BqoFf7z+nD72OhFE6kazjvqPGYOnElKYgoPv/Ewl86/lEL9X0yqiUiLzHDgRne/192/cPe88J/3AjcBI6OeUKQK6HlkT+YOmkvNpJo8/vbjjJ47moLCgqBjicRcpEWmKfBeKe+9F35fRErQ/fDuPD/keWon1ybz3UxGzB7BnsI9QccSialIi8zHhIbbL8kgYE354ohUbV3adGHhBQupm1KXKe9PYchzQ9hdsDvoWCIxE+nTZbcDWWbWCphB6MuYTYHzgXRKL0AiEnbGoWeweOhiMqZkMP2D6ewu3E1WvyxqJNUIOppI1EV0JuPu04AMoA7wIPAcoamNawMZ7q5JNUTK4NRDTmXZ8GU0rNmQ2R/Npu+0vuzasyvoWCJRF/GkZe6+2N1PJfRk2UFALXc/zd2XRD2dSBXW6eBOZI/IpnGtxiz4ZAG9pvZix+4dQccSiar9nk/G3QvdfZO+FyOy/44/6HhyRubQtE5TlqxbQs9nerItf1vQsUSiRvPJiASsQ9MO5IzIoXnd5mSvz6bHlB78lPdT0LFEokLzyYjEgaOaHEXuyFxa1m/JS1+8xFlPn8XWXVuDjiVSbppPRiROHNH4CFaMXMGhDQ7l9Y2v021SN77f8X3QsUTKJdIiU3Q+mS/MbKWZvVH0FYOMItVGmwPasGLUCtoe0Ja3vn6LrpO68t3274KOJbLfNJ+MSJxp1aAVuSNz6TKpC//99r+kT0xn6fClHFRXV6ql8ilzkTGzZOBJYL27b4xdJBFpUb8FuSNz6TqpK6u/W01aZhrLRyzn4HoHBx1NJCL783RZuxhlEZEiDqp7EDkjcji22bGs+X4NnSd05osfvwg6lkhE9HSZSBxrUqcJy4cv58TmJ/LpD5+SmpnKZz98FnQskTLT02Uica5x7cYsG76Mk1qcxPqt60nNTGXtlrVBxxIpEz1dJlIJNKzZkCXDlnD6Iafz5U9fkpqZykebPwo6lsg+6ekykUqifo36LBy6kHOnnkvO+hzSMtNYOnwpHZp2CDqaSKkiKjLuPipWQURk3+qm1OX5Ic/TJ6sPS9YtCT3ePGwpxx10XNDRREq0z8tlZjbEzBoVa2tlZknF2g42s+ujHVBEfql2cm3mDp7L2UeczeYdm0mfmM6bX70ZdCyREpXlnszTwOF7F8wsEfgMOLbYeocAt0UvWmTM7DAze8rMZhRpSzCzO8zsYTMbEVQ2kWirmVSTmQNm0rtdb37Y9QNdJ3XltQ2vBR1L5H+UpchYGdv2m5mNN7NNZraqWHuGma0xs7Vmdu2vbcPd17n76GLNvYGWwG5gQzQziwStRlINpp8/nf7t+/NT3k+c+fSZvPTFS0HHEvmFiCcti5FMQjNu/p/wGdM4oAfQHhhsZu3N7Bgzm1/s1bSU7bYDXnH3PwOXxTC/SCCSE5OZ2m8qQ44Zwrb8bXSf3J3sz7KDjiXyfyJ9uiwm3H2FmbUu1nwSsNbd1wGYWRbQ293vAnqWcdMbgPzw3wtKWsHMxgBjAJo1a0ZOTk6Zc2/bti2i9auS6tr3eO33hQdcyJZmW1j47UIyJmdw+9G389tGv43a9uO137GmfkeBu//qCygEfltkOTHcdkKx9U4GCva1vV/ZT2tgVZHl/sCTRZaHAY/8yucbA48BnwLXhdtqA08BDwOX7ytDx44dPRLZ2dkRrV+VVNe+x3O/CwoL/OK5Fztj8ZTbUnzemnlR23Y89zuW1O+yAd70Un6vlvVMZpGZ7SnWtqxYW6BnRe7+PXBpsbYdQPH7NCJVUoIl8FjPx0hJTGHcynH0fbYvz/Z/lvOOOi/oaFKNlaUw/D3mKUq2kdATa3u1DLeJSCkSLIGHezxMjcQa3P/a/Zw//Xye6fcMA44eEHQ0qab2WWTcPagisxI4wszaECoug4AhAWURqTTMjPvOuo8aSTW466W7GPzcYPIL8hl67NCgo0k1FBdPl5nZVOBVoJ2ZbTCz0e6+B7gCWAR8CExz99VB5hSpLMyMO7rcwdjUsRR6IcNnDWf8O+ODjiXVULw8XTa4lPYFwIIKjiNSJZgZt6TdQkpiCtcvv57Rc0ezu2A3l3S6JOhoUo3ExZmMiMTOdWdcx31n3gfApc9fysOvPxxwIqlOVGREqoGrT7uahzIeAuAPC//Afa/cF3AiqS5UZESqiStPvpLHznkMgL8u+St3rLgj4ERSHexXkQkP7zLMzK43s4PCbYebWb3oxhORaLqk0yWM7zUew7gx+0Zuyb5l75eZRWIiohv/ZlYXGA/0A/aEP78Q+Aa4E/gC+EuUM4pIFI06YRTJicmMmD2CW1fcSn5BPnd2vROzqI57KwJEfiZzP3Aa0A2oxy9HY15AsUEuRSQ+DT12KFP7TSXREvnHy//g6sVX64xGYiLSItMX+Ju7Z/O/A05+DhwalVQiEnMDjh7A9POnk5yQzAOvPcCVL1xJoRcGHUuqmEiLTC3g+1Leq0cpIx2LSHw676jzmDVwFjUSazBu5TgumXeJCo1EVaRFZiUwvJT3+gOvlC+OiFS0c448h7mD51IzqSZPvvMkF865kIJC/X9RoiPSInMT0NfMlgIXAQ6cbWZPA+cDt0Q5n4hUgLPansWCIQuonVybif+dyPDZw9lTWHzgdZHIRVRk3P1FoCtQA3iE0I3/vwOHAd3cfWXUE4pIhUhvk87CCxZSN6Uuz7z/DIOfG8zugt1Bx5JKLuLvybj7y+5+BlCf0PD79dz9dHd/OerpRKRCnXHoGSwZtoQGNRow44MZ9J/en7w9eUHHkkosoiJjZjeb2cEA7r7T3b8KTwyGmTU3s5tjEVJEKs4pLU9h2fBlHFDzAOaumct5z55HXoEKjeyfSM9kbiF09lKSg9E9GZEqoePBHckekc2BtQ/khbUvcMOqG9ixe0fQsaQSirTIGKGb/SVpCfxQvjgiEi+OO+g4ckbk0KxOM97a+hbnPHMO2/K3BR1LKpl9FhkzG2Fmy81sOaEC8+je5SKvV4DJQG6sA4tIxTm66dHkjszlwJQDyVmfQ8bkDH7K+ynoWFKJlOVMZgehL2B+T+hM5sciy3tfnwH3AGNiE1NEgtLuwHb86/h/cUj9Q3j5y5c58+kz2bpra9CxpJLY5wCZ7j4dmA5gZhOA29x9XayDiUj8aFGrBStGrSB9YjpvbHyDrpO6snjoYhrXbhx0NIlzkX5PZpQKjEj11Lpha3JH5tL2gLa8/fXbdJnUhe+2fxd0LIlzkT7C3H5fr1gFFZHgtWrQihWjVtCucTve+/Y90iam8c22b4KOJXEsovlkgFWU/nTZXon7mUVEKoGD6x1M7shcuk7qyurvVpOamcry4ctpUb9F0NEkDkX6CHM60KXYqx/wOKGh/ntHNZ2IxKVmdZuRPSKb45odx8fff0xqZipf/PhF0LEkDkV6Tya3hNdsd78MeAYYEJuYIhJvmtRpwvIRy+nYvCOf/vApnSd0Zt0PumUrvxTx2GW/IhudyYhUK41qNWLp8KWc0vIUPv/xc1IzU/nk+0+CjiVxJJpF5hxgaxS3JyKVQMOaDVk8dDG/a/U7Nvy0gdTMVD787sOgY0mciOjGv5lNK6E5BfgNcARwfTRCiUjlUq9GPRZesJBzp55L9vps0iamsXTYUo5pdkzQ0SRgkZ7JNCnhVQN4ETjX3e+ObjwRqSzqpNRh/pD5nNX2LDZt30T6xHTe+fqdoGNJwCI6k3H39FgFEZHKr3ZybeYMmkP/af15/pPn6TKpC4uHLua3LX4bdDQJSDTvyYiIUDOpJjMHzqTPb/qwdddWuj3djVe+fCXoWBKQfZ7JmNlK9v0FzP/j7ieVK5GIVHopiSlM6z+NobOGMm31NLpP7s7zQ56n86Gdg44mFawsl8tWE0GRCYqZ9SH0hFt94Cl3X2xmdYB/A/lAjrtPCTCiSLWSnJjMlL5TSE5IZsr7U+gxpQdzB82l62Fdg44mFagsozCPjHUIMxsP9AQ2uXuHIu0ZwIOEhqp50t3/8Ss5ZwOzzewA4D5gMdAXmOHu88zsWUBFRqQCJSUkMbHPRJITk8l8N5OeU3sye+Bsuh/ePehoUkH2656MmR1sZv3M7GIz62tmB5czRyaQUWwficA4oAfQHhgcHoTzGDObX+zVtMhHbwx/DkKzdX4Z/ntBOTOKyH5ITEjkqV5PcUnHS9i1Zxe9snoxb828oGNJBTH3sl8JC//ifxi4mF8OhFlAaPyyK929cL+CmLUG5u89kzGzU4Gx7t49vHwdgLvfVcrnDfgHsMTdl4bbhgE/uPt8M8ty90ElfG4M4cnWmjVr1jErK6vMmbdt20bdunXL3skqpLr2Xf3ef+7Ow58+zKyNs0i0RG4+6mY6N4nvezQ63mWTnp7+lrt3KvFNdy/zC7gd2AX8FWhF6DsyrcLLO4FbI9lesW23BlYVWe5P6BLZ3uVhwCO/8vk/AG8BjwGXhtvqABOAR4EL9pWhY8eOHons7OyI1q9Kqmvf1e/yKSws9KsXXe2MxRP/nuhZ72dFZbuxouNdNsCbXsrv1UiH+h8O3Oju9xVp+wK418w8/Iv+5gi3GRXu/hDwULG27cCoIPKIyP8yM+49815qJNbgzpfuZMjMIeQX5DPsuGFBR5MYifSeTFPgvVLeey/8frRsBA4pstwy3CYilZiZcXuX2/l72t8p9EJGzB7B+HfGBx1LYiTSIvMx8D/3NcIGAWvKF+cXVgJHmFkbM0sJb39uFLcvIgExM25OvZm7ut6F44yeO5pHVz4adCyJgUgvl90OZJlZK2AG8C2hs5fzCU1oVloB+lVmNhVIAw40sw3ALe7+lJldASwi9JDBeHdfvT/bF5H4dO3vriUlMYWrF1/N7xf8nvyCfP54yh+DjiVRFOnYZdPMbCvwd0LfX0kGdhO64Z7h7kv2J4S7Dy6lfQGwYH+2KSKVw59P/TMpiSlc+cKVXLXoKvIK8rjm9GuCjiVREumZDO6+GFhsZgnAgcBm38/HlkVEAK446QpSElO4dP6l/G3p38jbk8dNqTcFHUuiYL8HyHT3Qnff5O6FZtYwiplEpBoa03EME3pPwDBuzrmZm5bftPfrCVKJRVRkzOwyM7umyPLx4Xso35vZW2bWMuoJRaTaGHH8CCb3nUyiJXL7i7dz7dJrVWgquUjPZK4Efiqy/BDwFXBBeFulji0mIlIWQ44ZwtR+U0lKSOKeV+7hT4v+pEJTiUV6T6YV4ceUzawJcDrQ1d1zzCwfeCTK+USkGjr/6PNJSUzh/Onn8+DrD5JfkM8jZz9CgmkKrMom0iOWB6SE/54O7CA09TLAFqBhdGKJSHXX+ze9mT1oNjUSa/Dom49yybxLKNQzRpVOpEXmDeByMzua0BAyC9197+jGhxG6dCYiEhVnH3E28wbPo1ZSLZ5850lGzRlFQaEGVK9MIi0yVwNHA+8TGvLlhiLvDQRejlIuEREAzmx7JgsuWECd5DpM+u8khs4ayu6C3UHHkjKKqMi4+wfu3hZoArR294+LvP2X8EtEJKrSWqexaOgi6qXUI2tVFoOeG0R+QX7QsaQMIr6LFh5HrB/whJk9b2ZPmtnFwBp3/y7qCUVEgNNbnc6SYUtoUKMBMz+cSf9p/cnbkxd0LNmHSL8ncxTwCaGZJzsQmqysQ3h5rZm1j3pCEZGwk1uezPIRy2lUqxHzPp5Hn2f7sHP3zqBjya+I9EzmceBHoK27n+Luvdz9FOBwYCuhCcNERGLmxOYnkj0imwNrH8jCtQs5d+q5bM/fHnQsKUWkRaYTcLO7f1G0Mbx8C/DbaAUTESnNsc2OJWdEDs3qNGPZZ8s4+5mz+Tnv56BjSQkiLTLrgZqlvFeT0CyZIiIxd3TTo8kdmcvB9Q5mxecryJiSwY+7fgw6lhQTaZG5FrjdzE4u2mhmpwC3AX+LVjARkX1pd2A7VoxcQasGrXjly1c48+kz+WHnD0HHkiL2WWTMbKWZvWFmbxD6Xkx94BUz+9rM/mtmXxP6fkx94PrYxhUR+aW2jdqyYuQK2jRsw8qvVtJ1Ulc279gcdCwJK8vYZasBL7YsIhI3Dm14KCtGraDLxC688807dJnYhaXDl9K0TtOgo1V7+ywy7j6yrBszs+RypRER2U8t67ckd2QuXSZ14f1N75OWmcay4ctoXq950NGqtXIPaWohXc3sSeDbKGQSEdkvzes1J2dEDh2aduDDzR+SmpnKhp82BB2rWtvvImNmp5jZg8BGYDHQG5garWAiIvujWd1mZI/I5viDjueTLZ+QmpnK51s/DzpWtRXpN/6PMbM7zWwdoZv9Y4BmwJ+B5u5+eQwyiohE5MDaB7Js+DI6HdyJdT+so3NmZz7d8mnQsaqlsjxddpiZ3WBmq4B3CY3EvBoYDhwBGPCOu++JZVARkUg0qtWIpcOWcmrLU/nixy9IzUzl4+8/3vcHJarKciazFriV0LTLlwAHufu57j4F0FdsRSRuNajZgEVDF9H50M5s/HkjqZmpfPDdB0HHqlbKUmQ+J3S20gFIA04zs0inbRYRCUS9GvVYMGQBXdp04Ztt35CWmcZ7374XdKxqY59Fxt3bAKcBmUBXYB7wrZk9EV720j8tIhK8Oil1mD94Pt3bdue7Hd+RPjGdt79+O+hY1UKZbvy7+2vu/gegBXAWMJvQnDIzwqtcbGadYpJQRCQKaiXXYvag2fQ8sidbdm6h66SuvLHxjaBjVXmRzoxZ6O5L3X00oafKzgOmhf983cw+jEFGEZGoqJlUk+cGPEffo/qydddWuk3qxstfaNb4WNrv78m4+253n+Pug4GmwDBCE5qJiMStlMQUsvplMfDogfyc/zPdJ3cnd31u0LGqrHJ/4x/A3Xe4+zPu3isa2xMRiaXkxGQm953MsGOHsX33dnpM6cHSdUuDjlUlRaXIiIhUNkkJSUzoPYHRJ4xm556d9HymJy988kLQsaqcKlNkzKyPmT1hZs+a2VlF2uuY2Ztm1jPIfCISfxITEnn83Me5rNNl5BXk0efZPsxdMzfoWFVKXBQZMxtvZpvCowoUbc8wszVmttbMrv21bbj7bHe/GLgUGFjkrb8RejhBROR/JFgC484exx9P/iP5Bfn0m9aP5z54LuhYVUZcFBlC38HJKNpgZonAOKAH0B4YbGbtw+OnzS/2KjppxI3hz2FmZwIfAJsqohMiUjmZGQ90f4C/nvZX9hTuYeCMgWStygo6VpUQF9/cd/cVZta6WPNJwFp3XwdgZllAb3e/C/ifS19mZsA/gBfcfe+3rNKAOoSK1E4zW+DuhbHphYhUZmbG3d3upkZiDW5/8XYumHkB1xx5DWmkBR2tUouLIlOKFsCXRZY3ACf/yvpXAt2ABmZ2uLs/5u43AJjZSGBzSQXGzMYQGk2aZs2akZOTU+aA27Zti2j9qqS69l39rvq6JnRlY+uNTFg/gbvX3M2eqXs4p/k5QceqUFE93u4eFy+gNbCqyHJ/4Mkiy8OAR2KZoWPHjh6J7OzsiNavSqpr39Xv6uPul+52xuKMxce9MS7oOBUq0uMNvOml/F6Nl3syJdkIHFJkuWW4TUQk5q45/RoubxuaIuvyBZfzr9f+FWygSiqei8xK4Agza2NmKcAgQM8WikiF6d+yP+POHgfAnxb9ibtfujvgRJVPXBQZM5sKvAq0M7MNZjbaQ5OgXQEsAj4Eprn76iBzikj18/vf/p4nzn0Cw7h22bXclntb0JEqlbi48e+h8c9Kal8ALKjgOCIiv3DRiReRkpjCqDmjuDnnZvIK8rgt/TZCD7XKr4mLMxkRkXg3/LjhTOk7hURL5I4X7+CaJdfsfShJfoWKjIhIGQ3qMIhn+z9LUkIS9716H1ctvEqFZh9UZEREItCvfT+eG/AcKYkpPPTGQ1z2/GUU6jvepVKRERGJUK92vZgzaA41Emvwn7f+w0VzL6KgsCDoWHFJRUZEZD9kHJ7B/CHzqZVUiwnvTmDknJHsKdwTdKy4oyIjIrKfuh3WjRcueIE6yXWY/N5kLph5AbsLdgcdK66oyIiIlENq61QWDV1EvZR6TFs9jYEzBpJfkB90rLihIiMiUk6ntzqdpcOX0qBGA2Z9NIt+0/qxa8+uoGPFBRUZEZEoOKnFSSwfsZxGtRox/+P59M7qzc7dO4OOFTgVGRGRKDmx+Ylkj8imSe0mLP50Mec8cw7b87cHHStQKjIiIlF0bLNjyRmZw0F1DyJ7fTY9pvTg57yfg44VGBUZEZEoa9+kPbkjc2lRrwUvfvEiZ00+ix93/Rh0rECoyIiIxMCRjY8kd2QurRq04rUNr9Ht6W5s2bkl6FgVTkVGRCRG2jZqy4qRKzjsgMN486s36TqpK5t3bA46VoVSkRERiaFDGx5K7shcjmx8JO9+8y5pmWl8u+3boGNVGBUZEZEYa1m/JTkjcmjfpD2rv1tN2sQ0vvr5q6BjVQgVGRGRCtC8XnOyR2RzbLNj+WjzR6RmpvLlj18GHSvmVGRERCpI0zpNWT58OSccdAJrt6wlNTOV9VvXBx0rplRkREQqUOPajVk2fBkntTiJz7Z+RucJnfl0y6dBx4oZFRkRkQp2QK0DWDJsCacdchpf/vQlnTM7s2bzmqBjxYSKjIhIAOrXqM+ioYtIPTSVr37+itTMVD747oOgY0WdioyISEDqptRlwQUL6HZYN77d/i2pman895v/Bh0rqlRkREQCVDu5NnMHzaXH4T3YvGMzXSZ14e2v3w46VtSoyIiIBKxWci1mDZxFr3a92LJzC10mduH1Da8HHSsqVGREROJAjaQaTD9/Ov2O6sePeT9y5tNn8tIXLwUdq9xUZERE4kRKYgpZ/bMY1GEQP+f/TMbkDHLW5wQdq1xUZERE4khSQhKTz5vM8OOGs333ds6ecjZLPl0SdKz9piIjIhJnEhMSmdB7AhedcBE79+zk3KnnsuCTBUHH2i8qMiIicSjBEvjPuf/h8t9eTl5BHn2y+jDnozlBx4qYioyISJxKsAQe7vEwfzrlT+wu3E3/6f2Zvnp60LEikhR0gGgxsz7AOUB94Cl3X2xmrYCHgC3Ax+7+jwAjiohEzMz451n/JCUxhbtfvptBzw1id+FuhhwzJOhoZRIXZzJmNt7MNpnZqmLtGWa2xszWmtm1v7YNd5/t7hcDlwIDw83HADPc/ULghJiEFxGJMTPjrq53cVPnmyj0QobOHMrEdycGHatM4qLIAJlARtEGM0sExgE9gPbAYDNrb2bHmNn8Yq+mRT56Y/hzAK8Bo81sObAw5r0QEYkRM+PW9Fu5Lf02HGfUnFE88dYTQcfap7i4XObuK8ysdbHmk4C17r4OwMyygN7ufhfQs/g2zMyAfwAvuPveMRlGAbeEtz8DmBCrPoiIVIQbO99IjcQaXLP0GsbMH0N+QT6Xn3R50LFKZe4edAYAwkVmvrt3CC/3BzLc/aLw8jDgZHe/opTP/wEYAawE3nX3x8ysAzAW2Axsc/e/lPC5McAYgGbNmnXMysoqc+Zt27ZRt27dMq9flVTXvqvf1Us893vGhhmM+zR00eaywy5jwCEDorbtSPudnp7+lrt3KvFNd4+LF9AaWFVkuT/wZJHlYcAjsczQsWNHj0R2dnZE61cl1bXv6nf1Eu/9fnTlo85YnLH4nSvujNp2I+038KaX8ns1Xu7JlGQjcEiR5ZbhNhERAS7tdClP9XoKw7h++fX8Pefve/9THjfiucisBI4wszZmlgIMAuYGnElEJK5ceMKFTDpvEgmWwNjcsdy4/Ma4KjRxUWTMbCrwKtDOzDaY2Wh33wNcASwCPgSmufvqIHOKiMSjoccO5Zm+z5Boidz50p1cs+SauCk08fJ02eBS2hcAlXPAHhGRCjSww0CSE5MZNGMQ9716H3kFeTyY8SChB2+DExdnMiIiUn59j+rLcwOeIyUxhYffeJhL519KoRcGmklFRkSkCjm33bnMGTSHmkk1efztxxk9dzQFhQWB5VGRERGpYjIOz2D+4PnUSqpF5ruZDJ89nD2FewLJoiIjIlIFdT2sKwuHLqRuSl2eef8Zhjw3hN0Fuys8h4qMiEgV1fnQziweupj6Neoz/YPpnD/9fPL25FVoBhUZEZEq7NRDTmXpsKU0rNmQOWvm0HdaX3bt2VVh+1eRERGp4n7b4rcsH76cxrUas+CTBfSa2osdu3dUyL5VZEREqoETmp9A9ohsmtZpypJ1S+j5TE+252+P+X5VZEREqoljmh1DzogcmtdtTvb6bDKmZPBT3k8x3aeKjIhINXJUk6PIHZlLy/oteemLl+g+uTtbd22N2f5UZEREqpkjGh/BipErOLTBoby24TW6TerGlp1bYrIvFRkRkWqozQFtWDFqBW0PaMtbX79Fl4ld2Lxjc9T3oyIjIlJNtWrQityRubRr3I5aybWokVgj6vuIi1GYRUQkGC3qtyBnZA41k2pSr0a9qG9fRUZEpJo7qO5BMdu2LpeJiEjMqMiIiEjMqMiIiEjMqMiIiEjMqMiIiEjMqMiIiEjMqMiIiEjMqMiIiEjMqMiIiEjMqMiIiEjMqMiIiEjMqMiIiEjMmLsHnSFumNl3wOcRfORAIPoTMFQO1bXv6nf1on6XzaHu3qSkN1RkysHM3nT3TkHnCEJ17bv6Xb2o3+Wny2UiIhIzKjIiIhIzKjLl83jQAQJUXfuuflcv6nc56Z6MiIjEjM5kREQkZlRk9pOZZZjZGjNba2bXBp0nVszsEDPLNrMPzGy1mf0x3N7IzJaY2SfhPw8IOmssmFmimb1jZvPDy23M7PXwcX/WzFKCzhhtZtbQzGaY2Udm9qGZnVodjreZ/Sn8b3yVmU01s5pV9Xib2Xgz22Rmq4q0lXiMLeSh8M/gPTM7MZJ9qcjsBzNLBMYBPYD2wGAzax9sqpjZA1zt7u2BU4DLw329Fljm7kcAy8LLVdEfgQ+LLN8NPODuhwM/AKMDSRVbDwIL3f03wHGE+l+lj7eZtQD+AHRy9w5AIjCIqnu8M4GMYm2lHeMewBHh1xjg0Uh2pCKzf04C1rr7OnfPB7KA3gFnigl3/9rd3w7//WdCv3BaEOrvxPBqE4E+gQSMITNrCZwDPBleNqALMCO8SpXrt5k1ADoDTwG4e767b6UaHG8gCahlZklAbeBrqujxdvcVwJZizaUd497AJA95DWhoZs3Lui8Vmf3TAviyyPKGcFuVZmatgROA14Fm7v51+K1vgGZB5YqhfwHXAIXh5cbAVnffE16uise9DfAdMCF8mfBJM6tDFT/e7r4RuA/4glBx+RF4i6p/vIsq7RiX6/edioyUiZnVBZ4DrnL3n4q+56FHFKvUY4pm1hPY5O5vBZ2lgiUBJwKPuvsJwHaKXRqrosf7AEL/Y28DHAzU4X8vJ1Ub0TzGKjL7ZyNwSJHlluG2KsnMkgkVmCnuPjPc/O3eU+bwn5uCyhcjpwO9zGw9ocuhXQjdq2gYvpwCVfO4bwA2uPvr4eUZhIpOVT/e3YDP3P07d98NzCT0b6CqH++iSjvG5fp9pyKzf1YCR4SfPEkhdINwbsCZYiJ8H+Ip4EN3v7/IW3OBEeG/jwDmVHS2WHL369y9pbu3JnR8l7v7BUA20D+8WlXs9zfAl2bWLtzUFfiAKn68CV0mO8XMaof/ze/td5U+3sWUdoznAsPDT5mdAvxY5LLaPunLmPvJzM4mdM0+ERjv7ncEmyg2zOx3wIvA+/z/exPXE7ovMw1oRWjk6gHuXvxGYpVgZmnAX9y9p5kdRujMphHwDjDU3fMCjBd1ZnY8oYcdUoB1wChC/yGt0sfbzP4ODCT0ROU7wEWE7j1UueNtZlOBNEKjLX8L3ALMpoRjHC66jxC6fLgDGOXub5Z5XyoyIiISK7pcJiIiMaMiIyIiMaMiIyIiMaMiIyIiMaMiIyIiMaMiI4Ews7Fm5ma2qIT3ZphZTgVmSQtn6VBR+4yEmR1lZi+a2fZwztalrLfezO4rsjzAzEZWVM4i+00JH9/ji7W3DufvWdGZJDgqMhK0s8zst0GHiHP3Ag2BXsCphMbWKosBwMjYRPpVKYS+d3F8sfavCeV/qaIDSXCS9r2KSMxsITQ8xQ1UkdFtS2JmNd19Vzk28Rtgrrsvi1amSIW/kFejPP0If4nxteilkspAZzISJAfuIDRG2DGlrRS+9LK5hHY3syuKLK83s/vM7Foz+9rMfjSzf4aHwzg7PCHVz2Y2u5RJtw42s/nhy1JfmNmlJezzDDPLNbMdZva9mT1hZvWKvD8ynOskM8sxs53AX3+lb8eb2bLw9n4wsylm1iz8Xmszc6At8KfwdnNK21ax7WYC/YDU8OfczMYWeb+3mb1pZrvM7Bszuyc8Rt3e98ea2WYz+52ZrQR2AeebWR0ze8RCE/btMLPPzGycmdUvsvufw39OKLLv1iVdLrPQpHBjwz/vvPAxGlK8L+GsZ1po0qztZvaSmR1dbL3RFppcb2c4e27xdaTiqchI0KYDnxA6m4mGQYTm+xkF3AP8GbgfuA24CbgUSAXuKuGzTwHvAX2BBcCjxX4hng4sJTQMen/gKuBsYEIJ25oKzAu/P7+koGbWBMghNHfJEODKcLYlFhoTb+/lpW+AZ8J//30ZfgaE+5tNaCiUU8OvvfPiDCA0AOQbhC7B/Z3QZFTFfya1Cc0r8iShIUXeCLclEjpePQj9TLsQOo57dQn/eXuRfZd2ie/W8LYeD2d5GZhiZoOLrdeK0GXDO4DBQFPg2fAZFmbWGXgMeDqc60LgFaBBaT8gqSDurpdeFf4CxgKbw38fCRQAR4aXZwA5Ja1bbBsOXFFkeT2wFkgs0vYGobGo2hRpuwf4tshyWnhbjxfb/hLgtSLLLwLZxdbpEv5shyJ9ceCPZfgZ/APYCtQv0nZy+PODi/XrvjJs7xfrFf85htuM0LhUE4q1XwjsBBoX+Zk70Hsf+0wiNFqxA63CbXXDyyOLrds63N4zvNyI0FQCtxRbbwGwpshyZvgYHlGkrU94W78JL/8FeCvof9d6/e9LZzISDyYTGgX3uihsK8fdC4osrwXWu/tnxdqa2P/O1z6r2PJMoGP4kk5tQv8jn2ZmSXtfhG5i7wY6Fvvs82XIehKw2IvMz+OhIfbXA78rw+f3x5GEzgqK92M5UBMo+oSdAy8U34CZDbPQhGbbCPV97438IyPM0oHQmdH0Yu3PAkeGz/T2Wu/unxRZ/iD8Z8vwn+8CJ5jZA2bWuYRjKwFRkZHAeWjmwXuAoWZ2aDk3t7XYcn4pbUboKaiiis+RsonQ/9QPBA4gdJno34R+se595QHJ/HK+DQiNbLsvzUtZ71tC/8uPhQPDfy7gl/3YW4SL9uMHD00v/n/M7DxgEvAqcD5wCnBe+O2aEWbZO4Vv8Z/B3uWiP4OtxdbZm6smgLsvJXSJtDOhS5Cbw/eK6kSYSaJMT5dJvBgP3Aj8rYT3dlGsIJRy4768mpawvAfYTOiXmRO6jLSghM9+VWy5LMObf13CPiE07W2sZuTcOzz/GEL3a4oresZXUh/OB1539/+7N2RmqfuZZe99mqbA90Xa9077G9FUAu4+EZgYPgPqCzxA6CGEa3/1gxJTKjISF9w9z0JfJLyL0C/Y3UXe3gDUM7MWHpqLHeCsGMQ4j19eHjqP0HX+AmC7mb0GtHP3W6O0v9eBy8ysnrv/DGCh7wy1JjrfJcnnf88u1hB6bLy1uz+xH9usRejsragLStgvJey7uFWE5ic5n9ADAHsNAD529+/2Ix/hz/3HzPoC7fdnGxI9KjIST/5DaEK004DcIu0LCd2UHm9m/yQ0D/v/PF4cBT3M7I7wvvsCZxKa932va4BlZlZI6Kb6z4Tub5wD3ODuH0e4v/uBy4BFZnY3oRvm/yA0Qdxz5elI2EdAbzPrQ6hQf+XuX5nZ1cDT4ceOXyBUFA4jdDO9v7vv+JVtLgHGmdkNhIrk2YRmkfw/7p5vZp8BA8xsFaEz0feKb8hDE2L9C7jRzPYAbxL6uZ9N6AmyMrPQhGONCF8qA04g9KSezmICpnsyEjfCv9weKKF9M6HvfLQkNHvfUEKP/EbbRYTms58N9AQud/f/m1bb3V8idM2/CaFHZecRKjxfUrZ7ML8Q/h93OqFfwlOBcYSeYDuz+L2Q/fRvYDGhS5ErCV0iw92fJVQ8jyd0030moUej3+b/n4WU5j/AP4E/hj93KCUfi0sJ3f9ZGt73waVs72ZCZ6+XEXrUuzOh2SezytC/olYSOmt5DFgU3t5Y4MEItyNRppkxRUQkZnQmIyIiMaMiIyIiMaMiIyIiMaMiIyIiMaMiIyIiMaMiIyIiMaMiIyIiMaMiIyIiMaMiIyIiMfP/AM3WLzlTtVfBAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": { "image/png": { "height": 376, "width": 409 }, "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "s=time.time()\n", "#Array of iterations\n", "Niter = np.arange( 1, 100 )\n", "\n", "plt.figure( figsize=(6,6) )\n", "#plt.semilogy( Niter, 2.0**-Niter, color=\"green\", lw = 2 )\n", "plt.semilogy( Niter, np.exp2(-Niter), color=\"green\", lw = 2 )\n", "plt.grid(True)\n", "plt.xlabel(\"Number of Iterations\",size=15)\n", "plt.ylabel(\"Absolute Error $|p_n-p|$\",size=15)\n", "f=time.time()-s" ] }, { "cell_type": "markdown", "metadata": { "cell_id": "00019-b2770622-2ad4-402d-a96e-536b3ad56f04", "colab_type": "text", "deepnote_cell_type": "markdown", "id": "pd9PwTe1w7V-", "jp-MarkdownHeadingCollapsed": true, "tags": [] }, "source": [ "### Example 1" ] }, { "cell_type": "markdown", "metadata": { "cell_id": "00020-b74fc684-6ecc-47f5-9692-7bc288468de2", "colab_type": "text", "deepnote_cell_type": "markdown", "id": "vlEpT9Sjw7V_" }, "source": [ "Find the root of the function $f$\n", "\n", "$f(x) = x^3 - 2$\n", "\n", "for $20$ iterations, show the result and the relative error in each iteration." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1.2599210498948732" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "2**(1/3)" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "cell_id": "00021-b8212579-3028-497b-9892-c66c38ac1ec8", "colab": {}, "colab_type": "code", "deepnote_cell_type": "code", "id": "GpdA20yow7WB" }, "outputs": [], "source": [ "#Defining Bisection function\n", "def Bisection( f, a, b, Nmax, printer=False ):\n", " \"\"\"\n", " Find the root of function f between a and b\n", " \n", " f: function\n", " a,b: the initial interval\n", " Nmax: Nmax de interations\n", " printer: bool, to print internal steps\n", " \"\"\"\n", " #verifying the STEP1, a and b with different signs\n", " if f(a)*f(b)>0:\n", " raise Exception(\"Error, f(a) and f(b) should have opposite signs\")\n", " #Assigning the current extreme values, STEP2\n", " ai = a\n", " bi = b\n", " #Iterations\n", " n = 1\n", " while n<=Nmax:\n", " #Bisection, STEP3\n", " pi = (ai+bi)/2.0\n", " #Evaluating function in pi, STEP4 and STEP5\n", " if printer:\n", " print( f\"Value for {n} iterations: {pi}\" )\n", " #Condition A\n", " if f(pi)*f(ai)>0:\n", " ai = pi\n", " #Condition B\n", " else:\n", " bi = pi\n", " #Condition C: repeat the cycle\n", " n+=1\n", " #Final result\n", " return pi" ] }, { "cell_type": "markdown", "metadata": { "cell_id": "00022-67251c83-642f-4f07-9fce-0f5b8de62694", "colab_type": "text", "deepnote_cell_type": "markdown", "id": "Z1XWelBgw7WJ" }, "source": [ "Defining function" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "cell_id": "00023-2c8a0a16-ef64-4440-9965-c4baa27f7f20", "colab": {}, "colab_type": "code", "deepnote_cell_type": "code", "id": "0Vd7ahTbw7WK" }, "outputs": [], "source": [ "def function(x):\n", " return x**3.0 - 2.0" ] }, { "cell_type": "markdown", "metadata": { "cell_id": "00024-cda4b3e4-5d95-489e-b7cc-903a1a562024", "colab_type": "text", "deepnote_cell_type": "markdown", "id": "Cy9OT7XDw7WN" }, "source": [ "Finding the root of the function. The real root is $\\sqrt[3]{2}$, so `a` and `b` should enclose this value" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "cell_id": "00025-245deaf7-fc02-41aa-bb3c-198220f3026e", "colab": {}, "colab_type": "code", "deepnote_cell_type": "code", "id": "pvl1fZSKw7WO", "outputId": "e8accba7-9299-4533-c81c-5be42aeb1869" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Value for 1 iterations: 1.0\n", "Value for 2 iterations: 1.5\n", "Value for 3 iterations: 1.25\n", "Value for 4 iterations: 1.375\n", "Value for 5 iterations: 1.3125\n", "Value for 6 iterations: 1.28125\n", "Value for 7 iterations: 1.265625\n", "Value for 8 iterations: 1.2578125\n", "Value for 9 iterations: 1.26171875\n", "Value for 10 iterations: 1.259765625\n", "Value for 11 iterations: 1.2607421875\n", "Value for 12 iterations: 1.26025390625\n", "Value for 13 iterations: 1.260009765625\n", "Value for 14 iterations: 1.2598876953125\n", "Value for 15 iterations: 1.25994873046875\n", "Value for 16 iterations: 1.259918212890625\n", "Value for 17 iterations: 1.2599334716796875\n", "Value for 18 iterations: 1.2599258422851562\n", "Value for 19 iterations: 1.2599220275878906\n", "Value for 20 iterations: 1.2599201202392578\n", "Real value: 1.2599210498948732\n", "Absolute error 9.296556153781665e-07\n" ] } ], "source": [ "#Defining a and b\n", "a = 0.0\n", "b = 2.0\n", "Nmax = 20\n", "result = Bisection(function, a, b, Nmax, printer=True)\n", "print ( \"Real value:\", 2**(1/3))\n", "print ( \"Absolute error\", abs((2**(1/3)-result)) )" ] }, { "cell_type": "markdown", "metadata": { "cell_id": "00026-fdc43ff7-e829-4000-9ebf-f75e1ee04c29", "colab_type": "text", "deepnote_cell_type": "markdown", "id": "CmP1u-lPw7WQ" }, "source": [ "Using the error analysis, we can predict the produced error at $20$ iterations by computing:\n", "\n", "$$ \\left( \\frac{1}{2^{20}}\\right) \\approx 9.53674316\\times 10^{-7} $$\n", "\n", "This value is very close to the obtained relative error.\n", "\n", "If we were interested in a double precision, i.e. $\\epsilon \\sim 10^{-17}$, the number of required iterations would be:\n", "\n", "$$ 10^{-17} = \\left( \\frac{1}{2^{N}}\\right) \\longrightarrow N = \\frac{17}{\\log_{10}(2)} \\approx 56 $$" ] }, { "cell_type": "markdown", "metadata": { "cell_id": "00027-54ddeed1-55b9-4a2a-9ad7-d35b738d3e52", "colab_type": "text", "deepnote_cell_type": "markdown", "id": "p2lZUHGSw7WR" }, "source": [ "Now the same but using scipy implementation" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Value for 1 iterations: 1.0\n", "Value for 2 iterations: 1.5\n", "Value for 3 iterations: 1.25\n", "Value for 4 iterations: 1.375\n", "Value for 5 iterations: 1.3125\n", "Value for 6 iterations: 1.28125\n", "Value for 7 iterations: 1.265625\n", "Value for 8 iterations: 1.2578125\n", "Value for 9 iterations: 1.26171875\n", "Value for 10 iterations: 1.259765625\n", "Value for 11 iterations: 1.2607421875\n", "Value for 12 iterations: 1.26025390625\n", "Value for 13 iterations: 1.260009765625\n", "Value for 14 iterations: 1.2598876953125\n", "Value for 15 iterations: 1.25994873046875\n", "Value for 16 iterations: 1.259918212890625\n", "Value for 17 iterations: 1.2599334716796875\n", "Value for 18 iterations: 1.2599258422851562\n", "Value for 19 iterations: 1.2599220275878906\n", "Value for 20 iterations: 1.2599201202392578\n", "Value for 21 iterations: 1.2599210739135742\n", "Value for 22 iterations: 1.259920597076416\n", "Value for 23 iterations: 1.2599208354949951\n", "Value for 24 iterations: 1.2599209547042847\n", "Value for 25 iterations: 1.2599210143089294\n", "Value for 26 iterations: 1.2599210441112518\n", "Value for 27 iterations: 1.259921059012413\n", "Value for 28 iterations: 1.2599210515618324\n", "Value for 29 iterations: 1.2599210478365421\n", "Value for 30 iterations: 1.2599210496991873\n", "Value for 31 iterations: 1.2599210506305099\n", "Value for 32 iterations: 1.2599210501648486\n", "Value for 33 iterations: 1.259921049932018\n", "Value for 34 iterations: 1.2599210498156026\n", "Value for 35 iterations: 1.2599210498738103\n", "Value for 36 iterations: 1.259921049902914\n", "Value for 37 iterations: 1.2599210498883622\n", "Value for 38 iterations: 1.2599210498956381\n", "Value for 39 iterations: 1.2599210498920002\n", "Value for 40 iterations: 1.2599210498938191\n", "Value for 41 iterations: 1.2599210498947286\n", "Value for 42 iterations: 1.2599210498951834\n", "Value for 43 iterations: 1.259921049894956\n", "Value for 44 iterations: 1.2599210498948423\n", "Value for 45 iterations: 1.2599210498948992\n", "Value for 46 iterations: 1.2599210498948707\n", "Value for 47 iterations: 1.259921049894885\n", "Value for 48 iterations: 1.2599210498948779\n", "Value for 49 iterations: 1.2599210498948743\n", "Value for 50 iterations: 1.2599210498948725\n", "Value for 51 iterations: 1.2599210498948734\n", "Value for 52 iterations: 1.259921049894873\n", "Value for 53 iterations: 1.2599210498948732\n", "Value for 54 iterations: 1.259921049894873\n", "Value for 55 iterations: 1.259921049894873\n", "Value for 56 iterations: 1.259921049894873\n", "Real value: 1.2599210498948732\n", "Absolute error 2.220446049250313e-16\n" ] } ], "source": [ "#Defining a and b\n", "a = 0.0\n", "b = 2.0\n", "Nmax = 56\n", "result = Bisection(function, a, b, Nmax, printer=True)\n", "print ( \"Real value:\", 2**(1/3))\n", "print ( \"Absolute error\", abs((2**(1/3)-result)) )" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "cell_id": "00028-c539015e-ae21-43d8-89d6-8cc433cb4a58", "colab": {}, "colab_type": "code", "deepnote_cell_type": "code", "id": "LKwXO9vfw7WS" }, "outputs": [], "source": [ "from scipy import optimize" ] }, { "cell_type": "markdown", "metadata": { "cell_id": "00029-fb204e73-8727-4e9c-b600-2eedd82dbf10", "colab_type": "text", "deepnote_cell_type": "markdown", "id": "x92AU0K6w7WV" }, "source": [ "See help with `optimize?:`" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\u001b[0;31mSignature:\u001b[0m \u001b[0mBisection\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mb\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mNmax\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mprinter\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mFalse\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;31mDocstring:\u001b[0m\n", "Find the root of function f between a and b\n", "\n", "f: function\n", "a,b: the initial interval\n", "Nmax: Nmax de interations\n", "printer: bool, to print internal steps\n", "\u001b[0;31mFile:\u001b[0m ~/ComputationalMethods/material/\n", "\u001b[0;31mType:\u001b[0m function\n" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "Bisection?" ] }, { "cell_type": "code", "execution_count": 27, "metadata": { "cell_id": "00030-3deb7c62-edc6-42ae-9ad6-b130d0cbc0b5", "deepnote_cell_type": "code" }, "outputs": [ { "data": { "text/plain": [ "\u001b[0;31mSignature:\u001b[0m\n", "\u001b[0moptimize\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbisect\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\u001b[0m\n", "\u001b[0;34m\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n", "\u001b[0;34m\u001b[0m \u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n", "\u001b[0;34m\u001b[0m \u001b[0mb\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n", "\u001b[0;34m\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n", "\u001b[0;34m\u001b[0m \u001b[0mxtol\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m2e-12\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n", "\u001b[0;34m\u001b[0m \u001b[0mrtol\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m8.881784197001252e-16\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n", "\u001b[0;34m\u001b[0m \u001b[0mmaxiter\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m100\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n", "\u001b[0;34m\u001b[0m \u001b[0mfull_output\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mFalse\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n", "\u001b[0;34m\u001b[0m \u001b[0mdisp\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n", "\u001b[0;34m\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;31mSource:\u001b[0m \n", "\u001b[0;32mdef\u001b[0m \u001b[0mbisect\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mb\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n", "\u001b[0;34m\u001b[0m \u001b[0mxtol\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0m_xtol\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mrtol\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0m_rtol\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmaxiter\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0m_iter\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n", "\u001b[0;34m\u001b[0m \u001b[0mfull_output\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mFalse\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdisp\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\n", "\u001b[0;34m\u001b[0m \u001b[0;34m\"\"\"\u001b[0m\n", "\u001b[0;34m Find root of a function within an interval using bisection.\u001b[0m\n", "\u001b[0;34m\u001b[0m\n", "\u001b[0;34m Basic bisection routine to find a zero of the function `f` between the\u001b[0m\n", "\u001b[0;34m arguments `a` and `b`. `f(a)` and `f(b)` cannot have the same signs.\u001b[0m\n", "\u001b[0;34m Slow but sure.\u001b[0m\n", "\u001b[0;34m\u001b[0m\n", "\u001b[0;34m Parameters\u001b[0m\n", "\u001b[0;34m ----------\u001b[0m\n", "\u001b[0;34m f : function\u001b[0m\n", "\u001b[0;34m Python function returning a number. `f` must be continuous, and\u001b[0m\n", "\u001b[0;34m f(a) and f(b) must have opposite signs.\u001b[0m\n", "\u001b[0;34m a : scalar\u001b[0m\n", "\u001b[0;34m One end of the bracketing interval [a,b].\u001b[0m\n", "\u001b[0;34m b : scalar\u001b[0m\n", "\u001b[0;34m The other end of the bracketing interval [a,b].\u001b[0m\n", "\u001b[0;34m xtol : number, optional\u001b[0m\n", "\u001b[0;34m The computed root ``x0`` will satisfy ``np.allclose(x, x0,\u001b[0m\n", "\u001b[0;34m atol=xtol, rtol=rtol)``, where ``x`` is the exact root. The\u001b[0m\n", "\u001b[0;34m parameter must be nonnegative.\u001b[0m\n", "\u001b[0;34m rtol : number, optional\u001b[0m\n", "\u001b[0;34m The computed root ``x0`` will satisfy ``np.allclose(x, x0,\u001b[0m\n", "\u001b[0;34m atol=xtol, rtol=rtol)``, where ``x`` is the exact root. The\u001b[0m\n", "\u001b[0;34m parameter cannot be smaller than its default value of\u001b[0m\n", "\u001b[0;34m ``4*np.finfo(float).eps``.\u001b[0m\n", "\u001b[0;34m maxiter : int, optional\u001b[0m\n", "\u001b[0;34m if convergence is not achieved in `maxiter` iterations, an error is\u001b[0m\n", "\u001b[0;34m raised. Must be >= 0.\u001b[0m\n", "\u001b[0;34m args : tuple, optional\u001b[0m\n", "\u001b[0;34m containing extra arguments for the function `f`.\u001b[0m\n", "\u001b[0;34m `f` is called by ``apply(f, (x)+args)``.\u001b[0m\n", "\u001b[0;34m full_output : bool, optional\u001b[0m\n", "\u001b[0;34m If `full_output` is False, the root is returned. If `full_output` is\u001b[0m\n", "\u001b[0;34m True, the return value is ``(x, r)``, where x is the root, and r is\u001b[0m\n", "\u001b[0;34m a `RootResults` object.\u001b[0m\n", "\u001b[0;34m disp : bool, optional\u001b[0m\n", "\u001b[0;34m If True, raise RuntimeError if the algorithm didn't converge.\u001b[0m\n", "\u001b[0;34m Otherwise the convergence status is recorded in a `RootResults`\u001b[0m\n", "\u001b[0;34m return object.\u001b[0m\n", "\u001b[0;34m\u001b[0m\n", "\u001b[0;34m Returns\u001b[0m\n", "\u001b[0;34m -------\u001b[0m\n", "\u001b[0;34m x0 : float\u001b[0m\n", "\u001b[0;34m Zero of `f` between `a` and `b`.\u001b[0m\n", "\u001b[0;34m r : `RootResults` (present if ``full_output = True``)\u001b[0m\n", "\u001b[0;34m Object containing information about the convergence. In particular,\u001b[0m\n", "\u001b[0;34m ``r.converged`` is True if the routine converged.\u001b[0m\n", "\u001b[0;34m\u001b[0m\n", "\u001b[0;34m Examples\u001b[0m\n", "\u001b[0;34m --------\u001b[0m\n", "\u001b[0;34m\u001b[0m\n", "\u001b[0;34m >>> def f(x):\u001b[0m\n", "\u001b[0;34m ... return (x**2 - 1)\u001b[0m\n", "\u001b[0;34m\u001b[0m\n", "\u001b[0;34m >>> from scipy import optimize\u001b[0m\n", "\u001b[0;34m\u001b[0m\n", "\u001b[0;34m >>> root = optimize.bisect(f, 0, 2)\u001b[0m\n", "\u001b[0;34m >>> root\u001b[0m\n", "\u001b[0;34m 1.0\u001b[0m\n", "\u001b[0;34m\u001b[0m\n", "\u001b[0;34m >>> root = optimize.bisect(f, -2, 0)\u001b[0m\n", "\u001b[0;34m >>> root\u001b[0m\n", "\u001b[0;34m -1.0\u001b[0m\n", "\u001b[0;34m\u001b[0m\n", "\u001b[0;34m See Also\u001b[0m\n", "\u001b[0;34m --------\u001b[0m\n", "\u001b[0;34m brentq, brenth, bisect, newton\u001b[0m\n", "\u001b[0;34m fixed_point : scalar fixed-point finder\u001b[0m\n", "\u001b[0;34m fsolve : n-dimensional root-finding\u001b[0m\n", "\u001b[0;34m\u001b[0m\n", "\u001b[0;34m \"\"\"\u001b[0m\u001b[0;34m\u001b[0m\n", "\u001b[0;34m\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtuple\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\n", "\u001b[0;34m\u001b[0m \u001b[0margs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\n", "\u001b[0;34m\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mxtol\u001b[0m \u001b[0;34m<=\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\n", "\u001b[0;34m\u001b[0m \u001b[0;32mraise\u001b[0m \u001b[0mValueError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"xtol too small (%g <= 0)\"\u001b[0m \u001b[0;34m%\u001b[0m \u001b[0mxtol\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\n", "\u001b[0;34m\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mrtol\u001b[0m \u001b[0;34m<\u001b[0m \u001b[0m_rtol\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\n", "\u001b[0;34m\u001b[0m \u001b[0;32mraise\u001b[0m \u001b[0mValueError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"rtol too small (%g < %g)\"\u001b[0m \u001b[0;34m%\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mrtol\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0m_rtol\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\n", "\u001b[0;34m\u001b[0m \u001b[0mr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_zeros\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_bisect\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mb\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mxtol\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mrtol\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmaxiter\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfull_output\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdisp\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\n", "\u001b[0;34m\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mresults_c\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfull_output\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mr\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;31mFile:\u001b[0m /etc/anaconda3/lib/python3.6/site-packages/scipy/optimize/zeros.py\n", "\u001b[0;31mType:\u001b[0m function\n" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "optimize.bisect??" ] }, { "cell_type": "markdown", "metadata": { "cell_id": "00031-9ab1181b-ca31-446a-a0d0-72406b969751", "colab_type": "text", "deepnote_cell_type": "markdown", "id": "VqmyDCFIw7Wa" }, "source": [ "Check the help: See [here](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.bisect.html)" ] }, { "cell_type": "code", "execution_count": 30, "metadata": { "cell_id": "00032-bdb45b3c-df79-4fd2-80b6-21dce8364047", "deepnote_cell_type": "code" }, "outputs": [], "source": [ "def f(x,n,c=2):\n", " return x**n-c" ] }, { "cell_type": "code", "execution_count": 35, "metadata": { "cell_id": "00033-45328d77-7f76-4020-9989-5b17b34ab3c9", "deepnote_cell_type": "code" }, "outputs": [ { "data": { "text/plain": [ "1.2599210498938191" ] }, "execution_count": 35, "metadata": {}, "output_type": "execute_result" } ], "source": [ "optimize.bisect(f,0,2,args=(3,)) " ] }, { "cell_type": "code", "execution_count": 36, "metadata": { "cell_id": "00034-2564a11e-2fcf-4694-914e-0eb3dcd04c2f", "deepnote_cell_type": "code" }, "outputs": [ { "data": { "text/plain": [ "1.2599210498938191" ] }, "execution_count": 36, "metadata": {}, "output_type": "execute_result" } ], "source": [ "optimize.bisect(lambda x: f(x,3),0,2)" ] }, { "cell_type": "code", "execution_count": 26, "metadata": { "cell_id": "00035-36235759-01d1-457e-9607-c28caae798f1", "colab": { "base_uri": "https://localhost:8080/", "height": 1000 }, "colab_type": "code", "deepnote_cell_type": "code", "id": "s7ZzL1BEw7Wa", "outputId": "7b4af95b-ae70-4a02-e3ea-eb5ad062a87f" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Help on function bisect in module scipy.optimize.zeros:\n", "\n", "bisect(f, a, b, args=(), xtol=2e-12, rtol=8.881784197001252e-16, maxiter=100, full_output=False, disp=True)\n", " Find root of a function within an interval using bisection.\n", " \n", " Basic bisection routine to find a zero of the function `f` between the\n", " arguments `a` and `b`. `f(a)` and `f(b)` cannot have the same signs.\n", " Slow but sure.\n", " \n", " Parameters\n", " ----------\n", " f : function\n", " Python function returning a number. `f` must be continuous, and\n", " f(a) and f(b) must have opposite signs.\n", " a : scalar\n", " One end of the bracketing interval [a,b].\n", " b : scalar\n", " The other end of the bracketing interval [a,b].\n", " xtol : number, optional\n", " The computed root ``x0`` will satisfy ``np.allclose(x, x0,\n", " atol=xtol, rtol=rtol)``, where ``x`` is the exact root. The\n", " parameter must be nonnegative.\n", " rtol : number, optional\n", " The computed root ``x0`` will satisfy ``np.allclose(x, x0,\n", " atol=xtol, rtol=rtol)``, where ``x`` is the exact root. The\n", " parameter cannot be smaller than its default value of\n", " ``4*np.finfo(float).eps``.\n", " maxiter : int, optional\n", " If convergence is not achieved in `maxiter` iterations, an error is\n", " raised. Must be >= 0.\n", " args : tuple, optional\n", " Containing extra arguments for the function `f`.\n", " `f` is called by ``apply(f, (x)+args)``.\n", " full_output : bool, optional\n", " If `full_output` is False, the root is returned. If `full_output` is\n", " True, the return value is ``(x, r)``, where x is the root, and r is\n", " a `RootResults` object.\n", " disp : bool, optional\n", " If True, raise RuntimeError if the algorithm didn't converge.\n", " Otherwise, the convergence status is recorded in a `RootResults`\n", " return object.\n", " \n", " Returns\n", " -------\n", " x0 : float\n", " Zero of `f` between `a` and `b`.\n", " r : `RootResults` (present if ``full_output = True``)\n", " Object containing information about the convergence. In particular,\n", " ``r.converged`` is True if the routine converged.\n", " \n", " Examples\n", " --------\n", " \n", " >>> def f(x):\n", " ... return (x**2 - 1)\n", " \n", " >>> from scipy import optimize\n", " \n", " >>> root = optimize.bisect(f, 0, 2)\n", " >>> root\n", " 1.0\n", " \n", " >>> root = optimize.bisect(f, -2, 0)\n", " >>> root\n", " -1.0\n", " \n", " See Also\n", " --------\n", " brentq, brenth, bisect, newton\n", " fixed_point : scalar fixed-point finder\n", " fsolve : n-dimensional root-finding\n", "\n" ] } ], "source": [ "help(optimize.bisect)" ] }, { "cell_type": "markdown", "metadata": { "cell_id": "00036-c0648f60-b13b-4d60-b251-69523c1d13fb", "colab_type": "text", "deepnote_cell_type": "markdown", "id": "yC9nyRCwGJaw", "tags": [] }, "source": [ "## General arguments and options of a SciPy function" ] }, { "cell_type": "markdown", "metadata": { "cell_id": "00037-dfcd6c0d-7ec6-4552-8d7e-e190ec9d4259", "colab_type": "text", "deepnote_cell_type": "markdown", "id": "B-D9tsnEGTHD" }, "source": [ "Let define a general multiparameter \n", "* `f`\n", "with a mandatory parameter and two optional parameters, corresponding to a straight line with the slope and the intercept as parameters\n", "$$\n", "f(x)=mx+b\\,,\n", "$$" ] }, { "cell_type": "code", "execution_count": 44, "metadata": { "cell_id": "00038-f34d5918-9956-40d0-a4f7-7579a531398a", "colab": {}, "colab_type": "code", "deepnote_cell_type": "code", "id": "3Jxl8Cyx4J32" }, "outputs": [], "source": [ "def f(x,m=3,b=-2):\n", " return m*x+b" ] }, { "cell_type": "code", "execution_count": 45, "metadata": { "cell_id": "00039-615aaa9e-ef2b-44a3-948a-6a7f0c842e04", "colab": { "base_uri": "https://localhost:8080/", "height": 35 }, "colab_type": "code", "deepnote_cell_type": "code", "id": "WrYXrulI44Ij", "outputId": "b5bb125b-60da-4edc-fc0b-0249f80b66ca" }, "outputs": [ { "data": { "text/plain": [ "0.0" ] }, "execution_count": 45, "metadata": {}, "output_type": "execute_result" } ], "source": [ "f(0.66666666666666666,-3,2)" ] }, { "cell_type": "markdown", "metadata": { "cell_id": "00040-520dc43c-1a90-4944-b852-d236ef814bfd", "deepnote_cell_type": "markdown" }, "source": [ "We have at least three ways to use a multiparameter function, the first two with the implicit `lambda` function" ] }, { "cell_type": "code", "execution_count": 46, "metadata": { "cell_id": "00041-a764ed04-9dfd-4083-9143-ead7816985a0", "colab": { "base_uri": "https://localhost:8080/", "height": 35 }, "colab_type": "code", "deepnote_cell_type": "code", "id": "R6ihC-WN4ZUR", "outputId": "9d5a4737-95fe-4940-e9bc-f4dee796793c" }, "outputs": [ { "data": { "text/plain": [ "0.6666666666671972" ] }, "execution_count": 46, "metadata": {}, "output_type": "execute_result" } ], "source": [ "optimize.bisect(lambda x: -3*x+2,-10,10)" ] }, { "cell_type": "code", "execution_count": 47, "metadata": { "cell_id": "00042-de8c5d98-e2b0-4246-ae5c-3f99002ad753", "colab": { "base_uri": "https://localhost:8080/", "height": 35 }, "colab_type": "code", "deepnote_cell_type": "code", "id": "iZHq9BrU4xFh", "outputId": "1a105e92-22a5-47e6-873b-6f17d8f6a412" }, "outputs": [ { "data": { "text/plain": [ "0.6666666666671972" ] }, "execution_count": 47, "metadata": {}, "output_type": "execute_result" } ], "source": [ "optimize.bisect(lambda x: f(x,-3,2),-10,10,)" ] }, { "cell_type": "markdown", "metadata": { "cell_id": "00043-86e4bf79-9a0c-4739-bde6-7bf886e7c524", "deepnote_cell_type": "markdown" }, "source": [ "and the recommended way to used the designed implementation in Scipy through the `args` parameter which is common to many methods " ] }, { "cell_type": "code", "execution_count": 48, "metadata": { "cell_id": "00044-68c0ee1a-085a-4429-b1b2-4225957cb1b1", "colab": { "base_uri": "https://localhost:8080/", "height": 35 }, "colab_type": "code", "deepnote_cell_type": "code", "id": "ojivSEZh5YLx", "outputId": "ed0c6690-c165-43d1-8e8e-4251fe0ea33b" }, "outputs": [ { "data": { "text/plain": [ "0.6666666666671972" ] }, "execution_count": 48, "metadata": {}, "output_type": "execute_result" } ], "source": [ "optimize.bisect(f,-10,10,args=(-3,2))" ] }, { "cell_type": "markdown", "metadata": { "cell_id": "00045-d4261383-c2e7-4949-917a-8b3a5ac9bd1f", "colab_type": "text", "deepnote_cell_type": "markdown", "id": "7wLEoVWK8FqU" }, "source": [ "* `xtol`" ] }, { "cell_type": "code", "execution_count": 49, "metadata": { "cell_id": "00046-3410e4a9-1d40-4862-8921-54aad836246f", "colab": { "base_uri": "https://localhost:8080/", "height": 35 }, "colab_type": "code", "deepnote_cell_type": "code", "id": "1UyHY0m76fIk", "outputId": "2b172be2-43dd-4940-ba04-75659bda3d83" }, "outputs": [ { "data": { "text/plain": [ "0.6666666666666671" ] }, "execution_count": 49, "metadata": {}, "output_type": "execute_result" } ], "source": [ "optimize.bisect(f,-10,10,args=(-3,2),xtol=1E-17)" ] }, { "cell_type": "markdown", "metadata": { "cell_id": "00047-56cf8e40-248f-444e-86ad-1eb9bc97a9f2", "colab_type": "text", "deepnote_cell_type": "markdown", "id": "yjQIGQ5fGpIU" }, "source": [ "* `maxiter`" ] }, { "cell_type": "code", "execution_count": 50, "metadata": { "cell_id": "00048-aeba2fbd-eff7-4719-951b-56237027be40", "colab": { "base_uri": "https://localhost:8080/", "height": 35 }, "colab_type": "code", "deepnote_cell_type": "code", "id": "xyyTsH7L7R-E", "outputId": "bc87c14b-b6a1-4c22-b68f-bef351841a12" }, "outputs": [ { "data": { "text/plain": [ "0.6666666666666671" ] }, "execution_count": 50, "metadata": {}, "output_type": "execute_result" } ], "source": [ "optimize.bisect(f,-10,10,args=(-3,2),xtol=1E-17,maxiter=10000)" ] }, { "cell_type": "code", "execution_count": 51, "metadata": { "cell_id": "00049-323945d2-f445-4c73-a03b-7cf551a293c7", "deepnote_cell_type": "code" }, "outputs": [ { "data": { "text/plain": [ "0.6666666666666671" ] }, "execution_count": 51, "metadata": {}, "output_type": "execute_result" } ], "source": [ "optimize.bisect(f,-10,10,args=(-3,2),xtol=1E-17,rtol=1E-15,maxiter=10000)" ] }, { "cell_type": "markdown", "metadata": { "cell_id": "00050-c2807d89-4849-4822-a64c-29a46d6a4a7d", "colab_type": "text", "deepnote_cell_type": "markdown", "id": "FiXdFcKEGuIU" }, "source": [ "* `full_output`" ] }, { "cell_type": "code", "execution_count": 60, "metadata": { "cell_id": "00051-fc299f20-d410-4146-abb2-8300637e6ebd", "colab": {}, "colab_type": "code", "deepnote_cell_type": "code", "id": "tQnA1hbu7vYE" }, "outputs": [], "source": [ "x0,r=optimize.bisect(f,-10,10,args=(-3,2),xtol=1E-17,full_output=True)" ] }, { "cell_type": "code", "execution_count": 62, "metadata": { "cell_id": "00052-3c1b9419-4a40-42bc-bc2f-a9a1b669af4e", "deepnote_cell_type": "code" }, "outputs": [ { "data": { "text/plain": [ "0.6666666666666671" ] }, "execution_count": 62, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x0" ] }, { "cell_type": "code", "execution_count": 61, "metadata": { "cell_id": "00053-9299c6e3-46eb-45f5-8cbf-8876ec77ecc1", "colab": { "base_uri": "https://localhost:8080/", "height": 35 }, "colab_type": "code", "deepnote_cell_type": "code", "id": "YIE2XkN68Bfg", "outputId": "4f964c0f-7ebb-488a-fe1a-084f6629647f" }, "outputs": [ { "data": { "text/plain": [ "scipy.optimize.zeros.RootResults" ] }, "execution_count": 61, "metadata": {}, "output_type": "execute_result" } ], "source": [ "type(r)" ] }, { "cell_type": "markdown", "metadata": { "cell_id": "00054-ea09397a-105a-44c7-8649-b1a8501965fd", "colab_type": "text", "deepnote_cell_type": "markdown", "id": "QmSG5wGKw7We" }, "source": [ "Check the code" ] }, { "cell_type": "code", "execution_count": 67, "metadata": { "cell_id": "00055-5293ed38-736c-4ebe-a108-315ccf0f95b3", "colab": { "base_uri": "https://localhost:8080/", "height": 35 }, "colab_type": "code", "deepnote_cell_type": "code", "id": "jBi0uZeE8NGk", "outputId": "aacf42d1-5bb2-4028-bbdc-e4ae465ff554" }, "outputs": [ { "data": { "text/plain": [ "0.6666666666666671" ] }, "execution_count": 67, "metadata": {}, "output_type": "execute_result" } ], "source": [ "r.root" ] }, { "cell_type": "markdown", "metadata": { "cell_id": "00056-c437f4a1-6b36-4ba1-bd74-250c0e273414", "deepnote_cell_type": "markdown" }, "source": [ "The help of the function is write directly in the code of the Scipy method `bisect`" ] }, { "cell_type": "code", "execution_count": 70, "metadata": { "cell_id": "00057-23e6e725-6f8f-44ba-b401-6ecdc26281f2", "colab": {}, "colab_type": "code", "deepnote_cell_type": "code", "id": "csnji5Tsw7Wf" }, "outputs": [], "source": [ "optimize.bisect??" ] }, { "cell_type": "markdown", "metadata": { "cell_id": "00058-31c575e2-f31a-49bc-afe7-75efb1cd75d8", "colab_type": "text", "deepnote_cell_type": "markdown", "id": "6_7lOL6Ww7Wj" }, "source": [ "Check the internal referenced code" ] }, { "cell_type": "code", "execution_count": 68, "metadata": { "cell_id": "00059-b095ae37-c407-4598-ab47-e48d54e37419", "colab": {}, "colab_type": "code", "deepnote_cell_type": "code", "id": "Qh9fqUvgw7Wk" }, "outputs": [ { "data": { "text/plain": [ "\u001b[0;31mType:\u001b[0m module\n", "\u001b[0;31mString form:\u001b[0m \n", "\u001b[0;31mFile:\u001b[0m /usr/local/lib/python3.7/dist-packages/scipy/optimize/_zeros.cpython-37m-x86_64-linux-gnu.so\n" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "optimize._zeros??" ] }, { "cell_type": "markdown", "metadata": { "cell_id": "00060-7cfe2241-3f28-42a9-bddf-0e04f4537090", "colab_type": "text", "deepnote_cell_type": "markdown", "id": "IJIfqraOG_fz", "tags": [] }, "source": [ "### Example\n", "Find the roots of\n", "$$x^3=2$$" ] }, { "cell_type": "code", "execution_count": 71, "metadata": { "cell_id": "00061-002ecac9-390d-4f4c-8656-292e9927b2a5", "colab": {}, "colab_type": "code", "deepnote_cell_type": "code", "id": "O1E6o0IRw7Wo" }, "outputs": [], "source": [ "def function(x):\n", " return x**3.0 - 2.0" ] }, { "cell_type": "code", "execution_count": 72, "metadata": { "cell_id": "00062-e91c31ee-715e-4f6d-a409-1cc8fdd80883", "deepnote_cell_type": "code" }, "outputs": [], "source": [ "todo=optimize.bisect(function,a,b,full_output=True)" ] }, { "cell_type": "code", "execution_count": 73, "metadata": { "cell_id": "00063-b14fbdeb-bf0d-4ad2-9e13-9d2607a28e8a", "deepnote_cell_type": "code" }, "outputs": [ { "data": { "text/plain": [ "tuple" ] }, "execution_count": 73, "metadata": {}, "output_type": "execute_result" } ], "source": [ "type(todo)" ] }, { "cell_type": "code", "execution_count": 74, "metadata": { "cell_id": "00064-f1e070b4-f3fd-46c8-b9f7-cb183c9037ef", "deepnote_cell_type": "code" }, "outputs": [ { "data": { "text/plain": [ "2" ] }, "execution_count": 74, "metadata": {}, "output_type": "execute_result" } ], "source": [ "len(todo)" ] }, { "cell_type": "code", "execution_count": 75, "metadata": { "cell_id": "00065-242e95c6-c8a9-4dee-80f1-19430fab181f", "deepnote_cell_type": "code" }, "outputs": [ { "data": { "text/plain": [ "1.2599210498938191" ] }, "execution_count": 75, "metadata": {}, "output_type": "execute_result" } ], "source": [ "result=todo[0]\n", "result" ] }, { "cell_type": "code", "execution_count": 76, "metadata": { "cell_id": "00066-415dd707-ea15-4ba2-bda9-9ab1b47031f1", "deepnote_cell_type": "code" }, "outputs": [ { "data": { "text/plain": [ " converged: True\n", " flag: 'converged'\n", " function_calls: 42\n", " iterations: 40\n", " root: 1.2599210498938191" ] }, "execution_count": 76, "metadata": {}, "output_type": "execute_result" } ], "source": [ "r=todo[1]\n", "r" ] }, { "cell_type": "code", "execution_count": 138, "metadata": { "cell_id": "00067-a527be75-2a66-434d-b67c-faa705eb776c", "deepnote_cell_type": "code" }, "outputs": [ { "data": { "text/plain": [ "scipy.optimize.zeros.RootResults" ] }, "execution_count": 138, "metadata": {}, "output_type": "execute_result" } ], "source": [ "type(r)" ] }, { "cell_type": "code", "execution_count": 131, "metadata": { "cell_id": "00068-89541b5a-6813-45e9-83a2-9f331abb7af9", "deepnote_cell_type": "code" }, "outputs": [ { "data": { "text/plain": [ "40" ] }, "execution_count": 131, "metadata": {}, "output_type": "execute_result" } ], "source": [ "r.iterations" ] }, { "cell_type": "code", "execution_count": 78, "metadata": { "cell_id": "00069-9ef2c972-7d32-4f42-9324-82b31f965bd5", "deepnote_cell_type": "code" }, "outputs": [], "source": [ "todo=optimize.bisect(function,a,b,full_output=False)" ] }, { "cell_type": "code", "execution_count": 79, "metadata": { "cell_id": "00070-53f7d183-2b5c-4a51-aadd-797bcd559dd3", "deepnote_cell_type": "code" }, "outputs": [ { "data": { "text/plain": [ "float" ] }, "execution_count": 79, "metadata": {}, "output_type": "execute_result" } ], "source": [ "type(todo)" ] }, { "cell_type": "code", "execution_count": 80, "metadata": { "cell_id": "00071-25c9a171-1542-459f-90bf-525f7917f52d", "colab": {}, "colab_type": "code", "deepnote_cell_type": "code", "id": "kABWCrNcw7Wq" }, "outputs": [], "source": [ "a = 0.0\n", "b = 2.0\n", "\n", "result,r=optimize.bisect(function,a,b,full_output=True)" ] }, { "cell_type": "code", "execution_count": 81, "metadata": { "cell_id": "00072-4d7e471e-09f6-4637-b351-24166e731076", "colab": { "base_uri": "https://localhost:8080/", "height": 35 }, "colab_type": "code", "deepnote_cell_type": "code", "id": "nlxlkNYyw7Wt", "outputId": "1ca5ef81-e255-4df2-fe73-ebff00268deb" }, "outputs": [ { "data": { "text/plain": [ "1.2599210498938191" ] }, "execution_count": 81, "metadata": {}, "output_type": "execute_result" } ], "source": [ "result" ] }, { "cell_type": "code", "execution_count": 82, "metadata": { "cell_id": "00073-e7c8de5c-de2c-4abe-813d-1379bc7f3dd5", "colab": { "base_uri": "https://localhost:8080/", "height": 35 }, "colab_type": "code", "deepnote_cell_type": "code", "id": "L3HyS2OGw7Wv", "outputId": "d0094839-53ce-43cc-ab3b-c276f0275ff3" }, "outputs": [ { "data": { "text/plain": [ "1.2599210498948732" ] }, "execution_count": 82, "metadata": {}, "output_type": "execute_result" } ], "source": [ "2**(1/3.0)" ] }, { "cell_type": "markdown", "metadata": { "cell_id": "00074-de669b29-3f1d-47ba-b3e2-91eeb44f1a96", "colab_type": "text", "deepnote_cell_type": "markdown", "id": "TKlBRLXUw7Wy" }, "source": [ "Check `r` object by using the `` key next" ] }, { "cell_type": "code", "execution_count": 122, "metadata": { "cell_id": "00075-5ba7bdc2-e19e-4a3b-9ea6-6bcfa0d668e3", "colab": { "base_uri": "https://localhost:8080/", "height": 35 }, "colab_type": "code", "deepnote_cell_type": "code", "id": "TkfTsXumw7Wz", "outputId": "d618fc43-05eb-4e22-de7e-4db81ad77435" }, "outputs": [ { "data": { "text/plain": [ "40" ] }, "execution_count": 122, "metadata": {}, "output_type": "execute_result" } ], "source": [ "r.iterations" ] }, { "cell_type": "code", "execution_count": 123, "metadata": { "cell_id": "00076-9e9de48b-551d-47e0-bc7b-952db24db485", "deepnote_cell_type": "code" }, "outputs": [ { "data": { "text/plain": [ "42" ] }, "execution_count": 123, "metadata": {}, "output_type": "execute_result" } ], "source": [ "r.function_calls" ] }, { "cell_type": "code", "execution_count": 74, "metadata": { "cell_id": "00077-2d6988f0-e37b-4d62-b418-1b2964a9d409", "colab": { "base_uri": "https://localhost:8080/", "height": 35 }, "colab_type": "code", "deepnote_cell_type": "code", "id": "HpBNHlWDw7W3", "outputId": "25eb56bd-9ec6-4c4f-9a08-b3035ac7260c" }, "outputs": [ { "data": { "text/plain": [ "9.094947017729282e-13" ] }, "execution_count": 74, "metadata": {}, "output_type": "execute_result" } ], "source": [ "2**(-r.iterations)" ] }, { "cell_type": "code", "execution_count": 75, "metadata": { "cell_id": "00078-1beba11e-b7dc-43a4-a5c2-81c77a882086", "colab": { "base_uri": "https://localhost:8080/", "height": 35 }, "colab_type": "code", "deepnote_cell_type": "code", "id": "tvkJ0Jy5w7W6", "outputId": "c9d1ebf4-0f00-4992-c901-562c1ae6dfde" }, "outputs": [ { "data": { "text/plain": [ "1.0540457395791236e-12" ] }, "execution_count": 75, "metadata": {}, "output_type": "execute_result" } ], "source": [ "abs((2**(1/3.0)-result))" ] }, { "cell_type": "code", "execution_count": 76, "metadata": { "cell_id": "00079-d30fe040-6a03-44e0-8f2f-504565fb9a97", "colab": { "base_uri": "https://localhost:8080/", "height": 35 }, "colab_type": "code", "deepnote_cell_type": "code", "id": "4m6VCcpiw7W8", "outputId": "8229c9e2-8cb1-4a0c-a2e1-51f77ae36dc4" }, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 76, "metadata": {}, "output_type": "execute_result" } ], "source": [ "r.converged" ] }, { "cell_type": "code", "execution_count": 77, "metadata": { "cell_id": "00080-0cd89992-2a5e-4f13-900f-69d39f2c8d91", "colab": {}, "colab_type": "code", "deepnote_cell_type": "code", "id": "7bp9qun4w7W_" }, "outputs": [], "source": [ "result,r=optimize.bisect(function,a,b,xtol=1E-17,full_output=True)" ] }, { "cell_type": "code", "execution_count": 78, "metadata": { "cell_id": "00081-8788695c-7a8e-4f56-b96a-3606127a2750", "colab": { "base_uri": "https://localhost:8080/", "height": 35 }, "colab_type": "code", "deepnote_cell_type": "code", "id": "kmUCIBwlw7XB", "outputId": "70a1d062-9ec4-4082-afbc-4056962e0d46" }, "outputs": [ { "data": { "text/plain": [ "51" ] }, "execution_count": 78, "metadata": {}, "output_type": "execute_result" } ], "source": [ "r.iterations" ] }, { "cell_type": "code", "execution_count": 79, "metadata": { "cell_id": "00082-5dd094e4-636b-4826-b41e-f9f0abb5a298", "colab": { "base_uri": "https://localhost:8080/", "height": 35 }, "colab_type": "code", "deepnote_cell_type": "code", "id": "zQTyAChgw7XE", "outputId": "c1a5a72d-53f8-4f53-ef5c-26ef8fcd2305" }, "outputs": [ { "data": { "text/plain": [ "2.220446049250313e-16" ] }, "execution_count": 79, "metadata": {}, "output_type": "execute_result" } ], "source": [ "abs((2**(1/3.0)-result))" ] }, { "cell_type": "code", "execution_count": 80, "metadata": { "cell_id": "00083-059b3bf7-9af3-4004-8cfd-979ed9009fad", "colab": { "base_uri": "https://localhost:8080/", "height": 35 }, "colab_type": "code", "deepnote_cell_type": "code", "id": "eXvvDw-4EZr-", "outputId": "bb036e6a-67cb-4b17-b1db-8f1fe1e38021" }, "outputs": [ { "data": { "text/plain": [ "4.440892098500626e-16" ] }, "execution_count": 80, "metadata": {}, "output_type": "execute_result" } ], "source": [ "2**(-r.iterations)" ] }, { "cell_type": "markdown", "metadata": { "cell_id": "00084-2c87b9b5-238e-4d1e-833f-9fd2332a8f6c", "colab_type": "text", "deepnote_cell_type": "markdown", "id": "Mf8r3RfJw7XG" }, "source": [ "__ ACTIVITY __" ] }, { "cell_type": "markdown", "metadata": { "cell_id": "00085-e99f02c8-9220-4cc8-a9a4-42fa0b50a329", "colab_type": "text", "deepnote_cell_type": "markdown", "id": "gDKsdwSZw7XH" }, "source": [ "\n", "In an IPython notebook, use the scipy implementation and find the first solution to the equation\n", " \n", "$ 7 = \\sqrt{x^2+1}+e^x\\sin x $\n", " \n", "CLUE: Check graphically (with matplotlib) that the solution is within the interval $[0,2]$.\n", "" ] }, { "cell_type": "markdown", "metadata": { "cell_id": "00086-92f77970-8631-4dca-b7ab-ae74f1c42d56", "colab_type": "text", "deepnote_cell_type": "markdown", "id": "3Sr162oVw7XI" }, "source": [ "[solution](./solutions/sltn_bisect.ipynb)" ] }, { "cell_type": "code", "execution_count": 38, "metadata": { "cell_id": "00087-964a4748-9da4-4939-be15-755aa1713103", "colab": {}, "colab_type": "code", "deepnote_cell_type": "code", "id": "e5at_XXzHSlk" }, "outputs": [], "source": [ "a,b,c=(2,3,4)" ] }, { "cell_type": "code", "execution_count": 39, "metadata": { "cell_id": "00088-6ace9e02-7b44-451a-bc75-cfc036b3ca60", "colab": { "base_uri": "https://localhost:8080/", "height": 35 }, "colab_type": "code", "deepnote_cell_type": "code", "id": "geD3ecfDHXFE", "outputId": "e823e103-fbd6-461d-c9ed-82c11a9c5d28" }, "outputs": [ { "data": { "text/plain": [ "4" ] }, "execution_count": 39, "metadata": {}, "output_type": "execute_result" } ], "source": [ "c" ] }, { "cell_type": "code", "execution_count": 40, "metadata": { "cell_id": "00089-2a25631d-85d7-4206-958d-3d532c3a40d7", "colab": {}, "colab_type": "code", "deepnote_cell_type": "code", "id": "zrY4Rc19E_zC" }, "outputs": [], "source": [ "import numpy as np" ] }, { "cell_type": "code", "execution_count": 41, "metadata": { "cell_id": "00090-fd1017c0-95b8-4254-b587-6d468a2c4785", "colab": { "base_uri": "https://localhost:8080/", "height": 35 }, "colab_type": "code", "deepnote_cell_type": "code", "id": "SZhTQgvLFDmg", "outputId": "51f03fe3-2248-4cd6-832e-0f4ecebc3059" }, "outputs": [ { "data": { "text/plain": [ "2.718281828459045" ] }, "execution_count": 41, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Solo usar para constante\n", "np.e" ] }, { "cell_type": "code", "execution_count": 42, "metadata": { "cell_id": "00091-3c32f057-2093-4e45-ae31-3efe926363c4", "colab": { "base_uri": "https://localhost:8080/", "height": 35 }, "colab_type": "code", "deepnote_cell_type": "code", "id": "ela6ZuPcFGPQ", "outputId": "9237d95b-cf27-43b2-aab8-242954b8ec85" }, "outputs": [ { "data": { "text/plain": [ "2.718281828459045" ] }, "execution_count": 42, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.exp(1)" ] }, { "cell_type": "markdown", "metadata": { "cell_id": "00092-16a277e5-2710-41fa-8654-b258dc4899c3", "colab_type": "text", "deepnote_cell_type": "markdown", "id": "hBSyBUNrw7XJ", "jp-MarkdownHeadingCollapsed": true, "tags": [] }, "source": [ "### Example 2" ] }, { "cell_type": "markdown", "metadata": { "cell_id": "00093-5f93b55d-57e1-47c9-8ef7-ca48489465a2", "colab_type": "text", "deepnote_cell_type": "markdown", "id": "qpIK0RxFw7XK" }, "source": [ "In orbital mechanics, when solving the central-force problem it becomes necessary to solve the Kepler's equation. This is a transcendental equation that relates the orbital parameters of the trajectory.\n", "\n", "*Kepler equation:* $M = E - \\epsilon \\sin E$\n", "\n", "where $M$ is the mean anomaly, $E$ the eccentric anomaly and $\\epsilon$ the eccentricity. The mean anomaly can be computed with the expression\n", "\n", "$$M = n\\ t = \\sqrt{ \\frac{GM}{a^3} } t$$\n", "\n", "where $n$ is the mean motion, $G$ the gravitational constant, $M$ the mass of the central body and $a$ the semi-major axis. $t$ is the time where the position in the trajectory will be computed.\n", "\n", "The coordinates $x$ and $y$ as time functions can be recovered by means of the next expressions\n", "\n", "$$x(t) = a(\\cos E - \\epsilon)$$\n", "\n", "$$y(t) = b\\sin E$$\n", "\n", "where $b = a \\sqrt{1-\\epsilon^2}$ is the semi-minor axis of the orbit and the implicit time-dependence of the eccentric anomaly $E$ is computed through the Kepler's equation." ] }, { "cell_type": "markdown", "metadata": { "cell_id": "00094-d03c5ad9-4375-4392-813c-5324ba26d2cd", "colab_type": "text", "deepnote_cell_type": "markdown", "id": "FlKupbitw7XL" }, "source": [ "**Problem:**\n", "\n", "For a stallite orbiting the earth in a equatorial trajectory with eccentricity $\\epsilon = 0.5$ at a geostationary distance for the semi-major axis, tabulate the positions $x$ and $y$ within the orbital plane in intervals of $15$ min during $5$ hours.\n", "\n", "**Parameters:**\n", "\n", "- $\\epsilon = 0.5$\n", "\n", "- $a = 35900$ km\n", "\n", "- $G = 6.67384 \\times 10^{-11}$ m$^3$ kg$^{-1}$ s$^{-2}$\n", "\n", "- $M_{\\oplus} = 5.972\\times 10^{24}$ kg" ] }, { "cell_type": "code", "execution_count": 43, "metadata": { "cell_id": "00095-42d0cae2-815b-4a70-a287-262d37a0119a", "colab": {}, "colab_type": "code", "deepnote_cell_type": "code", "id": "URwGARM5w7XL" }, "outputs": [], "source": [ "#====================================================================\n", "#Parameters\n", "#====================================================================\n", "#Eccentricity\n", "eps = 0.5\n", "#Semi-major axis [m]\n", "a = 35900e3\n", "#Gravitational constant [m3kg-1s-2]\n", "GC = 6.67384e-11\n", "#Earth mass [kg]\n", "Me = 5.972e24\n", "\n", "#Semi-minor axis [m]\n", "b = a*(1-eps**2.0)**0.5\n", "#Mean motion\n", "n = ( GC*Me/a**3.0 )**0.5\n", "\n", "#Hour to Second\n", "HR2SC = 3600.\n", "#Initial time [hr]\n", "t0 = 0*HR2SC\n", "#Final time [hr]\n", "tf = 5*HR2SC\n", "#Time step [hr]\n", "tstep = 0.25*HR2SC\n", "#Number of maxim iterations\n", "Niter = 56\n", "#Root interval\n", "a0 = -10\n", "b0 = 10\n", "\n", "#====================================================================\n", "#Kepler Function\n", "#====================================================================\n", "def kepler( E ):\n", " func = E - eps*np.sin(E) - n*t\n", " return func\n", "\n", "#====================================================================\n", "#Position function\n", "#====================================================================\n", "def r(E):\n", " x = a*(np.cos(E)-eps)\n", " y = b*np.sin(E)\n", " return [x/1.e3, y/1.e3]" ] }, { "cell_type": "code", "execution_count": 45, "metadata": { "cell_id": "00096-922ad7be-f5db-47f8-ab75-c3ff9d860122", "colab": {}, "colab_type": "code", "deepnote_cell_type": "code", "id": "z4PZbsmIw7XQ", "outputId": "bb53b14d-8d19-459f-ecf6-4e9d8988266d" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "In 0.000000 hours, the satellite is located at (17950.000000,0.000000) km\n", "In 0.250000 hours, the satellite is located at (17454.741542,5146.426647) km\n", "In 0.500000 hours, the satellite is located at (16033.097675,10023.437750) km\n", "In 0.750000 hours, the satellite is located at (13848.847528,14430.262364) km\n", "In 1.000000 hours, the satellite is located at (11104.379909,18261.701894) km\n", "In 1.250000 hours, the satellite is located at (7989.437466,21493.410338) km\n", "In 1.500000 hours, the satellite is located at (4657.168469,24151.489610) km\n", "In 1.750000 hours, the satellite is located at (1221.066166,26286.121177) km\n", "In 2.000000 hours, the satellite is located at (-2238.747501,27954.872718) km\n", "In 2.250000 hours, the satellite is located at (-5667.264143,29214.008624) km\n", "In 2.500000 hours, the satellite is located at (-9027.373695,30114.739827) km\n", "In 2.750000 hours, the satellite is located at (-12294.392344,30702.085866) km\n", "In 3.000000 hours, the satellite is located at (-15452.164136,31014.965936) km\n", "In 3.250000 hours, the satellite is located at (-18490.361033,31086.789919) km\n", "In 3.500000 hours, the satellite is located at (-21402.633716,30946.195086) km\n", "In 3.750000 hours, the satellite is located at (-24185.347784,30617.769816) km\n", "In 4.000000 hours, the satellite is located at (-26836.717825,30122.702148) km\n", "In 4.250000 hours, the satellite is located at (-29356.211228,29479.336137) km\n", "In 4.500000 hours, the satellite is located at (-31744.135350,28703.638662) km\n", "In 4.750000 hours, the satellite is located at (-34001.350021,27809.586870) km\n" ] }, { "data": { "text/plain": [ "Text(0, 0.5, '$y$ coordinate [km]')" ] }, "execution_count": 45, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaIAAAEKCAYAAABQRFHsAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAAu0ElEQVR4nO3deZhU1bX38e+iGQQRmaRVZiLXiMYBUDFO4IBoYjBGE7UjOLZRjPpqNCr3OtzITYxRo9dERUUxtOI8xOCACg65QcEZRKVFRhUHUEQEBdb7x94diqaa7oKqOtVVv8/znKdO7TPU2jTVq88+++xt7o6IiEhSmiQdgIiIlDYlIhERSZQSkYiIJEqJSEREEqVEJCIiiVIiEhGRRDVNOoCkdezY0Xv06JF0GOv5+uuv2XzzzZMOIy9U1+JTKvWE0qlr7Xq+8sorn7n7Vtk4d8knoh49ejBt2rSkw1jP5MmTGThwYNJh5IXqWnxKpZ5QOnWtXU8zm5utc6tpTkREEqVEJCIiiVIiEhGRRCkRiYhIopSIpLBVVUGPHtCkSXitqko6IhHJMiUiyZ1NTSJVVVBZCXPngnt4razcuGSkhCZSsEq++7bkSE0SWb48vK9JIgAVFSGxrFgBX31Fy4UL4bXXYOlS+OqrsCxdChdeuPb4GsuXw5lnwkcfQYsWYWnefN3X2mUTJ8Jll4XPSxeLiCRKiUjSq6qCkSNh3jzo1g1Gjar/l/bKlWH/Dz6As85Kn0SGD4df/zokm1WrANgz09i++ALOPz/To9aP5dRT4YUXoEsX6Nx53dc2bdIftzH/LiKyQUpEsr66rmbWrIH99guJZs6c8Jq6fPhhuNLZkNWrYcmSsN68ObRpwzfNmtGyUyfYYouwtGkTXu+5JySs2tq0gZNPhm+/DcmvvtfXX08fyzffwM03p9/WuvX6CerDD+Huu8M5U/9dQMlIZBMoEcn6Lroo/dXMsGEbPq5Jk3CV0LMnTJ0KX3+9/j7bbgtvvhkSTfPmALxU15PpAweumxABWrWCv/41s1/8PXqEpFHbVlvBJZfAwoWwYMHa1wULYNkyeOedsGzI8uUwYgS0bw99+0J5ecPjEhFAiUhWrgyJYdo0eOWV8Dp/ft37b7NN+MXes+f6S5cu0KxZ2K/2VRWEJPLHP0KHDg2LrSbZbGpT2KhR6WO59tr053IPzX+piWnhwnCfKZ0vv4TDDgvr224bElLq0qULmGUWs0gpcfe8LMBmwMvAG8AM4PJY3hN4CagG7gGax/IW8X113N4j5VwXxfJ3gUNSyofEsmrgwobE1a9fPy9EkyZN2vSTjBvn3r27u1l4HTvW/dVX3UePdq+sdO/b171ZM/fwq7f+pWvXTfv8cePS7paVumYplg3q3j39v8sWW7jvt194Tbe9Y0f3wYPdL7zQp196qXt1tfuaNdmJqUDl5WdaIEqlrrXrCUzzbOWHbJ2o3g8CA1rH9WYxuQwA7gWOieU3AafH9TOAm+L6McA9cb1PTGYtYhJ7HyiLy/tAL6B53KdPfXEVbSIaN869Zcv6k4uZ+w47uB9/vPt117m/+KL7bbe5t2q17n6tWuXsF2Wj+SKPG7fhf5fVq93fe899/Hj3Cy5wP+gg9/bt0/+7t2zp3qRJ3v6N863R/EyzoFTqmstElLemuRj4svi2WVwcOAA4LpaPBS4DbgSGxnWA+4EbzMxi+Xh3Xwl8YGbVwB5xv2p3nw1gZuPjvm/nrlYFZuVKePFFePJJ+POf4bvv1t+naVM4+mjo3z8su+0W7tek2nvv0PVZvcPWVV9TYZMm0Lt3WH7xi1DmHvZ99VV49VU+nziRDnPmwKJF65+/5n7TjjvCLruoOU9KhoX8kKcPMysDXgG2A/4CXAVMcfft4vauwOPuvpOZTQeGuPuCuO19Qk/fy+Ix42L5bcDj8SOGuPspsfx4YE93PzNNHJVAJUB5eXm/8ePH56jGG2/ZsmW0bt16wzu503L+fNpPnUr7qVNp+8YblNU8K1PXIWY89+yzWYx00zWorkWipq77H3AAtoHv3sqOHfl8zz1ZvOeeLOnXj9WtWuUxyk1Xij/TYle7noMGDXrF3ftn5eTZurTKZAHaApOAfQhXMTXlXYHpcX060CVl2/tAR+AG4Jcp5bcBR8Xl1pTy44Eb6oul4Jvmat9HuOUW9wceCPd40t2z2Hln9/PPd+/UKX2TUPfuidWpLqXStOGeUte67je1bu2+7bbrljVr5n7gge5XX+3+zjvh/lKBK8mfaZHLZdNcIkP8uPsXhES0F9DWzGqaCLsAC+P6QkJiIm7fEvg8tbzWMXWVN17phrg59VT42c9g9OjwvmNHOPZYuOOO8JzLG2+EnmnXXBN6hqVq1So0JUnyRo1K//O56abQS++118I+P/xhePbqmWfgvPPg+9+H7bYLDww/8cTa0SJAwxhJ45WtjFbfAmwFtI3rLYEXgB8D97FuZ4Uz4voI1u2scG9c35F1OyvMJnRUaBrXe7K2s8KO9cVVqFdE/6qqct9yy/R/Nbdo4T5qlPu0aeEGeV0aSa+sUvmL0r1WXRv68/nsM/eqKveKCvcOHdb9v9CypfuPf+x+4onrd05JsPNDyf5Mi1hRdFYAtgHGxvtETWJieczM3gbGm9kVwGuEpjbi699iZ4TFhGSEu88ws3sJnRBWASPcfTWAmZ0JPElITGPcfUb+qpcFixbBvffCXXcxYMqUuvf79lu4+OL6z1dRoQ4GhayhP58OHeC448KyejW8/DL84x8wYUK4cnrssfTHLV8eOlbo/4AUuHz2mnsT2C1N+WzW9npLLV8BHF3HuUYB67UxufsEYMImB5tPS5fCww/DXXfB00+HXzTA6s02o6ysLP3oBN265TdGKRxlZbDXXmG54orQHPv443DKKen3nzcvjOnXVM+uS+HSNBD5ktp+3707nHsu/PznYUiY4cNDl2szOPxwuPtu/vngg2EcNN3nkQ3Zdtsw7l737um3u4c/XP7rv8L4gCIFSIkoH2p3Opg3Lwwvc9994WbzfvuFm9QffwyPPgrHHMOali1Dk8ro0eGXjFl4HT1aTS2yvnSdH5o1g623DlNmXHEF9OoFhx4KDz2U/hkzkYQoEeXaihVw9tnrDyIK0LZtSErPPQennZZ+DLaKivCX7Jo14VVJSNJJ90fL7beHprvnngvbmzcPPe2OPDJcJY0cGUZNF0mYElGufPhhaA7p1g0+/zz9Pl9+CV27pt8mkql0f7SYhSvucePCwK3XXgs77BCuvv/nf8JV0iGHwAMP6CpJEqNElE3u8K9/hed6uncPzSGffrp2ROra1OlA8qlDBzjnHJgxI0wIePzxYSinp56Co44KfxRddBG8/37YX88lSZ4oEW2M2l/QsWPhb3+DPfYIDyCOHx+S0tFHhy/87ber04EUDjPYZx+4885w5X7dddCnT3h84A9/CA/M7rQTnHTSug9TV1YqGUlOKBFlKt1oByecECaNmzYt/NV50UWh7f3ee8MXXp0OpFC1bx9GaZg+PQyYO2wYbLZZuGr69tt19615Lkkky/RwQaZGjkzf8aBZM7jxxvDQYcuW62/Xw6VSyMzCqOt77x1Gbm/fPv1+8+blNSwpDboiylRdX8RVq8LzHOmSkEhj0q5d3c8lmcH114cpR0SyRIkoU3V1MFDHAykm6Z5LatIk9Mg7+2z4j/+AMWPCH2Aim0iJKFN1jZqsjgdSTNLd17zzTnjkkdCRYd680AKw447hXuiaNUlHLI2YElGm1PFASkW655J+8hN4/fXQaed734P33guz0fbrFwZh9fxNtCnFQ4loY2i0AyllZWWhU87MmWE8xM6dQ3L60Y9g333h+eeTjlAaGSUiEdk4zZqFRxlmzYKrrw6PLvzzn7D//ux8wQXwyitJRyiNhBKRiGyali3DaPKzZ8Pll8MWW9B+6lTo3z+M2PD220lHKAVOiUhEsqNNG7jkEvjgA+b94hfhwdgHHoAf/CA89P3nP2vIIElLiUhEsqtDB2b/6ldhzLrTTw+JZ+xY+H//T0MGSVpKRCKSG9tuC3/9K7z7Lmy++frbNWSQREpEIpJbvXqlHxYLNGSQAEpEIpIPGxp5ZPRoPX9U4pSIRCT30o1IUlYWEtBpp8GQITB/fjKxSeKUiEQk99KNSDJ2bJi7q337MDnfTjvBHXfo6qgEKRGJSH6kG5HkF78Icx8NHQpLl8KJJ4ZhhD78MOloJY/ylojMrKuZTTKzt81shpmdHcsvM7OFZvZ6XA5LOeYiM6s2s3fN7JCU8iGxrNrMLkwp72lmL8Xye8yseb7qJyIbaeut4aGHwizHbdvCY4+Fq6OqKl0dlYh8XhGtAs5z9z7AAGCEmfWJ2651913jMgEgbjsG2BEYAvzVzMrMrAz4C3Ao0Ac4NuU8V8ZzbQcsAU7OV+VEZBOYwS9/GWaKPfRQWLIkvP/Zz+CTT5KOTnIsb4nI3T9y91fj+lfATKDzBg4ZCox395Xu/gFQDewRl2p3n+3u3wLjgaFmZsABwP3x+LHAETmpjIjkRufO8I9/wK23whZbhCulHXeE++5LOjLJIfMELn3NrAfwPLATcC5wArAUmEa4alpiZjcAU9x9XDzmNuDxeIoh7n5KLD8e2BO4LO6/XSzvCjzu7jul+fxKoBKgvLy83/jx43NT0U2wbNkyWrdunXQYeaG6Fp9s1LPFokVsf9VVtI+Dp34yaBCzzj6b77bcMhshZk2p/kwHDRr0irv3z8rJ3T2vC9AaeAU4Mr4vB8oIV2ejgDGx/AbglynH3QYcFZdbU8qPj/t2JFwp1ZR3BabXF0+/fv28EE2aNCnpEPJGdS0+WavnmjXuN97ovvnm7uBeXu7+8MPZOXeWlOrPFJjmWcoLee01Z2bNgAeAKnd/EMDdF7n7andfA9xCaHoDWBiTSY0usayu8s+BtmbWtFa5iDRWZvCrX8Gbb8L++8OiRXDEETBsWLiPJEUhn73mjHBVM9Pdr0kp3yZlt58C0+P6o8AxZtbCzHoCvYGXgalA79hDrjmhQ8OjMUNPIlwxAQwHHsllnUQkT3r1gmefheuuC9NO/O1voaxTJ43mXQSa1r9L1uxNaEZ7y8xej2UXE3q97Qo4MAc4DcDdZ5jZvcDbhB53I9x9NYCZnQk8SWjSG+PuM+L5fguMN7MrgNcIiU9EikGTJnDWWWEUhh/9CKqr126rGc0bNGNyI5S3ROTuLwKWZtOEDRwzinDfqHb5hHTHufts1jbtiUgx+o//gG+/Xb+8ZjRvJaJGRyMriEjjU9e4dBrNu1FSIhKRxmdDo3lPnJi/OCQrlIhEpPHZ0Gjehx0GY8YkE5dsFCUiEWl80o3mfccdcMEFsGoVnHwy/Od/aqy6RiKfveZERLKnoiJ9x4RevWDEiHDVNHs23H47tGiR//ikwXRFJCLF5bTT4O9/h9at4e674eCD4fPPk45KNkCJSESKz6GHwosvhkFUX3gBfvhDeP/9pKOSOigRiUhx2mUXmDIFdt4Z3nsPBgyAf/0r6agkDSUiESleXbqEK6IhQ+Czz2DQILj//vqPk7xSIhKR4tamTbhndNppsHIlHH00XHWVetQVECUiESl+TZvCjTfClVeG9xdcAGecEbp6S+Lq7b5tZu0bcJ417v7FpocjIpIjZiEB9egRppG46aYwWOo994TZYCUxDXmO6MO4pBuwtEYZsIExN0RECsTPfx7uHf3kJ/D447DvvmF68s6dk46sZDWkaW6mu/dy9551LYRJ6UREGocf/jD0qOvdG954A/bcE37/+3C1pPmN8q4hiWivLO0jIlI4ttsudOfeZx9YuBAuvjg01bmvnd9IySgv6k1E7r4iG/uIiBScDh3CaN21B1CFtfMbSc41eKw5M+sPjAS6x+MMcHffOUexiYjk3mabwTffpN+m+Y3yIpNBT6uA84G3gDW5CUdEJAHduoXmuHTlknOZPEf0qbs/6u4fuPvcmiVnkYmI5Eu6+Y0ATj89/7GUoEwS0aVmdquZHWtmR9YsOYtMRCRfas9v1LJlKL/hhtCRQXIqk0R0IrArMAQ4PC4/zkFMIiL5V1EBc+bAmjVh2oi994YFC8JI3l9+mXR0RS2Te0S7u/v2OYtERKRQtGwJjzwSktFbb8FPfxoeftUEezmRyRXR/5lZn5xFIiJSSDp0gCeegK23hkmT4MQTw9WSZF0miWgA8LqZvWtmb5rZW2b2ZkMPNrOuZjbJzN42sxlmdnYsb29mE81sVnxtF8vNzK43s+r4eX1TzjU87j/LzIanlPeLcVXHYzc0LJGIyIb16AETJqyd7fXCC5OOqChlkoiGAL2Bway9P3R4BsevAs5z9z6EpDYiXmFdCDzj7r2BZ+J7gEPj5/UGKoEb4d+DsF4K7AnsQehE0S4ecyNwaspxQzKIT0RkfbvtBg8+GEbwvuoq+N//TTqiopNJIuqY2m07dt3+QUMPdveP3P3VuP4VMBPoDAwFxsbdxgJHxPWhwJ0eTAHamtk2wCHARHdf7O5LgInAkLitjbtPcXcH7kw5l4jIxjv4YLjttrB+9tnwwAPJxlNkMumscIuZDXP36QBmdixwDvBYph9qZj2A3YCXgHJ3/yhu+hgoj+udgfkphy2IZRsqX5CmPN3nVxKusigvL2fy5MmZViHnli1bVpBx5YLqWnyKsp7dutHtlFPodeutrDn2WN7405/4cuedi7OuaeSynpkkoqOA+83sOGBfYBihmS4jZtYaeAA4x92Xpt7GcXc3s5xPm+juo4HRAP379/eBAwfm+iMzNnnyZAoxrlxQXYtP0dZz//2hWTOa3Hgju112Gfzzn0yG4qxrLbn8mTa4ac7dZwPHAA8CPwMGu3tGnevNrBkhCVW5+4OxeFFsViO+fhLLFwJdUw7vEss2VN4lTbmISHaYhXtEQ4fCkiUwZAjNP/ss6agavXoTUU3vuNhD7n6gPdATeCnDXnMG3EaY3+ialE2PAjU934YDj6SUD4u95wYAX8YmvCeBwWbWLnZSGAw8GbctNbMB8bOGpZxLRCQ7ysrgrrtgwACYN4+df/tbWLo06agatYY0zWVr9IS9geOBt8zs9Vh2MfAH4F4zOxmYC/w8bpsAHAZUA8sJIzvg7ovN7HfA1Ljff7v74rh+BnAH0BJ4PC4iItnVqhX8/e+w9960fu89OPLI0M27efOkI2uUGpKIHnL3vhvawcxerW8fd3+RuqcbPzDN/g6MqONcY4AxacqnATttKA4Rkazo2BGeeIJv+/Wj+TPPwEknwZ13hhleJSMNSUQ71NMEZ8CWWYpHRKTx6NmTN3//e/qfd16YzbVLF/jDH5KOqtFpSCL6fgP2Wb2pgYiINEbLtt8e7r8fDj8crrwSunaFEWkbc6QODZkqfG4DlgX1nUdEpGgNGQK33BLWzzwTOnUKTXQ9eoQrJdkgNWaKiGTDCSfAUUeF9U8/Bfcw62tlpZJRPZSIRESy5eWX1y9bvhxGjsx/LI2IEpGISLbMn5++fN68/MbRyDQ4EcUHS39pZpfE993MbI/chSYi0sh065ZZuQCZXRH9FdgLODa+/wr4S9YjEhFprEaNCg+71vab3+Q/lkYkk0S0p7uPAFYAxCkY9BixiEiNigoYPRq6dw/j0m22WSh/6qnQeUHSyiQRfWdmZYADmNlWgObNFRFJVVEBc+aEacVnzYI2bcJwQHfdlXRkBSuTRHQ98BDQycxGAS8Cv89JVCIixaBLF7gmjvF81lnw8cfJxlOgMpkGogq4gJB8PgKOcPd7cxWYiEhROOkkGDwYFi8OIy6oiW49mfSau9Ld33H3v7j7De4+08yuzGVwIiKNnlkYdaF1a3jwQbjvvqQjKjiZNM0dnKbs0GwFIiJStLp1g6uuCusjRoSRF+TfGjIx3ulm9hawfc0EeXH5AHgr9yGKiBSBykoYNAg++wx+/eukoykoDbkiugs4nDBj6uEpSz93r8hhbCIixaNJE7j11vCc0T33wEMPJR1RwWjI6Ntfuvscdz8WWAqUA92Bncxsv1wHKCJSNHr1Wjtf0emnw+efJxtPgciks8IpwPPAk8Dl8fWy3IQlIlKkRoyAffeFRYvgnHOSjqYgZNJZ4Wxgd2Cuuw8CdgO+yEVQIiJFq0kTuO22MOrCuHHw2GNJR5S4TBLRCndfAWBmLdz9HWD73IQlIlLEevcO49IBnHYafPFFouEkLZNEtMDM2gIPAxPN7BFgbi6CEhEpemefDXvtBR9+COeem3Q0icpkZIWfuvsX7n4Z8F/AbcDQXAUmIlLUyspgzBho0QJuvx2eeCLpiBKTSWeFFmZ2nJldDOwP7ApclMHxY8zsEzObnlJ2mZktNLPX43JYyraLzKzazN41s0NSyofEsmozuzClvKeZvRTL7zEzjQwuIoXt+9+Hyy8P66eeCkuXJhtPQjJpmnuEcAW0Cvg6ZWmoO4Ahacqvdfdd4zIBwMz6AMcAO8Zj/mpmZXH0778QRnToAxwb9wW4Mp5rO2AJcHIGsYmIJOO882D33WHBAjj//KSjSUTTDPbt4u7pEkmDuPvzZtajgbsPBca7+0rgAzOrBmpmg61299kAZjYeGGpmM4EDgOPiPmMJXctv3Nh4RUTyomnT0ETXt2+Yy+jnP4cDD0w6qrzK5Iro/8zsBzmI4cw4ZNAYM2sXyzoDqZO/L4hldZV3AL5w91W1ykVECt9OO8Ell4T1U06BZcuSjSfPMrki2gc4IY4xtxIwwN195034/BuB3xEm2/sdcDVw0iacr0HMrBKoBCgvL2fy5Mm5/siMLVu2rCDjygXVtfiUSj0he3W1AQPo27s3W8yaxapOnShbsYKVnTox+5RT+OSggzY90E2Uy59pJoko6yNtu/uimnUzuwWoebJrIdA1ZdcusYw6yj8H2ppZ03hVlLp/us8dDYwG6N+/vw8cOHDTKpIDkydPphDjygXVtfiUSj0hy3U98US4+GKafvMNAJstWkSfa6+lzw47hJlfE5TLn2km3bfnpls25cPNbJuUtz8FanrUPQocE3vq9QR6Ay8DU4HesYdcc0KHhkfd3YFJwFHx+OGEzhUiIo3HzTevX7Z8OYwcmf9Y8qjeKyIze9Hd9zGzrwhNaJb66u5tGvJBZnY3MBDoaGYLgEuBgWa2azzfHOA0wklnmNm9wNuEXnoj3H11PM+ZhHHuyoAx7j4jfsRvgfFmdgXwGuE5JxGRxmPevMzKi0S9icjd94mvW2zKB8XRu2urM1m4+yhgVJryCcCENOWzWduzTkSk8enWDeamaWjq1i3/seRRQ66INjj2hLtfk71wRERK2KhRYQK95cvXljVvvnZcuiLVkM4KNVdC2xNG3340vj+ccN9GRESyoaZDwsiRa6+MOneG446r+5gi0JCJ8S5398sJPdH6uvt57n4e0A8o7utFEZF8q6iAOXPCVdFWW8EHH8BzzyUdVU5l8kBrOfBtyvtvY5mIiGRby5Zw5plh/aqrko0lxzJJRHcCL8eBSi8DXiKMHyciIrlwxhkhIU2YADNm1L9/I9WgRGRmRkhEJxIGFF0CnOjuv89hbCIipa1jx/CQK8Cf/pRsLDnUoEQUHxid4O6vuvt1cXktx7GJiMi554bpxauqYGGdA8Y0apk0zb1qZrvnLBIREVnf974HRx4J330H11+fdDQ5kUki2hP4l5m9H0fLfsvM3sxVYCIiEtXMU3TTTUU5eV4mg54eUv8uIiKSdXvsAfvtB88/D7fcEibTKyKJDnoqIiINVHNV9Oc/h2a6IpJJ0xxmtouZnRmXXXIVlIiI1HLYYbDDDmFK8XvuSTqarGpwIjKzs4EqoFNcxpnZr3MVmIiIpGjSBH7zm7B+1VXgnmw8WZTJFdHJwJ7ufom7XwIMAE7NTVgiIrKeigrYemt4802YODHpaLImk0RkwOqU96tjmYiI5EOLFnDWWWG9iIb9ySQR3Q68FIf4uRyYgiafExHJr1/9CjbfHJ5+Gl4rjnEFMuk1dw1hiJ/FwGeEIX7+nKO4REQknXbt4NR4V6RIhv3JqNccoTnO47Im++GIiEi9zjkHyspC77l0M7o2MhvTa64j6jUnIpKc7t3hF7+A1avDc0WN3Mb0mrtUveZERBJW84DrLbfAkiXJxrKJ1GtORKQx2nVXOOgg+PpruPnmpKPZJBvba+4yQq+5MTmJSkRE6lfzgOt118HKlcnGsgk2ttfcYkKvuWtzFZiIiNRj8GDYeWf4+OMwX1EjlUlnhbHAbHe/3t2vB+aYWYOviMxsjJl9YmbTU8ram9lEM5sVX9vFcjOz682sOk450TflmOFx/1lmNjylvF+cmqI6HqtmQxEpbmZrr4pOOy0MA9SjR6NLSpk0ze3s7l/UvHH3JcBuGRx/BzCkVtmFwDPu3ht4Jr4HOBToHZdK4EYIiQu4lDA30h7ApTXJK+5zaspxtT9LRKT4uIeEtGpVWJ87FyorG1UyyiQRNUn5pV+TFBo8n5G7P09o0ks1FBgb18cCR6SU3+nBFKCtmW1DmBNporsvjolwIjAkbmvj7lPitOZ3ppxLRKR4XXLJ+gOgLl8OI0cmE89GyGRivKsJM7TeF98fDYzaxM8vd/eP4vrHQHlc7wzMT9lvQSzbUPmCNOUiIsVt3rzMygtQJlc0d5rZNOCAWHSku7+drUDc3c0sL+Oam1klocmP8vJyJk+enI+PzciyZcsKMq5cUF2LT6nUE5Kv64BOndhs0aL1yld06sSULMaV03q6e94WoAcwPeX9u8A2cX0b4N24fjNwbO39gGOBm1PKb45l2wDvpJSvs9+Gln79+nkhmjRpUtIh5I3qWnxKpZ7uBVDXcePcW7VyDw10YWnVKpRnUe16AtM8S7kh07Hmsu1RoKbn23DgkZTyYbH33ADgSw9NeE8Cg82sXbxfNRh4Mm5bamYDYm+5YSnnEhEpXhUVMHo0bLtteN+kSXjAtaIi2bgykLdEZGZ3A/8CtjezBWZ2MvAH4GAzmwUcFN8DTABmA9XALcAZAO6+GPgdMDUu/x3LiPvcGo95H3g8H/USEUlcRUWYQrxbN1izBvr0STqijDT4HpGZTQR+4+5vbMwHufuxdWw6MM2+Doyo4zxjSDOig7tPA3bamNhERBo9Mzj4YLjttjB7a9++9R9TIDK5Ivot8Gczuz12lxYRkUJy0EHh9emnk40jQ5kM8fOquw8CHgOeMLNLzaxl7kITEZGMHBgbmF54Ab75JtlYMpDRPaLYEeBdwigGvwZmmdnxuQhMREQytNVWsNtuYQDUF19MOpoGy2SsuX8CC4FrCQ+LngAMBPYws9G5CE5ERDJ08MHhdeLEZOPIQCYjK1QCb8eOBKl+bWYzsxiTiIhsrIMOgj/+sVElokzuEc1Ik4Rq/ChL8YiIyKbYZx9o0QJefx0+/TTpaBokK88RufvsbJxHREQ2UcuWsO++Yf2ZZ5KNpYGSHllBRESyrZHdJ1IiEhEpNjXPE02cuP4UEQVIiUhEpNjsuit07Ajz58OsWUlHUy8lIhGRYtOkydqHWxtB85wSkYhIMWpE94mUiEREilFNIpo0CVatSjaWeigRiYgUo27doHdvWLoUXn456Wg2SIlIRKRY1VwVFfho3EpEIiLFqpHcJ1IiEhEpVoMGhR50U6bAV18lHU2dlIhERIrVllvCHnuEzgqTJycdTZ2UiEREilkjuE+kRCQiUswawX0iJSIRkWI2YAC0bg0zZ8KCBUlHk5YSkYhIMWvWDPbfP6wXaPOcEpGISLFr1y68nngi9OgBVVWJhlNbQSQiM5tjZm+Z2etmNi2WtTeziWY2K762i+VmZtebWbWZvWlmfVPOMzzuP8vMhidVHxGRglFVBfffv/b93LlQWVlQyaggElE0yN13dff+8f2FwDPu3ht4Jr4HOBToHZdK4EYIiQu4FNgT2AO4tCZ5iYiUrJEjYcWKdcuWLw/lBaKQElFtQ4GxcX0scERK+Z0eTAHamtk2wCHARHdf7O5LgInAkDzHLCJSWObNy6w8AU2TDiBy4Ckzc+Bmdx8NlLv7R3H7x0B5XO8MzE85dkEsq6t8PWZWSbiaory8nMkF+KDXsmXLCjKuXFBdi0+p1BMKv64DOnVis0WL1itf0akTUzKIO5f1LJREtI+7LzSzTsBEM3sndaO7e0xSWRET3WiA/v37+8CBA7N16qyZPHkyhRhXLqiuxadU6gmNoK5XXx3uCS1fvrasVSs2u/rqjOLOZT0LomnO3RfG10+Ahwj3eBbFJjfi6ydx94VA15TDu8SyuspFREpXRQWMHh2mDgdo1Sq8r6hINq4UiSciM9vczLaoWQcGA9OBR4Ganm/DgUfi+qPAsNh7bgDwZWzCexIYbGbtYieFwbFMRKS0VVTA+PFhfffdCyoJQWE0zZUDD5kZhHjucvcnzGwqcK+ZnQzMBX4e958AHAZUA8uBEwHcfbGZ/Q6YGvf7b3dfnL9qiIgUsE6dwusnn2x4vwQknojcfTawS5ryz4ED05Q7MKKOc40BxmQ7RhGRRq889vdK03EhaYk3zYmISB506BDmJlq8GL77Lulo1qFEJCJSCsrKYKutwvqnnyYbSy1KRCIipaLmPlGBNc8pEYmIlIoCvU+kRCQiUiqUiEREJFEF2oVbiUhEpFToikhERBKlRCQiIomqSURqmhMRkUSo+7aIiCRKTXMiIpKomiuiTz+FNWuSjSWFEpGISKlo3hzatYPVq8OYcwVCiUhEpJQU4H0iJSIRkVJSgPeJlIhEREpJAXbhViISESklapoTEZFEqWlOREQSpUQkIiKJ0j0iERFJlO4RiYhIotQ0JyIiiUptmnNPNpao6BKRmQ0xs3fNrNrMLkw6HhGRgvLww2AGK1ZA9+5QVZV0RMWViMysDPgLcCjQBzjWzPokG5WISIGoqoLKyrVXQvPnh/cJJ6OiSkTAHkC1u89292+B8cDQhGMSESkMI0fC8uXrli1fHsoTZF4gbYTZYGZHAUPc/ZT4/nhgT3c/s9Z+lUAlQHl5eb/x48fnPdb6LFu2jNatWycdRl6orsWnVOoJjauu+x9wAJbmd76b8dyzz27w2Nr1HDRo0Cvu3j8bcTXNxkkaG3cfDYwG6N+/vw8cODDZgNKYPHkyhRhXLqiuxadU6gmNrK7dusHcuesVW7du9dYhl/Ustqa5hUDXlPddYpmIiIwaBa1arVvWqlUoT1CxJaKpQG8z62lmzYFjgEcTjklEpDBUVMDo0aG3nFl4HT06lCeoqJrm3H2VmZ0JPAmUAWPcfUbCYYmIFI6KisQTT21FlYgA3H0CMCHpOEREpGGKrWlOREQaGSUiERFJlBKRiIgkSolIREQSVVQjK2wMM/sUWP8Jr+R1BD5LOog8UV2LT6nUE0qnrrXr2d3dt8rGiUs+ERUqM5uWreEzCp3qWnxKpZ5QOnXNZT3VNCciIolSIhIRkUQpERWu0UkHkEeqa/EplXpC6dQ1Z/XUPSIREUmUrohERCRRSkQiIpIoJaI8MrPfmdmbZva6mT1lZtvGcjOz682sOm7vm3LMcDObFZfhKeX9zOyteMz1ZmaxvL2ZTYz7TzSzdgnU8yozeyfW5SEza5uy7aIY87tmdkhK+ZBYVm1mF6aU9zSzl2L5PXF6D8ysRXxfHbf3yGcdU+I72sxmmNkaM+tfa1tR1bWh6qpfITOzMWb2iZlNTylL+13K5vc138ysq5lNMrO34//bs2N5snV1dy15WoA2KetnATfF9cOAxwEDBgAvxfL2wOz42i6ut4vbXo77Wjz20Fj+R+DCuH4hcGUC9RwMNI3rV9bEAPQB3gBaAD2B9wnTdZTF9V5A87hPn3jMvcAxcf0m4PS4fkbKv98xwD0J/Ux3ALYHJgP9U8qLrq4N/Peos36FvAD7AX2B6Sllab9L2fy+JlDPbYC+cX0L4L34fzXRuuqKKI/cfWnK282Bmp4iQ4E7PZgCtDWzbYBDgInuvtjdlwATgSFxWxt3n+LhJ38ncETKucbG9bEp5Xnj7k+5+6r4dgphplwIsY1395Xu/gFQDewRl2p3n+3u3wLjgaHxL6kDgPvj8an1Sa3n/cCBSfyV6e4z3f3dNJuKrq4NlLZ+CcdUL3d/Hlhcq7iu71I2v6955e4fufurcf0rYCbQmYTrqkSUZ2Y2yszmAxXAJbG4MzA/ZbcFsWxD5QvSlAOUu/tHcf1joDyrFcjcSYS/iiDzenYAvkhJaqn1/PcxcfuXcf9CUUp1TVVX/Rqjur5L2fy+JiY28e4GvETCdS26ifGSZmZPA1un2TTS3R9x95HASDO7CDgTuDRXsbi7m1lO+ufXV8+4z0hgFVCVixjypSF1leKWy+9SEsysNfAAcI67L029wE6irkpEWebuBzVw1yrCTLKXAguBrinbusSyhcDAWuWTY3mXNPsDLDKzbdz9o3iZ/EmGVWiQ+uppZicAPwYOjJfoUHc9qaP8c0JTQNN4JZC6f825FphZU2DLuH/WZfAzTdUo65oFG6p3Y1PXdymb39e8M7NmhCRU5e4PxuJE66qmuTwys94pb4cC78T1R4FhsYfKAODLeJn8JDDYzNrFXiyDgSfjtqVmNiDeKxgGPJJyrpoeLMNTyvPGzIYAFwA/cfflKZseBY6JvcB6Ar0JNzanAr1jr7HmhBvyj8YENgk4Kh6fWp/Ueh4FPJuS8ApBKdU1Vdr6JRzTxqrru5TN72texc+/DZjp7tekbEq2rrnomaGlzh4rDwDTgTeBvwOdY7kBfyH0NnqLdXtfnUS40V0NnJhS3j+e633gBtaOktEBeAaYBTwNtE+gntWE9uPX43JTyraRMeZ3SelNQ+id817cNjKlvBfhF3g1cB/QIpZvFt9Xx+29EvqZ/pTQDr4SWBS/jEVZ1wz+TdLWr5AX4G7gI+C7+PM8ua7vUja/rwnUcx9CJ6k3U76fhyVdVw3xIyIiiVLTnIiIJEqJSEREEqVEJCIiiVIiEhGRRCkRiYhIopSIREQkUUpEIgkws8vM7Ddx/f828hxtzeyMjTiuh5l9Y2E6kh6WMvXBRpyrZTzPt2bWcWPPI6VNiUgkx+JT6XV+19z9hxt56raEKSI2xvvuvutGHvtv7v5NPM+Hm3ouKV1KRFJSLEwKdnBcv8LM/jfNPsPiJGBvmNnfYtm5ZjY9Luek7FtXeQ8Lk8PdSXjKvKuZjTSz98zsRcIcRjX7Lks5ZqaZ3WJh0rKnzKxl3Pawmb0SyyvjoX8AvhevSK6K+/3SzF6OZTebWVkG/za9zOw1M9s9xvKOmd0RY64ys4PM7J8WJkLbo6HnFalX0kNraNGSz4UwAdpkwjQc/wDKam3fkTA8Tcf4vj3QjzC8yeZAa2AGYfj8tOXxuB7AGmBAfF+zbyugDWFYlN/EbctSjlkF7Brf3wv8siaO+NqSkNg6xP1TJ3LbgTB0VLP4/q/AsDT/Bv8+rmadkBhfA3apFcsPCH+wvgKMIQz5MhR4uNY559T8m2nRkumi0belpLj783EwxnOBge6+utYuBwD3uftncf/FZnY88JC7fw1gZg8C+xJ+Kacrfy2ea66HycSI5Q95HATWzOoaCPQDd389rr9CSAgAZ5nZT+N6V8Igqh/XOvZAQsKbGqpISxo2+vpWhIEpj3T3t2vF8laMdwbwjLu7mb2VEpfIJlMikpJiZj8gTJf8uYcZKnPp6404ZmXK+mqgpZkNBA4C9nL35WY2mTAQam0GjHX3izL8zC+BeYQBMVMTUWosa1Ler0G/OySLdI9ISoaFeVaqCE1Ly+J0FbU9CxxtZh3iMe2BF4AjzKyVmW1OGHH7hQ2Up/N83LelmW0BHJ5B6FsCS2IS+j4wIJZ/BWyRst8zwFFm1qkmdjPr3oDzfxtjH2Zmx2UQl0hW6K8aKQlm1gp4EDjP3Wea2e+AK4EnUvdz9xlmNgp4zsxWA6+5+wlmdgdhCgaAW939tXjetOW1ufurZnYP8AahuWxqBuE/AfzKzGYSppSYEs/5eew8MB143N3PN7P/BJ6KvfS+A0YAc+v7AHf/2sx+DEyMnSfezCA+kU2iaSBESoyZ9QAec/edsnjOOYS5aj7L1jmldKhpTqT0rAa2NLPXN/VENQ+0As0I945EMqYrIhERSZSuiEREJFFKRCIikiglIhERSZQSkYiIJEqJSEREEqVEJCIiiVIiEhGRRCkRiYhIov4/QtIXIVY8j5QAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": { "image/png": { "height": 266, "width": 418 }, "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "#====================================================================\n", "#Solving for different times\n", "#====================================================================\n", "#Time array\n", "times = np.arange( t0, tf, tstep )\n", "\n", "rpos = []\n", "for t in times:\n", " #Finding the new eccentric anomaly\n", " E = optimize.bisect( kepler, a0, b0 )\n", " #Computing coordinates at this time\n", " ri = r(E)\n", " print (\"In %f hours, the satellite is located at (%f,%f) km\"%(t/HR2SC, ri[0], ri[1]) )\n", " rpos.append( ri )\n", "rpos = np.array(rpos) \n", "\n", "#Plotting\n", "plt.plot( rpos[:,0], rpos[:,1], \"o-\", color=\"red\", lw = 2 )\n", "plt.grid(True)\n", "plt.xlabel(\"$x$ coordinate [km]\")\n", "plt.ylabel(\"$y$ coordinate [km]\")" ] }, { "cell_type": "markdown", "metadata": { "cell_id": "00097-c6e00802-c115-4c54-bfb2-da772f8eff96", "colab_type": "text", "deepnote_cell_type": "markdown", "id": "RJlKLGjww7XU" }, "source": [ "- - - " ] }, { "cell_type": "markdown", "metadata": { "created_in_deepnote_cell": true, "deepnote_cell_type": "markdown", "tags": [] }, "source": [ "\n", "Created in deepnote.com \n", "Created in Deepnote" ] } ], "metadata": { "colab": { "include_colab_link": true, "name": "one-variable-equations.ipynb", "provenance": [] }, "deepnote": {}, "deepnote_execution_queue": [], "deepnote_notebook_id": "568ec0db-8924-4d77-814a-ffc3d2455285", "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.2" }, "toc": { "colors": { "hover_highlight": "#DAA520", "running_highlight": "#FF0000", "selected_highlight": "#FFD700" }, "moveMenuLeft": true, "nav_menu": { "height": "426px", "width": "252px" }, "navigate_menu": true, "number_sections": true, "sideBar": true, "threshold": 4, "toc_cell": false, "toc_section_display": "block", "toc_window_display": false } }, "nbformat": 4, "nbformat_minor": 4 }