Assignment-1.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problem 1 - Taylor Approximation\n",
    "Consider the univariate function\n",
    "\n",
    "$f(x) = x^2 + \\log(x)$\n",
    "\n",
    "Suppose we want to approximate the function values in the vicinity of x = 1 using Taylor polynomials. Lets call the linear approximation function (linear Taylor polynomial) as L(x) and the quadratic approximation function (quadratic Taylor polynomial) as Q(x):<br>\n",
    "1. Plot f(x), L(x) and Q(x) in the interval [0, 2] in the same graph.\n",
    "2. Let the error associated with the linear approximation be given by eL(x) = f(x) − L(x), and the error associated with the quadratic approximation be given by eQ(x) = f(x)−Q(x). Plot eL(x)/x−1 and eQ(x)/((x−1)^2) in the interval [0, 2] in the same graph.\n",
    "## Part 1 Solution\n",
    "$L(x) = f(x_0) + f'(x_0)(x-x_0)$ <br>\n",
    "\n",
    "Since, $f'(x) = 2x + \\frac{1}{x}$ <br>\n",
    "we have, $f'(1) = 2 + 1 = 3$ <br>\n",
    "and      $f(1) = 1^2 + \\log(1) = 1 + 0 = 1$. <br>\n",
    "Therefore, \n",
    "\n",
    "$L(x) = 1 + 3(x-1) = 3x - 2$ <br>\n",
    "\n",
    "$Q(x) = f(x_0) + f'(x_0)(x-x_0) + \\frac{f''(x_0)((x-x_0)^2)}{2}$ <br>\n",
    "\n",
    "Since $f''(x) = 2 - \\frac{1}{x^2}$ <br>\n",
    "we have $f''(1) = 2 - 1 = 1$ <br>\n",
    "Therefore, \n",
    "\n",
    "$Q(x) = 1 + 3(x-1) + \\frac{(x-1)^2}{2} = 3x - 2 + \\frac{(x^2 - 2x + 1)}{2} = \\frac{x^2}{2} + 2x - \\frac{3}{2}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAyjklEQVR4nO3deVxU9f7H8dcXhh1EWRRFFLdccEXctdTU1NxKWyxNs8wsy8rqtt2617K6ldf2RdOWn+2L+76muSsuKG6IoqAiguzLMDPf3x9DXjRRlIGZgc/z8eDhzFk/58z45vA9y1dprRFCCOG8XOxdgBBCiLKRIBdCCCcnQS6EEE5OglwIIZycBLkQQjg5gz1WGhQUpMPDw+2xaiGEcFq7du06r7UOvny4XYI8PDycnTt32mPVQgjhtJRSCVcaLk0rQgjh5CTIhRDCyUmQCyGEk7NLG/mVFBYWkpiYSH5+vr1LKReenp7UrVsXNzc3e5cihKhkHCbIExMT8fPzIzw8HKWUvcuxKa01qampJCYm0qBBA3uXI4SoZBymaSU/P5/AwMBKF+IASikCAwMr7V8bQgj7cpggBypliP+lMm+bEMK+HCrIhRCiskrPNfKvhQfIzC+0+bIlyIUQopytik2m74wNzN2awPb4NJsvX4K8mA8//JDmzZtz//33M3/+fKZOnXrV6Z999lnWrl1bQdUJIZxNRm4hz/y8h/Hf7iTQx50Fk7rRp0Utm6/HYa5acQSffvopq1evpm7dunTt2pWFCxdedfonnniC8ePH07t37wqqUAjhLFbHJvPSvBhSc4w82bsxk3o3wd1QPsfODhnk/150gNjTmTZdZos61XhtcESJ4x999FHi4+MZMGAAo0aNwsPDg6CgIACGDh3K8OHDeeCBB/jiiy/YsGED3333HfXr1yc1NZWzZ88SEhJi03qFEM4pPdfI1EWx/L47iWYhfswZ24GWof4AnMo8RVi1MJuvU5pWinz++efUqVOHdevWUbNmTSIjIy+OmzlzJlOnTmXjxo1Mnz6djz766OK4yMhINm3aZI+ShRAOZuWBs/SdsYEFe0/zZO/GLJzUnZah/mQUZDBt6zQGzR/EhsQNNl+vQx6RX+3IuSKcOXOG4OD/PSmyVq1aTJ06lV69ejFv3jwCAgIujqtZsyanT5+2R5lCCAeRlmO9ImXh3tM0r12Nr4qOwi3awryj83g/+n3SC9K5p+k9tAluY/P1O2SQ25uXlxcZGRmXDIuJiSEwMPBvoZ2fn4+Xl1dFlieEcBBaa5bEnOG1BdbLCp/ucxMTezbC3eDCgfMHeHPbm+w7v492NdvxRacvaBbQrFzqkCC/gubNmzN37tyL77dv386yZcvYvXs3t9xyC/369bt4q/2RI0e466677FWqEMJOzmXm888F+1lxIJnWdf35bkQnmoVU40L+Bd7a8SG/HfmNAM8ApnWfxuCGg8v1pkAJ8iu4+eabmTJlClprjEYj48eP56uvvqJOnTpMnz6dcePGsXbtWkwmE3FxcURFRdm7ZCFEBdFa88uuRN5YHEu+ycILA5rxcPcGoCz8cOgHPt79MTmFOYxqMYqJbSbi5+5X7jVJkBdz4sSJi6/79OnDmjVr6NOnD3v37r04fMiQIQwZMgSAxYsXM2LECAwG2Y1CVAWn0nJ5aV4MG4+ep0N4Dd4e3ppGwb7sSt7FW9ve4vCFw3QM6cgLHV+gSY0mFVaXJFAJXnrpJbZt23bVaUwmE1OmTKmgioQQ9mK2aL7efIL3VhzGRcHrQyO4v1N9UvLO8fyGqSw7vowQnxCm3zKdvvX7VvizlSTIS1CrVq2LR94lkbZxISq/w2ez+Mdv+9hzKp2eTYOZdkcrgv1cmXNgNjP3zcRsMTOh9QTGtRyHt5u3XWqUIBdCiCsoMJn5ZN0xPlsfh5+nGx/c25bBrWuzIWkD49e+w6msU/QO681zHZ6jrl9du9ZqsyBXSrkCO4EkrfUgWy1XCCEq2s4Tafzjt30cS8nhjnahvHJ7czLMSTy29jE2JW2ioX9Dvuj7BV3rdLV3qYBtj8gnAweBajZcphBCVJiMvELeWX6I77adJLS6F9+M60i7cA8+3/sBPxz8AS+DF89FPcfI5iNxc3GcbhttEuRKqbrA7cA04BlbLFMIISqK1prl+8/y2sIDnM8u4OHuDXjy1kasPLWIwfM+4kL+Be5scidPtHuCQK/AG1vJ2f2w7TPo82/wCbJp/bY6In8feB4o8YJJpdQjwCMA9erVs9FqbcvX15fs7Oy/DX///fcJCAjggQceKHHee++9l9dff50mTSrukiMhRNklpefx2oL9rD54jog61Zg9pgNGtzgeXj2ag2kHiawZyWd9PqNFYIvrX7jWELcGtnwE8evBzRtaDIMmfW26DWUOcqXUIOCc1nqXUqpnSdNprWcCMwGioqJ0WddbUUwmE3PmzCE6Ovqq002cOJF33nmHWbNmVVBlQoiyMJktfLMlgekrD6M1vDywObe1cefDPa+z4sQKQnxCePfmd7kt/Lbrv5zQVAD7foYtn0DKQfANgVtfg6gHwauGzbfFFkfk3YAhSqmBgCdQTSk1V2s96oaXuOwFOBtjg9KKCWkFA96+7tnWrl1LZGQkBoMBk8lEly5dePfdd+nZsycvvvgiLi4uTJs2jR49ejB27FhMJpPcICSEg4tJzODFefvYn5RJr6bBvHR7I5Ynfc+di75BoXiszWOMbTkWL8N1PkcpNw12zIbtMyHnHNRqCcM+g5YjwOBePhuDDYJca/0i8CJA0RH5s2UKcQezadMm2rdvD4DBYODrr79mxIgRfPTRRyxfvvziTUMuLi40btyYvXv3XpxeCOFYsgtMTF95mG82nyDQ14MPR7YBn2geXf8i5/LOMaDBAJ5p/wwhPtfZv0DqMevR957vwZQHjftAl0nQsCdUwM1BjnnoeANHzuXlzJkzNG/e/OL7iIgIRo8ezaBBg9iyZQvu7v/7LfvXI20lyIVwLH+dzPz3oliSs/IZ1ak+A6Ly+Xjv88ScjyEiMILpPafTtmbb61koJGyGLR/D4WXg6gat74Euj0PN5tee34ZsGuRa6/XAelsu0968vLzIz8+/ZFhMTAzVq1fn3LlzlwyXR9oK4XhOpeXy6oL9rDucQova1Xh9RCirzs5hwpplBHsF80a3NxjcaDAuqpT97JhNEDvfegR+Ohq8AuDm56DDw+Bn+/44S8Mxj8gdSPPmzYmLi7v4/vfffyctLY0NGzYwaNAgtm/fTvXq1QHrI21btmxpp0qFEMUZTRZmbYznwzVHcXVRPD8gnELfNby4/VuA67+tPj8Tor+FbZ9DxikIbAy3/xfajAR3+9ya/xcJ8mJyc3OpW/d/t9o+88wzDB8+nNGjRwNw/vx5XnjhBdasWUNYWBiTJk1i8uTJfPPNNyQnJ+Pl5SV9dwrhALYcS+WV+TEcS8mhf0RNurQ5wTeHniDlRAoDGgzg6cinqe1bu3QLSz9lDe/ob6EgE+p3gwHvwE39wcUxesuUIC/GYrFccXhgYCBHjx6lSZMmHDly5OLwJ5988uLr77//ngkTJpR7jUKIkp3LyufNJQeZv+c0YQFevHSnG6uT3+W93QdpHdyaGb1mlL6rtaRoa/v3gfnW9xHDrCcwQyOvNpddSJCXwttvv82ZM2euerNP9erVLx65CyEqlslsYe7WBKavPEKBycKYW3xIc/+djw6uo7ZPbd7u8TYDGwy89vXgFgscWW4N8IRN4O4HnSdCpwlQ3TFvZAQJ8lJp2rQpTZs2veo0Dz74YAVVI4QobldCGv+cf4DYM5l0aeJJg8abWZjwO+6u7kyOnMyo5qPwNHhefSHGXNj7A2z9FFLjoFpd6PcGRD4Anv4VsyFlIEEuhHBK57ML+M+yQ/yyK5EQfwMj+xxnQ8oPxCZkM7zJcB5r+xhBXtd4pkn2Odg+C3Z8CXlpUKcdDJ9tvY3e1Xni0XkqFUIIrM0o3207yfSVh8k1mhjQ6RzHzT+zOCmRbnW6MSVqyrW7WTt3yNp8su9nMBuh6UDr9d/1u1bIDTy2JkEuhHAaO06k8doCazNKuyYZGIIW82d6DI2rN+bzPp/TLbRbyTNrDcf/gM0fQ9wqMHhBu1HQ+TEIalxxG1EOJMiFEA7vXGY+by07xLzdSYQEZNOty0b2pf9BUEEQr3V5jTsa34Gri+uVZzYZ4cDv1gBPjgGfmtDrFYgaBz43+EhaByNBXkxiYiKPP/44sbGxmM1mBg4cyPTp0/Hw8GD37t18/PHHzJ49u8T5P/74Y7y9vRk3blwFVi1E5WU0Wfhq03E+XHOUQp1Np/a7OJK/nKPZbkxsM5GxEWNLvqEn7wLs/Mr6AKusMxDcHIZ8DK3uArdrnPx0Mo5xNbsD0Fpz5513MmzYMI4ePcrRo0fJy8vj+eefB+DNN9+85LrxKxk3bhwfffRRRZQrRKX3x5EU+n+wgbeWx1C/4Q4Cmv2XQ3lLGdJoCIvvWMxjbR+7coinHYelz8N/I2DNvyG4Kdz/Gzy2BSJHV7oQBwc9Iv/P9v9wKO2QTZfZLKAZ/+j4jxLHr127Fk9Pz4uXEbq6ujJjxgzq16/PG2+8wb59+2jTxnojweTJkwkMDOTVV19lxYoVTJs2jfXr1+Pt7U14eDjbt2+nY8eONq1fiKrixPkc3lgSy+qDZ6kdeoj6rVdwyphMt1rdeKb9M9xU46Yrz3hqO2z+EA4tAeUKrUZYT2CGtKrYDbADhwxyezhw4MDfnlpYrVo1wsPD+fLLLy95hspbb71Fhw4d6NGjB08++SRLly7FpehW3aioKDZu3ChBLsR1yi4w8cm6OGZvPI6bbxwN26wmxRhPc9/mvN3+dbrU6fL3mSxmOLjIegVK4g7rNd/dJkPHCVCtlLfgVwIOGeRXO3K2h7y8PIKDgy++9/b2ZtasWdx8883MmDGDRo0aXRxXs2ZNDh2y7V8TQlRmFotm3u4k/rP8EOeNJ6jXbB2pln0Y3GrzZsc3ub3h7X9/MmFBNuyea72BJz0BaoTDgHeh7X3g4WuX7bAnhwxye2jRogW//vrrJcMyMzM5e/YsDRo0uOQZK2B9lG1gYCCnT5++ZLg8ylaI0os+eYF/L4pl35kThISvx9dtG0aDL8+2fpZ7m92Lh6vHpTNknoZtX8CuryA/A8I6We/AbHY7lHTVShUgJzuL3HrrreTm5vLtt9ZHXJrNZqZMmcKkSZNo167dJY+yTUhIYPr06ezevZtly5Zd7CUI5FG2QpTGmYw8nvpxN3d+vpoE/RP+N03H6BHNmIgxLLtzGWMixlwa4mf2we8T4P1W1nbwhj3hodXw0EpoMaRKhzjIEflFSinmzZvH448/zuuvv05KSgr33HMPL7/8MgAZGRlkZWXh6+vLQw89xHvvvUedOnWYPXs2Y8eOZceOHXh6erJp0yb+9a9/2XdjhHBQeUYzX2w4xhcbDoP/nwQ2+4NCncvghoOZ1HbSpY+WtVggbrW1B/rjG8DNBzqMh86PWptSxEUS5MWEhYWxcOFCADZv3szIkSOJjo4mMjKScePG8dNPP/Hwww+zevXqi/O0b9+emBhrR9G7d+8mIiKCwMDKcZOBELZisWgW7j3N28sPcJ7N+Ddah5E0OtXpweTIyTQNKPZQusJ82PeTtf075RD41YE+/4b2Y8qlB/rKQIK8BF27diUhIeHi+4kTJ/LLL79cdZ7z58/z+uuvl3dpQjiVnSfSmLokltgLm6lWZxVermdpGtSKp9u/R4eQDv+bMOe8tQf6HbMgJ8V62eAdMyHijnLtgb4ycKgg11pf+3nBduLp6XnN54337du3xHFaa1uXJIRDO5WWy9vLD7H86CZ8a6/EK+wEdarV56nIGdxa79b//V8/f9R6+eDeH8GUD036WTtwaHCzUz7Ayh4cJsg9PT1JTU0lMDDQYcP8RmmtSU1NxdOz8t1RJsTlMvIK+WRdHN/s2IwheBne4YcJ9ArmsbavMazxMAwuhqIHWG20BviR5eDqAW3utd7AE3z1Z/+Lv3OYIK9bty6JiYmkpKTYu5Ry4enpeUl/oEJUNoVmC99vO8mM9ZvJ912Ge/gefN38GN/6aUY2G4mXwQvMhbDvF2uAn9kD3oFwyz+sJzF9g6+5DnFlDhPkbm5uNGjQwN5lCCGuk9aalbHJTFuxjWSXxbjX2YGPq4EHIh5mbMRY/D38rdd8b5tl7cQ4MwkCm8Cg961H4W5y30VZOUyQCyGcz+6TF3hj6S5icubjEbgZTxcLI5oM59E2jxLsHQzpJ2Hdf6w90BuzILwHDJoBjfs6TA/0lYEEuRDiuiWk5vDW8r2sPf0bHkEb8PAqYGD4ACa1m0RYtTBI2gVLnofYBdYTlhF3Wtu/67S1d+mVkgS5EKLUUrMLeH9NLL8c+QVDwFo8aubQI/QWJkc+QdPqja0nLn99BE5uBo9q1vDuNAH85fxQeZIgF0JcU67RxOw/j/FF9E/o6itxq5lB26Aonu34FG38G8Oe72HrvZAWD/714La3rM/+9vCzd+lVggS5EKJEJrOFn3ee5L1NP1HguxyX4PM09Y/guU7v0dk33Nr7zs451t54QtvDiK+g+RCn6oG+MpC9LYT4G601y/efYdr6X0hzX4RrYDL1fBryXMd/0su9FmrrJxDzi/Vywma3W2/gqddZbuCxEwlyIcQlNsed51+rfyeRebhWSyTEoy7PdfwPt5ndcPnjIzi21toDfeQD1h7oAxtde6GiXEmQCyEA2J+UwasrFnCw4GcM3gkEutXi6XavMiQvH8PyN+DcAfCtBb3/ae2B3jvA3iWLIhLkQlRxx8/n8NqKJezI+A6DTzz+fgE80epp7k5PwW3JK5B9FmpGwNBPrf1gGjyuvVBRoSTIhaiizmbk86+Vy9hwbi6uvkfw9avOhCYPMvr8KTwWvwyFudCoNwz71PqvtH87LAlyIaqYtBwjb65exbLEb3DxPYiXny/j6g7jweQ4vJdPBRcDtL7beg14rQh7lytKQYJciCoiK7+Q/6xdx/wTX6F89uPu68OYoJ6MPx2Dz4YPrZ029JgCHceDX4i9yxXXQYJciEou12jiv+v/4Ke42eATg5uPF6P82vHIyd34Hf8WAhrCwPesPdC7+9i7XHEDJMiFqKTyC818tHEjcw9/icV7LwZvD0a6N+LRhGiqHT8M9bpY78BsOqDKd17s7Moc5EqpMOBboBaggZla6w/KulwhxI0xmix8unkj38R+iclrDwYvd+6nFhMS9+Bvjrf2Ot/lCajb3t6lChuxxRG5CZiitY5WSvkBu5RSq7TWsTZYthCilArNFj7f/CdzYmdh8tyDwdONBwp9GX/6ENUNPhD1iLUH+ur17F2qsLEyB7nW+gxwpuh1llLqIBAKSJALUQFMZgszt27iy/1fUOi5Fzd3A2NzFQ+lHKe6b224tagHek9/e5cqyolN28iVUuFAO2DbFcY9AjwCUK+eHBEIUVZmi+aLrRuYvX8mRo99uLu78nBWAWPTEvGv1QqG/RMihoGrm71LFeXMZkGulPIFfgOe0lpnXj5eaz0TmAkQFRUlXcoLcYOsAb6e2ftnYfSIwdPNhQnp2YzOuIB/49vg9kkQ3l1u4KlCbBLkSik3rCH+ndb6d1ssUwhxKbNF89mWNXx1YBZGj1i8DYrxF9K5L6eQaq3vtT7AKvgme5cp7MAWV60oYDZwUGv937KXJIQozmzRfLhpOXMPzcLofhQ/A0xMS+cekzt+HSZBh4fAJ8jeZQo7ssUReTdgNBCjlNpTNOwlrfVSGyxbiCqr0GRmxqYl/HzkSwrcj1PD1cITqRnc5V4bn5vfgNb3gJunvcsUDsAWV638CUhjnBA2YjSZeWfDPBYem0We+2mCXcw8lJrBnQFt8Ro0GRr3kR7oxSXkzk4hHESu0cibf/zCqhNfkut+njrKxLPnsxka3g+PgZOhdmt7lygclAS5EHaWlZ/P1HX/x8akr8lxyyScQh5MMzIo4n7cO0+EanXsXaJwcBLkQthJam4O/149k+0p35NjyOcmi5GxmW4MiHoCQ+QY8PC1d4nCSUiQC1HBTmdcYOrK99mduZBcg4nWpgJGFwZxW/fnUC2GyAOsxHWTIBeighxNOcvbq94kJv8P8lwtdC7MZ6RrU3r1fwVVv7O9yxNOTIJciHIWnRjPB2teY79lD0YX6JlfwMjA7nQd/Ir1WeBClJEEuRDlZO2RaGZv/BcHXOMB6Jtn4r76d9DuruekB3phUxLkQtjYrzuX8vPu9zjkdg4PF83gXAP3t3yYZl0ngMHd3uWJSkiCXAgbsFgszFn/FYvjZnHMI4dqrmbuLqjOmM7PEdZ6mDzASpQrCXIhyiDPWMCny95m9bn5JLqbqOlqZoypPmP7TCWofgd7lyeqCAlyIW5AavYFPl70En/k/EmKG4RjZqJLO8bc8TY+AWH2Lk9UMRLkQlyHY6eP8vnKl9hkOUiWqyLCrHmwel9G3jMVg6efvcsTVZQEuRClsC12Pd9ufoOtrmcxuig6FLpxR9j9DL71GXmAlbA7CXIhrmLhn1/zW+zn7HbPxmCAbgX+jGz/DF2jhtu7NCEukiAX4jLGQiPfLJvGirMLOOxhxs9gYaC5HuN6TuWmRnICUzgeCXIhilzISOGLJS+yLncrp90UIS4WRrpEMn7wWwQH1LV3eUKUSIJcVHlHEvYxZ80r/MkxMlxdaGxx4XG/AYy9/TU8PbztXZ4Q1yRBLqqs9TsX8FP0dLa5p2FygXYF3gyt/yDDek7AxVWeQCichwS5qFIsZjPfrfovy07+SIyHEQ83TbfCYO7v9CKdW/ezd3lC3BAJclElpGdfYNbiV1iftYGT7hDgamGYbs7D/aZRP7SpvcsTokwkyEWldujkfuaseZXN+jAZri40BMZ79eSh29/Ax8ff3uUJYRMS5KJSWr5zPr/snkG0WypmBZFGTwaF3sedfZ6S9m9R6UiQi0qjwFTI7JUzWJ34M0c9CvAyWOhlDOLejs/Rse0ge5cnRLmRIBdOLzH9HJ8t+yebczdz3gC1XcyMtDRjTL83CA1tYe/yhCh3EuTCaf1xeBtzN01jt4qnwEXR2mThAb/e3DdiKh4+NexdnhAVRoJcOJVCk4nZm+ay8siXHHXPwF1peuS7M/SmB+jV6wnpgV5USRLkwimczbrAjFXvsO3CclINJmq6mLgvvxb3dnuBBi3727s8IexKglw4tD+O72H2xv9wwLwfowu0MRXyoKE1IwZNw6dWE3uXJ4RDkCAXDqfQXMhnW39j+aGZnDKk4GGx0DfXQv/QQdxy90soL2n/FqI4CXLhME5mJPPu+s/ZeX4h2QYjdbSJh7M9uaPdY9TrOEZ6oBeiBBLkwq601iw7upVZ2z8l3rQHi4LOhXkMMNdnQK+X8WrcS3qgF+IaJMiFXWQX5PD+1h9ZcWwu6a7n8bVYuCc7j/4B3Ykc/jIEy/NPhCgtCXJRofYmH+a9zXM4kL6SQhcTN5mMPJphpl+z+wju8QT4Btu7RCGcjgS5KHeF5kLmxizh/2L+jxTLEdy0pl9uLoML/ejY5XncIu8DNy97lymE05IgF+XmZEYi/936f2w4vZBCl2xqF5qZnJXFbd7NCOszBZr0kx7ohbABCXJhU2aLmaXH1vF59FxO5u1CATfn5XFXdi4d6/XHa+ATUKedvcsUolKRIBc2kZyTzGe7vmfJ8fnkk4a/CR7JymBovqZuu7G4dH4U/EPtXaYQlZJNglwp1R/4AHAFvtRav22L5QrHZraY2XBqE59Gz+VQxjZQFtrlWhidlUo310C8u70I7UaBh5+9SxWiUitzkCulXIFPgL5AIrBDKbVQax1b1mULx5Sck8z/7f+FX4/+Ro75PF4mA/dn53J/Vhp1glvjOvB1aD4EXOUPPiEqgi3+p3UE4rTW8QBKqR+BoYAEeSVispj4M+lPZu/9kT2pmwFNg1xPHs26QJ/cbNya3Y4aNgnqdZYbeISoYLYI8lDgVLH3iUCnyydSSj0CPAJQr149G6xWVISk7CR+Ofw7Px/6jSxTKgaTJ7dlezA5K55Q3HBpNwo6T4TARvYuVYgqq8L+9tVazwRmAkRFRemKWq+4foXmQtaeWsv3sb8QnbIdNATlBjEpR3FXzhEMPrVQN78IUePAO8De5QpR5dkiyJOAsGLv6xYNE04m7kIcvx/9nXlHF5JtyoDCarTJDOWf+YdpakxA12yB6vMptBoBBg97lyuEKGKLIN8BNFFKNcAa4PcC99lguaICZBuzWX5iOb8dmcf+1H2gXXHLasjI3Go8XbALL0s+NOoNXSahGvWW9m8hHFCZg1xrbVJKTQJWYL38cI7W+kCZKxPlRmvNruRdzIubx4oTKykw54OxFsEXOvAP80n6FawFFwOq9d3Q5XGoFWHvkoUQV2GTNnKt9VJgqS2WJcrP2ZyzLDq2iPlx8zmZdRJXPDGmtyQq05t/GXbQsOA3tFcNVI8p0HE8+IXYu2QhRCnIhb6VXIG5gHWn1jH/6Hy2nNmCRVvwNDXBcm4Y9xrTeNxtNdXNZ8G/Idz6HqrtfeDuY++yhRDXQYK8EtJas//8fhYcW8DS40vJMmbh6xqES8ateJ6rzzOe0dyh5+CusyGkC3SZDk0HSA/0QjgpCfJKJDknmcXxi1l4bCHxGfG4u7gT7NqBtMQW1MiGV2qspbPhG5RJo1oMhS6ToG57e5cthCgjCXInl2fKY+3JtSw6tuhi00kjv5aE6zHsP1if9q6H+aTa7zTy2A0mP+g0wfpTo769SxdC2IgEuROyaAvRydEsil/EihMryCnMIcSnNj2C7+bIsWYcPOjCA95b+Mb/WwLyEsA1FPq+Du3HgKe/vcsXQtiYBLkTOZFxgkXxi1gSv4Sk7CS8Dd70rHsrngUdWbnLhy0ZKUyutoK7q63E03gBqreFga9Ci6Hg6mbv8oUQ5USC3MGl5aex/PhyFscvJuZ8DC7KhU4hnRjT7FHiT4bz47pkahYk8Hr1tfTyXoOr0Qg3DYCuk6B+N7mBR4gqQILcAeWZ8vjj1B8sjl/MpqRNmLSJm2rcxJT2U2hToxe/bs/ktR9OEWVZyffVV9OKbVDoCe3ut97AE9TE3psghKhAEuQOwmwxs/3sdpbEL2H1ydXkFOZQ06smo1uMZlCjQVgKQvhs/THe2ruHIa5bWe+7ijr5R0EFQ8+XoMND4BNk780QQtiBBLkdaa2JTYtlSfwSlh9fTkpeCr5uvvSt35dBDQcRVSuKmKQs3lkYx/aDexjrvo4dvquoVpgCfs2g30fQ6m5w87T3pggh7EiC3A4SMhNYGr+UpceXciLzBG4ubvQI7cHtDW/n5ro342nwZOeJNMZ+tZPjcbFM9FjBp97rcbfkQd1boOsT0OhW6YFeCAFIkFeY5Jxklp9YztLjS4lNjUWh6BDSgbERY+lTvw/+HtbLArcfT+ODNXvIPbaVxz2W0dtjO8rFBdVyuLX9u3YbO2+JEMLRSJCXowv5F1iVsIqlx5cSnRyNRtMisAXPRj1L//D+1PKpdXHaXQkXmLEyFp/jK3nefRltPA6jPaqhop6AjhOkB3ohRIkkyG0sy5jF2pNrWXZiGVtPb8WszTTwb8DEthMZED6AcP/wS6aPSczg4xV7qBX/G2+5rSDM/SyW6vWhyzuotveDh699NkQI4TQkyG0gtzCXdafWseLECv5M+pNCSyGhvqGMiRjDwAYDuanGTajLrueOO5fN7GWbqHt0Lu8Y1uDvloM5tAN0eweXZoPkAVZCiFKTIL9BuYW5bEzayIoTK9iYuJF8cz41vWtyT9N7GNBgAK2CWv0tvAHOZuTzw6KlhB3+iqkum3A1aMxNb4fuT+Ia1tEOWyKEcHYS5Nchz5THxsSNrExYyYbEDeSZ8gj0DGRY42H0b9CfdjXb4aKufCVJVp6RFQu+p/bBL3laxWB088LU5kHcekzCJaBBBW+JEKIykSC/hr+OvFclrLoY3gGeAQxuOJj+DfoTWTMS16s0g5iN+exc9AVBMbMYwSnS3YJI7/gS1buPlx7ohRA2IUF+BbmFuWxI3MDKhJX8mfTnxfAe1HAQt4XfRvta7TG4XGPX5aZxcsVH+O6bQyedzglDA052/S/1bh4NBveK2RAhRJUgQV4k05jJH6f+YFXCKjaf3kyBuYBAz0CGNBpCv/r9aF+r/VWPvC86H0fOhg9wi/mJerqALS7tUN0m0an3nSi5gUcIUQ6qdJCn5aex7uQ6Vp9czdYzWzFZTNT0rsnwJsPpF96PtsFtSxfeWkPCZiybP0IdWY6bdmWB7kF++0e5a2BfPN3kChQhRPmpckF+Nucsa06uYc3JNexK3oVFWwj1DWVU81H0qd+HVkGtSjxh+TdmE8TOhy0fw+ndZCk/vjYNI6HhfTw9rDthAd7lui1CCAFVJMjj0+MvhveB1AMANK7emIdbPUzf+n1pWqPpFS8VLFF+JkR/A9u+gIxTnPcIY0bhODb79OWlu9ozuUWtay9DCCFspFIGuUVb2H9+P2tPrmXNyTWcyDwBQKugVjwV+RS31rv1b3dYlkr6Kdj2Oez6BoxZZNTswJtuo/g1M4IHujZkUb+m+HpUyl0qhHBglSZ1jGYj285sY92pdaw/tZ6UvBQMykD7kPbc1/w+eoX1IsQn5MYWnhRtbT45MB+AwubD+MI4gPf2e9Mw2IdfRrchsl4Nm22LEEJcD6cO8oyCDDYkbmD9qfVsOr2JnMIcvAxedA/tTu96vekR2uPiUwWvm8UCR5ZbAzxhE3hUg84T2R82kseXpHAqLZcJNzfk6b43yclMIYRdOV2Qn8o8ZT3qTlxPdHI0Zm0myCuI/uH96V2vN51qd8LD1ePGV2DMhb0/wNZPITUO/MOg3zTM7Ubz2ZZzzJh7lJBqnvz4SBc6NpAbeoQQ9udUQf7G1jf46fBPgPVk5YMtH6RXWC9aBrUs/ZUmJck+B9tnwY4vIS8N6kTCiDnQfCgpuWYmz93N5mOpDGlTh9eHtcTfS3qlF0I4BqcK8p5hPWng34Bb6t5CXb+6tlnouYOw5RPY9zOYjdB0oLUH+npdQCl2nEhj0vfRpOcW8u6I1oxoX/f6rnARQohy5lRB3j20O91Du5d9QVpD/DprgMetBoMXtBtl7YEnsFHRJJqvNx3njSUHCavhxdePd6R57WplX7cQQtiYUwV5mZmMsP9Xa4An7wefmtDrFYgaBz6BFycrNFt4dcEBfth+kr4tajH97jZU85SmFCGEY6oaQZ6bBru+gm0zIfssBDeHoZ9Aq7vAcOmJ0fRcIxPnRrMlPpXHejbi2X5NcXGRphQhhOOq3EGeFg9bP4Pdc6EwFxr2sgZ441vhCu3cp9PzGD17G6fS8vjv3W24M9JG7fBCCFGOKmeQn9wGWz6Cg4vBxWA98u7yOIS0LHGWuHPZPDB7G1n5Jv7voY50ahhY4rRCCOFIKk+QW8xwcJH1Bp7EHeBZHXo8Ax3GQ7XaV511X2I6Y+Zsx9XFhR8ndCaizg3eRCSEEHbg/EFekG1tOtn6KaQnQI0GMPA9aHsfuPtcc/Z9iencP2sb/t5uzH2oE+FB155HCCEcifMGeeZp69MHd30F+RkQ1glum2a9DryUPdDHns5k9Ozt+Hu78fOELtSp7lXORQshhO2VKciVUu8CgwEjcAx4UGudboO6SnZmn7X5ZP9voC3QfAh0fQLqRl3XYo4mZzF69ja83V35YXxnCXEhhNMq6xH5KuBFrbVJKfUf4EXgH2UvqwTLXoBtn4Gbj7Xtu/OjUCP8uhdzJiOPUbO34eKi+O7hTtIBhBDCqZUpyLXWK4u93QqMKFs519CoN/iFQPux4FX9hhaRU2Dioa93klNg5teJXWgY7GvTEoUQoqLZso18HPBTSSOVUo8AjwDUq1fvxtZwUz/rzw0yWzSTf9zDobOZzBnbgWYhcsu9EML5XTPIlVKrgSv1yPCy1npB0TQvAybgu5KWo7WeCcwEiIqK0jdUbRm9vewgqw8m8+8hEfRsWtMeJQghhM1dM8i11n2uNl4pNRYYBNyqtbZLQJfG8v1nmLXxOKM712dM13B7lyOEEDZT1qtW+gPPA7dorXNtU5LtnU7P4x+/xdC6rj//HNTC3uUIIYRNlbE3Bj4G/IBVSqk9SqnPbVCTTZktmqd+2oPJbOHDe9vhbijrJgshhGMp61UrjW1VSHn5ZF0c24+nMf2uNnLXphCiUqrUh6f7kzL4YM1RhrWtw52RofYuRwghykWlDXKLRfPqgv3U8Hbj30NbSvdsQohKq9IG+e+7k4g+mc4/+jeTjpKFEJVapQzyjLxC3l52kHb1qjNcOocQQlRyzvv0w6t4f/URUnOMfP1gR+mmTQhR6VW6I/K4c1l8uyWB+zrWo2WodBAhhKj8Kl2Qf7r+GO6uLkzp19TepQghRIWoVEGelJ7Hwj2nubdjGAE+7vYuRwghKkSlCvLZG4+jgYd7NLR3KUIIUWEqTZCn5xr5ccdJhrSpQ6j09iOEqEIqTZB/uyWBXKOZCbfI0bgQomqpFEGeZzTz9eYT9GoaLJ1FCCGqnEoR5Iv2niYtx8iEWxrZuxQhhKhwlSLIF8ecoV6AN50aBNi7FCGEqHBOH+TpuUY2x51nYKva8mAsIUSV5PRBvvJAMiaL5vZWte1dihBC2IXTB/mSmDOEBXjRMlROcgohqianDvL0XCObpFlFCFHFOXWQr4y1NqsMbCnNKkKIqsupg3xZzBnq1vCidV15yqEQoupy2iDPyCvkT2lWEUII5w3y1bHJFJo1A+VqFSFEFee0Qb4z4QL+Xm60kWYVIUQV57RBfuB0Bi1Dq0mzihCiynPKIC80Wzh0JouWdeRoXAghnDLIjyZnYzRbiJA+OYUQwjmD/MDpDAAi6sjdnEII4aRBnomPuysNAn3sXYoQQtidUwb5/qQMWtSphouLnOgUQginC3KzRRN7JpMIOdEphBCAEwb58fM55BrNtJQTnUIIAThhkMuJTiGEuJQTBnkm7gYXGtf0tXcpQgjhEJwuyPcnZdA8xA83V6crXQghyoVTpaHWmv1JGXIjkBBCFONUQZ54IY/MfJO0jwshRDFOFeT7k6wnOuUZK0II8T82CXKl1BSllFZKBdlieSU5cDoTVxdF0xC/8lyNEEI4lTIHuVIqDOgHnCx7OVcXFuDF8MhQPN1cy3tVQgjhNAw2WMYM4HlggQ2WdVX3dKjHPR3qlfdqhBDCqZTpiFwpNRRI0lrvLcW0jyildiqldqakpJRltUIIIYq55hG5Umo1EHKFUS8DL2FtVrkmrfVMYCZAVFSUvo4ahRBCXMU1g1xr3edKw5VSrYAGwN6i7tbqAtFKqY5a67M2rVIIIUSJbriNXGsdA9T8671S6gQQpbU+b4O6hBBClJJTXUcuhBDi72xx1QoAWutwWy1LCCFE6ckRuRBCODkJciGEcHJK64q/ElAplQIkXMcsQYAjnkR11LrAcWtz1LrAcWtz1LrAcWtz1LqgbLXV11oHXz7QLkF+vZRSO7XWUfau43KOWhc4bm2OWhc4bm2OWhc4bm2OWheUT23StCKEEE5OglwIIZycswT5THsXUAJHrQsctzZHrQsctzZHrQsctzZHrQvKoTanaCMXQghRMmc5IhdCCFECCXIhhHBydg9ypVR/pdRhpVScUuqFK4z3UEr9VDR+m1IqvNi4F4uGH1ZK3VbBdT2jlIpVSu1TSq1RStUvNs6slNpT9LOwgusaq5RKKbb+h4uNG6OUOlr0M8aWdZWythnF6jqilEovNq4899kcpdQ5pdT+EsYrpdSHRXXvU0pFFhtXbvusFHXdX1RPjFJqs1KqTbFxJ4qG71FK7bRlXaWsradSKqPYZ/ZqsXFX/R6Uc13PFatpf9H3KqBoXLntM6VUmFJqXVEmHFBKTb7CNOX3PdNa2+0HcAWOAQ0Bd2Av0OKyaR4DPi96fS/wU9HrFkXTe2B9nO4xwLUC6+oFeBe9nvhXXUXvs+24v8YCH19h3gAgvujfGkWva1RkbZdN/wQwp7z3WdGybwYigf0ljB8ILAMU0BnYVkH77Fp1df1rfcCAv+oqen8CCLLjPusJLC7r98DWdV027WBgbUXsM6A2EFn02g84coX/m+X2PbP3EXlHIE5rHa+1NgI/AkMvm2Yo8E3R61+BW5VSqmj4j1rrAq31cSCuaHkVUpfWep3WOrfo7Vasz2Mvb6XZXyW5DViltU7TWl8AVgH97VjbSOAHG66/RFrrDUDaVSYZCnyrrbYC1ZVStSnnfXaturTWm4vWCxX3Hftr3dfaZyUpy3fU1nVV5HfsjNY6uuh1FnAQCL1ssnL7ntk7yEOBU8XeJ/L3jb84jdbaBGQAgaWctzzrKu4hrL9p/+KprN3abVVKDbNRTddT1/CiP91+VdbOsa9n3vKujaJmqAbA2mKDy2uflUZJtZf3Prsel3/HNLBSKbVLKfWInWrqopTaq5RappSKKBrmEPtMKeWNNQx/Kza4QvaZsjb/tgO2XTaq3L5nNnuMbVWllBoFRAG3FBtcX2udpJRqCKxVSsVorY9VUEmLgB+01gVKqQlY/5rpXUHrLq17gV+11uZiw+y5zxyaUqoX1iDvXmxw96L9VRNYpZQ6VHS0WlGisX5m2UqpgcB8oEkFrv9aBgObtNbFj97LfZ8ppXyx/vJ4SmudactlX429j8iTgLBi7+sWDbviNEopA+APpJZy3vKsC6VUH6x9lw7RWhf8NVxrnVT0bzywHutv5wqpS2udWqyWL4H2pZ23vGsr5l4u+5O3HPdZaZRUe3nvs2tSSrXG+jkO1Vqn/jW82P46B8zDds2KpaK1ztRaZxe9Xgq4KaWCcIB9VuRq37Fy2WdKKTesIf6d1vr3K0xSft+z8mj4v44TBAasDfsN+N+JkYjLpnmcS092/lz0OoJLT3bGY7uTnaWpqx3WkzpNLhteA/Aoeh0EHMVGJ3tKWVftYq/vALbq/51QOV5UX42i1wEV+VkWTdcM60knVRH7rNg6win5xN3tXHoSantF7LNS1FUP67mfrpcN9wH8ir3eDPS3ZV2lqC3kr88QayCeLNp/pfoelFddReP9sbaj+1TUPiva9m+B968yTbl9z2z6wd/gDhiI9QzvMeDlomFTsR7lAngCvxR9obcDDYvN+3LRfIeBARVc12ogGdhT9LOwaHhXIKboCxwDPFTBdb0FHCha/zqgWbF5xxXtxzjgwYr+LIve/wt4+7L5ynuf/QCcAQqxtj8+BDwKPFo0XgGfFNUdg7Xv2XLfZ6Wo60vgQrHv2M6i4Q2L9tXeos/65XL4LK9V26Ri37OtFPtlc6XvQUXVVTTNWKwXQhSfr1z3GdZmLw3sK/Z5Dayo75ncoi+EEE7O3m3kQgghykiCXAghnJwEuRBCODkJciGEcHIS5EII4eQkyIUQwslJkAshhJP7f/xYbrJw6BwAAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def f(x):\n",
    "    return x**2 + np.log(x)\n",
    "\n",
    "def L(x):\n",
    "    return 3*x - 2\n",
    "\n",
    "def Q(x):\n",
    "    return x**2/2 + 2*x - 3/2\n",
    "\n",
    "x = np.linspace(0, 2, 100)\n",
    "while 0 in x:\n",
    "    x = np.delete(x, np.where(x == 0)) # remove 0 from x because log(0) is undefined\n",
    "\n",
    "plt.plot(x, f(x), label='f(x)')\n",
    "plt.plot(x, L(x), label='L(x)')\n",
    "plt.plot(x, Q(x), label='Q(x)')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 2 Solution\n",
    "$eL(x) = f(x) - L(x)$\n",
    "\n",
    "$eQ(x) = f(x) - Q(x)$\n",
    "\n",
    "To plot:\n",
    "\n",
    "$\\frac{eL(x)}{(x-1)}$ and $\\frac{eQ(x)}{(x-1)^2}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAArH0lEQVR4nO3deXxU9b3/8dc3k2Wy7yHLZIGwgwFCFBGvdUHF0sr1utWrUgpqvXaztv3d3vq7dvnZa29vbe+13t7WKg9rte5iq1YUt4tVQQMEQlhC2LJDEsgkIfvM9/fHmYQhhDAhM3NyJp/n4zGPMznnzDkfToY3X77ne85RWmuEEEJYV5jZBQghhBgbCXIhhLA4CXIhhLA4CXIhhLA4CXIhhLC4cDN2mpaWpgsKCszYtRBCWNaWLVuatdbpQ+ebEuQFBQWUlpaasWshhLAspdTh4eZL14oQQlicBLkQQlicBLkQQlicKX3kQghz9fX1UVtbS3d3t9mliGHY7XYcDgcRERE+rS9BLsQEVFtbS3x8PAUFBSilzC5HeNFa09LSQm1tLZMnT/bpM2PuWlFK5Sql3ldK7VJKVSilvjXWbQohAqu7u5vU1FQJ8XFIKUVqauqo/rfkjxZ5P/AdrfVWpVQ8sEUptUFrvcsP2xZCBIiE+Pg12t/NmFvkWusGrfVWz/t2YDeQM9btDue9PUf4zQdVgdi0EEJYll9HrSilCoAFwOZhlt2llCpVSpU2NTWd0/Y/3NfMo+9VIfdQF0KIk/wW5EqpOOBl4F6tddvQ5Vrrx7TWJVrrkvT0064w9UlOUjSdvS7auvrHWK0QYjwqKCigubn5tPmvvvoqP/nJT0b87He/+13ee++9U+Y999xz/PSnP/V5G0OtXr2ajIwM5s6de9Z9mckvQa6UisAI8We01q/4Y5vDyUqMBqCutStQuxBCjEM///nPueeee0Zc5xvf+AY/+9nPTpn35ptvsmzZMp+3MdSqVatYv369T/sy05hPdiqjV/4JYLfW+pdjL+nMspPsANS3djE7OyGQuxJiwvjxaxXsqj/tP9FjMjs7gR9+cc6I6zz99NM88sgj9Pb2smjRIn7zm98Mu15lZSVRUVGkpaUBsGLFCq6//npWrlzJ7373OzZu3MgzzzxDfn4+LS0tNDY2kpmZidaasrIyiouLfd7GUJdccgmHDh06bf7QfZnNH6NWlgC3A+VKqTLPvB9orf/qh22fIifJaJE3OKVFLoSV7d69m+eff56PPvqIiIgI7rnnnmGDFOCjjz6iuLh48OfHHnuMJUuWMHnyZB5++GE2bdo0uKy4uJiPPvqI66+/nm3btjFv3jyUUqPahq+892W2MQe51vpvQFDGMaXFRRFhU9S1ytVoQvjL2VrOgfDuu++yZcsWzj//fAC6urrIyMgYdt2Ghga8z6tNmjSJn/zkJ1x22WWsW7eOlJSUwWUZGRnU19cDsH79eq655ppRb8NX3vsym6XutRIWpshMtEuLXAiL01rz5S9/mbKyMsrKyti7dy8/+tGPhl03Ojr6tItjysvLSU1NPS1Iu7u7iY42/uf+9ttvc9VVV/m8jZqaGubPn8/8+fP57W9/e9Y/g/e+zGapIAfjhGe9nOwUwtKuuOIKXnrpJY4ePQrAsWPHOHx42FttM2vWLKqqTl4/8umnn/Lmm2+ybds2fvGLX3Dw4MHBZZWVlcydOxen00l/fz+pqak+byM3N3fwH5a77777rH+GgX2NB5YL8pykaOqla0UIS5s9ezYPPvggV111FUVFRVx55ZU0NDQAUFRUhMPhwOFwcN9993HJJZewbds2tNb09PRw5513snbtWrKzs3n44YdZvXo1Wmv6+vqoqqqipKSEDRs2sHTp0sH9+bqNoW655RYWL17M3r17cTgcPPHEEwCn7Gs8sNxNs7IS7TS2deNya2xhcomxEFZ18803c/PNN58yb7gRIgBLly7l3XffZenSpWzfvn1w/rXXXsu1114LwOuvv84NN9xAeHg469ev54477hhcLyYmxqdtDPXss88OO997X+OB5Vrk2UnRuNyao+3SKhdiovjBD35AZ2fniOv09/fzne98B4DHH3+cCy+8cNTb8JX3vsaD8fHPySgMDEGsb+0evEBICBHaJk2adMZW84Abb7xxzNvw1dn2FWyWa5FneV0UJIQQwoJBni0XBQkhxCksF+QJ9gjiosJl5IoQQnhYLsjBuOeKdK0IIYTBokEeTb10rQgRcpxOJytXrmTq1KkUFhZy6623cvz48cHlDQ0NfOELXxhxG6+//joPPPDAKfMaGhoGr/L0ZRtDPfroo0ydOhWl1Cm32R1uXwAffPABcXFxpwyBBCgrK2Px4sXMmTOHoqIinn/++VHVcSaWDHLj6k7pWhEi1KxZs4YpU6ZQVVXF/v37mTp1KqtWrRpc/stf/pI777xzxG0sX76c11577ZShhuvXr+fqq6/2eRtDLVmyhHfeeYf8/Pyz7mvnzp3cc889bNq0ifb2dn784x8PLouJieGpp56ioqKC9evXc++999La2jqqWoZjueGHADlJdo6d6KW7z4U9wmZ2OUJY25vfh8Zy/24z8zy4ZuT7dQ+9je19993Hli1bTmmlPvDAAxQWFrJ3715mzJjByy+/zIMPPgjAr371K8rLy1m7di3l5eXccsstfPrpp8TExHDppZfy+uuvc9NNNwFGkP/whz8E8Hkb3hYsWDDsn0Epdcq+6urqWLNmDa+++irTp0/nT3/6E7fffjtr165l9erVTJ8+ffCz2dnZZGRk0NTURFJS0uiO7xCWbZGDDEEUwqq8b2NbVlaGzWZj8+bNzJ8/H5vtZOPMZrOxYMECdu/ezcGDB0lOTiYqKgqAb33rW1RVVbFu3Tq+8pWv8Lvf/W4wgEtKSvjwww8BcLlc7N27l9mzZ49qG77y3ldOTg6bN28eDGybzcaf/vQnVq9efdrnPv30U3p7eyksLBzl0TudJVvkJ4cgdjMlPc7kaoSwuLO0nANhuNvYbtq06bSuC29Db0UbFhbGk08+SVFREV/96ldZsmTJ4DLvW8xu3ryZRYsWjXobvjqX29k2NDRw++2384c//IGwsLG3py0a5MZFQfLINyGsaeA2tg899NDgvKqqKq688krcbvdguLndbrZv305xcTEtLS2n3Yp23759xMXFjXg7W+/HvQ13O9vhtnH11Vdz5MgRSkpKePzxx0f8s4z2drZtbW0sX76cn/70p6fdRuBcWbJrJTNRru4UwsqGu41tREQECxYsGOy/BnjwwQe54ooryMvLY/r06afcVMvpdPLNb36TjRs30tLSwksvvTS4zPsWswM3ygJ83sZbb71FWVnZWUN86L7Opre3l+uuu46VK1dyww03+PQZX1gyyKPCbaTHR9EgI1eEsKQz3cZ27dq17Nu3j8LCQtLT09m0adPgQx5iY2MpLCwcvK/4t7/9bb72ta8xffp0nnjiCb7//e8P/sPw/vvvs3z5cpqamrDb7cTHx496G94eeeQRHA4HtbW1FBUVnTKscGBfvnjhhRfYuHEjTz755OBDLMrKys75OA7SWgf9tXDhQj1W1/76Q33b45vGvB0hJqJdu3aZXcJZ7dmzRxcWFuo33nhjcN4rr7yi77///hE/19jYqC+//HKttdZ//OMf9UMPPXTKcl+24SvvffnbcL8joFQPk6mW7CMH44Rn5ZF2s8sQQgTIjBkzTnmqD8B1111HS0vLiJ+rrq7m4YcfBuC22247bbkv2/CV977MZNkgz0qM5oO9TWitUUoeMCHEaFn1787QqyWHGhgJM5Zt+MqXfZ0LPczTikZiyT5ygMlpMXT1uWTkihDnwG6309LSMurAEIGntaalpQW73e7zZyzbIi9yJAGwo9aJI3l0A/iFmOgGTtw1NTWZXYoYht1ux+Fw+Ly+ZYN8ZlY8kbYwtte08vnzsswuRwhLiYiIYPLkyWaXIfzEsl0rUeE2ZmUnUFbTanYpQghhKssGOcB8RyLldU5cbunnE0JMXJYO8nm5SXT2utjf1GF2KUIIYRpLB/nACU/pXhFCTGSWDvIpabHER4WzXYJcCDGBWTrIw8IURbmJbK9tNbsUIYQwjaWDHGCeI4k9De1097nMLkUIIUxh/SDPTaLframobzO7FCGEMIXlg3x+bhKA9JMLISYsywf5pAQ7mQl2dkg/uRBigrJ8kAPMy01ke63T7DKEEMIUfglypdRapdRRpdROf2xvtIocSRxsPsHxE71m7F4IIUzlrxb5k8AyP21r1C4qTAXgg8rTH9EkhBChzi9BrrXeCBzzx7bOxTxHEhnxUbxdccSsEoQQwjRB6yNXSt2llCpVSpX6+x7IYWGKK2dP4oO9TTKeXAgx4QQtyLXWj2mtS7TWJenp6X7f/tVzMunqc/G3fc1+37YQQoxnITFqBeDCKanE28N5q6LR7FKEECKoQibII8PDuHxmBu/sPkK/y212OUIIETT+Gn74LPAJMEMpVauUWuOP7Y7WVbMzOd7ZR+nh42bsXgghTOGXZ3ZqrW/xx3bG6nMz0okMD+PtiiNcOCXV7HKEECIoQqZrBSAuKpyLp6bxVkUjWsvj34QQE0NIBTnA1XMmUdfaxQ65ZF8IMUGEXJAvm5tFdISNpzcdNrsUIYQIipAL8sToCK5fmMOft9fT0tFjdjlCCBFwIRfkAF9eXEBvv5vnPqsxuxQhhAi4kAzyaZPiuXhqGk9vOkyfjCkXQoS4kAxygFUXFdDg7JYbaQkhTKW1prqlk3Xbarl/XTk1xzr9vg+/jCMfjy6bmUFeSgxPfnyQ5UVZZpcjhJggevpd7KxrY+vh42w5fJwt1cdpajfO18VFhXPVnExyU2L8us+QDXJbmGLl4nwefGM35bVOznMkml2SECIEHW3rZmu1Edpbq1spr3XS6+nSzUuJYUlhKgsLUijJT2b6pHhsYcrvNYRskAPcWJLLr9+r4udv7eGPaxaZXY4QwuL6XW72NLYPBveWw8epPd4FGPd7KspJZNWSAorzkijOTyYj3h6UukI6yBOjI/jmFdP4f6/v4oO9R7l0RobZJQkhLOTYiV62Hj7O1mrjtb3GSZfnmQcZ8VEszE9m1UUFFOcnMyc7gahwmyl1hnSQA9x+YT5PfXKIf/vrbi6emka4LWTP7wohxqDf5WbvkXa2VreyzRPeh1qME5PhYYrZ2QncfH4uC/KSWJifTE5SNEr5v5vkXIR8kEeGh/HPy2ZyzzNbeWlLLV+6IM/skoQQ40BzRw9l1a2Dre0dtU46e43WdlpcJAvykvnSBXkU5yVzXk4i0ZHmtLZ9EfJBDnDN3EwW5ifz8IZKvjgvm9ioCfHHFkJ49Pa72d3Qxrbq42yraWVbdSvVx05tbd+40EFxfjLFeck4ksdPa9sXEyLRlFLcv3wW//Cbj/nPdyq5f/lss0sSQgSI1pp6Zzfbqo9TVt3KtppWyuuc9PYbI0kmJURRnJfMbRfmMT93/Le2fTEhghygOC+ZWxfl8fjfDnL5zEksLpT7lQsRCk709LOj1klZTStlNcfZVt3KUc+47ajwMM7LSWTlhfkU5yezIC+JrMRokyv2vwkT5AD3L5/FR1XNfPfF7bx579+RYI8wuyQhxCi43Jp9R9spq271BHcrlUfacXseP1CQGsOSqWnMz01ifm4Ss7ISiAwP/QEOEyrIYyLD+dXN87nht5/w47/s4uGb5pldkhDiDLTWNDi72V5zMrTL606ekEyMjmBebhJXzclkQV4S8x1JJMdGmly1OSZUkAMsyEvma5cW8sh7VVw2M50vFGWbXZIQAnB29VFe62R7rRHa22tOdpFE2BSzsxO5caGD+XlJzM9NpiA1xlInJANpwgU5wDeumMaHVc1878Ud5KfEyuX7QgRZd5+LXQ1t7KhpZXutk+01rRxoPjG4fEpaLBcVphpdJHnJzMqKN+1iGytQZjzbsqSkRJeWlgZ9v96Otndz3X9/TJ/LzZ+/viQkT4AIMR70u9xUHulgR60R2uV1rexpaKff07GdHh/F/Nwk5jkSKXIkMc+RRGKMnL8ajlJqi9a65LT5EzXIAfY0tnHD/3xCXkoML969WMaXCzFGbrfmQPMJyuta2V7jpLzOSUW9k+4+Y+hfgj2cIkcS5zkSmedIYl5uIpkJduki8dGZgnxCJ9fMzAR+/Y8LWPPkZ9z99BZ+v7IEe4T8900IX7jdmsPHOimvc1Je28qOWicV9W109PQDEB1hY25OArcuyqfI09qWfu3AmNBBDnDZjAx+dn0R//zyDu74Qym/X1li+YsDhPA3rTWHW4zQ3llntLTL65y0dxuhHRkexuysBP6hOIe5OUZruzA9Vu5tFCQTPsgBbirJJUwpvvfSdlY/+RlPrCohJlIOjZiY3G7NwZYT7KwzWtjltU521nuFti2MmVnxfHFeNkU5iZznSGT6pHgiJLRNI2nlccNCBxE2xbefL+PWxzfz2O0lpMdHmV2WEAHV73JT1dTBzro2KuqdVHimJzxjtSNtYczyhPZ5OYmcl2OE9kS4yMZKJMi9rJifQ1R4GPc+X8aKR//GYytLmJsjQxNFaOjqdbGnsY2KeuO1q97JnsZ2ejz3IImOsDE7O4HrFzqYm5PI3OxEpk2Kk5a2BUzoUStnsrPOyZ1PldLa2ccvbpwnz/wUlnPsRC+76tvY1eD0hHYb+5s6Bi9lT7CHMyc7kbk5CczJTmROdgJT0uMC8hgy4T8y/HCUjrZ3c/cft7C1upUvnZ/Lv35htgxPFOOOy6051HKC3Q1tnlc7u+rbaGzrHlwnM8HOnOwE5mQnMDvbCG6r3aY1aFz90N8Ffd3Q1wn93dDXNcK05+T6/Z6fvef393jW7Tam/T3wxf+C/MXnVJ4MPxyljHg7z921mF+9U8lv/3c/mw608Mub51Ocl2x2aWKCauvuY09DO3sajdDe1dBOZWP74KPHbGGKqelxLC5MZVZWPHOyE5mVlUBKKNx/RGtPSHYar95O6DthhGavZ15f18l5gz97vx+Ydp9839916jx337nXGG43XhHREB4F4QNTO0TGQEzKyfd+Ji1yH2w+0MJ9L2yn3tnFbYvy+e5VM+TKMxEwfS43B5tPsKexnb2NbextbGd3Qzt1rV2D6yRGRzArK55ZWQnMykpgdlYCUzPizL8OQmtw9UJPB/R2QO8Jz+ss7/s6T84bCOreDq/Q7gTtGl0tygaRsZ5gtUNEjPE+IvrU94PLPNOBMB5cNtLUK7CD8D8c6VoZo7buPn75diVPfXKI5JhI/vmamVxf7JA+RXHO3G5NXWsXlUfa2Xuknb2Nxmt/Uwd9LuPvZXiYYkp6LDMzE5iZFc/MTCO8/Xo1pNZGUPZ0QE879LR5pp5Xb4dnXofXzwPTjlOnvR3g7vd93xExxisy1hO6Q9/HQESsZxpz6jqDgRxzcrn3PFtEUMI1mCTI/aSi3sm/vrqTrdWtTMuI4ztXTefqOZnS3yjOaOB2rJVH2qk62uEJ7g6qjrQPDvMDyEq0MyMz3nhNimdmZgKFGbEj3yxKa6Pvtds5/KunDbrbhp8Ovm/3rbWrbBAVb7wi4yAqzmsaf/LnyFjPOp5Ajhx4H+NZHncyoMNkRMxoSJD7kdut+evOBn65oZIDTSeYm5PAP31uKlfPmSRXsk1gAy3sfUeNwDZC25gOXLYOxoN9p2UYgT1tUhwzM+OZlhpBgrsDuo6f+upu9Uyd0NXq+bnVE9Seqat35MLCwiEqAewJnmmiMY2K98zzBLM9AaISjUAeWD4Y0vFB6z4QZxbQIFdKLQP+C7ABj2utfzbS+lYP8gH9LjfrttXx6PtVHG7pxJEczVeWTOaGYof0oYewrl4XB5tPsL+pgwNNJ6hqMsL6QFMHPf0u4ukiWbUzOaaLOYn9FMb1kh/dRVZkF6lhJ7D3DgnrzmPGSbczUWFG+NqTIDrp5DQq4eTP9gTPNPFkUA+8j4iWAA4RAQtypZQNqASuBGqBz4BbtNa7zvSZUAnyAS63ZsOuIzzxtwN8dug4UeFhfP68LG4qyWXR5BTCpB/dclxuTX1rFwebOqg9cpSmhhqcLQ10Hj+C6mwilTZSlfHKjjhBhq2DZNqI6W/Fps/QR6xsxsiF6BTPNPnkKybFCOKBaXTSyWWR8dIFIYDADj+8AKjSWh/w7Og5YAVwxiAPNbYwxbK5mSybm8nOOifPfVbNn8vqWbetjqxEO58/L4vlRVnMdyRJqI8Xfd242htpbqjm2JEa2ptr6TnegLvjCOFdzcT1tpCqnFyAk0vUkCFpnv9suSITULFphMWmQawDYlIhNs2YDr7SjHCOSTFaydIyFgHgjxb5DcAyrfUdnp9vBxZprb8+ZL27gLsA8vLyFh4+fHhM+x3vunpdvL2rkde2N7Cxsolel5v0+CgunZ7O5TMzuGhqGonR0v3id24XdByBtgZor6frWC0dTTX0Hq9FtzUS0XWU2N5m4tztp39UK5xhCXREpNFnT0PFZxCVmEl8WjZxKZmouAyITTfCOTbN6DMWIohMvyBIa/0Y8BgYXSvB2q9ZoiNtrJifw4r5ObR19/Hu7iO8u/sob1U08uKWWpSCOdkJLJ6SSklBCgtyk8hIsJtd9vimNXS2QGs1OGvRzlq6mw/Tc6wG7awn4kQ9MT3NhHFyBEY0EK5tHCWJJp1MS1g63VGzccdNIjwxi9iUHJImOcjMKSA9I5vk8Ajkki9hNf4I8jog1+tnh2ee8EiwR3DdAgfXLXDQ73KztbqVj/c3s+lAC3/45DC///AgANmJdoocSczONi7ymJkZT05S9MTpjtEaTjTD8UP0HztE55ED9LYcRLfWENleS0x3AxHunsHVFaB0BMd0Kg06lQZm0KwuoicmEx2fRWSyg7j0XNIm5ZCbEsfklGgWxITAVY5CDOGPIP8MmKaUmowR4F8C/tEP2w1J4bYwLpicwgWTUwDjIbQV9U62VbeyraaVijon6ysaB9e3R4QxJS2Oyemx5KfEkJsSQ25yDJmJdrIS7Za5/4vbrWnv6ae1o5P2I4foa6rC3XIAW+shojuqie+qIaW3Abs27hESDiQALTqeOp1Onc6ggdm0RWXRG5uNTnQQlZpPcuokspJjyEmKZnZSNEkxETKmX0w4Y04BrXW/UurrwFsYww/Xaq0rxlzZBGGPsLEwP4WF+SmD8zp6+j2XZhtD2vY3dbCzzslbOxsHH1g7ID4qnLT4KNLiIkmNjSI5NoKE6AiSoiOJs4cTF2UjLiqC6Agb9ogw7BE2IsPDiLCFER6mCLcpwpRCKVAotNZowK01/S6Ny63pd7vp6XfT229Mu/tcdPcZ085eF529/ZzocdHR00dHTz9tXX2oE0dJPHGI1O7DpPfWku2qo4AGctVR8tXJro8uHUkNGey3ZVFqn0dnTA698bmEJRcQlT6ZlORkJiXYOT/BTmps5MT534kQoyAXBFmIy61pcHZRe7yLRmc3Dc5ujrR109TRQ0tHD80dvbR29uHs6h28xDuQFG5yVAvTVC3TVC2zwuuZGlbPZF1PHCcG1+tVUbRG53IiNp+exMnolCmEp00lJnMayZNyiYmSk75C+ML0k51i7GxhCkdyDI7kke+eprWmq89FR3c/HT1Ga7mrz+VpSbvodbnpc7np87S4tTZa4FprlKd1HqbUYIvdFhZGrKuNlPZ9JLbvJb5tH7Gte7EfrySs72RgEzcJ0qZD2iWQPgPSpkHqNCITcsiQcdBCBIwEeQhSShETGU5MZDgZo/2w1uCshYbt0LgDGnYY0zav89fRKTBpDhTcBhkzIX2WEdwxKWferhAiYCTIJ7qOo1C3Feq2QP1WqC+DzmZjmQqD1GmQtxgyz4PMuTBprtHylhOKQowbEuQTiavfaF3XfAq1n0Htp8aYbDBCO30WTF8G2fMha77R6g7ATfCFEP4lQR7K+rqMwD70EVR/ArWlxhNUAOKzIfd8uOAuyCmBrCLjVqNCCMuRIA8l/b1QVwoH/hcObjTeu3oBZXSLLLgV8i6E3EWQ6DC7WiGEn0iQW5nW0LIf9r8LVe/Cob95WtzK6B5Z9FXIv9gI7+gkk4sVQgSKBLnV9PfC4Y+g8i3Y9xYcO2DMT5kC82+BKZdCwcXG7U+FEBOCBLkV9LTDvg2w5w3Y97bxiK5wO0y+BC68B6YuhZTJZlcphDCJBPl41dMBleuhYh1UvWM8lzEmDWavgJnLYfLnZESJEAKQIB9fXP1w4H3Y8QLsed14snlcJhSvhNl/b/R1h43wIF4hxIQkQT4eNFVC2dOw/TnjoQj2JCi6Gc670bgYRy5vF0KMQILcLP09sPs1+OxxY4y3ssH0q2H+rTDtSnn6jBDCZxLkwdbeCJ/+Hrb+AU40QfJkWPpjmHcLxE8yuzohhAVJkAfLkV3wyaNQ/iK4+oxL4S+4A6ZcLl0nQogxkSAPtLqtsPEXsPcNiIiBhatg0d2QWmh2ZUKIECFBHij1ZfDeg1C1AeyJcOm/GPc1kVu9CiH8TILc31r2GwFe8YpxdeUVD8D5d4I9wezKhBAhSoLcX7qd8MG/w6e/A1sk/N13Yck3jda4EEIEkAT5WLndUPYMvPtjONEMxbfDZf9XRqAIIYJGgnwsmvfBX74J1R8bt4a99UXIXmB2VUKICUaC/Fy4+uDjR4yulAg7rPhv40IeefyZEMIEEuSjdewgvHyH8dCG2Svgmv+QbhQhhKkkyEdjx4vw+reN51vesBbmXm92RUIIIUHuk/4e+Ov3jMvqcy+E638PSXlmVyWEEIAE+dm1H4EXboeazfB334FLfwA2OWxCiPFDEmkk9dvg2X+E7la48UmYc53ZFQkhxGkkyM9k//vw3K0Qkwpr3obM88yuSAghhiVBPpyKV+GVOyF1Gtz+CsRnml2REEKckdw/daitT8GLqyC7GL7yhoS4EGLckxa5tx0vGldqTr0CbvqjPNxYCGEJEuQD9rwB674KBRfDzU9DRLTZFQkhhE+kawXgwAee7pQFcMuzEuJCCEuRIG+qhOduM05s3voiRMWbXZEQQozKxA7ybic8d4tx46tbX5Cn9wghLGlMQa6UulEpVaGUciulSvxVVFC4XfDynXD8ENz0FCQ6zK5ICCHOyVhb5DuBfwA2+qGW4Hr/32DfW3DNv0P+RWZXI4QQ52xMo1a01rsBlNXuw33ob/DhL2DBbVCyxuxqhBBiTILWR66UukspVaqUKm1qagrWbk/X0wGv3gMpU+Can8vDIIQQlnfWFrlS6h1guMsb79da/9nXHWmtHwMeAygpKdE+V+hvG/4VWqth9XqIjDWtDCGE8JezBrnWemkwCgmK/e9B6VpY/HXIu9DsaoQQwi8mzvDD3hPw529A2nS4/P+aXY0QQvjNWIcfXqeUqgUWA28opd7yT1kB8PGj0FYL1/5artwUQoSUsY5aWQes81MtgdN+BD76L5h1rXSpCCFCzsToWvngIXD1wNIfmV2JEEL4XegHedNe4x7jJWsgtdDsaoQQwu9CP8g3/NAYZvi5fza7EiGECIjQDvK6LVD5Jlx8L8Smml2NEEIERGgH+ab/gch4uOAusysRQoiACd0gb2uAinXG/VTkHuNCiBAWukFeuta4Ve0iaY0LIUJbaAZ5X7cR5NOXGTfHEkKIEBaaQV7xCnQ2w4V3m12JEEIEXOgFudbGSc70WTD5c2ZXI4QQARd6QV77GTTugEVflXuNCyEmhNAL8opXwRYJc683uxIhhAiK0ApyrWH3a1B4OdgTzK5GCCGCIrSCvGE7OKth1hfNrkQIIYImtIJ892ugbDD9GrMrEUKIoAm9IC9YIvdVEUJMKKET5E17oXmv8fAIIYSYQEInyHf/xZjOXG5uHUIIEWQhFOSvgeN8SMg2uxIhhAiq0Ajy44eNESsyWkUIMQGFRpBXvWNMZ0i3ihBi4gmNIK8thdh0eSanEGJCCo0gryuFnBK5t4oQYkKyfpB3HYfmSnCUmF2JEEKYwvpBXrfVmEqQCyEmKOsHeW0poCC72OxKhBDCFNYP8rpSSJ8pdzsUQkxY1g5yrY0WuWOh2ZUIIYRprB3kxw5A1zFjxIoQQkxQ1g7yui3GVE50CiEmMGsHee1nEBFrPGhZCCEmKIsHeSlkLwBbuNmVCCGEaawb5H3d0FguJzqFEBOedYO8cQe4+4xb1wohxAQ2piBXSv2HUmqPUmqHUmqdUirJT3Wd3cCJThmxIoSY4MbaIt8AzNVaFwGVwL+MvSQfHTsIUQmQkBW0XQohxHg0piDXWr+tte73/LgJcIy9JB85ayAxeLsTQojxyp995KuBN8+0UCl1l1KqVClV2tTUNPa9SZALIQTgQ5Arpd5RSu0c5rXCa537gX7gmTNtR2v9mNa6RGtdkp6ePvbKW2sgMXfs2xFCCIs76wBsrfXSkZYrpVYBXwCu0FprP9U1sp526G6VFrkQQuBDkI9EKbUM+D/A57TWnf4pyQfOOmOalBe0XQohxHg11j7yR4F4YINSqkwp9Vs/1HR2zhpjKi1yIYQYW4tcaz3VX4WMymCQSx+5EEJY88rO1hpQNojPNLsSIYQwnTWD3FkLCTkQZjO7EiGEMJ11gzxJulWEEAIsG+RyMZAQQgywXpC7+qGtXk50CiGEh/WCvKMRtEta5EII4WG9IG/1DD2UPnIhhACsGOTOWmMqXStCCAFYMsirjal0rQghBGDJIK+F6BSIjDW7EiGEGBesGeTSGhdCiEHWC/LWGrnroRBCeLFWkGstFwMJIcQQ1gry7lbo7ZARK0II4cVaQT449FBa5EIIMcCiQS4tciGEGGCtIJerOoUQ4jTWCnJnDdiiICbN7EqEEGLcsFaQp06FohshzFplCyFEII3pmZ1Bt/DLxksIIcQgadoKIYTFSZALIYTFSZALIYTFSZALIYTFSZALIYTFSZALIYTFSZALIYTFSZALIYTFKa118HeqVBNweBQfSQOaA1TOWIzXumD81jZe64LxW9t4rQvGb23jtS4YW235Wuv0oTNNCfLRUkqVaq1LzK5jqPFaF4zf2sZrXTB+axuvdcH4rW281gWBqU26VoQQwuIkyIUQwuKsEuSPmV3AGYzXumD81jZe64LxW9t4rQvGb23jtS4IQG2W6CMXQghxZlZpkQshhDgDCXIhhLA404NcKbVMKbVXKVWllPr+MMujlFLPe5ZvVkoVeC37F8/8vUqpq4Nc131KqV1KqR1KqXeVUvley1xKqTLP6y9BrmuVUqrJa/93eC37slJqn+fl9yd0+FDbr7zqqlRKtXotC+QxW6uUOqqU2nmG5Uop9Yin7h1KqWKvZQE7Zj7UdaunnnKl1MdKqXleyw555pcppUr9WZePtV2qlHJ6/c4e8Fo24vcgwHV9z6umnZ7vVYpnWcCOmVIqVyn1vicTKpRS3xpmncB9z7TWpr0AG7AfmAJEAtuB2UPWuQf4ref9l4DnPe9ne9aPAiZ7tmMLYl2XATGe9/80UJfn5w4Tj9cq4NFhPpsCHPBMkz3vk4NZ25D1vwGsDfQx82z7EqAY2HmG5Z8H3gQUcCGwOUjH7Gx1XTSwP+Cagbo8Px8C0kw8ZpcCr4/1e+Dvuoas+0XgvWAcMyALKPa8jwcqh/m7GbDvmdkt8guAKq31Aa11L/AcsGLIOiuAP3jevwRcoZRSnvnPaa17tNYHgSrP9oJSl9b6fa11p+fHTYDDT/seU10juBrYoLU+prU+DmwAlplY2y3As37c/xlprTcCx0ZYZQXwlDZsApKUUlkE+JidrS6t9cee/ULwvmMD+z7bMTuTsXxH/V1XML9jDVrrrZ737cBuIGfIagH7npkd5DlAjdfPtZz+hx9cR2vdDziBVB8/G8i6vK3B+Jd2gF0pVaqU2qSU+ns/1TSauq73/NftJaVU7ig/G+ja8HRDTQbe85odqGPmizPVHuhjNhpDv2MaeFsptUUpdZdJNS1WSm1XSr2plJrjmTcujplSKgYjDF/2mh2UY6aM7t8FwOYhiwL2PbPWw5fHIaXUbUAJ8Dmv2fla6zql1BTgPaVUudZ6f5BKeg14Vmvdo5T6Ksb/Zi4P0r599SXgJa21y2uemcdsXFNKXYYR5Bd7zb7Yc7wygA1KqT2e1mqwbMX4nXUopT4PvApMC+L+z+aLwEdaa+/We8CPmVIqDuMfj3u11m3+3PZIzG6R1wG5Xj87PPOGXUcpFQ4kAi0+fjaQdaGUWgrcD1yrte4ZmK+1rvNMDwAfYPzrHJS6tNYtXrU8Diz09bOBrs3LlxjyX94AHjNfnKn2QB+zs1JKFWH8HldorVsG5nsdr6PAOvzXregTrXWb1rrD8/6vQIRSKo1xcMw8RvqOBeSYKaUiMEL8Ga31K8OsErjvWSA6/kdxgiAco2N/MidPjMwZss7XOPVk5wue93M49WTnAfx3stOXuhZgnNSZNmR+MhDleZ8G7MNPJ3t8rCvL6/11wCZ98oTKQU99yZ73KcH8XXrWm4lx0kkF45h57aOAM5+4W86pJ6E+DcYx86GuPIxzPxcNmR8LxHu9/xhY5s+6fKgtc+B3iBGI1Z7j59P3IFB1eZYnYvSjxwbrmHn+7E8B/znCOgH7nvn1F3+OB+DzGGd49wP3e+b9BKOVC2AHXvR8oT8Fpnh99n7P5/YC1wS5rneAI0CZ5/UXz/yLgHLPF7gcWBPkuh4CKjz7fx+Y6fXZ1Z7jWAV8Jdi/S8/PPwJ+NuRzgT5mzwINQB9G/+Ma4G7gbs9yBfy3p+5yoCQYx8yHuh4Hjnt9x0o986d4jtV2z+/6/gD8Ls9W29e9vmeb8PrHZrjvQbDq8qyzCmMghPfnAnrMMLq9NLDD6/f1+WB9z+QSfSGEsDiz+8iFEEKMkQS5EEJYnAS5EEJYnAS5EEJYnAS5EEJYnAS5EEJYnAS5EEJY3P8Hs2jmkNtmyq4AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def f(x):\n",
    "    return x**2 + np.log(x)\n",
    "\n",
    "def L(x):\n",
    "    return 3*x - 2\n",
    "\n",
    "def Q(x):\n",
    "    return x**2/2 + 2*x - 3/2\n",
    "\n",
    "def eL(x):\n",
    "    return f(x) - L(x)\n",
    "\n",
    "def eQ(x):\n",
    "    return f(x) - Q(x)\n",
    "\n",
    "x = np.linspace(0, 2, 100)\n",
    "while 0 in x:\n",
    "    x = np.delete(x, np.where(x == 0)) # delete 0 from x because log(0) is undefined\n",
    "while 1 in x:\n",
    "    x = np.delete(x, np.where(x == 1)) # delete 1 from x because denominator will be 0 for eL and eQ\n",
    "\n",
    "eL = eL(x) / (x-1)\n",
    "eQ = eQ(x) / (x-1)**2\n",
    "\n",
    "plt.plot(x, eL, label='eL(x)/(x-1)')\n",
    "plt.plot(x, eQ, label='eQ(x)/(x-1)^2')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problem 2 - Combinational Descend Method\n",
    "Write a code to compute the unconstrained minimum of the following optimization problem by implementing a combination descent algorithm with initial point(1, 1)\n",
    "\n",
    "$f(x, y) = (x^2 - 3y^2)^2 + (\\sin(x^2 + y^2))^2$\n",
    "### Algorithm: Line Search Subroutine\n",
    "Input: \n",
    "- x: point in the domain of f \n",
    "- d: descent direction for f at x \n",
    "- β: reduction factor (0 < β < 1) \n",
    "- α: acceptable descent parameter (0 < α < 0.5)\n",
    "\n",
    "Output: \n",
    "- t: stepsize \n",
    "\n",
    "Algorithm: \n",
    "```python\n",
    "t := 1\n",
    "while f(x) − f(x + td) < α * Df(x)(td):\n",
    "    t := β * t\n",
    "return t as the stepsize\n",
    "```\n",
    "### Algorithm: Combinational Descend Method\n",
    "Input: \n",
    "- x: initial point in the domain of f \n",
    "- η: tolerance level for stopping criterion\n",
    "\n",
    "Output: \n",
    "- x: minimum point \n",
    "\n",
    "Algorithm: \n",
    "```python\n",
    "while norm(∇f(x)) > η:\n",
    "    if Hf(x) is positive definite:\n",
    "        d = [Hf(x)]^-1 * (-∇f(x))\n",
    "    else:\n",
    "        d = -∇f(x)\n",
    "    Line search: choose a stepsize t > 0 using the backtracking line search subroutine\n",
    "    x := x + td\n",
    "Return x as the minimum.\n",
    "```\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1.53499005 0.88622693]\n",
      "3.5244023582564334e-15\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import numdifftools as nd\n",
    "\n",
    "def f(var):\n",
    "    return (var[0]**2 - 3*var[1]**2)**2 + (np.sin(var[0]**2 + var[1]**2))**2 # function to be minimized\n",
    "\n",
    "def lineSearchSubroutine(f, x, d, beta=0.5, alpha=0.5):\n",
    "    t = 1\n",
    "    while f(x) - f(x + t * d) < -t*alpha*np.dot(nd.Gradient(f)(x), d): # line search subroutine to find the step size\n",
    "        t = beta*t\n",
    "    return t\n",
    "\n",
    "def combinationDescent(f, x, tolerance):\n",
    "    gradient = nd.Gradient(f)(x) # gradient of the function\n",
    "    while np.linalg.norm(gradient) > tolerance: # checking if norm of gradient is less than or equal to tolerance\n",
    "        hessian = nd.Hessian(f)(x) # hessian of the function\n",
    "        evals , evecs = np.linalg.eig(hessian) # eigen values and eigen vectors of the hessian\n",
    "        test_evals = np.all(evals > 0) # checking if all eigen values are positive\n",
    "        if test_evals:\n",
    "            d = -np.linalg.inv(hessian)@gradient # if all eigen values are positive then matrix multiplication of inverse of hessian and gradient\n",
    "        else :\n",
    "            d = -gradient # if all eigen values are not positive then d is equal to negative of gradient\n",
    "        t = lineSearchSubroutine(f, x, d) # finding the step size\n",
    "        x = x + t*d\n",
    "        gradient = nd.Gradient(f)(x) # updating the gradient of the function\n",
    "    return x\n",
    "\n",
    "x = np.array([1.0, 1.0])\n",
    "tolerance = 1e-6\n",
    "y = combinationDescent(f, x, tolerance)\n",
    "print(y)\n",
    "print(f(y))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problem 3\n",
    "Given the set S defined as:\n",
    "$$ S = {x \\in \\mathbb{R}^2: g_1(x_1,x_2) \\ge 0, g_2(x_1,x_2) \\ge 0, g_3(x_1,x_2) \\ge 0} $$\n",
    "where:<br>\n",
    "$g_1(x_1,x_2) = 2x_2 - x_1$<br>\n",
    "$g_2(x_1,x_2) = 2x_1 - x_2$<br>\n",
    "$g_3(x_1,x_2) = 1 - x_1 - x_2$<br>\n",
    "The associated potential function is defined as:\n",
    "$$ \\Psi(x) = -\\sum_{j=1}^{m} \\log g_j(x) $$\n",
    "The goal is to find the analytic center of the set S, which is the vector that minimizes the potential function, i.e.,\n",
    "$$ x_{analytic_center} = \\arg\\min_{x} \\Psi(x) $$\n",
    "### Solution\n",
    "$\\Psi(x_1, x_2) =$ $-\\Sigma \\log(g_i) =$\n",
    "$-\\log(g_1) - \\log(g_2) - \\log(g_3) =$\n",
    "$-\\log(g_1 \\cdot g_2 \\cdot g_3) =$\n",
    "$\\log\\left(\\frac{1}{g_1 \\cdot g_2 \\cdot g_3}\\right)$ <br>\n",
    "<br>\n",
    "This will be the new f(x). Beyond this, the process will be exactly the same as the previous problem(2)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.33333333 0.33333333]\n",
      "3.295836866004329\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import numdifftools as nd\n",
    "\n",
    "def g1(x):\n",
    "    return 2*x[1] - x[0]\n",
    "\n",
    "def g2(x):\n",
    "    return 2*x[0] - x[1]\n",
    "\n",
    "def g3(x):\n",
    "    return 1 - x[0] - x[1]\n",
    "\n",
    "def f(x):\n",
    "    return np.log(1/(g1(x)*g2(x)*g3(x))) # function to be minimized\n",
    "\n",
    "def lineSearchSubroutine(f, x, d, beta=0.5, alpha=0.5):\n",
    "    t = 1\n",
    "    while f(x) - f(x + t * d) < -t*alpha*np.dot(nd.Gradient(f)(x), d): # line search subroutine to find the step size\n",
    "        t = beta*t\n",
    "    return t\n",
    "\n",
    "def combinationDescent(f, x, tolerance):\n",
    "    gradient = nd.Gradient(f)(x) # gradient of the function\n",
    "    while np.linalg.norm(gradient) > tolerance: # checking if norm of gradient is less than or equal to tolerance\n",
    "        hessian = nd.Hessian(f)(x) # hessian of the function\n",
    "        evals , evecs = np.linalg.eig(hessian) # eigen values and eigen vectors of the hessian\n",
    "        test_evals = np.all(evals > 0) # checking if all eigen values are positive\n",
    "        if test_evals:\n",
    "            d = -np.linalg.inv(hessian)@gradient # if all eigen values are positive then matrix multiplication of inverse of hessian and gradient\n",
    "        else :\n",
    "            d = -gradient # if all eigen values are not positive then d is equal to negative of gradient\n",
    "        t = lineSearchSubroutine(f, x, d) # finding the step size\n",
    "        x = x + t*d\n",
    "        gradient = nd.Gradient(f)(x) # updating the gradient of the function\n",
    "    return x\n",
    "\n",
    "x = np.array([0.25, 0.25])\n",
    "tolerance = 1e-6 # tolerance value\n",
    "print(combinationDescent(f, x, tolerance))\n",
    "print(f(combinationDescent(f, x, tolerance)))"
   ]
  }
 ],
 "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.9.13"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "26a3cffd5f7a2724922a6aa05ca9c5d61dc993dd43746af46199d47f632ed425"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}