Skip to content

Instantly share code, notes, and snippets.

@adcroft
Last active December 1, 2023 15:55
Show Gist options
  • Save adcroft/85fc4ace76ba87d9635a96e6f058a826 to your computer and use it in GitHub Desktop.
Save adcroft/85fc4ace76ba87d9635a96e6f058a826 to your computer and use it in GitHub Desktop.
Interpolation with Lagrange cubic polynomial
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"import numpy"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"See https://en.wikipedia.org/wiki/Lagrange_polynomial"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"def lagrange_interp(x, y, q, e=1e-27):\n",
" \"\"\"Lagrange polynomial interpolation. Retruns f(q) which f(x) passes through four data\n",
" points at x[0..3], y[0..3].\n",
" \n",
" Optional argument \"e\" is a neglegible cell width to avoid division by zero for vanishing cells.\n",
" \"\"\"\n",
" # n - numerator, d - denominator\n",
" n0 = ( q - x[1] ) * ( q - x[2] ) * ( q - x[3] )\n",
" d0 = ( x[0] - x[1] - e ) * ( x[0] - x[2] - e ) * ( x[0] - x[3] - e )\n",
" n1 = ( q - x[0] ) * ( q - x[2] ) * ( q - x[3] )\n",
" d1 = ( x[1] - x[0] + e ) * ( x[1] - x[2] - e ) * ( x[1] - x[3] - e )\n",
" n2 = ( q - x[0] ) * ( q - x[1] ) * ( q - x[3] )\n",
" d2 = ( x[2] - x[0] + e ) * ( x[2] - x[1] + e ) * ( x[2] - x[3] - e )\n",
" n3 = ( q - x[0] ) * ( q - x[1] ) * ( q - x[2] )\n",
" d3 = ( x[3] - x[0] + e ) * ( x[3] - x[1] + e ) * ( x[3] - x[2] + e )\n",
" return ( (n0/d0)*y[0] + (n3/d3)*y[3] ) + ( (n1/d1)*y[1] + (n2/d2)*y[2] )"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAtdklEQVR4nO3df3xU9Z3v8ffk1wSVDIIQEhMhKIZfhUKCEkv8FQkllsoue9e712uxtW7TRRCy3GKwj9tHvd2N3XURUxCkplrKKrYdsLggS4ok0RpaAuFH5ac1khgT06hkYtAJSc79Y0xwyCRkhiTfzOT1fDzOg8f5zvcwn+/j64N5e873nGOzLMsSAACAIWGmCwAAAIMbYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAURGmC+iJtrY2ffDBBxo6dKhsNpvpcgAAQA9YlqXGxkbFx8crLKzr8x9BEUY++OADJSYmmi4DAAAEoKqqSgkJCV1+HhRhZOjQoZI8g4mJiTFcDQAA6AmXy6XExMSO3/GuBEUYab80ExMTQxgBACDIXGqJBQtYAQCAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGXFUby8vJks9m0bNmybvsVFxcrJSVF0dHRGjdunDZs2HA5XwsAAC7DU4WnlL/ntM/P8vec1lOFp/q1noDDyP79+7Vx40ZNnTq1234VFRXKyspSenq6ysvLtWrVKi1dulROpzPQrwYAAJchPMym1T4CSf6e01pdeErhYd0/vr23BfRumk8//VT33Xeffv7zn+snP/lJt303bNig6667TmvWrJEkTZw4UWVlZXryySe1cOHCQL4eAABchqUZ4yVJq784A7I0Y3xHEMmZc2PH5/0loDCyePFi3X333brrrrsuGUZKS0uVmZnp1TZ37lwVFBTo/PnzioyM7HSM2+2W2+3u2He5XIGUCQAAuvDlQLL29XfU3NpmJIhIAVym2bJliw4ePKi8vLwe9a+trVVsbKxXW2xsrFpaWlRfX+/zmLy8PDkcjo4tMTHR3zIBAMAlLM0Yr6jwMDW3tikqPMxIEJH8DCNVVVV65JFHtHnzZkVHR/f4uItfHWxZls/2drm5uWpoaOjYqqqq/CkTAAD0QP6e0x1BpLm1rctFrX3Nr8s0Bw4cUF1dnVJSUjraWltbVVJSorVr18rtdis8PNzrmNGjR6u2ttarra6uThERERoxYoTP77Hb7bLb7f6UBgAA/HDxGpH2fUkDe81IRkaGjh496tX27W9/WxMmTNDKlSs7BRFJSktL06uvvurVtnv3bqWmpvpcLwIAAPqWr8Wqvha19he/wsjQoUM1ZcoUr7Yrr7xSI0aM6GjPzc1VdXW1Nm3aJEnKzs7W2rVrlZOTo4ceekilpaUqKCjQSy+91EtDAAAA/mhts3wuVm3fb22z+rWegO6m6U5NTY0qKys79pOSkrRz504tX75c69atU3x8vPLz87mtFwAAQ5bPubHLz0wsYrVZ7atJBzCXyyWHw6GGhgbFxMSYLgcAAPRAT3+/eTcNAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKP8CiPr16/X1KlTFRMTo5iYGKWlpem1117rsn9RUZFsNlun7cSJE5ddOAAACA0R/nROSEjQE088oRtuuEGS9Mtf/lL33HOPysvLNXny5C6PO3nypGJiYjr2R44cGWC5AAAg1PgVRubPn++1/y//8i9av3699u3b120YGTVqlIYNGxZQgQAAILQFvGaktbVVW7ZsUVNTk9LS0rrtO336dMXFxSkjI0N79+695N/tdrvlcrm8NgAAEJr8DiNHjx7VVVddJbvdruzsbG3btk2TJk3y2TcuLk4bN26U0+nU1q1blZycrIyMDJWUlHT7HXl5eXI4HB1bYmKiv2UCAIAgYbMsy/LngObmZlVWVurs2bNyOp167rnnVFxc3GUgudj8+fNls9m0ffv2Lvu43W653e6OfZfLpcTERDU0NHitPQEAAAOXy+WSw+G45O+3X2tGJCkqKqpjAWtqaqr279+vp59+Ws8++2yPjp81a5Y2b97cbR+73S673e5vaQAAIAhd9nNGLMvyOotxKeXl5YqLi7vcrwUAACHCrzMjq1at0rx585SYmKjGxkZt2bJFRUVF2rVrlyQpNzdX1dXV2rRpkyRpzZo1Gjt2rCZPnqzm5mZt3rxZTqdTTqez90cCAACCkl9h5MMPP9T999+vmpoaORwOTZ06Vbt27dKcOXMkSTU1NaqsrOzo39zcrBUrVqi6ulpDhgzR5MmTtWPHDmVlZfXuKAAAQNDyewGrCT1dAAMAAAaOnv5+824aAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEb5FUbWr1+vqVOnKiYmRjExMUpLS9Nrr73W7THFxcVKSUlRdHS0xo0bpw0bNlxWwQAA854qPKX8Pad9fpa/57SeKjzVzxUhmPkVRhISEvTEE0+orKxMZWVluvPOO3XPPffo7bff9tm/oqJCWVlZSk9PV3l5uVatWqWlS5fK6XT2SvEAADPCw2xa7SOQ5O85rdWFpxQeZjNUGYKRzbIs63L+guHDh+vf//3f9eCDD3b6bOXKldq+fbuOHz/e0Zadna3Dhw+rtLS0x9/hcrnkcDjU0NCgmJiYyykXANBL2oNHzpwbtTRjfKd9oKe/3xGBfkFra6t+85vfqKmpSWlpaT77lJaWKjMz06tt7ty5Kigo0Pnz5xUZGenzOLfbLbfb3bHvcrkCLRMA0EeW3nG94hvK1bJ3hf5x703a3fJVgggC4ncYOXr0qNLS0vT555/rqquu0rZt2zRp0iSffWtraxUbG+vVFhsbq5aWFtXX1ysuLs7ncXl5efrxj3/sb2kAgP5Qf1o68mvpyMv6u7NnpAhpeGujisJnEEQQEL/vpklOTtahQ4e0b98+ff/739eiRYt07NixLvvbbN7XDduvCl3c/mW5ublqaGjo2KqqqvwtEwDQmxqqpbd+Jj17m7Q2VSr5N+nsGTWHX6mXW27XC9Y31Nza1uWiVqA7fp8ZiYqK0g033CBJSk1N1f79+/X000/r2Wef7dR39OjRqq2t9Wqrq6tTRESERowY0eV32O122e12f0sDAPSmxlrp2Hbp7a1S5ZfW+dnCpRvu0q6wdC07fK3+ac5X9OKX1oxI4gwJ/BLwmpF2lmV5re/4srS0NL366qtebbt371ZqamqX60UAAAY1VEvHX5WO/e6LAPKlexwSZ0lf+Ttp0gLl//Fsp8Wq7X8SSOAvv8LIqlWrNG/ePCUmJqqxsVFbtmxRUVGRdu3aJclzeaW6ulqbNm2S5LlzZu3atcrJydFDDz2k0tJSFRQU6KWXXur9kQAAAvPXU9KJ/5JO7JCqy7w/uzZVmvw30uQFkiOho7m17ROfi1Xb91vbLutGTQwyfoWRDz/8UPfff79qamrkcDg0depU7dq1S3PmzJEk1dTUqLKysqN/UlKSdu7cqeXLl2vdunWKj49Xfn6+Fi5c2LujAAD0XGuL9P6fpJM7pZOvSR+986UPbdJ1s6SJ86WJ35SGJfr8K5bPubHLv54zIvDXZT9npD/wnBEAuExNH0l/2SOd+m/pnd9Ln5+98FlYpDTuNmnC3VJyljR0tLEyEVr6/DkjAIABrK1Vqj7oCR7v/F6qPiCv9R9DrpbGZ0rJ86TrM6Ro/kcP5hBGACBUfHJGenev9JfXpXeLvc9+SNKoydKNmdL4uVLiTVJYuJEygYsRRgAgWDXVSxUlUkWxJ3x8UuH9ud3hufwyfo7n7IfjWjN1ApdAGAGAYHHuY+nMH6SKN6T33pTqLnpJqS1cSpgpXX+HNO4O6doUKZx/5jHw8V8pAAxUjbXSmbcubBeHD0mKnSKNTZfG3S6N/ZpkH9rvZQKXizACAAOBZXlusa3c98X2lvTxu537jZwgjZ3tCSBjZ0tXXtP/tQK9jDACACac/1z6oFyq+uOF7dxHF3Wyec58jLnli+1r0lUjjZQL9CXCCAD0NcuSzlZK7++X3i/zPHCs5ojUdt67X7hdunaGdF2aZ0u8SRoyzEjJQH8ijABAb3M3ep7xUV0mvX/AE0Ka6jr3u3KUdN3NUsJNnqeexk2TInhJKAYfwggAXI7W81LdMc9Dxd4/4Pnzryfk9YAxSQqLkEZ/xRM8EmZKiTOlYWMkm81I2cBAQhgBgJ6yLM+i0uqDntDxwUGp5rDU8nnnvo5Ez621Came8BE3TYoc0v81A0GAMAIAvliW5PrAEziqD3r+/KBc+ryhc99ohxQ/3fOG22tTPNvQ2P6vGQhShBEAkDwPFGsPHe1/fvph537hdiluqidwxM/w/Dl8nBQW1v81AyGCMAJg8HF/6rm80h48qg9IZ8907mcLl0ZN/OKsR4rnTpdRk6TwyP6vGQhhhBEAoa31vFR33BM4qg94wsdfj0tWW+e+I27wBI/4GZ7gMXqqFHVF/9cMDDKEEQChw7KkhirPszyqD3j+rDkstXzWue/QeE/guHaGJ3zET+eZHoAhhBEAwau5ybOotONhYvt9r/OwO6Rrp19YXBo/Q4qJ6/96AfhEGAEQHNqfYlr1pwuPT//wbclq9e4XFuF5hHpC6oW7W0bcwAJTYAAjjAAYmFpbpNojnpfGVe3zhJDGms79Yq69EDwSb+J5HkAQIowAGBjOf+ZZ53HmLenMH6Sq/dL5Ju8+YRGesJF48xdPMb1JciSYqRdAryGMADCj+ZznhXHvvenZqg9Irc3efaIdnuCReLPn3S3xM7i7BQhBhBEA/aOl2fPiuIoSz/b+/s7h46rR0pg0aczXPG+tHTWJtR7AIEAYAdA3LMvzArm/7JXeLfJcejl/zrvP0HgpKV0aO9sTQIaP48VxwCBEGAHQe5rqpb+8Lr2zR3p3b+fbbK+4Rkq69cJG+AAgwgiAy9HW5nnOx+ndnu2DcknWhc8jhnjOeoy73bNx2QWAD4QRAP5xN3rOfpz6b892rt7789gp0g0Z0vUZnkWnEXYzdQIIGoQRAJfWWCud3Cmd2ClVFHsvPI0aKl1/hzQ+U7rhLp5sCsBvhBEAvn1cIR1/VTq+3XPny5cNHyfdOE9K/rqUOEuKiDJTI4CQQBgBcMFHf5GOvSK9/Yrn6adfljBTSs6SJtwtXXMjC08B9BrCCDDYNbwv/dnp2WoOX2i3hXkWn078pjThG1x+AdBnCCPAYHTuY+nY76Qjv5Yq37rQbgv3PPdj8t94AsiV15irEcCg4VcYycvL09atW3XixAkNGTJEt9xyi376058qOTm5y2OKiop0xx13dGo/fvy4JkyY4H/FAALT0iy9UygdfslzF0zHIlSbNOYWacpCadI9BBAA/c6vMFJcXKzFixdr5syZamlp0WOPPabMzEwdO3ZMV155ZbfHnjx5UjExMR37I0eODKxiAP6p/bNUvlk6+mvp3EcX2mO/Ik39e08IcVxrrj4Ag55fYWTXrl1e+88//7xGjRqlAwcO6NZbb+322FGjRmnYsGF+FwggAJ+7pKO/kQ5ukmoOXWi/arQngEz7n1LsZGPlAcCXXdaakYaGBknS8OHDL9l3+vTp+vzzzzVp0iT98Ic/9Hnppp3b7Zbb7e7Yd7lcl1MmMHhUH5DKnvcsRm1/D0xYpDQhS/rq/5auv1MKZ6kYgIEl4H+VLMtSTk6OZs+erSlTpnTZLy4uThs3blRKSorcbrd+9atfKSMjQ0VFRV2eTcnLy9OPf/zjQEsDBpfmc9KffyvtL/A+C3JNspSySJp6L+tAAAxoNsuyrEt362zx4sXasWOH3nzzTSUkJPh17Pz582Wz2bR9+3afn/s6M5KYmKiGhgavdSfAoPbJGWn/z6WDv5I+P+tpC4+SJi2QUr/jeRQ7zwIBYJDL5ZLD4bjk73dAZ0aWLFmi7du3q6SkxO8gIkmzZs3S5s2bu/zcbrfLbud9FkAnliVVlkql6zyPZ7faPO3DxkgzH/RcirlyhNkaAcBPfoURy7K0ZMkSbdu2TUVFRUpKSgroS8vLyxUXxwOUgB5rbfE8GfWtn3lfihl3u3Rztue9MGHhhooDgMvjVxhZvHixXnzxRf3ud7/T0KFDVVtbK0lyOBwaMmSIJCk3N1fV1dXatGmTJGnNmjUaO3asJk+erObmZm3evFlOp1NOp7OXhwKEoOYmz2WYfeuks5WetohozzqQWd+XRk00Wx8A9AK/wsj69eslSbfffrtX+/PPP68HHnhAklRTU6PKysqOz5qbm7VixQpVV1dryJAhmjx5snbs2KGsrKzLqxwIZec+lv70c+mPG6TPPva0XXGNdNM/ei7HsCAVQAgJeAFrf+rpAhgg6H36V6l0rbT/Oan5U0/b1WOlW5ZIX71PihxitDwA8EefLmAF0Ms+rZP+8LTn9tyWzzxtsVOk2cs9d8fwbBAAIYx/4QCTmuqlP6yR/vTchRASP0O67QfSjV/n1lwAgwJhBDDhs7OeO2P2rZfON3nark2Rbs+VbriLEAJgUCGMAP3p/GfSH5+V3nzqwoPK4qZJd/xQGj+HEAJgUCKMAP2hrVU69KK091+lxg88bSMnSnc+Jk34BiEEwKBGGAH62ju/l3b/X6nubc++I1G64zHP23N5UBkAEEaAPlN3XNr9Q08YkaToYdKtK6SZD0mR0UZLA4CBhDAC9LZzH3sux5QVeN4dExYp3fw9Kf2fpSuGm64OAAYcwgjQW9papQMvSK//P+mzTzxtE74hzXlcGnG90dIAYCAjjAC9oepP0o5/lmqPePZHTZK+/oQ07jazdQFAECCMAJfj3MfS738kHfS8GFJ2h+cOmdQHeWoqAPQQ/1oCgbAs6fAWafdj0rmPPG1fvU+668fSVSPN1gYAQYYwAvjro79I/7VMqijx7I+aJN29WhqTZrQsAAhWhBGgp1pbpNKfSUVPSC2fSxFDpNtXSmkPS+GRpqsDgKBFGAF6ovbP0u/+Sao57Nkfd7v0jTXS8CSTVQFASCCMAN1pPS+9sVoq+TeprcXz4LKv50nT/oFHuANALyGMAF358G1pW/aF23UnfEO6+z+koaPN1gUAIYYwAlysrVUqXSu9/hOptdlzNuTu/5CmLORsCAD0AcII8GVnKz1nQ878wbN/49el+U9zNgQA+hBhBGh35Neep6i6XVLUVdLcf5VmfIuzIQDQxwgjwOcuTwg5+mvPfsJN0t8+Kw0fZ7YuABgkCCMY3KoPSL99UPqkQrKFSbetlNJX8Ch3AOhH/IuLwcmypH3PSIU/ktrOS45EaeFz0nWzTFcGAIMOYQSDz2efSK/8k3Ryp2d/4jelb+ZLQ642WxcADFKEEQwu1Qel3yzy3DUTHuVZpDrzuyxSBQCDCCMYHCxLOvC89NpKz7NDrh4r/Y9fSvFfNV0ZAAx6hBGEvvOfee6WOfSfnv3ku6UFz0hDhhktCwDgQRhBaDtbJb18n+cFd7YwKeNH0tce4bIMAAwghBGEroo3POtDzn0kXTFC+rvnpXG3ma4KAHARwghCj2VJ+5/zrA+xWqW4adK9m6Vh15muDADgA2EEoaWlWXrt/0gHXvDsf+XvPbftRg4xWhYAoGth/nTOy8vTzJkzNXToUI0aNUoLFizQyZMnL3lccXGxUlJSFB0drXHjxmnDhg0BFww8VXhK+XtOd/7g3Md6/2dzvwgiNmnO49LfbiSIAMAA51cYKS4u1uLFi7Vv3z4VFhaqpaVFmZmZampq6vKYiooKZWVlKT09XeXl5Vq1apWWLl0qp9N52cVjcAoPs2n1xYGk/rTO5qcroeGgmsOvlP7XyyxUBYAgYbMsywr04L/+9a8aNWqUiouLdeutt/rss3LlSm3fvl3Hjx/vaMvOztbhw4dVWlrao+9xuVxyOBxqaGhQTExMoOUihOTvOa3VhaeUM+dGLR1Xq8//8x8U3dKoBnu8HN9xSrGTTJcIAINeT3+/L2vNSENDgyRp+PDhXfYpLS1VZmamV9vcuXNVUFCg8+fPKzIystMxbrdbbre7Y9/lcl1OmQhBSzPGS5LO7HlOzSU/V7StVTVDpyrue1ulq0Yarg4A4A+/LtN8mWVZysnJ0ezZszVlypQu+9XW1io2NtarLTY2Vi0tLaqvr/d5TF5enhwOR8eWmJgYaJkIYUtnRCkv8ueKsrVqR1ua4pbuJogAQBAKOIw8/PDDOnLkiF566aVL9rVddN2+/crQxe3tcnNz1dDQ0LFVVVUFWiZCWP7BZj3W8qA2tN6jh5sXK7/kfdMlAQACENBlmiVLlmj79u0qKSlRQkJCt31Hjx6t2tpar7a6ujpFRERoxIgRPo+x2+2y2+2BlIZB4sKakX9UdsZ4NX+xL124hAMACA5+hRHLsrRkyRJt27ZNRUVFSkpKuuQxaWlpevXVV73adu/erdTUVJ/rRYBL8Vq8+kXwaP+TQAIAwcevyzSLFy/W5s2b9eKLL2ro0KGqra1VbW2tPvvss44+ubm5+ta3vtWxn52drTNnzignJ0fHjx/XL37xCxUUFGjFihW9NwoMKq1tllcQabc0Y7xy5tyo1raAbxADABjg1629Xa3xeP755/XAAw9Ikh544AG99957Kioq6vi8uLhYy5cv19tvv634+HitXLlS2dnZPS6SW3sBAAg+Pf39vqznjPQXwggAAMGnp7/fAd9NAwAA0BsIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwyu8wUlJSovnz5ys+Pl42m02vvPJKt/2Liopks9k6bSdOnAi0ZgAAEEIi/D2gqalJ06ZN07e//W0tXLiwx8edPHlSMTExHfsjR47096sBAEAI8juMzJs3T/PmzfP7i0aNGqVhw4b5fRwAAAht/bZmZPr06YqLi1NGRob27t3bbV+32y2Xy+W1AQCA0NTnYSQuLk4bN26U0+nU1q1blZycrIyMDJWUlHR5TF5enhwOR8eWmJjY12UCAABDbJZlWQEfbLNp27ZtWrBggV/HzZ8/XzabTdu3b/f5udvtltvt7th3uVxKTExUQ0OD17oTAAAwcLlcLjkcjkv+fhu5tXfWrFk6ffp0l5/b7XbFxMR4bQAAIDQZCSPl5eWKi4sz8dUAAGCA8ftumk8//VTvvPNOx35FRYUOHTqk4cOH67rrrlNubq6qq6u1adMmSdKaNWs0duxYTZ48Wc3Nzdq8ebOcTqecTmfvjQIAAAQtv8NIWVmZ7rjjjo79nJwcSdKiRYv0wgsvqKamRpWVlR2fNzc3a8WKFaqurtaQIUM0efJk7dixQ1lZWb1QPgAACHaXtYC1v/R0AQwAABg4BvQCVgAAgHaEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAY5XcYKSkp0fz58xUfHy+bzaZXXnnlkscUFxcrJSVF0dHRGjdunDZs2BBIrQAAIAT5HUaampo0bdo0rV27tkf9KyoqlJWVpfT0dJWXl2vVqlVaunSpnE6n38UCAIDQE+HvAfPmzdO8efN63H/Dhg267rrrtGbNGknSxIkTVVZWpieffFILFy709+sBAECI6fM1I6WlpcrMzPRqmzt3rsrKynT+/Hmfx7jdbrlcLq8NAACEpj4PI7W1tYqNjfVqi42NVUtLi+rr630ek5eXJ4fD0bElJib2dZkAAMCQfrmbxmazee1bluWzvV1ubq4aGho6tqqqqj6vEQAAmOH3mhF/jR49WrW1tV5tdXV1ioiI0IgRI3weY7fbZbfb+7o0AAAwAPT5mZG0tDQVFhZ6te3evVupqamKjIzs668HAAADnN9h5NNPP9WhQ4d06NAhSZ5bdw8dOqTKykpJnkss3/rWtzr6Z2dn68yZM8rJydHx48f1i1/8QgUFBVqxYkXvjAAAAAQ1vy/TlJWV6Y477ujYz8nJkSQtWrRIL7zwgmpqajqCiSQlJSVp586dWr58udatW6f4+Hjl5+dzWy8AAJAk2az21aQDmMvlksPhUENDg2JiYkyXAwAAeqCnv9+8mwYAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUQGFkWeeeUZJSUmKjo5WSkqK3njjjS77FhUVyWazddpOnDgRcNEAACB0+B1GXn75ZS1btkyPPfaYysvLlZ6ernnz5qmysrLb406ePKmampqObfz48QEXDQAAQoffYWT16tV68MEH9d3vflcTJ07UmjVrlJiYqPXr13d73KhRozR69OiOLTw8POCiAQBA6PArjDQ3N+vAgQPKzMz0as/MzNRbb73V7bHTp09XXFycMjIytHfv3m77ut1uuVwurw0AAIQmv8JIfX29WltbFRsb69UeGxur2tpan8fExcVp48aNcjqd2rp1q5KTk5WRkaGSkpIuvycvL08Oh6NjS0xM9KdMAAAQRCICOchms3ntW5bVqa1dcnKykpOTO/bT0tJUVVWlJ598UrfeeqvPY3Jzc5WTk9Ox73K5CCQAAIQov86MXHPNNQoPD+90FqSurq7T2ZLuzJo1S6dPn+7yc7vdrpiYGK8NAACEJr/CSFRUlFJSUlRYWOjVXlhYqFtuuaXHf095ebni4uL8+WoAABCi/L5Mk5OTo/vvv1+pqalKS0vTxo0bVVlZqezsbEmeSyzV1dXatGmTJGnNmjUaO3asJk+erObmZm3evFlOp1NOp7N3RwIAAIKS32Hk3nvv1UcffaTHH39cNTU1mjJlinbu3KkxY8ZIkmpqaryeOdLc3KwVK1aourpaQ4YM0eTJk7Vjxw5lZWX13igAAEDQslmWZZku4lJcLpccDocaGhpYPwIAQJDo6e8376YBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYFRAYeSZZ55RUlKSoqOjlZKSojfeeKPb/sXFxUpJSVF0dLTGjRunDRs2BFRsb3iq8JTy95z2+Vn+ntN6qvBUP1cEAMDg5ncYefnll7Vs2TI99thjKi8vV3p6uubNm6fKykqf/SsqKpSVlaX09HSVl5dr1apVWrp0qZxO52UXH4jwMJtW+wgk+XtOa3XhKYWH2YzUBQDAYGWzLMvy54Cbb75ZM2bM0Pr16zvaJk6cqAULFigvL69T/5UrV2r79u06fvx4R1t2drYOHz6s0tLSHn2ny+WSw+FQQ0ODYmJi/CnXp/bgkTPnRi3NGN9pHwAAXL6e/n5H+POXNjc368CBA3r00Ue92jMzM/XWW2/5PKa0tFSZmZlebXPnzlVBQYHOnz+vyMjITse43W653W6vwfSm9sCxuvCU1r7+jppb2wgiAAAY4tdlmvr6erW2tio2NtarPTY2VrW1tT6Pqa2t9dm/paVF9fX1Po/Jy8uTw+Ho2BITE/0ps0eWZoxXVHiYmlvbFBUeRhABAMCQgBaw2mze6yosy+rUdqn+vtrb5ebmqqGhoWOrqqoKpMxu5e853RFEmlvbulzUCgAA+pZfl2muueYahYeHdzoLUldX1+nsR7vRo0f77B8REaERI0b4PMZut8tut/tTml+6WjMiiTMkAAD0M7/OjERFRSklJUWFhYVe7YWFhbrlllt8HpOWltap/+7du5WamupzvUhf87VYdWnGeOXMudHnXTYAAKBv+XVmRJJycnJ0//33KzU1VWlpadq4caMqKyuVnZ0tyXOJpbq6Wps2bZLkuXNm7dq1ysnJ0UMPPaTS0lIVFBTopZde6t2R9FBrm+VzsWr7fmubXzcXAQCAy+R3GLn33nv10Ucf6fHHH1dNTY2mTJminTt3asyYMZKkmpoar2eOJCUlaefOnVq+fLnWrVun+Ph45efna+HChb03Cj8sn3Njl59xiQYAgP7n93NGTOjt54wAAIC+19Pfb95NAwAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIzy+3HwJrQ/JNblchmuBAAA9FT77/alHvYeFGGksbFRkpSYmGi4EgAA4K/GxkY5HI4uPw+Kd9O0tbXpgw8+0NChQ2Wz2Xrt73W5XEpMTFRVVVXIvvMm1MfI+IJfqI8x1Mcnhf4YGV/gLMtSY2Oj4uPjFRbW9cqQoDgzEhYWpoSEhD77+2NiYkLyP7AvC/UxMr7gF+pjDPXxSaE/RsYXmO7OiLRjASsAADCKMAIAAIwa1GHEbrfrRz/6kex2u+lS+kyoj5HxBb9QH2Ooj08K/TEyvr4XFAtYAQBA6BrUZ0YAAIB5hBEAAGAUYQQAABhFGAEAAEaFfBh55plnlJSUpOjoaKWkpOiNN97otn9xcbFSUlIUHR2tcePGacOGDf1UaeD8GWNRUZFsNlun7cSJE/1Ycc+VlJRo/vz5io+Pl81m0yuvvHLJY4JpDv0dX7DNX15enmbOnKmhQ4dq1KhRWrBggU6ePHnJ44JlDgMZX7DN4fr16zV16tSOB2KlpaXptdde6/aYYJk/yf/xBdv8XSwvL082m03Lli3rtl9/z2FIh5GXX35Zy5Yt02OPPaby8nKlp6dr3rx5qqys9Nm/oqJCWVlZSk9PV3l5uVatWqWlS5fK6XT2c+U95+8Y2508eVI1NTUd2/jx4/upYv80NTVp2rRpWrt2bY/6B9sc+ju+dsEyf8XFxVq8eLH27dunwsJCtbS0KDMzU01NTV0eE0xzGMj42gXLHCYkJOiJJ55QWVmZysrKdOedd+qee+7R22+/7bN/MM2f5P/42gXL/H3Z/v37tXHjRk2dOrXbfkbm0AphN910k5Wdne3VNmHCBOvRRx/12f8HP/iBNWHCBK+2733ve9asWbP6rMbL5e8Y9+7da0myPvnkk36orndJsrZt29Ztn2Ccw3Y9GV8wz59lWVZdXZ0lySouLu6yTzDPYU/GF+xzaFmWdfXVV1vPPfecz8+Cef7adTe+YJ2/xsZGa/z48VZhYaF12223WY888kiXfU3MYcieGWlubtaBAweUmZnp1Z6Zmam33nrL5zGlpaWd+s+dO1dlZWU6f/58n9UaqEDG2G769OmKi4tTRkaG9u7d25dl9qtgm8NABev8NTQ0SJKGDx/eZZ9gnsOejK9dMM5ha2urtmzZoqamJqWlpfnsE8zz15PxtQu2+Vu8eLHuvvtu3XXXXZfsa2IOQzaM1NfXq7W1VbGxsV7tsbGxqq2t9XlMbW2tz/4tLS2qr6/vs1oDFcgY4+LitHHjRjmdTm3dulXJycnKyMhQSUlJf5Tc54JtDv0VzPNnWZZycnI0e/ZsTZkypct+wTqHPR1fMM7h0aNHddVVV8lutys7O1vbtm3TpEmTfPYNxvnzZ3zBOH9btmzRwYMHlZeX16P+JuYwKN7aezlsNpvXvmVZndou1d9X+0DizxiTk5OVnJzcsZ+Wlqaqqio9+eSTuvXWW/u0zv4SjHPYU8E8fw8//LCOHDmiN99885J9g3EOezq+YJzD5ORkHTp0SGfPnpXT6dSiRYtUXFzc5Q92sM2fP+MLtvmrqqrSI488ot27dys6OrrHx/X3HIbsmZFrrrlG4eHhnc4Q1NXVdUp87UaPHu2zf0REhEaMGNFntQYqkDH6MmvWLJ0+fbq3yzMi2OawNwTD/C1ZskTbt2/X3r17lZCQ0G3fYJxDf8bny0Cfw6ioKN1www1KTU1VXl6epk2bpqefftpn32CcP3/G58tAnr8DBw6orq5OKSkpioiIUEREhIqLi5Wfn6+IiAi1trZ2OsbEHIZsGImKilJKSooKCwu92gsLC3XLLbf4PCYtLa1T/927dys1NVWRkZF9VmugAhmjL+Xl5YqLi+vt8owItjnsDQN5/izL0sMPP6ytW7fq9ddfV1JS0iWPCaY5DGR8vgzkOfTFsiy53W6fnwXT/HWlu/H5MpDnLyMjQ0ePHtWhQ4c6ttTUVN133306dOiQwsPDOx1jZA77bGnsALBlyxYrMjLSKigosI4dO2YtW7bMuvLKK6333nvPsizLevTRR63777+/o/+7775rXXHFFdby5cutY8eOWQUFBVZkZKT129/+1tQQLsnfMT711FPWtm3brFOnTll//vOfrUcffdSSZDmdTlND6FZjY6NVXl5ulZeXW5Ks1atXW+Xl5daZM2csywr+OfR3fME2f9///vcth8NhFRUVWTU1NR3buXPnOvoE8xwGMr5gm8Pc3FyrpKTEqqiosI4cOWKtWrXKCgsLs3bv3m1ZVnDPn2X5P75gmz9fLr6bZiDMYUiHEcuyrHXr1lljxoyxoqKirBkzZnjdcrdo0SLrtttu8+pfVFRkTZ8+3YqKirLGjh1rrV+/vp8r9p8/Y/zpT39qXX/99VZ0dLR19dVXW7Nnz7Z27NhhoOqeab+N7uJt0aJFlmUF/xz6O75gmz9fY5NkPf/88x19gnkOAxlfsM3hd77znY5/X0aOHGllZGR0/FBbVnDPn2X5P75gmz9fLg4jA2EObZb1xaoUAAAAA0J2zQgAAAgOhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABG/X/JYDQX7hMgPwAAAABJRU5ErkJggg==",
"text/plain": [
"<Figure size 640x480 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"x = [0., 1., 3., 4.] # Positiong of data, x[i]\n",
"y = [0., 2., 3., 4.] # Values of data, y[i]\n",
"q = numpy.linspace(x[1],x[2],100) # Sample interval between x[1]..x[2]\n",
"plt.plot(x,y,'x');\n",
"plt.plot(q, lagrange_interp(x,y,q));"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAABCuklEQVR4nO3de1xUZf4H8M+ZGWa4o4giCgKiIkreQBHQTE1aMzdzS63WS2mr29Wo/EluF+1Ctaa1W7CamlmWVl62LSypvKBkCuIVb3kDEURQGS4yw8w8vz9QCrnIIHDmDJ/36zWvbc48B77PPul8es5zniMJIQSIiIiIZKKSuwAiIiJq3RhGiIiISFYMI0RERCQrhhEiIiKSFcMIERERyYphhIiIiGTFMEJERESyYhghIiIiWWnkLqAhLBYLzp8/Dzc3N0iSJHc5RERE1ABCCBQXF6NTp05Qqeqe/1BEGDl//jz8/PzkLoOIiIgaITs7G76+vnV+rogw4ubmBqCyM+7u7jJXQ0RERA2h1+vh5+dX9T1eF0WEkeuXZtzd3RlGiIiIFOZmSyy4gJWIiIhkxTBCREREsmIYISIiIlkxjBAREZGsGEaIiIhIVgwjREREJCuGESIiIpIVwwgRERHJimGEiIiIZMUwQkRE1NpsiQe2vVP7Z9veqfy8BTGMEBERtTYqNbDljZqBZNs7lcdV6hYtRxHPpiEiIqImNGxO5f9ueeP399eDyPB5v3/eQhhGiIiIWqM/BpLt/wTMRlmCCMDLNERERK2WZegLqIAGMBsh1FpZggjAMEJERNRqqVL+CQeYYJIcIJmNdS9qbe46ZPmtREREJK8/rBHRvFJQeYmmtkWtLYBrRoiIiFqZC/9bAO/0d6uvEaltUWsLYRghIiJqRUoNJmw6kIPCivvh5/IgJvzxw+sBxGJu0ZoYRoiIiFqRV745jK+L/wwfD0ds6uVds4FS7qZJSEhAYGAgHB0dERYWhpSUlHrbf/jhhwgJCYGTkxOCg4OxatWqRhVLREREjbch4xy+Tj8HlQS8N7Ef2jhr5S4JQCNmRtauXYvZs2cjISEB0dHRWLJkCUaPHo3MzEx06dKlRvvExETExcXho48+wsCBA7F792489thjaNu2LcaOHdsknSAiIqL6/ZZfjBfXHwIAPDWiOyK6tpO5ot9JQghhzQkREREYMGAAEhMTq46FhIRg3LhxiI+vuZd9VFQUoqOj8c9//rPq2OzZs5GWloYdO3Y06Hfq9Xp4eHigqKgI7u7u1pRLRETU6pUZTRj34U4cv1CC6G7tsOrRCKhVUrP/3oZ+f1t1mcZoNCI9PR0xMTHVjsfExCA1NbXWcwwGAxwdHasdc3Jywu7du1FRUVHnOXq9vtqLiIiIGic58wKOXyhBezcd3pvYv0WCiDWsCiMFBQUwm83w9q6+4MXb2xt5eXm1nnPXXXdh2bJlSE9PhxACaWlpWLFiBSoqKlBQUFDrOfHx8fDw8Kh6+fn5WVMmERER/cG9/Trjw4cG4F+T+qO9m07ucmpo1AJWSaqeqIQQNY5d99JLL2H06NEYPHgwHBwccO+992LatGkAALW69qcCxsXFoaioqOqVnZ3dmDKJiIjomjF9fBAZZDvrRP7IqjDi5eUFtVpdYxYkPz+/xmzJdU5OTlixYgXKyspw5swZZGVlISAgAG5ubvDy8qr1HJ1OB3d392ovIiIiargSgwn/9/UB5OvL5S7lpqwKI1qtFmFhYUhOTq52PDk5GVFRUfWe6+DgAF9fX6jVaqxZswb33HMPVCruRk9ERNTUhBD4v68PYG1aNh79ZA+svFelxVl9a29sbCwmT56M8PBwREZGYunSpcjKysKsWbMAVF5iycnJqdpL5Pjx49i9ezciIiJw+fJlLFq0CIcOHcInn3zStD0hIiIiAMCylNP47mAuHNQS5v85tM6lFLbC6jAyceJEFBYWYsGCBcjNzUVoaCiSkpLg7+8PAMjNzUVWVlZVe7PZjHfffRfHjh2Dg4MDhg8fjtTUVAQEBDRZJ4iIiKjSLycL8db3RwEAL93TC2H+bWWu6Oas3mdEDtxnhIiI6Obyispxz79TUFBixPj+nfHuhL6yzoo0yz4jREREZJuMJgseX52OghIjenZ0wxv33Wbzl2euYxghIiKyA5fLjCgxmODuqMGSyWFw0ta+fYYt4lN7iYiI7IC3uyPWPx6NExeK4d/ORe5yrMKZESIiIgUrLv/90SquOg36d7H9Bas3YhghIiJSqOxLZRi+cBuWbj9p83uJ1IdhhIiISIHKjCY8tioNBSUGfLP/PIxmi9wlNRrDCBERkcJYLALPf7UfR/OK4eWqw9LJ4dBplLNg9UYMI0RERAqzKPk4kg7mwUEt4T9/HYBObZzkLumWMIwQEREpyLr0c/hgy28AgPjxfRAe4ClzRbeOYYSIiEghzhSUYu76AwCAx+8Iwv1hvjJX1DS4zwgREZFCBHi5YN7dIUg7exnPxwTLXU6T4bNpiIiIFEYIoYit3vlsGiIiIjtgMJnxxneZKCr7fXMzJQQRazCMEBER2SiLReC5L/fjo5TTmLZyt6I3NqsPwwgREZGNit90BN8eyIVGJeH5mGC7mxG5jmGEiIjIBq3YcRofpZwGAPzzgT6I7uYlc0XNh2GEiIjIxmw6mIvXvssEAMz5UzDu628ft/DWhWGEiIjIhvx6qhDPrN0HIYC/Du6Cvw8LkrukZsd9RoiIiGxIezcd2rvq0KuTO14d29tu14n8EcMIERGRDena3hXrH4+Ch5MDNOrWcQGjdfSSiIjIhuUVlSP1t4Kq997ujnB0UO5TeK3FMEJERCSjK2VGTFnxK6Z+vBvJmRfkLkcWDCNEREQyKS6vwNSP9+D4hRJ4umjRs6Ob3CXJgmGEiIhIBmVGEx5duQf7s6+grbMDPp0eAT9PZ7nLkgXDCBERUQsrrzDjsVVp2HPmMtwcNfh0egR6eLfOWRGAYYSIiKhFGU0WPL56L3b+VggXrRqfPDoIoZ095C5LVgwjRERELUitkuDpooVOo8LyaQMxoEtbuUuSnSQU8AhAvV4PDw8PFBUVwd3dXe5yiIiIbonFInAivwTBdr5gtaHf35wZISIiamYVZgtW7jwNk9kCAFCpJLsPItbgDqxERETNqMJswTNrMpB0MA8Hc/R4d0JfuUuyOQwjREREzaTCbMHTX2Rg06E8aNUqjOnTUe6SbBLDCBERUTMwmMx4+osM/HD4ArRqFZZMDsPwnh3kLssmMYwQERE1satGM2Z+lo7txy9Cq1FhyV8ZROrDMEJERNSEhBD4++rKIOLkoMZHU8IxpLuX3GXZNN5NQ0RE1IQkScL0IYHwctXi0+mDGEQaoFFhJCEhAYGBgXB0dERYWBhSUlLqbb969Wr07dsXzs7O8PHxwSOPPILCwsJGFUxERGSL/rht19Du7bF9znCEB3jKWJFyWB1G1q5di9mzZ2PevHnIyMjA0KFDMXr0aGRlZdXafseOHZgyZQqmT5+Ow4cP46uvvsKePXswY8aMWy6eiIjIFmQVluG+hFT8ll9SdcxZy5UQDWV1GFm0aBGmT5+OGTNmICQkBO+99x78/PyQmJhYa/tdu3YhICAATz/9NAIDAzFkyBDMnDkTaWlpt1w8ERGR3A6fL8L4xFTsy76Cf2w8KHc5imRVGDEajUhPT0dMTEy14zExMUhNTa31nKioKJw7dw5JSUkQQuDChQv4+uuvMWbMmDp/j8FggF6vr/YiIiKyNbtOFWLSkl0oKDEgxMcd/5rUX+6SFMmqMFJQUACz2Qxvb+9qx729vZGXl1frOVFRUVi9ejUmTpwIrVaLjh07ok2bNvj3v/9d5++Jj4+Hh4dH1cvPz8+aMomIiJrdD4fzMGXFbhQbTBgU6Ik1fxuMDu6OcpelSI1awCpJUrX3Qogax67LzMzE008/jZdffhnp6en4/vvvcfr0acyaNavOnx8XF4eioqKqV3Z2dmPKJCIiahYrd57GrM/SYTRZMKqXN1Y9OggeTg5yl6VYVq2u8fLyglqtrjELkp+fX2O25Lr4+HhER0fjhRdeAAD06dMHLi4uGDp0KF5//XX4+PjUOEen00Gn01lTGhERUYswmS1IOpgHIYAHB3XBa/f2hkbNnTJuhVX/72m1WoSFhSE5Obna8eTkZERFRdV6TllZGVSq6r9GrVYDqH4bFBERkRJo1CosnRKGN+4LxZv3hTKINAGr/x+MjY3FsmXLsGLFChw5cgTPPvsssrKyqi67xMXFYcqUKVXtx44di/Xr1yMxMRGnTp3Czp078fTTT2PQoEHo1KlT0/WEiIiomRSUGPDZrrNV79s4a/FwhH+dSxTIOlbfBD1x4kQUFhZiwYIFyM3NRWhoKJKSkuDv7w8AyM3NrbbnyLRp01BcXIwPPvgAzz33HNq0aYMRI0bg7bffbrpeEBERNZOjeXpMX5mGnCtXoZIkPBTRRe6S7I4kFHCtRK/Xw8PDA0VFRXB3d5e7HCIiaiV+zLyAZ9ZkoNRoRkA7ZyybOhDdOrjKXZZiNPT7m9vDERER3UAIgaXbT+Gt749CCCCyazsk/nUA2jhr5S7NLjGMEBER/UF5hRnzNhzCur3nAAAPRXTB/D/3hgMXqjYbhhEiIqI/2Jt1GeszzkElAS/f0wtTowK4ULWZMYwQERH9QVSQF+bdHYKeHd0xpLuX3OW0CpxzIiKiVk0IgU9/OYPsS2VVx2YM7cog0oIYRoiIqNUqMZjwzJp9eOm/hzHz03SUV5jlLqlV4mUaIiJqlY7m6fH46r04dbEUapWEv4T5Qqfhf6PLgWGEiIhana/SsvHSfw+hvMKCju6O+OCh/ggP8JS7rFaLYYSIiFqN8gozXtp4CF+lV962e3uP9lg8oS/aufLhrHJiGCEiolZDrZJwPL8EKgl49s4eeGJ4N6hUvG1XbgwjRERk1ywWAYsQ0KhVcFCr8K9J/XD+Sjkig9rJXRpdw5U6RERktwpKDHj0kz14N/l41TH/di4MIjaGYYSIiOzS5sN5uGvxdmw9dhErd55BfnG53CVRHXiZhoiI7EpxeQXm/y8TX19bpNqzoxvem9QPHdwcZa6M6sIwQkREdmPXqUI89+V+5Fy5CkkCZt4ehGdHdYdOo5a7NKoHwwgREdmF4vIKPLYqDcXlJvh5OmHRhH4YyL1DFIFhhIiI7IKbowPm3R2C/eeuYN6YXnDV8StOKbiAlYiIFElfXoG56w5g67H8qmOTBnVB/Pg+DCIKw9EiIiJFEUJg06E8vPrNYeQXG5ByogBbnr8DWj5XRrEYRoiISDGyL5XhlW8O4+ejlbMhgV4uePsvfRhEFI5hhIiIbF6F2YKPd57G4uQTuFphhoNawt/v6IbH7wiCowPvlFE6hhEiIrJ5e05fwptJRwEAgwI98eZ9t6FbB1eZq6KmwjBCREQ2yWAyV+0PEtXNCxPD/RDm3xYPhPtCkvhwO3vCMEJERDalvMKM5TtOY9UvZ/C/p4ZU7Zz69v19ZK6MmgvDCBER2QQhBP53IBdvbzqKnCtXAQBf7snGkyO6y1wZNTeGESIikl1G1mW89m0m9mZdAQD4eDji//7UE/f26yRvYdQiGEaIiEg2QgjEfrkfGzJyAADOWjX+PiwIM4Z2hZOWd8m0FgwjREQkG0mS4OmihSQB9w/wxfN3BcPbnU/XbW0YRoiIqMUUlVXgo5RTiOntjT6+bQAATw7vhgfCfdGzo7u8xZFsGEaIiKjZ6csr8PGOM1i24xSKy01IP3sZnz8WAUmS0NZFi7YuWrlLJBkxjBARUbMpMZjwSeoZLN1+CkVXKwAAPTu6YVp0gLyFkU1hGCEiIuttiQdUamDYnJqfbXsHsJix2vkh/POHY7hSVhlCgtq74NlRPXB3qA9UKm5aRr9jGCEiIuup1MCWNyr/+Y+BZNs7lceHz4NKknClrAJdvVzw5IhuuLdfZ6gZQqgWDCNERGS96wHkWiDJvu1JnPjqZYzI/QgYPg8YNgfjTWa4OzrgT6EdGUKoXgwjRETUOMPmIPVSHr49+B/c9ct7GFFeiuUOD2La0BegBqDTqDGmj4/cVZICMIwQEZFVTGYzEnd/h9VHP0Op6gjg5oICjQpReUb0nPg6OAlC1lI15qSEhAQEBgbC0dERYWFhSElJqbPttGnTIElSjVfv3r0bXTQREbW80opSfHH0Cwz5fDSWHp+HUtURSAIYXVKKp66UwgEViM5ZwSfqktWsnhlZu3YtZs+ejYSEBERHR2PJkiUYPXo0MjMz0aVLlxrt33//fbz11ltV700mE/r27YsHHnjg1ionIqIWse10Jn46twHJ2d+ipKKk8qDFESNER/zf+Z3oNHRu5RqS64tXgdrvsiGqgySEENacEBERgQEDBiAxMbHqWEhICMaNG4f4+Pibnr9x40aMHz8ep0+fhr+/f4N+p16vh4eHB4qKiuDuzh36iIiaW3mFER/8+g3WnfgKJarMquP+7v6YFDwJY86fRtsdC6sWq1b5w900DCTU0O9vq2ZGjEYj0tPTMXfu3GrHY2JikJqa2qCfsXz5ctx55531BhGDwQCDwVD1Xq/XW1MmERE1Utq5k1i861Mc0P8IqIsAFSCEhPbqflgw/DFEd46GSlIBF+JrDxzX31vMLV88KZZVYaSgoABmsxne3t7Vjnt7eyMvL++m5+fm5mLTpk34/PPP620XHx+P+fPnW1MaERE1ktFsxJbsLXhz+ycotByCJAlADcDsit5uI/F85DSE+3arftLwuLp/IGdEyEqNupvmxsVJQogGLVhauXIl2rRpg3HjxtXbLi4uDrGxsVXv9Xo9/Pz8GlMqERHVQgiBb4+m4ZD+RySdTsIVwxUAgCQBrpaeGNv1Pjwz+C9w0enkLZRaBavCiJeXF9RqdY1ZkPz8/BqzJTcSQmDFihWYPHkytNr6H4ik0+mg4x8AIqImdygvGx/s+RK/5ifDpMmpOt7BqQNG+o3ByM73IKJLDxkrpNbIqjCi1WoRFhaG5ORk3HfffVXHk5OTce+999Z77rZt2/Dbb79h+vTpjauUiIgaJbf4Mj74dSN+zt6EYulo5WUYDSAsavR0j8QzEQ8hslMkNCpuPUXysPrfvNjYWEyePBnh4eGIjIzE0qVLkZWVhVmzZgGovMSSk5ODVatWVTtv+fLliIiIQGhoaNNUTkREdbpquoqUcylYf/xb7MhJgaQyASpAAuBoDsKwTn/CM4MfgF+bdnKXSmR9GJk4cSIKCwuxYMEC5ObmIjQ0FElJSVV3x+Tm5iIrK6vaOUVFRVi3bh3ef//9pqmaiIhquFxWimXpSdhzcSvOXN2Dq6arAABJBahN3ujvORx/D5+AQX7dZa6UqDqr9xmRA/cZISKq3YWSIixLS8KPWcm4aN4PSWWs+qyza2fEBMQgyvtODOoUCpWqUZtuEzVas+wzQkRE8iu8WoiEX7/BD2d+whVxuPISDCpnQCRTW/Rwi0Zs5AOI9O3PrdlJERhGiIhsnMViwbbTh3G0eBdSz2/H/ov7ISAAqfJWXLWpPXq6R2FS73vw556DOANCisMwQkRkg4rKy7DmwFZsOrUFp0v3wKIprPZ5kHtPtEU/PBQ6BiOD+jCAkKIxjBAR2YgjF89gadp3SLuwC5fFYUiqisoPrt2G6+/SB1P6jMEwv2Ho6NJR3mKJmhDDCBGRTPJL9NievQsn9OlIPZ+KM/ozlR9cu/wCswf8HPtjlP9wTO53J7xcuICf7BPDCBFRCymvMOJ/R/cg6bftyLyyB6XSKUjS7w+UU0tquEvd0MN9IMb3HIk/dR/Ayy/UKjCMEBE1EyEEThWdwoe/bEJa/m5cNh8B1OWVH17bgExtbofxISMQ3Skag3wGwU3rJmvNRHJgGCEiaiIWiwU7zh7F5lM7YdScwO683Sgs/8PCUzUAsxPaqXthQIcIPNB7OAb7BfP2W2r1GEaIiBrJYrFg+5nD+Pb4TuzNT8NF0xFAra/WxlHtCF/nXvBS98bd3YdiTI+B0Gr4Vy/RH/FPBBFRA5ktZpy4cgLpF9KxITMFx4r2A+qS3xuoK+96cUEQxvYYitHdhuI2r9ugVdf/pHKi1o5hhIioDvryq/j22G78fHoXjlzZD5PDaZSZbgwfGrgiCD08+mJkQCTuDYlEGycX+YomUiCGESKia84XX8LXh1KQei4Np4oPoUw6U7XVOiQAJsDFwQX92vdDr7b94OXQE2N7RsBN5yRr3URKxzBCRK2SxWLBnnMnceTyfpwtO4x9+fvw25Xffm+grswfMLuinToYoe364d6QoRge2BcaFf/qJGpK/BNFRK1CWYUBm46n46dTvyLz0gEUmo8B6uIa7dTmDvDRhaB/h/4Y0yMKkX7B3OuDqJkxjBCRXSoyFGH/xf3Yk7sXXx5KQSlO/769OlC53kOo0UYViHEh0ejfoT/6deiHdk7t5CuaqJViGCEixbNYLNh97gS+O/4L9uZn4LL5GIotOb83UF2/5OKMNqruCG57G+7wH4R7ggdxsSmRDWAYISLFqbBUIOloOn44+Uu9l1z83f3Rr30/uEvdMdBnAIb694JGrZahYiKqD8MIEdm8iyV6JJ3YhRLpBDIuZOBAwQFcNV39vcG1/T2chD8CXEMxuFMYJvcfig4u7eUrmogajGGEiGzO6UsXsD5zB1Jz9uBs6SGUq7IhSZZqbXQqVziau6Jn2z4Y7h+Be3oOgoejs0wVE9GtYBghItkVXC1AWl4avjy0FRn5e2HS5P7+4bVbbCWTJ6J8wzEiIAL9O/RHUJsgqCTe5UJkDxhGiKjFHb14Dl8f3o5d53fDpD2BnNKs3z+89reS2tQRnZ16Y1DHcPy5ZzT6dwqUp1gianYMI0TU7M5cyseXh7ZiR84uZJUdgFlz4fcPKwAJEoI9g9GrbT+4ih4Y32sogtp1lK9gImpRDCNE1OSumq5i74W92JW7C1vO7sSZ4t8gSaLyQw0ghAStxReBLn1wd/do3N/7dnjoPOQtmohkwzBCRLfMZDbj+xPp+O+xrTh4aTeuqk7CAlPV55IEaEw+6OJ8G4b4RmLibcPQpQ3vdCGiSgwjRNQoJwvz8On+zdiZk4q8iv2A+trTbK+tKe3o0hGDfQZjsM9ghLQZgK6ePvIVS0Q2jWGEiBrEIiw4XHAYKTkpSDmXgkMFh4Hrl17UgLDo0FYKwYD2g/CXXiMwNKAXJEmSt2giUgSGESKqU36JHp9k/ICfsrYg17gPFtUfdjmVAAeTL3q4h2N00DDc33soXHQ6+YolIsViGCGiajIvZGNFxnf4JW8bisQRSCpz5QcqwEntjCG+0RjaeSgGd4yCj5u3vMUSkV1gGCEinNWfxY9nf8Saw0nIMxyvPChVLjxVmbwQ5DIIY7qNxKTbhnH2g4iaHMMIUSu15eQBfHLgf8g370Z26alqn+nMgejnGY2H+9yNYQG9oVJxp1Miaj4MI0StSMqZTCzP2ID9l7ZV23JdI2kwsONA3OE3An09o9Hb20/GKomotWEYIbJzWfrzeHPbF9hz8ScY1dmVBzWAEGq0QS/EBNyFZyLHcdMxIpINwwiRHbp8tQhbzv2Eb099i7S8NAiIyttvhQru6IXhnUdh1sB74demndylEhExjBDZC6PJhOXpP+Dr4xuRb04HpIqqz/yceyPEbRgeHzgeQe14BwwR2RaGESKF23HmCD7Y8zkyi7dAqIsqD0qAj5M/JoSMw92Bd6OTayd5iyQiqgfDCJEClZvKsWzvRnx6+EuUqU5UHlQDMDsjyGkoHun7AMb2HMi7YIhIERr1N1VCQgICAwPh6OiIsLAwpKSk1NveYDBg3rx58Pf3h06nQ1BQEFasWNGogolas+OXfsPbu9/GiK9GYMmRN1CmOgEhJHiI2/DXwJfwy8PbsPHBhbi3VwSDCBEphtUzI2vXrsXs2bORkJCA6OhoLFmyBKNHj0ZmZia6dOlS6zkTJkzAhQsXsHz5cnTr1g35+fkwmUy1tiWi6q5WGPCvXRux4bcvUao6XnW8k0snBDmNwKywiejTMUC+AomIbpEkhBDWnBAREYEBAwYgMTGx6lhISAjGjRuH+Pj4Gu2///57TJo0CadOnYKnp2ejitTr9fDw8EBRURHc3d0b9TOIlOZEQS7eSPkYey8nVa0FEULCMN9hmNRzAqI6RUGtUstcJRFR3Rr6/W3VzIjRaER6ejrmzp1b7XhMTAxSU1NrPeebb75BeHg43nnnHXz66adwcXHBn//8Z7z22mtwcnKq9RyDwQCDwVCtM0StxQ8n0rF493KcM6ZWPhdGDcDshtvcY/BC1FT07xQod4lERE3KqjBSUFAAs9kMb+/qtwZ6e3sjLy+v1nNOnTqFHTt2wNHRERs2bEBBQQEef/xxXLp0qc51I/Hx8Zg/f741pREpmkVYsCNnBz45/Al25+0GAEgqQGcOwNiACXgu+gG46hxlrpKIqHk06m4aSZKqvRdC1Dh2ncVigSRJWL16NTw8Knd4XLRoEe6//358+OGHtc6OxMXFITY2tuq9Xq+Hnx+3pyb7U2Y04K2Uz/HrpfU4X3YGAKCS1GgnhWNWv2l4IDS6zj9bRET2wqow4uXlBbVaXWMWJD8/v8ZsyXU+Pj7o3LlzVRABKteYCCFw7tw5dO/evcY5Op0OOj4ZlOxYYVkxXvppGXZcXFe1HsTFwQX3d78fD4c8DB9XH5krJCJqOVbd+6fVahEWFobk5ORqx5OTkxEVFVXrOdHR0Th//jxKSkqqjh0/fhwqlQq+vr6NKJlIuc7rL+HRjfG4Y80opFxaURlEzO4Y3GYKvh+/Gc8PfJ5BhIhaHasv08TGxmLy5MkIDw9HZGQkli5diqysLMyaNQtA5SWWnJwcrFq1CgDw0EMP4bXXXsMjjzyC+fPno6CgAC+88AIeffTROhewEtmbIkMR/i/5A+y8uBFQlQNqQGXywl2+D+IfwybD3ZF/Foio9bI6jEycOBGFhYVYsGABcnNzERoaiqSkJPj7+wMAcnNzkZWVVdXe1dUVycnJeOqppxAeHo527dphwoQJeP3115uuF0Q2Sm/UY9XhVVh9ZDVKKkoAFaA2+WB84BTMGToBjg5auUskIpKd1fuMyIH7jJDSXCgpQtyPH+JA8TcwWEoBAD3a9kCk5yTMjhwPjZr7gxCR/WuWfUaIqH56Qxnikpdie/4aQF0ZQoI8uuGJ/o9jZJeRUEncop2I6EYMI0RNwGgy4fVtn2LjmY8hNJevrQlpj790nY642yfBgTMhRER1YhghukUrM37A+3sXwaQ5X/knyuSBu32nYv6IqVwTQkTUAAwjRI104vIJLExbiNTzqZV/ksxOGNzuAbwz6nG0dXaRuzwiIsVgGCGy0qlLF/DWrsX4tWATLMICB5UDBnqOxYtRT8K/bXu5yyMiUhyGEaIGulphwAs/JGJb/mpAXQ4AGOU/Cs+GPQs/Nz6ugIiosRhGiBrg4/TN+Nf+f8KkzgPUgIPZDy9GzMX9vW+XuzQiIsVjGCGqx6G8LDyTPB/5lt2AGoDZBaM7P4rXRz4KrYZ/fIiImgL/NiWqhcliwqpDq7E4/d+AygAhJHTVjcKHo+Pg18ZL7vKIiOwKwwjRDQ4XHMb8X+bjyKUjgArQmQPxatTLuKdnuNylERHZJYYRomsulugx67s3cKL8ewhY4KZ1wzP9n8VfunP7diKi5sQwQgTgoz2b8O+Db0GoLwEARgfcjTmDXoCXEy/JEBE1N4YRatVy9Jfw2LevILtiK6AGJFNbzAydgyci7pG7NCKiVoNhhFqtD3f9D0sy34ZQF0EICUG6GHz0l5fRwZVPhiYiakkMI9TqlBhL8PxPr2FnflLVA+2eH/ASJvcfLndpREStEsMItSp78vbgHzv+gfOl5wEhobvjaCx/4CW0dXaVuzQiolaLYYRahWLDVTye9Cb26/8LAYHOrp0xP/I1RHQaKHdpREStHsMI2b2tpw8hdssLqFCfAwCM7z4ecwbOgYsDn6xLRGQLGEbIblksFrz44zJ8m/MfSOoKwOyMqT3m4Pmov8hdGhER/QHDCNmlc0WXMPm/z6FApEFSAa6WXlh+z7vo1cFX7tKIiOgGDCNkd7458gv+kToXQnMJQqgR1XYyEu6ZzV1UiYhsFMMIKc+WeEClBobNqXZYCIFP//cIFl3eC6ERUJna4ZXBb2F878EyFUpERA3BMELKo1IDW96o/OdrgaSw7Ape++4h/FSWDQCI6DAcC6JfRSd3T7mqJCKiBmIYIeW5PiNyLZB873MnXk+ZiSJNOTRQYU7EXEwKngRJkmQskoiIGophhJTpWiD5Zs97WNDuExg0KribdEgc+zH6dLhN5uKIiMgaDCOkSCazBS9dicAPXp6okCRElpXjxYk/IcCzg9ylERGRlVRyF0BkrcISA6as2A2v3aswt/AyZl4pRuKFfAQcXCl3aURE1AicGSFFuVJmxJ8/2InxxavxnMPXOOH3NLpPeA3Y9k6NRa1ERKQMDCOkKG2ctXjdMwnDy79G4cDn0X3MS5Uf3LColYGEiEg5GEbI5gkhcLXCDGdt5b+uw7p5ojxoLtqNjKve8HoAsZhbuEIiIroVkhBCyF3Ezej1enh4eKCoqAju7u5yl0MtqLzCjLnrDiC3qByfTo+AVsNlTkREStHQ72/OjJDNKiwxYMaqNGRkXYFaJSH97GVEBrWTuywiImpiDCNkk05dLMEjK/fgbGEZPJwckPjwAAYRIiI7xTBCNif97CXM+CQNl8sq4OfphI+nDUK3Dq5yl0VERM2EYYRsSnLmBTzx+V4YTRb09fXAsqkD0d5NJ3dZRETUjBq1GjAhIQGBgYFwdHREWFgYUlJS6my7detWSJJU43X06NFGF032y7+dM3QaFe4M8cYXfxvMIEJE1ApYPTOydu1azJ49GwkJCYiOjsaSJUswevRoZGZmokuXLnWed+zYsWoradu3b9+4ismu9fB2w4bHoxHo5QK1ig+6IyJqDayeGVm0aBGmT5+OGTNmICQkBO+99x78/PyQmJhY73kdOnRAx44dq15qtbrRRZP9MJkteGnjIew6VVh1rFsHVwYRIqJWxKowYjQakZ6ejpiYmGrHY2JikJqaWu+5/fv3h4+PD0aOHIktW7bU29ZgMECv11d7kf0przDj8dV78emus5j1WTr05RVyl0RERDKwKowUFBTAbDbD29u72nFvb2/k5eXVeo6Pjw+WLl2KdevWYf369QgODsbIkSOxffv2On9PfHw8PDw8ql5+fn7WlEkKUGIwYeqK3diceQFajQpvje8Dd0cHucsiIiIZNOpuGkmqPoUuhKhx7Lrg4GAEBwdXvY+MjER2djYWLlyI22+/vdZz4uLiEBsbW/Ver9czkNiRK2VGTP14D/ZnX4GrToOPpoRzDxEiolbMqjDi5eUFtVpdYxYkPz+/xmxJfQYPHozPPvuszs91Oh10Ot5FYY8KSgz467JfcTSvGG2cHfDpoxG4zddD7rKIiEhGVl2m0Wq1CAsLQ3JycrXjycnJiIqKavDPycjIgI+PjzW/muxE4taTOJpXDC9XHdb+LZJBhIiIrL9MExsbi8mTJyM8PByRkZFYunQpsrKyMGvWLACVl1hycnKwatUqAMB7772HgIAA9O7dG0ajEZ999hnWrVuHdevWNW1PSBHm/CkYxeUVmDUsCF3bc1dVIiJqRBiZOHEiCgsLsWDBAuTm5iI0NBRJSUnw9/cHAOTm5iIrK6uqvdFoxPPPP4+cnBw4OTmhd+/e+O6773D33Xc3XS/IphWWGODpooUkSdBp1Hjn/r5yl0RERDZEEkIIuYu4mYY+gphsz9nCUjy4dBdienfEK2N71bnQmYiI7E9Dv78btR08UUNcDyLni8qx/cRF6MtNcpdEREQ2iGGEmsXZwlJMuhZEgtq7YM1jg+HhxH1EiIioJj61l5pc9qUyTFq6C7nXgsgXfxuMDm6OcpdFREQ2ijMj1KRyi67ioWUMIkRE1HAMI9Sk9mdfQc7lq/Bv54zPH2MQISKim+NlGmpSfwr1QeJfwxDa2QPe7gwiRER0cwwjdMuKrlbAYDJXzYLc1bujzBUREZGS8DIN3ZKrRjOmr9yDiUt24fyVq3KXQ0RECsQwQo1WYbbg8dXpSDt7GQUlBlwpq5C7JCIiUiCGEWoUi0Xgha/2Y8uxi3B0UGHFtIHo1Ym74xIRkfUYRshqQggs+DYTG/edh0YlIfHhMAwM8JS7LCIiUiiGEbLav3/+DStTzwAAFj7QF8N7dpC3ICIiUjSGEbJKicGEdXvPAQBeHdsL4/p3lrkiIiJSOt7aS1Zx1Wmw7u9RSM68gAcHdZG7HCIisgOcGaEGuWo0V/2zl6uOQYSIiJoMwwjdVFZhGYYv3Iqv0rLlLoWIiOwQwwjV61KpEVM/3o08fTlWpp5Bhdkid0lERGRnGEaoTuUVZjy2Kg2nC0rRuY0TPp42EA5q/itDRERNi98sVCshBOZ8fQDpZy/DzVGDlY8MRAc++I6IiJoBwwjV6v2fTuCb/ZWbmv3nr2Ho7u0md0lERGSnGEaohrQzl/DejycAAK+NC0V0Ny+ZKyIiInvGfUaohjD/togd1QOlBhNv4SUiombHMEI1SJKEp0d2hxBC7lKIiKgV4GUaAgCUGU14+/ujKDOaqo5JkiRjRURE1FowjBCEEHj+q/1I3HoSsz7bK3c5RETUyjCMED74+TckHcyDg1rC0yO6yV0OERG1Mgwjrdzmw3l4N/k4AGDBvaEID/CUuSIiImptGEZaseMXivHs2n0AgCmR/rxzhoiIZMEw0koVXa3A31alodRoxuCunnjpnl5yl0RERK0Uw0grlX2pDMXlJnRu44SEh8P4zBkiIpIN9xlppUI7e+Dbp4fgSlkFPF20cpdDREStGMNIK1NhtlTNgvh4OMHHw0nmioiIqLXj3Hwrkn2pDMMXbsX3h3LlLoWIiKgKw0grUV5hxt9Xp+Pc5atI3HoSZgu3eiciItvAMNJKzP9fJg7l6OHpokXCX8OgVnGrdyIisg0MI63Af/fl4IvdWZAk4P1J/dC5DdeJEBGR7WhUGElISEBgYCAcHR0RFhaGlJSUBp23c+dOaDQa9OvXrzG/lhrh5MUSxK0/CAB4akR3DO3eXuaKiIiIqrM6jKxduxazZ8/GvHnzkJGRgaFDh2L06NHIysqq97yioiJMmTIFI0eObHSxZJ2rRjOeWL0XZUYzIru2wzMju8tdEhERUQ1Wh5FFixZh+vTpmDFjBkJCQvDee+/Bz88PiYmJ9Z43c+ZMPPTQQ4iMjGx0sWQdtUpCVJAX2rvp8P6D/bhOhIiIbJJVYcRoNCI9PR0xMTHVjsfExCA1NbXO8z7++GOcPHkSr7zySoN+j8FggF6vr/Yi62k1Krw8theSn70dHdwc5S6HiIioVlaFkYKCApjNZnh7e1c77u3tjby8vFrPOXHiBObOnYvVq1dDo2nYHmvx8fHw8PCoevn5+VlTZqt3sdgAk9lS9b6NM3dYJSIi29WoBaySVH26XwhR4xgAmM1mPPTQQ5g/fz569OjR4J8fFxeHoqKiqld2dnZjymyVDCYzpn28G5OW7sL5K1flLoeIiOimrNoO3svLC2q1usYsSH5+fo3ZEgAoLi5GWloaMjIy8OSTTwIALBYLhBDQaDTYvHkzRowYUeM8nU4HnU5nTWl0zT+/P4bD5/Vo6+wAVS0BkYiIyNZYNTOi1WoRFhaG5OTkaseTk5MRFRVVo727uzsOHjyIffv2Vb1mzZqF4OBg7Nu3DxEREbdWPVWz9Vg+lu04DQD45/190dGD60SIiMj2Wf2gvNjYWEyePBnh4eGIjIzE0qVLkZWVhVmzZgGovMSSk5ODVatWQaVSITQ0tNr5HTp0gKOjY43jdGsuFhvw/Ff7AQBTIv1xZ6+aM1VERES2yOowMnHiRBQWFmLBggXIzc1FaGgokpKS4O/vDwDIzc296Z4j1LQsFoHnvtqPghIjenZ0w4t3h8hdEhERUYNJQgibf2KaXq+Hh4cHioqK4O7uLnc5Nmf5jtN47dtM6DQqfPvUEHT3dpO7JCIiogZ/f1s9M0K2Z0g3L4T4uOOhiC4MIkREpDgMI3YguKMb/vtENBzUvHuGiIiUh0/tVbC8ovKqf9ZqVLXu9UJERGTrGEYUatepQgx952e89+NxKGDZDxERUZ0YRhRIX16B577cjwqzwPkrVzkjQkREisYwokCvfnMYOVeuoounM14e21vucoiIiG4Jw4jC/HA4D+v35kAlAYsm9IWrjmuQiYhI2RhGFORSqRHzNhwEAPzt9iCEB3jKXBEREdGtYxhRkFe+OYyCEiO6d3DFs6O6y10OERFRk2AYUQghBMK6tIGrToOFD/SFTqOWuyQiIqImwQUHCiFJEqZFB2J8mC/cHR3kLoeIiKjJcGbExgkhUF5hrnrPIEJERPaGYcTGfXcwFzGLt+PXU4Vyl0JERNQsGEZs2KVSI17572FkXSrDzpMMI0REZJ8YRmzYa99morDUiGBvNzw5vJvc5RARETULhhEbteVYPjZkVG5u9vb9faDVcKiIiMg+8RvOBpUYTJi3vnJzs0ejA9HPr428BRERETUjhhEb9M73R3G+qBx+nk6IjekhdzlERETNimHExpgtlU/iBYC3xveBs5ZbwRARkX3jN52NUaskfDQlHBnZVzCgS1u5yyEiImp2nBmxQZIkMYgQEVGrwTBiI37LL0Hc+gO4XGqUuxQiIqIWxcs0NkAIgXkbDuLX05dQZjTj/Un95S6JiIioxXBmxAZ8nX4Ov56+BEcHFZ6PCZa7HCIiohbFMCKzS6VGvJl0BADw7J094OfpLHNFRERELYthRGZvJh3B5bIK9OzohkeHBMpdDhERUYtjGJHRLycL8XX6OUgS8MZ9t8FBzeEgIqLWh99+MlqUfAwA8HBEF4T581ZeIiJqnRhGZPTRlHDMGBKIF+7qKXcpREREsuGtvTJq46zFP+7pJXcZREREsuLMiAzSz16CEELuMoiIiGwCw0gL2378Iv6S+AseWbkHZgsDCREREcNICzKYzHj1m8MAgK5erlCrJJkrIiIikh/DSAtalnIapwpK0d5Nh9mjustdDhERkU1gGGkhOVeu4t8/nwAAzLs7BO6ODjJXREREZBsYRlrIm98dQXmFBYMCPXFvv05yl0NERGQzGhVGEhISEBgYCEdHR4SFhSElJaXOtjt27EB0dDTatWsHJycn9OzZE4sXL250wUq061QhvjuYC5UEzP9zb0gS14oQERFdZ/U+I2vXrsXs2bORkJCA6OhoLFmyBKNHj0ZmZia6dOlSo72LiwuefPJJ9OnTBy4uLtixYwdmzpwJFxcX/O1vf2uSTtg6lSSha3sXRAW1Q4iPu9zlEBER2RRJWLnhRUREBAYMGIDExMSqYyEhIRg3bhzi4+Mb9DPGjx8PFxcXfPrppw1qr9fr4eHhgaKiIri7K/PLvMJsgcFkgauO+8wREVHr0NDvb6su0xiNRqSnpyMmJqba8ZiYGKSmpjboZ2RkZCA1NRXDhg2rs43BYIBer6/2UjoHtYpBhIiIqBZWhZGCggKYzWZ4e3tXO+7t7Y28vLx6z/X19YVOp0N4eDieeOIJzJgxo8628fHx8PDwqHr5+flZU6bNWPjDMXy0/RSMJovcpRAREdmsRv2n+o0LMIUQN12UmZKSgpKSEuzatQtz585Ft27d8OCDD9baNi4uDrGxsVXv9Xq94gLJsbxiJG47CbNF4DZfDwzu2k7ukoiIiGySVWHEy8sLarW6xixIfn5+jdmSGwUGBgIAbrvtNly4cAGvvvpqnWFEp9NBp9NZU5pNEULgtW8zYbYI3NXbm0GEiIioHlZdptFqtQgLC0NycnK148nJyYiKimrwzxFCwGAwWPOrFWXLsXzs+K0AWrUK8+7mU3mJiIjqY/VlmtjYWEyePBnh4eGIjIzE0qVLkZWVhVmzZgGovMSSk5ODVatWAQA+/PBDdOnSBT179gRQue/IwoUL8dRTTzVhN2xHhdmCN747AgB4ZEgAurRzlrkiIiIi22Z1GJk4cSIKCwuxYMEC5ObmIjQ0FElJSfD39wcA5ObmIisrq6q9xWJBXFwcTp8+DY1Gg6CgILz11luYOXNm0/XChqzZnYWTF0vh6aLFE8O7yV0OERGRzbN6nxE5KGWfEaPJgui3f8bFYgNeu7c3JkcGyF0SERGRbJplnxGqn1ajwuoZEXg4ogseHFRzN1oiIiKqibtwNbEe3m54477b5C6DiIhIMTgz0kRyi67KXQIREZEiMYw0gfSzlzDk7S14aeMhKGAJDhERkU1hGLlFQgi8mXQUZouA0WS56U60REREVB3DyC3anHkB6Wcvw9FBhdiYHnKXQ0REpDgMI7fAZLbgne+PAgCmDwmEt7ujzBUREREpD8PILfgy7RxOXixFW2cHzBwWJHc5REREisQw0khlRhMW/3gcAPDUiO5wd3SQuSIiIiJlYhhppCO5epRXmOHn6YSHB3ODMyIiosbipmeNFObviZQ5w5F96Sp0GrXc5RARESkWw8gtaOOsRRtnrdxlEBERKRov01jp3OUybDmWz83NiIiImgjDiJUWJ5/AIx/vQfymo3KXQkREZBcYRqxw4kIxNmScAwCMuc1H5mqIiIjsA8OIFd7dfBwWAdzV2xt9/drIXQ4REZFdYBhpoP3ZV/D94TxIEvB8TLDc5RAREdkNhpEGWrj5GADgvv6d0d3bTeZqiIiI7AfDSAOknixAyokCOKglPHsnH4ZHRETUlLjPSAMFtXdBVJAX/Dyd5S6FiIjIrjCMNEBUkBc2PzsM5RVmuUshIiKyO7xM00BqlQQXHbMbERFRU2MYqcfWY/lYseM0Z0SIiIiaEf9Tvw4Wi8Bbm47iaF4xyowmPDmiu9wlERER2SXOjNTh+8N5OJpXDDedBpMHB8hdDhERkd1iGKmF2SKwOPk4AODRIYHwcHaQuSIiIiL7xTBSi+8O5uJEfgncHTV4dEig3OUQERHZNYaRG5gtAu/9WDkr8tjQrvBw4qwIERFRc2IYucE3+3Nw6mIp2jg7YFp0gNzlEBER2T3eTXODEB93jOrljf5d2sDNkbMiREREzY1h5AY9O7rjoynhEELIXQoREVGrwMs0dZAkSe4SiIiIWgWGkWuSDubiHxsPIufKVblLISIialV4mQaVu60uTj6OE/klaO/qiGfu5G6rRERELYUzI6jcbfVEfgncHDW8g4aIiKiFNSqMJCQkIDAwEI6OjggLC0NKSkqdbdevX49Ro0ahffv2cHd3R2RkJH744YdGF3zLtsQD296pemuxCPzrpxMAgP/4/QSPXQvlqoyIiKhVsjqMrF27FrNnz8a8efOQkZGBoUOHYvTo0cjKyqq1/fbt2zFq1CgkJSUhPT0dw4cPx9ixY5GRkXHLxTeKSg1seaMqkPx45AKO5hXjOe1GRGcvqfyciIiIWowkrLyHNSIiAgMGDEBiYmLVsZCQEIwbNw7x8fEN+hm9e/fGxIkT8fLLLzeovV6vh4eHB4qKiuDu7m5NubXb9g6w5Q2IO17E2INRGJ63Es85fA0MnwcMm3PrP5+IiIga/P1t1QJWo9GI9PR0zJ07t9rxmJgYpKamNuhnWCwWFBcXw9PTs842BoMBBoOh6r1er7emzJu7FjikLW9gndBA52BCWfRcODOIEBERtTirLtMUFBTAbDbD29u72nFvb2/k5eU16Ge8++67KC0txYQJE+psEx8fDw8Pj6qXn5+fNWU2zLA5EGotdJIJJskBzqPimv53EBER0U01agHrjRuCCSEatEnYF198gVdffRVr165Fhw4d6mwXFxeHoqKiqld2dnZjyqzftncgmY2AWguNqKi2qJWIiIhajlWXaby8vKBWq2vMguTn59eYLbnR2rVrMX36dHz11Ve48847622r0+mg0+msKc0619aMVK0Ruf4e4JoRIiKiFmbVzIhWq0VYWBiSk5OrHU9OTkZUVFSd533xxReYNm0aPv/8c4wZM6ZxlTaVG4MIUPm/w+dVu8uGiIiIWobVO7DGxsZi8uTJCA8PR2RkJJYuXYqsrCzMmjULQOUllpycHKxatQpAZRCZMmUK3n//fQwePLhqVsXJyQkeHh5N2JUGsphrv2vm+nuLueVrIiIiasWsvrUXqNz07J133kFubi5CQ0OxePFi3H777QCAadOm4cyZM9i6dSsA4I477sC2bdtq/IypU6di5cqVDfp9TX5rLxERETW7hn5/NyqMtDSGESIiIuVp6Pc3n01DREREsmIYISIiIlkxjBAREZGsGEaIiIhIVgwjREREJCuGESIiIpIVwwgRERHJimGEiIiIZMUwQkRERLKy+tk0cri+Saxer5e5EiIiImqo69/bN9vsXRFhpLi4GADg5+cncyVERERkreLi4nofjquIZ9NYLBacP38ebm5ukCSpyX6uXq+Hn58fsrOz7faZN/beR/ZP+ey9j/beP8D++8j+NZ4QAsXFxejUqRNUqrpXhihiZkSlUsHX17fZfr67u7td/gv2R/beR/ZP+ey9j/beP8D++8j+NU59MyLXcQErERERyYphhIiIiGTVqsOITqfDK6+8Ap1OJ3cpzcbe+8j+KZ+999He+wfYfx/Zv+aniAWsREREZL9a9cwIERERyY9hhIiIiGTFMEJERESyYhghIiIiWdl9GElISEBgYCAcHR0RFhaGlJSUettv27YNYWFhcHR0RNeuXfGf//ynhSptPGv6uHXrVkiSVON19OjRFqy44bZv346xY8eiU6dOkCQJGzduvOk5ShpDa/untPGLj4/HwIED4ebmhg4dOmDcuHE4duzYTc9Tyhg2pn9KG8PExET06dOnakOsyMhIbNq0qd5zlDJ+gPX9U9r43Sg+Ph6SJGH27Nn1tmvpMbTrMLJ27VrMnj0b8+bNQ0ZGBoYOHYrRo0cjKyur1vanT5/G3XffjaFDhyIjIwMvvvginn76aaxbt66FK284a/t43bFjx5Cbm1v16t69ewtVbJ3S0lL07dsXH3zwQYPaK20Mre3fdUoZv23btuGJJ57Arl27kJycDJPJhJiYGJSWltZ5jpLGsDH9u04pY+jr64u33noLaWlpSEtLw4gRI3Dvvffi8OHDtbZX0vgB1vfvOqWM3x/t2bMHS5cuRZ8+feptJ8sYCjs2aNAgMWvWrGrHevbsKebOnVtr+zlz5oiePXtWOzZz5kwxePDgZqvxVlnbxy1btggA4vLlyy1QXdMCIDZs2FBvGyWO4XUN6Z+Sx08IIfLz8wUAsW3btjrbKHkMG9I/pY+hEEK0bdtWLFu2rNbPlDx+19XXP6WOX3FxsejevbtITk4Ww4YNE88880ydbeUYQ7udGTEajUhPT0dMTEy14zExMUhNTa31nF9++aVG+7vuugtpaWmoqKhotlobqzF9vK5///7w8fHByJEjsWXLluYss0UpbQwbS6njV1RUBADw9PSss42Sx7Ah/btOiWNoNpuxZs0alJaWIjIystY2Sh6/hvTvOqWN3xNPPIExY8bgzjvvvGlbOcbQbsNIQUEBzGYzvL29qx339vZGXl5erefk5eXV2t5kMqGgoKDZam2sxvTRx8cHS5cuxbp167B+/XoEBwdj5MiR2L59e0uU3OyUNobWUvL4CSEQGxuLIUOGIDQ0tM52Sh3DhvZPiWN48OBBuLq6QqfTYdasWdiwYQN69epVa1sljp81/VPi+K1ZswZ79+5FfHx8g9rLMYaKeGrvrZAkqdp7IUSNYzdrX9txW2JNH4ODgxEcHFz1PjIyEtnZ2Vi4cCFuv/32Zq2zpShxDBtKyeP35JNP4sCBA9ixY8dN2ypxDBvaPyWOYXBwMPbt24crV65g3bp1mDp1KrZt21bnF7bSxs+a/ilt/LKzs/HMM89g8+bNcHR0bPB5LT2Gdjsz4uXlBbVaXWOGID8/v0biu65jx461ttdoNGjXrl2z1dpYjeljbQYPHowTJ040dXmyUNoYNgUljN9TTz2Fb775Blu2bIGvr2+9bZU4htb0rza2PoZarRbdunVDeHg44uPj0bdvX7z//vu1tlXi+FnTv9rY8vilp6cjPz8fYWFh0Gg00Gg02LZtG/71r39Bo9HAbDbXOEeOMbTbMKLVahEWFobk5ORqx5OTkxEVFVXrOZGRkTXab968GeHh4XBwcGi2WhurMX2sTUZGBnx8fJq6PFkobQybgi2PnxACTz75JNavX4+ff/4ZgYGBNz1HSWPYmP7VxpbHsDZCCBgMhlo/U9L41aW+/tXGlsdv5MiROHjwIPbt21f1Cg8Px8MPP4x9+/ZBrVbXOEeWMWy2pbE2YM2aNcLBwUEsX75cZGZmitmzZwsXFxdx5swZIYQQc+fOFZMnT65qf+rUKeHs7CyeffZZkZmZKZYvXy4cHBzE119/LVcXbsraPi5evFhs2LBBHD9+XBw6dEjMnTtXABDr1q2Tqwv1Ki4uFhkZGSIjI0MAEIsWLRIZGRni7NmzQgjlj6G1/VPa+P39738XHh4eYuvWrSI3N7fqVVZWVtVGyWPYmP4pbQzj4uLE9u3bxenTp8WBAwfEiy++KFQqldi8ebMQQtnjJ4T1/VPa+NXmxrtpbGEM7TqMCCHEhx9+KPz9/YVWqxUDBgyodsvd1KlTxbBhw6q137p1q+jfv7/QarUiICBAJCYmtnDF1rOmj2+//bYICgoSjo6Oom3btmLIkCHiu+++k6Hqhrl+G92Nr6lTpwohlD+G1vZPaeNXW98AiI8//riqjZLHsDH9U9oYPvroo1V/v7Rv316MHDmy6otaCGWPnxDW909p41ebG8OILYyhJMS1VSlEREREMrDbNSNERESkDAwjREREJCuGESIiIpIVwwgRERHJimGEiIiIZMUwQkRERLJiGCEiIiJZMYwQERGRrBhGiIiISFYMI0RERCQrhhEiIiKSFcMIERERyer/Ac/vgpNCVmupAAAAAElFTkSuQmCC",
"text/plain": [
"<Figure size 640x480 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Generate y from a cubic (to illustrate equivalence to cubic polynomial interpolation)\n",
"def cubic_poly(x):\n",
" x = numpy.array(x)\n",
" return 0.2 + 0.5*x - 0.2*(x**2) + 0.03*(x**3)\n",
"y = cubic_poly(x)\n",
"plt.plot(numpy.linspace(x[0],x[3],100), cubic_poly( numpy.linspace(x[0],x[3],100) ), '--', label='cubic');\n",
"plt.plot(x, y, 'x', label='data');\n",
"plt.plot(q, lagrange_interp(x,y,q), label='Lagrange interpolation');"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAiIAAAGdCAYAAAAvwBgXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAw0UlEQVR4nO3de5yWY+LH8c/MVFOqGSq1pYkklhKrWqZ1iJSNjRxDWix27SaH7C7ht8vKjmXXMSKsc5tDSk6RdECiYlYOq9BuUUrSTMep5rl/f9wVKekwz1zP4fN+vZ5Xz3X31P3d52Wb7+u+r+u6c6IoipAkSQogN3QASZKUvSwikiQpGIuIJEkKxiIiSZKCsYhIkqRgLCKSJCkYi4gkSQrGIiJJkoKpETrA5iQSCebOnUv9+vXJyckJHUeSJG2BKIpYsmQJzZo1Izd389c8UrqIzJ07l6KiotAxJEnSNpgzZw7Nmzff7GdSuojUr18fiP+HFBQUBE4jSZK2RHl5OUVFRet/jm9OSheRdbdjCgoKLCKSJKWZLZlW4WRVSZIUjEVEkiQFYxGRJEnBWEQkSVIwFhFJkhSMRUSSJAVjEZEkScFYRCRJUjAWEUmSFEy1FZGSkhJycnK4+OKLq+uUkiTpO24eM4Pbxs7c5O/dNnYmN4+ZUa15qqWITJkyhSFDhtCuXbvqOJ0kSfoeebk53LSJMnLb2JncNGYGebnV+7T7pBeRpUuX0rt3b+655x522mmnZJ9OkiRtxoVdWtO/654blJF1JaR/1z25sEvras2T9Ife9e3bl2OOOYYjjzySgQMHbvazFRUVVFRUrB+Xl5cnO54kSVnnwi6tqVvxJQXj+tN3XAeeW9MhSAmBJBeRYcOG8fbbbzNlypQt+nxJSQnXXHNNMiNJkpTdVi2HNwZxzju3QI1ldEx8xNi8DkFKCCTx1sycOXO46KKLeOSRR6hdu/YW/ZkBAwZQVla2/jVnzpxkxZMkKbtEEUx/EgZ1hHHXweplTEu05veJfqys5HsnsCZb0q6ITJs2jQULFtC+ffv1xyorK5k4cSKDBg2ioqKCvLy8Df5Mfn4++fn5yYokSVJ2mlsKz/8BPnsLgCX5Tbhiycm0PuJMnjxyz/VzRIDMmSPSpUsXpk+fvsGxs88+mx//+MdcdtllG5UQSZJUxZYvgrF/gWkPABHU3IE3mv2Ssz46iL5d264vHet+DVFGklZE6tevT9u2bTc4VrduXRo2bLjRcUmSVIUSlTDtfhh7LaxcHB9rexJ0u5bJby6lb4ucjcrGunFlIqrWqElfNSNJkqrR52/Ds5fAvNJ43LgNHH0D7HYwAJd0/f4/mnGrZr5r/Pjx1Xk6SZKyx4rF8Mq1MOU+IIL8QjjiSuhwDuSl7nWH1E0mSZJ+2LrVMC9eAcsWxMfa9YKu10L9JmGzbQGLiCRJ6WrRLHiuP3zySjxu2BqO+QfsfljYXFvBIiJJUrqpXAOT74BxJbBmBeTlw6F/gJ9dCDXSaxsMi4gkSelkbimMugC+WLtFxm6HwC9ugUZ7hEy1zSwikiSlg9UrYcLf4PVbIaqE2jtCt4HwkzMgp3qfmFuVLCKSJKW62W/GV0EWxhuO0eZ46H4D1GscNlcVsIhIkpSqVi2Pd0Z98y4ggnpN4JibYO9fhE5WZSwikiSlotlvwsjzYdGn8Xj/3nDUdVBnp7C5qphFRJKkVLJ6JYwbCJMGARHUbwbH3g6tjwydLCksIpIkpYrPp8GI87+ZC7J/bzjqr1Bnx6CxkskiIklSaJWrYeLfYeKN8YqYek2gx62wV/fQyZLOIiJJUkgLZ8JTv4a5b8fjNifEu6Pu0CBsrmpiEZEkKYQogin3wkv/F++OWrswXhGz70mhk1Uri4gkSdVt6QJ4ui/MfCke794ZjrsTCncJGisEi4gkSdVpxosw8newfGH8jJiu18BPfwO5uaGTBWERkSSpOqxeAWP+BG8NiceN28CJ90KTfcLmCswiIklSsi34EJ78FSz4IB4f9Dvo8meoWTtsrhRgEZEkKVmiCKY9AKMHxBNS6zaGnoMzdnOybWERkSQpGVYshmcuhA+ejsetusDxd2XEg+qqkkVEkqSqNmdKfCumbDbk1ohvwxRfkLUTUjfHIiJJUlWJInhjELx8NSTWwE67wUn/hF3ah06WsiwikiRVheWLYORvYcboeNzmeOhxG9QuCJsrxVlEJEnaXrPfjG/FlH8W7w3y8xLo8CvIyQmdLOVZRCRJ2lZRBG/cAS//Ob4V06AVnPwANG0XOlnasIhIkrQtVpbFO6T+59l43OYEOPY2yK8fNleasYhIkrS15r0Lj/8Svp4FuTXjWzEdz/VWzDawiEiStDXeeQSe7Q+VFVDYAk55wFUx28EiIknSllhTAS/8Md4pFaD1UfEGZTs0CBor3VlEJEn6IYvnxLdi5r4N5MARV8LBl7pBWRWwiEiStDmfjIuX5q5YBHV2ip+Yu4fPiqkqFhFJkjYliuD1W2HsNRAloOn+cMpDsNOuoZNlFIuIJEnfVbEURl0A74+Ix/ufAcf8A2rWDpsrA1lEJEn6tq8+gcfOgAUfxA+s6/436HCOS3OTxCIiSdI6M1+G4b+KNyur1wROfhB2LQ6dKqNZRCRJiiKYdFv81NwoAc07wikPQ0HT0MkynkVEkpTdVq+AUf1g+hPx+Cd94vkgNfLD5soSFhFJUvYq+wyG9YZ5pZCTF88Hcav2amURkSRlp9lvxpNSly2AHRrG80FaHhI6VdaxiEiSsk/pUHjmIqhcBU32hVMfdX+QQCwikqTskaiEl/8Mk26Px3v3gJ53QX69sLmymEVEkpQdVpbD8HNg5kvx+NA/QucBPi8mMIuIJCnzff1fGNoLvvwP1KgNPe+EtieGTiUsIpKkTPe/N+Cx3rD8K6jfFE4dCrscEDqV1rKISJIyV+m/4JkL40mpTfeD04ZBQbPQqfQtFhFJUuZJJOCVa+G1m+Lx3j3g+LuhVt2wubQRi4gkKbOsXgEjfgMfPB2PD7kUDr/KSakpyiIiScocSxfAv06Fz6dBbk049nbY/7TQqbQZFhFJUmZY8CE8egqUzYY6O0GvR2G3n4VOpR9gEZEkpb+Px8ITZ0FFOTRoBb2fgIatQqfSFrCISJLS27QH4dlLIKqEXQ+GXg/DDg1Cp9IWsohIktJTIgHjBsKr/4jH7U6N54TUqBU2l7aKRUSSlH7WVMDI38F7T8bjwy6HzpdDTk7YXNpqFhFJUnpZvgiG9YbZkyC3xtqVMaeHTqVtZBGRJKWPr/8Lj5wEX82E/ELo9RDs3jl0Km0Hi4gkKT18/nb84LplC6Cgebwypsk+oVNpO1lEJEmpb8aL8fLc1cuhyb5xCSloGjqVqoBFRJKU2qY9sHZ5bgJ2PxxOeQhqF4ROpSqS1I33Bw8eTLt27SgoKKCgoIDi4mJeeOGFZJ5SkpSGbh4zg9vGztzwYBTBK9fBMxfFJWT/3vGVEEtIRklqEWnevDnXX389U6dOZerUqRxxxBEcd9xxvP/++8k8rSQpzeTl5nDTt8tI5WoYdQFMvAGAN4vOhePugLyaAVMqGZJ6a6ZHjx4bjK+77joGDx7M5MmTadOmTTJPLUlKIxd2aQ3ATWNmUKNyBb/78lqY+RKVUQ7jWw+gyxmXBU6oZKm2OSKVlZU88cQTLFu2jOLi4k1+pqKigoqKivXj8vLy6oonSQrswi6tqbPqazq+dibkfsqKqBYvtymhxynnho6mJEp6EZk+fTrFxcWsXLmSevXqMWLECPbZZ9PLrUpKSrjmmmuSHUmSlIoWzeK8mb+B3E9ZFNXj/Mo/8rglJOPlRFEUJfMEq1atYvbs2SxevJjhw4dz7733MmHChE2WkU1dESkqKqKsrIyCAicnSVLGmvfveKOyZQv4LGrEOWuu4KPKH9G/657rb9sofZSXl1NYWLhFP7+TfkWkVq1a7LHHHgB06NCBKVOmcOutt3L33Xdv9Nn8/Hzy8/OTHUmSlEo+nRBv2b5qCR8mWjDpoLt58ehO3DZ2JjeNmQFgGclg1b6PSBRFG1z1kCRlsfeeghG/gcpVTE7sTenPBnP+UT8BNpzA+u2xMktSi8gVV1xB9+7dKSoqYsmSJQwbNozx48czevToZJ5WkpQO3roHnv8DEDGz4RFM+/G19O3adoOPrCsflYmkziJQQEktIvPnz6dPnz7MmzePwsJC2rVrx+jRo+natWsyTytJSmVRBONLYMLf4nGHc2h99I20zs3b5Me9EpLZklpE7rvvvmT+9ZKkdJOojK+CTF3786HzADjsMsjJCZtLwfisGUlS9VhTAU/9Gj4YCeTA0TfCT88LnUqBWUQkSclXsSReGTNrAuTWhBOGQNsTQqdSCrCISJKSa9lX8OhJMPdtqFkXTn0UWh0eOpVShEVEkpQ8ZZ/Dw8fDwo+gTgPo/SQ0bx86lVKIRUSSlBwLZ8YlpGwOFOwCfUbAznuFTqUUYxGRJFW9uaXwyImwfCE03AP6jIQdi0KnUgqyiEiSqtb/JsHQXlBRDk33g97Dod7OoVMpRVlEJElVZ8ZL8HgfWLMSdv0ZnDYMavvQUn0/i4gkqWpMfzJ+bkxiDbQ+Ck55EGrWCZ1KKS43dABJUgaY+k8Yfm5cQvY9OV6iawnRFrCISJK2z2u3wLOXABF0OAeOHwJ5NUOnUprw1owkadtEEbxyLbz6j3h8cH/o8iefG6OtYhGRJG29RAJGXwZvDYnHR14NB18SNJLSk0VEkrR1KtfAqH7w76FADhzzd+h4buhUSlMWEUnSlluzCoafAx+Ogpw86DkY9usVOpXSmEVEkrRlVi2P9wj5+GXIqwUn3Q97/yJ0KqU5i4gk6YdVLIGhp8L/XoOaO6x9gu4RoVMpA1hEJEmbt3wRPHoSfD4N8gvg9Mdh1+LQqZQhLCKSpO+39Et4uCfMfw/qNIA+T0Gzn4ROpQxiEZEkbVr5XHjwWPhqJtRrEj9Bt8k+oVMpw1hEJEkb+/p/8NCx8PV/oaA5nDkKGrYKnUoZyCIiSdrQwo/jElL+OezUMi4hO7YInUoZyiIiSfrG/A/goeNg2QJotBf88mkoaBo6lTKYRUSSFJtbCg8fDysWwY/2jeeE1G0UOpUynEVEkgSfTYWHT4CKMtilPZwxHOrsFDqVsoBFRJKy3X9fh6GnwKql0KI43iekdkHoVMoSFhFJymafjIN/nQZrVkDLQ+G0YVCrbuhUyiK5oQNIkgKZ8RIM7RWXkD26xldCLCGqZhYRScpG/3kOhp0OlRWw1zHxs2Nq1gmdSlnIIiJJ2eb9EfD4LyGxGvbpCac8CDXyQ6dSlrKISFI2efdxePJXkFgD+54CJ94HeTVDp1IWs4hIUrZ45xF46tcQJWD/M+D4uyDPNQsKyyIiSdlg6v3wdF8ggvZnw7G3Q25e6FSSRUSSMt5b98CzF8fvDzwffnEz5PrPv1KD/yVKUiZ74054/vfx++IL4OfXQ05O2EzSt1hEJClTvX4rvDggfn9wf+g20BKilOMsJUnKRBP/Dq9cG78/7HLofLklRCnJIiJJmWb832D8X+P3R1wFh/4hbB5pMywikpQpogjG/RUm3hCPj7waDr4kaCTph1hEJCkTRBGMvQZeuzkedxsInfqFzSRtAYuIJKW7KIIxf4JJt8Xjn18PB/02bCZpC1lEJCmdRRG8eCVMviMed78RDvx12EzSVrCISFK6iiIYPQDeHByPj7kJOp4TNpO0lSwikpSOoghe+CO8NSQe/+IW6HB20EjStrCISFK6SSTghT/AlHuBnPi5MQf0CZ1K2iYWEUlKJ4kEPH8pTP0nkAM974T9Tw+dStpmFhFJSheJBDx3CUx7gLiEDIb9TwudStouFhFJSgeJBDx7Ebz9EOTkQs+7YL9eoVNJ280iIkmpLpGAZ/rBO4/EJeT4u6HdKaFTSVXCIiJJqSyRgFH9oHRdCRkC7U4OnUqqMhYRSUpVicq1JeTRuISccA/se1LoVFKVsohIUipKVMLTF8C/h0JOHpx4D7Q9MXQqqcpZRCQp1SQq4em+8O9/rS0h90LbE0KnkpLCIiJJqeS7JeSk+6DN8aFTSUljEZGkVJGohJG/g3eHWUKUNSwikpQKNioh/4Q2PUOnkpIuN5l/eUlJCR07dqR+/fo0btyYnj178tFHHyXzlJKUfr5bQk6+3xKirJHUIjJhwgT69u3L5MmTGTNmDGvWrKFbt24sW7YsmaeVpPSxqRKyz3GhU0nVJieKoqi6Tvbll1/SuHFjJkyYwKGHHvqDny8vL6ewsJCysjIKCgqqIaEkVSNLiDLU1vz8rtY5ImVlZQA0aNBgk79fUVFBRUXF+nF5eXm15JKkamcJkYAk35r5tiiK6N+/PwcffDBt27bd5GdKSkooLCxc/yoqKqqueJJUfdYt0bWESNV3a6Zv374899xzvPbaazRv3nyTn9nUFZGioiJvzUjKHBvtE+LqGGWelLs1069fP0aNGsXEiRO/t4QA5Ofnk5+fXx2RJKn6rd+23RIirZPUIhJFEf369WPEiBGMHz+eli1bJvN0kpS61j3Abt2zY066zxIikeQi0rdvX4YOHcrTTz9N/fr1+eKLLwAoLCykTp06yTy1JKWORAJGXbj2Kbprnx3jjqkSkOQ5Ijk5OZs8fv/993PWWWf94J93+a6ktJdIwDP94J1HfIquskbKzBGpxi1KJCn1JBLwzIVrS0iuJUTahGpbvitJWSWRgGcvhncejkvICZYQaVMsIpJU1daVkLcfjEvI8UNg35NCp5JSkkVEkqpSIgHP9f9WCbkb2p0cOpWUsiwiklRVogievxSm3R+XkJ53QbtTQqeSUppFRJKqQhTBc5fC1H8COdBzMOzXK3QqKeVZRCRpe0URPP8HmHof35SQU0OnktKCRUSStkcUwQt/hCn3ADlw3B2w/2mhU0lpwyIiSdsqimD05fDWENaXkJ/0Dp1KSisWEUnaFlEEowfAm3fF42Nvt4RI28AiIklbK4rgpavgzcHxuMdtcECfsJmkNGURkaStsa6EvDEoHve4FdqfGTaTlMYsIpK0paIIxvzfNyXkFzdD+7OCRpLSnUVEkrZEFMHLf4ZJt8fjY26CDr8Km0nKABYRSfohUQRjr4HXb43HR/8dOp4TNpOUISwikrQ5UQSvXAuv3RyPu98IPz0vbCYpg1hEJOn7RBGMuw5e/Uc8/vnf4MBfh80kZRiLiCR9n/ElMPHG+P1RJXDQ+WHzSBnIIiJJmzL+epjwt/j9UX+F4t+FzSNlKIuIJH3XhBviqyEA3QZCcd+weaQMZhGRpG+bcGM8LwSg67XQqV/YPFKGs4hI0joTb4RxA+P3R14NP7swaBwpG1hEJAnilTGvrC0hXf4MB18SNo+UJSwikvTqTTD2L/H7Ln+CQ/qHzSNlEYuIpOz22i3xrqkAR1wFh1waNI6UbSwikrLX67fGz48BOPxKOPQPYfNIWcgiIik7vX4bjPlT/L7zFXDYH8PmkbKURURS9pl0O4z5v/h95wHQ+bKweaQsZhGRlF0mDYKXrorfH3Y5dL48bB4py1lEJGWPN+6Al66M3x92GRw+IGweSRYRSVnijTvgxSvi94f+Mb4lIyk4i4ikzPfGnd8qIX+Aw6+AnJywmSQBFhFJme6NO+HFtVc/Dvl9vEzXEiKlDIuIpMw1+a5vlZBL4w3LLCFSSrGISMpMk++C0WuX5R5yKRzxf5YQKQVZRCRlHkuIlDYsIpIyiyVESisWEUmZwxIipR2LiKTMMHmwJURKQxYRSenvjTth9Nqt2i0hUlqpETqAJG2X7+4T4hJdKa1YRCSlrw22bf+Dm5VJacgiIik9fffZMW7bLqUl54hISj+TbreESBnCIiIpvbx+K7x0VfzeEiKlPW/NSEofr90ML18dvz/scjh8QNA4krafRURSenj1Jhh7Tfy+8wDofHnYPJKqhEVEUuqb+Hd45dr4/eFXwmF/DJtHUpWxiEhKbRNugHHXxe+PuCpepispY1hEJKWu8dfD+JL4fZc/xbumSsooFhFJqSeKYNxfYeIN8fjIa+Dgi4NGkpQcFhFJqSWK4JWB8Orf43G3gdCpX9hMkpLGIiIpdURRvDz39Vvi8VF/heK+IRNJSjKLiKTUEEXxRmVvDIrHP/8bHHR+2EySks4iIim8KILRA+DNwfH46L/DT88Lm0lStbCISAoriuD5P8CUe+LxL26BDmcHjSSp+lhEJIWTSMDzl8LUfwI5cOztcECf0KkkVaOkPvRu4sSJ9OjRg2bNmpGTk8PIkSOTeTpJ6SSRgGcu/KaE9LzTEiJloaQWkWXLlrHffvsxaNCgZJ5GUrpJVMLTv4N3HoacXDhhCOx/euhUkgJI6q2Z7t27071792SeQlK6qVwDI34D7z0JOXlw4r3Q9oTQqSQFklJzRCoqKqioqFg/Li8vD5hGUpWrXA3Dz4UPRkJuDTjpftjn2NCpJAWU1FszW6ukpITCwsL1r6KiotCRJFWVNRXw+JlrS0hNOOVhS4ik1CoiAwYMoKysbP1rzpw5oSNJqgqrV8JjZ8BHz0FePpw6FH58dOhUklJASt2ayc/PJz8/P3QMSVVp1XIYdhp8Oh5q1IHT/gWtDg+dSlKKSKkiIinDVCyFf50K/30VataF3o/DbgeHTiUphSS1iCxdupSPP/54/XjWrFmUlpbSoEEDWrRokcxTSwptZRk8ejLMeRNq1YcznoQWB4VOJSnFJLWITJ06lcMP/+YSbP/+/QE488wzeeCBB5J5akkhLV8Ej5wAc9+B2oVwxgho3j50KkkpKKlFpHPnzkRRlMxTSEo1yxbCQz1h/nSo0wB+ORKa7hc6laQU5RwRSVVnyRfw4LGw8COo2xjOHAWN9w6dSlIKs4hIqhpln8UlZNEnUL8ZnPkMNNojdCpJKc4iImn7LZoVl5Cy2VDYIr4S0qBl6FSS0oBFRNL2+XIGPHQsLJkHDVrFJaSweehUktKERUTStvviPXjoOFi+EHbeG375NNRvEjqVpDRiEZG0bT5/Gx4+HlYuhh+1gz4joW7D0KkkpRmLiKSt979J8OgpsGoJNO8IvZ+EOjuGTiUpDVlEJG2dj8fCsN6wZgXsdkj87Jj8+qFTSUpTFhFJW+7DZ+HJs6FyFbTuBqc8BDXrhE4lKY1ZRCRtmXefgBG/gagS9jkOTrgXatQKnUpSmssNHUBSGph6Pzx1XlxC9jsdTvynJURSlbCISNq8SbfDsxcDEXQ8F467A/K8mCqpaviviaRNiyIY91eYeEM8PvgS6PJnyMkJm0tSRrGISNpYIgEvXgFvDo7HXf4Eh1waNpOkjGQRkbShyjXwzEVQ+kg8Pvrv8NPzwmaSlLEsIpK+saYChp8DHz4DObnxfJD9Tw+dSlIGs4hIiq1aFm9U9uk4yKsFJ/0T9u4ROpWkDGcRkQQrvoahvWDOm1CzLpz6KLQ6PHQqSVnAIiJluyXz4ZETYP57ULsQeg+Hoo6hU0nKEhYRKZt9/V94qCd8PQvqNYE+I6BJm9CpJGURi4iUreZ/AA8fD0u/gB1bQJ+R0LBV6FSSsoxFRMpGc6bAoyfBysXQeB844ykoaBo6laQsZBGRss3HY+GxM2D1cmjeEU5/HHZoEDqVpCxlEZGyyXvD4anfQGI1tDoCej0CteqGTiUpi/nQOylbvHUPPHlOXELaHA+nPWYJkRScV0SkTBdFMP56mHB9PO54LnS/AXLzwuaSJCwiUmZLVMILf4Qp98bjzgPgsMt8gq6klGERkTLVmgp46tfwwUggB46+0YfXSUo5FhEpE60sh8d6w6yJkFsTThgCbU8InUqSNmIRkTLNkvnw6InwxXSoVR9OfQR27xw6lSRtkkVEyiRffRLvlrr4f1B3Z+j9JDTbP3QqSfpeFhEpU3z+Njx6MixfCDu1hD5PQYPdQ6eSpM2yiEiZYObL8PgvYfUyaLpffCWkXuPQqSTpB1lEpHRXOhRG9YPEGtj9cOj1MOTXD51KkraIRURKV1EEr/4DXrk2HrfrBccOghq1wuaSpK1gEZHS0Xc3KvvZxXDk1W5UJintWESkdLNqOQw/Bz56HsiBn18PB50fOpUkbROLiJROli2Eob3g86lQozaccA/sc2zoVJK0zSwiUrr46hN45ET4ehbU2QlOGwYtDgqdSpK2i0VESgdzpsC/esHyr2DHFnDGU9CodehUkrTdLCJSqvtgFDx1HqxZCU33h9Mfh/pNQqeSpCphEZFSVRTB5DvhxSuBCPb8OZx4H+TXC51MkqqMRURKRYlKGD0A3ro7Hnc8F37+N8jz/7KSMov/qkmppmJpfCvmo+fjcbeBUHyBe4RIykgWESmVlM+DoafAF+/Gy3OPvxva9AydSpKSxiIipYovpsd7hJR/Djs0ipfnFnUMnUqSksoiIqWCGS/Bk2fDqqXQaC/o/TjstFvoVJKUdBYRKaQogreGwOjLIUpAy0PhlIehzo6hk0lStbCISKFUroHRl33z4LqfnAHH3OzTcyVlFYuIFMKKxfGtmE9eAXKg6zXQ6UJXxkjKOhYRqbotmhVPSl34EdTcIX5w3d6/CJ1KkoKwiEjV6b+vwWN9YMUiqN8MTh8GTfcLnUqSgrGISNVl2gPw3KWQWAPNfgKn/gsKmoZOJUlBWUSkZKtcAy9dCW/eFY/bngjH3QE164TNJUkpwCIiJdOKxfDkr+CTsfH48Kvg0N87KVWS1rKISMny5QwYdhp89XE8KfX4u2Cf40KnkqSUklsdJ7nzzjtp2bIltWvXpn379rz66qvVcVopnBkvwr1d4hJSWAS/Gm0JkaRNSHoReeyxx7j44ou58soreeeddzjkkEPo3r07s2fPTvappeoXRfDqTfHy3IpyaNEJzhvnyhhJ+h45URRFyTzBgQceyAEHHMDgwYPXH9t7773p2bMnJSUlm/2z5eXlFBYWUlZWRkFBQTJjSttv1XIYdQG8Nzwetz8but/gTqmSss7W/PxO6hyRVatWMW3aNC6//PINjnfr1o1JkyZt9PmKigoqKirWj8vLy5MZT6o6X/8Xhp0B86dDbo24gHQ8J3QqSUp5Sb01s3DhQiorK2nSpMkGx5s0acIXX3yx0edLSkooLCxc/yoqKkpmPKlqfDIOhnSOS0jdneGXoywhkrSFqmWyas53lipGUbTRMYABAwZQVla2/jVnzpzqiCdtmyiC12+DR06AFV9DswPg1+Nht5+FTiZJaSOpt2YaNWpEXl7eRlc/FixYsNFVEoD8/Hzy8/OTGUmqGhVLYVQ/eP+peLx/bzjmJqhZO2wuSUozSb0iUqtWLdq3b8+YMWM2OD5mzBg6deqUzFNLybPwY7j3yLiE5NaA7jeu3SnVEiJJWyvpG5r179+fPn360KFDB4qLixkyZAizZ8/m/PPPT/appar3n+dgxPnx0tx6TeDkB2HX4tCpJCltJb2I9OrVi6+++oq//OUvzJs3j7Zt2/L888+z6667JvvUUtVJVMIrA+G1m+Jxi2I4+QGo/6OgsSQp3SV9H5Ht4T4iSglLF8Dwc2DWxHh84G+h27WQVzNsLklKUSmzj4iU9v73Bjx5NiyZBzXrwrG3wb4nhU4lSRnDIiJtShTBG3fAmD9BVAmN9oJTHoLGPw6dTJIyikVE+q4Vi+HpvvCfZ+Nx25Ogx62QXy9oLEnKRBYR6ds+nwZPnAWLZ0NeLTjqr9DxXNjEBnySpO1nEZEgvhXz5t3w0lWQWA077hqvitnlgNDJJCmjWUSkFV/Hu6R++Ew83vtYOPZ2qLNj0FiSlA0sIspuc96CJ8+BstmQWzO+FfPT87wVI0nVxCKi7JRIwOs3wyvXxatidmoJJ/3TWzGSVM0sIso+S+bDiF/Dp+Pj8b4nxw+sq+2meZJU3Swiyi4fjY6X5i5fCDV3gKP/Dvuf7q0YSQrEIqLssHoFvPR/MOWeeNxkXzjpPth5r7C5JCnLWUSU+b54D4afC19+GI8P6gtH/hlq5IfNJUmyiCiDJRLw5mB4+RqorIC6jeH4wbDHkaGTSZLWsogoM5V9BiPOh/++Go9bHwXH3QH1dg6bS5K0AYuIMs+7T8Bzl0JFWTwh9ajroP3ZTkiVpBRkEVHmWL4oLiDvPxWPd2kPxw+BRnuEzSVJ+l4WEWWGj16AURfCsgWQkweH/REO+T3k+Z+4JKUy/5VWeltZBqMHQOmj8bjRXvGE1F3ah80lSdoiFhGlr49fhlEXQflnQA50ugAOvwpq1g6dTJK0hSwiSj8rvoYXr4LSR+Jxg92h52BocVDYXJKkrWYRUXr5z/Pw7CWw9AsgBw76LRxxFdSqGzqZJGkbWESUHpZ+CaMvh/eejMcN94j3BfEqiCSlNYuIUlsUwb//BS9eEd+SycmFTv2g8wCoWSd0OknSdrKIKHUtmgXPXgyfjo/HTdrCsbe5IkaSMohFRKmncjW8MQjG/w3WrIAateGwy+IrIXk1Q6eTJFUhi4hSy/8mwbP9v3lS7m6HQI9boWGrsLkkSUlhEVFqWPYVvPwneGftktwdGkK3gbDfaT4jRpIymEVEYSUq4e2HYOxfYMWi+NgBZ8KRV8MODYJGkyQln0VE4Xw2FZ7/Pcx9Jx43bgO/uMkluZKURSwiqn7LFsLLV8M7D8fj/AI4/AroeJ4PqZOkLOO/+qo+a1bBW3fDhBugojw+tt/p0PUaqNc4bDZJUhAWESVfFMGM0fDilbDok/hY0/2g+43Q4sCw2SRJQVlElFxfvAcvXQWfjovH9ZpAlz/FV0Jyc8NmkyQFZxFRcpR9DuOug9KhQAR5taC4LxxyKeTXD51OkpQiLCKqWivL4fVb4I07411RAdqcEF8FadAyaDRJUuqxiKhqrF4BU+6FV2/6Zj+QFp3iTcma+2wYSdKmWUS0fSrXQOkj8XNhlsyNjzVsHa+E2etod0WVJG2WRUTbJlEJ7w2H8dd/sxKmoDl0vjzelt39QCRJW8CfFto6iUp4fwRM+BssnBEf26EhHPJ76PArqFk7bD5JUlqxiGjLrCsgE2+EL/8TH6u9I3TqBwf+xpUwkqRtYhHR5q1ZBe8+Bq/d/M0tmNqFUHwBHHg+1C4Im0+SlNYsItq0VcvjZ8G8fhuUfxYfq7MTHPjb+ApInR2DxpMkZQaLiDa09Et4a0i8FHfdMtx6TeIrIB1+Bfn1wuaTJGUUi4hiC2fCG3fAv/8Fa1bGx3bcNZ4D8pM+TkKVJCWFRSSbJRLwySvw5mD4+OVvjjc7AH52Iex9LOTmhcsnScp4FpFstLIM/v0YvHU3fPXx2oM5sOfP4ysgu3ZyIzJJUrWwiGSTz9+Gqf+MNyJbvTw+ll8Q33r56bnQYPew+SRJWccikulWlsXFY9oDMO/f3xxvtBf89Lx4F1QnoEqSArGIZKJEAmZNgNJH4cNnvpl8mlcL9ukJHc6GFsXefpEkBWcRySTz34fpT8C7T3yz9wfAzj+Gn5wB+50OdRuGyydJ0ndYRNLdolnxrZfpT8KXH35zvHYh7Hsy7H96vArGqx+SpBRkEUlHC/4T33L5cBR88e43x/NqQetusO9JsGd39/6QJKU8i0g6qFwDn70FM0bDRy9889RbgJw8aHlIfPXjx79w63VJUlqxiKSqJfPh0/Ew80X4eCysXPzN7+XWhFaHxxuO7XW08z4kSWnLIpIqVpbD/ybFq10+HQ8LPtjw9+vsBHt0hT2PgtZd4zkgkiSlOYtIKOXzYPYbMHsyzJ4Ur3iJEt/6QA40bQetusTlo3lHt1uXJGWcpBaR6667jueee47S0lJq1arF4sWLk3m6H3TzmBm8OesrOrVqxIVdWm/we7eNncmkTxZyYMuGXNJ1z6o98bKF8WZic9+Gz9+Bue/Akrkbf26n3WD3zvFrt0O95SJJynhJLSKrVq3i5JNPpri4mPvuuy+Zp9oiebk5TP50EZM/jR9vv66M3DZ2JjeNiSeAdmrVaNv+8iiCZV/CV5/Ez29Z8CEseB/mfwDLFmz8+ZxcaNI23lisxUHxrwVNt+3ckiSlqaQWkWuuuQaABx54IJmn2WLrisdNY2asLx7rxgD9u+650ZUSIC4ZFUtg+cL46kb551D2+dpfP4Ov/xvv57FqyfecOQcatIz389jlgPjXpu2gVt0q/l8oSVJ6Sak5IhUVFVRUVKwfl5eXV/k5LuzSmp2Wz2L1m/eRN66SPBKU1EjQbpd6tFlcB4Yth1VLYdWy+LXia1j+FVSu2oK/PQd2LIIGraDx3mtfbaDxjy0dkiRtQkoVkZKSkvVXUZKpz965MG30hgfnr31tTs0dYIdGUNAMCneBgrWvHVtAw1bxHI8a+UlKLUlS5tnqInL11Vf/YFmYMmUKHTp02OowAwYMoH///uvH5eXlFBUVbfXf80Me/E8OS9ccyxrySES5rCGPTq0b87M9fxRfuahVL/615g7xstkdGsavWjtUeRZJkrLZVheRCy64gFNPPXWzn9ltt922KUx+fj75+cm9onDb2Jnc9NpK4FT6d92TPODWMTO44z/Qv2hPLuy0iTkikiQpKba6iDRq1IhGjbZxZUlg314d892Jqd+ewLrJCauSJKnKJXWOyOzZs1m0aBGzZ8+msrKS0tJSAPbYYw/q1auXzFNvUmUi4qDdG2y0j8i695M+WUhlIqr2XJIkZaucKIqS9pP3rLPO4sEHH9zo+Lhx4+jcufMP/vny8nIKCwspKyujoKAgCQklSVJV25qf30ktItvLIiJJUvrZmp/fudWUSZIkaSMWEUmSFIxFRJIkBWMRkSRJwVhEJElSMBYRSZIUjEVEkiQFYxGRJEnBWEQkSVIwSX3WzPZat+lreXl54CSSJGlLrfu5vSWbt6d0EVmyZAkARUVFgZNIkqSttWTJEgoLCzf7mZR+1kwikWDu3LnUr1+fnJycKv27y8vLKSoqYs6cOT7HJon8nquH33P18HuuHn7P1SdZ33UURSxZsoRmzZqRm7v5WSApfUUkNzeX5s2bJ/UcBQUF/odeDfyeq4ffc/Xwe64efs/VJxnf9Q9dCVnHyaqSJCkYi4gkSQoma4tIfn4+f/7zn8nPzw8dJaP5PVcPv+fq4fdcPfyeq08qfNcpPVlVkiRltqy9IiJJksKziEiSpGAsIpIkKRiLiCRJCiYri8idd95Jy5YtqV27Nu3bt+fVV18NHSnjTJw4kR49etCsWTNycnIYOXJk6EgZqaSkhI4dO1K/fn0aN25Mz549+eijj0LHyjiDBw+mXbt26zd9Ki4u5oUXXggdK+OVlJSQk5PDxRdfHDpKRrn66qvJycnZ4PWjH/0oWJ6sKyKPPfYYF198MVdeeSXvvPMOhxxyCN27d2f27Nmho2WUZcuWsd9++zFo0KDQUTLahAkT6Nu3L5MnT2bMmDGsWbOGbt26sWzZstDRMkrz5s25/vrrmTp1KlOnTuWII47guOOO4/333w8dLWNNmTKFIUOG0K5du9BRMlKbNm2YN2/e+tf06dODZcm65bsHHnggBxxwAIMHD15/bO+996Znz56UlJQETJa5cnJyGDFiBD179gwdJeN9+eWXNG7cmAkTJnDooYeGjpPRGjRowI033sg555wTOkrGWbp0KQcccAB33nknAwcOZP/99+eWW24JHStjXH311YwcOZLS0tLQUYAsuyKyatUqpk2bRrdu3TY43q1bNyZNmhQolVR1ysrKgPiHpJKjsrKSYcOGsWzZMoqLi0PHyUh9+/blmGOO4cgjjwwdJWPNnDmTZs2a0bJlS0499VQ+/fTTYFlS+qF3VW3hwoVUVlbSpEmTDY43adKEL774IlAqqWpEUUT//v05+OCDadu2beg4GWf69OkUFxezcuVK6tWrx4gRI9hnn31Cx8o4w4YN4+2332bKlCmho2SsAw88kIceeog999yT+fPnM3DgQDp16sT7779Pw4YNqz1PVhWRdXJycjYYR1G00TEp3VxwwQW8++67vPbaa6GjZKS99tqL0tJSFi9ezPDhwznzzDOZMGGCZaQKzZkzh4suuoiXXnqJ2rVrh46Tsbp3777+/b777ktxcTGtWrXiwQcfpH///tWeJ6uKSKNGjcjLy9vo6seCBQs2ukoipZN+/foxatQoJk6cSPPmzUPHyUi1atVijz32AKBDhw5MmTKFW2+9lbvvvjtwsswxbdo0FixYQPv27dcfq6ysZOLEiQwaNIiKigry8vICJsxMdevWZd9992XmzJlBzp9Vc0Rq1apF+/btGTNmzAbHx4wZQ6dOnQKlkrZdFEVccMEFPPXUU7zyyiu0bNkydKSsEUURFRUVoWNklC5dujB9+nRKS0vXvzp06EDv3r0pLS21hCRJRUUFH374IU2bNg1y/qy6IgLQv39/+vTpQ4cOHSguLmbIkCHMnj2b888/P3S0jLJ06VI+/vjj9eNZs2ZRWlpKgwYNaNGiRcBkmaVv374MHTqUp59+mvr166+/2ldYWEidOnUCp8scV1xxBd27d6eoqIglS5YwbNgwxo8fz+jRo0NHyyj169ffaH5T3bp1adiwofOeqtDvf/97evToQYsWLViwYAEDBw6kvLycM888M0ierCsivXr14quvvuIvf/kL8+bNo23btjz//PPsuuuuoaNllKlTp3L44YevH6+773jmmWfywAMPBEqVedYtQ+/cufMGx++//37OOuus6g+UoebPn0+fPn2YN28ehYWFtGvXjtGjR9O1a9fQ0aSt9tlnn3HaaaexcOFCdt55Zw466CAmT54c7Odg1u0jIkmSUkdWzRGRJEmpxSIiSZKCsYhIkqRgLCKSJCkYi4gkSQrGIiJJkoKxiEiSpGAsIpIkKRiLiCRJCsYiIkmSgrGISJKkYCwikiQpmP8HGWCGrHwHVLAAAAAASUVORK5CYII=",
"text/plain": [
"<Figure size 640x480 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Vanishing cells?\n",
"x = [0., 1.e-2, 4., 5.] # Positiong of data, x[i]\n",
"y = [-1., -1., 3., 4.] # Values of data, y[i]\n",
"q = numpy.linspace(x[0],x[3],100) # Sample interval between x[1]..x[2]\n",
"plt.plot(x,y,'x');\n",
"plt.plot(q, lagrange_interp(x,y,q));"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAiIAAAGdCAYAAAAvwBgXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAvHElEQVR4nO3deXhV5aG28TsJJCAkwYCASEBEsSqiMqhBUZRBqbVSax1qqVp7PLbgUOpXRc+pQ9F4ag9VarVSrdo6YK3icFQUKZMiyiCK1go4ATKDJGFKQrK/PxYEEUSG7P3u4f5d17r2u1Z2WI9bLvfjGt6VFYvFYkiSJAWQHTqAJEnKXBYRSZIUjEVEkiQFYxGRJEnBWEQkSVIwFhFJkhSMRUSSJAVjEZEkScE0CB1gZ2pra1m8eDH5+flkZWWFjiNJknZBLBajoqKCNm3akJ2982MeSV1EFi9eTHFxcegYkiRpDyxcuJC2bdvu9D1JXUTy8/OB6B+koKAgcBpJkrQrysvLKS4urvse35mkLiJbTscUFBRYRCRJSjG7clmFF6tKkqRgLCKSJCkYi4gkSQrGIiJJkoKxiEiSpGAsIpIkKRiLiCRJCsYiIkmSgrGISJKkYBJWREpLS8nKyuLqq69O1C4lSdJX/H7cXEaOn7fDn40cP4/fj5ub0DwJKSLTp09n1KhRdOnSJRG7kyRJXyMnO4sROygjI8fPY8S4ueRkJ/Zp93EvImvXruXCCy/kz3/+M/vuu2+8dydJknbiyj6HMLRfp23KyJYSMrRfJ67sc0hC88T9oXeDBw/mjDPOoG/fvgwfPnyn762srKSysrJuvby8PN7xJEnKOFf2OYS86jIOnXgpl03ozyubjmZov0MTXkIgzkVk9OjRzJo1i+nTp+/S+0tLS7n55pvjGUmSJAH/uc9EyHmH/WrXMDGna5ASAnE8NbNw4UKuuuoqHnnkERo1arRLvzNs2DDKysrqloULF8YrniRJmWtTJeum3APAg7Xfoaom9rUXsMZb3I6IzJw5k+XLl9OtW7e6bTU1NUyePJm7776byspKcnJytvmdvLw88vLy4hVJkiQB4564m37Vq1ib25Lf/ffNtJv4KSM23y2TNteI9OnThzlz5myz7ZJLLuFb3/oW11577XYlRJIkxd/IV+dy2od/gWxoevIQyGlYVz5ClJG4FZH8/Hw6d+68zbYmTZrQvHnz7bZLkqTEaLt6KodmL4LcfOh2cd32LeWjpjaW0Dxxv2tGkiQlj7M3PB0Nuv4YGhVu87O0u2vmqyZOnJjI3UmSpC9b8i58MgmycuD4y0OnAXzWjCRJmeONP0avRwyEZu2CRtnCIiJJUiYo+xze+0c0LhkSNsuXWEQkScoEb90HtZug/YlwQNfQaepYRCRJSncby2HGQ9G45xVBo3yVRUSSpHQ362GoLIMWneCQ/qHTbMMiIklSOttUBW9E07nT80rITq6v/uRKI0mS6td7/4CKxdC0NXQ5N3Sa7VhEJElKV7W18PrIaHz8z6BB8j3PzSIiSVK6mj8OVnwQTefe/ZLQaXbIIiJJUrp6/a7otfsl203nniwsIpIkpaOF0+Gz1yG7YXRaJklZRCRJSkdTNx8N6XIuFLQJm2UnLCKSJKWblfPhg/+Lxkk2gdlXWUQkSUo3U0cCMTjkNGh5WOg0O2URkSQpnZQvgXcej8Yn/iJsll1gEZEkKZ1M+yPUVEG7EmhfEjrNN7KISJKULjZ8ATMejMYpcDQELCKSJKWP6fdD1VpoeXjSPdzu61hEJElKB1XrYdqfovGJv4CsrLB5dpFFRJKkdPD2I7B+JTRrB0ecHTrNLrOISJKU6mqqYeofonHPKyGnQdg8u8EiIklSqnvvaShbAPu0gGN+FDrNbrGISJKUympr4fU7o/HxP4OGjYPG2V0WEUmSUtncl2D5vyA3H3r8NHSa3WYRkSQpVcViMPl30fjYn0LjZkHj7AmLiCRJqerjCbB4FjRoDMcPDp1mj1hEJElKVVNGRK/dLoKm+4XNsocsIpIkpaIFb8KnUyC7IfS8InSaPWYRkSQpFU3ZfG3IUedDYduwWfaCRUSSpFSz5B2Y9wpkZafMw+2+jkVEkqRUM+V/o9cjzobmHcNm2UsWEUmSUsmKD+Ffz0XjXr8Mm6UeWEQkSUolr/0eiMGhZ0Crw0On2WsWEUmSUsXqj+Hdv0fjk1L/aAhYRCRJSh1TRkCsBg7uCwd0C52mXlhEJElKBWsWwDuPR+OTrw2bpR5ZRCRJSgWv/R5qN0GHk6H42NBp6o1FRJKkZFf2Obz9SDROo6MhYBGRJCn5vX4X1FRB+xPgwBNCp6lXFhFJkpJZxVKY9XA0PvlXYbPEgUVEkqRkNvUPsGkjtD02uj4kzVhEJElKVmtXwIy/ROOTr4WsrLB54sAiIklSspo6EqrXQ5uucHCf0GniwiIiSVIyWrsCpt8fjXsPS8ujIWARkSQpOU29a+vRkEP6hU4TNxYRSZKSzdoV8Fb6Hw0Bi4gkScln6l2waUPaHw0Bi4gkScklg46GgEVEkqTksuVoyAHd0v5oCFhEJElKHhl2NAQsIpIkJY/X79x6NOTgvqHTJIRFRJKkZFCxdOu8ISdflxFHQ8AiIklScpgyYvMzZXpkxLUhW1hEJEkKrexzmPlgND7lhow5GgIWEUmSwpvyO6ipgvYnwEG9Q6dJKIuIJEkhffEZzPpbNM6woyFgEZEkKazJv4Xa6uhIyIEnhE6TcHEtIvfeey9dunShoKCAgoICSkpKeOmll+K5S0mSUseqj2D249H4lP8KmyWQuBaRtm3bcvvttzNjxgxmzJjBqaeeyllnncX7778fz91KkpQaJv0PxGrgkP5Q3CN0miCyYrFYLJE7LCoq4o477uDSSy/9xveWl5dTWFhIWVkZBQUFCUgnSVKCrPgQ7jkeYrVw2URoc0zoRPVmd76/GyQoEzU1NTz55JOsW7eOkpKSHb6nsrKSysrKuvXy8vJExZMkKbH+OTwqId/6TlqVkN0V94tV58yZQ9OmTcnLy+Pyyy9nzJgxHH744Tt8b2lpKYWFhXVLcXFxvONJkpR4i9+GD54DsqI7ZTJY3E/NVFVVsWDBAtasWcNTTz3F/fffz6RJk3ZYRnZ0RKS4uNhTM5Kk9PLI92H+q9DlPDh7VOg09W53Ts0k/BqRvn370rFjR+67775vfK/XiEiS0s5nU+HBAZDdAIZMh6KDQieqd7vz/Z3weURisdg2Rz0kScoYsRiMvyUad/1xWpaQ3RXXi1Wvv/56BgwYQHFxMRUVFYwePZqJEycyduzYeO5WkqTkNH88LHgDGjSCk/5f6DRJIa5FZNmyZQwaNIglS5ZQWFhIly5dGDt2LP36Zc5TBSVJAqC2FsbfHI17/BQK2oTNkyTiWkQeeOCBeP7xkiSljg+eg6XvQm5TOHFo6DRJw2fNSJIUbzWbYMKt0bhkCDRpHjZPErGISJIUb7MfhZVzoXERlAwOnSapWEQkSYqn6g0w8fZofNI10MjpKL7MIiJJUjy9NQoqFkNhMXT/5uesZRqLiCRJ8bJhDUwZEY1PuR4aNgoaJxlZRCRJipfX74SNa6Dl4dF07tqORUSSpHgoXwLT/hSN+/wasnPC5klSFhFJkuJh0u2waQMUHw+dTg+dJmlZRCRJqm8r58Osv0XjvjdBVlbQOMnMIiJJUn0bfxPEaqIjIe1LQqdJahYRSZLq04I34YPnISs7OhqinbKISJJUX2IxGPfraHz0hdDysLB5UoBFRJKk+vLvF2DhNGjQOJo3RN/IIiJJUn2oqYZXb4zGJYOhoE3YPCnCIiJJUn2Y9VdYNR/2aQ4nXBU6TcqwiEiStLcqK2BiaTQ++TofbLcbLCKSJO2tqXfDuhVQdBB0uzh0mpRiEZEkaW+UL4GpI6NxnxuhQW7YPCnGIiJJ0t6YMByq10PbHnD4WaHTpByLiCRJe2rpHHj70Wh82m1O5b4HLCKSJO2JWAxe+S8gBkd8D4qPDZ0oJVlEJEnaE/PGwccTISfXqdz3gkVEkqTdVbNp89EQ4Lj/hH0PDBonlVlEJEnaXbMehpUfQuMi6HVN6DQpzSIiSdLu2FgOE26Lxr2vg8bNgsZJdRYRSZJ2x2sjYP1KaH4wdP9J6DQpzyIiSdKuWv0JvPHHaNzvN5DTMGyeNGARkSRpV437NdRUwUG94dABodOkBYuIJEm74tPX4IPnICvbycvqkUVEkqRvUlsDY4dF426XQKsjwuZJIxYRSZK+yexHYem7kFcIp1wfOk1asYhIkrQzG8th/G+ice9roUmLsHnSjEVEkqSdeW0ErFse3a7b4z9Cp0k7FhFJkr7O6o+33q7b/1ZokBs2TxqyiEiS9HVeviG6XbfjqdDptNBp0pJFRJKkHZn3Knz4ImQ3gNNv93bdOLGISJL0VZuqYOx10fi4y2G/Q8PmSWMWEUmSvuqtUbBqHjTZD07+Veg0ac0iIknSl1Usg4m3R+M+N0KjwrB50pxFRJKkLxt/C1RVQJtj4OgLQ6dJexYRSZK2WDQTZj8SjQfcAdl+Tcabn7AkSRA9T+bFX0bjoy6A4h5h82QIi4gkSQCz/gqL34a8Auh7c+g0GcMiIknS+tUwfnP5OOV6yG8VNk8GsYhIkjT+FtjwBbQ8wufJJJhFRJKU2T6fBTMfisbfvgNyGgSNk2ksIpKkzFVbCy9eA8Sgy3lw4AmhE2Uci4gkKXO9/Tf4fCbk5kO/W0KnyUgWEUlSZlq/Gl69KRqfcj3ktw4aJ1NZRCRJmenVG2HDamh5OBzrBaqhWEQkSZlnwZvRvCEAZ4yAnIZh82Qwi4gkKbPUbIIXhkbjo38E7UvC5slwFhFJUmZ56z5Y9h403tcLVJOARUSSlDnKPocJt0XjvjdDk+Zh88giIknKIC8Pg6q10PZYOGZQ6DTCIiJJyhTzXoV/PQtZOfCd30O2X4HJIK7/FkpLS+nRowf5+fm0bNmSgQMH8uGHH8Zzl5Ikba9q/dYLVI//GbTuHDaP6sS1iEyaNInBgwczbdo0xo0bx6ZNm+jfvz/r1q2L524lSdrWpP+BNZ9BQVvofV3oNPqSuD7ZZ+zYsdusP/jgg7Rs2ZKZM2dy0kknxXPXkiRFlr4HU/8Qjc/4HeTlh82jbST0EYNlZWUAFBUV7fDnlZWVVFZW1q2Xl5cnJJckKU3V1sDzV0GsBg47Ew4dEDqRviJhV+rEYjGGDh3KiSeeSOfOOz43V1paSmFhYd1SXFycqHiSpHQ04y/w+YzooXYDfhs6jXYgYUVkyJAhvPvuuzz++ONf+55hw4ZRVlZWtyxcuDBR8SRJ6aZ8Mbx6czTueyMUtAmbRzuUkFMzV1xxBc899xyTJ0+mbdu2X/u+vLw88vLyEhFJkpTuXvoVVFXAAd2h+09Cp9HXiGsRicViXHHFFYwZM4aJEyfSoUOHeO5OkqTIv1+AD56P5gw58y7IzgmdSF8jrkVk8ODBPPbYYzz77LPk5+ezdOlSAAoLC2ncuHE8dy1JylQby+CFX0bjnkOcMyTJZcVisVjc/vCsrB1uf/DBB7n44ou/8ffLy8spLCykrKyMgoKCek4nSUpLz18NMx+EooPgZ1Ohof/jm2i78/0d91MzkiQlzKevRSUE4MyRlpAU4ET7kqT0UL0RnrsyGne9CDr0CptHu8QiIklKD5P+B1Z/BE1bQ79bQqfRLrKISJJS39I58Ppd0fiM30HjZkHjaNdZRCRJqa1mEzw7ZPM07t+NpnJXyrCISJJS29S7YMlsaNQMvn1H6DTaTRYRSVLqWv5vmHh7ND79dshvHTaPdptFRJKUmmpr4NnBUFMFh/SHo84PnUh7wCIiSUpN0+6JnqybVwDfuRO+ZhJNJTeLiCQp9aycD/8cHo1PuxUKDwibR3vMIiJJSi21tfDcENi0EQ46BY4ZFDqR9oJFRJKUWt66Dxa8AblN4bsjPSWT4iwikqTUsXIevHpTNO53CzRrFzSO9p5FRJKUGmpr4JmfbT0l0/0noROpHlhEJEmpYepIWDQ9ukvmrLs9JZMmLCKSpOS37F8w4bZofHopFLYNm0f1xiIiSUpuNdXwzOXRxGWdToejLwydSPXIIiJJSm5TRsCSd6JnyZx5l6dk0oxFRJKUvD6fBZN/G42//TufJZOGLCKSpORUvQHG/CfUboLDvgtHnhM6keLAIiJJSk6v3gQr50LTVp6SSWMWEUlS8vloArz5p2h81h9hn6KweRQ3FhFJUnLZ8AU88/No3P0ncEi/sHkUVxYRSVJyefH/QcViKDoI+g8PnUZxZhGRJCWP956COU9CVjZ8bxTkNgmdSHFmEZEkJYc1C+H5X0TjXr+E4h5h8yghLCKSpPBqa2DM5VBZBgd0g5OvDZ1ICWIRkSSF9/pd8Nlr0LAJnP1nyGkYOpESxCIiSQrr81kw4dZo/O3fQvOOYfMooSwikqRwqtbBUz+NZk89/CwfaJeBLCKSpHDGDoPVH0F+G/jOnc6emoEsIpKkMN4fA7MeBrLg7PucPTVDWUQkSYm3ZgE8d1U0PvFq6HBS0DgKxyIiSUqsmk3RdSGVZdC2B5xyQ+hECsgiIklKrImlsPBNyCuA79/vrboZziIiSUqcTybDlP+NxmfeCfseGDKNkoBFRJKUGOtWwdOXATE4ZhB0/n7oREoCFhFJUvzV1sIzl0PFEmjRCQb8T+hEShIWEUlS/E0dCfNegQaN4JwHfaqu6lhEJEnxteBNGH9LND79dmjdOWweJRWLiCQpftavhn9cArEa6HwOdLs4dCIlGYuIJCk+amthzOVQ/jkUdYzuknEKd32FRUSSFB9v3A3zXoacPDj3YcjLD51IScgiIkmqfwumwfibo/GA26H1kWHzKGlZRCRJ9WvtCnjyYqjdFM0V0u2S0ImUxCwikqT6U1sDT/1k83whh8KZI70uRDtlEZEk1Z8Jt0XTuDdsAuf+FfKahk6kJGcRkSTVj7kvw5TfRePvjoSW3wqbRynBIiJJ2ntffLb5OTJAj/+AI88Jm0cpwyIiSdo71Rvg74Ng4xo4oBucdmvoREohFhFJ0p6LxeCFX8KSd2Cf5vCDh6FBXuhUSiEWEUnSnpvxAMx+FLKy4Zy/QLPi0ImUYiwikqQ9s+BNeOm6aNz3Jjiod8g0SlEWEUnS7qtYBn//MdRWw+FnQc8rQydSirKISJJ2T001PHkRrF0K+30Lzvqjk5Zpj1lEJEm7Z+x1sOANyCuA8x7xYXbaKxYRSdKum/kQTL8fyIKzR0GLQ0InUoqLaxGZPHkyZ555Jm3atCErK4tnnnkmnruTJMXTgmnwwjXR+NQb4NABYfMoLcS1iKxbt46jjjqKu+++O567kSTFW9nn8MSgrRen9romdCKliQbx/MMHDBjAgAE2ZklKadUb4IkLYd1yaNUZzrrHi1NVb+JaRHZXZWUllZWVdevl5eUB00iSiMXg+atg8dvQuAjOf9Qn6qpeJdXFqqWlpRQWFtYtxcXO0CdJQb02At59ArJy4AcPwb4Hhk6kNJNURWTYsGGUlZXVLQsXLgwdSZIy1wfPw/hbovG3fwsHnRw2j9JSUp2aycvLIy/PhyVJUnBL3oWnL4vGx14GPX4aNo/SVlIdEZEkJYGKZfD4BVC9Hg46BU4rDZ1IaSyuR0TWrl3L/Pnz69Y/+eQTZs+eTVFREe3atYvnriVJe6J6Y3SHTPkiaH4I/OBByEmqg+dKM3H92zVjxgxOOeWUuvWhQ4cCcNFFF/HQQw/Fc9eSpN1VWwvP/AwWTYdGzeCHT0DjfUOnUpqLaxHp3bs3sVgsnruQJNWXCbfC+09DdgM472/QvGPoRMoAXiMiSYK3H4Upv4vGZ46EDieFzaOMYRGRpEz3yWR4/spo3OsaOObCsHmUUSwikpTJVsyFJ34EtZvgiLPhlBtCJ1KGsYhIUqZauxwePQc2lkHbY2HgvZDt14ISy79xkpSJqtbBY+fCms+iadsveBwaNgqdShnIIiJJmaZmEzx5ydYH2f3oaWjSInQqZSiLiCRlklgMXvwlzHsZGjSCH/7d23QVlEVEkjLJlP+FmQ8BWfD9B6C4R+hEynAWEUnKFLMfh3/+Jhp/+w447Dth80hYRCQpM8x9GZ4dHI1PuAqO/Y+weaTNLCKSlO4WvgV/vwhiNXDUBdDnptCJpDoWEUlKZ8v/DY/+ADZtgIP7wXf/4FwhSir+bZSkdFW2CB45GzaugbY94NyHIadh6FTSNiwikpSO1q2Cv50N5Z9Di07Rbbq5TUKnkrZjEZGkdLOxHB79Pqz8EPLbRBOW7VMUOpW0QxYRSUon1Rtg9A+jWVP3aQ4/fgaaFYdOJX0ti4gkpYua6mjq9k+nQG4+/Ogp2O/Q0KmknbKISFI6qK2N5gmZ+9LmqdtHQ5tjQqeSvpFFRJJSXSwGL14D7z4B2Q3g3L/CgSeGTiXtEouIJKWyWAxe+S+Y8QCQBd+7DzqdFjqVtMssIpKUyibcCm/cHY2/+wc48pyweaTdZBGRpFQ1+Xcw+Y5oPOAO6DoobB5pD1hEJCkVvXHP1ifp9rsFjrssbB5pD1lEJCnVvPVneHlYNO59ffQ0XSlFWUQkKZVMvz+6QwbgxF/Ayb8Km0faSxYRSUoVM/4CL/wyGp9wFfS5EbKywmaS9pJFRJJSwcyH4P9+EY1LhkDfmy0hSgsWEUlKdjMfhuc3Xwdy/GDoP9wSorRhEZGkZDb9AXj+ymh83OVw2q2WEKUVi4gkJatpf4IXhkbj4wfD6bdbQpR2LCKSlIym/gHGXhuNT7jKIyFKWw1CB5AkfcWUETD+5mh80v+DU26whChtWUQkKVnEYjDhNpj822i99/XQ+9qwmaQ4s4hIUjKIxeDl62HaPdF6nxuh19CwmaQEsIhIUmi1NfB/V8Osv0brA+7w2THKGBYRSQqpphrG/Ce89xRkZcN374ZjLgydSkoYi4gkhVK1Hv5xCcwdC9kN4fv3wxEDQ6eSEsoiIkkhbFgDj58PC96ABo3g3L9Bp/6hU0kJZxGRpESrWAqPfB+WvQd5hfDD0dC+Z+hUUhAWEUlKpNUfw9++B198Ck1bwY+egtZHhk4lBWMRkaREWfIuPHoOrF0G+x4Ig8ZA0UGhU0lBWUQkKRE+mgBPDIKqCmjVOToSkt86dCopOIuIJMXbO0/Asz+H2k1wYC84/1FoVBg6lZQUfOidJMVLLAav3QljLotKyBFnR0dCLCFSHY+ISFI81NbA2OvgrVHReskQ6PcbyPb//6Qvs4hIUn2rXAtPXRpNVEYWnHYblPw8dCopKVlEJKk+lS+Bx8+DJe9EE5V97z5nS5V2wiIiSfVl2fvw6LlQvgj2aQEXjIbiHqFTSUnNIiJJ9WHeq/DkxdHtuS06wQ//DkUdQqeSkp5FRJL2RiwG0+6FV26AWG10e+55f4PG+4ZOJqUEi4gk7alNVfDiNTDr4Wj9mEFwxghokBs2l5RCLCKStCfWr4a//xg+nQJZ2dB/OBz/c8jKCp1MSikWEUnaXcv+BaN/CF98Arn5cM4D0Om00KmklGQRkaTd8a/nYMzlUL0OmrWDC56AVoeHTiWlLIuIJO2K2lqYWAqTfxutdzgJfvAw7FMUNpeU4iwikvRNNpbD05fB3Jei9eN/Hk3XnuN/QqW9lZCHHtxzzz106NCBRo0a0a1bN6ZMmZKI3UrS3lv2LxjVOyohOXkw8E9weqklRKoncS8iTzzxBFdffTU33HADb7/9Nr169WLAgAEsWLAg3ruWpL0z5x9wfx9Y/REUHAA/eQmOviB0KimtZMVisVg8d3DcccfRtWtX7r333rpthx12GAMHDqS0tHSnv1teXk5hYSFlZWUUFBTEM6YkbVVTDa/8N7y5+b9bB/WG7z8ATVoEjSWlit35/o7rscWqqipmzpzJddddt832/v37M3Xq1O3eX1lZSWVlZd16eXl5PONJ0vbKF8OTl8DCadF6r1/CKTdAdk7YXFKaiuupmZUrV1JTU0OrVq222d6qVSuWLl263ftLS0spLCysW4qLi+MZT5K2Nf9V+NOJUQnJK4DzH4M+v7aESHGUkItVs74y02AsFttuG8CwYcMoKyurWxYuXJiIeJIyXc0mGH8LPPJ9WL8KWneByybCt84InUxKe3E9NdOiRQtycnK2O/qxfPny7Y6SAOTl5ZGXlxfPSJK0rfIl8NRP4bPXovXul8Jpt0HDRmFzSRkirkdEcnNz6datG+PGjdtm+7hx4+jZs2c8dy1J32zuy/CnE6ISkpsP5/wFvjPCEiIlUNxvhB86dCiDBg2ie/fulJSUMGrUKBYsWMDll18e711L0o5tqoRxN269K6b1kXDOQ9Di4KCxpEwU9yJy3nnnsWrVKm655RaWLFlC586defHFF2nfvn28dy1J21s5D/5xCSydE60f/3PoexM08LSwFELc5xHZG84jIqnexGIw8yF4+XqoXg/7NIeB9/rUXCkOkmYeEUlKCmtXwHNDYO7YaL3DSfC9UVCwf9hckiwiktLch2OjErJuBeTkQp8bo9Mx2QmZvUDSN7CISEpPlRXwyn9Fp2MAWh4BZ4+C1p2DxpK0LYuIpPTzyWR4djCs2fxwzZIhcOp/e1uulIQsIpLSR9V6ePUmeOu+aL1ZOzjrj9E1IZKSkkVEUnr4bCo8OwRWfxStd7sY+g+HvPygsSTtnEVEUmqrrIBXb4bpf47W89vAWX+Ag/uGzSVpl1hEJKWu+ePh+augbPMDMrv+GPr9Bho3CxpL0q6ziEhKPetWRXfEvPNYtN6sPXx3JBzUO2gsSbvPIiIpdcRi8M7oaHbUDauBLDjucujz35DbJHQ6SXvAIiIpNaycDy/8Iro1F6J5Qc68E4qPDRpL0t6xiEhKbtUb4LU74bXfQ00lNGgMva+N5gbJaRg6naS9ZBGRlLzmvgwv/Qq++DRa73gqnDECijoEjSWp/lhEJCWfLz6DscPgwxei9fw2cPptcPhAyMoKGk1S/bKISEoeVevh9Tvh9btg00bIbgDH/wxOvtaJyaQ0ZRGRFF4sBu89BeN+DeWfR9sO7AXfvgNaHhY2m6S4sohICmvx7Og0zIKp0XphO+j/Gzj8LE/DSBnAIiIpjLJFMP438O7oaL1BY+g1FHpeAQ0bh80mKWEsIpISq7IiuhX3jT9G14EAdDkP+vwaCtuGzSYp4SwikhJjUxXMfBAm3wHrVkTb2p8QPSH3gK5hs0kKxiIiKb5qa+G9f8A/h8Oaz6JtRR2h3y3wrTO8DkTKcBYRSfERi8G8V6LrQJbNibY1bRXditv1x86KKgmwiEiqb7EYfDwRJtwKi6ZH2/IK4ISrojlBfDidpC+xiEiqP5++HhWQz16P1hs0hmP/A078BexTFDabpKRkEZG0d2Ix+HQKTPpt9AqQkwvdfwInDoX8VmHzSUpqFhFJe2bLKZhJv906GVl2QzjmR3DSNd6KK2mXWEQk7Z7aWvjwxWgukM9nRNtycqMLUE+4GpoVB40nKbVYRCTtmppqmPMkvHYnrPww2tagEXS7OLoQtaBNyHSSUpRFRNLObSyHWX+FafdC+aJoW14B9PhpdBdM05Zh80lKaRYRSTtWtgje/BPMfBgqy6NtTfaD438OPS6FRoVh80lKCxYRSdtaNCMqIO+PgdpN0bYWnaBkcPRMGB9IJ6keWUQkRc+B+dczUQH5fObW7Qf2ip6Ge3A/yM4OFk9S+rKISJmsbBHMfCi6BmTtsmhbTi4c+QM49jJoc3TIdJIygEVEyjS1tfDxP2H6X2DuSxCrjbbn7x9d+9H1Ymi6X9CIkjKHRUTKFOVLYPaj8Pbf4ItPt24/sFc0C+phZ/ogOkkJZxGR0llNNcwbF516mffy1qMfeYVw9AVRAdnv0LAZJWU0i4iUjpbOgdmPw5y/w7oVW7e3K4lmQD38LJ+CKykpWESkdFGxDN57Ct55LCoiW+zTIjr6ccyPYb9O4fJJ0g5YRKRUtrEcPng+OvLxyeStp16yG8Khp8PRF8LBfb32Q1LSsohIqaZqHcx9OZpwbN4rsGnj1p+17QFHngtHngP7FIXLKEm7yCIipYLKtTB/HLz/TFRCNm3Y+rMWnbaWj6IOwSJK0p6wiEjJav3qqHR88Dx8NH7bIx/7HghHnA1HDITWXSArK1RKSdorFhEpmaz+BOaOhQ9fhM+mbn3WC0Tl47DvQuezYf+jLR+S0oJFRAqpphoWvhVd6zF3LKz497Y/b3lENNHYYWdCqyMsH5LSjkVESrTyxTD/1WiisY8nQmX51p9l5UD7nnDoAOh0OjTvGCymJCWCRUSKt8oK+PR1+HgCfDQBVn647c/3aQ4d+8Ah/eGQvtB43zA5JSkAi4hU36o3wMI34ZMp8OkU+Hzmttd6ZGVDm2Pg4H5wSL9onJ0TLq8kBWQRkfZWZUVUPD57Axa8AYumQ03Vtu9p1h46ngodT4EOJ3nUQ5I2s4hIuyMWg7JFUfFY+Fb0uvTdrTOabpG/f/RU2w69olfn95CkHbKISDtTWQGLZ0enVz6fGR3tqFiy/fuatYP2J0QPlWt/QnSRqXe4SNI3sohIW1SujR4Wt+SdaFn89ubbaWPbvi8rB/bvAm2PheJjo/JReECQyJKU6iwiyjyxWHQL7bL3ouKx7D1Y+h6sms92pQOgsBgO6AoHdIuWNl0hd5+Ex5akdGQRUfqKxWDdCljxYXRkY/m/YPkH0evGsh3/Tn4b2P8oaHP05teukN8qobElKZNYRJT6qjfCF59GRzRWzYdV82DFXFg5Fzau2fHvZOVED4tr3RlabV5aH2npkKQEs4go+cVisOELWPNZVDhWfxK9fvEJrP4Uyhayw1MqAGTBvu1hv8Og5WHQ8nBo+S1ofgg0bJSwfwRJ0o5ZRBRe9cboTpTyz6HscyhfFL2WLYI1C6Klet3O/4y8guhOleYHQ1FH2K8TtDg0WrdwSFLSimsRufXWW3nhhReYPXs2ubm5rFmzJp67+0a/HzeXNz9ZBUDPji24ss8hdT8bOX4eUz9aCcBxHZrzi36dgmRMG7W10WmRdSuiZe1yWLts87IcKpZuXhZHRzt2RdPW0W2yRR2iJ9Huu/m1+cHQpIW3y0pSCoprEamqquIHP/gBJSUlPPDAA/Hc1S7Jyc5i2serAeper+xzCCPHz2PEuLl17+vZsUWQfEmpthaqKmBjeXSB58ayqDhsXBO9bvgC1q+GDauj1/WromXdSojV7Pp+GjSCgjZQcAAUto2WggOi4tGsfbTukQ1JSjtxLSI333wzAA899FA8d7PLthwB2VI6Royby7SPVzH1o1V17xnar9M2R0qSXiwWTSe+qTJaaja/Vm/YvG1DdOqjen20rXpd9Fq1HqrWRtur1kXjyrXRBF7bLOV8/fUXu6BRs+hoRdNW0LTltq/5raO7VAr2j97nEQ1JyjhJdY1IZWUllZWVdevl5eU7efeeubLPIey7/lOq37w/2vAZ9N38KXRt14yjq5rBS5vfHNvRF3DsS9u3jL/mNRaLpv7e4VITHW2I1UBtTfRQtC+Ptyw1m6C2GmqqN69XR8Xjy0si5ORCo8JoabxvVBwa7wuNm0VPj21cBPsUReN9iqBJy2jcIDcx+SRJKSmpikhpaWndUZR4GnRYNswcu/0PFm9eUllOLuTkRacxGjTe/NoIGjaGhvtES+4+0XpuPuQ2+dLSFPLyIW/za27+1vLhaRFJUhzsdhG56aabvrEsTJ8+ne7du+92mGHDhjF06NC69fLycoqLi3f7z/kmD/8bKjadtd324zoU0ePA5ttu3OHpgqwvbd8y/vIr0aPeyYpes7a85mx+zYbszevZOdu+5jSA7K8sOQ0hu+HmnzWEBnnRtpzcresNGm1ez67Xz0qSpHja7SIyZMgQzj///J2+58ADD9yjMHl5eeTl5e3R7+6qkePnMeK1jcB5APTs2HzrNSLzYOiBKXaNiCRJKWy3i0iLFi1o0SI17yr56t0xWy5M/fL2La+WEUmS4i+u14gsWLCA1atXs2DBAmpqapg9ezYABx98ME2bNo3nrneopjbG8QcVAdvOI7Lldcs8IjW1e3GXiCRJ2mVZsdgObw2pFxdffDEPP/zwdtsnTJhA7969v/H3y8vLKSwspKysjIKCgjgklCRJ9W13vr/jWkT2lkVEkqTUszvf395iIUmSgrGISJKkYCwikiQpGIuIJEkKxiIiSZKCsYhIkqRgLCKSJCkYi4gkSQrGIiJJkoKJ67Nm9taWSV/Ly8sDJ5EkSbtqy/f2rkzentRFpKKiAoDi4uLASSRJ0u6qqKigsLBwp+9J6mfN1NbWsnjxYvLz88nKyqrXP7u8vJzi4mIWLlzoc2ziyM85MfycE8PPOTH8nBMnXp91LBajoqKCNm3akJ2986tAkvqISHZ2Nm3bto3rPgoKCvyLngB+zonh55wYfs6J4eecOPH4rL/pSMgWXqwqSZKCsYhIkqRgMraI5OXlceONN5KXlxc6Slrzc04MP+fE8HNODD/nxEmGzzqpL1aVJEnpLWOPiEiSpPAsIpIkKRiLiCRJCsYiIkmSgsnIInLPPffQoUMHGjVqRLdu3ZgyZUroSGln8uTJnHnmmbRp04asrCyeeeaZ0JHSUmlpKT169CA/P5+WLVsycOBAPvzww9Cx0s69995Lly5d6iZ9Kikp4aWXXgodK+2VlpaSlZXF1VdfHTpKWrnpppvIysraZmndunWwPBlXRJ544gmuvvpqbrjhBt5++2169erFgAEDWLBgQehoaWXdunUcddRR3H333aGjpLVJkyYxePBgpk2bxrhx49i0aRP9+/dn3bp1oaOllbZt23L77bczY8YMZsyYwamnnspZZ53F+++/Hzpa2po+fTqjRo2iS5cuoaOkpSOOOIIlS5bULXPmzAmWJeNu3z3uuOPo2rUr9957b922ww47jIEDB1JaWhowWfrKyspizJgxDBw4MHSUtLdixQpatmzJpEmTOOmkk0LHSWtFRUXccccdXHrppaGjpJ21a9fStWtX7rnnHoYPH87RRx/NnXfeGTpW2rjpppt45plnmD17dugoQIYdEamqqmLmzJn0799/m+39+/dn6tSpgVJJ9aesrAyIviQVHzU1NYwePZp169ZRUlISOk5aGjx4MGeccQZ9+/YNHSVtzZs3jzZt2tChQwfOP/98Pv7442BZkvqhd/Vt5cqV1NTU0KpVq222t2rViqVLlwZKJdWPWCzG0KFDOfHEE+ncuXPoOGlnzpw5lJSUsHHjRpo2bcqYMWM4/PDDQ8dKO6NHj2bWrFlMnz49dJS0ddxxx/HXv/6VTp06sWzZMoYPH07Pnj15//33ad68ecLzZFQR2SIrK2ub9Vgstt02KdUMGTKEd999l9deey10lLR06KGHMnv2bNasWcNTTz3FRRddxKRJkywj9WjhwoVcddVVvPLKKzRq1Ch0nLQ1YMCAuvGRRx5JSUkJHTt25OGHH2bo0KEJz5NRRaRFixbk5ORsd/Rj+fLl2x0lkVLJFVdcwXPPPcfkyZNp27Zt6DhpKTc3l4MPPhiA7t27M336dO666y7uu+++wMnSx8yZM1m+fDndunWr21ZTU8PkyZO5++67qaysJCcnJ2DC9NSkSROOPPJI5s2bF2T/GXWNSG5uLt26dWPcuHHbbB83bhw9e/YMlErac7FYjCFDhvD000/zz3/+kw4dOoSOlDFisRiVlZWhY6SVPn36MGfOHGbPnl23dO/enQsvvJDZs2dbQuKksrKSDz74gP333z/I/jPqiAjA0KFDGTRoEN27d6ekpIRRo0axYMECLr/88tDR0sratWuZP39+3fonn3zC7NmzKSoqol27dgGTpZfBgwfz2GOP8eyzz5Kfn193tK+wsJDGjRsHTpc+rr/+egYMGEBxcTEVFRWMHj2aiRMnMnbs2NDR0kp+fv521zc1adKE5s2be91TPbrmmms488wzadeuHcuXL2f48OGUl5dz0UUXBcmTcUXkvPPOY9WqVdxyyy0sWbKEzp078+KLL9K+ffvQ0dLKjBkzOOWUU+rWt5x3vOiii3jooYcCpUo/W25D79279zbbH3zwQS6++OLEB0pTy5YtY9CgQSxZsoTCwkK6dOnC2LFj6devX+ho0m5btGgRF1xwAStXrmS//fbj+OOPZ9q0acG+BzNuHhFJkpQ8MuoaEUmSlFwsIpIkKRiLiCRJCsYiIkmSgrGISJKkYCwikiQpGIuIJEkKxiIiSZKCsYhIkqRgLCKSJCkYi4gkSQrGIiJJkoL5/2lDqFkJeA+KAAAAAElFTkSuQmCC",
"text/plain": [
"<Figure size 640x480 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Vanishing cells?\n",
"x = [0., 1.e-2, 2.e-2, 5.] # Positiong of data, x[i]\n",
"y = [-1., -1., -1., 4.] # Values of data, y[i]\n",
"q = numpy.linspace(x[0],x[3],100) # Sample interval between x[1]..x[2]\n",
"plt.plot(x,y,'x');\n",
"plt.plot(q, lagrange_interp(x,y,q));"
]
}
],
"metadata": {
"gist": {
"data": {
"description": "Interpolation with Lagrange cubic polynomial",
"public": true
},
"id": ""
},
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.11.5"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment