M368K Workbook - Fixed Point algorithms for nonlinear systems.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Iterative methods for solving nonlinear systems of equations\n",
    "by [Richard Tsai](https://www.researchgate.net/profile/Richard_Tsai3)\n",
    "\n",
    "* Direct generalization of the **fixed-point iterations** for a single equation\n",
    "$$ x_{k+1}:=g(x_k) $$ with $x_k$ being a vector and $g(x_k)$ also. Here, \n",
    "$$g:\\mathbb{R}^d \\mapsto \\mathbb{R}^d.$$\n",
    "\n",
    "* Fixed point of $g$ corresponds to the solution of $f(x)=0$.\n",
    "* Generalization of  **Jacobi** and **Gauss-Seidel** iterations to solve nonlinear systems:\n",
    "\n",
    "* **Newton's method**. Caution: need to solve a linear system at every iteration. \n",
    "\n",
    "    Generally difficult and costly for large systems. \n",
    "    \n",
    "    Try use **quasi-Newton method**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fixed point iterations for systems of equations\n",
    "\n",
    "Learning objectives:\n",
    "* The existence of a fixed point of $g$. \n",
    "* The convergence theory of the fixed point iterations.\n",
    "* Theorem 10.6 in Burden and Faires covers both issues.\n",
    "* Know how to use Gauss-Seildel iterations to solve nonlinear systems.\n",
    "* Know Newton's method and the quasi-Newton's method, called Broyden's method.\n",
    "* Know how to derive the Jacobian matrix of a given nonlinear equation.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gauss-Seidel iterations for solving nonlinear systems"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example: solving the steady-state viscous Burger's equation\n",
    "\n",
    "$$ uu_x = \\epsilon u_{xx},\\,\\,\\,\\,\\,0<x<1, $$ \n",
    "with boundary conditions\n",
    "$$u(0)=1,\\,\\,u(1)=-1. $$\n",
    "\n",
    "Discretize on a uniform grid with step size $h$.\n",
    "\n",
    "$u(x_j) \\approx u_j $\n",
    "\n",
    "$u_x(x_j) \\approx (u_{j+1}-u_{j-1})/2h $\n",
    "\n",
    "$u_{xx}(x_j) \\approx (u_{j+1}-2u_j+u_{j-1})/h^2 $\n",
    "\n",
    "We are lead to a system involving the unknowns $u_1,u_2,\\cdots, u_{n-1}$ with the boundary conditions $u_0=1$ and $u_n=-1$:\n",
    "\n",
    "$$ u_j \\frac{u_{j+1}-u_{j-1}}{2h} = \\epsilon\\frac{u_{j+1}-2u_j+u_{j-1}}{h^2}, j=1,2,\\cdots,n-1.\\,\\,\\,\\,\\, (P1)$$\n",
    "\n",
    "The above equations can be organized into:\n",
    "$$f(u_1,u_2,\\cdots,u_{n-1})=(f_1(u_1,u_2,\\cdots,u_{n-1}), f_2(u_1,u_2,\\cdots,u_{n-1}),\\cdots, f_{n-1}(u_1,u_2,\\cdots,u_{n-1}))^\\top = 0,$$\n",
    "where\n",
    "$$ f_j(u_1,u_2,\\cdots,u_{n-1}):=u_j \\frac{u_{j+1}-u_{j-1}}{2h} - \\epsilon\\frac{u_{j+1}-2u_j+u_{j-1}}{h^2}.$$\n",
    "\n",
    "In the following, we use Gauss-Seidel iterations to solve this system.\n",
    "\n",
    "Let us first re-arrange the system into a more convenient form: for $k=1,2,\\cdotsfor $j=1,2,\\cdots,n-1$\n",
    " \n",
    "$$ \\left(2\\frac{\\epsilon}{h^2}+\\frac{u_{j+1}^{(k)}-u_{j-1}^{(k+1)}}{2h}\\right)u_j^{(k+1)} = \\epsilon\\frac{u_{j+1}^{(k)}+u_{j-1}^{(k+1)}}{h^2} $$\n",
    "\n",
    "We solve for $u_j^{(k+1)}$ by an explicit formula. \n",
    "\n",
    "** In general, you might want to use Newton's method for solve the resulting nonlinear system of equations **\n",
    "\n",
    "### Homework:\n",
    "Derive the Jacobian of the nonlinear system (P1) and solve it by Newton's method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "GaussSeidel_Burgers! (generic function with 1 method)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function GaussSeidel_Burgers!(u, h, ϵ, nitrs)\n",
    "    \n",
    "    half_h=h/2.0\n",
    "    \n",
    "    for t=1:nitrs\n",
    "    for j=2:length(u)-1\n",
    "            coeff = (2ϵ+ half_h*(u[j+1]-u[j-1]))\n",
    "            \n",
    "            u[j] = ϵ*(u[j+1]+u[j-1])/coeff\n",
    "    end\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "k= 1: 0.000420187392,\t 7.986026111337\n",
      "k= 2: 0.001799711694,\t 5.547493974707\n",
      "k= 3: 0.001214294211,\t 4.130644042113\n",
      "k= 4: 0.001201154148,\t 3.091776803341\n",
      "k="
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjsAAAGgCAYAAABMn6ZGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzs3XlcVmX+//HXuW82QbgVEUTFHRVDU9EEyxYr0rKyVaei5aeW07cpy2bKbKayKafNlmnMyaymchora9rMUSvTEtfEfV8CFcT1vkGU9f79cQRFRAG5OXDf7+fjcR7ncDj3dT73TDO8u851nctwu91uRERERLyUzeoCRERERDxJYUdERES8msKOiIiIeDWFHREREfFqCjsiIiLi1RR2RERExKsp7IiIiIhXU9gRERERr6awIyIiIl5NYUdERES8msKOiIiIeDU/qwuwQklJCXv27CE0NBTDMKwuR0RERKrA7XaTk5NDy5Ytsdmq3l/jk2Fnz549xMTEWF2GiIiI1EBGRgatW7eu8vU+GXZCQ0MB8z+ssLAwi6sRERGRqnC5XMTExJT9Ha8qnww7pY+uwsLCFHZEREQamOoOQdEAZREREfFqCjsiIiLi1RR2RERExKsp7IiIiIhXU9gRERERr6awIyIiIl5NYUdERES8msKOiIiIeDWFHREREfFqCjsiIiLi1TwadhYsWMC1115Ly5YtMQyD//73v2f9zE8//URCQgJBQUF06NCBKVOmVLhm8uTJtG/fnqCgIBISEli4cKEnyhcREREv4NGwc+TIEc4//3zefPPNKl2/Y8cOrr76agYMGMDKlSt54oknePDBB5k5c2bZNTNmzGDMmDGMHz+elStXMmDAAAYPHkx6erqnvoaIiIg0YIbb7XbXyY0Mgy+++IKhQ4dWes1jjz3GV199xYYNG8rOjR49mlWrVpGamgpAv3796N27N2+99VbZNXFxcQwdOpSJEyeett38/Hzy8/PLfi5dNdXpdNbqQqAbP5/N4fc+wlyf7PgiZYb53UuXLDN/ZZQenPj55NPHfzh5oTPj+C9P/N5W2gRGaROGgYGBYYDNMDBsBjbDwHa8Btvxnw0D7Daj3LkT523Y7QZ2w8DPZpjX2U7JxKdbgO3Uc2f7ubauqU+12O3g51d+8/eveO7U3wcHV9yCguDU/9xFRHycy+XC4XBU++93vVr1PDU1leTk5HLnrrrqKqZNm0ZhYSFut5sVK1bw+OOPl7smOTmZRYsWVdruxIkTeeaZZzxS88mcy1aS+M10j99HfESjRmbwCQmBpk2hWTNzCw8/cRwZCa1bQ0yMuQ8MtLpqEZF6p16FnaysLKKiosqdi4qKoqioiP379+N2uykuLj7tNVlZWZW2O27cOB555JGyn0t7dmpbkwv7kTr8PnCDG3DjxnC7cUPZOUpKf2MqPVf2sxsMSnCXXu8uPX+iA87tdoP7xGdKP4fbjfv4+ZIS87jE7abEzfFjcJeY59wnnS+7rgRK3G6KS0rglP4+49QTp2Gc0kl4us/YgEA/G0F+dgL9bQTZbQT62wgO8CM4wE5woJ0QfzvBAXb87baTvtwpqnKurq8pKYGiInMrLDxxfOp28u8KCuDoUcjLM7eTeiA5etTcDhyAqj6mjYqC2Fjo1u3EFh8P0dFV+7yIiBeqV2EHKPfoBk78kTcMo9zxqdeceu5kgYGBBNbBv/F2GTIQhgz0+H08ze12U1jsJr+omGOFJeQXFZNfVMKxQnOfX1jCsaJijhUUk5tfxJH8Io4UFJNz7PhxfpF5vqAI19EiDuUVcOhIAUcKiqtVR2igH62aNqJdsxDaRYTQIcLct48IIaJxwBn/O2+wiovLh5+8PMjNhYMHzdBTui893rsXMjLMMHTsmPnz3r3w88/l223fHi666MQWF3f6R3UiIl6oXoWdFi1aVOihyc7Oxs/Pj2bNmuF2u7Hb7ae95tTeHqk5wzAI8DMI8LMRGlR77R4rLOZwXmFZ+Dl4fL8vt4C9zmNkuo6R5TxKpvMYOceKyMkvYmNWDhuzciq01TjQj06RjYmLDqNbdChdo8Po2iKU0CD/2ivYCnY7NG5sbtXhdpsB6LffYNMmWL/+xLZlC+zYYW4ffmhe36oVXHst3HgjDBxo3ldExEvVq7CTlJTE119/Xe7cnDlz6NOnD/7+5h+xhIQE5s6dyw033FB2zdy5c7n++uvrtFapviB/Oy0cdlo4zp6gcvOLyHIeJePgUXbsP8LOA0fYsd/cdh8+Sm5+EWkZh0nLOFz2GcOA2MjG9IppSu+2TejVpimdmjfGZvOBHgzDgIgIc0tIKP87lwsWLzZ7e37+2TzevRumTDG3li3hjjtg9GizB0hExMt4dDZWbm4uW7duBaBXr15MmjSJyy67jPDwcNq0acO4cePYvXs3H3zwAWBOPY+Pj+e+++5j1KhRpKamMnr0aD7++GNuuukmwJx6npKSwpQpU0hKSuLtt99m6tSprFu3jrZt21aprpqO5pb6Ib+omPQDeWzMymFDpuv4lkOW61iFa5sG+5PUsRn9O0ZwUacI2jYL9s7HX9Vx7Bj8+CP897/w2Wfm4zAwZ39dfz088QT06WNtjSIip1HTv98eDTvz58/nsssuq3D+rrvu4v333+fuu+9m586dzJ8/v+x3P/30Ew8//DDr1q2jZcuWPPbYY4wePbrc5ydPnsyLL75IZmYm8fHxvPrqq1x88cVVrkthxzvty8knLeMwK9MP8Wv6IVbvcpJ3yjihVk0acUVcJFed14K+7cNPDIL2VQUFMGsWvPUWzJlz4vwNN8Bzz5lje0RE6ol6GXbqK4Ud31BYXMLqXYf5ZesBft66n5XphygsPvGPu6ORP5d3jeTq7tFc0qW5gs/69fC3v8FHH5ljgPz94Y9/hPHjzSnwIiIWU9ipBoUd35RXUMSirQeYsz6LeRuyOXikoOx3EY0DuL5nK27q3ZpuLX38n4n16+Gxx+Cbb8yfO3WCjz/Woy0RsZzCTjUo7EhxiZvlOw8ye10WX6/aw/7cE8GnW3QYdya1ZWivVgT5++gsJbcbvvwS/vAH2LXL7OV54QUYM0ZT1kXEMgo71aCwIycrLC5hweZ9zPx1F/PWZ1NQXAJAeEgAt/drQ0piWyLDanEOfkNy6BCMHAmff27+fNdd8PbbEBBgbV0i4pMUdqpBYUcqczivgM9W7OK9X3ay+/BRAPztBkN7tuLBy2OJCffBsStuN0yeDA89ZL708PLLYeZMcDisrkxEfIzCTjUo7MjZFBWXMGf9Xt79eQfLfzsEgJ/NYFjfGB4Y2IloRyOLK7TAd9/Brbeab3ROTDRnb4WGWl2ViPgQhZ1qUNiR6ljx2yFem7eZhVv2AxDgZ+OOfm35w8BONA3xscc5v/4KV15pvptnwAAzAIWEWF2ViPgIhZ1qUNiRmliy/QCvzNnM0p3mS/iaBvvz2KCu3Nonxjfe0lxqxQpziQmXCwYNMmdtabkJEakDNf377eMvFhGpun4dmjHjvkQ++H8X0CUqlEN5hTz++RpueGsRa3Y5rS6v7iQkwOzZ5rt3Zs82p6mLiNRj6tlRz47UQGFxCR+k/sarczeTm1+EYUBKYlseH9yV4IB6teSc53z6qTmGB+D9982ZWiIiHqSeHZE65G+3MeKi9vww9hKG9myJ2w0fpP7GkDd+Lrc4qVe75Rb4y1/M4/vug7Vrra1HRKQS6tlRz47UgoVb9vHHT1eT5TqG3Wbw4MBY/u+yjvh5+xIUJSVw7bXm+lrnnw9LlkBgoNVViYiXUs+OiIUGxDZn9pgBDOkRTXGJm1fnbeamKalkHMyzujTPstlg2jSIiIBVq+Cpp6yuSESkAoUdkVrSJDiAv/+uF68P70lokB+rMg5z7Zs/88vW/VaX5lktWsDUqebxiy9Caqq19YiInEJhR6QWGYbB9T1bMXvMxfRo7eBwXiEp05bwzsLtePUT46FDzQHKbjc88ID5pmURkXpCYUfEA1o1acQn9yVxU+/WlLjhr99uYMyMNI4WeHEIePFFcwmJX3+Fd96xuhoRkTIKOyIeEuRv5+VbevD0td2w2wy+TNvDrf9MZX9uvtWleUZkJDz7rHn8xBNw4IC19YiIHKewI+JBhmFw94XtmT6yH+EhAazZ7eTmtxZ578Dl3/8eunc3l5OYMMHqakREAIUdkTqR2KEZM3/fn9ZNG7HzQB43vbWIjVkuq8uqfX5+8Oqr5vE//wm7dllbj4gICjsidaZ9RAgzf9+fLlGhZOfkc+uUVJYdX2fLqwwcCJdcAvn5MHGi1dWIiCjsiNSlqLAgPrkviT5tm+I6VsQd7yxhweZ9VpdVuwwDnnnGPJ46FdLTra1HRHyewo5IHXME+/PhiH4M7BpJflEJoz5YTuo2LxvMe8klZg9PYSH89a9WVyMiPk5hR8QCjQLsTLkjgcuPB54R/1rGcm97pFXau/Ovf0F2trW1iIhPU9gRsUiAn41/3N6bAbER5BUUc897y1i9y4sWEb3oIujXDwoKYMoUq6sRER+msCNioSB/O2+n9OGC9uHk5BeRMm0p6/d40Sythx4y95MnmwOWRUQsoLAjYrFGAXbevbsvvds0wXm0kDvfXeI97+G5+WZo2RL27oVPPrG6GhHxUQo7IvVA40A/3rvnArpFh7E/t4B73l+GM6/Q6rLOnb+/uVYWmO/f8eb1wUSk3lLYEaknHI38effuvrQIC2Jrdi6jP1pBQVGJ1WWdu3vvhaAgWLkSli+3uhoR8UEKOyL1SAtHEO/e3ZeQADup2w/w+OerG/5q6c2awY03msfvvWdtLSLikxR2ROqZbi3D+MftvbHbDD7/dTdvfL/V6pLO3T33mPuPP4Zjx6ytRUR8jsKOSD10aZdInr0+HoBX523mq1V7LK7oHA0cCDExcPgwfPml1dWIiI9R2BGpp27r14b7Lu4AwGOfrWZTVo7FFZ0Dmw3uuss81qMsEaljdRJ2Jk+eTPv27QkKCiIhIYGFCxdWeu2ll16KYRgVtmuuuabsmrvvvrvC7xMTE+viq4jUqT8N6spFnSI4WljM6I9W4DrWgGdo3X23uZ8zR6uhi0id8njYmTFjBmPGjGH8+PGsXLmSAQMGMHjwYNIrWRzw888/JzMzs2xbu3YtdrudW265pdx1gwYNKnfdrFmzPP1VROqc3Wbwxu960apJI3bsP8LYT1ZRUtJAByx37AgDBpjTz2fMsLoaEfEhHg87kyZNYsSIEYwcOZK4uDhee+01YmJieOutt057fXh4OC1atCjb5s6dS3BwcIWwExgYWO668PDwSmvIz8/H5XKV20QaivCQACbf3psAu4256/cyZcE2q0uquWHDzP1nn1lbh4j4FI+GnYKCAlasWEFycnK588nJySxatKhKbUybNo3hw4cTEhJS7vz8+fOJjIykc+fOjBo1iuwzLDQ4ceJEHA5H2RYTE1P9LyNiofNjmvDM9ecB8PL/NvHzlv0WV1RDN9wAhgGLF0NGhtXViIiP8GjY2b9/P8XFxURFRZU7HxUVRVZW1lk/v3TpUtauXcvIkSPLnR88eDDTp0/nhx9+4JVXXmHZsmUMHDiQ/ErW3hk3bhxOp7Nsy9D/yUoDNLxvDLf2aU2JGx78z0qyXQ1wCnfLlnDhhebx559bW4uI+Iw6GaBsGEa5n91ud4VzpzNt2jTi4+O54IILyp0fNmwY11xzDfHx8Vx77bV89913bN68mW+//fa07QQGBhIWFlZuE2loDMNgwvXxxEWHcfBIAWM/baDjd26+2dzrUZaI1BGPhp2IiAjsdnuFXpzs7OwKvT2nysvL4z//+U+FXp3TiY6Opm3btmzZsuWc6hWp74L87bwxvCeBfjYWbtnPu7/ssLqk6it9m/Ivv8CeBv7+IBFpEDwadgICAkhISGDu3Lnlzs+dO5f+/fuf8bOffPIJ+fn53HHHHWe9z4EDB8jIyCA6Ovqc6hVpCGKjQnlySDcAXpy9iXV7nBZXVE0xMZCYaM7K0qMsEakDHn+M9cgjj/DOO+/w7rvvsmHDBh5++GHS09MZPXo0AHfeeSfjxo2r8Llp06YxdOhQmjVrVu58bm4ujz76KKmpqezcuZP58+dz7bXXEhERwQ033ODpryNSL9zRrw1XxEVRUFzCQ/9J42hBsdUlVc9NN5n7r7+2tg4R8QkeDzvDhg3jtddeY8KECfTs2ZMFCxYwa9Ys2rZtC0B6ejqZmZnlPrN582Z+/vlnRowYUaE9u93OmjVruP766+ncuTN33XUXnTt3JjU1ldDQUE9/HZF6wTAMXripO81DA9manctzs9ZbXVL1lL4k9Kef4MgRa2sREa9nuBv8ksrV53K5cDgcOJ1ODVaWBm3B5n3c+e5SAN65sw9XdDvzWLh6w+2GDh1g506zd2fIEKsrEpEGoKZ/v7U2lkgDdnHn5oy4qD0AT3yxBmdeA1lOwjBg8GDzWG8/FxEPU9gRaeD+eFUXOkSEkJ2Tz4RvGtDjrKuvNvezZpk9PSIiHqKwI9LABfnbeemWHhgGzPx1Fz9s3Gt1SVVz2WUQGAi//QYbN1pdjYh4MYUdES+Q0DacERcef5z1+VqcRxvA46yQELj0UvNYj7JExIMUdkS8xNjkLrRrFkyW6xjPfdtAHmed/ChLRMRDFHZEvESjADsv3nw+hgGfLN/FT5v3WV3S2V11lbn/+Wc4etTaWkTEaynsiHiRC9qHc3f/dgCMm7ma3Pwiaws6m86dIToaCgogNdXqakTESynsiHiZP17VhZjwRuxxHuPVuZutLufMDAMGDjSPf/zR2lpExGsp7Ih4meAAP569Ph6A937Zwdrd9XztrMsuM/cKOyLiIQo7Il7o0i6RDOkRTYnbfNlgcUk9fo9NadhZskRLR4iIRyjsiHipvwzpRmiQH6t3OfkwdafV5VSufXto0waKisyByiIitUxhR8RLRYYF8adBXQF4ec5mspzHLK6oEoahR1ki4lEKOyJe7PYL2tAzpgm5+UU88/U6q8upnMKOiHiQwo6IF7PZDCbe2B27zeC7tVn1dymJ0rCzfDm4XNbWIiJeR2FHxMvFRYeVrYz+9FfrOVZYbHFFp9GmDbRrByUl5kBlEZFapLAj4gMevDyWyNBA0g/mMe3nHVaXc3oXXmjuFy2ytg4R8ToKOyI+oHGgH09cHQfAmz9sZc/herg0Q//+5l5hR0RqmcKOiI+4vmdL+rZrytHCYp6ftcHqcioqDTupqVBcDx+1iUiDpbAj4iMMw+Dp687DZsA3qzNJ3XbA6pLKi4+Hxo0hJwfW1eOZYyLS4CjsiPiQ81o6uK1fGwCe+XodRcUlFld0Ej8/SEw0j/UoS0RqkcKOiI8Ze2UXmgT7szErh48W/2Z1OeVp3I6IeIDCjoiPaRoSwKPJXQCYNHczB3LzLa7oJAo7IuIBCjsiPuh3F7ThvJZhuI4V8dL/NlldzgmJiebyEdu2wd56+gJEEWlwFHZEfJDdZvDMdecBMGN5Bqt3Hba4ouMcDnOgMqh3R0RqjcKOiI/q0y6cG3q1wu2Gp75aR0mJ2+qSTP36mftly6ytQ0S8hsKOiA97fHBXQgLsrEw/zBcrd1tdjqlvX3O/fLm1dYiI11DYEfFhUWFBPDAwFoAX/7eRI/lFFlcE9Olj7pcvB3c96W0SkQZNYUfEx/2/i9rRJjyYva58pvy0zepyzDE7gYFw6JA5UFlE5Bwp7Ij4uEA/O09c3RWAtxdsZ9ehPGsLCgiA8883j/UoS0RqgcKOiHDVeS1I7BBOflEJf/tuo9XlnBi3o0HKIlILFHZEBMMw+MuQE+tmLdt50NqCTh63IyJyjuok7EyePJn27dsTFBREQkICCxcurPTa999/H8MwKmzHjh2rcZsicnbdWoYxrK+5btaEr9dbOxW9tGdnxQqtgC4i58zjYWfGjBmMGTOG8ePHs3LlSgYMGMDgwYNJT0+v9DNhYWFkZmaW24KCgs6pTRE5u7HJnQkN9GPNbief/brLukK6doWQEDhyBDbWg8dqItKgeTzsTJo0iREjRjBy5Eji4uJ47bXXiImJ4a233qr0M4Zh0KJFi3LbubSZn5+Py+Uqt4lIRRGNA/nD5Z0AeOl/m8i1aiq63Q69e5vHepQlIufIo2GnoKCAFStWkJycXO58cnIyi87wKvjc3Fzatm1L69atGTJkCCtXrjynNidOnIjD4SjbYmJizuFbiXi3u/u3p12zYPbl5DP5x63WFVI6bkeDlEXkHHk07Ozfv5/i4mKioqLKnY+KiiIrK+u0n+natSvvv/8+X331FR9//DFBQUFceOGFbNmypcZtjhs3DqfTWbZlZGTUwrcT8U4BfjaeuDoOgHd+3kHGQYumoickmPu0NGvuLyJeo04GKBuGUe5nt9td4VypxMRE7rjjDs4//3wGDBjAJ598QufOnfn73/9e4zYDAwMJCwsrt4lI5a7sFsWFnZpRUFTCxO82WFNEz57mftUqKCmxpgYR8QoeDTsRERHY7fYKPS7Z2dkVemYqY7PZ6Nu3b1nPTm20KSJnZhgGfx7SDZsBs9ZksXj7gbovoksX803KubmwfXvd319EvIZHw05AQAAJCQnMnTu33Pm5c+fSv3//KrXhdrtJS0sjOjq61toUkbPr2iKM2/qdmIpeXNdT0f38oHt381iPskTkHHj8MdYjjzzCO++8w7vvvsuGDRt4+OGHSU9PZ/To0QDceeedjBs3ruz6Z555hv/9739s376dtLQ0RowYQVpaWtn1VWlTRGrHw1d0JjTIj/WZLj5dbsFYt9JHWQo7InIO/Dx9g2HDhnHgwAEmTJhAZmYm8fHxzJo1i7Zt2wKQnp6OzXYicx0+fJh7772XrKwsHA4HvXr1YsGCBVxwwQVVblNEakezxoE8dHksf/12Ay/P2cQ1PaIJDfKvuwIUdkSkFhhut9vC16Raw+Vy4XA4cDqdGqwschYFRSUMem0B2/cf4b5LOjBucFzd3fyXX+Cii6BVK9hl4UsORaReqOnfb62NJSJnFOBnY/w1ZsB57+ed/HbgSN3dvEcPc797N+zbV3f3FRGvorAjImc1sGskA2IjKCgu4blv63AqemgodDLf6MyqVXV3XxHxKgo7InJWpVPR7TaDOev3smjr/rq7ucbtiMg5UtgRkSrpHBXK7aVT0b+pw6noCjsico4UdkSkyh6+ojOORv5szMrhP8vS6+ampWHnpDXyRESqQ2FHRKqsaUgAY66IBeCVOZtxHSv0/E1LBylv3gz5+Z6/n4h4HYUdEamWOxLb0rF5CAePFPD377d4/oatW4PDAUVFZuAREakmhR0RqRZ/u40nh3QD4P1FO9mx38NT0Q0D4uPN4zVrPHsvEfFKCjsiUm2XdYnk0i7NKSx289y36z1/w9Kws3at5+8lIl5HYUdEauTJa8yp6PM2ZLNwi4df+Fe6IKjCjojUgMKOiNRIp8jGpCSa69E9+816iopLPHczPcYSkXOgsCMiNTbmiliaBPuzeW8uHy/14FT00rCzcyfk5HjuPiLilRR2RKTGmgQH8MiVnQGYNHczzjwPTUVv1gyio83j9XUwRkhEvIrCjoick9suaENsZGMO5RXyuienoutRlojUkMKOiJwTP7uNPx+fiv5B6k62Zud65kYapCwiNaSwIyLn7OLOzbm8ayRFJR6ciq7p5yJSQwo7IlIrxl8Th5/N4MdN+5i/Kbv2b6DHWCJSQwo7IlIrOjRvzF392wHw1283UFjbU9G7dTPfppydbW4iIlWksCMitebBy2MJDwlga3YuHy3+rXYbDwmBdu3M440ba7dtEfFqCjsiUmscjfwZm3x8KvqczWS7jtXuDeLizP2GDbXbroh4NYUdEalVw/u24fzWDnLyi3huVi2Hkq5dzb3CjohUg8KOiNQqu83gr0O7YxjwZdoeFm3dX3uNl/bs6DGWiFSDwo6I1LrurR1l62Y9+eVaCopqabCyHmOJSA0o7IiIR4xN7kJE40C27zvC1IXba6fR0sdY6emQ66GXF4qI11HYERGPcDTy58lrzJ6Yv/+whYyDeefeaLNm0Ly5ebx587m3JyI+QWFHRDzm+p4tSewQzrHCEp75el3tNKpHWSJSTQo7IuIxhmHw16Hx+NsN5m3IZvbazHNvVDOyRKSaFHZExKM6RYYy+pKOAPz5y3U4jxaeW4OakSUi1aSwIyIe93+XdaJD8xD25eTzt+/OsUdGj7FEpJoUdkTE44L87fztxh4AfLw0g8XbD9S8sdLHWFu2QFFRLVQnIt5OYUdE6sQF7cO5rV8bAJ74fA3HCotr1lBMDAQHQ2EhbK+lKe0i4tXqJOxMnjyZ9u3bExQUREJCAgsXLqz02qlTpzJgwACaNm1K06ZNueKKK1i6dGm5a+6++24Mwyi3JSYmevpriMg5enxwVyJDA9m+/whv/rC1Zo3YbBqkLCLV4vGwM2PGDMaMGcP48eNZuXIlAwYMYPDgwaSnp5/2+vnz5/O73/2OH3/8kdTUVNq0aUNycjK7d+8ud92gQYPIzMws22bNmuXpryIi5ygsyJ8J158HwJSftrEh01WzhhR2RKQaPB52Jk2axIgRIxg5ciRxcXG89tprxMTE8NZbb532+unTp3P//ffTs2dPunbtytSpUykpKeH7778vd11gYCAtWrQo28LDwyutIT8/H5fLVW4TEWsMio/mqvOiKCpx88fPVlFYXIOlJDQjS0SqwaNhp6CggBUrVpCcnFzufHJyMosWLapSG3l5eRQWFlYIM/PnzycyMpLOnTszatQosrOzK21j4sSJOByOsi0mJqb6X0ZEas2z18fjaOTP2t0u3pq/rfoNaEaWiFSDR8PO/v37KS4uJioqqtz5qKgosrKyqtTG448/TqtWrbjiiivKzg0ePJjp06fzww8/8Morr7Bs2TIGDhxIfn7+adsYN24cTqezbMvIyKj5lxKRcxYZFlT2OOuN77ewbo+zeg2c/BjL7a7l6kTE2/jVxU0Mwyj3s9vtrnDudF588UU+/vhj5s+fT1BQUNn5YcOGlR3Hx8fTp08f2rZty7fffsuNN95YoZ3AwEACAwPP4RuISG277vyWzFqTyf/W7WXsJ6v46oGLCPCr4r8ny0IFAAAgAElEQVR/xcaC3Q45ObBnD7Rq5dliRaRB82jPTkREBHa7vUIvTnZ2doXenlO9/PLLPP/888yZM4cePXqc8dro6Gjatm3Lli1bzrlmEakbhmHw3A3dCQ8JYGNWDm/+UI3//QYEQEfzrcwatyMiZ+PRsBMQEEBCQgJz584td37u3Ln079+/0s+99NJLPPvss8yePZs+ffqc9T4HDhwgIyOD6Ojoc65ZROpORONAnr0+HoB/zN/G6l2Hq/5hzcgSkSry+GysRx55hHfeeYd3332XDRs28PDDD5Oens7o0aMBuPPOOxk3blzZ9S+++CJPPvkk7777Lu3atSMrK4usrCxyc3MByM3N5dFHHyU1NZWdO3cyf/58rr32WiIiIrjhhhs8/XVEpJZd0yOaIT2iKS5xM/aTVVV/2aAGKYtIFXl8zM6wYcM4cOAAEyZMIDMzk/j4eGbNmkXbtm0BSE9Px2Y7kbkmT55MQUEBN998c7l2nnrqKZ5++mnsdjtr1qzhgw8+4PDhw0RHR3PZZZcxY8YMQkNDPf11RMQDJlwfz+LtB9mSncsLszfy1LXnnf1Dmn4uIlVkuN2+N5XB5XLhcDhwOp2EhYVZXY6IAD9uyuae95YB8P49fbm0S+SZP7BkCSQmQnS0OUhZRLxeTf9+a20sEakXLusSyd392wHw6Ker2Z97+ldJlCkds5OZCc5qTl0XEZ+isCMi9cbjg7vSJSqU/bn5PPbZas7Y8exwQMuW5rEeZYnIGSjsiEi9EeRv5/Xf9STAz8b3G7P5aMnp19AroxlZIlIFCjsiUq90bRHG44PMEPPXb9azZW9O5RdrRpaIVIHCjojUO3f3b8fFnZuTX1TCA/9eydGCSqajl/bs6DGWiJyBwo6I1Ds2m8HLt/QgonEgm/bm8PRX605/YZcu5l5vTxeRM1DYEZF6KTI0iDeG98RmwIzlGXz+666KF3XubO63boWiorotUEQaDIUdEam3+neK4KHLzUAz/ou1FcfvxMRAUBAUFsJvv1lQoYg0BAo7IlKvPTCwExd2asbRwmLun/4reQUn9eDYbOYK6ACbNllToIjUewo7IlKv2W0Grw3rRfPQQLZk5/KXL08Zv1M6bmfz5rovTkQaBIUdEan3mocG8vrx8TufrdjFx0tPev9O6bgd9eyISCUUdkSkQejfMYKxyWYvzlNfriMt47D5i9Kwo54dEamEwo6INBj3X9qRq86LoqC4hN9/tMJcP0uPsUTkLBR2RKTBMAyDl285nw7NQ8h0HuOBf/9KUcdO5i937YIjR6wtUETqJYUdEWlQQoP8eTslgZAAO4u3H+SFpdkQEWH+Ui8XFJHTUNgRkQanU2QoL99yPgBTF+7gYKt25i80SFlETkNhR0QapMHdoxl9SUcAfnQ3MU9q3I6InIbCjog0WH+8qguXdmnO1qatADi2Tqufi0hFCjsi0mDZbQZv/K4XR9uZPTzpi9M4VljJCuki4rMUdkSkQQsL8mfUyEEAtMhKZ9zM1bjdbourEpH6RGFHRBq8Vn264zYMwvKPsPCXdfxzwXarSxKRekRhR0QavqAgjLZtAWh/cDcvzN7IvPV7LS5KROoLhR0R8Q7H36R8m+Mobjc89J+VbN6bY3FRIlIfKOyIiHc4vkbWdcFHSOwQzpGCYkb+azkHjxRYXJiIWE1hR0S8w/GeHfvWLbx1ewJtwoNJP5jH6I9WUFBUYnFxImIlhR0R8Q6lq59v2kTTkACm3dWH0EA/lu44yJP/XaMZWiI+TGFHRLxD6ern27ZBURGxUaG8cVsvbAZ8snwXUxdqhpaIr1LYERHv0Lo1BAVBYSHs3AnAZV0iefKabgBM/E4ztER8lcKOiHgHmw1iY83jk9bIuufCdtzWr03ZDK0NmS6LChQRqyjsiIj3KH2UdVLYMQyDZ647j/4dm5XN0NqXk29RgSJiBYUdEfEeJw1SPpm/3cbk23vTPiKE3YePcu+Hy7WGlogPqZOwM3nyZNq3b09QUBAJCQksXLjwjNfPnDmTbt26ERgYSLdu3fjiiy/K/d7tdvP000/TsmVLGjVqxKWXXsq6des8+RVEpCE4Tc9OqSbB5gwtRyN/VqYf5jGtoSXiMzwedmbMmMGYMWMYP348K1euZMCAAQwePJj09PTTXp+amsqwYcNISUlh1apVpKSkcOutt7JkyZKya1588UUmTZrEm2++ybJly2jRogVXXnklOTl6W6qIT6ukZ6dUh+aNeev23vjZDL5M28ObP2ytw+JExCqG28P/atOvXz969+7NW2+9VXYuLi6OoUOHMnHixArXDxs2DJfLxXfffVd2btCgQTRt2pSPP/4Yt9tNy5YtGTNmDI899hgA+fn5REVF8cILL3DfffdVaDM/P5/8/BPP6F0uFzExMTidTsLCwmrz64qIlQ4ehGbNzOOcHGjc+LSXTV/yG+O/WAvAP27rzTU9ouuqQhE5By6XC4fDUe2/3x7t2SkoKGDFihUkJyeXO5+cnMyiRYtO+5nU1NQK11911VVl1+/YsYOsrKxy1wQGBnLJJZdU2ubEiRNxOBxlW0xMzLl8LRGpr8LDISLCPN5aea/N7f3acs+F7QAY+2kaq3cdroPiRMQqHg07+/fvp7i4mKioqHLno6KiyMrKOu1nsrKyznh96b46bY4bNw6n01m2ZWRk1Oj7iEgDcJZHWaXGXx3HpV2ac6ywhJH/Wk6W81gdFCciVqiTAcqGYZT72e12VzhX3eur02ZgYCBhYWHlNhHxUmcYpHwyP7uNN37Xi9jIxmTn5DPyg2XkFRTVQYEiUtc8GnYiIiKw2+0Velyys7Mr9MyUatGixRmvb9GiBUC12hQRH1LFnh2AsCB/pt3Vl/CQANbudvHIjFWUlGiGloi38WjYCQgIICEhgblz55Y7P3fuXPr373/azyQlJVW4fs6cOWXXt2/fnhYtWpS7pqCggJ9++qnSNkXEh5SGnbP07JRq0yyYKXck4G83mL0ui0lzq/Y5EWk4/Dx9g0ceeYSUlBT69OlDUlISb7/9Nunp6YwePRqAO++8k1atWpXNzHrooYe4+OKLeeGFF7j++uv58ssvmTdvHj///DNgPr4aM2YMzz//PLGxscTGxvL8888THBzMbbfd5umvIyL1XeljrE2bwO2GMzwyL3VB+3Am3tiDRz9dxZs/bqVjZAg39Grt4UJFpK54POwMGzaMAwcOMGHCBDIzM4mPj2fWrFm0bdsWgPT0dGy2Ex1M/fv35z//+Q9PPvkkf/7zn+nYsSMzZsygX79+Zdf86U9/4ujRo9x///0cOnSIfv36MWfOHEJDQz39dUSkvuvY0Qw4LhdkZ0MVH2/fnNCardm5TPlpG4/NXEPH5o3p0bqJh4sVkbrg8ffs1Ec1nacvIg1Ehw6wYwcsWAADBlT5YyUlbkZ9sJzvN2bT0hHE13+4iGaNAz1YqIhUR718z46IiCWqMUj5ZDabwaRhPWkfEcIe5zH+79+/UlRc4oECRaQuKeyIiPep4vTz03E08uftlARCAuws3n6Qv323sZaLE5G6prAjIt6nhj07pWKjQnnl1vMBeOfnHXyZtru2KhMRCyjsiIj3qeb089MZFB/N/13WEYDHZq5mQ6arNioTEQso7IiI9yl9jLVtGxTV/K3Ij1zZhUs6m0tK3D/9V3KOFdZSgSJSlxR2RMT7tG4NjRpBYSHs3FnjZuw2g9eG9aSlI4gd+48w7vM1+OAEVpEGT2FHRLyPzQaxsebxOTzKAmgaEsDfb+uNn83gm9WZfLQkvRYKFJG6pLAjIt7pHAcpnyyhbVMeG9QVgGe/Xs/a3c5zblNE6o7Cjoh4p1oYpHyykQPac0VcFAXF5vgdl8bviDQYCjsi4p1OXiOrFhiGwSu3nE+rJo1IP5jHY5+t1vgdkQZCYUdEvFMt9+wAOIL9+cftvfG3G3y3NotPlmfUWtsi4jkKOyLinUrDzu7dkJtba832jGnCI1eavUbPfL2eHfuP1FrbIuIZCjsi4p3CwyEiwjzesqVWm7734g4kdggnr6CYMTPSKNT6WSL1msKOiHgvDzzKAvP9O5Nu7UlYkB+rMg7zxve1G6ZEpHYp7IiI96rlQcona9mkEc/f2B2Af/y4lWU7D9b6PUSkdijsiIj38lDPTqkhPVpyY+9WlLjh4Rlpmo4uUk8p7IiI9yrt2fFQ2AF45rrziAlvxK5DR3numw0eu4+I1JzCjoh4r5Pfouyhd+KEBvkz6daeGAbMWJ7BT5v3eeQ+IlJzCjsi4r06dQLDAJcLsrM9dpu+7cK5K6kdAONmrtbq6CL1jMKOiHivwEBo18489sAg5ZP9aVAX2oQHs8d5jOdnbfTovUSkehR2RMS7eXiQcqngAD9euKkHAB8vTefnLfs9ej8RqTqFHRHxbh6cfn6qpI7NuDOpLQCPzVxNbn6Rx+8pImensCMi3q2OenZKPTaoK62bNmL34aO88J0eZ4nUBwo7IuLd6mD6+clCAv148fjjrA8X/8aK3/SyQRGrKeyIiHcr7dnZtg2K6uaxUv9OEdyS0BqA8V+s1dpZIhZT2BER79a6NTRqBIWFsHNnnd123NVxNAn2Z2NWDu//Unf3FZGKFHZExLvZbBAbax7XwSDlUuEhATwxOA6ASXM3s/vw0Tq7t4iUp7AjIt6vjgcpl7o5oTV92zXlaGExT3+1rk7vLSInKOyIiPer40HKpWw2g+du6I6fzWDu+r3MWZdVp/cXEZPCjoh4v5PXyKrrW0eFMuriDgA8/dU6jujdOyJ1TmFHRLyfRY+xSj04MJbWTRuxx3mMN3/cakkNIr7Mo2Hn0KFDpKSk4HA4cDgcpKSkcPjw4UqvP3jwIH/4wx/o0qULwcHBtGnThgcffBCn01nuOsMwKmxTpkzx5FcRkYasNOzs3g25uXV++0YBdp669jwApi3cwW8HjtR5DSK+zKNh57bbbiMtLY3Zs2cze/Zs0tLSSElJqfT6PXv2sGfPHl5++WXWrFnD+++/z+zZsxkxYkSFa9977z0yMzPLtrvuusuTX0VEGrLwcIiIMI+3bLGkhCviIhkQG0FBcQnPfbvBkhpEfJWfpxresGEDs2fPZvHixfTr1w+AqVOnkpSUxKZNm+hSOmDwJPHx8cycObPs544dO/Lcc89xxx13UFRUhJ/fiXKbNGlCixYtqlRLfn4++fn5ZT+7XK6afi0Raai6dIH9+81xO7161fntDcPgL0O6Mej1hcxZv5eft+znotiIOq9DxBd5rGcnNTUVh8NRFnQAEhMTcTgcLFq0qMrtOJ1OwsLCygUdgAceeICIiAj69u3LlClTKCmp/A2lEydOLHuU5nA4iImJqf4XEpGGzeJxOwCxUaGkJJoLhT7z9TqK9GZlkTrhsbCTlZVFZGRkhfORkZFkZVVt+uWBAwd49tlnue+++8qdf/bZZ/n000+ZN28ew4cPZ+zYsTz//POVtjNu3DicTmfZlpGRUb0vIyINn0XTz0/18BWdaRrsz5bsXKYvSbe0FhFfUe2w8/TTT592gPDJ2/LlywGz2/ZUbrf7tOdP5XK5uOaaa+jWrRtPPfVUud89+eSTJCUl0bNnT8aOHcuECRN46aWXKm0rMDCQsLCwcpuI+BgLp5+fzBHsz9hkM3hNmruZQ0cKLK1HxBdUe8zOAw88wPDhw894Tbt27Vi9ejV79+6t8Lt9+/YRFRV1xs/n5OQwaNAgGjduzBdffIG/v/8Zr09MTMTlcrF3796zti0iPqq0Z2fTJnC7oQr/0uUpv7ugDR8t/o2NWTlMmruZZ4fGW1aLiC+odtiJiIggIuLsg+qSkpJwOp0sXbqUCy64AIAlS5bgdDrp379/pZ9zuVxcddVVBAYG8tVXXxEUFHTWe61cuZKgoCCaNGlS9S8iIr6lUyew2yEnx5yC3rq1ZaXYbQZPX3cew99ezL+XpnNX/7Z0igy1rB4Rb+exMTtxcXEMGjSIUaNGsXjxYhYvXsyoUaMYMmRI2Uys3bt307VrV5YuXQqYPTrJyckcOXKEadOm4XK5yMrKIisri+LiYgC+/vprpk6dytq1a9m2bRvvvPMO48eP59577yUwMNBTX0dEGrqAgBMLgq5fb20tQGKHZiR3i6K4xM3fvrP20ZqIt/Poe3amT59O9+7dSU5OJjk5mR49evDhhx+W/b6wsJBNmzaRl5cHwIoVK1iyZAlr1qyhU6dOREdHl22lg4r9/f2ZPHkySUlJ9OjRg9dff50JEybwyiuvePKriIg36NbN3G+oH++5+dOgrthtBvM27GXJ9gNWlyPitQy32+22uoi65nK5cDgcZdPaRcRHPPkkPPcc3Hsv/POfVlcDwPgv1jB9STrnxzThv/f3r9IEDhFfVdO/31obS0R8R2nPTj14jFXqoStiCQ6wsyrjMN+uybS6HBGvpLAjIr7j5LBTTzq1I0ODuO/ijgC8OHsTBUV60aBIbVPYERHf0aWLOeX84EHYt8/qasqMHNCe5qGBpB/MY/qS36wuR8TrKOyIiO9o1AjatzeP69GjrJBAPx6+wnzp4Rvfb8F5tNDiikS8i8KOiPiWejhuB+DWPq3p2DyEQ3mFvL1gm9XliHgVhR0R8S31bPp5KT+7jT8N6grAe7/sZH9uvsUViXgPhR0R8S1xcea+nvXsACR3i+L81g7yCoqZ/KN6d0Rqi8KOiPiWevoYC8zFk0sXCf1oyW9kOo9aXJGId1DYERHfUtqzk5UFhw5ZW8tpDIiN4IL24RQUlfDG91utLkfEKyjsiIhvCQ2FmBjzuJ727vzxKrN359PlGfx24IjFFYk0fAo7IuJ74uPN/Zo11tZRib7twrmkc3OKSty8Nm+L1eWINHgKOyLie0rDztq11tZxBo8eH7vz37TdbN6bY3E1Ig2bwo6I+J7u3c19PQ473Vs7GHReC9xumDRns9XliDRoCjsi4ntOfoxVT9bIOp1HkjtjGDB7XRZrdjmtLkekwVLYERHfExcHNpu5RlZWltXVVKpzVChDe7YC4OU5myyuRqThUtgREd8TFASxseZxPR2kXGrMFbH42Qx+2ryPZTsPWl2OSIOksCMivqkBDFIGaNsshFv7mlPlX/rfJtz1+LGbSH2lsCMivqkBDFIu9YeBnQjws7F0x0EWbtlvdTkiDY7Cjoj4pnr+rp2TRTsacUe/tgC8Om+zendEqklhR0R8U2nYWbcOSkqsraUKRl/agSB/GyvTDzN/8z6ryxFpUBR2RMQ3deoEgYFw9Cjs2GF1NWcVGRpESqLZu/PaXPXuiFSHwo6I+Ca7/cQK6A3gURbAfZd0pJG/nVW7nPywMdvqckQaDIUdEfFdpYOUV6+2to4qimgcyF392wEauyNSHQo7IuK7evY092lp1tZRDfde3IGQADtrd7uYu36v1eWINAgKOyLiuxpg2AkPCeDuC9sB8Oq8LZSUqHdH5GwUdkTEd51/vrnfsQMOH7a2lmoYNaADjQP92JDp4n/r6u9yFyL1hcKOiPiu8HBo08Y8biDjdgCaBAfw/y5qD8Br6t0ROSuFHRHxbQ3wURbAiIvaExrkx6a9Ocxam2l1OSL1msKOiPi2Bhp2HI38GXlRB8Ds3SlW745IpRR2RMS3NdCwA3DPRe1wNPJna3Yu36zeY3U5IvWWwo6I+LbSsLNuHRQUWFtLNYUF+XPvxWbvzuvztlBUXP+XvRCxgkfDzqFDh0hJScHhcOBwOEhJSeHwWWY8XHrppRiGUW4bPnz4ObcrInJa7dpBWJgZdDZutLqaarurfzuaBvuzff8Rvlql3h2R0/Fo2LnttttIS0tj9uzZzJ49m7S0NFJSUs76uVGjRpGZmVm2/fOf/6yVdkVEKjCMBv0oq3GgH/de3BGAN75X747I6fh5quENGzYwe/ZsFi9eTL9+/QCYOnUqSUlJbNq0iS5dulT62eDgYFq0aFFr7ebn55Ofn1/2s8vlOpevJiLepmdPWLAAVq6EO++0uppquzOpLe8s3M7OA3l8sXI3t/SJsbokkXrFYz07qampOByOskACkJiYiMPhYNGiRWf87PTp04mIiOC8887j0UcfJScn55zanThxYtkjL4fDQUyM/o9ARE7Sq5e5//VXa+uooZBAP+67xBy788YPWyhU745IOR4LO1lZWURGRlY4HxkZSVZW5W/8vP322/n444+ZP38+f/7zn5k5cyY33njjObU7btw4nE5n2ZaRkVGDbyQiXqtPH3O/YgUUF1tbSw2lJLYjonEgGQePMnPFLqvLEalXqh12nn766QoDiE/dli9fDoBhGBU+73a7T3u+1KhRo7jiiiuIj49n+PDhfPbZZ8ybN49fT/o3ruq2GxgYSFhYWLlNRKRM164QHAxHjsCmTVZXUyONAuz8/lJz7M7ff9hKQZF6d0RKVTvsPPDAA2zYsOGMW3x8PC1atGDv3oor8u7bt4+oqKgq36937974+/uzZcsWgFprV0SkjJ8f9O5tHi9bZm0t5+D2fm2IDA1k9+GjfLpCPdgipaoddiIiIujatesZt6CgIJKSknA6nSxdurTss0uWLMHpdNK/f/8q32/dunUUFhYSHR0NUGvtioiUU/oo63jPdEMU5G/n/uO9O2/+sJX8oob5SE6ktnlszE5cXByDBg1i1KhRLF68mMWLFzNq1CiGDBlSNmNq9+7ddO3atSy4bNu2jQkTJrB8+XJ27tzJrFmzuOWWW+jVqxcXXnhhldsVEam2vn3NfQPu2QEYfkEbWoQFkek8xoxl6t0RAQ+/Z2f69Ol0796d5ORkkpOT6dGjBx9++GHZ7wsLC9m0aRN5eXkABAQE8P3333PVVVfRpUsXHnzwQZKTk5k3bx52u73K7YqIVFtp2ElLg8JCa2s5B0H+dv5vYCcA/vHjVo4VqndHxHC73T63epzL5cLhcOB0OjVYWURMJSUQHg5OpzkFvXQ6egOUX1TMZS/NZ4/zGE9d2417LmxvdUkitaKmf7+1NpaICIDN5hXjdgAC/ew8MDAWgMnzt6l3R3yewo6ISKnSsNPAx+0A3JzQmtZNG7EvJ5+PFv9mdTkillLYEREp5SWDlAEC/Gw8eLx3Z8pP28grKLK4IhHrKOyIiJQqDTtr1pgvGGzgbujdijbhwezPLeDDVPXuiO9S2BERKRUTA61amUtGeEHvjr/dxoOXm707/1ywnSP56t0R36SwIyJSyjCg9OWkZ1mwuKEY2rMl7SNCOHikgH+l7rS6HBFLKOyIiJzs+AtMvSXs+NltPHi5+d6dtxdsx3Ws4b5DSKSmFHZERE5W2rOTmmq+e8cLXHd+KzpFNuZwXiH//Gmb1eWI1DmFHRGRk/XsCY0awcGDsHmz1dXUCrvN4E9XmcvpTPt5B9muYxZXJFK3FHZERE7m7w8XXGAe//KLtbXUoiu7RZHQtinHCkt47fstVpcjUqcUdkRETuVlg5QBDMPgsUFdAZixLIPt+3Itrkik7ijsiIicygvDDsAF7cO5vGskxSVuXp6zyepyROqMwo6IyKmSksz9xo1w4IC1tdSyPw3qimHArDVZpGUctrockTqhsCMicqpmzSAuzjxeuNDaWmpZlxah3NirNQAvfLcRt9ttcUUinqewIyJyOpddZu5/+MHaOjzg4StjCbDbSN1+gJ8277O6HBGPU9gRETmd0rDz44/W1uEBrZsGc2dSWwCen7WBomLveJ+QSGUUdkRETufSS8392rWwz/t6P/4wMJYmwf5s3pvLx8syrC5HxKMUdkRETiciArp3N4/nz7e0FE9wBPvzyJWdAZg0ZxPOo1pGQryXwo6ISGUGDjT3XjhuB+C2C9oQG9mYQ3mF/F0vGhQvprAjIlIZLx63A+YioU8O6QbA+4t26kWD4rUUdkREKnPxxWAYsGkT7NljdTUecUnn5lzapTlFJW6en7XR6nJEPEJhR0SkMk2bQq9e5rGXPsoCePKaOOw2g3kb9vLL1v1WlyNS6xR2RETO5IorzP3//mdtHR7UKTKUlERzKvqz36zXVHTxOgo7IiJncvXV5n72bCjx3hDw0OXmVPSNWTm8v2in1eWI1CqFHRGRM+nfH8LCYP9+WL7c6mo8pmlIAI8fXxX91bmbyXQetbgikdqjsCMicib+/pCcbB7PmmVtLR52a58YerdpwpGCYiZ8vd7qckRqjcKOiMjZDB5s7r/7zto6PMxmM3juhu7YbQbfrc3ih417rS5JpFYo7IiInM2gQeZ+2TLIzra2Fg+Liw5jxEXtAfjLl+s4WlBscUUi505hR0TkbFq2NKegu91ePSur1EOXx9LSEcSuQ0d580e9WVkaPoUdEZGqKJ2V9c031tZRB0IC/Xj6uvMAeHvBdrbszbG4IpFzo7AjIlIV111n7mfNgqPeP1Mp+bwWXBEXSWGxmz/NXE1xidvqkkRqzKNh59ChQ6SkpOBwOHA4HKSkpHD48OFKr9+5cyeGYZx2+/TTT8uuO93vp0yZ4smvIiK+rm9fiImB3FyYM8fqaurEhOvjCQ30Y2X6Yd5esN3qckRqzKNh57bbbiMtLY3Zs2cze/Zs0tLSSElJqfT6mJgYMjMzy23PPPMMISEhDC6dDXHce++9V+66u+66y5NfRUR8nWHAzTebx599Zm0tdaRlk0b85VpzodBX525mU5YeZ0nDZLjdbo/0TW7YsIFu3bqxePFi+vXrB8DixYtJSkpi48aNdOnSpUrt9OrVi969ezNt2rQTRRsGX3zxBUOHDq1SG/n5+eTn55f97HK5iImJwel0EhYWVo1vJSI+bdEiuPBC8yWD2dkQGGh1RR7ndrsZ+a/lfL8xm/hWYXxx/4X42zUCQqzhcrlwOBzV/vvtsX9iU1NTcTgcZUEHIDExEYfDwaJFi6rUxooVK0hLS2PEiBEVfvfAAw8QERFB3759mTJlCiVneI37xIkTyx6lORwOYmJiqv+FREQSE82ZWS4XzJtndTV1wp9rkJMAABcqSURBVDAMJt7YHUcjf9budvGPH7daXZJItXks7GRlZREZGVnhfGRkJFlZWVVqY9q0acTFxdG/f/9y55999lk+/fRT5s2bx/Dhwxk7dizPP/98pe2MGzcOp9NZtmVkZFTvy4iIANhscNNN5vFJ4wi9XWRYEBOuN2dnvfnDVtbudlpckUj1VDvsPP3005UOIi7dlh9fP8YwjAqfd7vdpz1/qqNHj/Lvf//7tL06Tz75JElJSfTs2ZOxY8cyYcIEXnrppUrbCgwMJCwsrNwmIlIjpeN2/vtfOHbM2lrq0HXnt2RwfAuKStyM/WQVxwr1skFpOPyq+4EHHniA4cOHn/Gadu3asXr1avburfiq8X379hEVFXXW+3z22Wfk5eVx5513nvXaxMREXC4Xe/furVLbIiI1dtFF5qysjAz48ksYNszqiuqEYRj8dWg8y3YeZNPeHJ75ej0Tb+xudVkiVVLtnp2IiAi6du16xi0oKIikpCScTidLly4t++ySJUtwOp0VHkudzrRp07juuuto3rz5Wa9duXIlQUFBNGnSpLpfR0Skemw2KJ39+d571tZSx5o1DuS1Yb0wDPh4aTr/Xbnb6pJEqsRjY3bi4uIYNGgQo0aNYvHixSxevJhRo0YxZMiQsplYu3fvpmvXruUCEcDWrVtZsGABI0eOrNDu119/zdSpU1m7di3btm3jnXfeYfz48dx7770E+sDMCBGpB+6+29zPmQO7dllaSl27KDaCPwyMBeCJL9awNVvT0aX+8+j8wenTp9O9e3eSk5NJTk6mR48efPjhh2W/LywsZNOmTeTl5ZX73LvvvkurVq1ITk6u0Ka/vz+TJ08mKSmJHj168PrrrzNhwgReeeUVT34VEZETOnaEiy8218r64AOrq6lzD10eS/+OzcgrKOb+6b+SV1BkdUkiZ+Sx9+zUZzWdpy8iUub99+GeeyA2FjZtMl866EP25eRz9RsL2ZeTz80JrXn5lvOtLkl8QL17z46IiFe7+WYICYEtW+Cnn6yups41Dw3kjeG9sBnw2Ypd/GdputUliVRKYUdEpCYaN+b/t3fvYVHV+R/A32fuglxUYhiVlUurYFlrkIiXrF1/aJnm7/ltrusla8tL4ir69DP2MW+bkKiPv25YK7mrv621fulatstj2uYdySRo1yBUEMwlBAxhBAWG+f7+ODKAYDIjM4eZeb+e5zxz5sz3jJ/5ODbvzhUtt7957TVla1FIfGQ/LPuPwQCAlz46jaxzVQpXRNQ5hh0iIkctXiw/fvwxUOydN8pMfORuTLm/PyxWgQXv5uBcxVWlSyLqgGGHiMhR0dHAhAnygcpvvql0NYqQJAkbfnkfYgb1Qe11C36z/Utcvtpw+xWJXIhhh4joTiQlyY/vvCPfM8sLGbRqbJ0dg5/09cGFH+ox7885vMIy9SgMO0REdyIhAYiKAsxmICND6WoU06+3Hn98+kH4GTTIKa3G8l3/hNXqdSf7Ug/FsENEdCdUKuCFF+T5DRuAujpl61HQ3cG98YdZMdCoJOz9ugxrP/kGXnh1E+qBGHaIiO7UU08BERFARQWwZYvS1Shq1N1B2PTk/ZAkYMeJUmzaX6h0SUQMO0REd0yrBVatkufT0uRdWl5s6vABWDf1XgBA+sEi/M+BM9zCQ4pi2CEi6g4zZ8pXU758GXj9daWrUdzMuEFY8Vg0AOC1f5zFZgYeUhDDDhFRd9BogNWr5fn164F/847gcx+KwEuT5MDzxufnkJpZwMBDimDYISLqLr/+NTByJHD1KrB8udLV9AjPjY3A6slDAQAZR8/jhQ//CUuzVeGqyNsw7BARdReVSr64oCQBf/kLcOSI0hX1CM+MDsfGX94HtUrC7q8u4rn/PYXa601Kl0VehGGHiKg7xcQA8+bJ84mJQAOvJgwAT8aG4g+zYmDQqnCosBL/mX4cxZW8tQS5BsMOEVF3S0kBgoKA06dbj+MhjB9qxIfzR8EUYEBRZR2eSD+OQ4UVSpdFXoBhh4iou/Xr13o15Q0bgKNHla2nBxk2MAAfLxqNmEF9YL5xL62tR4p44DI5FcMOEZEzTJ0K/OY38k1Cn3oKqK5WuqIeI9jPgL/MjcOvYkNhFUBq5rd4dscpVJivK10aeSiGHSIiZ3n1VSA8HCgpkc/UsliUrqjH0GvUWP9fw7B2yj3QqVX4/NsKTHz1KPadLle6NPJADDtERM7i5wf89a+Ajw/w6afAiy8qXVGPIkkS5owKw97fjka0yR8/1DViwbs5eOHDr2Hm2VrUjRh2iIic6Wc/A7Zvl+c3b/b6e2d1JirEHx8ljsKCcZGQJGBXzkWM33wYH+X+m8fyULdg2CEicrYnnwTWrJHnExOBP/5R0XJ6Ir1GjeRHo/DBvHj8pK8PLtU2IOmDPDz59gmc/neN0uWRm5OEF8bm2tpaBAQEoKamBv7+/kqXQ0TeQAhg6VLgtdfkiw7+6U/AnDlKV9UjXW9qxjtHi5F+sAjXmpohScD0B0Ox+Bc/hSmgl9LlkYIc/f1m2GHYISJXEQJYuBB4+235+csvAytWyOGHOvi+5hpeyfwWe78uAwDo1CpMe3AgFoyLxMA+PgpXR0pg2LEDww4RKcZqBZKTgY0b5edPPw2kp8sHMVOnTp7/AZs+LcTJkh8AABqVhF/GDMT8cZEID/JVuDpyJYYdOzDsEJHitmwBfvtbOfzccw/wwQfyI91SdvFlvP6Ps8gqumxbNvanQZgZ9xP8ItoIrZqHoXo6hh07MOwQUY/w2WfArFnApUuAwQC89BLwwguAXq90ZT1aTukPSD9YhIOFFWj5BTP66/Gr2FA8MXwAIu/qrWyB5DQMO3Zg2CGiHuPSJXlX1r598vOoKHkX16RJPJbnNr77oR47T17A/536DlVXG23Lo03+ePw+EyYNMyGMu7k8CsOOHRh2iKhHEQLYuVM+W6vixo0xR44EVq4EJk4EVNw982MaLVbs+6Ycf/3qIo6drYLF2vqzNsToh4cGB+GhwXfhwbC+MGjVClZKd4phxw4MO0TUI1VXA2lpwOuvA9euycsGDwaefx6YMQMIDla2Pjdwpb4R+7+5hE/+WYasostobhN89BoVRoT3xYiwvogJ64OfhQbCR6dRsFqyF8OOHRh2iKhHKy8HNm2S75xeWysvU6uBhAR599bEiUBkpLI1uoHqukYcO1eFI2cqceRsJS7VNrR7Xa2SEG3yw/0DAxFt8sfQ/v6ICvFjAOrBGHbswLBDRG7BbAb+/Gf5dhNfftn+tbvvlkPPmDHAAw/I4Ye7u25JCIEzl64iq6gKOaXV+Kq0GmU1He+yLknAoL4+iLirNwb180F4kC/C+slT/0ADNDzjS1E9MuykpKTg73//O/Ly8qDT6XDlypXbriOEwNq1a7F161ZUV1cjLi4O6enpuKfNKZnV1dVYvHgx9u7dCwCYMmUK3njjDQQGBnapLoYdInI7BQXAxx/LNxQ9dqzjHdT9/IDhw4H775eDUGSkPIWH8+yuWyi7cg1fXajGN2W1yC+rRcH3tagwN9xyvFYtYUBgLxj9DTD6GxASYECwn942389Xh0AfHfwNGoYiJ+mRYWf16tUIDAzExYsXsW3bti6FnbS0NKSkpGD79u0YPHgw1q1bhyNHjqCwsBB+fn4AgEcffRQXL17E1q1bAQDz5s1DWFgYPvnkky7VxbBDRG6tthY4eBA4cAA4dQr4+mvgesetFADkTRVBQfLxPkZj+8fgYMDfXw5KbSd/f8DXVw5JXra1qOpqA86Um3H+ch1KL9fjfFUdSi/XoeRyPRot1i6/j79Bg0AfHQJ9tAjopZXne2nR26CBj1YNH70GPjr1jUkDX50avXRq+Oo1MGjU0GokaNUq6DQq6NQqaNUqqFU8O69Hhp0W27dvR1JS0m3DjhAC/fv3R1JSEl588UUAQENDA4xGI9LS0jB//nwUFBRg6NChyM7ORlxcHAAgOzsb8fHx+PbbbzFkyJAO79vQ0ICGhta0Xltbi9DQUIYdIvIMFou85eerr4DTp4Giotapru7O3lutlkOPTtf62HZerwe0WnmcStX66Oh08+n2d/q8O95DkiCEQF1jM65et6DeYsW1RgvqGppR39SM+kZ5ut7UjKbmrgcie0mSBJUEqCUJKpUElSRBrWqZB1QSIEECJEAFyfYxVFLb+TZjJAkS0GG+5dO3dkFq15IO7bmx0s2d1z40BjH/vaDbPj/geNjpUUdhnT9/HuXl5UhISLAt0+v1GDduHLKysjB//nycOHECAQEBtqADACNHjkRAQACysrI6DTuvvPIK1q5d65LPQETkchoNMGyYPLUlBFBZKR/wfOmSfFp7RUXrfGWlvJXIbG4/1de3vkdzs/y87TIvJAHofWOirvmisRHo5rDjqB4VdsrLywEARqOx3XKj0YjS0lLbmOBOTr8MDg62rX+z3/3ud1i2bJntecuWHSIijyZJrbur7GGxyOGmsVGeGhp+fL6hQb7txe2m5ubbv97WzTse7H3eHe/hyJ/ZTYQQsAqg2SrQbBWwCvmx+caj1YrWeSFPAGAV8rpCAAI3Hm8ss96o2SoAgbbjbszf9JFEazEQ7Wrr+Hrb5wKAYdTI7m6Jw+wOO2vWrLntVpIvv/wSsbGxDhcl3bSNTAjRbtnNr3c2pi29Xg89D9AjIuoajUY+bocUJQFQ35joztgddhYtWoTp06f/6JiwsDCHigkJCQEgb70xmUy25RUVFbatPSEhIbh06VKHdSsrKztsESIiIiKyO+wEBQUhKCjIGbUgPDwcISEhOHDgAIYPHw4AaGxsxOHDh5GWlgYAiI+PR01NDU6ePIkRI0YAAL744gvU1NRg1KhRTqmLiIiI3JdTzym8cOEC8vLycOHCBTQ3NyMvLw95eXm4evWqbUxUVBT27NkDQN49lZSUhNTUVOzZswenT5/G008/DR8fH8yYMQMAEB0djYkTJ2Lu3LnIzs5GdnY25s6di8cff7zTg5OJiIjIuzn1AOVVq1Zhx44dtuctW2sOHjyIhx9+GABQWFiImpoa25jly5fj2rVrWLhwoe2igvv377ddYwcA3nvvPSxevNh21taUKVPw5ptvOvOjEBERkZvi7SJ4EB4REZFbcPT327sujUlEREReh2GHiIiIPBrDDhEREXk0hh0iIiLyaAw7RERE5NEYdoiIiMijMewQERGRR2PYISIiIo/m1Cso91Qt11Gsra1VuBIiIiLqqpbfbXuvh+yVYcdsNgMAQkNDFa6EiIiI7GU2mxEQENDl8V55uwir1YqysjL4+flBkqRufe/a2lqEhobiu+++460onIh9dg322TXYZ9dgn13HWb0WQsBsNqN///5Qqbp+JI5XbtlRqVQYOHCgU/8Mf39//mNyAfbZNdhn12CfXYN9dh1n9NqeLToteIAyEREReTSGHSIiIvJo6jVr1qxRughPo1ar8fDDD0Oj8cq9hC7DPrsG++wa7LNrsM+u05N67ZUHKBMREZH34G4sIiIi8mgMO0REROTRGHaIiIjIozHsEBERkUdj2CEiIiKPxrDjgC1btiA8PBwGgwExMTE4evToj47fvXs3hg4dCr1ej6FDh2LPnj0uqtS92dPnjIwMjB07Fn369EGfPn0wfvx4nDx50oXVui97v88t3n//fUiShKlTpzq5Qs9gb5+vXLmCxMREmEwmGAwGREdHIzMz00XVui97+/zqq69iyJAh6NWrF0JDQ7F06VJcv37dRdW6pyNHjmDy5Mno378/JEnCRx99dNt1Dh8+jJiYGBgMBkRERODtt992QaVtCLLL+++/L7RarcjIyBD5+fliyZIlwtfXV5SWlnY6PisrS6jVapGamioKCgpEamqq0Gg0Ijs728WVuxd7+zxjxgyRnp4ucnNzRUFBgXjmmWdEQECAuHjxoosrdy/29rlFSUmJGDBggBg7dqx44oknXFSt+7K3zw0NDSI2NlY89thj4tixY6KkpEQcPXpU5OXlubhy92Jvn999912h1+vFe++9J86fPy8+/fRTYTKZRFJSkosrdy+ZmZlixYoVYvfu3QKA2LNnz4+OLy4uFj4+PmLJkiUiPz9fZGRkCK1WK3bt2uWiioVg2LHTiBEjxIIFC9oti4qKEsnJyZ2OnzZtmpg4cWK7ZRMmTBDTp093Wo2ewN4+38xisQg/Pz+xY8cOZ5TnMRzps8ViEaNHjxbvvPOOmDNnDsNOF9jb57feektERESIxsZGV5TnMeztc2Jiovj5z3/ebtmyZcvEmDFjnFajp+lK2Fm+fLmIiopqt2z+/Pli5MiRziytHe7GskNjYyNycnKQkJDQbnlCQgKysrI6XefEiRMdxk+YMOGW48mxPt+svr4eTU1N6Nu3rzNK9AiO9vn3v/897rrrLjz77LPOLtEjONLnvXv3Ij4+HomJiTAajbj33nuRmpqK5uZmV5Tslhzp85gxY5CTk2Pb5V1cXIzMzExMmjTJ6fV6k1v9Dp46dQpNTU0uqUH5azi7kaqqKjQ3N8NoNLZbbjQaUV5e3uk65eXldo0nx/p8s+TkZAwYMADjx493RokewZE+Hz9+HNu2bUNeXp4rSvQIjvS5uLgYn3/+OWbOnInMzEycPXsWiYmJsFgsWLVqlSvKdjuO9Hn69OmorKzEmDFjIISAxWLB888/j+TkZFeU7DVu9TtosVhQVVUFk8nk9BoYdhwgSVK750KIDsvuZDzJHO3bhg0bsHPnThw6dAgGg8FZ5XmMrvbZbDZj1qxZyMjIQFBQkKvK8xj2fJ+tViuCg4OxdetWqNVqxMTEoKysDBs3bmTYuQ17+nzo0CGkpKRgy5YtiIuLw7lz57BkyRKYTCasXLnSFeV6jc7+Xjpb7iwMO3YICgqCWq3u8H8JFRUVHVJri5CQELvGk2N9brFp0yakpqbis88+w3333efMMt2evX0uKipCSUkJJk+ebFtmtVoBABqNBoWFhYiMjHRu0W7Ike+zyWSCVquFWq22LYuOjkZ5eTkaGxuh0+mcWrM7cqTPK1euxOzZs/Hcc88BAIYNG4a6ujrMmzcPK1asgErFIz26w61+BzUaDfr16+eSGvg3aQedToeYmBgcOHCg3fIDBw5g1KhRna4THx/fYfz+/ftvOZ4c6zMAbNy4ES+//DL27duH2NhYZ5fp9uztc1RUFP71r38hLy/PNk2ZMgWPPPII8vLyEBoa6qrS3Yoj3+fRo0fj3LlztjAJAGfOnIHJZGLQuQVH+lxfX98h0KjVagj55B2n1eptbvU7GBsbC61W65oiXHYotIdoObVx27ZtIj8/XyQlJQlfX19RUlIihBBi9uzZ7Y78P378uFCr1WL9+vWioKBArF+/nqeed4G9fU5LSxM6nU7s2rVLfP/997bJbDYr9RHcgr19vhnPxuoae/t84cIF0bt3b7Fo0SJRWFgo/va3v4ng4GCxbt06pT6CW7C3z6tXrxZ+fn5i586dori4WOzfv19ERkaKadOmKfUR3ILZbBa5ubkiNzdXABCbN28Wubm5tlP8k5OTxezZs23jW049X7p0qcjPzxfbtm3jqefuID09XQwaNEjodDrxwAMPiMOHD9teGzdunJgzZ0678R9++KEYMmSI0Gq1IioqSuzevdvFFbsne/o8aNAgAaDDtHr1atcX7mbs/T63xbDTdfb2OSsrS8TFxQm9Xi8iIiJESkqKsFgsLq7a/djT56amJrFmzRoRGRkpDAaDCA0NFQsXLhTV1dUKVO4+Dh482Ol/b1t6O2fOHDFu3Lh26xw6dEgMHz5c6HQ6ERYWJt566y2X1iwJwW11RERE5Ll4zA4RERF5NIYdIiIi8mgMO0REROTRGHaIiIjIozHsEBERkUdj2CEiIiKPxrBDREREHo1hh4iIiDwaww4RERF5NIYdIiIi8mgMO0REROTR/h+AL4fLSMBh6gAAAABJRU5ErkJggg==",
      "text/plain": [
       "PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x13f99c790>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " 5: 0.000970486365,\t 2.405141480821\n",
      "k= 6: 0.000706687190,\t 0.942200691795\n",
      "k= 7: 0.000110332384,\t 0.116033574832\n",
      "k= 8: 0.000001637596,\t 0.001522042857\n",
      "k= 9: 0.000000000183,\t 0.000000222835\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "1-element Array{PyCall.PyObject,1}:\n",
       " PyObject <matplotlib.lines.Line2D object at 0x13f9b6dd0>"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x=linspace(0,1,201)\n",
    "u=zeros(201)\n",
    "\n",
    "h=x[2]-x[1]\n",
    "ϵ=0.0275\n",
    "\n",
    "#specify the boundary conditions\n",
    "u[1]=1.0\n",
    "u[end]=-1.0\n",
    "\n",
    "GaussSeidel_Burgers!(u,h,ϵ,750)\n",
    "\n",
    "plot(x, u)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Newton's method for systems\n",
    "$f:\\mathbb{R}^d \\mapsto \\mathbb{R}^d$; i.e. $f$ maps a vector of $d$ component, x, to a vector of $d$ component, denoted by $f(x_k)$. \n",
    "\n",
    "$f(x_*)=0$. Taylor expansion gives\n",
    "$$f(x_*) = 0=f(x_k)+f^\\prime(x_k)(x_*-x_k)+\\ldots$$\n",
    "\n",
    "Replace $x_*$ by $x_{k+1}$ yields Newton's method:\n",
    "$$ f^\\prime(x_k)(x_{k+1}-x_k))=-f(x_k), $$\n",
    "where  the derivative $f^\\prime(x)$ is a $d\\times d$ matrix.\n",
    "\n",
    "Here, the $g$ function that defines the fixedpoint iterations is: \n",
    "$$ g(x) := x-f^\\prime(x)^{-1} f(x) $$\n",
    "\n",
    "Typical algorithm would do:\n",
    "* Form $A_k:= f^\\prime(x_k)$ (this could be already expensive to compute\n",
    "* Solve $A_k z = -f(x_k)$, and obtain $z$\n",
    "* Update $x_{k+1} := z+x_k$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Newton! (generic function with 1 method)"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function Newton!(x, F::Function, JF::Function, h, ϵ; nitrs=10, tol=1e-12)\n",
    "   \n",
    "    for k in 1:nitrs\n",
    "        \n",
    "        f=F(x)\n",
    "        \n",
    "               \n",
    "        if norm(f)<tol\n",
    "           break \n",
    "        end\n",
    "        \n",
    "        #Solve BΔx=-F(x)  for Δx\n",
    "        Δx=-JF(x)\\f\n",
    "\n",
    "        #Update x\n",
    "        x.+=Δx\n",
    "        \n",
    "        @printf \"k=%2d: %1.12f,\\t %1.12f\\n\" k norm(f) norm(Δx)\n",
    "\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Worked example: Problem 5.22 in Gander-Gander-Kwok\n",
    "\n",
    "* radius of pipe: $r=4$\n",
    "\n",
    "* rope length: $L=30$\n",
    "\n",
    "Find the distance of the pipe to the ceiling above: $h$.\n",
    "\n",
    "## Setting up equations: \n",
    "\n",
    "* Need to variables to describe the system. $h$ and $\\theta$\n",
    "\n",
    "* radius of pipe = r  = $(h+r)\\cos\\theta$\n",
    "\n",
    "* length of rope = L = (length of rope wrapping around the pipe) + (legth of rope in the air) $$L = 2\\left( r(\\pi - \\theta) + (h+r)\\sin\\theta \\right)$$\n",
    "\n",
    "* Obtain two equations: \n",
    "$$f(h,\\theta):=\\left(\\begin{array}{c}\n",
    "f_{1}(h,\\theta)\\\\\n",
    "f_{2}(h,\\theta)\n",
    "\\end{array}\\right)=\\left(\\begin{array}{c}\n",
    "(h+4)\\cos\\theta-4\\\\\n",
    "(h+4)\\sin\\theta+4(\\pi-\\theta)-15\n",
    "\\end{array}\\right)=\\left(\\begin{array}{c}\n",
    "0\\\\\n",
    "0\n",
    "\\end{array}\\right)\n",
    "$$\n",
    "\n",
    "* Derivative of $f$: $$f^{\\prime}(h,\\theta)=\\left(\\begin{array}{cc}\n",
    "\\frac{\\partial f_{1}}{\\partial h} & \\frac{\\partial f_{1}}{\\partial\\theta}\\\\\n",
    "\\frac{\\partial f_{2}}{\\partial h} & \\frac{\\partial f_{2}}{\\partial\\theta}\n",
    "\\end{array}\\right)=\\left(\\begin{array}{cc}\n",
    "\\cos\\theta & -(h+4)\\sin\\theta\\\\\n",
    "\\sin\\theta & (h+4)\\cos\\theta-4\n",
    "\\end{array}\\right)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "f(h,θ)  = [ (h+4.0)*cos(θ)-4.0; (h+4.0)*sin(θ)+4.0*(π-θ)-15.0]\n",
    "\n",
    "f′(h,θ) = [cos(θ) -(h+4.0)*sin(θ); sin(θ)  (h+4.0)*cos(θ)-4.0] \n",
    "\n",
    "Newton_prob522(x) = x - f′(x[1],x[2])\\f(x[1],x[2]) #this is not usually how one should run Newton's method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#x is an array used to store the approximations computed along the iterations\n",
    "#solution at the k-th iteration is x[:,k]\n",
    "#function assumes that x[:,1] is the assigned initial guess\n",
    "function FixedPointItrs!(x::Array{Float64,2}, g)\n",
    "    for k=1:size(x,2)-1\n",
    "        x[:,k+1]=g(x[:,k])\n",
    "    end\n",
    "end\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x0=[ 9.0; π/3.0]\n",
    "x=zeros(2, 7)\n",
    "x[:,1]=x0\n",
    "\n",
    "FixedPointItrs!(x,Newton_prob522)\n",
    "\n",
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "println(\"Verify that the equations are solved (f should be essentially 0): f = \", f(x[1, end], x[2,end]))\n",
    "\n",
    "h, θ=x[:, end]\n",
    "\n",
    "if norm(f(x[1, end], x[2,end]))<1e-10\n",
    "    println(\"\\nFound a solution: h=\",round(h,10), \", θ=\", round(θ,10))\n",
    "else\n",
    "    println(\"\\nNewton method fails.\")\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example: solving the steady-state viscous Burger's equation\n",
    "\n",
    "$$ uu_x = \\epsilon u_{xx},\\,\\,\\,\\,\\,0<x<1, $$ \n",
    "with boundary conditions\n",
    "$$u(0)=1,\\,\\,u(1)=-1. $$\n",
    "\n",
    "We are lead to a system involving the unknowns $u_1,u_2,\\cdots, u_{n-1}$ with the boundary conditions $u_0=1$ and $u_n=-1$:\n",
    "\n",
    "$$ u_j \\frac{u_{j+1}-u_{j-1}}{2h} = \\epsilon\\frac{u_{j+1}-2u_j+u_{j-1}}{h^2}, j=1,2,\\cdots,n-1.\\,\\,\\,\\,\\, (P1)$$\n",
    "\n",
    "The above equations can be organized into:\n",
    "$$f(u_1,u_2,\\cdots,u_{n-1})=(f_1(u_1,u_2,\\cdots,u_{n-1}), f_2(u_1,u_2,\\cdots,u_{n-1}),\\cdots, f_{n-1}(u_1,u_2,\\cdots,u_{n-1}))^\\top = 0,$$\n",
    "where\n",
    "$$ f_j(u_1,u_2,\\cdots,u_{n-1}):=u_j \\frac{u_{j+1}-u_{j-1}}{2h} - \\epsilon\\frac{u_{j+1}-2u_j+u_{j-1}}{h^2}.$$\n",
    "\n",
    "In the following, we use Newton's method to solve this nonlinear system, with initial conditions computed by a few Gauss-Seidel iterations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "NewtonB! (generic function with 1 method)"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function FBurgers(U, h, ϵ)\n",
    "    \n",
    "    F=copy(U)\n",
    "    \n",
    "    for j=2:length(U)-1\n",
    "        F[j]=h*U[j]*(U[j+1]-U[j-1])/2-ϵ*(U[j+1]-2U[j]+U[j-1])\n",
    "    end\n",
    "    \n",
    "    F[2:end-1]\n",
    "end\n",
    "\n",
    "function JFBurgers(U, h, ϵ) \n",
    "    \n",
    "    JF=zeros(length(U), length(U))\n",
    "    \n",
    "    for j=2:length(U)-1\n",
    "        JF[j,j-1]=-h*U[j]/2-ϵ\n",
    "        JF[j,j]  = h*(U[j+1]-U[j-1])/2+2ϵ \n",
    "        JF[j,j+1]= h*U[j]/2-ϵ    \n",
    "    end\n",
    "    \n",
    "    JF[2:end-1,2:end-1]\n",
    "end\n",
    "\n",
    "function NewtonB!(x, F::Function, JF::Function, h, ϵ; nitrs=10, tol=1e-12)\n",
    "   \n",
    "    for k in 1:nitrs\n",
    "        \n",
    "        f=F(x,h,ϵ)\n",
    "        \n",
    "               \n",
    "        if norm(f)<tol\n",
    "           break \n",
    "        end\n",
    "        \n",
    "        #Solve BΔx=-F(x)  for Δx\n",
    "        Δx=-JF(x,h,ϵ)\\f\n",
    "\n",
    "        #Update x\n",
    "        x[2:end-1].+=Δx\n",
    "        \n",
    "        @printf \"k=%2d: %1.12f,\\t %1.12f\\n\" k norm(f) norm(Δx)\n",
    "\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjsAAAGgCAYAAABMn6ZGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzs3Xl8VOXd///XmUkyIdtADFmAkATZggGFICRY3O8IitbaBWuN2ltRvr2tpWprKfYWscpPW73thlRFqdYqtdbWrVRcECphEUH2gAiEJQuEZCYhMFlmfn+cZEjIQgKZTGbm/Xw8zmNmzrnOmc94P27y7nWu61yGx+PxICIiIhKkLP4uQERERMSXFHZEREQkqCnsiIiISFBT2BEREZGgprAjIiIiQU1hR0RERIKawo6IiIgENYUdERERCWoKOyIiIhLUFHZEREQkqCnsiIiISFAL83cB/uB2uzl06BCxsbEYhuHvckRERKQTPB4PVVVVDBgwAIul8/01IRl2Dh06RGpqqr/LEBERkTOwf/9+Bg0a1On2IRl2YmNjAfM/VlxcnJ+rERERkc5wOp2kpqZ6/453VkiGnaZbV3FxcQo7IiIiAaarQ1A0QFlERESCmsKOiIiIBDWFHREREQlqCjsiIiIS1BR2REREJKgp7IiIiEhQU9gRERGRoKawIyIiIkFNYUdERESCmsKOiIiIBDWfhp0VK1Zw7bXXMmDAAAzD4B//+Mdpz/nkk0/Izs4mMjKSIUOGsHDhwlZtFixYQEZGBpGRkWRnZ7Ny5UpflC8iIiJBwKdh59ixY5x//vn8/ve/71T7PXv2cPXVVzN58mQ2bNjAz3/+c+655x7eeOMNb5slS5Ywa9Ys5syZw4YNG5g8eTJTp06lqKjIVz9DREREApjh8Xg8PfJFhsGbb77J9ddf326bBx54gLfeeovt27d7982cOZMvvviCgoICACZOnMi4ceN45plnvG0yMzO5/vrrmT9/fpvXdblcuFwu7+emVVMdDke3LgS64+9LqXzxz97PTeuUeWi5YJlhtNzX9M5jGC1bGga0tdZZiwXQjDb2mXs9RosWjR/ayLdtLajWiX1GY82tDxmt33bhO1pfqnU7w2Lx/ucxDMN8b1iwNDa1NO2zGBgY5n4DLBjmvsbjlsbzsZhtLIaBxTAIsxhYrBasFoMww8BiNbAaBlarBcupdXf2t1mtEB4OYWGtX9vaFx5ubtHRLbeoKPNaIiIhxul0Yrfbu/z3u1etel5QUEBeXl6LfVdddRWLFi2irq4Oj8fD+vXr+dnPftaiTV5eHqtWrWr3uvPnz+fhhx/2Sc3NOdZtJOedV3z+PSJERprBJzYWEhJabykpkJ4OaWmQmgoREf6uWETEb3pV2CkpKSEpKanFvqSkJOrr6zly5Agej4eGhoY225SUlLR73dmzZ3Pvvfd6Pzf17HS3fhdNoODGu8wP3v6yNjrO3Cf3eZodN5o19Xj3eVq1o63OuFa7PG23O7Wtt01b12xe58l6Tr2s0Xi0vf3Njxm0Pr+tOg2Pp42f1LIej6fxv0tjTU2fzfce8z9Bi88e3I0nN+2n8Xs8jddwe8DjduP2eHB7oMHjwd34f682O9naqr2N/5aGByyeBsLcbsLc9YS5G5pt9YS53dhwE4mbCNxEehqIwENEQx0RtScIP1GDpabm5PedOGFu5eWwd28blTX/cgMGDoRRoyArC847D8aPN1/VQyQiIaBXhR1ovKXQTNNdNsMwWrw/tc2p+5qz2WzYbLZurrS14dMuh2mX+/x7pGd5PB5qG9zU1rupa/A0vrpxNb7W1rupbXBzoq6BY64Gamrrqak1X1t+bqDaVU/1iXoqampxHK+joqaWE3XuzhZCtLuWc6MMhkYbDOkDaeENpHGcNM9x4qorMcrL4fBhOHAA9u0ztxMnzM8HDsD775+8Xlwc5ObCJZfAtdea4aeD/z8SEQlUvSrsJCcnt+qhKSsrIywsjHPOOQePx4PVam2zzam9PSLdxTAMbGFWbGG+6QU5UddAZU0dlcdrqThWx5FqF6XOE5Q6T1DidFHqOEFp1QlKHCc4Vm9jkws2eYeghQORQD8SYtLJHBbHqEviOH9QXy5Mj6d/TIQZfnbvhq1bYcsW2LQJ1q0DpxP+/W9z+/nPYcgQ+PrXIT8fxo71yW8VEfGHXhV2cnNzefvtt1vse//99xk/fjzh4eEAZGdns2zZMr7xjW942yxbtoyvf/3rPVqrSHeJDLeSbLeSbI/ssJ3H46HU6WLPkWPsKz/GnvJj7DtSw66yKvYcOcaR6lpW7jrCyl1HvOdkJEQzIT2eCRmpXDx9HP3vaOzhrK+HzZvh009h6VL44AP46iv4v/8zt+xs+OEP4aabzEHSIiIBzKezsaqrq/nyyy8BGDt2LE899RSXXXYZ8fHxDB48mNmzZ3Pw4EFeeuklwJx6npWVxV133cWMGTMoKChg5syZvPrqq3zzm98EzKnn+fn5LFy4kNzcXJ599lmee+45tm7dSlpaWqfqOtPR3CK91fHaBgpLq9h2yMnWQw7W76ugsLSqxXAow4DzB/XlipGJTB2dzNDE2JMHq6vNW1yvvQb/+AfU1Zn7Bw+G2bPhjjvMGWIiIn50pn+/fRp2li9fzmWXXdZq/6233srixYu57bbb2Lt3L8uXL/ce++STT/jxj3/M1q1bGTBgAA888AAzZ85scf6CBQt44oknKC4uJisri//7v//j4osv7nRdCjsSChw1dXy27yhr9x7l0y+PsOWgs8XxrIFxXH/BQK67YACJsc16lY4cgUWLzB6e0tLGxlnwm9/A5RqTJiL+0yvDTm+lsCOhqNR5go93lPH+tlJW7DxMfeMsszCLwbQxKfz31zIYM6jvyROOH4fnn4e5c+HoUXPfnXfCU0+Z095FRHqYwk4XKOxIqDt6rJZ3Nx3ijc8PsnF/pXf/+LR+/PCKYVw8LOHkDMejR+F//xf+8Afz89Ch8PrrcMEFfqhcREKZwk4XKOyInLT5gIMXP93D25sOUddg/nOQMySeB6aMZOzgficbfvQR3HqrOYU9Ohr+8he47jo/VS0ioUhhpwsUdkRaK3Oe4I8rvuLlgn3UNpjP/rlmdAr/e+0okuIax/RUVMD06bBsmTni+emn4Z57/Fi1iIQShZ0uUNgRad/ByuM8vWwnb3x+ALcHYm1h/HTqSG6eONi8tVVXZwachQvNE377W3OauoiIj53p32+frnouIoFnYN8+/Orb5/PODydzfmpfqlz1/OIfW/j+4nUcrnKZz91ZsADmzDFPaB58RER6IYUdEWnTqAFx/P3/TeKha0dhC7OwvPAwU3+zgoLd5eYtrEcegZ/8xGz8gx/Au+/6t2ARkXYo7IhIu6wWg+9flMFbd3+NEUmxHKmuJX/RGl5evc8MPI8/bk5H93jMpy3v2OHvkkVEWlHYEZHTGpEcyz/vvojrzh9AvdvDL/6xhXlvb8PtAX73O5g82Vxr6+tfB4fD3+WKiLSgsCMinRIZbuU3N17AT64aAcALn+7hp29sot4aBn/7m7m0xM6dMGuWnysVEWlJYUdEOs0wDP7nsqE8+e3zsRjwt/UH+NFrG6k/J8F87o5hwOLF8NZb/i5VRMRLYUdEuuyb2YNY8L1sIqwW3t1czM/f3Ixn0iS4/36zwZ13Qnm5f4sUEWmksCMiZ2RKVjK//e5YLAb89bMDzP/XDpg3D0aNMhcQbQo+IiJ+prAjImdsSlYy/98NYwB4dsVX/HljqbliOsCf/gSff+7H6kRETAo7InJWvnNhqnfQ8ty3trI2aTh873vmdPQf/9h8FRHxI4UdETlrP7j0XKaNSaHe7eH//Xk9JT97CPr0gRUr4O9/93d5IhLiFHZE5KwZhsET3xpDZkoc5cdquefTI7jvu888OHs2NDT4t0ARCWkKOyLSLaIiwlh48ziiIqys3XOU5yd9C+LjYdcueP11f5cnIiFMYUdEuk3aOdHMve48AJ5YeZDS799lHnjsMY3dERG/UdgRkW717exBXD06mXq3hx/0uwhPbCxs3gzvvOPv0kQkRCnsiEi3MgyDR68fTXx0BOurYOO075oHHn1UvTsi4hcKOyLS7fpFR/DzqzMBuDvpEtyRkbBmDaxd6+fKRCQUKeyIiE98c9xAJmbEc9BmpyD7CnPnwoX+LUpEQpLCjoj4hGEYPPqN0YRbDZ7MuNTc+dprUFHh17pEJPQo7IiIzwxNjOF7E9P4fMBI9gw4F06cgJde8ndZIhJiFHZExKd+ePlQYiLDWXRenrnjj3/UQGUR6VEKOyLiU+fE2Ljr4iH847zLqImIhO3b4T//8XdZIhJCFHZExOdun5xBn4R+vDPia+aOv/zFvwWJSEhR2BERn4uKCOOui4fwVuYlAHhefx3q6vxclYiECoUdEekR350wmB0jszkc1RejvBw+/NDfJYlIiFDYEZEeEW0L45aLh/LeyIsA8Lz6mp8rEpFQobAjIj3m1tx0Phh9KQANf/+7ORVdRMTHeiTsLFiwgIyMDCIjI8nOzmblypXttr300ksxDKPVds0113jb3Hbbba2O5+Tk9MRPEZGzYI8KZ9S3pnIoNoGw6ir417/8XZKIhACfh50lS5Ywa9Ys5syZw4YNG5g8eTJTp06lqKiozfZ///vfKS4u9m5btmzBarXy7W9/u0W7KVOmtGj33nvv+fqniEg3uHlSBu9kXgxA1cuv+rkaEQkFPg87Tz31FLfffjt33HEHmZmZPP3006SmpvLMM8+02T4+Pp7k5GTvtmzZMqKiolqFHZvN1qJdfHx8uzW4XC6cTmeLTUT8IzU+iqNXTgUgbNm/ob7ezxWJSLDzadipra1l/fr15OXltdifl5fHqlWrOnWNRYsWceONNxIdHd1i//Lly0lMTGT48OHMmDGDsrKydq8xf/587Ha7d0tNTe36jxGRbnPRzddQERlLn2onNZ/oAYMi4ls+DTtHjhyhoaGBpKSkFvuTkpIoKSk57flr165ly5Yt3HHHHS32T506lVdeeYWPPvqIJ598knXr1nH55ZfjcrnavM7s2bNxOBzebf/+/Wf+o0TkrF00IpnPRk0E4KvFS/xcjYgEu7Ce+BLDMFp89ng8rfa1ZdGiRWRlZTFhwoQW+6dPn+59n5WVxfjx40lLS+Pdd9/lhhtuaHUdm82GzWY7w+pFpLtZLAbh110Ln39AzAf/9nc5IhLkfNqzk5CQgNVqbdWLU1ZW1qq351Q1NTW89tprrXp12pKSkkJaWhq7du06q3pFpOeMvWM69YaF9JI97Fq9yd/liEgQ82nYiYiIIDs7m2XLlrXYv2zZMiZNmtThuX/9619xuVzcfPPNp/2e8vJy9u/fT0pKylnVKyI9xz4wia9GXADAl7qVJSI+5PPZWPfeey/PP/88L7zwAtu3b+fHP/4xRUVFzJw5E4BbbrmF2bNntzpv0aJFXH/99Zxzzjkt9ldXV3P//fdTUFDA3r17Wb58Oddeey0JCQl84xvf8PXPEZFuZEybBoD9w/epa3D7uRoRCVY+H7Mzffp0ysvLmTdvHsXFxWRlZfHee++RlpYGQFFRERZLy8y1c+dO/vOf//D++++3up7VamXz5s289NJLVFZWkpKSwmWXXcaSJUuIjY319c8RkW6Unv8t+PU8xu3dxIrNB7jigsH+LklEgpDh8Xg8/i6ipzmdTux2Ow6Hg7i4OH+XIxK6PB6q4/sTU1nOk3Oe5b5fzvB3RSLSi53p32+tjSUi/mMYNFxyKQARnyzHeaLOv/WISFBS2BERv4q72nzo6IS9X/DBtlI/VyMiwUhhR0T8yrjiCgDGHirkg8/2+LkaEQlGCjsi4l9DhlA3MJUIdz3Hl6+gSreyRKSbKeyIiH8ZBmFXXg7AhK828uH29te5ExE5Ewo7IuJ3xuVm2Mkt2sR7m4v9XI2IBBuFHRHxv8suA2B0yZd8tmkv1a56PxckIsFEYUdE/C81FU9GBlaPm/P2b2fFzsP+rkhEgojCjoj0CkbjennjDu7gox0atyMi3UdhR0R6h8awk31wOx/vKMPtDrmHu4uIjyjsiEjvkJsLwNjiQo5Wn+CLA5V+LkhEgoXCjoj0DqNHQ3Q0sa4ahh0p4mPdyhKRbqKwIyK9Q1gYTJgAQPbBHXyosCMi3URhR0R6j6ZByod2sPWQkxLHCT8XJCLBQGFHRHqPxrAzqbQQgI8L1bsjImdPYUdEeo+cHAAGlu2nX42D/+w64ueCRCQYKOyISO8RHw8jRwLmKuif7j6iKegictYUdkSkd5k4EYALD++msqaObcVOPxckIoFOYUdEepfsbAAmOfYB8J8vdStLRM6Owo6I9C6NYWfYAXOQssbtiMjZUtgRkd7lggvAYiGq/DCJVeWs3XuUE3UN/q5KRAKYwo6I9C5RUZCZCcBk5z5q692s31fh56JEJJAp7IhI79N4KyvvxEFA43ZE5Owo7IhI79MYdsaU7gZglcKOiJwFhR0R6X0aw07iri0AbDnk5Jir3p8ViUgAU9gRkd6ncZCytaSEMZYaGtweNhRV+rsqEQlQCjsi0vtER3ufpHxtfTEAa/eU+7MiEQlgCjsi0js13srKqdwLwJo9R/1YjIgEMoUdEemdxo0DIOPALgA27q/EVa/n7YhI1ynsiEjvNGYMANE7t5MQE4Gr3s3mAw4/FyUigUhhR0R6p9GjATB27+ZryZGAbmWJyJnpkbCzYMECMjIyiIyMJDs7m5UrV7bbdvHixRiG0Wo7ceLEGV9TRAJQ//6QkgLAlQ2HAVirsCMiZ8DnYWfJkiXMmjWLOXPmsGHDBiZPnszUqVMpKipq95y4uDiKi4tbbJGRkWd1TREJQI23ssY59gOwfl8FDW6PPysSkQDk87Dz1FNPcfvtt3PHHXeQmZnJ008/TWpqKs8880y75xiGQXJycovtbK7pcrlwOp0tNhEJAI1hJ3nfTmJsYVS76tlVVuXnokQk0Pg07NTW1rJ+/Xry8vJa7M/Ly2PVqlXtnlddXU1aWhqDBg1i2rRpbNiw4ayuOX/+fOx2u3dLTU09i18lIj2mcdyOZfNmzk+1A/D5Pj1cUES6xqdh58iRIzQ0NJCUlNRif1JSEiUlJW2eM3LkSBYvXsxbb73Fq6++SmRkJBdddBG7du0642vOnj0bh8Ph3fbv398Nv05EfK6xZ4fNmxk7qC8AG4q0ArqIdE1YT3yJYRgtPns8nlb7muTk5JCTk+P9fNFFFzFu3Dh+97vf8dvf/vaMrmmz2bDZbGdavoj4y8iREBYGlZXk2mr4PfC5wo6IdJFPe3YSEhKwWq2telzKyspa9cy0x2KxcOGFF3p7drrjmiISIGw277IRY46aExB2Hz5GZU2tP6sSkQDj07ATERFBdnY2y5Yta7F/2bJlTJo0qVPX8Hg8bNy4kZTGKajdcU0RCSCNt7Jid24nIyEaMJ+mLCLSWT6/jXXvvfeSn5/P+PHjyc3N5dlnn6WoqIiZM2cCcMsttzBw4EDmz58PwMMPP0xOTg7Dhg3D6XTy29/+lo0bN/KHP/yh09cUkSDSOEiZTZsYe/1U9hw5xudFlVw6ItG/dYlIwPB52Jk+fTrl5eXMmzeP4uJisrKyeO+990hLSwOgqKgIi+VkB1NlZSV33nknJSUl2O12xo4dy4oVK5gwYUKnrykiQaRpkPKmTYz9UT/+vuGgBimLSJcYHo8n5J7Q5XQ6sdvtOBwO4uLi/F2OiHSkqAjS0iAsjC27ipm2cA2xkWF88b95WCxtT0oQkeB0pn+/tTaWiPRuqakQEwP19YysLqVPuJWqE/V8ebja35WJSIBQ2BGR3s0wIDMTgLCdhYweZD5c8AsNUhaRTlLYEZHerzHssG0bYwaaYWfzQYcfCxKRQKKwIyK936hR5uu2bd6enU0HFHZEpHMUdkSk92sWdsY0LhuxrdhJXYPbj0WJSKBQ2BGR3q8p7OzYQXpfG7GRYdTWuyks0QroInJ6Cjsi0vulp0NkJLhcGPv2MWaQxu2ISOcp7IhI72e1etfIYts2Rg80b2Vp3I6IdIbCjogEhuYzsrw9O5p+LiKnp7AjIoGh+YysxunnhSVVnKhr8GNRIhIIFHZEJDA0CzuD+vUhPjqCugYPOzRIWUROQ2FHRAJDU9jZvh0DvL07mw7oVpaIdExhR0QCw7nnQng4HDsG+/d7w862Q04/FyYivZ3CjogEhvBwM/AA7NzJqAHmisfbihV2RKRjCjsiEjiGDzdfd+5kVIoZdnaUVFGvJymLSAcUdkQkcDQLO4Pjo4iOsFJb7+arI8f8W5eI9GoKOyISOJqFHYvFILOxd2frIT1cUETap7AjIoGjKewUFgKcHLejQcoi0gGFHREJHE1hZ+9ecLk4T4OURaQTFHZEJHAkJ0NMDLjd8NVXjEo5Of3c4/H4uTgR6a0UdkQkcBgGjBhhvt+5k2FJMVgtBhU1dZQ4T/i3NhHptRR2RCSwNBukHBluZWj/GEDjdkSkfQo7IhJYmoUd0CBlETk9hR0RCSynhp0UDVIWkY4p7IhIYGlv+rnCjoi0Q2FHRALLsGHma2kpOBzeBwvuK6+h6kSdHwsTkd5KYUdEAovdDklJ5vtdu4iPjiDFHgmY62SJiJxKYUdEAk+z6efQbNyOBimLSBsUdkQk8GhGloh0gcKOiASeU8JO07IRW4u1IKiItKawIyKB55SwMzLZDDs7S6upb3D7qyoR6aV6JOwsWLCAjIwMIiMjyc7OZuXKle22fe6555g8eTL9+vWjX79+XHnllaxdu7ZFm9tuuw3DMFpsOTk5vv4ZItJbNA87Hg+D46PoE26ltt7NvqM1/q1NRHodn4edJUuWMGvWLObMmcOGDRuYPHkyU6dOpaioqM32y5cv57vf/S4ff/wxBQUFDB48mLy8PA4ePNii3ZQpUyguLvZu7733nq9/ioj0FkOGgMUCVVVQUoLFYjA8yVw2olAzskTkFD4PO0899RS33347d9xxB5mZmTz99NOkpqbyzDPPtNn+lVde4Qc/+AEXXHABI0eO5LnnnsPtdvPhhx+2aGez2UhOTvZu8fHx7dbgcrlwOp0tNhEJYDYbpKeb7xtvZQ1PigU0/VxEWvNp2KmtrWX9+vXk5eW12J+Xl8eqVas6dY2amhrq6upahZnly5eTmJjI8OHDmTFjBmVlZe1eY/78+djtdu+Wmpra9R8jIr3LKeN2RiSbYWenwo6InMKnYefIkSM0NDSQ1PQAsEZJSUmUlJR06ho/+9nPGDhwIFdeeaV339SpU3nllVf46KOPePLJJ1m3bh2XX345LperzWvMnj0bh8Ph3fbv33/mP0pEeod2BikXlirsiEhLYT3xJYZhtPjs8Xha7WvLE088wauvvsry5cuJjIz07p8+fbr3fVZWFuPHjyctLY13332XG264odV1bDYbNpvtLH6BiPQ6Q4ear7t3Ayd7dvaWH+N4bQN9Iqz+qkxEehmf9uwkJCRgtVpb9eKUlZW16u051a9//Wsee+wx3n//fcaMGdNh25SUFNLS0ti1a9dZ1ywiAeLcc83XxrCTEBNBfHQEHg98WVbtx8JEpLfxadiJiIggOzubZcuWtdi/bNkyJk2a1O55v/rVr3jkkUdYunQp48ePP+33lJeXs3//flJSUs66ZhEJEM17dhp7i0d4BylrEoKInOTz2Vj33nsvzz//PC+88ALbt2/nxz/+MUVFRcycOROAW265hdmzZ3vbP/HEEzz44IO88MILpKenU1JSQklJCdXV5v9Sq66u5v7776egoIC9e/eyfPlyrr32WhISEvjGN77h658jIr1Fero5/fzYMXMFdE7eytL0cxFpzudjdqZPn055eTnz5s2juLiYrKws3nvvPdLS0gAoKirCYjmZuRYsWEBtbS3f+ta3WlznoYceYu7cuVitVjZv3sxLL71EZWUlKSkpXHbZZSxZsoTY2Fhf/xwR6S0iImDwYNi71+zdSU4+GXY0SFlEmjE8Ho/H30X0NKfTid1ux+FwEBcX5+9yRORMXXklfPghLF4Mt97K50UV3LBgFYmxNtbOufK0p4tIYDnTv99aG0tEAtcpg5SbHixYVuWi4litv6oSkV5GYUdEAlfTIOUvvwQgxhZGanwfQLeyROQkhR0RCVyn9OwA3hlZGqQsIk0UdkQkcJ3SswMnZ2RpjSwRaaKwIyKBa8gQ8/XoUaioAGBE47IRO3UbS0QaKeyISOCKiYHkZPN907IRSScXBA3ByaYi0gaFHREJbKeM2xnSP5pwq0GVq56Dlcf9WJiI9BYKOyIS2E4ZtxNutXBu/xhAt7JExKSwIyKBrY0ZWcMab2XtKtWCoCKisCMiga75gqCNhiU29ewo7IiIwo6IBLqmnp1m08+HJ5lh58sy3cYSEYUdEQl0TT07hw5BTY25K7HxNlZZNW63ZmSJhDqFHREJbPHx0Lev+f6rrwBIPyeKcKtBTW0DhxyakSUS6hR2RCTwnTJIOcxqYUiCeStLg5RFRGFHRAJfG8tGDG0ct7NL43ZEQp7CjogEvjamnw9vHLejGVkiorAjIoGvjZ6dYd6eHYUdkVCnsCMiga+tnp2m6eelWiNLJNQp7IhI4GsKO/v2QV0dAGnnmGtkHatt4JDjhB+LExF/U9gRkcCXkgKRkdDQAEVFgLlGVkZCNKA1skRCncKOiAQ+iwWGDDHft1g2whyk/KUGKYuENIUdEQkObS4IqtXPRURhR0SCRVthp9myESISuhR2RCQ4NIWdxiUjoPmCoNWakSUSwhR2RCQ4tDFmJ+2caMIsBtWueoo1I0skZCnsiEhwaH4bq7EXJyJMM7JERGFHRIJFejoYBhw7BmVl3t3Dmt3KEpHQpLAjIsHBZoPUVPN9G4OU1bMjEroUdkQkeHQw/VwzskRCl8KOiASPNtfIOvlgQc3IEglNCjsiEjzamH6e3jgjq8pVT4lTM7JEQlGPhJ0FCxaQkZFBZGQk2dnZrFy5ssP2b7zxBqNGjcJmszFq1CjefPPNFsc9Hg9z585lwIAB9OnTh0svvZStW7f68ieISCBoY/p5RJiFdO+MLN3KEglFPg87S5YsYdasWcyZM4cNGzYwefJkpk6dSlGneYhAAAAgAElEQVTjYn2nKigoYPr06eTn5/PFF1+Qn5/Pd77zHdasWeNt88QTT/DUU0/x+9//nnXr1pGcnMx//dd/UVWlAYgiIa2N21gAwxIbx+1okLJISDI8Pr6JPXHiRMaNG8czzzzj3ZeZmcn111/P/PnzW7WfPn06TqeTf/3rX959U6ZMoV+/frz66qt4PB4GDBjArFmzeOCBBwBwuVwkJSXx+OOPc9ddd7W6psvlwuVyeT87nU5SU1NxOBzExcV1588VEX+qrIR+/cz3VVUQY4acp5bt5Lcf7mL6+FQe/9YYPxYoImfD6XRit9u7/Pfbpz07tbW1rF+/nry8vBb78/LyWLVqVZvnFBQUtGp/1VVXedvv2bOHkpKSFm1sNhuXXHJJu9ecP38+drvdu6U2TU8VkeDSty/Ex5vvm43b8fbslKlnRyQU+TTsHDlyhIaGBpKSklrsT0pKoqSkpM1zSkpKOmzf9NqVa86ePRuHw+Hd9u/ff0a/R0QCQAczsnZpRpZISOqRAcqGYbT47PF4Wu3ravuuXNNmsxEXF9diE5Eg1UbYSU+Iwto4I6vU6WrnRBEJVj4NOwkJCVit1lY9LmVlZa16ZpokJyd32D45ORmgS9cUkRDSRtixhVlJPycK0JOURUKRT8NOREQE2dnZLFu2rMX+ZcuWMWnSpDbPyc3NbdX+/fff97bPyMggOTm5RZva2lo++eSTdq8pIiGkjWftwMllI/QkZZHQE+brL7j33nvJz89n/Pjx5Obm8uyzz1JUVMTMmTMBuOWWWxg4cKB3ZtaPfvQjLr74Yh5//HG+/vWv889//pMPPviA//znP4B5+2rWrFk89thjDBs2jGHDhvHYY48RFRXFTTfd5OufIyK9XRvP2gEYnhTD0q2afi4SinwedqZPn055eTnz5s2juLiYrKws3nvvPdLS0gAoKirCYjnZwTRp0iRee+01HnzwQX7xi19w7rnnsmTJEiZOnOht89Of/pTjx4/zgx/8gIqKCiZOnMj7779PbGysr3+OiPR2TT07+/ZBfT2Emf/MDU1Sz45IqPL5c3Z6ozOdpy8iAcDthqgocLnM3p3Gnp4dJU6mPL2S2MgwNj2U1+EkCRHpnXrlc3ZERHqcxdLmrayMhGgsBlSdqKesSjOyREKJwo6IBJ92Z2Q1rZGlcTsioURhR0SCT3trZCU1rZGlcTsioURhR0SCz2mnn6tnRySUKOyISPBpZ/q5enZEQpPCjogEn+a3sZpNOG3q2dlZWqU1skRCiMKOiASfjAwwDKiuhsOHvbuH9DdnZDk1I0skpCjsiEjwsdlg0CDzfbNbWZHhmpElEooUdkQkOLUzI2t4UtOtLI3bEQkVCjsiEpzaDTvmIOWdJerZEQkVCjsiEpzaCzvJjT07mn4uEjIUdkQkODVNPz/lWTtNt7F2lVZrRpZIiFDYEZHg1E7PTvo50YRZDKpd9RxynPBDYSLS0xR2RCQ4NYWdkhI4dsy7OyLMwpD+mpElEkoUdkQkOPXrZ27QetmIphlZGqQsEhIUdkQkeLVzK2uEpp+LhBSFHREJXqebfq7bWCIhQWFHRILXaR4s+GVZNW63ZmSJBDuFHREJXu2EnbRzookIs3C8roEDFcf9UJiI9CSFHREJXu08a8dqMTi3v3krq1C3skSCnsKOiASvpp6dvXuhvr7FoREatyMSMhR2RCR4DRxoroBeXw/797c45J1+rrAjEvQUdkQkeFkskJFhvtf0c5GQpbAjIsHtNDOydh+upr7B3dNViUgPUtgRkeDWTtgZ1K8PfcKt1Na72Xe0xg+FiUhPUdgRkeDWTtixWAyGNQ5S3qVxOyJBTWFHRIJbO9PPAYYlmreyCks0bkckmCnsiEhwa96z42n5tOQRyY3Tz8vUsyMSzBR2RCS4ZWSAYUBVFRw50uJQ0/Rz3cYSCW4KOyIS3CIjzeftQLszsr46fIzaes3IEglWCjsiEvzaGaQ8wB5JjC2MereHveXH/FCYiPQEn4adiooK8vPzsdvt2O128vPzqaysbLf90aNH+eEPf8iIESOIiopi8ODB3HPPPTgcjhbtDMNotS1cuNCXP0VEAlk7YccwTs7IKizRrSyRYBXmy4vfdNNNHDhwgKVLlwJw5513kp+fz9tvv91m+0OHDnHo0CF+/etfM2rUKPbt28fMmTM5dOgQf/vb31q0ffHFF5kyZYr3s91u990PEZHA1k7YAfNJyhuKKrVshEgQ81nY2b59O0uXLmX16tVMnDgRgOeee47c3FwKCwsZMWJEq3OysrJ44403vJ/PPfdcHn30UW6++Wbq6+sJCztZbt++fUlOTu5ULS6XC5fL5f3sdDrP9GeJSCDqKOwkm+N2dqhnRyRo+ew2VkFBAXa73Rt0AHJycrDb7axatarT13E4HMTFxbUIOgB33303CQkJXHjhhSxcuBC3u/3BhfPnz/feSrPb7aSmpnb9B4lI4Bo2zHzdtavVoZHJcQBsL9b/CBIJVj4LOyUlJSQmJrban5iYSElJSaeuUV5eziOPPMJdd93VYv8jjzzC66+/zgcffMCNN97Ifffdx2OPPdbudWbPno3D4fBu+09Z/VhEglxT2Ckrg1PGDWammD07ByqO4zxR19OViUgP6HLYmTt3bpsDhJtvn332GWAO/juVx+Npc/+pnE4n11xzDaNGjeKhhx5qcezBBx8kNzeXCy64gPvuu4958+bxq1/9qt1r2Ww24uLiWmwiEkJiYyElxXy/c2eLQ32jIkixRwIapCwSrLo8Zufuu+/mxhtv7LBNeno6mzZtorS0tNWxw4cPk5SU1OH5VVVVTJkyhZiYGN58803Cw8M7bJ+Tk4PT6aS0tPS01xaREDV8OBQXm2FnwoQWhzJT4ih2nGBHsZML0+P9VKCI+EqXw05CQgIJCQmnbZebm4vD4WDt2rVMaPyHZc2aNTgcDiZNmtTueU6nk6uuugqbzcZbb71FZGTkab9rw4YNREZG0rdv387/EBEJLcOHwyeftOrZARiZHMtHO8rYrp4dkaDks9lYmZmZTJkyhRkzZvDHP/4RMKeeT5s2zTsT6+DBg1xxxRW89NJLTJgwgaqqKvLy8qipqeHPf/4zTqfTO3Oqf//+WK1W3n77bUpKSsjNzaVPnz58/PHHzJkzhzvvvBObzearnyMiga5pBmhbYSdFg5RFgplPn7PzyiuvcM8995CXlwfAddddx+9//3vv8bq6OgoLC6mpqQFg/fr1rFmzBoChQ4e2uNaePXtIT08nPDycBQsWcO+99+J2uxkyZAjz5s3jf/7nf3z5U0Qk0A0fbr62EXZGpTStfl6F2+3BYjn9uEIRCRyGx3PKMsAhwOl0YrfbvdPaRSQEFBbCyJEQHW0uCtpsokR9g5tRD/2b2no3y++/lPSEaD8WKiLtOdO/31obS0RCQ0YGWK1w7Jg5ULmZMKuFEUlNDxfUrSyRYKOwIyKhISLCDDzQ7iBlgG3FGqQsEmwUdkQkdDSN2yksbHUos3GQ8g4NUhYJOgo7IhI6OhikPDJFa2SJBCuFHREJHR2EnczGNbKKjtZQpWUjRIKKwo6IhI4OnrXTLzqC5DjzIaY7S9W7IxJMFHZEJHQ09ex89RXUte69abqVpUHKIsFFYUdEQseAARAVBfX1sHdvq8MapCwSnBR2RCR0WCwwbJj5voPp51o2QiS4KOyISGjpYPr5qKaenZIqGtwh93B5kaClsCMioaWDGVlD+sfQJ9xKTW0De45U93BhIuIrCjsiElo6CDtWi8GoAWbvzpaDupUlEiwUdkQktHQQdgBGD7QDsPmgo6cqEhEfU9gRkdDSFHYOHoTq1reqshR2RIKOwo6IhJb4eEhIMN+30buTNdC8jbXtkBO3BimLBAWFHREJPZmZ5uv27a0ODe0fQ2S4hWpXPXvKj/VwYSLiCwo7IhJ6Ro0yX7dta3UozGrxPlxwi25liQQFhR0RCT0dhB1oNkj5gMKOSDBQ2BGR0HOasNM0SHnLIYUdkWCgsCMioacp7OzeDS5Xq8NNPTtbD2qQskgwUNgRkdCTkgJxcdDQALt2tTo8LDEGW5iFKlc9+47W+KFAEelOCjsiEnoMo9ODlPW8HZHAp7AjIqHptON2NCNLJFgo7IhIaNKMLJGQobAjIqGpCzOyPB4NUhYJZAo7IhKamsLOzp1QX9/q8PCkWCLCLFSdqGdfuQYpiwQyhR0RCU2pqRAVBXV15hT0U4Q3G6S8SeN2RAKawo6IhCaL5eQaWVu3ttnkgkHmrayNRZU9VZWI+IDCjoiErvPOM1/bCTtjB/cDYMP+ip6qSER8QGFHRELX6NHm66ZNbR4eO7gvYD5J2VXf0FNViUg3U9gRkdA1Zoz5unlzm4cHx0cRHx1BbYObbYecPViYiHQnn4adiooK8vPzsdvt2O128vPzqazs+N73pZdeimEYLbYbb7zxrK8rItJKU9jZtQtqWs+4MgyDsalm784GjdsRCVg+DTs33XQTGzduZOnSpSxdupSNGzeSn59/2vNmzJhBcXGxd/vjH//YLdcVEWkhKQn69we3u93n7TTdytqwX2FHJFCF+erC27dvZ+nSpaxevZqJEycC8Nxzz5Gbm0thYSEjRoxo99yoqCiSk5O77boulwtXs5WNnU51R4sI5hpZY8bAhx+a43bGj2/VxDtIuUiDlEUClc96dgoKCrDb7d5AApCTk4PdbmfVqlUdnvvKK6+QkJDAeeedx/33309VVdVZXXf+/PneW152u53U1NSz/HUiEjSabmW1M0h5zCA7hgEHKo5TVnWiBwsTke7is7BTUlJCYmJiq/2JiYmUlJS0e973vvc9Xn31VZYvX84vfvEL3njjDW644Yazuu7s2bNxOBzebf/+/Wfwi0QkKJ1mRlZsZDjDE2MBPW9HJFB1+TbW3Llzefjhhztss27dOsAc3Hcqj8fT5v4mM2bM8L7Pyspi2LBhjB8/ns8//5xx48ad0XVtNhs2m63DmkUkRDXv2fF4zFtbpxg7uC+FpVV8XlRJ3nlt32IXkd6ry2Hn7rvvbjU76lTp6els2rSJ0tLSVscOHz5MUlJSp79v3LhxhIeHs2vXLsaNG0dycnK3XFdEBDDXyLJYoLwcSkogJaVVk3Fp/Xht3X4+23vUDwWKyNnqcthJSEggISHhtO1yc3NxOBysXbuWCRMmALBmzRocDgeTJk3q9Pdt3bqVuro6Uhr/Aequ64qIANCnDwwfDjt2mL07bYSdCenxAGw64OBEXQOR4daerlJEzoLPxuxkZmYyZcoUZsyYwerVq1m9ejUzZsxg2rRp3hlTBw8eZOTIkaxduxaA3bt3M2/ePD777DP27t3Le++9x7e//W3Gjh3LRRdd1Onrioh0yWkGKaedE0X/WBu1DW6+0BR0kYDj0+fsvPLKK4wePZq8vDzy8vIYM2YML7/8svd4XV0dhYWF1DQ+zCsiIoIPP/yQq666ihEjRnDPPfeQl5fHBx98gNVq7fR1RUS6pGmQ8hdftHnYMAxv78463coSCTiGx+Px+LuInuZ0OrHb7TgcDuLi4vxdjoj427vvwrRp5iro7TxccPGne5j79jYmD0vg5dsnttlGRHzrTP9+a20sEZHsbPN1xw6orm6zyYUZZs/O5/sqqG9w91RlItINFHZERJKTYcAAc+r5xo1tNhmZHEesLYxjtQ1sL65qs42I9E4KOyIicLJ3Z/36Ng9bLQbZ6ebSEWs1bkckoCjsiIjAacMOwIVNg5T3KOyIBBKFHRER6FTYmdg4bmfNnnLc7pCb2yESsBR2RESg5SDlY8fabHJ+al+iI6xU1NSxrdjZg8WJyNlQ2BERAfPJySkp4Ha3O0g53GphQmPvzqrdR3qyOhE5Cwo7IiJNOnEr66Kh5nI5n35Z3hMViUg3UNgREWnShbCzds9Rauv1vB2RQKCwIyLSpBNhZ0RSLOdER3C8roENRRU9VJiInA2FHRGRJuPHm6/btoGz7QHIFovBpKZbWbt1K0skECjsiIg0SUmB9HTzScpr17bb7KJzzwHg0y81SFkkECjsiIg0N2mS+bpqVbtNmsbtbNxfSdWJup6oSkTOgsKOiEhzubnma0FBu01S46PISIimwe3hP7vUuyPS2ynsiIg019SzU1BgPnOnHZeNSATgox1lPVGViJwFhR0RkebGjIGoKHA4zKcpt+PykWbY+bjwsJaOEOnlFHZERJoLC4MJE8z3HYzbmZART3SElSPVLrYccvRQcSJyJhR2RERO1fxWVjsiwix8bZg5UPnjHYd7oioROUMKOyIip2oapNxBzw6cvJX1UaHG7Yj0Zgo7IiKnyskxX3fsgKNH223WNEh504FKjlS7eqIyETkDCjsiIqdKSICRI833K1e22ywxLpKsgXF4PPDh9tIeKk5EukphR0SkLZddZr5+9FGHzaaclwzAv7aU+LoiETlDCjsiIm25/HLz9XRhJysFMJeOcBzX05RFeiOFHRGRtlx6qfm6ZQuUtT8AeWhiDMOTYqhr8OhWlkgvpbAjItKWhATzAYMAy5d32HRqY+/Oe5t1K0ukN1LYERFpTydvZV092gw7K3Yd1sKgIr2Qwo6ISHs6GXaGJ8UwpH80tfVurZUl0gsp7IiItOfii8FigV274MCBdpsZhsHVjbey3tp4qKeqE5FOUtgREWmP3Q7Z2eb7Dz/ssOn1YwcA8MnOw3rAoEgvo7AjItKRvDzz9V//6rDZ0MRYzh9kp97tUe+OSC+jsCMi0pFp08zXpUuhruPBxzeMGwTA3ze0f8tLRHqeT8NORUUF+fn52O127HY7+fn5VFZWttt+7969GIbR5vb6669727V1fOHChb78KSISqi68EPr3B4cDPv20w6bXnj+AcKvBloNOCkuqeqhAETkdn4adm266iY0bN7J06VKWLl3Kxo0byc/Pb7d9amoqxcXFLbaHH36Y6Ohopk6d2qLtiy++2KLdrbfe6sufIiKhymqFq68237/zTodN46MjvIuD/v1z9e6I9BY+Czvbt29n6dKlPP/88+Tm5pKbm8tzzz3HO++8Q2FhYZvnWK1WkpOTW2xvvvkm06dPJyYmpkXbvn37tmjXp0+fdmtxuVw4nc4Wm4hIpzXdyjpN2AH4ZrZ5K+uNzw9SW+/2ZVUi0kk+CzsFBQXY7XYmTpzo3ZeTk4PdbmfVqlWdusb69evZuHEjt99+e6tjd999NwkJCVx44YUsXLgQt7v9f1Tmz5/vvZVmt9tJTU3t+g8SkdCVlwdhYVBYaE5D78DlIxNJirNxpNrF0q16orJIb+CzsFNSUkJiYmKr/YmJiZSUdO4fgEWLFpGZmcmkSZNa7H/kkUd4/fXX+eCDD7jxxhu57777eOyxx9q9zuzZs3E4HN5t//79XfsxIhLa4uLgkkvM9+++22HTcKuF704YDMDLBXt9W5eIdEqXw87cuXPbHUTctH322WeAOZD4VB6Pp839pzp+/Dh/+ctf2uzVefDBB8nNzeWCCy7gvvvuY968efzqV79q91o2m424uLgWm4hIlzTdynrzzdM2/e6EwYRZDNbtrWB7sW6bi/hbl8PO3Xffzfbt2zvcsrKySE5OprS09QrAhw8fJikp6bTf87e//Y2amhpuueWW07bNycnB6XS2+X0iIt3ihhvM15Ur4eDBDpsmxUVy1XnJALxUsM/XlYnIaYR19YSEhAQSEhJO2y43NxeHw8HatWuZMGECAGvWrMHhcLS6LdWWRYsWcd1119G/f//Ttt2wYQORkZH07dv39D9ARORMDB4MF11kTj9//XWYNavD5vm5aby7uZh/bDjIA1NG0DcqoocKFZFT+WzMTmZmJlOmTGHGjBmsXr2a1atXM2PGDKZNm8aIESMAOHjwICNHjmTt2rUtzv3yyy9ZsWIFd9xxR6vrvv322zz33HNs2bKF3bt38/zzzzNnzhzuvPNObDabr36OiAh897vm66uvnrbpxIx4RibHcryugT+tUu+OiD/59Dk7r7zyCqNHjyYvL4+8vDzGjBnDyy+/7D1eV1dHYWEhNTU1Lc574YUXGDhwIHlNj2lvJjw8nAULFpCbm8uYMWP4zW9+w7x583jyySd9+VNEROBb3zIXBl27Fr76qsOmhmHwg8uGAvDiqj0cc9X3RIUi0gbD4/F4/F1ET3M6ndjtdhwOhwYri0jX/Nd/wQcfwGOPwezZHTatb3BzxVOfsK+8hgevyeSOyUN6qEiR4HSmf7+1NpaISFfceKP5+pe/wGn+t2KY1cLMS84F4PmVe3DVN/i6OhFpg8KOiEhX3HADREbCli3m7azTNR83kKQ4GyXOE/x1nZ7xJeIPCjsiIl3Rrx985zvm+z/+8bTNbWFW/qdx7M5vPvxSY3dE/EBhR0Skq2bONF9few0qKk7b/MYLB5N2ThRHql08t7Ljgc0i0v0UdkREuionB8aMgePHodkM0/ZEhFn4yVXmIzeeW/EVh6tcvq5QRJpR2BER6SrDgLvuMt8vXHjagcoAV2elMGaQnWO1DTy1bKePCxSR5hR2RETOxM03Q3Q0bN9uTkU/DYvF4MFrRgHw6toi1u87/e0vEekeCjsiImciLg6aFiqeP79Tp0zIiOdb2YMAmPPmZuoa3L6qTkSaUdgRETlT998P4eHw8cdQUNCpU35+dSZ9o8LZUVLFi5/u8XGBIgIKOyIiZy41FW65xXz/2GOdOiU+OoKfT80E4Mn3d7KztMpX1YlII4UdEZGz8cAD5npZ77wDX3zRqVO+PX4Qlwzvj6vezT2vbtCTlUV8TGFHRORsDBsG06eb73/6007NzDIMg199ewzx0RHsKKniV0sLfVykSGhT2BEROVuPPAIREfD++/Cvf3XqlMTYSJ745hgAnv/PHpZuKfZlhSIhTWFHRORsnXsu/OhH5vt774W6uk6dduWoJP77ogzztL9+wY4Sp68qFAlpCjsiIt1hzhzo3x8KC2HBgk6f9vOrR/K1oQnU1DZwx58+o7xaT1cW6W4KOyIi3cFuh1/+0nw/Zw581bk1sMKsFn5/01gGx0dxoOI4t764FueJzvUMiUjnKOyIiHSXO+6ASy6BY8fg+98Hd+ceGtg3KoLF37+Qc6Ij2HLQyR2LP+N4rWZoiXQXhR0Rke5iscCLL5rLSKxYAb/5TadPHdI/hj/99wRiI8NYu/co31+8lir18Ih0C4UdEZHulJEBTz5pvv/Zzzr9ZGWArIF2XrztQmJsYaz+6ig3PbdGY3hEuoHCjohId7vzTrjhBqitNV8PHuz0qePT43l1Rg7x0RFsPujgm8+sYpeesixyVhR2RES6m2HAn/4EWVlQUgLXX2+O4+mk0YPs/PWuXAb27cPe8hqu/8On/HtriQ8LFgluCjsiIr4QEwP//CfEx8Nnn8G110JNTadPH5oYw1t3X0TOkHiO1TZw18vr+d9/btHAZZEzoLAjIuIrQ4bAu++awefjj+Eb34ATJzp9+jkxNl6+fSK3f8188OBLBfuY9ruVrNt71FcViwQlhR0REV/KyTGXkIiKMpeTuOIKOHy406eHWy38YtooXvrvCSTG2th9+BjfXljAT17/giMavCzSKQo7IiK+9rWvmYGnb19YtQomToQtW7p0iYuH9+ffsy7mxgtTAXh9/QEmP/4xjy/dQWVNrS+qFgkahsfTiSV6g4zT6cRut+NwOIiLi/N3OSISKnbsgGnTYPdusNngiSfg7rvN5/N0wfp9R5n39ja+OOAAINYWxq2T0rk5J41ke6QvKhfpFc7077fCjsKOiPSkI0cgPx+WLjU/X3YZ/O53cN55XbqMx+Phw+1lPLlsJ9uLzQVErRaDq85LIj8nnYkZ8VgsRndXL+JXCjtdoLAjIn7l8cAzz8B995kDlq1WmDkTZs+GgQO7dCm328O/t5bw4qq9rN1zcuByij2Sa0ancO35AxgzyI5hKPhI4FPY6QKFHRHpFb76Cu6/H9580/wcEQG33go//CGMHt3ly+0ocfJSwT7e3niIKle9d39yXCQXD0/gkuGJXDT0HPpGRXTXLxDpUQo7XaCwIyK9ykcfwdy5sHLlyX0TJ8JNN8HXvw5paV263Im6Bj7ZeZi3vzjEh9vLOF7X8tk85/aPZtzgfoxL68e4wf04t380YVbNV5HeT2GnCxR2RKRXWrkSnn4a3noL6k/2zHD++XDNNeasrtxcc1ZXJ52oa2DtnqOs2HmYT3YeZldZdas2EVYL5ybGMCIphuHJsQztH0NqfBSp8VHE2MK645eJdIteGXYeffRR3n33XTZu3EhERASVlZWnPcfj8fDwww/z7LPPUlFRwcSJE/nDH/7Aec0G71VUVHDPPffw1ltvAXDdddfxu9/9jr6d/AdAYUdEerXSUnjlFfjHP+DTT8HtPnnMMGDUKBgzBjIzzfejRsG555q3wU7j6LFaNhRV8HlRBZ/vq+SLA5XUdPBU5n5R4Wbw6RdFsj2ShBgb/WNtJMRE0D/WRv8YG/HREeoZkh7RK8POQw89RN++fTlw4ACLFi3qVNh5/PHHefTRR1m8eDHDhw/nl7/8JStWrKCwsJDY2FgApk6dyoEDB3j22WcBuPPOO0lPT+ftt9/uVF0KOyISMA4fNp/C/PHH5jN6vvyy/baJiTBokLkNHAgJCeZyFc03ux2io82HHEZF4bZFcsBZS2FpFTtLqygsqWLPkWPsr6ihsqauUyUaBsTYwoiLDCc2Moy4PuHERYYT18fcFxcZRp+IMCLDLfQJt9InwkpkuNX7vk+4+Tky3EKE1UKY1UK41SDcaiHMYmC1GBpgLUAvDTtNFi9ezKxZs04bdjweDwMGDGDWrFk88MADALhcLpKSknj88ce566672L59O6NGjWL16tVMnDgRgNWrV5Obm8uOHTsYMWJEq+u6XC5crpNPGnU6naSmpirsiEjgKS2Fdetg27aWWxcWGm3FZvOGH2w2CA+HsKhwSFUAAAucSURBVDAarFZqDSsuLJzwGJzAiguD4x4Lxz0GxxugtsFDg2HgwQDDwAN4Gj97DPBg4DaajrXc72ne3ntO47WaMQzzCbiGxcBiGFgMzFcLWDD3GQZggPnS9Nl8Na/W+N44ebz1fryhqnkFHsNo8bm93OVtZZy6v/WJJ9+2fbFTv8Po6OCZ8nGADJ98Edk/mdmt1zzTsNOrbsbu2bOHkpIS8vLyvPtsNhuXXHIJq1at4q677qKgoAC73e4NOgA5OTnY7XZWrVrVZtiZP38+Dz/8cI/8BhERn0pKMh9MOG3ayX0eD5SXw4EDsH+/+XrwIBw9am4VFSffV1bC8ePm1sTlMreKihZfZQX6NG4iXbXG5YJuDjtnqleFnZKSEgCSkpJa7E9KSmLfvn3eNomJia3OTUxM9J5/qtmzZ3Pvvfd6Pzf17IiIBAXDMG9ZJSTABRd07hy323zGT03Nye3YMaitNQdH19WZrx29r683g5bbbb6eup3hfrfHg9vtwe0Bt8dDg7tpHzQ0O97QuM/t8Zg9RN7Xtt/j8ZyyHzx4Gl9Padv4ubmTR6DprQcwGq/Rsm3Ld23dQ+nMOYanxbd2nx6Ym2SblOPz7+isLoeduXPnnraXZN26dYwfP/6Mizr13qzH42mxr617t6e2ac5ms2Gz2c64HhGRoGOxnLx11ctY0MKN0r26HHbuvvtubrzxxg7bpKenn1ExycnJgNl7k5KS4t1fVlbm7e1JTk6mtLS01bmHDx9u1SMkIiIi0uWwk5CQQEJCgi9qISMjg+TkZJYtW8bYsWMBqK2t5ZNPPuHxxx8HIDc3F4fDwdq1a5kwYQIAa9asweFwMGnSJJ/UJSIiIoHLpz2FRUVFbNy4kaKiIhoaGti4cSMbN26kuvrkQ61GjhzJm42PSjcMg1mzZvHYY4/x5ptvsmXLFm677TaioqK46aabAMjMzGTKlCnMmDGD1atXs3r1ambMmMG0adPaHJwsIiIioe3/b+9+Q5ra/ziAf+bcXJlXLNOmDUsf+IcKapKZihaiUWg9EqlEouiPi9QeiCD+oXJYRvQgNVo+EBKFlCJilEZqpkUkGwSKZWZ/MCmhcmSl08/vwcVxl3qvZ7+74z3H9wv2wON39d7b2fm0nXPm1gOUS0pKqK6uzvH1zKs1bW1tlJSURERE/f399O3bN8eagoIC+vHjB+Xk5DguKtjS0uK4xg4RUX19PZ06dcpx1lZ6ejpduXLFnQ8FAAAAJAofF4Hr7AAAAEiCq/tvHPAOAAAAsoZhBwAAAGQNww4AAADIGoYdAAAAkDUMOwAAACBrGHYAAABA1jDsAAAAgKxh2AEAAABZc+sVlP+rZq6jODY2tshJAAAAYKFm9ttCr4e8JIcdm81GREQ6nW6RkwAAAIBQNpuNfH19F7x+SX5cxPT0NA0PD5OPjw8pFIp/9c8eGxsjnU5H79+/x0dRuBF6Fgd6Fgd6Fgd6Fo+7umZmstlsFBQURB4eCz8SZ0m+suPh4UFr165169/xxx9/4JdJBOhZHOhZHOhZHOhZPO7oWsgrOjNwgDIAAADIGoYdAAAAkDVlWVlZ2WKHkBulUklJSUnk6bkk3yUUDXoWB3oWB3oWB3oWz3+p6yV5gDIAAAAsHXgbCwAAAGQNww4AAADIGoYdAAAAkDUMOwAAACBrGHYAAABA1jDsuKC6uprWr19PGo2G9Ho9dXZ2/u365uZmioqKIi8vL4qKiqJbt26JlFTahPRsMpkoISGB/Pz8yM/Pj5KTk+nZs2cippUuoc/nGY2NjaRQKGjfvn1uTigPQnv++vUrGQwG0mq1pNFoKDIyksxms0hppUtoz5cvX6bw8HBatmwZ6XQ6ys/Pp58/f4qUVpoePXpEaWlpFBQURAqFgm7fvv2P9+no6CC9Xk8ajYZCQ0Pp6tWrIiT9CwZBGhsbWaVSsclk4t7eXs7NzWVvb29++/btnOu7u7tZqVSy0Wjkvr4+NhqN7OnpyU+fPhU5ubQI7Xn//v1cVVXFFouF+/r6+NChQ+zr68sfPnwQObm0CO15xtDQEAcHB3NCQgLv3btXpLTSJbTnX79+cXR0NO/evZsfP37MQ0ND3NnZyVarVeTk0iK05xs3brCXlxfX19fzmzdv+P79+6zVajkvL0/k5NJiNpu5qKiIm5ubmYj41q1bf7t+cHCQly9fzrm5udzb28smk4lVKhU3NTWJlJgZw45AW7du5ePHjztti4iI4MLCwjnXZ2Rk8K5du5y2paamcmZmptsyyoHQnn9nt9vZx8eH6+rq3BFPNlzp2W63c1xcHF+/fp2zs7Mx7CyA0J5ramo4NDSUJyYmxIgnG0J7NhgMvHPnTqdtp0+f5vj4eLdllJuFDDsFBQUcERHhtO3YsWO8bds2d0ZzgrexBJiYmKCenh5KSUlx2p6SkkLd3d1z3ufJkyez1qemps67Hlzr+Xfj4+M0OTlJK1eudEdEWXC15zNnztDq1avp8OHD7o4oC670fOfOHYqNjSWDwUCBgYG0YcMGMhqNNDU1JUZkSXKl5/j4eOrp6XG85T04OEhms5n27Nnj9rxLyXz7wefPn9Pk5KQoGRb/Gs4SMjo6SlNTUxQYGOi0PTAwkEZGRua8z8jIiKD14FrPvyssLKTg4GBKTk52R0RZcKXnrq4uqq2tJavVKkZEWXCl58HBQXr48CEdOHCAzGYzvXr1igwGA9ntdiopKREjtuS40nNmZiZ9/vyZ4uPjiZnJbrfTiRMnqLCwUIzIS8Z8+0G73U6jo6Ok1WrdngHDjgsUCoXT18w8a9v/sx7+5GpvFy5coIaGBmpvbyeNRuOueLKx0J5tNhsdPHiQTCYT+fv7ixVPNoQ8n6enpykgIICuXbtGSqWS9Ho9DQ8PU2VlJYadfyCk5/b2diovL6fq6mqKiYmhgYEBys3NJa1WS8XFxWLEXTLm+rnMtd1dMOwI4O/vT0qlctb/Ej59+jRrap2xZs0aQevBtZ5nXLx4kYxGIz148IA2bdrkzpiSJ7Tn169f09DQEKWlpTm2TU9PExGRp6cn9ff3U1hYmHtDS5Arz2etVksqlYqUSqVjW2RkJI2MjNDExASp1Wq3ZpYiV3ouLi6mrKwsOnLkCBERbdy4kb5//05Hjx6loqIi8vDAkR7/hvn2g56enrRq1SpRMuAnKYBarSa9Xk+tra1O21tbW2n79u1z3ic2NnbW+paWlnnXg2s9ExFVVlbS2bNn6d69exQdHe3umJIntOeIiAh68eIFWa1Wxy09PZ127NhBVquVdDqdWNElxZXnc1xcHA0MDDiGSSKily9fklarxaAzD1d6Hh8fnzXQKJVK4j9P3nFb1qVmvv1gdHQ0qVQqcUKIdii0TMyc2lhbW8u9vb2cl5fH3t7ePDQ0xMzMWVlZTkf+d3V1sVKp5IqKCu7r6+OKigqcer4AQns+f/48q9Vqbmpq4o8fPzpuNpttsR6CJAjt+Xc4G2thhPb87t07XrFiBZ88eZL7+/v57t27HBAQwOfOnVushyAJQnsuLS1lHx8fbmho4MHBQW5paeGwsDDOyMhYrIcgCTabjS0WC1ssFiYivnTpElssFscp/oWFhZyVleVYP3PqeX5+Pvf29nJtbS1OPZeCqqoqDgkJYbVazVu2bOGOjg7H9xITEzk7O9tp/c2bNzk8PJxVKhVHRERwc3OzyImlSUjPISEhTESzbqWlpeIHlxihz+e/wrCzcEJ77u7u5piYGPby8uLQ0FAuLy9nu90ucmrpEdLz5OQkl5WVcVhYGGs0GtbpdJyTk8NfvnxZhOTS0dbWNue/tzPdZmdnc2JiotN92tvbefPmzaxWq3ndunVcU1MjamYFM16rAwAAAPnCMTsAAAAgaxh2AAAAQNYw7AAAAICsYdgBAAAAWcOwAwAAALKGYQcAAABkDcMOAAAAyBqGHQAAAJA1DDsAAAAgaxh2AAAAQNYw7AAAAICs/Q9TM4cgrSICRwAAAABJRU5ErkJggg==",
      "text/plain": [
       "PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x140a0b3d0>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "k= 1: 0.000041625814,\t 0.982912862001\n",
      "k= 2: 0.000011372225,\t 0.665794128298\n",
      "k= 3: 0.000006671248,\t 0.039292389065\n",
      "k= 4: 0.000000030047,\t 0.000111411420\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "1-element Array{PyCall.PyObject,1}:\n",
       " PyObject <matplotlib.lines.Line2D object at 0x140a51590>"
      ]
     },
     "execution_count": 73,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x=linspace(0,1,401)\n",
    "u=zeros(401)\n",
    "\n",
    "h=x[2]-x[1]\n",
    "ϵ=0.025\n",
    "\n",
    "#specify the boundary conditions\n",
    "u[1]=1.0\n",
    "u[end]=-1.0\n",
    "\n",
    "GaussSeidel_Burgers!(u,h,ϵ,4*4*250)\n",
    "\n",
    "#u[2:end-1]=2.0*(rand(length(u)-2)-0.5)\n",
    "\n",
    "plot(x, u)\n",
    "\n",
    "NewtonB!(u, FBurgers, JFBurgers, h, ϵ, nitrs=30, tol=1e-12)\n",
    "\n",
    "plot(x, u,\"r\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gripes about Newton's method\n",
    "\n",
    "* Sensitive to initial guesses\n",
    "* Cost for evaluation of the Jacobian matrix $f^\\prime(x)=(f_1(x),f_2(x),\\cdots,f_N(x))^\\top$ (we also used the notation $Jf(x)$): \n",
    "\n",
    "    For $x\\in \\mathbb{R}^N$, the cost for evaluating $f^\\prime(x)$ ranges from $O(N^2)$, when each component of $f$ depends only on a few components of $x$, to $O(N^3)$, when each component of $f$ depends on all the component of $x$.\n",
    "\n",
    "* Cost for evaluation of $(f^\\prime)^{-1})(x)$ is $O(N^3)$ in general, using Gaussian elimination based direct methods.\n",
    "* The need to do the above two items at every update"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A convergence theorem for Newton's method\n",
    "\n",
    "**Theorem**: Let $f:\\mathbb{R}^n\\mapsto \\mathbb{R}^n$ be a $C^1$ function. Let $x^*$ be a solution to $f(x)=0$. Let $||x||$ denotes a vector norm of $x$ and $||A||$ be the associated operator norm for the matrix $A$. \n",
    "\n",
    "Suppose $(f')^{-1}(x)$ exists and that there exist positive constants $C,L$ and $R$ such that\n",
    "\n",
    "$$ ||f'(x^*)||\\le C,\\,\\,\\,\\, ||f'(x)-f'(y)||\\le L||x-y||,\\,\\,\\,\\forall x,y\\in B_R(x^*). $$ \n",
    "\n",
    "There exists $r>0$ such that for any $x_0\\in B_r(x^*)$, the Newton's method yields\n",
    "\n",
    "$$||x^{(k+1)}-x^k||\\le CL\\cdot||x^(k)-x^*||^2.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Broyden's method -- a quasi-Newton method\n",
    "\n",
    "\n",
    "**Newton's method is too computationally costly because of the evaluation of the Jacobian matrix and inverting it at every step. We can make it cheaper by approximating the Jacobina matrix somehow**\n",
    "\n",
    "* A generalization of the secant method\n",
    "* **Low-rank update of the Jacobian matrix**\n",
    "\n",
    "Notation: $ J(x_n):=\\frac{\\partial f}{\\partial x}(x_n) $\n",
    "\n",
    "(Ignore Burden-Faires. See Wikipedia's description of [Broyden's method](https://en.wikipedia.org/wiki/Broyden%27s_method) )\n",
    "\n",
    "**IDEA 1**: (**Secant method**) Approximate $J(x_n)$ by a matrix $A_n$ along the direction $x_n-x_{n-1}$.\n",
    "\n",
    "\n",
    "$$y_n:=f(x_n)-f(x_{n-1})=:\\Delta f_n,$$\n",
    "\n",
    "$$s_n:=x_n-x_{n-1}=:\\Delta x_n,$$\n",
    "\n",
    "Find $A_n$ such that $$A_n\\Delta x_n = \\Delta f_n\\approx J(x_n)\\Delta x_n.$$\n",
    "\n",
    "**IDEA 2**: Relate $A_n$ to $A_{n-1}$, $x_n$, $x_{n-1}$, $f(x_n)$ and $f(x_{n-1})$.\n",
    "\n",
    "Solve\n",
    "\n",
    "$$ \\min_{A}||A-A_{n-1}||_F\\,\\,s.t.\\,\\, A\\Delta x_n = \\Delta f_n,$$\n",
    "\n",
    "and derive a **low-rank update** formula:\n",
    "\n",
    "$$A_n = A_{n-1}+\\frac{\\Delta f_n-A_{n-1}\\Delta x_n}{|\\Delta x_n|}\\left(\\frac{\\Delta x_n}{|\\Delta x_n|}\\right)^\\top.$$\n",
    "\n",
    "In other words, we have the form $A_n=A_{n-1}+uv^\\top$. **What's low-rank? $uv^\\top$ is a rank-1 matrix.**\n",
    "\n",
    "\n",
    "**IDEA 3**: Sherman-Morrison-Woodbury formula for:\n",
    "\n",
    "**updating $A_n^{-1}$ from $A_{n-1}^{-1}$, $x_n$, $x_{n-1}$, $f(x_n)$ and $f(x_{n-1})$**\n",
    "\n",
    "($B=A_{n-1}$)\n",
    "\n",
    "$$ (B+uv^\\top)^{-1} = \\left(I-wv^\\top\\right)B^{-1},\\,\\,\\, w:=\\frac{(B^{-1}u)}{1+v^\\top B^{-1} u}$$\n",
    "\n",
    "This means\n",
    "\n",
    "$$ A_n^{-1} := \\left(I +  \\frac{ \\Delta x - A_{n-1}^{-1}\\Delta f}{ \\Delta x^\\top A_{n-1}^{-1}\\Delta f  }\\Delta x^\\top\\right)A_{n-1}^{-1}$$\n",
    "\n",
    "(Ignore Burden-Faire. See GGK p155 for more in-depth exposition.)\n",
    "\n",
    "**Superlinear convergence**: $$\\lim_{n\\rightarrow\\infty} \\frac{e_n}{e_{n-1}} =0. $$\n",
    "\n",
    "Learning objectives:\n",
    "* Know the three ideas used in Broyden's method.\n",
    "* Know Sherman-Morrison-Woodbury formula and how to prove it.\n",
    "* Know the computational complexity of Broyden's method, in particular, compared to the straight forward Newton's method.\n",
    "* Know how fast Broyden's method converges."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Algorithm 1\n",
    "\n",
    "(Without low-rank update of the inverse of the \"Jacobian\" matrix $A_n$)\n",
    "\n",
    "$$A_n = A_{n-1}+\\frac{\\Delta f_n-A_{n-1}\\Delta x_n}{|\\Delta x_n|}\\left(\\frac{\\Delta x_n}{|\\Delta x_n|}\\right)^\\top.$$\n",
    "\n",
    "Inputs: $A_0, x_0, f(x)$\n",
    "\n",
    "For $n=0,1,2,\\cdots$\n",
    "1. Solve $A_n\\Delta x =-f(x_n)$\n",
    "2. $x_{n+1}:=x_n+\\Delta x$\n",
    "3. $\\Delta f=f(x_{n+1})-f(x_n)$\n",
    "4. Rank-1 update of the Jacobina matrix: $ A_{n+1}:=A_n+\\frac{\\left(\\Delta f-A_n\\Delta x\\right)\\Delta x}{\\Delta x^\\top \\Delta x}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "function Broyden_QuasiNewton!(x, F::Function, B::Array{Float64}; nitrs=10, tol=1e-12)\n",
    "   \n",
    "    for k in 1:nitrs\n",
    "        f=F(x)\n",
    "\n",
    "        #Solve BΔx=-F(x)  for Δx\n",
    "        Δx=-B\\f\n",
    "\n",
    "        #Update x\n",
    "        x.+=Δx\n",
    "        \n",
    "        #Evaluate F at new x\n",
    "        fnew=F(x)\n",
    "        \n",
    "        if norm(f)>tol\n",
    "            \n",
    "            Δf=fnew.-f\n",
    "            \n",
    "            #rank-1 update of B\n",
    "            B.+= (Δf-B*Δx)*Δx'./(Δx'Δx)\n",
    "        else\n",
    "           break \n",
    "        end\n",
    "        \n",
    "        @printf \"k=%2d: %1.12f,\\t %1.12f\\n\" k norm(f) norm(Δx)\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Algorithm 2\n",
    "\n",
    "(Using low rank update on the inverse of the \"Jacobian\" matrix $A_{n}^{-1}$)\n",
    "\n",
    "$$ A_n^{-1} := \\left(I +  \\frac{ \\Delta x - A_{n-1}^{-1}\\Delta f}{ \\Delta x^\\top A_{n-1}^{-1}\\Delta f  }\\Delta x^\\top\\right)A_{n-1}^{-1}$$\n",
    "\n",
    "Inputs: $A_0, x_0, f(x)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "function Broyden_QuasiNewton2!(x, F::Function, B::Array{Float64}; nitrs=10, tol=1e-12)\n",
    "   \n",
    "    Binv=B^(-1)\n",
    "    for k in 1:nitrs\n",
    "        f=F(x)\n",
    "        \n",
    "        if norm(f)<tol\n",
    "            break\n",
    "        end\n",
    "\n",
    "        #Solve BΔx=-F(x)  for Δx\n",
    "        Δx=-Binv*f\n",
    "\n",
    "        #Update x\n",
    "        x.+=Δx\n",
    "        \n",
    "        #Evaluate F at new x\n",
    "        fnew=F(x)\n",
    "                    \n",
    "        Δf=fnew.-f\n",
    "            \n",
    "        #rank-1 update of Binv\n",
    "        BinvΔf=Binv*Δf\n",
    "        \n",
    "        Binv.+= (Δx-BinvΔf)*Δx'/(Δx'BinvΔf)*Binv\n",
    "        \n",
    "        @printf \"k=%2d: %1.12f,\\t %1.12f\\n\" k norm(f) norm(Δx)\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### A simple example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "F(x)=[ x[1]^2+x[2]^2-2.0; exp(x[1]-1.0)+x[2]^3-2.0]\n",
    "\n",
    "JF(x)=[ 2x[1] 2x[2]; exp(x[1]-1.0) 3x[2]^2]\n",
    "\n",
    "x0=[1.5; 2.0]\n",
    "\n",
    "x=copy(x0)\n",
    "A=JF(x)\n",
    "\n",
    "println(\"Newton's method\")\n",
    "@time Newton!(x, F, JF)\n",
    "\n",
    "println()\n",
    "x.=x0\n",
    "println(\"Broyden's method 1\")\n",
    "@time Broyden_QuasiNewton!(x, F, A)\n",
    "\n",
    "println()\n",
    "x.=x0\n",
    "A=JF(x)\n",
    "println(\"Broyden's method 2\")\n",
    "@time Broyden_QuasiNewton2!(x, F, A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this problem, we see that Algorithm 1 and Algorithm 2 take longer time than Newton's method. When do you expect that they become faster?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Julia 0.6.2",
   "language": "julia",
   "name": "julia-0.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}