10. Решение нелинейных уравнений.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import linalg\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Решение нелинейных уравнений\n",
    "\n",
    "Рассмотрим задачу численного нахождения корня $x^*$ нелинейного уравнения\n",
    "$$\n",
    "f(x) = 0, \\qquad a\\le x \\le b, \\qquad f(a)\\cdot f(b)<0.\n",
    "$$\n",
    "Предполагается, что функция $f(x)$ непрерывна на отрезке $[a,b]$ и в концах этого отрезка принимает значения разных знаков.\n",
    "\n",
    "## 1. Некоторые методы\n",
    "\n",
    "### 1.1. Метод деления пополам\n",
    "\n",
    "Метод *деления пополам* (или *дихотомии*, или *бисекции*) состоит в следующем.\n",
    "На первой итерации находим середину отрезка $[a,b]$\n",
    "$$\n",
    "m = \\frac{1}{2}(a+b)\n",
    "$$\n",
    "Если $f(a)\\cdot f(m) < 0$, то полагаем $b=m$, иначе полагаем $a=m$.\n",
    "В обоих случаях повторяем итерации для нового отрезка $[a,b]$, пока величина отрезка не станет меньше заданного значения.\n",
    "\n",
    "Метод деления пополам имеет *линейную скорость сходимости* (также в этом случае говорят, что последовательность сходится со скоростью геометрической прогрессии). Если $x_k$ - середина отрезка на $k$-й итерации, то\n",
    "$$\n",
    "|x_{k+1} - x^*| \\le \\frac{1}{2} |x_k - x^*| \\le \\frac{|b-a|}{2^k}\n",
    "$$\n",
    "\n",
    "\n",
    "### 1.2. Метод Ньютона\n",
    "\n",
    "*Метод Ньютона* (*метод касательных*) на каждой итерации $k$ аппроксимирует в текущей точке $x_k$ функцию $f(x)$ касательной к этой функции в точке $x_k$. Следующим приближением к $x^*$ является точка $x_{k+1}$ пересечения касательной \n",
    "$$\n",
    "y = f(x_k) + (x-x_k) f'(x_k)\n",
    "$$\n",
    "с осью $Ox$:\n",
    "$$\n",
    "x_{k+1} = x_k - \\frac{f(x_k)}{f'(x_k)}.\n",
    "$$\n",
    "\n",
    "<img src='http://www.math.by/newton1.png'>\n",
    "\n",
    "Если $f′(x)\\ne 0$ и $f(x)$ имеет непрерывную вторую производную, то имеет место *квадратичная сходимость* метода Ньютона:\n",
    "$$\n",
    "|x_{k+1} - x^*| \\le C |x_k - x^*|^2, \\quad где \\quad C = |f''(ξ)/(2f'(η))| \n",
    "$$\n",
    "Итак, метод Ньютона работает быстрее, чем метод деления пополам, но требует вычисления производной.\n",
    "\n",
    "\n",
    "### 1.3. Метод секущих\n",
    "\n",
    "Заменим в методе Ньютона производную разделенной разностью: \n",
    "$$\n",
    "f'(x_k)\\approx\\frac{f(x_k) - f(x_{k-1})}{x_k - x_{k-1}}\n",
    "$$\n",
    "получим *метод секущих* (или *метод хорд*):\n",
    "$$\n",
    "x_{k+1} = x_k - f(x_k)\\frac{x_k - x_{k-1}}{f(x_k) - f(x_{k-1})}\n",
    "$$\n",
    "\n",
    "<img src=https://studfile.net/html/2706/298/html_a70QXs6983.h1q0/img-y3y3j8.png>\n",
    "\n",
    "Если $f′(x)\\ne 0$ и $f(x)$ имеет непрерывную вторую производную, то имеет место _сверхлинейная сходимость_:\n",
    "$$\n",
    "|x_{k+1} - x^*| \\le C |x_k - x^*|^r, \\quad где \\quad r = \\frac{1+\\sqrt{5}}{2} = 1.618\\dots, \\quad C = \\sqrt[r]{|f''(ξ)/(2f'(η))|}\n",
    "$$\n",
    "— медленне, чем метод Ньютона, но не нужно вычислять производные.\n",
    "\n",
    "### 1.3. Обратная квадратичная интерполяция\n",
    "\n",
    "Рассмотрим квадратичную функцию\n",
    "$$\n",
    "g(y) = \n",
    " x_{k-2}\\frac{(y-f_{k-1})(y-f_k)}{(f_{k-2}-f_{k-1})(f_{k-2}-f_k)}\n",
    "+x_{k-1}\\frac{(y-f_{k-2})(y-f_k)}{(f_{k-1}-f_{k-2})(f_{k-2}-f_k)}\n",
    "+x_k    \\frac{(y-f_{k-2})(y-f_{k-1})}{(f_k-f_{k-2})(f_k-f_{k-1})}\n",
    "$$\n",
    "Имеем $x_k=g(f_k)$, $x_{k-1}=g(f_{k-1})$, $x_{k-2}=g(f_{k-2})$.\n",
    "\n",
    "В качестве очередного приближения берем\n",
    "$$\n",
    "x_{k+1} = g(0)\n",
    "$$\n",
    "Метод сходится чуть быстрее, чем метод секущих:\n",
    "$$\n",
    "|x_{k+1} - x^*| \\le C |x_k - x^*|^r, \\quad где \\quad r = 1.839\\dots\n",
    "$$\n",
    "но необходимы три начальные точки\n",
    "\n",
    "### 1.4. Метод Брента\n",
    "\n",
    "— сочетает методы деления пополам, секущих и обратной квадратичной интерполяции\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Методы из модуля `optimize`\n",
    "\n",
    "В модуле `optimize` реализован ряд методов для решения нелинейного уравнения.\n",
    "Наиболее общей является функция `minimize_scalar()` В ней используемый метод можно указать явным образом задав параметр `method`. Но можно также воспользоваться \"специальными\" методами. В частности,\n",
    "  * `bisect(f, a, b)` - реализует метод деления пополам; \n",
    "  * `brentq(f, a, b)` - модификация метода Брента;\n",
    "  * `newton(fun, x0, fprime)` - метод Ньютона, где `x0` - начальное приближение, `fun` - функция, вычисляющая $f(x)$, `fprime` - функция, вычисляющая $f'(x)$; если `fprime` не задано, то используется метод секущих.\n",
    "  \n",
    "В качестве дополнительных параметров можно задавать, в частности,\n",
    "  * `xtol` - абсолютная ошибка,\n",
    "  * `rtol` - относительная ошибка,\n",
    "  * `maxiter` - максимальное количество итераций.\n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import bisect, newton, brentq "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Примеры и эксперимент\n",
    "\n",
    "Давайте поэкспериментируем с этими методами. В качестве тестового возьмем уравнение\n",
    "$$\n",
    "\\sin(x) = 0,\\qquad  3\\le x \\le 4.\n",
    "$$\n",
    "Вначале зададим функцию, вычисляющую $f(x)$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(x):\n",
    "    return np.sin(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "И ее производную (понадобится для метода Ньютона):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fprime(x):\n",
    "    return np.cos(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Например, вызовем метод `brentq`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.141592653589793"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "brentq(f, 3, 4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Можно получить дополнительную информацию (сошелся ли метод, количество вычислений функции, количество итераций):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3.141592653589793,       converged: True\n",
       "            flag: 'converged'\n",
       "  function_calls: 7\n",
       "      iterations: 6\n",
       "            root: 3.141592653589793)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x0, r = brentq(f, 3, 4, full_output = True)\n",
    "x0, r"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3.141592653589793,       converged: True\n",
       "            flag: 'converged'\n",
       "  function_calls: 6\n",
       "      iterations: 3\n",
       "            root: 3.141592653589793)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x0, r = newton(f, 3, fprime=fprime, full_output=True)\n",
    "x0, r"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Нам будет интересно в каких точках метод вычисляет $f(x)$ на каждой итерации, поэтому добавим сохранение этих точек в списке `x_all`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f_ext(x):\n",
    "    x_all.append(x)\n",
    "    return np.sin(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.1. Метод деления пополам"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[3.0,\n",
       " 4.0,\n",
       " 3.5,\n",
       " 3.25,\n",
       " 3.125,\n",
       " 3.1875,\n",
       " 3.15625,\n",
       " 3.140625,\n",
       " 3.1484375,\n",
       " 3.14453125,\n",
       " 3.142578125,\n",
       " 3.1416015625,\n",
       " 3.14111328125,\n",
       " 3.141357421875,\n",
       " 3.1414794921875,\n",
       " 3.14154052734375,\n",
       " 3.141571044921875,\n",
       " 3.1415863037109375,\n",
       " 3.1415939331054688,\n",
       " 3.141590118408203,\n",
       " 3.141592025756836,\n",
       " 3.1415929794311523,\n",
       " 3.141592502593994,\n",
       " 3.1415927410125732,\n",
       " 3.1415926218032837,\n",
       " 3.1415926814079285,\n",
       " 3.141592651605606,\n",
       " 3.1415926665067673,\n",
       " 3.1415926590561867,\n",
       " 3.1415926553308964,\n",
       " 3.1415926534682512,\n",
       " 3.141592654399574,\n",
       " 3.1415926539339125,\n",
       " 3.141592653701082,\n",
       " 3.1415926535846666,\n",
       " 3.141592653642874,\n",
       " 3.1415926536137704,\n",
       " 3.1415926535992185,\n",
       " 3.1415926535919425,\n",
       " 3.1415926535883045,\n",
       " 3.1415926535901235,\n",
       " 3.1415926535901235]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_all = []\n",
    "x0 = bisect(f_ext, 3, 4)\n",
    "x_all.append(x0)\n",
    "x_bisect = x_all\n",
    "x_bisect"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Обратите внимание, что в методе деления пополам очередная цифра в ответе генерируется примерно за одно и то же счисло итераций (4 итераций). **Это по сути дела и означает линейную скорость сходимости.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.2. Метод секущих"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[3.0,\n",
       " 3.0004000000000004,\n",
       " 3.142542483087832,\n",
       " 3.1415894937885818,\n",
       " 3.1415926535902665,\n",
       " 3.1415926535897936]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_all = []\n",
    "x0 = newton(f_ext, 3)\n",
    "x_all.append(x0)\n",
    "x_secant = x_all\n",
    "x_secant"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.3. Метод Ньютона"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[3.0, 3.142546543074278, 3.141592653300477, 3.141592653589793]"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_all = []\n",
    "x0 = newton(f_ext, 3, fprime=fprime)\n",
    "x_all.append(x0)\n",
    "x_newton = x_all\n",
    "x_newton"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.4. Метод Брента"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[3.0,\n",
       " 4.0,\n",
       " 3.157162792479947,\n",
       " 3.14154625558915,\n",
       " 3.1415926554589646,\n",
       " 3.141592653589793,\n",
       " 3.1415926535907945,\n",
       " 3.141592653589793]"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_all = []\n",
    "x0 = brentq(f_ext, 3, 4)\n",
    "x_all.append(x0)\n",
    "x_brent = x_all\n",
    "x_brent"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Теперь для каждого метода построим график зависимости абсолютной ошибки от номера итерации."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAD4CAYAAAAKA1qZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3dd3zU9f3A8dc7l8skAbIgkDDCDhASCFOWIEOluAWkoHUgKmodbbG1rXXU9qe/aq1UREBwVKH8HKgoRfaGsGRD2GEmYSSQnXx+f1yShpCES3KXS3Lv5+NxD7jv3X3vfV/l3vdZ748YY1BKKeWePFwdgFJKKdfRJKCUUm5Mk4BSSrkxTQJKKeXGNAkopZQb83R1AJUREhJiWrVq5eowlFKqTtmyZUuKMSa0rMfqVBJo1aoVCQkJrg5DKaXqFBE5Vt5j2h2klFJuTJOAUkq5MU0CSinlxurUmIBSyn3l5uaSlJREVlaWq0OptXx8fIiIiMBqtdr9Gk0CSqk6ISkpiYCAAFq1aoWIuDqcWscYQ2pqKklJSbRu3dru12l3kFKqTsjKyiI4OFgTQDlEhODg4Eq3lDQJKKXqDE0AFavK9dEkAKz910I2/We9q8NQSqka5/ZJYNuilTR85QVO/eGP/GbBT6Rl5bo6JKVULWWxWIiNjaVbt250796ddevWAXDq1Cnuvvtuh73PihUris/tbG6dBC6cSeXK719AjKHDxROsWL2D4X9bxfJ951wdmlKqFvL19WX79u3s2LGD119/nRdeeAGAZs2asWDBAoe9jyaBGlBQUMCayc/QMOMSBb/9EwAzW10m0NeTX8zZzHPzd3ApQ1sFSqmypaWl0bhxYwCOHj1Kly5dANi9eze9evUiNjaWmJgYDh48CMAnn3xSfPzRRx8lPz8fgB9++IHu3bvTrVs3hg4dytGjR5k+fTpvvfUWsbGxrF692qmfw22niP7nL9Npu28zR+5+kFsm3sPhrz7Dd+MqvvlkEu8uS+SfKw6x+mAyr93RlWHRTVwdrlKqhD99s5s9p9Ices7oZoH88WedK3xOZmYmsbGxZGVlcfr0aZYtW3bNc6ZPn87TTz/N+PHjycnJIT8/n7179zJv3jzWrl2L1Wrl8ccf59NPP+Xmm2/mkUceYdWqVbRu3Zrz588TFBTE5MmTadCgAc8//7xDP2NZ3DIJ7F2zhfBP3uNQ667c/NKzAAQOH0Hy22/jkXyO54Z3YETnpjz/7x088lECd8Q15893dMXXy+LiyJVSrlTUHQSwfv16Jk6cyK5du656Tt++fXnttddISkrizjvvpF27dixdupQtW7bQs2dPwJZMwsLC2LBhAwMHDiye1x8UFFSzHwg3TAKXL6Zx+vnn8fX2o9f7f8fiaftiDxgxnOS33yZ9yRKCJk6kS/OGLJzSn2nLE3ln2UGOpV5h9gM9aeTn5eJPoJS63i/2mtC3b19SUlJITk6+6vh9991H7969+e677xgxYgQzZ87EGMP999/P66+/ftVzFy5c6PJpr243JrB0ygs0uXgWj9/9ibAW4cXHvVu3xrtDB9IW/6f4mJenB88Ma89747uz61Qad09fz6mLma4IWylVy+zbt4/8/HyCg4OvOn748GGioqJ46qmnGD16ND/99BNDhw5lwYIFnDtnm3Ry/vx5jh07Rt++fVm5ciVHjhwpPg4QEBBAenp6jXwOt0oCy/75Ce0TlnF42F30vnvkNY8HjBhO5tat5J69enbQyC7hfPRgL85eyuKu99Zx4GzN/MdRStUuRWMCsbGxjBkzhrlz52KxXN1NPG/ePLp06UJsbCz79u1j4sSJREdH8+qrrzJ8+HBiYmIYNmwYp0+fJjQ0lBkzZnDnnXfSrVs3xowZA8DPfvYzvvzyyxoZGBZjjFPfwJHi4+NNVTeVOfrTAVLvu5fk0EiGfr8Aq4/3Nc/JPnyYw7fcSpMXXyTo5+OveXzv6TTun72J7LwCZj8QT4+WNd9/p5S72rt3L506dXJ1GLVeWddJRLYYY+LLer5LWwIiEiUis0TEcRNsy3Dy+DaOTBmHRfLp8csbse6ZD9s+Kbx9artdPI53VBTe7dqS/sMPZZ6nU3gg//dYP4L9vbjvg438uOesM8NWSimnq3ISEJHZInJORHaVOj5SRPaLSKKITK3oHMaYw8aYh6oag722fv4hTc9lkNn/IqHb/gQLp8DXTxTeHrfdvrXNEgoYPoKMLVvIKzXYUyQyyI9/T+5Lx/BAHv1kC/M3n3B2+Eop5TTVmR00B3gX+KjogIhYgGnAMCAJ2CwiCwEL8Hqp1z9ojKmRpbnDn3qFMR6r6Nv7Lvp2fRQo7AIr6gpb9Qbs+j/IyyZw5AhSpk0j/ccfaTxuXJnnC27gzb8e7s1jn27l1//3E/nGMK5Xi5r4KEop5VBVbgkYY1YB50sd7gUkFv7CzwE+B24zxuw0xowqdbMrAYjIJBFJEJGE0lOx7OXt05CwmF5sOL8bGkVCoxa2W+OWtluHWyA3A5I249W2LV5RUaT9sLjCc/p7ezLr/ngGtg/lpYW7STyng8VKqbrH0WMCzYGS/SNJhcfKJCLBIjIdiBORF8p6jjFmhjEm3hgTHxoaWuXA+jTrQ+LFRJIzykgkrW4AscDhFYgIgSNHkLF5M3mpqRWe02rx4M17YvD39uTpz7eTk1dQ5fiUUsoVHJ0Eylr1UO70I2NMqjFmsjGmjTGmdHeRQ/UO7w3AxjMbr33QpyE07w6HVwIQMGIEFBSQ/uPS6543LMCHv94Vw+5TafxtyQGHxqyUUs7m6CSQBESWuB8BnHLwe1RJx8YdaejdkA2nNpT9hKjBcHILZF3Cu317vFq2JH1x2bOEShsW3YRxvVrw/qpDbDhcceuhpPyCujM9Vyll89prr9G5c2diYmKIjY1l48Yyflg62Jw5czh1yjlfpY5OApuBdiLSWkS8gLHAQge/R5VYPCz0atqLDac3UObaiKjBYPLh6FpEhICRI7mycRN5Fy7Ydf7fj+pEq2B/W/XRzIqrjxpj+HjDMbr8cTF/+X4fBZoMlKoT1q9fz7fffsvWrVv56aef+PHHH4mMjLz+C6upViYBEfkMWA90EJEkEXnIGJMHTAEWA3uB+caY3Y4Jtfr6hPfhbMZZjqUdu/bBiJ7g6QtHbF1CgSOGQ34+6T/+aNe5/bw8eWtMLGfSsvjD17vKfV5GTh7PzNvO77/aRZNAb6avPMSTn20jKze/Sp9JKVVzTp8+TUhICN7etsWmISEhNGvWjC1btjBo0CB69OjBiBEjOH36NACJiYncdNNNxZvQHDp0iMuXLzN06FC6d+9O165d+frrrwFbOepOnTrxyCOP0LlzZ4YPH05mZiYLFiwgISGB8ePHExsbS2amY0vXuM2KYYDjace59ctb+V3v3zG249hrn/DxnZB2Ep7YiDGGQyNG4tWiBS1mfmD3e/xj6UH+d8kB/j42lttirx4TP5R8mcc+2cLBc5d59qb2PHFjW2auOcyfF+2je4tGfDAxnuAG165kVkqVWgn7/VQ4s9Oxb9C0K9z8lwqfcvnyZfr3709GRgY33XQTY8aMoV+/fgwaNIivv/6a0NBQ5s2bx+LFi5k9eza9e/dm6tSp3HHHHWRlZVFQUICXlxcZGRkEBgaSkpJCnz59OHjwIMeOHaNt27YkJCQQGxvLvffey+jRo/n5z3/O4MGDefPNN4mPL3PR71Xq1IrhmhYZEEkz/2ZsOF3BuEDyPkg7bZslNGI4VzZsIP/iRbvf47HBbejRsjEvfrWLkyWKzX370ylG/2MNKZdz+OjBXjw5tB0eHsKkgW345/ju7D6Vxp3vreNQ8uXqfUillNM0aNCALVu2MGPGDEJDQxkzZgzvv/8+u3btYtiwYcTGxvLqq6+SlJREeno6J0+e5I477gDAx8cHPz8/jDH89re/JSYmhptuuomTJ09y9qyt+kDr1q2JjY0FoEePHhw9etTpn8mtSkmLCH2a9WHJsSXkF+Rj8Si1P0DUINufR1ZBtzEEDB9B6gczSV+6jEZ33WnXe3haPHjr3lhueWc1z87bztwHe/HXH/bx4dqjxLVoxLT7utOske9Vr7mlazhNG/rwyNwE7vznOmZM6EHvqOBy3kEpdb1f7M5ksVgYPHgwgwcPpmvXrkybNo3OnTuzfv36q56Xllb2pjeffvopycnJbNmyBavVSqtWrcjKygIo7mYqeh9Hd/2Uxa1aAgC9m/YmPSedvef3Xvtgk67gGwSHVwDg06Uz1ubNSftPxQvHSmsR7MdLozuz8ch5Br+xgg/XHuWBfq2YN6nvNQmgSPcWjfny8RsIbuDFhFmb+Hr7ycp+NKWUk+3fv794u0iA7du306lTJ5KTk4uTQG5uLrt37yYwMJCIiAi++uorALKzs8nIyODSpUuEhYVhtVpZvnw5x46VMUZZijNLS7tdEugV3gug7C4hDw9ba+DwCjDGNktoxAiurFtP/uXKddPc1b05o2LCScvK5R/j4nhpdGe8PCu+3C2C/fjisX7EtWjE059vZ9aaI5V6T6WUc12+fJn777+f6OhoYmJi2LNnDy+//DILFizgN7/5Dd26dSM2NrZ4k/iPP/6Yd955h5iYGPr168eZM2cYP348CQkJxMfH8+mnn9KxY8frvu8DDzzA5MmTdWC4ugPDRe5aeBeNfRozc/jMax9M+BC+/SVMSYCQdlxZt47jDz5E5KyZNLjhhkq9T36BISMnjwAfa6Vel52Xz+OfbGVNYgqrf3MjYQE+lXq9UvWRlpK2jw4M26F3eG+2nd1GVl7WtQ9GDbb9WdQlFBMDHh5kbtte6fexeEilEwCAt6eF34+KJje/gJmrtTWglHIet0wCfcL7kFOQw7Zz2659MKg1NGpZnAQsDRrg3b49mVu31miMrUL8Gd2tGZ9sOMb5Kzk1+t5KKffhlkkgvkk8nuLJxtPlLPeOGgRHVkOBbQGXb1wsmTt2YPJrdkHXEze2JTM3n9mVHBvYfeoSx1MznBSVUqo+ccsk4Gf1IyY0puL1AtmX4JStC8gvLo6CK1fITkyssRgB2jUJ4OYuTZm77uh1S1EUOZR8mTv/uY6b3lrJtOWJ5OZrZVOlVPncMgmArUtoT+oeLmVfuvbB1kXrBVYA4BsXB0DmtjK6j5zsiRvbkp6dx9x1R6/73Lz8Ap6bvwMfq4UhHcJ4Y/F+fvaPNew4Yf9iN6WUe3HfJNCsDwbD5jObr33QP8S2ZqBwXMAaEYElNMQlSaBzs4bc1CmM2WuPcDk7r8Lnzlh9mO0nLvLybZ2ZPqEHMyb04GJGLnf8cy2vfLuHjJyKX6+Ucj9umwS6hHTBz9Ovgi6hQXB8I+RmIiL4xcaRsbXmkwDAlCHtuJiRyycbyl9Usv9MOm8vOcjNXZoyulszAIZ3bsqSZwcyvndLZq05wvC3VrHyQNV2Z1NK2aoOPPfcc8X333zzTV566aVqnXP79u0sWrSompFVndsmAauHlfim8RUkgRshPxuO2x73jYsj98SJcjegd6bYyEYMaBfCzNWHycy5dnA6N7+AZ+dvJ8DHk1dv74LIf/f2CfCx8srtXfj35L54e3pw/+xNPDt/O1eu06pQSl3L29ubL774gpSUFIedU5OAC/UJ78OxtGOcvnz62gdb9gUPa3GXkG+crahTxvbKrxdwhCeHtCPlcg6fbTp+zWPTliey+1Qar93RpdwqpD1bBbHo6QE8NaQtX207yV3vrSPpgs4gUqoyPD09mTRpEm+99dY1jyUnJ3PXXXfRs2dPevbsydq1awHo2rUrFy9exBhDcHAwH330EQATJkzgP//5D3/4wx+YN28esbGxzJs3j/Pnz3P77bcTExNDnz59+OmnnwB46aWXePDBBxk8eDBRUVG88847jvlMDjlLHdUnvA9gKyFxR7s7rn7Qyx8ie/130VjnzojVSua27QQOG1bDkUKv1kH0bh3E+6sOMb5PC7w9bcXvdp28xLvLErktthkju4RXeA5vTwvPDu9Aj1ZBTPnXVm57dy3vT+hBfKugmvgISjnMXzf9lX3n9zn0nB2DOvKbXr+57vOeeOIJYmJi+PWvf33V8aeffppnnnmG/v37c/z4cUaMGMHevXu54YYbWLt2LS1btiQqKorVq1czceJENmzYwHvvvcfLL79MQkIC7777LgBPPvkkcXFxfPXVVyxbtoyJEyeyvfDH5759+1i+fDnp6el06NCBxx57DKu18gtSS3LrlkDbRm0J9gmueKro6R2QcR4PLy98unRxyeBwkSeHtONsWjb/TkgCbOUlnpu/gyB/L/40urPd5xnUPpQvH7+BAB9Pxn2wgfkJJ5wVslL1TmBgIBMnTrzml/iPP/7IlClTiI2NZfTo0aSlpZGens6AAQNYtWoVq1at4rHHHmPnzp2cPHmSoKAgGjRocM3516xZw4QJEwAYMmQIqampXLpkm8V466234u3tTUhICGFhYcUlqKvDpS0BEekEPA2EAEuNMe/V8PvTO7w3G0/bNpEp2ZcO2KaKLn8Njq6G6Nvw7R7HhY8+piA7Gw/vmt/85Ya2wcS1aMR7Kw4xpmckf//xIPvPpjP7gXga+XlV6lxtwxrw1RM3MOVf2/j1gp84cCadF27phMVDrv9ipVzMnl/szvTLX/6S7t2784tf/KL4WEFBAevXr8fX9+pKwQMHDmTatGkcP36c1157jS+//JIFCxYwYMCAMs9dVj23ou+m0qWm8/KqP7ZXne0lZ4vIORHZVer4SBHZLyKJIjK1onMYY/YaYyYD9wLX3zLHCfqE9yE1K5XEi2UsBGveHbwCiruE/OLiMLm5ZO3eU7NBFhIRnhzSlpMXM3n5mz1MX3mIe3pEMKRjkyqdr5GfF3N+0ZMH+rVi5pojPDhnM2lZ9i1KU8qdBQUFce+99zJr1qziY8OHDy/u0gGKu3AiIyNJSUnh4MGDREVF0b9/f958883iJFC6TPTAgQP59NNPAVixYgUhISEEBgY67bNUpztoDjCy5AERsQDTgJuBaGCciESLSFcR+bbULazwNaOBNcDSasRSZSXHBa5hsUKr/nDYtu+wb+GOP67sErqxQxidmwXy8YZjNA304fc/i67W+TwtHrw0ujN/vqMraxNTuGPa2qt2RFNKle255567apbQO++8Q0JCAjExMURHRzN9+vTix3r37k379u0BGDBgACdPnqR///4A3HjjjezZs6d4YPill14qPs/UqVOZO3eucz+IMabKN6AVsKvE/b7A4hL3XwBesPNc35VzfBKQACS0aNHCOMOIBSPMs8ufLfvB9f805o+Bxlw8YYwx5uCw4ebElClOicNey/aeNdG//96sOnDOoeddl5hi2v9ukfnj17scel6lHGHPnj2uDqFOKOs6AQmmnO9eRw8MNwdKjjImFR4rk4gMFpF3ROR9oMyJssaYGcaYeGNMfGhoqGOjLRTuH05KZjnzflv0tf1ZuF7ALy6WjG3by+y3qyk3dgxj2x+GM6CdY69H3zbB9GsTzCpdUKaU23B0EihrVLHcb0tjzApjzFPGmEeNMdMcHIvdgn2DOZ91vuwHm3QBqz+csFUc9Y3rTn5KCrknXDuj5nq7lFXVwPahHE65olVIlXITjv4mSQIiS9yPAE45+D0cLsgniNSs1LIftHhCRPxVK4fBteMCzjSova11sfJg5VoDxhgW7TzNufQyNupRStVajk4Cm4F2ItJaRLyAscBCB7+HwwX7BJOek05Ofjmbt7ToA2d3QXY63m3b4NGgARn1NAm0DvEnMsiXlfsrlwRWHkjm8U+3MviNFfxtyQHSdZaRUnVCdaaIfgasBzqISJKIPGSMyQOmAIuBvcB8Y8xux4TqPEG+thWz5XYJtegDpgCSNiMWC77dulVpu8m6QEQY1D6UdYdSyMmzfy+C7346TYC3Jzd2COOdpQcZ/MYK5qw9UqlzKKVqXpWTgDFmnDEm3BhjNcZEGGNmFR5fZIxpb4xpY4x5zXGhOk+wTzBA+V1CET1BPGxVRbF1CWUfOED+5cs1FWKNGtQ+jIycfBKOlZMUS8nNL+A/e84yLLoJ08Z356snbqB9kwBe+mYPQ/+2gq+3n6SgwHUD6Uqp8rl12YgiQT6FLYHMcr70vAOgSWc4UTQuEAvGkLl9R02FWKP6tgnGahG7y06vTUzhUmYut3S11S6KjWzEvx7pzdwHe9HA28rTn2/nZ++uYXUlxxmUqm0sFguxsbF069aN7t27s27dOoed+6uvvmLPnppfiKpJANvsIKigJQC2qaInNkN+Hr7duoGHR70dHG7g7Ul8yyC7xwUW7bR1BQ1oH1J8rKhb6bsn+/P2mFguZeYyYdYmJs7exL4zac4KXSmn8vX1Zfv27ezYsYPXX3+dF1544Zrn5FdxL3JNAi5U1B1U7pgAQGRvyL0CZ3dhadAA7/bt620SABjUIZR9Z9I5m1bxbJ+irqCbopsUVzYtycNDuD2uOUufG8SLt3Zi+/EL3PL31fxmwU/XPbdStVlaWhqNGzcGbOUdbrzxRu677z66du0KwCeffEKvXr2IjY3l0UcfLU4ODRo04He/+x3dunWjT58+nD17lnXr1rFw4UJ+9atfERsby6FDh2rsc7h1KekiflY/fD19Sc2sqCVgKy/BiY3QLBbfuFjSFn6Dyc9HLNd++dV1g9qH8pfv97HqQDL3xEeW+7x1h1K5mJHLzV2aVng+b08LDw+I4u4eEby7LJG564+ycMcpHhkYxaMDo/D31v8Vlf3O/PnPZO91bClp704dafrb31b4nMzMTGJjY8nKyuL06dMsW7as+LFNmzaxa9cuWrduzd69e5k3bx5r167FarXy+OOP8+mnnzJx4kSuXLlCnz59eO211/j1r3/NBx98wIsvvsjo0aMZNWoUd999t0M/1/VoS6BQkE9QxS2BhhEQGAHH1wO2YnIFV66QnVhG4bl6oGPTAMICvK87LrDop9M08PZkYHv7Vi838vPixVHR/PjsIIZ0KpxJ9OYKPtt0nHwdPFa1XFF30L59+/jhhx+YOHFicfWAXr160bp1awCWLl3Kli1b6NmzJ7GxsSxdupTDhw8D4OXlxahRowDo0aMHR48edclnKaI/vwoF+wRX3BIAW2vg2Dow5r+LxrZuxadDhxqIsGaJCAPbh7Jkz1nyC0yZJaZz8wtYvOcMQzuF4WOtXGuoZbA/0+7rzkP9L/Dn7/bywhc7+Xj9MV4a3ZlerXWTG1Wx6/1irwl9+/YlJSWF5MItZ/39/YsfM8Zw//338/rrr1/zOqvVWlwa2lHloKtDWwKFgnwrWDVcpEUfSD8Fl05gjYjAEhpSbxeNga1L6FJmLjuSLpb5+IbDtq6gollBVdG9RWP+Pbkv794Xx4WMHO59fz1PfbaNM5d0vEDVbvv27SM/P5/g4OBrHhs6dCgLFizg3LlzAJw/f55jx45VeL7SJaVriiaBQsE+FdQPKhLZ2/bn8Q2ICH6xcfV20RhA/7YheAjlzhJatPM0/l6W4lITVSUijIppxtLnBvHUkLb8sPsMQ/53BdOWJ5KVW7WZFko5Q9GYQGxsLGPGjGHu3LlYyhgTjI6O5tVXX2X48OHExMQwbNgwTp8uYy/zEsaOHcsbb7xBXFycDgy7QpBPEBeyLlBgCvCQcnJjk862TWaOb4CYe/GNiyN9yRLykpPxdFKFU1dq7O9Ft8hGrDyQzDPD2l/1WF5+AYt3n2VopyaV7goqj5+XJ88O78A98ZG8+t0e3li8n/kJJ/j9rdEM7RR27c5vStWw8qZ/Dh48mMGDB191bMyYMYwZM+aa514uscj07rvvLh4IvuGGG3SKqCsF+waTb/K5lH2p/Cd5WCCyZ4mKorZNZjK219/WwKD2oexIusiFK1fXVdpw+Dznr+RUqyuoPJFBfrw/IZ6PH+qF1eLBwx8l8MhHWzinU0qVcjhNAoWKS0dcb3A4sg+c3Q2ZF/Hp3Bk8Pcnauavi19Rhg9qHYgysTrx6v4XvCruCBndwXgtoQLtQvn96AL+9pSOrDyZz099W8sXWJJfu5aBUfaNJoFBx6YjrjQu06A0YSErAw8sLa7Nm5Ca5dm8BZ4qJaEQjP+tV4wK2rqAzDHFgV1B5rBYPJg1sw6KnB9CuSQDPzt/BIx8laKvATekPgIpV5fpoEihkV+kIgObxIJbiOkJeEc3JSTrp7PBcxuIhDGgXysoDycVF4DYeKewKus4CMUdqE9qA+Y/25cVbO7H6YIq2CtyQj48Pqamp+t+8HMYYUlNT8fHxqdTrdGC4kN0tAe8G0LRr8SYz1ohIspYscXZ4LjWofSjf7DjF3jNpdG7WkO92nsbXamFwh7AajcPiITw8IIohHcP41YKfeHb+DhbtPM2f7+hKWGDl/sdXdU9ERARJSUnF8/LVtXx8fIiIiKjUazQJFGro3RCLWK4/JgC29QJb5kJ+LtaICPIvXCD/8hUsDfyv/9o6aGA7W2G4VQdS6NAkgMW7zjCkUxi+Xq4plxFV2Cr4cO0R3li8nxFvr2L2Az2Ja9HYJfGommG1WotX5CrH0e6gQh7icf3SEUVa9IG8TDjzE16RtqybezLJyRG6TligD53CA1l54Bybjpwn9UoOtzphVlBlFLUKFj09gAAfK+NnbmTNwZTrv1ApdRWXJgERGSwiq0VkuogMdmUsULjXsD0tgcjCYnLHN2ItbHrlJtXfJAC2LqGEoxeYn3ACX6uFG2u4K6g8bUIbsGByX1oE+fHgnM18v7PiBTlKqatVZ3vJ2SJyTkR2lTo+UkT2i0iiiEy9zmkMcBnwwbZJvUsF+wZff2AYIDAcGrWA4+uLk0DOifo7QwhsSSCvwPDV9lMM6ei6rqCyhAX6MG9SX7o0D+SJf23l803HXR2SUnVGdVoCc4CRJQ+IiAWYBtwMRAPjRCRaRLqKyLelbmHAamPMzcBvgD9VIxaHsLs7CGytgRMbsTRsiIe/P7n1eIYQQI+WjfEv/OJ3xgKx6mroZ+WTh3vTv10oU7/YyfSVNbfsXqm6rDp7DK8CSn9j9gISjTGHjTE5wOfAbcaYncaYUaVu54wxRbuQXwC8y3ofEZkkIgkikuDsWQFFlUTtmoLWog9cPotcPIY1IrI/xRMAACAASURBVKLedwd5eXpwQ9sQfKwe3NixdpbI8PPyZObEeEbFhPOX7/fxl+/36XRCpa7D0bODmgMl+0WSgN7lPVlE7gRGAI2Ad8t6jjFmBjADID4+3qn/ooN8g8jKzyIzLxM/q1/FTy7aZOb4BqyREeS4uCZ4Tfj9qGjOpGXh51V7J5V5eXrw97FxNPS1Mn3lIS5m5PDaHV3LLIWtlHJ8EijrX1q5X9zGmC+ALxwcQ5WVLB1x3SQQ2gm8G8KJDXg1j+DKmrUYY+p1kbPIID8ig65zXWoBi4fw6u1daOznxbvLE7FaPHjl9i6uDkupWsnRs4OSgJJ7EUYApxz8Hk5TtGDMrsFhDw+I7GWbIRQZicnKIj9FpyjWFiLC8yM68OjAKD7ecIxPNlRcy10pd+XoJLAZaCcirUXECxgLLHTweziN3aUjirToDcl7sYY2BCCnno8L1EW/HtmRIR3DeGnhbtYd0iStVGnVmSL6GbAe6CAiSSLykDEmD5gCLAb2AvONMbsdE6rz2V06okjhegEvT1vSqO+Dw3WRxUP4+9hYWoX48/inWzmWesXVISlVq1RndtA4Y0y4McZqjIkwxswqPL7IGNPeGNPGGPOa40J1PrvLSRdp3gMQrHm2eemaBGqnAB8rMyfGYww8PDeB9KxcV4ekVK2hZSNKsFqsBHoF2t8S8PIDv2A8clLxDA3V7qBarFWIP++N787hlCv88vPt5BeUP9EsOT2bl7/ZQ/dXlvDNjjozpKVUlWgSKMXu0hFF/EPhSrJtrcAJTQK1Wb+2Ibz0s2iW7jvH/yzed83jFzNy+OsP+xj4P8uZu/4o/t4Wnpm3nR/3nK35YJWqIbV3wreL2F06ooh/CGSkYo2IIXPLFucFphxiQt9W7DuTzvsrD9OhSQB3do8gPSuX2WuOMnP1YS7n5DG6WzOeHtqO0ABvxs/cyOP/2sqcB3rSr22Iq8NXyuE0CZQS5BNE4sVE+1/gHwJnduIVGUHad99hcnMRq9V5Aapqe2l0Zw4lX2bqFzvZfyadeQknuJiRy4jOTXhmWHs6Ng0sfu7cX/Ri7IwNPPxRAp883JvuWq5a1TPaHVRKUekIu/mFwJUUrM0joKCA3NNaxbK2s1o8+Of4HjQJ9Ob9VYfpFtGIhVNu4P0J8VclAIDG/l58/FAvwgK8eWD2JnafuuSiqJVyDk0CpQT5BpGWk0Zuvp0zSPxDIesi1ma2rRbrezXR+iLI34sFk/vx7ZP9mftgL2IiGpX73LBAHz55uDcNvD2ZOGsTh5Iv12CkSjmXJoFSiqaJ2j1DyN/2fK9gWzmF+l5NtD5pEuhDl+YN7XpuRGM/Pnm4NyLw85kbOXE+w8nRKVUzNAmUUvkkYKuo6ekvYLXqWoF6LCq0AR8/1JuMnHzGz9zI2bQsV4ekVLVpEiil0qUj/GwzRiQrFWuzcHKStDuoPusUHsicX/Qk9XI2w99axT+WHiRNF5+pOkyTQCmVLh1R2BLgSipeEZHaHeQG4lo0ZsFj/ejZqjH/u+QA/f+yjL//eJBLmZoMVN2jSaCU4paAvTOE/AvnjhcvGNOWgDvoFB7IzPt78u2T/ekdFcxbPx6g/1+X8daSA5oMVJ2iSaAUP08/vC3e9rcEfBqBWCAjBWtEc/IvXiT/ss4ecRddmjfkg4nxfPdUf/q1CebvSw/S/y/L+NuSA1zOznN1eEpdlyaBUkSkcmsFPDzALxiuJOMVadtKQQeH3U/nZg15f0I8i54aQP92Ibyz9CBD3lzBF1uTKKigTpFSrqZJoAyVLx0RCldSbQvG0CTgzqKbBfLez3vw5eP9CG/ky7Pzd3DX9HXsOHHR1aEpVSZNAmUI8gmyvzsIbGsFriTjFWlLAlpNVMW1aMyXj/XjzXu6kXQhk9umreX5f+/gXLpOK1W1iyaBMgT7VrJ0hH8oZKTg0bAhHg0aaDVRBYCHh3B3jwiWPTeIRwdF8fX2kwx5cyXvrzxEdl6+q8NTCnBxEhCRASIyXURmisg6V8ZSUpBPEBeyLlBgCux7gX8oXElBRLBGRmp3kLpKgI+VF27uxH+eGUTv1kG8/v0+bvjLMt5YvI+kC7ryWLlWdbaXnC0i50RkV6njI0Vkv4gkisjUis5hjFltjJkMfAvMrWosjhbsE0yeySMtO82+F/iFQHYa5GXjFdFcu4NUmVqH+DPrgZ786+HexLVozHsrDjHwf5bz8NzNrNh/TgeQlUtUp5T0HOBd4KOiAyJiAaYBw4AkYLOILAQswOulXv+gMeZc4d/vAx6uRiwOVXLBWCOf8guLFSteK5CCNSKSy6tWY4xBRJwYpaqr+rUNoV/bEE5ezOSzjcf5fPNxftx7jpbBfozv3YJ7ekTS2N/L1WEqN1GdPYZXAaVHT3sBicaYw8aYHOBz4DZjzE5jzKhSt3MAItICuGSMKfNnt4hMEpEEEUlITk6uariVUunSEUVJoHCtgMnOJq+GYlV1V/NGvjw/ogPrpg7lnXFxNAnw4c+L9jHwjeUs3n3G1eEpN+HoMYHmQMkls0mFxyryEPBheQ8aY2YYY+KNMfGhoaEOCPH6iloC9ieBotIRyXhFFE0T1fIRyj5enh6M7taM+ZP78v3TA4gK8efRj7fw2nd7yM23c1xKqSpydBIoq/+jwo5OY8wfjTG1ZlAYqlA6wq+oOygVa/GCMS0foSqvU3gg8yf3ZWLflnyw+gjjZmzgzCWdVqqcx9FJIAmILHE/Ajjl4PdwuoZeDfEQj0oUkStRP6i5reGjg8Oqqrw9Lbx8Wxf+PjaWPafTuPWd1axNTHF1WKqecnQS2Ay0E5HWIuIFjAUWOvg9nM7iYaGxd2P7WwI+DcHDalsr4O2NZ1iYdgeparsttjkLp9xAY38vfj5rI/9YelBnECmHq84U0c+A9UAHEUkSkYeMMXnAFGAxsBeYb4zZ7ZhQa1alSkeI2FoDV2yDwVpNVDlK27AAvn7iBm7r1oz/XXKAX8zZzKKdp9ly7DwnzmeQlauLzlT1VHmKqDFmXDnHFwGLqhxRLVHp0hGFG84DeEVGcGXzZidFptyNv7cnb42JJb5VEC9/u4eVB66eedbIz0qTAB/CAr35WUwz7u0ZWc6ZlLpWddYJ1GvBvsGcOFeJX/P+/00C1uYR5C38BpOTg3jpfG9VfSLCz/u0ZHRsM05eyORsWhbn0rI5m5bF2XTb3w+eu8wLX+6kd1QQLYP9XR2yqiM0CZSj8kXkQuD8YQDbDCFjyD19Gq+WLZ0UoXJHgT5WAsOtdAoPvOaxc2lZDPif5by7LJE37unmguhUXaQF5MoR7BNMZl4mGbl21nbxD4UM2xiCV0ThDCEtJKdqUFigD+N7t+SLbSc5lnrF1eGoOkKTQDkqvdewXzDkXIbcTKwRuq+Aco3Jg6Lw9BDeXZZYqdclnrvMifNazM4daRIoR+VLRxStGk7BMywMsVrJPalJQNWsqrQGTpzP4PZpaxn6t5X8c0UiebpK2a1oEihHsI8tCZzPtHfB2H9LR4jFgrVZM+0OUi5R1Br4hx2tgfwCwzPztiPA4Pah/M8P+7n9n2vZfeqS8wNVtYImgXJUvYic7fm6r4BylaLWwJfbTnI0peLWwHsrEkk4doGXb+/MjInxTP95d85cyua2d9fy5uL9ug7BDWgSKEdjn8ZAJcYESpSOALBGNNcFY8pliscGlpffGthx4iJv/3iQn3Vrxu2xtskMI7uE8+OzA7k9rjnvLk/k1ndWs+VYJWbJqTpHk0A5vC3eBFgDqlBErmjBWCT5ly6Rn57upAiVKt/1WgMZOXn8ct52wgK8efW2LlftfdHIz4s37+nG3Ad7kZVbwN3T1/Onb3aTk6djBfWRJoEKVKp0hHcAWLz/2xJorjOElGtNHlx+a+CVb/dyNPUK/3tvLA39rGW+flD7UBY/M5AJfVry4dqj3D97E5cyc50dtqphmgQqUKkFY0X1g4rHBGxJQKuJKlcJC/Dh532ubQ0s2XOWzzYdZ9LAKPq2Ca7wHA28PXn5ti787d5uJBw7z13vrdOppPWMJoEKBPsG298dBFcVkdPNZVRt8GipmULn0rP4zf/9RHR4IM8Oa2/3ee7sHsFHD/bmXFoWd/xzLdtPXHRWyKqGaRKoQHWKyFkaNsQjIEAHh5VLFbUGvtp+kiMpV/jVv3/iSnYe74yLxdvTUqlz9W0TzBeP98PHamHsjPW6BWY9oUmgAsE+wVzMvkhugZ39oP6hxUkAbF1CObpgTLlYUWtgwqyNrDyQzO9u7UTbsIAqnattWABfPn4DHZoGMvmTLcxacwRjdI+DukyTQAWKSkdczLKz6esfAhn/TQJezSPI1QVjysWKWgNJFzIZ3CGUCX2qV9QwNMCbzx/pw/DoJrzy7R5eWribfN3sps7SKqIVKLlgLNTPjk3u/UMgNwNyroCXP55hYVzZuNHJUSp1fVNubIvV4sFD/VtfNR20qny9LPxzfA9eX7SXmWuOEOBj5fkRHRwQqappLm0JiEi0iMwXkfdE5G5XxlKW4iJy9paO8Lt6wRge2tBStUNjfy+m3tyR0ABvh53T4iG8OCqaUTHhzF57hPNXchx2blVzqrO95GwROSciu0odHyki+0UkUUSmXuc0NwP/MMY8BkysaizOUvUicpWYUaRUHff00HZk5uYzc/VhV4eiqqA6P1XnACNLHhARCzAN25d7NDCu8Nd+VxH5ttQtDPgYGCsibwAVT1h2gUqXky5VOkIpd9CuSQC3dA1n7rqjXNDWQJ1T5SRgjFkFlP527AUkGmMOG2NygM+B24wxO40xo0rdzhXengCmAimUQUQmiUiCiCQkJ9fsl2sDawO8PLzsXytQXESuzI+iVL311JB2XMnJZ9aaI3a/JvVyNpM+SmDBliSdYeRCju60bg6UnBifVHisTCLSSkRmAB8Bb5T1HGPMDGNMvDEmPjTUjsFZBxKRypWOKD0moJSb6NA0gFu7hjNn3VEuZly/NWCM4YUvdvKfPWd5/t87mDh7k65EdhFHJ4Gyph2Um+KNMUeNMZOMMeONMWscHItDBPkE2Z8EvPzB0/eqtQJKuYsnh7blcnYes+1oDfx7SxL/2XOWF27uyMu3dWbrsQuMeHsVc9Ye0emmNczRSSAJiCxxPwI45eD3qFHBvsH2zw4SuWbBGNrMVW6iY9NAbu7SlA/XHuVSRvkLLI+nZvCnhbvpExXEIwOimNi3FYufGUjPVkG89M0e7pm+joNntfpuTXF0EtgMtBOR1iLiBYwFFjr4PWpUpVoCAP7B/x0TqP50bKXqlKeGtiM9O49Za8tuDeQXGJ6dvx0PEd68pxseHrZ/JBGN/Zjzi5787d5uHE65wq3vrOEfSw9q+eoaUJ0pop8B64EOIpIkIg8ZY/KAKcBiYC8w3xiz2zGhukawTzDns87bP3DlH6pjAsptdQoPZETnJny49kiZZaenrzxUvJNZRGO/qx4TEe7sHsGPzw5ieOcm/O+SA9wzfR3n0rJqKny3VJ3ZQeOMMeHGGKsxJsIYM6vw+CJjTHtjTBtjzGuOC9U1gnyCyCvIIy0nzb4X+IXoOgHl1p4a2o70rDw+LNUa2HXyEm8tOcCtXcOLdzIrS0gDb969rzvvje/OwXOXuX3aWvaetvPfn6o0XdJ6HV4WLwDyCvLse0FROWkdC1BuqnOzhgyPbsLsNUdIy7K1BrJy83lm3naC/L149fYudpWuuLlrOPMf7Uu+Mdz93jqW7zvn7NDdkiYBR/MPgfxsyLns6kiUcpmnhrYjLSuPOWuPAvDXH/Zx8Nxl3rinG439vew+T5fmDfn6if60CvHnobmb+Wj9UafE6840CThacemIwnEBbREoN9SleUNu6tSEWWuO8MOu07btKfu2ZFD7yq/1adrQh/mP9mVIxzD+8PVurVrqYJoEHK14wViqQ6o1KlVXPT20HZcyc3ns0620CfVn6s2dqnwuf29P3p8Qz0P9WzNn3VEmfZTAlWw7u2hVhTQJOJrWD1IKgK4RDbmpUxgWEd4aE4uvV+V2MivN4iH8flQ0r9zehRUHkrln+nrSs3Tj++rS/QQcTesHKVXs7bFxJF3IoGPTQIedc0KfljRr6MNDcxP4fNMJHhkY5bBzuyNtCdjJlF/94mpaP0ipYg28PR2aAIoM7dSEXq2DmLv+qI4PVJMmgeuQyi779fIDq7/WD1LKyX7RrxVJFzJZuvesq0Op0zQJOIN/iCYBpZxsWHQTmjX0Yc66o5V6XXJ6Np9vOk5mTr5zAqtjNAk4Q9GCMdApoko5iafFgwl9W7HuUCr7z9hfcO6FL3Yy9YudDHtrJUv2aCtCk4Az+IcWDgzrFFGlnGlsz0i8PT3sbg2sS0zhx71nubtHBH5eFh75KIEH52zmWOoV5wZai2kScAbtDlKqRjT29+KOuOZ8uS3pupvZ5BcYXvluL80b+fLq7V347qkBvHhrJzYeTmXYW6t4a8kBsnLdr4tIk4Az+GkSUKqm3N+vFVm5BczbfKLC5y3YcoK9p9OYenNHfKwWrBYPHh4QxbLnBzOyc1P+vvQgw95a6XYDzZoErqNKq379Q6EgF/KyHR+QUuoqncID6RMVxEfrj5GXX/b+A5ez83hj8QF6tGzMqJjwqx5rEujDO+Pi+NcjvfH2tPDQ3ASm/GtrhRvj1CeaBJyhaMFYXqZr41DKTTzQrzUnL2by496yK41OX3GIlMvZ/H5UdLk/7Pq1CWHRUwN4fnh7fth1hpF/X8W6Q/W/RV9jSUBEokRklogsqOhYvVCUBHI1CShVE27qFEbzRr7MWXftjmYnL2bywerD3B7bjNjIRhWex8vTgylD2vHF4/3wtVoYP3Mjry/aS3Ze/R0rsCsJiMhsETknIrtKHR8pIvtFJFFEplZ0DmPMYWPMQ9c7Vi/4lUgCOkVUKafztHgwsW9LNhw+f80GNH/9fh8AvxrZ0e7zxUQ04tun+nNfrxa8v+owd0yrv/se29sSmAOMLHlARCzANOBmIBoYJyLRItJVRL4tdQtzaNS1XVE5ae0OUqrGjOkZiY/Vg7klpotuPX6BhTtOMWlgFM0b+VbqfH5enrx2R1dmToznbFoWo/6xhrnrjtq/1WwdYVcSMMasAs6XOtwLSCz8NZ8DfA7cZozZaYwZVermXlsC6ZiAUjWukZ8Xd8RF8OW2k1y4koMxhle+3UNYgDeTB7Wp8nlvim7C978cQN82wfxx4W4mfbyFnLyyB6DrouqMCTQHSs7JSio8ViYRCRaR6UCciLxQ3rEyXjdJRBJEJCE52XVF2SqV/T29wTtQxwSUqmEP9GtFdl4Bn28+wTc/nWbb8Ys8P6ID/t7VK5gcFuDDhw/05MVbO7Fkz1l+vWAHBfWkcF11rkxZQ+zlXhVjTCow+XrHynjdDGAGQHx8fN256n7BmgSUqmEdmgbQr00wH68/iojQuVkgd3ePcMi5RYSHB0SRnVfAG4v30yTQhxduqfpGObVFdVoCSUBkifsRwKnqhVOP+IdqElDKBR7o14pTl7I4eTGTF2+NxsPDseVbHh/chgl9WvL+qsPMXnPtbKS6pjotgc1AOxFpDZwExgL3OSSq+sA/BHJT0KUYStWsoZ2a0CbUn45NA+nbJtjh5xcRXhrdmXPpWbzy3R7CAr0ZFdPM4e9TU+ydIvoZsB7oICJJIvKQMSYPmAIsBvYC840xu50Xah3jHwK5GfZuRaOUchCLh/DtkwN4e2ysU9/j72Pj6NGiMc/O21GnF5XZOztonDEm3BhjNcZEGGNmFR5fZIxpb4xpY4x5zbmh1jF+IZCXpesElHIBXy9bbSBn8rFamHl/PC2C/Xj0oy3XrE+oK7Svwk52by9ZxD8UTP2ZRqaUulYjPy/mPtgLf29PHvhwEycv1r1xQE0C11GlAnLw37UC2iGkVL3WvJEvcx7sSUZOPvfP3sTaxJQ6VWaiepNnVfmKkoB2BylV73VsGsgHE+N5cM5mxs/ciK/VQp+oIAa2D2VQ+1Bah/hX/Qelk2kScJai0hHaElDKLfSJCmbz725i/aFUVh1MZtWBZJbvty1wjWjsy8D2oYztGUlMRMVF7GqaJgFn8dOWgFLuxt/bk5uim3BTdBMAjqdmsPJgMiv3J/P1tpN8te0k8yb1pWtEQxdH+l86JuAsfoXzk6uYBAqysijIyHBgQEqpmtYi2I8JfVoy8/54lv9qMI39vPjFnM2cOF97/m1rEnAWTy/brYrdQQf63cD+7j0cG5NSymXCAnyY+2BPcvMLuP/DTdfdE7mmaBKwU5XKx1r9qtwSMNoKUKreaRsWwIwJPUg6n8kjHyXUio3tNQlch5RZJ89OVl90YFgpVVLvqGDevLcbm49e4Ll/u74aqQ4MO5PVF0z93I1IKVV1o7s14/TFTF7/fh/NGvrwu1ujXRaLJgFnsvqhLQGlVFkmDYwq3P/4CM0b+fLADa1dEod2BzmTbyPbmMBl122Go5SqnUSEP/6sMzd1asKfvt3D4t1nXBKHtgScqWjB2JmfoO1Q18ailKp1LB7CP8bFMfaDDTz68Ra8PMv/XX5vfASv3t7V4TFoEnAmTQJKqevw9bLw4QM9+Xj9MTIrmC0UG+mcBWaaBOxU6SqiAFYfEA84s9PxASml6o0gfy+evqmdS95bxwSuo1pTRAvPwOmfHBKLUko5Wo0lARGJEpFZIrKgxLFOIjJdRBaIyGM1FUuNEg9ITYScK66ORCmlrmHv9pKzReSciOwqdXykiOwXkUQRmVrROYwxh40xD5U6ttcYMxm4F4ivbPB1gngABs7qzptKqdrH3pbAHGBkyQMiYgGmATcD0cA4EYkWka4i8m2pW1h5JxaR0cAaYGmVPkFtV1RD/Ix2CSmlah+7BoaNMatEpFWpw72ARGPMYQAR+Ry4zRjzOjDK3gCMMQuBhSLyHfCv0o+LyCRgEkCLFi3sPW2tYNtEQsC3sY4LKKVqpeqMCTQHTpS4n1R4rEwiEiwi04E4EXmh8NhgEXlHRN4HFpX1OmPMDGNMvDEmPjQ0tKyn1H5Nu2pLQClVK1VnimhZ02bKnUdpjEkFJpc6tgJYUY0YnM4hW8I1jYFNH0B+Hlh0Vq5SqvaoTksgCYgscT8COFW9cOqppjGQnw0pB1wdiVJKXaU6SWAz0E5EWouIFzAWWOiYsOqZ8Bjbn7poTClVy9g7RfQzYD3QQUSSROQhY0weMAVYDOwF5htjdB5kWYLbgaePjgsopWode2cHjSvn+CLKGdBVhYyxjQOERWsSUErVOlo2wplKDio37WqbJlrF7SaVUsoZNAnYqUp7DJcUHgNZF+FSkmMCUkopB9AkcB3VLyBXqGnR4LB2CSmlag9NAjWlSWe0oqhSqrbRJFBTvPwhuK1OE1VK1SqaBGpSeIx2BymlahVNAs5WckC5aVe4dAIyzrsuHqWUKkGTgJ2qtL1k6UHlprpyWClVu2gSqEmaBJRStYwmgZrUIBQCwnVcQClVa2gSqGlNY7QloJSqNTQJ1LSmXSF5P+RmujoSpZTSJFDjwmPA5MO5Pa6ORCmlNAk42zVzipp2tf2pXUJKqVpAk4CdqjRFtKytKRu1Au9ALR+hlKoVNAlch0P2GC7JwwOadNGWgFKqVqixJCAiUSIyS0QWlDg2WERWi8h0ERlcU7G4XHgMnN0FBfmujkQp5ebs3V5ytoicE5FdpY6PFJH9IpIoIlMrOocx5rAx5qHSh4HLgA+2jevdQ9OukJsB5w+7OhKllJuztyUwBxhZ8oCIWIBpwM1ANDBORKJFpKuIfFvqFlbOeVcbY24GfgP8qWofoQ4qWjl8eodr41BKuT179xheJSKtSh3uBSQaYw4DiMjnwG3GmNeBUXaet6DwrxcA77KeIyKTgEkALVq0sOe0tV9oR/Cw2sYFut7t6miUUm6sOmMCzYETJe4nFR4rk4gEi8h0IE5EXig8dqeIvA98DLxb1uuMMTOMMfHGmPjQ0NBqhFs9Vd5esqzXeXpBWEctH6GUcjm7WgLlKGvaTLnflMaYVGByqWNfAF9UIwanq9b2khXNLGoaAwcW25KEo2cgKaWUnarTEkgCIkvcjwBOVS8cN9I0BjJSIP2MqyNRSrmx6iSBzUA7EWktIl7AWGChY8JyA8Urh7VLSCnlOvZOEf0MWA90EJEkEXnIGJMHTAEWA3uB+caY3c4LtZ5p2hXaj7TtPayUUi5i7+ygceUcXwQscmhE7sInEO6b5+oolFJuTstGKKWUG9MkYKeq7TFM2VNElVKqltAkcB3VKiCnMz+VUrWcJgGllHJjmgSUUsqNaRJQSik3pklAKaXcmCYBe1V1ko/ODlJK1WJS5eqYLiAiycCxKr48BEhxYDj1mV4r++h1so9eJ/s48zq1NMaUWYa5TiWB6hCRBGNMvKvjqAv0WtlHr5N99DrZx1XXSbuDlFLKjWkSUEopN+ZOSWCGqwOoQ/Ra2Uevk330OtnHJdfJbcYElFJKXcudWgJKKaVK0SSglFJuzC2SgIiMFJH9IpIoIlNdHU9tIiKzReSciOwqcSxIRJaIyMHCPxu7MkZXE5FIEVkuIntFZLeIPF14XK9TKSLiIyKbRGRH4bX6U+Hx1iKysfBazSvcktbtiYhFRLaJyLeF92v8OtX7JCAiFmAacDMQDYwTkWjXRlWrzAFGljo2FVhqjGkHLC28787ygOeMMZ2APsAThf8P6XW6VjYwxBjTDYgFRopIH+CvwFuF1+oC8JALY6xNnsa2PW+RGr9O9T4JAL2ARGPMYWNMDvA5cJuLY6o1jDGrgPOlDt8GzC38+1zg9hoNqpYxxpw2xmwt/Hs6tn+0zdHrdA1jc7nwrrXwZoAhwILC43qtABGJAG4FZhbeF1xwndwhCTQHTpS4n1R4TJWviTHmNNi+AIEwF8dTa4hIKyAO2IhepzIVdnFsB84BS4BDwEVjTF7hU/TfoM3bwK+BgsL7wbjgOrlDEihrG6C6hQAAAYZJREFUfy+dF6sqTUQaAP8H/NIYk+bqeGorY0y+MSYWiMDWEu9U1tNqNqraRURGAeeMMVtKHi7jqU6/Tp7OfoNaIAmILHE/AjjloljqirMiEm6MOS0i4dh+0bk1EbFiSwCfGmO+KDys16kCxpiLIrIC2zhKIxHxLPyVq/8G4QZgtIjcAvgAgdhaBjV+ndyhJbAZaFc46u4FjAUWujim2m4hcH/h3+8HvnZhLC5X2Fc7C9hrjPlbiYf0OpUiIqEi0qjw777ATdjGUJYDdxc+ze2vlTHmBWNMhDGmFbbvpGXGmPG44Dq5xYrhwmz7NmABZhtjXnNxSLWGiHwGDMZWxvYs8EfgK2A+0AI4DtxjjCk9eOw2RKQ/sBrYyX/7b3+LbVxAr1MJIhKDbUDTgu1H5nxjzMsiEoVtUkYQsA34uTEm23WR1h4iMhh43hgzyhXXyS2SgFJKqbK5Q3eQUkqpcmgSUEopN6ZJQCml3JgmAaWUcmOaBJRSyo1pElBKKTemSUAppdzY/wNVOvCIf4UEtAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.semilogy(np.abs(np.array(x_bisect) - np.pi), label='Bisect')\n",
    "plt.semilogy(np.abs(np.array(x_secant) - np.pi), label='Secant')\n",
    "plt.semilogy(np.abs(np.array(x_newton) - np.pi), label='Newton')\n",
    "plt.semilogy(np.abs(np.array(x_brent)  - np.pi), label='Brent')\n",
    "plt.legend()\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Обратите внимание, что для метода деления пополам график ошибки в логарифмическом масштабе (по вертикальной оси) - практически прямая (**это и означает сходимость с линейной скоростью или со скоростью геометрической прогрессии!**) Для остальных методов - сходимость еще более быстрая!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Выводы\n",
    "\n",
    "1. Метод деления пополам - простой и надежный метод нахождения корня нелинейного уравнения, но не самый быстрый! Скорость сходимости у него - линейная.\n",
    "2. У метода Ньютона - скорость сходимости квадратичная, но он требует вычисления производной.\n",
    "3. Метод секущих работает быстрее метода деления пополам, но медленнее метода Ньютона.\n",
    "4. Метод Брента и разные его модификации - это надежный метод, не требующий вычисления производных. Этот метод работает не медленне метода деления пополам. Как правило, это первый метод, который нужно попробовать для решения вашей задачи"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}