{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "23776d43",
   "metadata": {},
   "source": [
    "# Multivariate Root Finding\n",
    "\n",
    "## Table of Contents\n",
    "\n",
    "1. [Fixed-point iteration](#first-bullet)\n",
    "2. [Newton's method](#second-bullet)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1b42eb2a",
   "metadata": {},
   "source": [
    "## 1. Fixed-point iteration <a class=\"anchor\" id=\"first-bullet\"></a>\n",
    "\n",
    "We want to rewrite $F(x)=0$ as $x = G(x)$ to obtain:\n",
    "\n",
    "$$\n",
    "x_{k+1} = G(x_k)\n",
    "$$\n",
    "\n",
    "We want $G$ to be a contraction, and satisfies if $||G(x) - G(y)|| \\leq L||x-y||$, then $||x_k - \\alpha||\\leq L^k ||x_0 - \\alpha||$. So $G$ will converge to a fixed point $\\alpha$.\n",
    "\n",
    "For multivariate case, we have to use **Jacobian matrix** $J_G \\in \\mathbb{R}^{n \\times n}$, which is \n",
    "\n",
    "$$\n",
    "(J_G)_{ij} = {\\partial G_i \\over \\partial x_j}, \\quad i,j = 1, ..., n\n",
    "$$\n",
    "\n",
    "If $||J_G(\\alpha) < 1||$, then there is some neighborhood of $\\alpha$ for which the fixed-point iteration converges to $\\alpha$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4c0a1a25",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c02e6ff1",
   "metadata": {},
   "source": [
    "The following program implements fixed-point iteration to solve a system of nonlinear equations:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "x_1^2 + x_2^2 - 1 & = 0\\\\\n",
    "5x_1^2 + 21 x_2^2 - 9 &= 0 \n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "It can be rearrange to $x_1 = \\sqrt{1-x_2^2} , x_2 = \\sqrt{(9-5x_1^2)/21}$. So we define \n",
    "\n",
    "$$\n",
    "G_1(x_1, x_2) = \\sqrt{1-x_2^2} ,\\quad G_2(x_1, x_2) = \\sqrt{(9-5x_1^2)/21}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "6ec63022",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   0                  0                  0                  -1                  -9\n",
      "   1                  1     0.654653670708      0.428571428571                   5\n",
      "   2     0.755928946018     0.436435780472     -0.238095238095      -2.14285714286\n",
      "   3     0.899735410842     0.540848413886      0.102040816327       1.19047619048\n",
      "   4     0.841120082507     0.485620906056    -0.0566893424036     -0.510204081633\n",
      "   5     0.874169511937     0.510022643476     0.0242954324587      0.283446712018\n",
      "   6     0.860160975133      0.49661417054    -0.0134974624771     -0.121477162293\n",
      "   7     0.867971408296     0.502404479634    0.00578462677588     0.0674873123853\n",
      "   8     0.864632719045     0.499195933091   -0.00321368154216    -0.0289231338794\n",
      "   9     0.866489134603     0.500573542753     0.0013772920895     0.0160684077108\n"
     ]
    }
   ],
   "source": [
    "from math import sqrt\n",
    "\n",
    "# Define f(x_1, x_2)\n",
    "def f(x1,x2):\n",
    "    return (x1*x1+x2*x2-1,5*x1*x1+21*x2*x2-9)\n",
    "\n",
    "# Store results\n",
    "ks = []\n",
    "x1s = []\n",
    "x2s = []\n",
    "f1s = []\n",
    "f2s = []\n",
    "\n",
    "# Print out solution\n",
    "def print_sol(k,x1,x2):\n",
    "    (f1,f2)=f(x1,x2)\n",
    "    ks.append(k)\n",
    "    x1s.append(x1)\n",
    "    x2s.append(x2)\n",
    "    f1s.append(f1)\n",
    "    f2s.append(f2)\n",
    "    print(\"%4d %18.12g %18.12g %19.12g %19.12g\" \\\n",
    "          %(k,x1,x2,f1,f2))\n",
    "\n",
    "# Define starting position\n",
    "(x1,x2)=(0.,0.)\n",
    "print_sol(0,x1,x2)\n",
    "\n",
    "# Perform fixed-point iteration\n",
    "for k in range(1,10):\n",
    "    (x1,x2)=(sqrt(1-x2*x2),sqrt((9.-5.*x1*x1)/21.))\n",
    "    print_sol(k,x1,x2)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "1e15e086",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA9sAAAIlCAYAAADFZusBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAxOAAAMTgF/d4wjAABuaElEQVR4nO3dd5ycZb3//9e1Pb2RkEYSIPTQAol0QUGRYgEEEaRIEetXxeOxnaN4POd4fh45niMiKggqxQKIDex00WQ3ofeQwm4a6XXrXL8/7tndyWZTNpnde2b29Xw85nHPzD1z35+ZvSeZ91zlDjFGJEmSJElS/pSlXYAkSZIkSaXGsC1JkiRJUp4ZtiVJkiRJyjPDtiRJkiRJeWbYliRJkiQpzwzbkiRJkiTlmWFbkiRJkqQ8M2xL6pdCCCeHEGIIoaKXtl+R3f7JvbH9nP0sCCFcmedtTsnWPjWf290dvfE686Wnx1II4bYQwu29XZe2L4RQHUL4eQhhXQghpljHFsdDIR/rhSqE8FAI4Wtp1yFJXfXKl0xJKhQhhIeAN3e5+x/AScC4GGNrnxeVXzOADTv74OyX+tYY42XbedjrwDjgjexzrgS+FGOcsutl7nR9U4FXgL1jjAtyVvXodRa4/5d2AX1lJ4+3tFwAnAwcB6xIt5QtlNKxnlfZH7RagFNijA/lrDoHaE6lKEnaDsO2pP7gW8B/5dxujjE2A0vTKSd/Yoxv9MI22+iF9yaEUJV933usN15nWmKMa9OuQQDsA7wQY3w27UJypX2shxAqSX4gSa21v6dijKvSrkGSumM3ckn9wcYY49Kcy6rcrr8hhPIQwt9CCD9sf0II4agQQlMI4fjs7dEhhDtCCGtCCCuy10flPH5YCOGeEMLmEMLLwNt2VFQI4SshhMdCCP8UQlgeQlgdQviPEELIecx+IYQ/Zre7PITwjdzuyrldTnO6f787hDArhLAx271yUvv+gIuAS7OP6/bLdG438mw3+B8Ak9uf0941PoSwTwjhNyGEDSGExSGEG0IIA7vU9s/Z92UT8IkQwp4hhLtDCEtDCOtDCI+EEI7I2f0r2eX87L6+0vV1Zm/PDCE8kf0bvR5C+GyX1xBDCJeFEP4cQtgUQqgLIRy2nb/FqBDCL0IIq7Lv21MhhGNz1l8SQng1u79nQgjv2MZ2jgwhtIUQxnW5/7chhG9nr3fXbfja7P43hhBeCCG8pcvzrw0hLAshrA0hfDN7/N22ndfzzhDC37Pv8eIQwo0hhEE9ff0hhAkhhNYQwsFdHv+T9v2HEKZnj+ON2WP44RDC8O0dbyGEU0MItdnj+uUQwkdz1rUff+fkPObP2RrfG0KYl93P/4TQ+Vnp5jUNCiHcnH3shuxxuGf73wD4MnBSdl8PbWMbD4UQ/r8Qwvey7+WCEML7ujzmHdljoil7jFzSzWvp9jO5jX3u9Gc65zmfCCG8lj3WZ4ec4SshhANDCPeH5N+tNdnre+esb/+38PQQwnPAZmAU3dje3y1nWy9k1/82JJ//BV3ez691eU7u660OIfw4JJ/pjSH53OZ+Fl7NLh/M1nxbd9sNO/fv5jY/cyH5t+33IRlisC6E8I9QQMNqJBUPw7akfi/bknspcEEI4ewQQjXwY+DbMcbHsw+7O7s8kaTr6XAgd9ztt4BDgLdkt3XdTu7+cOBY4BTgKuBj2ecTQigHfgU0ATOz918CfLbbLXX6CvDP2ecMBP4ne/9/A/cAPyfpJj6uuyd38TfgWqA+5zl/CyFUAX8gCcdHAe8i6f76zS7P/yzwe2Aa8FNgAPAIcFr2ec8Dvw4h1GQf3x5wZ2b39d9dCwohDAHuB54Djsju48shhPd3eei/At/OPmYxcOt2Xue/AUNIhhccRvL3a87u7zjgh8D/Zdf9ErgvhDCl60ZijHNJAsF5OfUOz77en21n/58FfpOt9VHg9ux7TAjhVOA/gS+SvC9VwNnb2RZADfDvJMfX+0iOry9v5/Hdvv4YYwPwZ+ADOa9nEPAeks8IJJ+Dx4FDgROAO7L3d3u8hRAOAO4Fvkvymfk0yd/vgi41/SvJsXcsMBn4BXAxybF2MfAR4KztvKb/IRlC8q7s65oA/CS77v+RfGafyNZ1zna28yHgReBI4Dbg1hDCmOxrmQLcl70clt3mD0P2R7ocX6H7z+TO2ubzQwgfzL6ej5B8zn4M3J9zfA4m+ffrhOylmeSz2NWXSf4NOhRY13Xljv5uIYRhJJ+NB0neq98Cn+vh66wAXiY5vg8Hfg38qv39Bo7JLs8l+bttNSSjB/9ubvMzB9wALCP5N+1oks9+poevRZIgxujFixcvJXsBHiL5crkh53IRSWCOQEXOYz8OLCFpyX0BqMnefxJJt+rcx47PPn8iMJRkHOHpOetPz64/eTu1fQXYBIzIue9rQG3ONjYDI3PWXwO8kXN7AXBl9vqU7D7Pz1l/IbAi5/btwG07eM/atzM1e/tKYEGXx1zSXmfOfceRfMEtz6nt1h3sqzz7Nzkpe3tqdt9Tujwu93VeAzR0+Xt8HZidczsCn825fWz2vsHbqOM3wL9sY91PgZ93ue/vwDey17c4loCvAo/lPPYyknHwIXv7NuD2Lq/txpzb47Lbm5a9/Ysujy8HFu3o79il3vcBr21n/fZe//uy+yvL+dvn3l4PnLiN5251vJH8cPHfXe77AvDn7RzHnyMJO2Ny7vs98M1t7HcIyWfyjJz7Dsxu95Ccz9pDO3jfHgLuz7ldAWwEzso57mZ1c7z8ogefye6Oh558pl9rryfnvj+SzLPQ3Wsam93mpC7H75t38F7s6O/2YZIf5XI/lz8l59+O7Pv5tS7b6Hi929jvi8AlOe//Vv+u5m6Xnf93c3ufuWeAD+zs58uLFy9etnWxZVtSf/ADktaL9suvt/G4G0haJa8ALo0xNmbvPxQYDazJdkfdQNL6Asm4z31IvgTOytlW7nVCCM+1PzeEcFPOqldjjKu7PO+A7PUDgFfiluMRnwD2CCGM3M7rfSbn+lJgVLa1J58OBQ7PeU0bgD+RtLpOyHnc3NwnhRAqQ9JV/oUQwhpgLUlL3V492PcBQF3ccnK7J+h839p1fR8AxtC9HwBfCCE8GkL412wrXu7+/t7l8d3tr91PgeNCCO2v6XySsB638fgd1bofUNe+MiY9MZ7czrYIIRwcQvhlCGFRCGE9Sav+9t7j7b3++0jC68nZ2x8gCYftLX03AH8MIdwXQvhoCGGP7dVGcux8rMux868kn6Ncue/JMpKwtLzLfaO3sY/2z2TH3y3G+CKwhm3/3balo47sMbeCzr/Nzh4bu/uZ7Pb5IYTBwN7Az7q8n6eQfT9DMsTlOyGEV0II6+jsit31eJjL9u3o73YAMKfL53JW143sSAjhMyGEp0MypGEDyfHf038fdubfze195m4Ebg4h/CFbT0/2L0kdnCBNUn+wOsb4au4dofuhnnsA+5O0oO1D5xfFwSRfUM/s5jkNJC1mkLSMbMsZQGX2em4Xze09Z5vjUXegpZvt7+q2tmUwSXfwD3WzbknO9U1d1v0zSbfOTwAvAY0k73MlO29nX0t370O3PzLHGH8dQtiHpPvqGcAXQwiXxBh/1oP9tW/r+ezY1/eGEG4FTmX7Xbi3qDXGGLPHZ3utge0fJ935NfA0SS+O5SS9M76/nZq3+fpjjI0hhJ8Bl4QQXiIZKvHxnOd+PoRwB0mX7g8AXw0hHBNjfKWbXUFy7FxP0lKaq+uZAbr+/Vq6rI8krfzdyefx3t1+c/82Pd3Grnwmt/X89nH47ycZVpFrfXb5TZLu158E5pN893uKrT9zXT+rXe3o77Yzx2mGrV93Rx0hhItJAvzHSX5Q2kjSNb1X/33o+pmLMX43hPAHks/D2cB1IYTTY4yP9qAOSTJsS1KOG0laO34O3BBCeDDGuIzki+kkYF2XljUAQgjzSL5wziQZxwzJWL8OMcaF29jnfiGE4THGNTnPeyl7/cXs+pE5rTTHkrTw7ersuy30/N/+FrYONU8B7wTqc3oA7IxjSLrY3gOQbTEa0WVfdLO/XC+SBNmKnFa0Y7P377IY4xKSQPr9EMKNJD8K/Cy73WO6PPxYkh8btuWnJKeWWgu8HmOcvRulvUwyvh3oGJN6BPDX7h6cbVneFzgvxvhk9r7zd7ST7bx+SFrG/0zS/bY220qc+9xngWeBr2d/aHgP8P/R/fH2FHBA1x/A8qz9M3kMyfh+QggHksy1sFvHSRcvkvz4kGu3j8UeWE7SKjspxvirbTzmGODmGOPvAEIIJ+7ivnb0d3sJOCeEUJ7tfQFd/h0kOZ3g2PYbIYTRubeztf41xvij7PrBJP/2tmsjCew7+vdht//djDG+Bvwv8L8hhPtJuu8btiX1iN3IJQkIIVxIMoP4B2OM3wdm09kS+EeSEH5vCOHE7Ey1p4UQvg8QY1wH3An8TwjhTSGEY0gmnNoZbSTdFQ8OIZxD0uL7nZz9zgduCyFMC8kM2NeRTMK0qxYCR4ZkluMddffNfc6eIYSjQwh7hOTUQHeQjIX/WQhhRkhmLj87hLDVhGZdzANOD8kM1tOBH5G0brdbmt3u20IyA/zAbrZxB1ANfDckMy1fSNIS9q2dfD1bCSFcF0I4K/u3PRo4ns4fPf6PJER8LISwfwjhqyQTQN24nU3+jCRofJLtT4y2M74HnB9C+GC2e/f1JKFxW62Iq7OXq7Kv5wK674HQYQevnxjjP0jGnX+BzonRCCEMCCH8X/ZzMTmEcAZJOGp/bnfH2zeAs0IIX8se94eEZOb4j/TkTdmeGON6khbYb2Vrm04yNvpPMcbn87UfksnCDg8hfDV7bHyMZHK8b+VxH9uUHZrwH8C/hRAuDyHsm/2cfi50zq49Dzgv+16fQPL+74od/d3uJJm/4n9DCAeEEK4G3t5lG4+QfJbeEkI4FLiZZJ6HdvNIhmCcGEI4hOTfh47vqtnX+zrwlhDCmGwY72q3/90MyUz3p2aP2xNJJr97aUfPk6SuDNuS+r0QwliScaefijEuyt59JclpgT6QHZt6OsmXrXtJumt+m2T8Z7tPZtc/TBIGdzZsPwXUknwJ/SHJl/fbALL7fRfJDN6zSb54/pikxXBX/QBYRTIL+M6ez/cRkpbaP2efc3w2zJxMEoz/lH0dX2PLLuTd+RrJF+HHSGaq/j6wsn1ljLEJ+CeSrqTL6Gbm9ey+zyAZQ/oUSQi4LsZ4506+nu60ksye/TzwO5Ku7V/K7u9vwAdJ/sbPkrTavjvGuGBbG8u2/s0hmR16t8J2jPHPJCH369m6Wkn+Fk3beHwbSffxt5Ecqx8ieT+3Z5uvP8ePSQJ+7kzWbSTjXO8iaYG/AfhqTivrVsdbjLGOZHb2N5OMRX8MuJyk1TyfriVpifwNyTHcQM6s6vmQ7bHybpJj4lmSY+SK7DHTJ2KM3yb5nHyWZGLH35D0smnIPuRakq7VdSSftx0dC9vaz3b/btneOe/JPuYpOns35LqZ5N/Qe0mOs7tIWufb3QT8haQ3wp9I/n5PddnGZ0mO7yUkx1vXOvPx72YlyXv1YrbGO7vblyTtSPvMqJKkPhaS8xCfGmM8Ie1aVDxCMrj0JeAHMcZdbaXclf3+HzAxxri902RJHUJy/uwvxRinpF2LJKXBMduSJBW4EMJngAdIWpY/QtJV+xd9tO/BJOc8vhR4b1/sU5KkUmA3ckmSCt9JJF2hZwFHA6dtrxt7nt1A0m39jhjjH/ton5IkFT27kUuSJEmSlGe2bEuSJEmSlGeGbUmSJEmS8sywLUmSJElSnpXMbOTV1dVx9OjRaZchSZIkSeonGhoammOM1d2tK5mwPXr0aOrr69MuQ5IkSZLUT4QQ3tjWOruRS5IkSZKUZ4ZtSZIkSZLyrGS6kUuSJEmS8ieTyRBjTLuMVIUQKCvbtTZqw7YkSZIkqUNzczOLFi2ipaUl7VIKQmVlJZMmTaKqqqpHzzNsS5IkSZI6LFq0iCFDhjBq1ChCCGmXk6oYIytXrmTRokVMnTq1R881bEuSJEmSgKTreEtLC6NGjaKiwrgIMGrUKFatWkUmk+lRl3InSJMkSZIkAXSM0e7vLdq52t+Lno5fN2xLkiRJkpRnhm1JkiRJkvLMsC1JkiRJKnh33HEHhx56KNOnT+fJJ5/cYt0555zDE0880eNtxhg58cQTmT9/fp6q7GTYliRJkiQVtEwmw7XXXssDDzzAnDlzOOKIIzrWzZo1izVr1nDsscf2eLshBD71qU9x3XXX5bHahGFbkiRJklSwVq1axYEHHsj69es566yz+O53v7vF+u9973tcdNFFHbevuOIKPv/5zwOwcOFCDjjgAB577LFtbv/ss8/m/vvvZ/369Xmt27ncJUmSJEnduvJHs1m4clOvbX/yqIHcfOmM7T5m5MiRfPWrX+W+++7jpz/96VbrH3roIT7zmc903P7a177G4Ycfzvnnn88HPvABbrjhBk444YRtbr+yspJp06bx+OOPc/rpp+/6i+nClm1JkiRJUkGrra3lqKOO6nZdfX09Y8eO7bg9btw4rrnmGo4//nj+/d//ndNOOw2A1atXc8UVV7DXXntttY2xY8dSX1+f15pt2ZYkSZIkdWtHrc59pa6ujn/5l3/pdt3AgQPZvHkzI0aMAGDFihX8+te/ZujQoVsE6xEjRnDLLbdw6qmnbrWNxsZGBgwYkNeaDduSJEmSpIIVY2Tu3LlMnz692/WHHXYYL774IuPHj2ft2rWcccYZfP7zn2fz5s1ce+21PPjggzvcxwsvvMDhhx+e17rtRi5JkiRJKlivvPIKo0aNYvjw4d2uP++883jggQfYtGkTZ511FldffTUXXHABl1xyCStXruS+++7b7vYXLFgAwLRp0/Jat2FbkiRJklSw9t9/f+bNm7fN9ZdddhkPPPAAMUYeffRRrrzySgDKysp4+umnefe7393x2GuuuYYXX3yRa665pmObN9100xYTrOWL3chV+lbOg0GjoWZo2pVIkiRJyrMhQ4bwrW99i/nz5++wdfqmm27a6r7x48dz+eWX570uw7ZKW9MGuOlE2PcUeN8daVcjSZIkqRd0N+nZzvrEJz6Rx0o62Y1cpW3xXGjZCC/+DlYvTLsaSZIkSf2EYVulraE2eyVC7S2pliJJkiSp/zBsq7Q11EEog+GTYc5PoKUx7YokSZIk9QOGbZW2+joYczC86RrYvAqeuzftiiRJkiT1A4Ztla51i2H9YphwFBzxfqgcCLO+n3ZVkiRJkvoBw7ZKV312vPaEo2DAcDjs/GTCtPq6VMuSJEmSVPoM2ypdDdlQPfHoZDnjqmQ5+wfp1CNJkiSp3zBsq3Q11EHVYBh9YHJ77DSYdCw8ey9sXJFubZIkSZJKmmFbpSnTlnQZH38klJV33j/zKmhrgjk/Tq82SZIkST12xx13cOihhzJ9+nSefPLJLdadc845PPHEEz3eZoyRE088kfnz5+epyk6GbZWmN16E5g0wYfqW9x94NgzeE2pvTQK5JEmSpIKXyWS49tpreeCBB5gzZw5HHHFEx7pZs2axZs0ajj322B5vN4TApz71Ka677ro8VpswbKs0tY/XnnD0lvdXVMFRl8HaRfDyH/q8LEmSJEk9s2rVKg488EDWr1/PWWedxXe/+90t1n/ve9/joosu6rh9xRVX8PnPfx6AhQsXcsABB/DYY49tc/tnn302999/P+vXr89r3YZtlab2mcgnHr31uqMuh7IKTwMmSZIkFYGRI0fy1a9+lbPPPpsnn3ySD3/4w1usf+ihhzjuuOM6bn/ta1/jlltuYe7cuZx55pnccMMNnHDCCdvcfmVlJdOmTePxxx/Pa90Ved2aVCga6mDIeBg6fut1Q8fBgWfB8/fBildhj6l9Xp4kSZJUFO58H6zO/3jmDiP2hvf/dIcPq62t5aijjup2XX19PWPHju24PW7cOK655hqOP/547rrrLk477TQAHnzwQW6//XYaGxsZMWIEN9xwQ8dzxo4dS319/W6+mC3Zsq3S07QBlj+/9XjtXDPbTwN2c9/UJEmSJGmX1dXVbTNsDxw4kM2bN3fcXrFiBb/+9a8ZOnQoe+21V8f9p5xyCrfccgt33HEHCxcu3KLbeGNjIwMGDMhrzbZsq/QseQpipvsu5O0mHw+jD4In74C3fAmqB/ddfZIkSVKx2IlW594WY2Tu3LlMn959Y9phhx3Giy++yPjx41m7di1nnHEGn//859m8eTPXXnstDz744BaP/+1vf8vBBx/MkCFDOu574YUXOPzww/Naty3bKj0N2fHaXSdHyxVC0rrdtA6e+Xnf1CVJkiSpx1555RVGjRrF8OHDu11/3nnn8cADD7Bp0ybOOussrr76ai644AIuueQSVq5cyX333dfx2Jtvvpna2lr+67/+q+O+BQsWADBt2rS81m3YVumpr4VQlpxje3sOuwCqh8KsmyHGvqlNkiRJUo/sv//+zJs3b5vrL7vsMh544AFijDz66KNceeWVAJSVlfH000/z7ne/G4B77rmH6667jqVLl3LNNdfwxhtvAHDTTTfxmc98Ju91241cpaehDkYfuOOu4dWD4fALYdb3YNETMPm47T9ekiRJUsEZMmQI3/rWt5g/f/52W6fPPfdczj333K3uHz9+PJdffnne67JlW6Vl3RJY1wATup88YSszkl+9PA2YJEmSVLxOPfXUXe4G/olPfIKysvxHY8O2SktDXbLc3uRouUbvD/ucDC/8BtYv7bWyJEmSJPUvhm2Vlp2ZHK2rGVdBphXqbuuVkiRJkiT1P4ZtlZb6WqgcmIzZ3ln7nw7D9oLaW6GtpfdqkyRJktRvGLZVOjJtsPjJZBby8h7M/VdeAUdfDhuWJt3JJUmSJGk3GbZVOla8DM3rd35ytFzTL4XyKph9c/7rkiRJkopECAGA6KlxO7S/F+3vzc7y1F8qHfXZ8do7OzlarkF7wCHvgad/Bsuegz0PyW9tkiRJUhEoKyujsrKSlStXMmrUqB4HzFITY2TlypVUVlb2eMZyw7ZKR8fkaLvQsg0w8+okbM/6AZz9rbyVJUmSJBWTSZMmsWjRIlatWpV2KQWhsrKSSZMm9fh5hm2VjoY6GDwWhk7YtedPOArGHQFP/xxOuw5qhuW1PEmSJKkYVFVVMXXqVDKZTL/vTh5C2OVzcBds2A4hfBn4CnBojPHZlMtRoWveBMuehwPeAbva1SUEmHkV/Oqj8ORdcMw1+a1RkiRJKiK7GjKVKMh3L4QwHTgGWJR2LSoSS56E2LbrXcjbTTsXBoyA2T+ATCYvpUmSJEnqfwoubIcQqoHvAB8B+nefBe28+t0cr92ucgAc+QFY+SrMf2i3y5IkSZLUPxVc2Aa+CtweY5yfdiEqIg11QEjOsb27ZlyRbGuWpwGTJEmStGsKKmyHEI4FZgA37sRjPx1CqG+/bNiwofcLVOFqqIPRB0LN0N3f1ogpsP/b4eUHYI0jGSRJkiT1XEGFbeDNwIHA/BDCAmAi8IcQwju6PjDGeH2McWL7ZfDgwX1cqgrG+mWw9nWYuJtdyHPNuApiBmp/mL9tSpIkSeo3Cipsxxi/HmMcH2OcEmOcAtQDb48xPpByaSpku3t+7e7s+xYYuQ/M+TG0NOZvu5IkSZL6hYIK29IuaahLlhOOzt82y8pgxpWwaSU898v8bVeSJElSv1DQYTvbwu05trV99bVQORDGHJzf7R7xfqgYkJwGTJIkSZJ6oKDDtrRDmQwsngvjjoDyivxue8AIOOy9Sct5w5z8bluSJElSSTNsq7iteBma1sGE6b2z/RlXJcvZngZMkiRJ0s4zbKu4tY/XnpjH8dq5xh0Gex0Dz9wNG1f2zj4kSZIklRzDtopbx0zkvRS2AWZeBW1NMPcnvbcPSZIkSSXFsK3iVl8Lg/eEYRN7bx8HvRMGjYHaWyDT1nv7kSRJklQyDNsqXs2bYNlzyfm1Q+i9/VRUwVGXwZpF8Mofe28/kiRJkkqGYVvFa+nTENuSsN3bjr4cQjnM8jRgkiRJknbMsK3iVZ8dr91bk6PlGjoeDjwT5v0FVs7r/f1JkiRJKmqGbRWvhlogwPgj+2Z/M69Olp4GTJIkSdIOGLZVvOrrYI/9oWZY3+xvygkw+iCYewc0b+ybfUqSJEkqSoZtFacNy2Htor7pQt4uBJhxBTSthWd+0Xf7lSRJklR0DNsqTg11ybIvJkfLdfj7oGoIzLoZYuzbfUuSJEkqGoZtFae+nBwtV/UQOOJCWPYMLPp73+5bkiRJUtEwbKs4NdRCRQ2MObjv9z3jymQ529OASZIkSeqeYVvFJ5OBhrkw7ggor+z7/Y8+APY+CZ7/Faxf2vf7lyRJklTwDNsqPitfTSYp6+su5LlmXg2ZVqj7UXo1SJIkSSpYhm0Vn4bseO2+nhwt1/7vgKEToe5WaGtJrw5JkiRJBcmwreJTXwBhu7wCjr4M1i+BF3+XXh2SJEmSCpJhW8WnoQ4GjYbhk9KtY/plUF4Fs5woTZIkSdKWDNsqLi2bYdmzMOFoCCHdWgaPhoPfDQsfg2XPp1uLJEmSpIJi2FZxWfJ0MjFZml3Ic828KlnOvjndOiRJkiQVFMO2ikv75GgTCyRsT5wBYw+Dp34KjWvTrkaSJElSgTBsq7g01CXL8dPTraNdCMlpwFo2JoFbkiRJkjBsq9jU18Ie+8OA4WlX0mnauVAzPOlKHmPa1UiSJEkqAIZtFY+NK2DNwsIZr92uaiAceTGseBnmP5x2NZIkSZIKgGFbxaMQzq+9LTOuAIKnAZMkSZIEGLZVTNrHa088Ot06ujNyH9jvNHjpfljzetrVSJIkSUqZYVvFo6EWyqthz2lpV9K9GVdBzEDdrWlXIkmSJCllhm0Vh0wmadkedziUV6ZdTfemngojpkDdj6C1Ke1qJEmSJKXIsK3isGpech7rQuxC3q6sDGZcCZtWwHP3pV2NJEmSpBQZtlUc2sdrF+LkaLmOuAgqamC2E6VJkiRJ/ZlhW8WhfSbyQm7ZBhg4Eg49D+pnw+K5aVcjSZIkKSWGbRWHhloYOAqGT067kh2bcVWynHVzunVIkiRJSo1hW4WvpRGWPgsTjoYQ0q5mx8YfARNnwrN3w6ZVaVcjSZIkKQWGbRW+pc9ApqXwu5DnmnkVtDbC3NvTrkSSJElSCgzbKnwN2fHahT45Wq6D3wWDRsPsmyHTlnY1kiRJkvqYYVuFr31ytAnT062jJyqqYfqlsGYhvPrntKuRJEmS1McM2yp8DbUwaioMGJF2JT1z9OUQymCWpwGTJEmS+hvDtgrbxpWwekEyOVqxGTYRDjwTXv0TrJyXdjWSJEmS+pBhW4WtoS5ZFtPkaLnaTwNW+8N065AkSZLUpwzbKmwNRTheO9feJ8EeB8Dcn0DzprSrkSRJktRHDNsqbPW1UF4Nex6adiW7JoTkNGCNa5PzbkuSJEnqFwzbKlwxJt3Ixx0GFVVpV7PrDrsAqgbDrO8nr0mSJElSyTNsq3Cteg0a1xTn5Gi5aobC4e+Dpc/A67PSrkaSJElSHzBsq3B1nF/7qHTryIcZVybL2Z4GTJIkSeoPDNsqXO2To00sgbA95iCYciI8dx9sWJ52NZIkSZJ6mWFbhauhDgaOghF7p11Jfsy8CjItUPejtCuRJEmS1MsM2ypMrU3JGOcJRyUzepeCA86EIeOTc263taZdjSRJkqReZNhWYVr6DLQ1l8Z47XblFXD0B2H9Ynjpd2lXI0mSJKkXGbZVmDomRyvymci7OupSKKuEWU6UJkmSJJUyw7YKU0NdspwwPd068m3wGDj4XbDgUVj+YtrVSJIkSeolhm0VpoZaGLkvDByZdiX5N/PqZDn75nTrkCRJktRrDNsqPJtWwarXSmu8dq69ZsLYQ+Gpu6BxXdrVSJIkSeoFhm0VnvYu5BNLbLx2uxBgxlXQvAGe/lna1UiSJEnqBQUXtkMINSGE+0IIL4cQngwh/D6EMCXtutSHOsZrl2jYBjj0vVAzLJkoLca0q5EkSZKUZwUXtrO+DxwQYzwC+G32tvqL+loor4Kx09KupPdUDYQjPwArXoL5j6RdjSRJkqQ8K7iwHWNsjDHeH2NHc9/fgX3SrEl9KMakZXvsoVBRnXY1vevoDybL2Z4GTJIkSSo1BRe2u/EJ4DdpF6E+sno+bF5V2l3I243aF6aeCi/eD2vr065GkiRJUh4VdNgOIXwB2A/4YjfrPh1CqG+/bNiwoe8LVP7Vl/jkaF3NvBpiG9TemnYlkiRJkvKoYMN2COEzwDnAO2KMm7qujzFeH2Oc2H4ZPHhw3xep/GuoTZaletqvrqaeCsMnw5wfQWtT2tVIkiRJypOCDNshhE8DFwKnxRjXpFyO+lJ9LQwYASP7yTD9snKYcQVsfAOe/3Xa1UiSJEnKk4IL2yGEicA3geHAg9nTf/0j3arUJ1qbYenTSat2CGlX03eO/ABU1DhRmiRJklRCKtIuoKsYYz3Qj5KWOix7Btqa+8fkaLkGjoRp58GTt8OSp2Dc4WlXJEmSJGk3FVzLtvqx/jY5Wq6ZVybLWbZuS5IkSaXAsK3C0T452vjp6daRhvFHJi36z9wNm1alXY0kSZKk3WTYVuFoqIMRe8OgUWlXko6ZV0PrZnjyjrQrkSRJkrSbDNsqDJtXw8pX+2cX8naHvBsG7gGzb4FMJu1qJEmSJO0Gw7YKQ0N2vHZ/mxwtV0U1TL8EVs+HeX9JuxpJkiRJu8GwrcLQPjnahKPSrSNtR38QQhnM+n7alUiSJEnaDYZtFYaGOiirhLGHpl1JuobvBQecAa/8CVbNT7saSZIkSbvIsK30xZjMRD72UKisSbua9M24EohQe0valUiSJEnaRYZtpW/1Ati0sn9PjpZrn5Nh1H4w5yfQvCntaiRJkiTtAsO20tfgeO0thAAzr4LGNfDsPWlXI0mSJGkXGLaVPmci39rh74PKQTD7B0k3e0mSJElFxbCt9NXXQs1wGLVv2pUUjpphcPgFsOSp5P2RJEmSVFQM20pXa3MSKCcclXSfVqcZVyVLTwMmSZIkFR3DttK17Floa3K8dnf2PBgmnwDP3wcb3ki7GkmSJEk9YNhWutrHazsTefdmXgltzTDnR2lXIkmSJKkHDNtKlzORb9+BZ8GQcVD7Q2hrTbsaSZIkSTvJsK101dfCiCkwaI+0KylM5ZVw1OWwrgFefiDtaiRJkiTtJMO20rN5Nax8xVbtHTnqUiirgFk/SLsSSZIkSTvJsK30LJ6bLD2/9vYNGQsHvwvmPwxvvJR2NZIkSZJ2gmFb6al3crSd1n4asNk3p1uHCtvGlXD3FfCP70FbS9rVSJIk9WuGbaWnoTbpHj32sLQrKXyTjoE9p8GTd0HT+rSrUSFq3gh3ng/P3g0PfBa+exy8/AeIMe3KJEmS+iXDttIRYzI52p7ToLIm7WoKXwgw40poXg9P/yztalRo2lrgF5clP2Ad93E4+Quwtj4J37efA8ueT7tCSZKkfsewrXSsWQSbVtiFvCcOOx+qhyUTpdlaqXYxwm8+Ca/8EQ5/P5z2b3DyP8PH6+DwC2HeX+Gm4+G3n4KNK9KuVpIkqd8wbCsdDbXJ0snRdl7VIDjyInjjRVjwWNrVqFD89Wvw5O0w9TR45/8lvSAAho6H99wEV/0VJs5MztX+f0fC4/8HrU3p1ixJktQPGLaVDidH2zUzrkyWsz0NmEh6OTz63zB+Opz/o+S87F1NOAo++Hs471aoGQ5/+hf4zpvghd/YQ0KSJKkXGbaVjobapEv0yH3TrqS4jNoX9n0rvPBbWNuQdjVK0/O/gvv/KfkMXfSLpOfDtoQA086Bj82Gt34ZNr4BP7sYfnQ2LHmq72qWJEnqRwzb6nttLckX/AnTocxDsMdmXgWxDepuS7sSpWXB43DPVTBoNHzgXhi0x849r7IGTvw0fHwOTL8kGY7wvTfDrz4K65f1bs2SJEn9jElHfW/Zc9DaaBfyXbXf22DYpCRstzanXY362rLn4K4LobwKLr4bRkzp+TaG7Anv/DZ86BGYcgLMvR2+PR0e+W9o2Zz3kiVJkvojw7b6npOj7Z6ycphxBWxcDi/8Ou1q1JfWvA63nwctm+CCn8C4w3dve+MOg0t/A++7EwaPgb/+G9wwE569x/HckiRJu8mwrb7XPjnahKPSraOYHfkBKK9OJshS/7BpFdx+LqxfnMwyvu8p+dluCHDgmfCRf8Db/h0a18LdH4Qfvr3zsypJkqQeM2yr7zXUwfBJMHh02pUUr0GjYNq58PrfYekzaVej3tayGe56H6x4Cd7+H3DoefnfR0UVHPcx+MTcZNb7+lq4+S1w79VOxidJkrQLDNvqW41rYcXLdiHPh5lXJUtbt0tbW2vS0vz6P+C4j8OxH+3d/Q0aBWd+Ez78eDLz/dM/g28fBQ/+JzRv7N19S5IklRDDtvpWwxwg2oU8HyZMT97Hp38Om1enXY16Q4xw/7Xw0v1w6Plw6lf7bt9jDkpmOr/o7qQnysNfh28fDU/9FDKZvqtDkiSpSBm21bfaJ0dzJvL8mHEVtG6GJ+9MuxL1hof/K5l1fp9T4F3fSedUefudlrRyv+MbybH2yw/BzW+FRX/v+1okSZKKiGFbfathDpRV7P4sykoc8h4YOApm32xrY6mpvRUe+s/ks3LBT5Ix1Wkpr4Q3XZ2M5z7mI7D06WQCtV9cBqsXpleXJElSATNsq+/EmEy6tOchUDkg7WpKQ2UNTL8EVr0G8/6adjXKlxd/B7/7dHIO7YvuhuohaVeUGDACTv/PZOby/d8Bz/0SbpgBf74OmtanXZ0kSVJBMWyr76x9PTk3tOO18+voD0Iog9lOlFYSFv09mRBtwEi4+N7k/NeFZo+p8P6fwiW/glFT4bHr4f+mw5wfQ6Yt7eokSZIKgmFbfac+O17bmcjza/gk2P90ePkPsHpB2tVodyx/Ee68AEI5XPQLGLVv2hVt3z4nwzWPwlnfgpiBX38cvv9mmP9I2pVJkiSlzrCtvtNQlyydHC3/Zl4FRJh9S9qVaFetbYDbz4XmDXDBj5PZ5otBWTkcfTl8Yg4c///gjZfgR2fDTy+ClfPSrk6SJCk1hm31nYY6qB4Ko/ZLu5LSs/fJSXfeuT+Bls1pV6Oe2rwG7jgP1tUns45PPTXtinquZhic9lX46Cw4+F3w4m/hO2+CP3wxeX2SJEn9jGFbfaOtBRY/CeOPTOf0RaWurAxmXJmcb/vZe9OuRj3R0gg/fT8sfz4Jq4e/L+2Kds/IveH8H8Nl98OeB8MTN8C3pycz5re1pl2dJElSnzH1qG8sfz45R69dyHvP4RdC5SCY9f1k5ncVvkwb3HsVLHw8OaXWcZ9Iu6L8mXI8XPUQvOtGKKuE310LN50Ar/4l7cokSZL6hGFbfaN9vLaTo/WeAcPhsPNhyZOd77cKV4zwwD/DC7+GQ86Bt/07hJB2VflVVgZHXgQfr4OTPgur58Pt58Ad74U3Xk67OkmSpF5l2FbfqG8P23132q9MJvLAM0t43/ef4AO3/IPr//QyD720nLWbWvqshj4386pkOcvTgBW8R7+ZnK5tyonwnptKe3hF9WB4yxfhY7Vw6HvhlT/CjcfA/Z+FTavSrk6SJKlXhFgi3U0nTpwY6+vr0y5D2/KdN0HzRvjUs72+q5a2DL9+cjE3PvQq897YSHVFEmKaWjMdj5k6ZjDTJw1n+qQRTJ88gqmjB1NWViKtij98BzTUwqdfgEF7pF2NujP3dvjVR2HPQ+Hy3yWTi/Unr8+G338uOU5rhsPJn0vmHCivTLsySZKkHgkhNMQYJ3a7zrCtXte4Dr4+KZmh+Pwf9d5uWtr4Re3r3PTwazSs2cyQ6gouOW4ylx+/N0NrKnlhyTrmLFrNnEVrmLNwNQ1rOmftHlJTwRF7dYbvI/YazrABRfrF/9l74e7L4a3/Cidem3Y16urlP8BdF8LQCXDln2DI2LQrSkeM8Mzd8OevJLOwj5qadKXf/+2l151ekiSVLMO20vXaw/Djd8LbvgbHfTzvm1/f2MId/1jEzY/OZ8WGJkYOquKKE/bm4mMmbzcwL1/XuEX4frphLc3Z1u8QYOrowdnwnYTwfYul9butBf5nWtJK+P+eSs6DrMJQXwu3nQWVA+CKP8IengaP5k3wxHfgseuhZRPsczK8/T9gz0PSrkySJGmHDNtK16PfhL98FS7/PUw+Nm+bXbWxmdsen89tf1vAusZWxg2r4eqT9uF9MyYxoKrnAbO5NcPzS9YxZ+Fq5ixazdxFa7Zo/R5aU8ERk0Z0dD8/YtJwhtYUaOv3g/8JD38d3ncnHHhm2tUIYMUrcMvboLURLv2NM/N3tW4J/PXf4Mk7IJTB9EvhlC/C4NFpV6ZCtXk1LHkKBo6CEVOgekjaFUmS+iHDttJ11/vh5d/D5+uhauBub27p2kZ+8Ohr3PmPRWxuaWPKqIF8+OR9ec+RE6mqyO8kU8vWNXaE7zmL1vBMl9bv/cZkW7+zLeD77FEgrd/rlsC3psGUE+CSX6VdjdYtSYL2uga48K6kq7S6t3gu/P4LsOhvUD0UTvoMvOkaqKhOuzKlraUR6mfBaw8ll8VzIXbOxcHAPZLzvI+YAiOyy5F7J9cH71nakxBKklJj2FZ6YoRvHgCDx8A1j+3Wphau3MhND7/GPXX1NLdlOHDsED56ylTOOHQc5X0UcJtbMzy3eG3S9XzRauYuXM3itY0d64cNqMwZ+z2cI/YazpC0Wr9/cRk898tkBmi7K6encS3ceiYsewbe9R048uK0Kyp8MSanRPvjv8CahUloOu2rcNA7Hc/dn2QysPTpznC96O/Qmu1tVDMc9j4RJh0LTeth1fzk1HKrF8CGZVtvq6ImG8KnJOF7i1A+2R9zJEm7zLCt9Kyth/85BI66HM7+1i5t4qWl67nxoVf5zVOLyUSYPmk4H3vLVE45YAyhAL54L12bHfudbQF/tmEdzW2drd/7jxnC9MnDOTLbAr7v6EF9U/eCx+G2M5JWwXf8V+/vT1trbYLbz4UFj8Jb/iVppdXOa2mEf9wEj/w3NK+Hyccn47nHH5F2Zeotq+Z3huv5j8Dm7Knhyqth0jHJmP59ToZxh297PormjbB6YRK+V2UDePv1NYsg0/X0jwGGjs9pDZ+SvZ4N5QNG+COPJGmbDNtKz3P3wS8u3aUWvbmLVnPjQ/P40/NJK8WJ++3BR06eyjH7jCyIkL0tTa1tPLc4Gfs9N9sCvqRL6/eR7acdmzSCw/ca1jut3zHCd4+Hta8npwGrHpz/fWjbMhm454NJ74IZV8EZ3/AL+67asBwe/HeY8+PkuD7iInjrv/TfmdxLycYVMP/hzoC9ZlF2RUh+VNn7zUm4nnRMMrHg7sq0JcM5Vi/YsjW8/Xrj2q2fUz00p0v6lC1bxodOhPKK3a9LklS0DNtKzx+/BH/7Nnzk7zDmoB0+PMbIE/NW8p2HXuXxV1cC8PZD9uQjJ0/l8L2G93KxvWfJ2s3MWbgmO/Z7Nc91af0+YM8h2Zbv4UyfPIJ99shT63ftD+G3n4Izr4cZV+z+9rRzYkzOI/2Pm5Kuz++9zVnh82Hps/CHLyThrHIQnPgpOPZj+Qlh6hvNG2HhEzD/oSRcL32mc93IfTpbrqecCANH9n19m1dv3Rq+ekFyWVsPdPnOVFYBw/bqHBveNZT7I6cklTzDttLzw3ckX6Y+t3C7YSOTifzlxeV858FXefL1NZSXBd55+Hg+fPK+7L9n6c0w29TaxrMN65ibDd9zFq5h6brO1u/hAys5Mue834fvNZzB1bvQetK0Aa4/CIZNhA//zZbVvvLYt+DPX066PV98L1TW9MpuGlvaWLK2kSVrNrM4u1y2vpGKsjIGVZczuLqSwdXlDK6pYHB1JYOqyxmSXQ6uqWBIdSU1lWUF3VNkKzEmEy7+4Yuwal4SdE79Ckw71+O7ELW1JhOZtbdcv/6Pzm7cA/foDNf7vBmGT0qvzp3R2gRrXt+6Nbz9euvmrZ8zaPQ2xolPSXpmeMxKUtErqrAdQtgP+BGwB7AGuCzG+PyOnmfYLkBtrfD1vZJTHF36m24f0tqW4XfPLOHGB+fx0rL1VJWX8d6jJ/Khk/Zl0qjdn7m8mCxes7kjeM9ZtJrnFq+lpS35fJYF2H/PIUyfnJ35fNJw9t7Z1u8H/jlpYb3sfphyfC+/CvHUT+GXH4IxB8PlD8CA4bu0mZa2DMvWNbJkbSOL12zuWC5e08iStcntVRubd7vcsgCDqyuSS00Fg7LXh9RUMKiqIhvUk8ug7P3t13OfN7i6guqKPgzurc0w++bkFHeNa2HiTDj9Pz2lWtpihBUvw2vZruELHoWmdcm6yoHJD1DtAXvMwaUzQ3iMycRsHSF8wZYt4xuXb/2cigGdk7Z1bRkfPslJ2/qj1qbk37PGtbB5Tfb6ms77mjdAeVVybFTU5Fyytyu73O76mMoB9rKSekGxhe2/Aj+OMd4WQjgPuDbGuMOTMxu2C9DSZ+CmE+CET8OpX95iVVNrG/fOaeCmh+excOUmBlaVc/Exk7nyhL0ZM7R3WgGLTWNLWzLzeU7382XrmjrWjxhY2dn1fFLS+j2ou9bvFa/ADUfDwe+G83/Udy+gP3r1z3DnBTB4LFz5p2TSpW5kMpEVG5poyAnRS9YmIbo9TC9f30R3/zyXlwXGDq1h3LAaxg0fwPhhNYwfPoBx2eWeQ2uIMbK+qZUNja1sbGplfVOy3NB+adzy+sbm7u5r6/HLrygLSWCvyob13DBe1Rnmh7SH9Zqc61sE/nKqK3byC+GmVfDQf8LsWyC2waHnJ//eDOv2/zz1hnVLthx3vX5Jcn8oT3782OfkZOz1xBlQUbVLu2jLRJpbMzS3Zmhqa+u43pqJBKCsLFBRFigLgfKynEsIlGWv564vC/Rtj46mDZ3d0btO3LZmEWRauzwhwNAJ2RA+eetTmTlpW2Fqa+kMxrkheXPO9a0uOetaG3ewgzwoq8gJ5ANygnk2jOfe3q313fwY4PwGKlFFE7ZDCGOAl4E9YoytIfmfcAlwTIxxwfaea9guQLW3wm8/CRfcAQedBcCm5lbu/McifvDoayxb18SwAZVcfvwULj12CiMG7dqXsP4ixsjitVue9/v5Lq3fB4wd2hG+p08ewZRRA5MvlD9+d9LC9MlnYei4dF9IqWqYA7edRSyvZO37f0d9+V5bBeglaxpZvHYzy9Y1dvzdcoUAewyuZnw2RI8bNoDxw5PluOE1jB82gNFDqvvkVHeZTGRjcysbm9rY0NTC+sbO6xua2tjQ2MLG5jbWN7ayoamFjU1t2cckgb095G9obGVzS8+De1V5WUfw7ugOX13B4JrO67khfWzzQg597huMWvIwmYoa1k//CPG4TzBoyDAqy0uk9bSPtGUiLW0ZmrKBtrkt0xFum1sztG5eTU3D3xmy+DGGL32Cwete7XjumkH7UD/iTSwcPpOFg49gAwO32kZTl+21r2vJLpu67rMtQ1sm/99VygJUlJVRVkZHKK8oaw/j3Yf2LQJ7WaA8dxvZdTu1jezt8rJAOW2MbHuDkU2LGdncwIimBoY3LWZ4Yz3DNtdT3bZhq9qbK4awYdBebBo4kU2DJrFp8F40DZlE4+BJNA8aT3lFOeVlZdl9Jq+vojz3x4bkdsePEV1q7b727GvN+aEixkiMyUj2TIxk2m/H5Hb7/cl9kUyXZefzOreVyVkmf/Zkmemyfut9tK/vZp852yBuuc/cfcS2FspbNlDevJaK5vVUtqyjvGkdFS3rqGxJbld2XF9PVcs6KlvXU9V+PbPzYbmlrIbG8iE0VQymsWIIjeVDaCwfTGP5YDaXD2FzWXZZPojG8iFsKhtEUxhIWWylkmYqM01UxmYqYzMVmWaqYjMVsYmKTDMVsZnKTOft5DFNVGQf23E7e195Jtleec5zymPXH4B2XyZU0FpWRVtZFW1l1bSVVyfLLrcz5dnbufeVVdFWXkOmvJq28ioy2XWZ8uqO623Z67G8JruNGjIVVcTyaihLJqJt/40qhED7/6QhQCB0rut4XNjiOVtGps4b7ffnru68b0eP2/rftty72p8fM5EQ2wixDTKtlMU2iBlCbCVk2giZ1uR6bCPk3F8Wt1xH+30dl+zzc26XdTyuNef+Nsoynetz7+/YT/bxZbl15Gyn/XnJ7cwW2+q6rn17zZf+gWETD9zWIVUwiilsHwX8JMZ4cM59s4DPxBgf6fLYTwOfbr89bNiwCWvWrOmrUrUzfvVRmHs7XPsSa8tH8aMnFnDr4/NZvamF0UOquerEvXn/mybv2lhkAUnr97MNa7fofr58fWfr98hBVRy513DOHfQUZzx3Lc0nfJaqU7+YYsXFb31jy5at0Ws207T8VT7y2oepyWzmsrYv8kRL9+c1HzGwstsAndsqXVVResGwtS3Dxua2jiCeG8o35LTAb8hthc9pac99XlNrZrv7enPZU3yx4nb2L2tgaRzBN1ou4HdlJzG4pqojnLe3rm/V0t6lW3x7t/ncdfn+oSOTiTRvJ9g2t7XR3Bq7hM/Olt2tgmnu7bZu7usafHci2FbSypHhFY4vf5YTyp7l8DCPipD8HZbGETyWOZTH2qbxt8whLGfETr3u8rJAVXkZVRXZS3mXZdfrFWVUd7OuorwMYqQtRtoy0JbJ0JZJwlRbJnt/W7LMZG+3ZrLXM7Hjca051zMZaM1kaIt0PK59W+3baMu5f4ttdDyOjsftnsgwNjI5LGNyWMZeYTmTw3ImlyXXx7GKsrDlPlpiOfVxDxbFPVkUx7Awu2y/vYnd7z1WFpLAUEBfIQEIZBjCZoaGTQxlI0PDJoaxkaFhI0PZtMX9ye2NW9weEroZd78NjbGSdQxiXRzIOgayLg5iHQNZm12u67LMvX89A2mhsL/7lJGhmmaqaUkuoYWajtvNXW53Xd9CdWimpv2xXdeH5DE1NFPVcbu5877Q8x9od6Q1ltFIFU1U0kQljbHzehNVNMXsspv1zVRQToYKMpTTRnl22X67IrTfbqMs53G5y7KQoYK2re7v2N4213fuo9hkYqCVMtoop5Vy2iijlXIylCX3x/It1nfcT+f9Yy+5lSlTdzzBctqKLWz/OMZ4SM59s0m6kj+y7Wfasl2QbjyWts1r+cbB93L73xeyoamViSMGcM2b9+W8oyZSU+m4oXyLMdKwZjNzFq3JnnpsNc8tXkcm08Yj1Z+kilauHHErh00Z3XHqscntrd/aasKxJFBv2TK9vmnLX/v3YC33VH2ZieENrhv4eRaMPqWjVXrc8BomZLt4jxs2gAFVHvO7q6Uts0VQ35gN4bnXNzU2sX/93ZxY/wMGta1lftX+/HDw1dTGAzuf29jacUaAnhhQWd5lHHt5x+RzHV2duwu22wi7rb3QYttVWSAnoJZT3V2YLc8NtTCldQEHbp7D/hvrmLLxSaqyLXZN5YNYMmIGy/Y4hhVjjqVp2L5UVZZvOyRvsd3yjut90TujELS3nLaH8tb2gJ4T3lu7hPa2nEDfHvwzHT8mdG4nk4nElkYqN9ZTs24RNRsWMXDj6wzMLgdtqqci07RVTZsqR7F2wATW1kxgdXVyWVU9gVWV41lbPjL5kaHLDwqdPzS0/6ARKQtJCziB7PWkHXDr+5JW8faWxI77Q+hoUSwLyVCAECPVmU0MyGxgQNt6ato2bHlpXUdN23qqWzdQ3b5sXU916zqqWtdT1bqR0HXG+G3IhApaKofSUjmElqphtFQOobVyKK1V7cvk0lY1lLbqobRVZpdVw2irHkqoqNmi/tDldeW+/vbbWy67vifZ96PrY9ly2+3f29t/7IjEjqbSLe7ruN7ZetrRsprT6rrVYzrWd91Gl+d02c+2ttv+Odjyx5m49XZzH5NpI2SaKGttIrQ1UtbWRGhtIrQ1EVobk3VtTZS1NhLastc7Lo2E1ibKMk2UtTVTln1+xyXTRPm2bmeaKG+fwHE3ZEIFMZQTy8qTZSjvcl9Fx/2dj6nIWVdGLMt9THKdUE4m+xhy7m/fDh3Pz9lWWVnHbTr20/78zu2Ts2/Kkm1QtuW+c28n26jI3s4+p6wc2p8XygnZ+QFyv2Lm9ibY+l626F0AcNC4oUWRF4opbI8BXgFG2Y28uDUsW8647+7PHzMzuab5/7HfmMF85JR9Ofuw8UlLhPpMY0sbzzSsJfPI9bzptW/zhfJPc+fGzkmkRg2q4shJw7Pjv5Pzfg+sKuxf3HdFS1uGpWsbt+ravTMTjlVXlHUE5vHDO1umJw5qZeZDl1Cz4hk461tw9OV9+6K0fZvXwCPfgH98L5kB++B3w2nXJWNfSeaO2NjU1hHS28eud9e6vqGxlQ3NW7bC565rD807E2wrywNVFUlA7Tb45tzOXV+5vZbe3Q22axZ1jrl+7WHYtCK5v6wyOcf1Pm+GfU6BcUc47rJYxAjrl3Y/Tnz1Atj4xtbPaZ+0batTme0Nw/fa8aRtMULLph2MVV6z9Vjl3EvcyR/BQjnUDNv6MmB4zu1tXR+WjC/2h2Z1lckkY+dbG5MJ61oboa05Gete1h46Kzpvd70eyjyu+qGiCdsAIYSHgNtyJkj7TIzxmB09z7BdGF5dvoHvPjSPZU/9kdsrv8atAz/IuDP+mbcdvCdl/aQlo2BtXAHXH0SccBT1776XOYtWM3dR0vX8+cXrOsJCeVngwLFDsuO+k/Hfk0YWdut3JhN5Y0PTFrN27+yEYxVlgT2H1myza/e4YTWMHFS19etvbYY7z4fXHoSTPw8nf65vXqx6buU8+OO/wEu/g/JqOPYjycSNNUPzsvkYk27eybjYIvkxcdMqmP9IEq7nPwyrXutcN/bQzhnDJx0LVYNSKlK9qml956RtXU9jtvb17idtGzYxCeDDJkLL5u4n+drqedsSks/g9kLxFsG5y6VqsKFGUkEotrB9AHAbMApYB1waY3xuR88zbKfr2Ya1fOfBV/n9c0uJEf5j9J94//pbiZf9jjDlhLTLU7tfXgNP3QXXPA5jp3Xcvbk5af1Oxn4nk6+t2NDZ/XCPwVUcsVdn+D584vA+6xIdY2T1ppZuA3T7qbCWrWvstjtuCDB6cHXHrN15m3Ask0lO7/XMz+Goy5JWbb/0Fb7XHk7Oz73smeT8x2/5Ehz5gf5xKpyWzbDo752t10ueoqOT5/BJSav1PifD3ifBoD3Sq1OFoa0V1tV3fxqz1Qs6T+dWPbSbIDx8J1qch0HVkNI59Zukfq2owvauMmynY9b8VXznwVd5+OWkO9pbDxzDR07Zl6Oe+Di89AB8/nVbRQpJfR3c/JYkIJ79v9t8WIyR+tWbtwjfzy9Z1zHZT3lZ4KBxQzrGfU+fNIK9Rg7Ypdbv9gnHGtYkY6K7du1esnYzjS3ddyvsOuFYbhfvccNqem/CsT9+Cf72bTjgTDj/x3arLSaZtmTixr/+W9KNds9p8Pb/SLpJl5JMGyx5sjNcL/oHtGV/QBswIjkVV3vr9ci9UytTRShGaN7oOZslKcuwrbyKMfLQy29w44OvMnvBakKAMw8dx0dOnsrB47PdMr95IAzcAz78WLrFamvfPxneeAk+/ULS0rCTNje38XT9mmTytUXJ5GsrNnSOcd5jcDVHtp92bNJwDps4nBDYbtfu7iYcaze4uiInRHcG6PFpTzj2txvgj1+Evd4El/wq+cKp4tO4Dh67Hp64MQmhB5wJb/s3GLVv2pXtmhiTruCvPZjtGv5I0qUXkvPbTj6uM1zveagtipIk5YlhW3nRlon84bmlfOfBV3lu8ToqywPnHDmRD715H/YZPbjzgWsb4H8OhumXwjv/L72C1b25d8CvPgKnfx2O+fAubybGyOursq3f2csLS9Z3tH6XBdjWRMvVFWVbhOatxksPr2FoTeUu19Zrnrkb7rkC9jgAPvh7GDgy7Yq0u1YvgD99GZ6/L5kMbObV8OZ/Slp/C92G5UnX+PbW63XZ/wNDGYyfnp3U7GSYOBMqd/90T5IkaWuGbe2WlrYM981t4LsPz+O1NzZSU1nGhTMncdWJ+zB+eDetes//Cn5+Cbzz2zD9kr4vWNvXshmuPzgJEx+rzWsL16bmVp6uT8Z+P/36WqoqyraacGz88AGMGFhZ0BOudWveg3DHe5Oxvlf8MZmZV6Vj4RPwh8/D4rkwYCSc8gU46vLCGiLQtB4W/q0zYC/Pmc5k1H6dLddTTuhRrxVJkrTrDNvaJY0tbfxs9ut8/5HXaFizmSHVFVxy3GQ+ePzejBq8ndN//Olf4fH/hQ8/AXse3HcFa+e1/40uvhemvjXtagrfkqfg1jOSU3588Pce16Uqk4GnfwZ/uQ7WL4HRB8Lb/h32OzWdetpaoKGus+W6fnbnTM+D98wZd/3mZHZoSZLU5wzb6pF1jS3c/veF/PCx+azY0MyoQVV88IS9+cCxk3eua++tZyYT83xukZOnFKrVC+B/j4AD3gEX3pV2NYVt1Xy45W3J+NcP/BKmHJ92ReptzRuTH6Me/z9o3QxTT4O3/zuMPqB39xsjLH+h83RcCx6D5g3JuqrBSYt1e+v16AOdAV+SpAKwvbBdQP3jlLZVG5u59fH53Pa3BaxvbGX8sBq+cvbBXDBj0s5PRJVpS7phTphu0C5kI6bA/qcnM8avXggjJqddUWHauAJuPyeZtfr8Hxu0+4uqQUk38umXwF++mrR2z/srHP3B5Jzqg0blb19r6zu7hc9/GDYsS+4vq0jGWrePu55wFJQX4DwGkiRpmwzbYsnazfzgkfncNWsRm1va2HuPQfzLmQfz7iMn9Py0SctfgJaNyRdDFbaZV8LLD0DtD+G069KupvA0bUjGaK96Dc78Jhz8zrQrUl8bNhHO+T7M/BD8/nMw+wfJudXf/M8w4yqoqOr5NjevgQWPdgbsla90rhtzCEw7LwnXk4+D6sHb2IgkSSoGhu1+bMGKjdz08DzumVNPS1vkoHFD+egp+/KOaeMoL9vF7okNdcly4tH5K1S9Y5+3wMh9YM6Pk9Y6Zyvu1NYCv7gMFs+Bk/4JZlyZdkVK08Sjkknxnr0H/vwV+MMXYPYt8LavJUMxttedu6UR6md1jrtePBdi9rzxQyfCERcn4Xrvk2DInr3/WiRJUp8xbPdDLyxZx3cfmsdvn15MJsJRk0fwsVOmcvIBo3d/huiG2mQ5wbBd8MrKkta5P3wenvslHHFh2hUVhhjh15+AV/8ER14Mp3wx7YpUCEKAQ8+DA8+EJ74Dj14PP70wCclv/08YOy15XCYDS5/uDNeL/p6M+waoGZY8f5+TYZ9Tkh+7HHctSVLJcoK0fmTOotXc+OCr/PmF5QCcuN8efPSUqbxp75H5Ow3TjcfB5tVw7Qv52Z561+Y1cP1ByWRLVz+YdjWF4c/XwWPXJ2PaL7ijsE79pMKxfin89d+S89aHAIddkJxWb/4jsHlV8pjyaph0TOeM4eOOcC4LSZJKjBOk9WMxRh5/dSXfefBVnnhtJQCnHzKWj5yyL4dNHJ7fnTVtgDdeSFpuVBwGDIdD3wtzfpQMAejvY+3/8b0kaE84Gs77oUFb2zZkLLzrO9neIV+Ep+4CAow7PJlYbZ+Tk6BdOSDtSiVJUkr8JlmiMpnIn19YxncemsdTr6+hvCxwzpET+PDJ+7LfnkN6Z6dLnkzGItqFvLjMvCoJ27Nuhvf047D93C/hgX+GUVPh/T9PZqSWdmT8EXDZb2HZczB0PAwcmXZFkiSpQBi2S0xrW4bfPr2EGx96lZeXbaCqooyLj5nEh07al71GDuzdnddnx2s7OVpxGXsoTDo2mfzpbV/L72mNisX8R+Heq2HwGLj43v75HmjXhdA5ZluSJCnLsF0iGlvauGdOPTc9PI/XV21mUFU5HzppH644YW/GDO2jWaYbaiGUJeMSVVxmXAmLroC5P4YTPpV2NX1r6bPw0/dDRQ1cfI/nHJckSVJeGLaL3MamVu78xyJ+8OhrLF/fxPCBlXzq1P259LjJDB+4C+eA3R31dTD6IM8NW4wOeicMGgOzfwjHfaL/TOK0ZhHcfi60NiZBe+yhaVckSZKkEmHYLlJrNjXzo78t5Na/zWfNphbGDKnmi2ccxPvfNIlB1Sn8WdcthvWLYb9T+37f2n0VVXD05fDwf8HLf4ADz0i7ot63aRX85BzYsCyZDG3vk9KuSJIkSSXEsF1klq9r5JbH5nP73xeysbmNvUYO4J/efgDnTp9ITWWKrZENdcnSydGK11GXwSP/DbN/UPphu3kT3Hk+rHwFTv8vmHZO2hVJkiSpxBi2i8TrqzbxvUfm8fPaeppbM+w3ZjAfPWUqZx02jorysrTL65wcrb+fOqqYDR0PB50Fz/8KVrwKe0xNu6Le0dYKd18O9bPh+E/CMdekXZEkSZJKkGG7wL26fD03PjiPXz21mLZM5PCJw/jIKVM57aA9KSsLaZfXqaEOKgfBmIPSrkS7Y+bVSdiuvQVO/8+0q8m/GOG3n4SXfw+HXwinfiXtiiRJklSiDNsF6un6Ndz44Dz+8PxSYoRj9xnFR0+ZyvFTRxFCAYVsgEwbLJ4L44/sPxNrlarJxyeT3M29A97ypdI71/SD/wFzfwL7vhXe+e3klE2SJElSLzBsF5AYI/+Yv4rvPPgqj76yAoBTDxrDh0+eylGTR6Rc3Xa88RI0b4CJdiEveiHAzCvhd9fC0z9PJk0rFbNvgUf+v+RHofN/DOWVaVckSZKkEmbYLgAxRh566Q2+8+Cr1C5cTVmAsw8fz0dO3peDxg1Nu7wda3C8dkk57H3w5+tg9s3JpGml0Pr7wm/g/s/AyH3g/b/w9HSSJEnqdYbtFLVlIg88u4TvPDiPF5aso7I88L4Ze/GhN+/L3nsUUffdjsnRnIm8JFQPTsYzz/oeLHoCJh+XdkW7Z+ETcPcVMHAUXHwvDB6ddkWSJEnqBwzbKWhuzXDf3Aa++/A85q/YSE1lGR88fm+uOmlvxg0bkHZ5PdcwB4aMg2ET0q5E+TLjyiRsz/pBcYft5S/AXRckXcYvuhtG7p12RZIkSeonDNt9aHNzGz+dvYjvP/IaS9Y2MqSmgo+dMpXLj5/CqMHVaZe3a5o3wvLn4IASPy9zfzN6f9j7zfDCr2H9UhgyNu2Kem5tPdx+bnJO7Yt+DuOPSLsiSZIk9SOG7T6yamMzp13/MCs3NjNqUBWfPf0ALj5mMkNrinySpsVPQsw4XrsUzbwa5j8MdbfByZ9Lu5qe2bw6CdrrGuCcm2Hft6RdkSRJkvoZw3YfGTmoirceNIaDxw3lghmTGFBVIqfIap8cbaLjtUvO/qfD0IlQeyuceG3xzN7dshnuuhDeeBHe9jU47L1pVyRJkqR+qCztAvqT/++8w7ns+L1LJ2gDNNQBITmdkkpLeUVy6q8NS+HF36Zdzc7JtME9VyYTux37MTju42lXJEmSpH7KsK3dU18HYw6C6iFpV6LeMP1SKK+CWTenXcmOxZicH/zF38K08+C0f0u7IkmSJPVjhm3tuvVLYV09TJiediXqLYNHwyHvgYWPwbLn0q5m+x75BtTdCvucDO/+LpT5z5skSZLS47dR7TrPr90/zLgqWc4u4Nbtuh/Bg/8OYw+D838CFVVpVyRJkqR+zrCtXddQlyydHK20TTwaxh0BT/0MGtemXc3WXrwffvtJGD45OZd2zdC0K5IkSZIM29oNDbVQORBGH5R2JepNIcDMq6BlIzx5V9rVbGnRP+Duy2HACPjAL2HInmlXJEmSJAGGbe2qTBs0zE1aPMs9g1zJm3ZuEmhn35xMRFYI3ngJ7roAQhlc9AsYtW/aFUmSJEkdDNvaNStehub1MPGotCtRX6gcAEdeDCtfgdceSrsaWLcEbj8XmtYnY7QneBxKkiSpsBi2tWvax2s7OVr/cfQVQIBZP0i3js1rkqC99nV45w2w36np1iNJkiR1w7CtXdM+E7mTo/UfI/eG/d4GLz8Aa15Pp4aWRvjpRbD8OTj1K3DEhenUIUmSJO2AYVu7pqEWBu8JQyekXYn60syrIGag9od9v+9MG/zy6uSc32+6Bo7/ZN/XIEmSJO0kw7Z6rnkTLHs+6UIeQtrVqC/t+1YYsTfM+VHSytxXYoTffw6e/xUc8h54+3967EmSJKmgGbbVc0uegtjm5Gj9UVkZzLgSNq2E5+/ru/0+9j8w6/sw5UR4z/eSOiRJkqQC5jdW9VxDdry2M0D3T0deBBUD+m6itLl3wF+ugz2nwfvugIrqvtmvJEmStBsM2+q5+logwPjpaVeiNAwYAYe9N/nRpWFO7+7rlT/Brz8OwybBRXdDzbDe3Z8kSZKUJ4Zt9VxDHYw+AGqGpl2J0jLjqmQ5++be20d9Hfz8kuQ4u/geGDqu9/YlSZIk5ZlhWz2zfllyfmPPr92/jTsM9noTPHsPbFqV/+2veBXufG8yMdr7fw6j98//PiRJkqReZNhWzzTUJcsJdiHv92ZcBa2NMPcn+d3u+mVw+3tg8xp4722w18z8bl+SJEnqA4Zt9Uz75GgTbdnu9w5+Fwwak3Qlz7TlZ5uN6+COc2HNIjj7W3DA6fnZriRJktTHDNvqmfraZCbqMYekXYnSVlEFR12aBONX/rT722ttgp9dDEufgVO+BNMv2f1tSpIkSSnpUdgOIQwKIbwthPDBEMKFIYSDe6swFaBMBhbPhfFHQHlF2tWoEBx1OYRymL2bpwHLZOC+D8P8h+HoK+Ckz+SnPkmSJCklOxW2QwiTQwi3Aq8B/wy8FTgPuDuE8FwI4fJerFGFYuUr0LTO82ur07AJcOCZ8OqfYeW8XdtGjPDHLyaTrR10NpzxDQghv3VKkiRJfWxnW7Z/AvwKmBBjfGuM8aIY47kxxoOBM4H9Qwgf67UqVRjqs+O1DdvKNbP9NGC37Nrz//Zt+PuNMOk4OOdmKCvPX22SJElSSnaqL3CM8aTtrFsAfD5fBamAOTmaujPlRBh9IDx5O7zlS1A1cOef+/TP4U//AqMPggvvhMqa3qtTkiRJ6kM9HbN9XG8VoiLQUJfMPj1sr7QrUSEJAWZcCY1r4Zlf7PzzXv1LMk576AS4+B4YMKL3apQkSZL6WE9nI/9hCOEtuXeEEN6Zx3pUqFo2w7Lnki7kjqdVV4e/D6qGwKwfJGOwd2TxXPj5JVA1GC6+Nxn7LUmSJJWQnobttwP/G0I4K4RwRAjhr8BX81VMCOE/QggvhBCeCiHM6hrslaIlT0GmFSY6XlvdqB6SBO5lz8Dr/9j+Y1e9Bne8Nzme3v8zGHNg39QoSZIk9aEehe0Y40LgA8DPgPuB24Aj81jPo8D0GOPhwFXAPSEEB3EWgo7J0RyvrW2YcWWynLWd04BteAN+cg5sWgnn3gKTjumb2iRJkqQ+1tMx2/8B/B74b2AZsDnGnekzunNijA/EGDdnbz4DlAN75Gv72g0NdclywvR061DhGnMg7H0SPP8rWL9s6/VNG+CO82D1fDjzm3DQWX1foyRJktRHetqNfAxwRIzxy8BpwGd68RzblwPzYoz1vbR99URDLeyxP9QMS7sSFbIZV0GmBeb8aMv7W5vh5x+AJU/Cmz8HR38wlfIkSZKkvtLTbuRXxhiXZq+vAN5GEop3Sgjh0RDCim1c9sp53FuBLwPv2862Ph1CqG+/bNiwoScvRT2x4Q1Ys8gu5NqxA85IZhevvRXaWpP7Mhn49cdg3l9h+qVw8ufSrVGSJEnqAzsVtkMIb+ru/hjjWuDtIYQBIYSDd7SdGOOJMcY9tnF5PbuvNwO3AmfHGF/azraujzFObL8MHjx4Z16KdkXH+bWdHE07UF4BR18O6xfDS79L7vvLV+Dpn8H+74Azr3c2e0mSJPULO9uy/dkQwp9DCJeHEA4IIQwLIYwNIZwMfAX4BzB2d4sJIZwE/AR4V4zxqd3dnvKkY7y2LdvaCdMvhbLKZKK0J26Ex/8XJs6E836YhHFJkiSpH9ipb74xxnNDCDOADwFfBCYCG4GngXuA42KM+ejHfQtQDdwaOlu/PhBjfCYP29auqq+FihrY85C0K1ExGDwGDnk3PPMLWPBoMtb//T+DqoFpVyZJkiT1mZ1uZooxzgZm92ItxBj3683taxdkMtAwB8YdDuWVaVejYjHz6iRsDxkHF98DA0emXZEkSZLUp+zTqe1b+So0rbULuXpmr5nwvjthz2kwfFLa1UiSJEl9bmcnSDustwtRgWofr+3kaOqpA8+EEZPTrkKSJElKxc5OkPZDgBDCh3qxFhWi9pnIJxi2JUmSJGln7WzYrgohjCeZIE39SX0tDNwDhttCKUmSJEk7a2fHbP8XUAcMDiHcmb1eB8yJMa7rreKUspbNsOxZmHqq50aWJEmSpB7YqZbtGOMdwARgMfArknNq/yuwMITwSu+Vp1QtfQYyrU6OJkmSJEk91JNTf2VCCMfEGFcDP2u/P4Tg6bpKVX37eO3p6dYhSZIkSUVmZ8dsA5AN2l3vs2W7VDk5miRJkiTtkh6FbfUz9bUwaj8YMDztSiRJkiSpqBi21b2NK2DNQpjoeG1JkiRJ6inDtrrXUJcs7UIuSZIkST1m2Fb36h2vLUmSJEm7yrCt7jXUQnk17Dkt7UokSZIkqegYtrW1GJNu5OMOh4qqtKuRJEmSpKJj2NbWVs6DxrV2IZckSZKkXWTY1tbaz6/tTOSSJEmStEsM29qak6NJkiRJ0m4xbGtrDXUwcBSMmJJ2JZIkSZJUlAzb2lJLIyx9JmnVDiHtaiRJkiSpKBm2taWlz0CmBSY4XluSJEmSdpVhW1vqmBzN8dqSJEmStKsM29pSQ12yHD893TokSZIkqYgZtrWl+loYuS8MHJl2JZIkSZJUtAzb6rRxJaye7/m1JUmSJGk3GbbVafGcZOnkaJIkSZK0Wwzb6lTv5GiSJEmSlA+GbXVqqIXyKthzWtqVSJIkSVJRM2wrEWMyE/nYw6CiOu1qJEmSJKmoGbaVWPUabF7t5GiSJEmSlAeGbSXaz6/t5GiSJEmStNsM20q0T442YXq6dUiSJElSCTBsK9FQCwNGwsh90q5EkiRJkoqeYVvQ2gRLn4EJR0EIaVcjSZIkSUXPsC1Y+iy0NTs5miRJkiTliWFbSRdySFq2JUmSJEm7zbCtnMnRDNuSJEmSlA+GbSWn/Rq5DwwcmXYlkiRJklQSDNv93aZVsGqe59eWJEmSpDwybPd3DXOSpV3IJUmSJClvDNv9XfvkaM5ELkmSJEl5Y9ju7xrqoLwKxh6adiWSJEmSVDIM2/1ZjMlM5GMPhYrqtKuRJEmSpJJh2O7PVs+Hzascry1JkiRJeWbY7s/q65KlM5FLkiRJUl4ZtvuzhmzYdnI0SZIkScorw3Z/1lALNcNh5D5pVyJJkiRJJcWw3V+1NsOSp5Px2iGkXY0kSZIklRTDdn+17Bloa7ILuSRJkiT1AsN2f9UwJ1k6OZokSZIk5Z1hu7+qr02WnvZLkiRJkvLOsN1fNdTCiCkwaFTalUiSJElSyTFs90ebV8PKV+1CLkmSJEm9xLDdH7WP13ZyNEmSJEnqFYbt/qihLlk6XluSJEmSekVBhu0QwskhhLYQwsfSrqUk1ddCWSWMPSztSiRJkiSpJBVc2A4hDAH+C3gg7VpKUozJ5Ghjp0FlTdrVSJIkSVJJKriwDVwPfANYkXYhJWnNQti00snRJEmSJKkXFVTYDiG8AxgeY7w77VpKlufXliRJkqReV9GXOwshPAoctI3VRwJfB07byW19Gvh0++1hw4btdn39QvvkaM5ELkmSJEm9pk/DdozxxG2tCyGcAIwDZoUQAPYAzg4hjI4xfrmbbV1P0uUcgIkTJ8b8V1yC6muhZhiM3DftSiRJkiSpZPVp2N6eGONjwJj22yGE24DaGOMNqRVVatpaYMlTMOV4KCuoEQSSJEmSVFJMXP3Jsmehrcnx2pIkSZLUywqmZburGONladdQcjomR3O8tiRJkiT1Jlu2+xMnR5MkSZKkPmHY7k8a6mD4ZBi0R9qVSJIkSVJJM2z3F5vXwIqXHa8tSZIkSX3AsN1fLJ6TLO1CLkmSJEm9zrDdX9Rnx2s7OZokSZIk9TrDdn/RUAdlFTDusLQrkSRJkqSSZ9juD2KEhlrY8xCoHJB2NZIkSZJU8gzb/cGaRbDxDbuQS5IkSVIfMWz3Bw21ydLJ0SRJkiSpTxi2+4OG7EzktmxLkiRJUp8wbPcH9bVQPQxGTU27EkmSJEnqFwzbpa6tBZY8CROOhDL/3JIkSZLUF0xfpW7Zc9DaaBdySZIkSepDhu1S11CXLJ0cTZIkSZL6jGG71LWH7QlHpVuHJEmSJPUjhu1SV18LwybB4DFpVyJJkiRJ/YZhu5Q1roUVL8NEW7UlSZIkqS8ZtkvZ4rlAdHI0SZIkSepjhu1SVl+bLB2vLUmSJEl9yrBdyhrqIJTDuMPTrkSSJEmS+hXDdqmKMWnZ3vMQqBqYdjWSJEmS1K8YtkvV2nrYuNzza0uSJElSCgzbparB8dqSJEmSlBbDdqnqmBzNlm1JkiRJ6muG7VLVUAfVQ2GP/dOuRJIkSZL6HcN2KWprhcVPwvgjocw/sSRJkiT1NZNYKVr+PLRudry2JEmSJKXEsF2K2idHcyZySZIkSUqFYbsU1dclSydHkyRJkqRUGLZLUUMdDNsLhuyZdiWSJEmS1C8ZtktN4zp440WYMD3tSiRJkiSp3zJsl5rFc4FoF3JJkiRJSpFhu9Q4OZokSZIkpc6wXWoa5kAoh3GHp12JJEmSJPVbhu1SEiPU18KYg6FqUNrVSJIkSVK/ZdguJesaYMNSmHhU2pVIkiRJUr9m2C4l9dnx2k6OJkmSJEmpMmyXkoa6ZOnkaJIkSZKUKsN2KWmog6rBsMf+aVciSZIkSf2aYbtUtLUm59gefySUladdjSRJkiT1a4btUvHGC9CyyS7kkiRJklQADNulon28tpOjSZIkSVLqDNulomMmck/7JUmSJElpM2yXioY6GDoBho5LuxJJkiRJ6vcM26WgaT0sf8FWbUmSJEkqEIbtUrD4SSA6OZokSZIkFQjDdilocLy2JEmSJBUSw3YpqK+FUAbjjki7EkmSJEkShu3S0FAHYw6G6sFpVyJJkiRJwrBd/NYthvVL7EIuSZIkSQXEsF3sPL+2JEmSJBUcw3axa58czZnIJUmSJKlgGLaLXX0dVA2G0QemXYkkSZIkKcuwXcwybbB4Low/EsrK065GkiRJkpRVcGE7hPCREMILIYRnQwhPhxBq0q6pYL3xIrRshAnT065EkiRJkpSjIu0CcoUQ3gVcBBwTY1wbQhgDtKRcVuHqmBzN8dqSJEmSVEgKKmwD/wR8Oca4FiDGuDzlegqbk6NJkiRJUkEqtG7kBwNHhxAeDyHUhhA+kXZBBa1hDgwZD0PHp12JJEmSJClHn7ZshxAeBQ7axuojSerZFzgJGAY8HEJ4NcZ4fzfb+jTw6fbbw4YNy3/BhaxpAyx/Hg44I+1KJEmSJEld9GnYjjGeuL31IYRFwF0xxjZgVQjhAWAmsFXYjjFeD1zffnvixIkxz+UWtiVPQszYhVySJEmSClChdSO/EzgdIDsL+ZuBp1KtqFA5OZokSZIkFaxCC9v/A4wNITwP1AEPxBh/mXJNhamhDkJZco5tSZIkSVJBKajZyGOMm4FL0q6jKDTUwegDoXpw2pVIkiRJkrootJZt7Yx1S2BdA0w4Ku1KJEmSJEndMGwXI8+vLUmSJEkFzbBdjBrqkqUt25IkSZJUkAzbxai+FioHwuhtnbJckiRJkpQmw3axybTB4rnJLOTlBTW/nSRJkiQpy7BdbN54CZo32IVckiRJkgqYYbvYOF5bkiRJkgqeYbvYOBO5JEmSJBU8w3axqa+DwWNh6IS0K5EkSZIkbYNhu5g0b4TlzyWt2iGkXY0kSZIkaRsM28VkyVMQMzBhetqVSJIkSZK2w7BdTOqz47UnOF5bkiRJkgqZYbuYNNQCITnHtiRJkiSpYBm2i0l9HYw+EGqGpl2JJEmSJGk7DNvFYv1SWFfv+bUlSZIkqQgYtotFQ12ynGjYliRJkqRCZ9guFk6OJkmSJElFw7BdLBpqoWIAjDk47UokSZIkSTtg2C4GmQw0zIXxR0B5RdrVSJIkSZJ2wLBdDFa8DM3rnRxNkiRJkoqEYbsYNGTHa090vLYkSZIkFQPDdjFwcjRJkiRJKiqG7WLQUAeDxsCwiWlXIkmSJEnaCYbtQte8CZY9l3QhDyHtaiRJkiRJO8GwXeiWPAWxzcnRJEmSJKmIGLYLnZOjSZIkSVLRMWwXuoY6IMD4I9OuRJIkSZK0kwzbha6+DvbYH2qGpV2JJEmSJGknGbYL2YblsHaRXcglSZIkqcgYtgtZx/m1nRxNkiRJkoqJYbuQNdQlS8O2JEmSJBUVw3Yha6iFihrY85C0K5EkSZIk9YBhu1BlMtAwB8YdAeWVaVcjSZIkSeoBw3ahWvkKNK1zcjRJkiRJKkKG7ULVMV57erp1SJIkSZJ6zLBdqDpmIrdlW5IkSZKKjWG7UDXUwqDRMHxS2pVIkiRJknrIsF2IWjbDsueSVu0Q0q5GkiRJktRDhu1CtORpyLR6fm1JkiRJKlKG7ULUkB2vPdGwLUmSJEnFyLBdiNonRxvvTOSSJEmSVIwM24WooRb22B8GDE+7EkmSJEnSLjBsF5oNb8CaRY7XliRJkqQiZtguNA11ydKwLUmSJElFy7BdaDomRzs63TokSZIkSbvMsF1o6muhvBr2nJZ2JZIkSZKkXWTYLiSZDCyeA+MOh/LKtKuRJEmSJO0iw3YhWTUPGtfahVySJEmSipxhu5C0n1/bydEkSZIkqagZtguJk6NJkiRJUkkwbBeShjoYOAqGT067EkmSJEnSbjBsF4qWRlj6LEw4GkJIuxpJkiRJ0m4wbBeKpU9DpsUu5JIkSZJUAgzbhaKhLlk6OZokSZIkFb2CCtshhH1DCH8JITwZQngxhPDNEEJB1dhrOmYin55uHZIkSZKk3VZoQfa/gV/FGI8AjgDeBpyeZkF9pqEWRk2FASPSrkSSJEmStJsKLWwDDMsuBwCVwJIUa+kbG1fA6gXJ5GiSJEmSpKJXaGH7k8B7QwiLgcXAj2OMc7t7YAjh0yGE+vbLhg0b+rLO/GqYkyydHE2SJEmSSkKfhu0QwqMhhBXbuOwFfAj4SYxxPDAZeH8I4S3dbSvGeH2McWL7ZfDgwX35UvKrwfHakiRJklRKKvpyZzHGE7e3PoTwCWCf7GOXhxAeAN4M/LUPyktPfS2UV8Oeh6ZdiSRJkiQpDwqtG/lrwDsAQgiDgLcAz6ZaUW+LMTnt17jDoKIq7WokSZIkSXlQaGH7UuDqEMLTQC3wR+DudEvqZateg8Y1To4mSZIkSSWkT7uR70h2MrTj066jT3WcX/uodOuQJEmSJOVNobVs9z/tk6NNNGxLkiRJUqkwbKetvhYGjoIRe6ddiSRJkiQpTwzbaWptgqXPJF3IQ0i7GkmSJElSnhi207T0Gci0OF5bkiRJkkqMYTtNHZOjORO5JEmSJJUSw3aa2idHmzA93TokSZIkSXll2E5TQx2M3BcGjky7EkmSJElSHhm207JpFax6zfHakiRJklSCDNtpaahLlhMdry1JkiRJpcawnRYnR5MkSZKkkmXYTktDHZRXwdhpaVciSZIkScozw3YaYkzC9thDoaI67WokSZIkSXlm2E7Dqtdg8yq7kEuSJElSiTJsp8HJ0SRJkiSppBm209Aetj3tlyRJkiSVJMN2GuprYcAIGLlP2pVIkiRJknqBYbuvtTbB0qeTVu0Q0q5GkiRJktQLDNt9bemz0Nbs5GiSJEmSVMIM233N8dqSJEmSVPIM232toTZZGrYlSZIkqWQZtvtafS2M2BsGjUq7EkmSJElSLzFs96VNq2DVPM+vLUmSJEklzrDdlxbPSZZ2IZckSZKkkmbY7kv17ZOj2bItSZIkSaXMsN2XGmqhrBLGHpp2JZIkSZKkXmTY7isxJpOjjT0UKmvSrkaSJEmS1Isq0i6g32hthEPeAyOmpF2JJEmSJKmXGbb7SuUAOOv6tKuQJEmSJPUBu5FLkiRJkpRnhm1JkiRJkvLMsC1JkiRJUp4ZtiVJkiRJyjPDtiRJkiRJeWbYliRJkiQpzwzbkiRJkiTlmWFbkiRJkqQ8M2xLkiRJkpRnhm1JkiRJkvLMsC1JkiRJUp4ZtiVJkiRJyjPDtiRJkiRJeWbYliRJkiQpzwzbkiRJkiTlmWFbkiRJkqQ8M2xLkiRJkpRnhm1JkiRJkvLMsC1JkiRJUp4ZtiVJkiRJyrMQY0y7hrwIITQBb6Rdx04YDGxIuwipj3i8qz/xeFd/4vGu/sJjXTsyOsZY3d2KkgnbxSKEUB9jnJh2HVJf8HhXf+Lxrv7E4139hce6dofdyCVJkiRJyjPDtiRJkiRJeWbY7nvXp12A1Ic83tWfeLyrP/F4V3/hsa5d5phtSZIkSZLyzJZtSZIkSZLyzLAtSZIkSVKeGbb7SAhhvxDC30IIL4cQZoUQDk67Jqk3hBBqQgj3ZY/1J0MIvw8hTEm7Lqm3hRC+HEKIIYRpadci9ZYQQnUI4YYQwishhOdCCLenXZPUW0IIbw8h1IUQ5oYQng0hXJp2TSouFWkX0I98D/h+jPG2EMJ5wC3AsSnXJPWW7wMPxBhjCOFj2dtvS7kmqdeEEKYDxwCL0q5F6mVfBzLA/tl/48elXZDUG0IIAbgTOCXG+HS24eDFEMK9Mcb16VanYmHLdh8IIYwBpgPtv/7eA+xta59KUYyxMcZ4f+ycffHvwD5p1iT1phBCNfAd4COAs46qZIUQBgGXA19o/zc+xrgk3aqkXjc8uxwKrASa0itFxcaw3Tf2AhbHGFsBsv9BLQImpVqV1Dc+Afwm7SKkXvRV4PYY4/y0C5F62b4kYeNLIYTaEMKjIYS3pl2U1Buy39fPB+4NISwEHgMujTE2p1uZiolhu+90be0IqVQh9aEQwheA/YAvpl2L1BtCCMcCM4Ab065F6gOVJD2Vno8xHg18DPhpCGF0umVJ+RdCqAA+D7wrxjgZeCvwoxDCyHQrUzExbPeN14GJ2Q9t+xiQvXBsn0pYCOEzwDnAO2KMm9KuR+olbwYOBOaHEBYAE4E/hBDekWpVUu9YSDJe+w6AGONTwHzgkDSLknrJEcD4GOPjADHG2cBi4PA0i1JxMWz3gRjjcmAucHH2rnOBBTHGBakVJfWiEMKngQuB02KMa1IuR+o1McavxxjHxxinxBinAPXA22OMD6RcmpR3McYVwF+AtwOEECYDewMvpVmX1EvaG8sOAAghTCUZSvFyqlWpqITOOYzUm7If1NuAUcA6kjEfz6ValNQLQggTSf6Deg1on62zKcb4pvSqkvpGtnX7rBjjs2nXIvWGEMI+wA9Jvs+0AdfFGH+ZblVS7wghXAh8gaRHRwD+I8b403SrUjExbEuSJEmSlGd2I5ckSZIkKc8M25IkSZIk5ZlhW5IkSZKkPDNsS5IkSZKUZ4ZtSZIkSZLyzLAtSZIkSVKeGbYlSZIkScozw7YkSf1QCOHAEEJ9CGGf7O1/CiHcH0IIadcmSVIpCDHGtGuQJEkpCCFcCFwLfAa4FZgRY1yRblWSJJWGirQLkCRJ6Ygx3hVCOAX4A/BWg7YkSfljN3JJkvqpEEIFMA1YBUxIuRxJkkqKYVuSpP7r68BLwEnAN0MIU1OuR5KkkmE3ckmS+qEQwlnA6cDMGOOmEMJngF+EEI6NMTamXJ4kSUXPCdIkSZIkScozu5FLkiRJkpRnhm1JkiRJkvLMsC1JkiRJUp4ZtiVJkiRJyjPDtiRJkiRJeWbYliRJkiQpzwzbkiRJkiTlmWFbkiRJkqQ8M2xLkiRJkpRn/z/ddnHhsbHaLwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1200x640 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "plt.figure(figsize=(15, 8), dpi=80)\n",
    "plt.plot(ks,f1s, label = r\"$f\\ (x_1)$\")\n",
    "plt.plot(ks,f2s, label = r\"$f\\ (x_2)$\")\n",
    "plt.ylabel(r\"$f\\ (x)$\")\n",
    "plt.xlabel('x')\n",
    "plt.title('Fixed-point iteration solving a system of nonlinear equations')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "88569896",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a02bae99",
   "metadata": {},
   "source": [
    "## 2. Newton's method  <a class=\"anchor\" id=\"first-bullet\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a4ce2e33",
   "metadata": {},
   "source": [
    "The generalization of Newton's method is\n",
    "\n",
    "$$\n",
    "x_{k+1} = x_k - J_F(x_k)^{-1} F(x_k), \\quad k = 0,1,2, ...\n",
    "$$\n",
    "\n",
    "Rewrite it to the standard form for a linear system:\n",
    "\n",
    "$$\n",
    "J_F(x_j) \\Delta x_k = - F(x_k), \\quad k = 0,1,2,...\n",
    "$$\n",
    "\n",
    "where $\\Delta x_k = x_{k+1}- x_k$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "87870d3a",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e28b4862",
   "metadata": {},
   "source": [
    "If $x_0$ is sufficiently close to $\\alpha$, then Newton's method converges quadratically in $n$ dimension. This relies on **Taylor's Theorem**.\n",
    "\n",
    "The **multivariate Taylor theorem** is: \n",
    "\n",
    "Let $\\phi(s) = F(x+ s \\delta)$, then 1D Taylor Theorem yields\n",
    "\n",
    "$$\n",
    "\\phi(1) = \\phi(0) + \\sum^k_{l=1} {\\phi^{(l)} (0) \\over l!} + \\phi^{(k+1)}(\\eta), \\quad \\eta \\in (0,1)\n",
    "$$\n",
    "\n",
    "We have\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\phi(0) &= F(x)\\\\\n",
    "\\phi(1) &= F(x+\\delta)\\\\\n",
    "\\phi'(s) &= {\\partial F(x + s\\delta) \\over \\partial x_1} \\delta_1 + {\\partial F(x + s\\delta) \\over \\partial x_2} \\delta_2 + ... + {\\partial F(x + s\\delta) \\over \\partial x_n} \\delta_n\\\\\n",
    "\\phi''(s) &= {\\partial^2 F(x + s\\delta) \\over \\partial x_1^2} \\delta_1^2 + ... + {\\partial^2 F(x + s\\delta) \\over \\partial x_1x_n} \\delta_1 \\delta_n + ... + {\\partial^2 F(x + s\\delta) \\over \\partial x_1x_n} \\delta_1\\delta_n + ... + {\\partial^2 F(x + s\\delta) \\over \\partial x_n^2}\\delta_n^2\\\\\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "where $x$ is a vector, $s$ is a scalar (telling how far along the vector $\\delta$), and $\\delta$ is a vector (considered as a direction).\n",
    "\n",
    "\n",
    "Hence, we have\n",
    "\n",
    "$$\n",
    "F(x + \\delta) = F(x) + \\sum^k_{l=1} {U_l(\\delta)\\over l!} + E_k\n",
    "$$\n",
    "\n",
    "where $U_l(x)$ is a direction derivative operator, evaluated at $x$.\n",
    "\n",
    "$$\n",
    "U_l(x) = \\left[\\left({\\partial \\over \\partial x_1} \\delta_1 + ... + {\\partial \\over \\partial x_n} \\delta_n\\right)^l F\\right]_{(x)}, \\quad l = 1,2,...,k,\n",
    "$$\n",
    "\n",
    "and the error, evaluated at $(x+ \\eta \\delta)$.\n",
    "\n",
    "$$\n",
    "E_l = (U_{k+1})_{(x+ \\eta \\delta)}, \\quad \\eta \\in (0,1)\n",
    "$$\n",
    "\n",
    "Let $A$ be an upper bound on the absolute value of all derivative of order $k+1$, we can derive the upper bound of $|E_k|$:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "|E_k| &\\leq {1\\over (k+1)!} |(A, ..., A)^T(||\\delta||^{k+1}_\\infty, ..., ||\\delta||^{k+1}_\\infty)|\\\\\n",
    "&= {1\\over (k+1)!} A ||\\delta||^{k+1}_\\infty |(1,...,1)^T (1,...,1)|\\\\\n",
    "&= {n^{k+1} \\over (k+1)!} A ||\\delta||_\\infty^{k+1}\n",
    "\\end{align}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d27f6295",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7f0c8c32",
   "metadata": {},
   "source": [
    "To analyse Newton's method, we only need to use the first order Taylor expansion:\n",
    "\n",
    "$$\n",
    "F(x+\\delta) = F(x) + \\nabla F(x)^T \\delta + E_1\n",
    "$$\n",
    "\n",
    "For each $F_i$, we have\n",
    "\n",
    "$$\n",
    "F_i(x+\\delta) = F_i(x) + \\nabla F_i(x)^T \\delta + E_{i,1}\n",
    "$$\n",
    "\n",
    "So for $F$ we have\n",
    "\n",
    "$$\n",
    "F(x+\\delta) = F(x) + J_F(x) \\delta + E_F\n",
    "$$\n",
    "\n",
    "where $||E_F||_\\infty \\leq \\max_{1 \\leq i \\leq n}|E_{i,1}| \\leq {1\\over 2} n^2 \\left(\\max_{1 \\leq i, j, l \\leq n} |{\\partial^2 F_i \\over \\partial x_j \\partial x_l}|\\right) ||\\delta||^2_{\\infty}$.\n",
    "\n",
    "\n",
    "Let's Taylor expand $F(\\alpha)$ at $F(x_k)$,\n",
    "\n",
    "$$\n",
    "0 = F(\\alpha) = F(x_k) + J_F(x_k)[\\alpha - x_k] + E_F\n",
    "$$\n",
    "\n",
    "So that\n",
    "\n",
    "$$\n",
    "x_k - \\alpha = [J_F(x_k)]^{-1} F(x_k) + [J_F(x_k)]^{-1} E_F\n",
    "$$\n",
    "\n",
    "Also the Newton iteration can be rewritten as\n",
    "\n",
    "$$\n",
    "J_F(x_k)[x_{k+1} - \\alpha] = J_F(x_k) [x_k - \\alpha] - F(x_k)\n",
    "$$\n",
    "\n",
    "Combine these two equations we have\n",
    "\n",
    "$$\n",
    "x_{k+1} - \\alpha = [J_F(x_k)]^{-1} E_F\n",
    "$$\n",
    "\n",
    "so that $||x_{k+1} - \\alpha||_{\\infty} \\leq const. ||x_k - \\alpha||^2_\\infty$, i.e. **quadratic convergence**."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9f81e521",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "72e3e22d",
   "metadata": {},
   "source": [
    "The following program implements a Newton's method for the two-point Gaussian quadrature rule.\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "F_1(x_1, x_2, w_1, w_2) &= w_1 + w_2 -2 = 0\\\\\n",
    "F_2(x_1, x_2, w_1, w_2) &= w_1x_1 + w_2x_2 = 0\\\\\n",
    "F_3(x_1, x_2, w_1, w_2) &= w_1x_1^2 + w_2 x_2^2 - 2/3 = 0\\\\\n",
    "F_4(x_1, x_2, w_1, w_2) &= w_1x_1^3 + w_2x_2^3 = 0\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "The Jacobian of the system is required\n",
    "\n",
    "$$\n",
    "J_F(x_1, x_2, w_1, w_2) = \\begin{bmatrix}0 & 0& 1 & 1 \\\\ w_1 & w_2 & x_1 & x_2 \\\\ 2w_1x_1 & 2w_2x_2 & x_1^2 & x_2^2\\\\ 3w_1x_1^2 & 3w_2x_2^2 & x_1^3 & x_2^3 \\end{bmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "fee8d85c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-1.  1.  1.  1.]\n",
      "[-0.66666667  0.66666667  1.          1.        ]\n",
      "[-0.58333333  0.58333333  1.          1.        ]\n",
      "[-0.57738095  0.57738095  1.          1.        ]\n",
      "[-0.57735027  0.57735027  1.          1.        ]\n",
      "[-0.57735027  0.57735027  1.          1.        ]\n",
      "\n",
      "Solution is [-0.57735027  0.57735027  1.          1.        ]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy.optimize import fsolve\n",
    "\n",
    "# Function f(x_1,x_2,w_1,w_2)\n",
    "def f(x):\n",
    "    print(x)\n",
    "    (x1,x2,w1,w2)=(x[0],x[1],x[2],x[3])\n",
    "    x1s=x1*x1;x2s=x2*x2\n",
    "    f=np.array([w1+w2-2,\n",
    "                w1*x1+w2*x2, \n",
    "                w1*x1s+w2*x2s-2./3,\n",
    "                w1*x1s*x1+w2*x2s*x2])\n",
    "    return f\n",
    "\n",
    "# Analytical Jacobian\n",
    "def jacobian(x):\n",
    "    (x1,x2,w1,w2)=(x[0],x[1],x[2],x[3])\n",
    "    x1s=x1*x1;x2s=x2*x2\n",
    "    return np.array([[0,0,1,1],\n",
    "                     [w1,w2,x1,x2],\n",
    "                     [2*w1*x1,2*w2*x2,x1s,x2s],\n",
    "                     [3*w1*x1s,3*w2*x2s,x1s*x1,x2s*x2]])\n",
    "\n",
    "# Initial condition\n",
    "x0=np.array([-1.,1.,1.,1.])\n",
    "\n",
    "# Newton's method\n",
    "if True:\n",
    "    xk=x0\n",
    "    fk=f(xk)\n",
    "    err=np.linalg.norm(fk)\n",
    "    while err>1e-16:\n",
    "        xk-=np.linalg.solve(jacobian(xk),fk)\n",
    "        fk=f(xk)\n",
    "        err=np.linalg.norm(fk)\n",
    "\n",
    "# Print solution\n",
    "print(\"\\nSolution is\",xk)"
   ]
  }
 ],
 "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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}