{ "cells": [ { "cell_type": "markdown", "id": "organizational-budget", "metadata": {}, "source": [ "# Chebyshev Interpolation Example" ] }, { "cell_type": "markdown", "id": "flush-petroleum", "metadata": {}, "source": [ "This notebook demonstrates how to use Chebyshev expansions for interpolation when the samples fall on the proper grid. The samples have to be on this grid, however." ] }, { "cell_type": "code", "execution_count": 1, "id": "sharing-cradle", "metadata": {}, "outputs": [], "source": [ "from orthopoly.chebyshev import *\n", "from numpy import *\n", "import matplotlib.pyplot as plt\n", "%matplotlib inline" ] }, { "cell_type": "markdown", "id": "radical-future", "metadata": {}, "source": [ "First, I define a function to sample from and interpolate. In practice, the samples could measurements or anything else. The function here is simple and fast to evaluate, so you wouldn't need an interpolating function, but just imagine the function is cumbersome and slow or the measurements are very costly to get! We'll look at the range from $x=-6$ to $x=3$." ] }, { "cell_type": "code", "execution_count": 2, "id": "informed-evening", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEICAYAAABVv+9nAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAmrklEQVR4nO3deXxU9b3/8dcnG0kgCQkQAiQkhCCLsi8igguIUrTue91FtHXrdr1X2yq917a3vS6trf4sYtWqdaPgvivKoiKrrAoYQMIOAcKShCzf3x8ZClhIRsjJmTnzfj4e83Amk5nz5vuI75x855zvMeccIiISPHF+BxAREW+o4EVEAkoFLyISUCp4EZGAUsGLiASUCl5EJKA8LXgza2lmE83sSzNbamYneLk9ERHZz+s9+D8BbzvnugG9gaX1ffOTTz7pgCO6bd68+YhfG7SbxkLjofGIqbE4rIT6njwaZpYBnARcA+Cc2wvsre81q1atOuLtVVZWHvFrg0ZjcTCNx8E0HvsFfSy83IPvBGwGnjCzeWY2wcyae7g9ERE5gGd78KH37gfc6pybaWZ/Av4L+NWB32RmY4GxACNHjqSkpOSINlZaWnp0aQNEY3EwjcfBNB77BWEscnNzD/uclwVfApQ452aGHk+kruAP4pwbD4wHGDdunKsvbEOO5rVBo7E4mMbjYBqP/YI8Fp5N0TjnNgBrzKxr6EsjgCVebU9ERA7m5R48wK3As2aWBBQD13q8PRERCfG04J1z84EBXm5DREQOTWeyiogElNdTNFHHOUfp7r2s3V7Oropqdu+tYc/eavZW15IYHxe6GRkpieRkJNM2PZnkxHi/Y4uI/JuYLPidFVWsKS1nzbY9rCndQ8m2ctaU7mHNtrr7e/bWAGAGqYnxpDZLICk+jqqa2tDNsauy+l/v1zI1kS7ZLejeLp1euS0ZXJhFbmaqX/88EREgoAW/u7Ka9TsqKNm2hzXbyikJlfe+Ut++pwqA5knx5GWlkpuZSl5WKkOKWpOXmUJeViodMlNIa5aAmR1yG5XVNWzeWcnGsgrWbq9g+cadLF1fxusL1lO6ey/5rVI5vUdbvtezHX1yWxIXd+j3ERHxSiAK/rdvLmXeyk1sr/yaDTsq2Bnau06KjyM3M4UOodLu2aEleVkp5IUKPTM18bAF3pBmCfHkZtb9cuifv//rtbWOZZt2Mm3ZFt5evIHHpq2koFUqlx/fkQv755HVPKkx/skiIg0KRMEnJ8ZzTJsUjsnNpm16Mu0yksnJSKZNi2ZNvuccF2d0y0mnW046N5xUyPod5bw0u4QnZqzi/neXcfnxHbnp5M60TU9u0lwiEnsCUfA/HXkMJSUlEXlGWruMFG4b0YUfndKZtxZt4C8fruDZmd9wzZACbhleRHpyot8RRSSgdJhkE0mIj+P7vdvz1u3DePDiPry5cD3D7/uIF2evwbl6V/wUETkiKvgmFhdnnNmrHe//9GSuPqGAX728iKufmMW67eV+RxORgFHB+yQ5MZ5bR3ThzduHsbOiijMenMprX6zzO5aIBIgK3med27Rg4k1DuOmUztz+/DzueWURldU1fscSkQBQwUeA+Djj5lOLeHbMYN5YuIFLx3/G5p3BvtKMiHhPBR9BTujcijdvG0p1jeO8R2awfONOvyOJSBRTwUeY7PRkXrhxMN3bpXP+//uEWaui/4ozIuIPFXwESk1K4NEr+nN27/Zc9fjnzFixxe9IIhKFVPARKj7OuPfc47hicEeufXIWH3650e9IIhJlVPARzMy4a3R3bjypkJuemas9eRH5TlTwEc7M+OnIY7j6hHzGPDWbOas1Jy8i4VHBR4F9e/Ln9evANU/M4ssNZX5HEpEooIKPEmbGveccx9Ci1lz3xCw2llX4HUlEIpwKPorExRkPXtKHnIxkrntyFrsPuKqUiMi3qeCjTHJiPI9dNYBdldXc/vx8amu1EqWIHJoKPgq1atGMCVcN4NOvt/DQh8v9jiMiEUoFH6W6tE3jgUv68NAHy3lviY6RF5F/p4KPYmccm8Mtpxbxkxfms3rrbr/jiEiEUcFHudtPO4ZeuRnc9tw89lbX+h1HRCKICj7KxYeOrFmzrZz73/3K7zgiEkFU8AHQNj2Z+y/qzfhpxXy8bLPfcUQkQqjgA+LUbtlcf2InfvbifLburvI7johEABV8gNwxqhvtMlL43YdrcE7Hx4vEOk8L3sxWmdlCM5tvZrO93JZAUkIcD17Sm7klu5g4p8TvOCLis6bYgz/VOdfHOTegCbYV84qy07j++Bz++/UlbNih9WpEYpmmaALokj5tKGzTgrsmL9RUjUgM87rgHfCumc0xs7Eeb0tCEuKM+y7sxfTlW5g0d63fcUTEJwkev/9Q59xaM8sG3jOzL51zUw/8hlDxjwUYOXIkJSVHNndcWqoLYexTWlpKVhZcOyibe15dROfmVbRukeh3LN/oZ+NgGo/9gjAWubm5h33O04J3zq0N/XeTmU0GBgFTv/U944HxAOPGjXP1hW3I0bw2aHJzc/mPdu2Z8c0MJszdxl8u7+d3JF/pZ+NgGo/9gjwWnk3RmFlzM0vbdx84HVjk1fbk3yXEx3HvuT15Y+F6pi3XCVAiscbLOfi2wHQz+wL4HHjDOfe2h9uTQ+iT15LLB3Xk7lcWU1ld43ccEWlCnhW8c67YOdc7dDvWOfcbr7Yl9bvjjG6UlVfx14+L/Y4iIk1Ih0nGgIzURO4a3Z2/TFmhZYVFYogKPkac368DfXJbcs+ri3VsvEiMUMHHCDPjf849jmnLt/DRV/rAVSQWqOBjSNecNC4f1JF731hCVY0uDiISdCr4GPOTkcewaWcl/5j5jd9RRMRjKvgYk9U8iduGd+HB95exY4/WjRcJMhV8DLpqSD4ZKYn86YPlfkcREQ+p4GNQs4R47hrdnb9/uorizbv8jiMiHlHBx6jTe7RlQEEm//vWl35HERGPqOBjlJlx5/e68+6SjcxZvc3vOCLiARV8DOud15LRPXP4/dtf6uQnkQBSwce4n53elTmrt/HRMp38JBI0KvgY17lNCy4ekMsf3v6K2lrtxYsEiQpeuH3EMRRv3sVrC9b5HUVEGpEKXsjJSOaaEwu4/91l7K3WEgYiQaGCFwB+dHIR2/bs5cXZa/yOIiKNRAUvQN2a8WOGFvLIlBW68pNIQKjg5V+uHVrArspqXpxd4ncUEWkEKnj5l/TkRMYM0168SFCo4OUg15xYwG7txYsEQoMFb2bnm9lyM9thZmVmttPMypoinDS99OREbtBevEgghLMH/wfgbOdchnMu3TmX5pxL9zqY+OeaEwvYs7eGF2fpiBqRaBZOwW90zi31PIlEjLTkRG4Y1omHp3ytvXiRKBZOwc82sxfM7LLQdM35Zna+58nEV1cPKaC8qoYXtBcvErXCKfh0YA9wOvD90O0sL0OJ//btxT8y5WsqqrQXLxKNEhr6BufctU0RRCLP1UMKeGzaSibOKeGKwfl+xxGR7yico2hyzWyymW0K3f5pZrlNEU78lZacyNVDCvjr1K+prtEaNSLRJpwpmieAV4H2odtroa9JDLh2SAFbd+3l9QXr/Y4iIt9ROAXfxjn3hHOuOnR7EmjjcS6JEJnNk7hsUEce+WiF1osXiTLhFPxWM7vCzOJDtyuArV4Hk8hxw7BCVm7ZzQdfbvI7ioh8B+EU/HXAxcAGYD1wIRD2B6+hXwrzzOz1I4sofsvJSObC/rk8PGWFrt0qEkUaLHjn3Grn3NnOuTbOuWzn3LnOuW++wzZuB3SiVJS78aTOLCjZzqfF+uNNJFoc9jBJM7vDOfcHM/sz8G+7bc652xp689DRNmcCvwF+ejRBxV8FrZtzZq/2PDLla4Z0bu13HBEJQ33Hwe/b6559FO//R+AOIO1w32BmY4GxACNHjqSk5MhWMSwtLT2i1wWRV2NxfvcWXPv8Mt6bs4zubVM92YYX9LNxMI3HfkEYi9zcwx+1ftiCd869Frq7xzn30oHPmdlFDW3UzM4CNjnn5pjZKfVsZzwwHmDcuHGuvrANOZrXBo0XY5GbC8Pnb2fikjL+2v+YRn9/L+ln42Aaj/2CPBbhfMh6Z5hf+7YTgbPNbBXwPDDczJ75DtkkAv3olM68s3gjKzbt9DuKiDSgvjn47wGjgQ5m9tABT6UD1Q29sXPuTkK/CEJ78D93zl1xNGHFfwMKsuifn8mEaSv53wt6+R1HROpR3x78Ourm3yuAOQfcXgXO8D6aRKqxJxUyae5aNu2s8DuKiNSjvjn4L4AvzGwysNs5VwN1x7UDzb7LRpxzHwEfHXlMiSSndW9Lh8wU/v7Jan5+Rle/44jIYYQzB/8ukHLA4xTgfW/iSDSIjzPGDOvE05+tZndlg7N1IuKTcAo+2Tm3a9+D0P3oOUZOPHFBv1wS4oyXZuuCICKRKpyC321m/fY9MLP+QLl3kSQaJCfGc9UJBUyYvlJLCYtEqHAK/sfAS2Y2zcymAy8At3iaSqLClSfks2VXJW8t2uB3FBE5hHCu6DTLzLoB+z5N+8o5V+VtLIkGWc2TuKh/HuOnFnNWr3aYmd+RROQA4ezBAwwEegH9gMvM7CrvIkk0GTOsE4vW7eCz4ug/5VskaBrcgzezp4HOwHxg39WXHfB372JJtMhv1ZxRx+bw2LRiTujcyu84InKABgseGAD0cFoIXA5j7EmFnPfIJyzfuJMubQ+7rpyINLFwpmgWATleB5Ho1bdjJoMKsnhsWrHfUUTkAOEUfGtgiZm9Y2av7rt5HUyiyw0nFfLyvHVsKtPyBSKRIpwpmnFeh5DoN6JbNrlZKTz16Sr+44xufscREcI7TPLjpggi0S0uzrh+aCf+8PZX3HxqEalJ4ew7iIiXGpyiMbOdZlYWulWYWY2ZlTVFOIkuF/TLJT7OmDjnyK7KJSKNK5yLbqc559Kdc+nULTR2AfCI58kk6iQnxnPF4Hwen76SmloddCXit3BPdALA1XkZrQcvh3Hl4HzW76jgvSUb/Y4iEvPCOdHp/AMexlF3XLwOlZBDapPWjPP6dODx6cWMOk5H14r4KZw9+O8fcDsD2Amc42UoiW5jhnVi1qptzF+z3e8oIjGtvmuynu+cm+Scu9bMMp1z25oymESvLm3TOKVrGx6bVszDl/dr+AUi4on69uB/ecD9D7wOIsFyw7BC3lq4njWle/yOIhKz6it4O8x9kQYN6dyKrjnpPDFjld9RRGJWfQWfYmZ9Q1dwSg7d77fv1lQBJTqZGTcM68QLs75hR7kuHyDih/oKfj3wAHAfsCF0//7Q7T7vo0m0O6tXe1okJ/D859/4HUUkJh32Q1bn3KlNGUSCJykhjquHFPDkJ6u4bmgnEuO/02kXInKU9H+ceOoHg/LZUV7FmwvX+x1FJOao4MVTGamJXDwgj8emFaNrxog0LRW8eO66EzuxZF2Zrtsq0sTCWU3SzOwKM7s79LijmQ3yPpoERcdWqZxxbA4TdMUnkSYVzh78I8AJwGWhxzuBhz1LJIE0ZlghH3y5iRWbdvkdRSRmhFPwxzvnbia0wFhoyYIkT1NJ4PTPz6Rvx5Y8Pn2l31FEYkY4BV9lZvGAAzCzNkBtQy8ys2Qz+9zMvjCzxWb266PMKlHuhmGFTJpbwtZdlX5HEYkJ4RT8Q8BkINvMfgNMB34bxusqgeHOud5AH2CUmQ0+0qAS/c44Nofs9GY8/dlqv6OIxIRwruj0LHAH8Dvqzm491zn3Uhivc865fROuiaGbjpOLYfFxxnUnduLpT1dTUVXjdxyRwDtswZtZ1r4bsAl4DvgHsDH0tQaZWbyZzQ+9/j3n3MxGyCxR7OIBeVTV1PLyvLV+RxEJvPqu6DSHuj3uQ60k6YDCht7cOVcD9DGzlsBkMzvOObfowO8xs7HAWICRI0dSUnJkF2wuLdUx1vtE+lh8v0cmj05ZxpB2Rpx5v1BppI9HU9N47BeEscjNzT3sc/WtRdOpsQI457ab2RRgFLDoW8+NB8YDjBs3ztUXtiFH89qgieSxuPWM1rzw+w/5enczTu2W3STbjOTx8IPGY78gj0VYZ7Ka2flm9oCZ3W9m54b5mjahPXfMLAUYCXx5pEElOHIykvl+7/Y8phOfRDwVzpmsjwA3AQup2/u+yczCOdGpHTDFzBYAs6ibg3/9aMJKcIwZ1olPvt7K4nU7/I4iElj1zcHvMxzo7kIrRZnZU8Dihl7knFsA9D26eBJUx7bPYEjnVkyYtpIHL+njdxyRQApnimYF0PGAx3mhr4kclRuGFfLaF+tYv6Pc7ygigRROwacBS83sIzP7CFgCpJvZq2b2qqfpJNBOPqYNBa2b89QnOvFJxAvhTNHc7XkKiUlxccaYoZ347ZtLuXV4Ec2bhfPjKCLhCudM1o+dcx8D86j7oHUhsPCAr4scsXP7diAxPo4XZ6/xO4pI4IRzFM1YM9sALABmU3cC1Gyvg0lsSE6M58oT8vnbjJXU1GolC5HGFM4c/H8AxznnCpxzhc65Ts65Bs9iFQnXlYPz2VRWyTuLN/gdRSRQwin4r4E9XgeR2NWqRTPO75erE59EGlk4n2rdCXxiZjOpWwIYAOfcbZ6lkphz/dBOnPbAx8xZXUr//LDWshORBoSzB/9X4EPgM+rm3/fdRBpNUXYLRnTL5rGpuuKTSGMJZw8+0Tn3U8+TSMwbM6yQyyd8xuqtu8lv1dzvOCJRL5w9+LdCR9K0+9Ya8SKNanBhFse2T+eJGav8jiISCOEU/GWE5uHZPz2jwySl0ZkZY4YW8uLsNezYU+V3HJGoF86JTp0OcdNhkuKJM3u1IyMlkWc/1/IFIkcrrHPDzew4oAeQvO9rzrm/exVKYldifBzXDCngbzNWMmZoIUkJYV2yQEQOIZwzWe8B/hy6nQr8ATjb41wSwy4d1JFdFdW89sU6v6OIRLVwdo8uBEYAG5xz1wK9gQxPU0lMy0hJ5JKBHXlsWjGhyxCIyBEIp+DLnXO1QLWZpQObqFsTXsQz155YwLKNO5mxYqvfUUSiVjgFPzt0bdXHqDuCZi7wqZehRPKyUvlez3ZMmK7lC0SOVDhH0fzIObfdOfcodRfOvjo0VSPiqTFDO/HRV5tZvnGn31FEolI4H7Jev+++c24VsDj0wauIp/p2zGRAfiYTpmn5ApEjEc4UzQgzezN0Juux1K1Jk+ZxLhGgbvmCyfPWsmlnhd9RRKJOOFM0lwNPUXclpzeBHzvnfu51MBGAkT3akpuZwuPTtRcv8l2FM0XTBbgd+CewGrjSzFK9DiYCEB9n/OjUIp75dDXbdu/1O45IVAlniuY14FfOuRuBk4HlwCxPU4kc4Jw+7clsnsQTn6zyO4pIVAmn4Ac55z4AcHXuB87zNpbIfonxcfzwlM48OWMlZRVahEwkXIcteDO7A8A5V2ZmF33r6Wu8DCXybRf2zyU1KYGnP9UiZCLhqm8P/tID7t/5redGeZBF5LCaJcRz48mFPD59JXv2VvsdRyQq1Ffwdpj7h3os4rlLB3YkzuAfM7/xO4pIVKiv4N1h7h/qsYjnUpLiGTOskPFTi6moqvE7jkjEq6/ge5tZmZntBHqF7u973LOJ8okc5IrB+eytqeWl2Wv8jiIS8Q5b8M65eOdcunMuzTmXELq/73FiQ29sZnlmNsXMlpjZYjO7vXGjSyxq0SyB607sxKMfF7O3utbvOCIRzcvL5VQDP3PO9QAGAzebWQ8Ptycx4uohBZRVVPHSHO3Fi9THs4J3zq13zs0N3d8JLAU6eLU9iR0ZKYmMGVrIXz5cobl4kXo0yQUvzawA6AvMbIrtSfBdN7SA8qoanv9cR9SIHE5YF90+GmbWgrp1bH7snCs7xPNjgbEAI0eOpKSk5Ii2U1paejQxAyVWxuKS3q348wfLGNo+nuTEw++rxMp4hEvjsV8QxiI3N/ewz3la8GaWSF25P+ucm3So73HOjQfGA4wbN87VF7YhR/PaoImFsbitTQ4TF0xhyppqbjipsN7vjYXx+C40HvsFeSw8m6IxMwMeB5Y65x7wajsSu5o3S+CHp3Tm/338NbsrdXaryLd5OQd/InAlMNzM5oduoz3cnsSgKwbnkxBnPKmVJkX+jWdTNM656WhJA/FYcmI8twwv4r53vuIHx3ekZWqS35FEIkaTHEUj4qVLB3Yks3kSj3z0td9RRCKKCl6iXlJCHD8/vStPfrKKtdvL/Y4jEjFU8BIIZ/ZsR7ecNB54d5nfUUQihgpeAiEuzvivUd2YNK+Epev/7XQLkZikgpfAGFLUmpO6tOH3b3/pdxSRiKCCl0D5z1Hd+HjZZj75eovfUUR8p4KXQOnRPp0L+uXyP68vpaZW16WR2KaCl8C544yufLN1N8/P0kJkEttU8BI42enJ3DK8C/e98xU79lT5HUcagXOOiqoayiqq2Lqrkg07KijdvZeqGl30pT6eryYp4ofrhhbw/Kxv+OMHy7i+b4bfcSRM5XtrWLJ+BwtKdrB0fRlrt5ezfnsF63dUUH6Ytf9TEuNpmZpIbmYK+a2ak5+VSlF2C3rltaR9RjJ1y2LFJhW8BFKzhHh+eWYPbnpmDsPzjyHACwZGteqaWuav2c6Urzbx0VebWbq+jFoHha2b06N9Or1zWzLq2GTaZaSQ1SKJpPg4khLiSIyPo7K6hrLyanZWVLFtTxVrSveweutuPvxqE3+dWsyuympat2hG79wMBnXKYliXNnTLSSMuLnYKXwUvgXVa92yGdG7FQ9PWMrRn55jek4skzjlmrdrGxDlreGfxRsoqquid25KRPdryyzN7cGyHdNKTG7zsc71qax3FW3Yxf80O5q/Zxouz1/C7t76kdYskhha15vRjczila5tG+hdFLhW8BJaZcfdZPRj1p6m8tWgDo3u28ztSTNu8s5LnPv+GiXNKWLNtD0OLWnP3WT04pWsbWrVo1qjbioszirLTKMpO48L+dX++rdtezvQVW/j4q838/KUvqHWOQXlpnD/QGN49+6h/qUQiFbwEWpe2aVzeN5txry5maJfWgfyfONKt2LSLx6cX88+5a2mXkcylg/I4r28H2mWkNGmO9i1TuHhAHhcPyKOiqoapyzbzz8+L+dUri7hjYi2n9cjmvL65nHxMG5ISgnH8iQpeAu+agW2ZumoXf3j7S+49t6ffcWLG0vVlPPDeMt5bspH++Zn8+bK+nNa9LfERMAeenBjP6cfm0COjmuyc9kxbvplJ89Zy8z/m0jwpnrN7t+e8frn0zs2I6qk9FbwEXrOEOH5zbk+u/NtMzuvbgf75WX5HCrRVW3bz4PvLePWLdQzr0oZ//vCEiB7zpIQ4RnRvy4jubdlRXsVbC9czae5annp4BoVtmnN+3w6c27cDuZmpfkf9zlTwEhOGdmnNuX06cOekhbx+67DA/AkeSbbv2cuD7y3jmZnf0CevJc/dMJjBha38jvWdZKQkcumgjlw6qCNrSvcwed5aJs4p4b53lzGoIIvv92nP6ONyGv0zA6+o4CVm/PLM7ox44GMem1bMzacW+R0nMGpqHS/MWsP/vfMlGSmJ/PWK/ozonh3VUxsAeVmp3DaiC7cOL2Lemu28On8dD32wnHGvLmZI51ac3bs9px+bQ0ZK5H6uo4KXmNGqRTN+eWYP7pq0kOHdsuneLt3vSFFvYckO7py8gOLNu7l1eBeuG1pAs4R4v2M1KjOjX8dM+nXM5Fdn9WBm8VZeW7COe99Yyi8mL+LEolac0jWbU7q2Ib9Vc7/jHkQFLzHlgn4deGfxBn7ywnxeueXEwJVRUynfW8Mf31/GhOkrGXVcDhOuGkhORrLfsTwXH2cMKWrNkKLW/Prs45i+YjPvLdnEXz/+mnteXUyn1s0ZWJBJ77yW9MlrSde2aSTE+zcdqIKXmGJm/O78noz641T++P5y/nNUN78jRZ1Pv97KnZMWUF5Vw6NX9Gdkj7Z+R/JFUkIcw7u1ZXi3tjjnWL5pF1OXbWbuN9t4ZMrXrN1eTlJ8HLlZKeRnpdKqRTPSkhNISYxnb3UtldW1bCyroGRbOe1bJjPh6oGNnlEFLzGndYtm/O78Xtz49GxGdMtmQEHkHuERSXZVVvObN5by/KxvuHRgR+4c3U3nFYSYGce0TeOYtmn/+tqmsgoWry8LLaGwh22797KmtJyKqhqSEuJIio+jXUYyAwoy6XLA6xqTCl5i0sgebbmwfy4/ffEL3rx9GC2a6X+F+sxZXcpPXvgCh+MfYwZzQufoOjrGD9npyWSn+zttpWPFJGb96qweOBx3v7wI53RxkEOpqqnlgXe/4qJHP+X4Tlm8dftJKvcoot0WiVlpyYn8+bJ+XPToJxxfmMUlAzv6HSmiFG/exU9emM83pXt45Af9GHWc1vKJNtqDl5jWJ68lvxjdnbtfWczS9WV+x4kIzjme+Ww1Zz40nYzUJN758Ukq9yilPXiJeVcPKeDzVaX88Jk5vHLzUDJSY/eDw807K/nPfy5gxoot3DW6O1edkB/1JyzFMu3BS8wzM35/QS+SEuK45bm5VMfoZeDeX7KRUX+cysayCt64bShXDylQuUc5FbwIdfPxj101gIVrd/DbN7/0O06T2rO3mrsmL2Ts07O5aEAek390IkXZ3hy2J01LUzQiIfmtmvPI5f248m+fU9imOVcMzvc7kucWlGznx8/Pp7K6luduGMzxUbY4mNTPsz14M/ubmW0ys0VebUOksQ0pas3vzuvJ3a8s4u1FG/yO45maWsfDU1Zw/iOf0DM3gzdvH6ZyDyAv9+CfBP4C/N3DbYg0uosH5rFpZwW3PT+Pp68bFLji21C2l5++/hlL15dx30W9ObdvB78jiUc824N3zk0FSr16fxEv3XxqEZcNzOP6p2YzZ/U2v+M0CuccL89by9XPfwXAWz8epnIPON/n4M1sLDAWYOTIkZSUlBzR+5SW6nfJPhqLgx3peFzfL4PSHTu5csJnPHB2Ice1i6ylYL+LrburuO/jEj5ZVcblPTMYc2Iu7C6lZHds/6wE4f+V3Nzcwz7ne8E758YD4wHGjRvn6gvbkKN5bdBoLA52pOPxpyty+cXLC/n566uYcPWAqLtCkXOOyfPW8uvXlpOXlcLrtw4jraZMPx8HCPJY+F7wIpEsLs74zbk9aZYQz1WPf84Dl/TmrF7t/Y4Vlg07KvjF5IVMW76F20/rwo0nFZIQH0dJic7YjRUqeJEGxMUZ93y/B+0ykrntuXms217ODcMKI/YkoOqaWp76dDUPvreMztkteOO2oZ4tRyuRzbOCN7PngFOA1mZWAtzjnHvcq+2JeMnMuPHkzuRkJHPHxAV8UbKD31/QK+KWGZ5ZvJW7X1nMuh3l3DGqKz84Pp/4uMj8RSTe8+yn0zl3mVfvLeKXc/p0oEt2Gjc9M4dzH57BHy/pw3EdMvyORfHmXdz/3jLeXLieC/rl8uwNx9O6RTO/Y4nPtFSByHfUo306r90ylK45aZz78AweeG8Ze6v9Wb9mY1kFd01eyMgHp7KprIKJNw3hvot6q9wF0By8yBHJSE3k4cv78fqCdfzq5UW8vmAdvxjdneHdsptkbn7llt08Pr2YiXNKKGjVnMeu6s+pXZtm2xI9VPAiR+GsXu0Z0rk1f3p/GWOfnsPAgkxuPrWIoUWtG71sa2sdM1eW8rcZK3l/6UZ6dcjg/y7szeie7TTPLoekghc5SlnNk/j1Ocdx5QkF/OmD5VzzxCy6ZLfgB8d3ZHTPdrQ6iukS5xxL1+/k7UXrmTRvLeu2lzOie1teGHsCAwsytccu9VLBizSSouwW/PmyvtxxRlee+mQVD324gnGvLWFwYRYnFLZiYEEW3XLS672gSEVVDcWbd/NFyXbmrt7G9BVbWL+jgm45aVw5OJ9z+nQgJ8PfCzlL9FDBizSyvKxUfnlWD+4c3Z2ZxVt5d8lG3lq0gfvfW4ZzkJmaSE5GCi2axZOcGE9ldS179lazeWclG8sqAWiXkUzfji354SmdGd4tm9zMVJ//VRKNVPAiHomPM4YUtWZIUWsAdlZUsWrLHlZu3c3GHRXs2VtDeVUNyYlxpCbFk9W8GfmtUslvlUp2mvbS5eip4EWaSFpyIj1zM+iZ6/9x8xIbdBy8iEhAqeBFRAJKBS8iElAqeBGRgFLBi4gElApeRCSgVPAiIgGlghcRCShzzvmd4V/MbAJQcoQv7w/MacQ40UxjcTCNx8E0HvsFYSxWOeeePNQTEVXwR8PMZjvnBvidIxJoLA6m8TiYxmO/oI+FpmhERAJKBS8iElBBKvjxfgeIIBqLg2k8Dqbx2C/QYxGYOXgRETlYkPbgRUTkAIEreDO71cy+NLPFZvYHv/P4zcx+ZmbOzFr7ncVPZvZ/oZ+LBWY22cxa+p2pqZnZKDP7ysxWmNl/+Z3HT2aWZ2ZTzGxJqCtu9zuTFwJV8GZ2KnAO0Ns5dyxwn8+RfGVmecDpwDd+Z4kA7wHHOed6AcuAO33O06TMLB54GPge0AO4zMx6+JvKV9XAz5xzPYDBwM1BHI9AFTzwQ+B/nXOVAM65TT7n8duDwB1AzH/Q4px71zlXHXr4GZDrZx4fDAJWOOeKnXN7geep2xmKSc659c65uaH7O4GlQAd/UzW+oBX8McAwM5tpZh+b2UC/A/nFzM4B1jrnvvA7SwS6DnjL7xBNrAOw5oDHJQSw0I6EmRUAfYGZPkdpdFF3TVYzex/IOcRTv6Du35NF3Z9cA4EXzazQBfRQoQbG4i7qpmdiRn3j4Zx7JfQ9v6Duz/NnmzKbRCYzawH8E/ixc67M7zyNLeoK3jl32uGeM7MfApNChf65mdUCrYHNTZWvKR1uLMysJ9AJ+MLMoG46Yq6ZDXLObWjCiE2qvp8NADO7BjgLGBHUX/r1WAvkHfA4N/S1mGVmidSV+7POuUl+5/FC0KZoXgZOBTCzY4AkYIufgfzgnFvonMt2zhU45wqo+3O8X5DLvSFmNoq6zyPOds7t8TuPD2YBXcysk5klAZcCr/qcyTdWt+fzOLDUOfeA33m8ErSC/xtQaGaLqPsQ6eoY3FOTQ/sLkAa8Z2bzzexRvwM1pdAHzLcA71D3geKLzrnF/qby1YnAlcDw0M/DfDMb7XeoxqYzWUVEAipoe/AiIhKighcRCSgVvIhIQKngRUQCSgUvIhJQKngRkYBSwYuIBJQKXqQeoTXDR4bu32tmf/Y7k0i4om4tGpEmdg/w32aWTd2Kg2f7nEckbDqTVaQBZvYx0AI4JbR2uEhU0BSNSD1CK3O2A/aq3CXaqOBFDsPM2lG3bvw5wK7QipQiUUMFL3IIZpYKTKLuup1Lgf+hbj5eJGpoDl5EJKC0By8iElAqeBGRgFLBi4gElApeRCSgVPAiIgGlghcRCSgVvIhIQKngRUQC6v8DLwso2ZddntkAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "fun = lambda x: sin(1.1*x) + exp(x/2) - x\n", "xa = -6\n", "xb = 3\n", "\n", "x = linspace(-6, 3, 250)\n", "plt.plot(x, fun(x))\n", "plt.ylabel('Example Function')\n", "plt.xlabel('$x$')\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "wrapped-conjunction", "metadata": {}, "source": [ "Next, I setup the grid and get a function that will compute the coefficients for the interpolating chebyshev expansion. To start, I'll use 6 interpolation nodes." ] }, { "cell_type": "code", "execution_count": 3, "id": "killing-phase", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWoAAAEGCAYAAABM7t/CAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAA7C0lEQVR4nO3dd3hUVf7H8fdJZtJIQhJCCCmQ0CGQhNB7E1CpFpQiCkoRRV0LCitKVv2tu7IrK3YRVhCwIAooriDSpDcBaVIDCTUktBBS5/z+uBBA0oCpyff1PHmWzNw7883Z8ZOTc885V2mtEUII4bzcHF2AEEKI4klQCyGEk5OgFkIIJydBLYQQTk6CWgghnJxNgvqzzz7TwC19paam3vK5Ze1L2kLaQ9qjXLVFkWwS1ElJSbd8bnZ2tvUKcXHSFteT9rietMdVZb0tZOhDCCGcXKmCWikVoJT6Rim1Rym1WynVytaFCSGEMJhKedw7wE9a6/uVUh6Ajw1rEkIIcY0Sg1opVRFoDwwB0FrnADm2LUuIsiM3N5eUlBSysrKs+rp5eXlcuHDBqq/pqlypLby8vIiIiMBsNpf6nNL0qKOBVOC/Sqk4YDPwjNb64q2VKUT5kpKSgp+fH1FRUSilrPa6OTk5eHh4WO31XJmrtIXWmrS0NFJSUoiOji71eaUJahOQADyltV6vlHoHGAu8cu1BSqkRwAiArl27kpKSUuoirpWenn5L55VF0hbXc9X2yMjIICwsjNzcXKu+bl5enlVfz5W5Ulv4+flx/PjxGzIyIiKiyHNKE9QpQIrWev3l77/BCOrraK0/AT4BSExM1MW9aUlu59yyRtrieq7YHhcuXMDT09Mmr+0KvUh7caW2MJlMN/VZLjGotdYnlFLJSqm6Wus/gC7Artuo0Xryc+FcMpw7ChdPwZUtW70qgn84BESCp59jaxRCiNtU2lkfTwGzLs/4OAgMtV1JJTi1G3bOg8Or4ehmyM0EFPgEgZvJCOuss5CfYzweXAcimkHNTlC7G3j5O6x0IcqzpKQkevbsyY4dO0p9zuTJk/nwww9JSEhg1qxZVqtjzZo1DBw4EIBNmzYxY8YMJk+ebJXXt4VSBbXWeivQ1LalFCM3C7bOhE2fwcnfIbQR1OgILUdBlRjwCwPTNX/2WCyQmQZp+40wT14H3/8F8rOhVldo+ijU7Axust5HCGf2wQcfsGTJEqsOeSUlJTF79uyCoG7atClNmzou3krDuZMqPxfWfQST4+GX16D2HfDEenh8FXR7A+r1gMCo60MajAD2rQzVW0Hr0fDgTHjxADw4C0ye8MWD8F4T2DzdeA8hyriZM2fSvHlz4uPjGTlyJPn5+WzcuJHY2FiysrK4ePEiMTEx7Nixg4yMDLp06UJCQgKNGjVi/vz5gBFw9erVY8iQIdSpU4dBgwaxZMkS2rRpQ+3atdmwYQMAiYmJDB48mFatWlG7dm2mTJlyQz35+fmMGTOGZs2aERsby8cff3zDMY8//jgHDx7krrvuYtKkSSQmJvKvf/2r4PmGDRuSlJREUlISsbGxDB8+nJiYGLp168alS5cA2L9/P3fccQdxcXEkJCRw4MABxo4dy6+//kp8fDyTJk1i+fLl9OzZEzAuWPft25fY2FhatmzJ9u3bC36mRx99lI4dO1KjRg27975LO/Rhf0fWww/PwvkUaPssNH2sxGELrTWnM3JIOZNJyplLl78yScvIISsvn6xcP7LzhhEUdB93Zy3krh/GkvnTm/waPoIT1XtTrVIF6ob6Ur1SBczuzv07TLimvHwLx89ZZz51bm4uZrMx26FqRS9MRXxmd+/ezVdffcXq1asxm8088cQTzJo1i4cffpjevXszfvx4Ll26xEMPPUTDhg3Jy8vju+++w9/fn9OnT9OyZUt69+4NGME3Z84cpk2bRrNmzZg9ezarVq1iwYIF/P3vf2fevHkAbN++nXXr1nHx4kUaN25Mjx49rqtp6tSpVKxYkY0bN5KdnU2bNm3o1q3bdVPWPvroI3766SeWLVtGcHAwiYmJRbbF/v37+fLLL5kyZQoPPPAAc+fO5aGHHmLQoEGMHTuWe+65h6ysLCwWC//4xz/417/+xQ8//ADA8uXLC15nwoQJNG7cmHnz5rF06VIefvhhtm7dCsCePXtYtmwZFy5coG7duowaNeqm5kLfDucLaksuLHoZ1r4PsQ/Cw/ON3jFgsWhOZ2STfOYSR89euiGQj565RHaeBYBAHzMRgT6EB3hT2c8Tbw93vExueJrd0TqEw3m1+Tj7EeKOzKTH4TfZnzKHV/OGsCkrArO7okFYRVrVqESrmpVoHhWEt4e7I1tFlBHHz2XR7q1lVn/dX1/sRGRQ4QuGf/nlFzZv3kyzZs0AuHTpEiEhIQC8+uqrNGvWDC8vr4Jeotaav/71r6xcuRI3NzeOHj3KyZMnAYiOjqZRo0YAxMTE0KVLF5RSNGrU6LrN2Pr06YO3tzfe3t506tSJDRs2EB8fX/D84sWL2b59O9988w0A586dY9++fTc1t/haUVFRBa/fpEkTkpKSuHDhAkePHuWee+4BjIUmJVm1ahVz584FoHPnzqSlpXH+/HkAevTogaenJ56enoSEhHDy5Em7zUJyqqA+e/wgft89TE5mCmuafcBWz6acXHyclDMHOXrmEilnL5FzOYgrVfAgPNCbiEBv6oX6cUf9ECICvQkP8CE80Btfz9L+aC3hzAs0WPQyc/4Yy8X2T7E5agSbj2ay9sBppq46iMnNjS71Q+gZG0aX+iHS2xa3rGpFL359sZNVXsvoUZsLXrcoWmseeeQR3nzzzRueS0tLIyMjg9zcXLKysqhQoQKzZs0iNTWVzZs3YzabiYqKKlhVee00Qzc3t4Lv3dzcrpvL/OeFPX/+XmvNu+++S/fu3Uv985pMJiwWS8H31670vLYud3f3gqEPa/rze9hz7rbzBLXFwrkpvTiZV4Hn8/8PvSOUUP/TVKnoRYMwf7o1qEJEoI8RxoHe+HhYsfTAKOg/C7X7e3y/f4YOR5bR4d5PoGtrMnPyWP5HKgu3H+eZL38jwMfMgObVGNSiOpX9bDM3VpRdJne3Inu+N6u0q/G6dOlCnz59ePbZZwkJCSE9PZ0LFy5QvXp1Ro4cyeuvv86hQ4d46aWXeO+99zh37hwhISGYzWaWLVvG4cOHb7q2+fPnM27cOC5evMjy5cv5xz/+QU7O1Z0nunfvzocffkjnzp0xm83s3buX8PBwKlSoUORrRkVFFQxXbNmyhUOHDhVbg5+fHxEREcybN4++ffuSnZ1Nfn4+fn5+RS43b9euHbNmzeKVV15h+fLlBAcH4+/v+JlizhPUbm6YBs7GO78CK+vUtOpS21Kr3wsiW8CCp+HjDnD3RHyaPMLdjapyd6OqnMvMZc7mZD5fd5iPVhxgcMvqjGhfUwJbOLUGDRrwxhtv0K1bNywWC2azmffff58VK1ZgNpsZOHAg+fn5tG7dmqVLlzJo0CB69epFo0aNaNq0KfXq1bvp94yNjaVTp06cPn2aV155hbCwsOuGRoYNG0ZSUhIJCQloralcuXLB+HZR7rvvPmbMmEFMTAwtWrSgTp06Jdbx+eefM3LkSF599VXMZjNz5swhNjYWd3d34uLiGDJkCI0bNy44/spFw9jYWHx8fJg+ffpN/+y2oLQu9sYCtyQxMVEXN/BfnJSUFMevPtMaNkyBReOg8UNw11vGbJHL8i2aH7Yf451f9nHiXBajO9fisbbReJqsO47tFG3hRFy1PXbv3k39+vWt/rrOur9FYmIivr6+vPDCC3Z7T2dti6IU8Zkosncqg62FUQpajIAhC+GP/8F/74ILJwqedndT9IkP5+dnOzC+RwM+/fUQ3SetZPX+0w4sWghRVklQF6daSxi5EpQbfNoVUv+47ml3N8XAFtVY9nxH2tYO5qGp63l1/g4yc1xngxghrC0xMdGuvenyQIK6JH6h8PACqBoLU7vC4TU3HFLRx8wbfRsx87EW/LL7FD0mr2L38fMOKFYIURZJUJeGhw88MAMa9YMZfWDPwkIPa1MrmP/9pR31Qv2454PVfLP51rZ6FUKIa0lQl5abO9z9L2g/Br5+GHZ8W+hh/l5mPhiUwAvd6jJ27nbe+GEXFov1L9gKIcoP55me5wqUgg4vgskL5g4z9gmJe7CQwxTD2tWgflV/Hv98M8fOXeLtB+LxMsvqRiHEzZMe9a1o8zTc+Q+Y9zhs/aLow2oFM2dUK347cpZBn67nbKbcalI4jyFDhhQs4S6NpKQkGjZseMvvN2zYMHbtKn4r+3nz5pV4jDWU5me5ssveFZs2beLpp5+2dWmFkqC+VS1GGPOr5z8Ju+YXeVi9UH++e6INF7PzGDhlPWcuSliL8unTTz+lQYMGxR5zK0Ftq6Xcfw7qpk2bOmzPagnq29F8OHQeD988BvuXFHlYaEUvZg9viQYGfrqedAlrYWczZswgNjaWuLg4Bg8eXPD4ypUrad26NTVq1Liudz1x4sSCLUgnTJhQ8HheXh6DBg2ifv363H///WRmZrJ06VL69u1bcMzPP/9csBHStTp27MimTZsA8PX15eWXXyYuLo6WLVty8uRJ1qxZw4IFCxgzZgzx8fEcOHCAAwcOcOedd9KkSRPatWvHnj17AOOvgccff5wWLVrw4osv8vrrrxe6tarWmjFjxtCwYUMaNWrEV199dUNdSUlJtGvXjoSEBBISElizxpjZ5VTboWqtrf41YcIEfauSk5Nv+VyH+TlR69eraH1oVbGHpWdk67v+s1J3n7RCp2Vkl/iyLtkWNuSq7bFr166r3+Tlap2eZJWv7JP7rn6fl1vk++/YsUPXrl1bp6amaq21TktL01pr/cgjj+j7779f5+fn6507d+qaNWtqrbVetGiRHj58uLZYLDo/P1/36NFDr1ixQh86dEgDetUq43M+dOhQPXHiRG2xWHTdunX1qVOntNZaDxgwQC9YsOCGOjp06KA3btyotdYaKDhmzJgx+vXXXy+oac6cOQXndO7cWe/du1drrfW6det0p06dCo7r0aOHzsvL01prPX78eB0bG6szMzN1amqqjoiI0EePHtXffPONvuOOO3ReXp4+ceKEjoyM1MeOHdOHDh3SMTExWmutL168qC9duqS11nrv3r26SZMmWmutly1bpnv06FFQy7Xfjx49WicmJmqttf7ll190XFyc1lrrCRMm6FatWumsrCydmpqqg4KCdE5OTvGfiauKzFS5mGgNXV6F7Asw+0F49CcILXzsK7CCB7OGtWDgp+t59LONzB7ewrqbSwnnd/4ovBNrlZe6bsH0M9shsHqhxy1dupR+/foRHBwMQFBQUMFzffv2xc3NjQYNGhRsZbp48WIWL15csAdGRkYG+/bto1q1akRGRtKmTRsAHnroISZPnswLL7zA4MGDmTlzJkOHDmXt2rXMmDGj+No9PAp6p02aNOHnn3++4ZiMjAzWrFlDv379Ch7Lzs4u+He/fv1wd796gb6wrVVXrVrFgAEDcHd3p0qVKnTo0KHghglX5ObmMnr0aLZu3Yq7uzt79+4ttnaw/3aokhLWoJQxXn0xFWb1g2FLoGJ4oYcGVvBg+qPNuO/DNTw5awufPNxUtk0tT/zDjVC1gpzcXDyubFzvX/jnrSTXbt2pL+/7o7Vm3LhxjBw58rpjk5KSity+dOjQofTq1QsvLy/69euHyVR8tJjN5oJzi9oy1GKxEBAQULBx/5/9eae9krZWLcqkSZOoUqUK27Ztw2KxlGrf6uLYYjtUSQhrcXODez6GgGow+wHIKnplYoifF9OHNmdbyjle/u73gv9ARDngbjJ6vtb4Cqh29d/uRQdj586dmTNnDmlpaYAxvlqc7t27M23aNDIyMgA4evQop06dAuDIkSOsXbsWgNmzZ9O2bVsAwsLCCAsL44033mDo0Fu/9/W1W5D6+/sTHR3NnDlzAOMXyLZt24o8d/78+WRlZZGWlsby5ctp1qwZ7dq146uvviI/P5/U1FRWrlxJ8+bNrzvv3LlzVK1aFTc3Nz7//HPy8/NvqOXPrmyHCthlO1QJamsye8GALyAvy1gUU8z9GGtU9mXakGZ8v+047y7db8ciRXkTExPDyy+/TIcOHYiLi+O5554r9vhu3boxcOBAWrVqRaNGjbj//vsLAqtu3bq8//771K9fnzNnzjBq1KiC8wYNGkRkZORt7RTYv39/Jk6cSOPGjTlw4ACzZs1i6tSpxMXFERMTU3D/xsJc2Vq1ZcuWBVur3nPPPQUXUTt37sxbb71FaGjodec98cQTTJ8+nbi4OPbs2VPQU792O9RJkyZdd05iYiKbN28mNjaWsWPH2n471OIGsG/1q9xdTPyztANa/7OG1vNHa22xFHvoTzuO6+ixP+hFO47f8FyZaAsrctX2KOLC0W3Lzi75grQ9Pfnkk/rTTz91yHuPHz9eT5w40SHvfStu9mKi9KhtIaiG0bPe9qWxr3UxuseE8nSX2jz71Vb2niz8zywhnF2TJk3Yvn07Dz30kKNLKZPkYqKtRDaHnpOMu8WE1IPo9kUe+nTn2uw+fp7hMzYx/8k2BPi4zgboQgBs3rzZoe//yiuvuNSNA26W9KhtqfFD0HwEfP0InEkq8jA3N8W/H4jH0+TGc19vk4uLZZD8fyquuJXPggS1rXV7w5hX/eUgyM4o8jBfTxMfDEpg7YE0pq4q/qadwrV4eXmRlpYmYS3QWpOWlnbTUwBl6MPW3E3Qbzp80hHmjTL2tS5ifmetED9e79uQsXO30zQqiGAH3N9XWF9ERAQpKSmkpqZa9XXz8vJKnK9cXrhSW3h5ed30AhjX+MlcnU+QcXFxShdY+z60Hl3kofc3iWDN/tOMnr2FKffXtGORwlbMZjPR0dFWf11XvdmvLZT1tpChD3upEmNcXFwyAY6sL/bQ1/s2xMPkxr+Xyx1ihBClDGqlVJJS6nel1Fal1CZbF1VmxQ+A+IEwZwhcLPqO5RU8TfznwXiW7T/Lwu3H7VefEMIp3UyPupPWOl5r3dRm1ZQHd70FFSrBt8PBkl/kYbERATzctArj5/3OqQtZdixQCOFsZOjD3szexsXFlE2wcmKxhz7StArhgd6Mmyv7gQhRnpU2qDWwWCm1WSk1wpYFlQuVakKf92HFP+HgiiIPM7kr3n4gnl/3n2aO3NFciHKrtLM+2mqtjyqlQoCflVJ7tNYrrz3gcoCPAOjatSspKbcWLCXt7FVm+CcQUL8/3t8M4+S932LxCrzhkPT0dIKC4NFmVXjt+53U888jyMfsgGKdQ7n5bJSStMdVZaEtipu1Uqqg1lofvfy/p5RS3wHNgZV/OuYT4BOAxMREfTtTZcryNJvr3PM2TOlM2Ma/Q//Zhc6vjoiI4PmqYSw/tIopm8/y7oDGDijUeZSbz0YpSXtcVZbbosShD6VUBaWU35V/A92AHbYurFwwe8P90+DAUtg0tejD3N34x32xLNx+jGV/nLJjgUIIZ1CaMeoqwCql1DZgA7BQa/2TbcsqR0LqQ/f/g0Uvw8mi774cHxnAw62iGP/dDjJzbHPXZSGEcyoxqLXWB7XWcZe/YrTW/2ePwsqVpo9BrTtg7mOQe6nIw17oXheL1vxnyT47FieEcDSZnucMlILe78Kls/Dzq0Ue5utp4rU+Dfn014PsPHbOfvUJIRxKgtpZ+ATBvR/Dxk9h3413ZL6ia4MqdK5Xhb8t2CVzq4UoJySonUl0e2j5BMwfDZlFTzd6pWd9tiaf5XtZXi5EuSBB7Ww6vwLeAbDw+SIPqV6pAsPbR/P3hbvlwqIQ5YAEtbMxe8E9H8PuBXgf+LHIw57oWAuAD5YdsFdlQggHkaB2RmHx0OElAte8DucLH96o4Gli3N31+GTlQQ6nXbRvfUIIu5KgdlZtnyPPvzosGA1FXDTsHRdGXGRF3li4287FCSHsSYLaWbmbSO/wJiStgs3/LfQQpRQTesWwZPdJ1h1Ms3OBQgh7kaB2YnkB0dD1NWPVYlrhY9ENwyvSNz6cv/+4G4tFpusJURZJUDu7ZsMhohnMfxIslkIPeb5bHfacuMD324/ZuTghhD1IUDs7Nzdj1eLx7cZimEJEBPowtE0UExf9QXZe0XeNEUK4JglqVxBYHbr+DZYkwpmkQg95omMtMrLz+HztYbuWJoSwPQlqV9H0MWPa3oKnCp0FUtHbzFOda/Pu0v2cy8y1f31CCJuRoHYVV4ZAkjfC5s8KPWRwy+pU9DbzwfL99q1NCGFTEtSupFJN6PIKLH4Fzt14qzMPkxvPda3D9LVJnDovdy4XoqyQoHY1LR43bjbw/TOFDoH0igujWpAP7y+TXrUQZYUEtatxczfuYH7oV9g6+4an3d0Uz3Wty+wNR0g5k+mAAoUQ1iZB7Yoq14GOY2HRuEL3AukeU4V6of68+4v0qoUoCySoXVXrpyEwGhY+d8MQiFKK57vV4ZstKRw6LRs2CeHqJKhdlbsJ+n5g3A1mx9wbnu5QpzIJ1QL4z5K9DihOCGFNEtSurEoMtH0Wfhp7wx1hjF51XRZsO8YfJy44qEAhhDVIULu6ds+DV4AxZe9PWtaoRNtawbz98x/2r0sIYTUS1K7O7AW9J8PWWXBwxQ1PP9u1Dot2nmT38fMOKE4IYQ0S1GVB9dbQ5BFjbnXupeueSqgWSLvawby3VGaACOGqJKjLijv+ZoT0in/e8NTTXWrz447j7D0pY9VCuCIJ6rLCOwDungirJ8OJ3697qllUEC2ig6RXLYSLkqAuSxr0hrp3GTvsWa7fl/rpLrX5YfsxDqRmOKg4IcStkqAua+6eaNy2a/1H1z3cqkYlEqoFyh4gQrggCeqyxj8M7pgAS9+AM1dvIqCU4ukutZm/9RiH02S1ohCupNRBrZRyV0r9ppT6wZYFCSto8iiExt6wvLxd7WAahlfkg2WF3yhXCOGcbqZH/Qyw21aFCCtyc4Ne7xjzqn//puBhpRTPdKnF3C0pJKfLznpCuIpSBbVSKgLoARR+d1XhfELqGasW/7S8vFPdEOpV9ePDFdKrFsJVlLZH/R/gRcBiu1KE1bV7DnwqwaKXCx5SSjG6Uy2+2ZQid4ERwkWYSjpAKdUTOKW13qyU6ljMcSOAEQBdu3YlJeXGW0WVRnp6eskHlRPWaAuPluOp/MPDnA7rQnZ4SwAaVNSE+pl556ftjGoddtvvYS/y2bietMdVZaEtIiIiinyuxKAG2gC9lVJ3A16Av1Jqptb6oWsP0lp/AnwCkJiYqIt709spuLy57baIiIATQ6m87g14Yi2YvQEY3QVe/2EXL/VujL+X2QqV2od8Nq4n7XFVWW6LEoc+tNbjtNYRWusooD+w9M8hLZxclwmXl5e/VfBQn8ZhVPA0MXPd4WJOFEI4A5lHXR54B8Ddb8GayXBiBwCeJneGtYtm2qoksnLziz9fCOFQNxXUWuvlWuuetipG2FD93lC7m7HD3uXl5QOaVyM338I3m2/teoIQwj6kR11eKGUsL0/dAxunAlDB08QjraP4ZOVB8vJlQo8QzkqCujypGAFdXoVf/gbnjF70kNZRnLqQxY87Tji4OCFEUSSoy5tmw6ByXfhxDGhNUAUP+jerxofLD6D/dDdzIYRzkKAub9zcoddk2LcYdn8PwLB20ew7eYEVe1MdXJwQojAS1OVRaENo/ZTRq846R0SgD73jw/hwuSwrF8IZSVCXVx1eMha/LPkbAI93qMn6Q+lsOXLGwYUJIf5Mgrq8MntDr//ApmlwZB11qvjRpV4IU1YedHRlQog/kaAuz2p0hLgBxtzqvGyGt6/BTztPyI0FhHAyEtTlXbc34GIqrH6HFtFBNAqvyNRVhxxdlRDiGhLU5V2FStD9TVg5EZW2n+HtavD1pmTOXMxxdGVCiMskqAXEPgDV28D3z3BXTAiVKnjKZk1COBEJamEsL+85CY5uwbR9No+1jWb6WtmsSQhnIUEtDEHR0HEsLB7Pgw08ycmz8N1vRx1dlRACCWpxrVZPQsVqVFg6nkEtqzPl14NYLLKsXAhHk6AWV7mbofc7sPM7RobuJzk9k6V7Tjm6KiHKPQlqcb3wJtB8JAHLxvFAo0A++VUWwAjhaBLU4kadXwZt4XnPb9lwKJ2tyWcdXZEQ5ZoEtbiRpx/0fJug7Z/ySNQZpkivWgiHkqAWhavTHRr0YUz2Byz+PYUjaZmOrkiIckuCWhTtzn9SITOZlwKWMW21LCsXwlEkqEXR/Kqgur3OkJwvWLVxM2czZVm5EI4gQS2K1/hh3MLjed00jVmyrFwIh5CgFsVzc8Ot1zs0ZwfHVs2UZeVCOIAEtShZ5brktXmW5yzT+GnDLkdXI0S5I0EtSsWz4wtor0C8lk+QZeVC2JkEtSgdkyeq12S65Sxl26/zHV2NEOWKBLUotUoxHVkf2JPQleMg95KjyxGi3JCgFjclqO+buOdlcnLhG44uRYhyQ4Ja3JS6UZF8GTya4K0fwkm5sCiEPZQY1EopL6XUBqXUNqXUTqXU3+xRmHBeCd2HsMwSR/Z3o8FicXQ5QpR5pelRZwOdtdZxQDxwp1KqpU2rEk6tTe1gpgc+Bad2w8Ypji5HiDKvxKDWhozL35ovf8n8rHJMKcW9HZvzz7z+6CWJkC77gAhhS6Uao1ZKuSultgKngJ+11uttWpVwej1jw/if590cq1AfFjwlQyBC2JCpNAdprfOBeKVUAPCdUqqh1nrHtccopUYAIwC6du1KSkrKLRWUnp5+S+eVRc7eFvc0qsSo34YwL+NFzi15m4sN+tv0/Zy9PexN2uOqstAWERERRT5XqqC+Qmt9Vim1DLgT2PGn5z4BPgFITEzUxb1pSW7n3LLGmdtiVHAVPt+cytYGz5Kw8W0Cm/WDwOo2fU9nbg9HkPa4qiy3RWlmfVS+3JNGKeUNdAX22Lgu4QL8vMwMaFGNccnN0WFxsGA0aLl8IYS1lWaMuiqwTCm1HdiIMUb9g23LEq5iSOsoDpy+xPpGr0HyRtj8X0eXJESZU+LQh9Z6O9DYDrUIFxQW4E2vuDAm/5ZFyzsmwOJXoNYdEFDN0aUJUWbIykRx24a3q8GaA2nsiOgPoY2MWSAyBCKE1UhQi9vWIMyfdrWDmbIqCfq8D0fWwZbpji5LiDJDglpYxfB2Nfhh+3GOuodBl1dh0Xg4m+zosoQoEySohVW0qx1M7RBfpq06BC0ehyoNZAhECCuRoBZWoZRieLsafLnhCOeyLdD3Q0heDxs/dXRpQrg8CWphNb3iwvD1MvHFhiNQqSZ0e92YBXJ6n6NLE8KlSVALq/EwuTG0TTT/XX2InDwLNH0MotrAtyMgP9fR5QnhsiSohVUNbFGNi9n5fL/tGCgFvd+DM4fg1387ujQhXJYEtbAqfy8z/ZtFMuXXg2itwb8q9JwEK96Co5sdXZ4QLkmCWljd0LbR7DuVwcp9p40HYu6BhvcZQyA5mY4tTggXJEEtrC48wJuesVX5eMWBqw/ePdG4c/mSCY4rTAgXJUEtbOLxDjVZcyCNrclnjQe8A6DvB7BhChxY6sjShHA5EtTCJupX9adzvRA+WLb/6oM1OhqLYb4bBRfTHFabEK5GglrYzJOdarJ410n2nbxw9cE7EsGnEsx/UlYtClFKEtTCZppUD6J5dBAfLr9mrNrsBfdPg4PLjWEQIUSJJKiFTT3RsSbztx0jOf2a2R4h9eDON2HxeDixo+iThRCABLWwsQ51KlMv1I9PVh68/okmQ6BON/jmUZmyJ0QJJKiFTSmleKJjLb7elEzqhexrn4BekyHnIiwa57gChXABEtTC5u5sGEp4gDfTVh+6/gmfILhvCmyZAbvmO6Y4IVyABLWwOXc3xeMdavL52sOcu/SnzZmqt4b2Lxp7V5894pgChXByEtTCLvo2DsfPy8TMdYdvfLL9GKjSEOYMgbwcu9cmhLOToBZ24WFyY3i7GkxddYhLOfnXP+luMqbsnU02ZoIIIa4jQS3spn/zSBQwa30hvWq/ULh/KmycAju+tXttQjgzCWphNz4eJka0r8FHKw7e2KsGiG4PnV42xqvlrjBCFJCgFnY1uFV1LFoze0MRFw7bPmdcYPxqsDF1TwghQS3s62qv+gBZuYX0qt3c4J6PjZBe+LzsByIEEtTCAQa3rE6+RTNrfRG9ap8g6PcZ/P4NbJlu19qEcEYS1MLuKniW0KsGiGhi7Afy4xhI3mjfAoVwMhLUwiGu9KpnF9WrBmg2DGIfgK8egvPH7VecEE6mxKBWSkUqpZYppXYppXYqpZ6xR2GibLvSq/6wuF61UtDjbQiINMI6L7vw44Qo40rTo84DntdaNwBaAk8qpRrYtixRHpSqV23yhAdnwvmjBK5+TS4uinKpxKDWWh/XWm+5/O8LwG4g3NaFibLv2l51ofOqr/ALhQdn4XPwR1j/sf0KFMJJ3NQYtVIqCmgMrLdJNaLcGdyyOlprpq9NKv7AiCacaTMBFv3VuDuMEOWIqbQHKqV8gbnAX7TW5wt5fgQwAqBr166kpKTcUkHp6em3dF5ZVF7aYnBCZd5buo/24e74exX9kUwPbo85ZhAVvhzMqV4zyQusaccqnU95+XyURlloi4iIiCKfK1VQK6XMGCE9S2td6EYMWutPgE8AEhMTdXFvWpLbObesKQ9tMSo0jDm/p/PD/ixevLNescf63fsf+Oo0ob88BcN+Ad/K9inSSZWHz0dpleW2KM2sDwVMBXZrrd+2fUmivPEwufFc1zpMW32IU+ezij/YzR3u+xS8A+GL/pB7yT5FCuFApRmjbgMMBjorpbZe/rrbxnWJcqZPfDjVgyoweWkpNmPyqAADv4KMk/DtCLBYbF+gEA5Umlkfq7TWSmsdq7WOv/z1oz2KE+WHu5tiTPe6fLkhmaTTpdiMyS8UBn5tXFhc8qrN6xPCkWRlonAaXeqHEBcZwNs/7y3dCVUawAPTYd2HsGGKbYsTwoEkqIXTUErx0p31WLDtGDuOnivdSTU7G3cz/9+LsGOubQsUwkEkqIVTaR4dRJd6Ifzfwt3o0q5CbDwIur4G346E/UtsW6AQDiBBLZzOX3vUZ2NSOj/vOln6k1o/Ba1HGzcckN32RBkjQS2cTs3KvjzUsjp//3E3OXk3MaOjywRodD/Muh9O7bZdgULYmQS1cEp/uaM2ZzJzmVHS0vJrKQU9/2Pce3FGHzi931blCWFXEtTCKQX4ePBMl9pM/mUfZy7mlP7EKwtiwhrD9J6QdsB2RQphJxLUwmkNblWdYD9P/rOklNP1rjB5wgMzILQRfCZhLVxfqTdlEsLezO5ujO9Rn+EzNtO/eTX8buZkkyc88Dl8NQim94IhCyEo2lalimLkWzSpF7I5du4SJ85lcSErl8ycfC7l5uOuFB4mNzxN7gT6mAn286SyrydhAd54mKQfeYUEtXBqnetVoVPdyoyft4NJPSNv7mSzFzw4C74cYPSsH1kAlcr3jnu2prVm78kM1h9KY8fRc+w6fp69JzLIyTcuClf0NuPvbcLHbMLLwx2LRZOTZyE7L58zmbmcu5QLGCtVoyr5UCvElzpV/GhcLYCEaoEE+Hg48sdzGAlq4fQm9Iqh66QV/G+3LyMjbyGs+8+Grx+GaXfC4G+NIRFhNZk5eSzbk8qinSdYc+A0pzNyiAzyJjYigLsaVuX5rv5EBvkQFuCFj0fxkZOdl8+p89kcTstk/6kL7DuVwdoDaXy88iA5eRZqhfjSPDqIjnUq06ZWMBU8y0eElY+fUri0yCAfnupcmw9W7OeBNvUIrHCTvSqzt9Gznvc4fNYDBs6Bai1sU2w5kW/RrNybypzNySzdcwo3pehUL4QXu9ejda1KRAT63NLreprciQzyITLIh7a1gwsez87LZ+ex82xOOsOaA6d56ovf0BqaRQfSqW4I8cFQdjc5laAWLmJ4uxp8vSGJtxbt4c17Y2/+BUwecO8UWPg8fN7XuA9jrS5Wr7OsO3Uhi683JvPFhmROXciiW4NQ3unfmA51KuNldrfZ+3qa3EmoFkhCtUCGt69BVm4+6w6msfyPVP67OomjZy+RUO04PWLD6NGoKqEVva5/AYsFss/DpTOQdRb0lfn5Crwqgk8QeFYEN+ccF5egFi7Bw+TGcx3C+cv8g/RrGklCtcCbfxE3d+g5yfgPc/aDcM9HxgIZUaIjaZl8tPIA32xKoUpFTx5qWZ1+TSMI9vV0SD1eZnc61g2hY90QJvRqwM9b9rHheB7TV/7B9z8uoGfwSdoGpFHD7RQe55Pg7BGw5BX/om4mCIyCyvUgpD5ENINqLY3Pi4NJUAuX0STCjz5xYYyb+zvfP9X21mYFKAVd/wYVguHb4ZB+ENqPMR4XNziYmsG7S/ezYNsxYsL8eXdgY7rWr4Kbm5O0V/YF1OE1tDr0A93StvNy3k6URy7pWWFsTw5nVV4IpuCe1GrekIQGdfDxDwYvfyOUwehZZ52HzDS4eArSDkLqbkhaBasngyUXwptAgz7GV0A1h/yYEtTCpbzSswHdJq3kvWX7ea5rnVt/odZPGb2nuZfDutc7xpQ+AcDpjGz+s2QvX2xIpkm1QD4b2oy2tYJRzvALLfUP2L0A9v8CKRtBuWOu0hjq9UR1eRXCGhPkE0Q7i8Z8MI35W4/y73UnyFl1gjvqa3rHm+lYtzKepstDNZ5+UDHc+Heta94nNwuOboK9P8H6T2DxeKjRCZoPhzp3Gn+h2Ykq9Q5lNyExMVEnJibe0rkpKSll+t5nN0Pa4npX2uOH7cf4y5dbmT+6DTFht/ln6dHN8MUAqFTLGLf2CbJOsXZgi8/HpZx8Pv31IB+tOEBYgDfj7q5Hp7ohjg1oreH4NiOcd38Pp/dCcF2o0x1qdoJqrUg5mVZsW2Tn5bP8j1QWbD3Gkt0n8TC5cVfDUHrHhdOqZiXcS/oLQWvjl8LGqbDzW6gYAR3GGkNn1gvsIouQHrVwOT0aVeWHbccZM2c780e3wex+GxeAwpsYN8md/SBM6WQskql6CxcrXZzWmkU7T/La9zvJtWjG92xAvyYRmG6nbW9XRipsmw1bZkDafqgaD7EPQv1eULnuTb2Up8md7jGhdI8J5UJWLj/vOsn8rcd45L8bCKrgQc/YqtzdqCqNIwMK/5mVgsjmxle3N2DNZPj+Gfj133DXP41fGDYkQS1cjlKK1/s2pOukFXyw7ADP3FH79l4wIBIeWwzzn4SpXY1hkLj+1inWBSSnZzJhwU5W7k3lsbbRPN2ltuPmJ1sscHCpEc57fgTfKtD4IYgfYAxVWYGfl5l7EyK4NyGCtIxsfvz9OPO3HuOzNUn4eppoWyuY9nUq0zw6iOhKFW4cj/etDN1eh1ajYfmb8Pk9EHMP3PmmcYs4G5CgFi6psp8nr/VpyLNfbaV9nWAa38oskGt5+kK/z2Dt+zDvCePP3O5/L9Pj1rn5Fj5ZeZB3l+6jUXhFFj7djrqhN7VQ33rOH4PfZsKWz+HCMWMMeMAXxh18bDgWXMnXk8GtohjcKoq0jGxW7T/Nir2pvP3zXlIvZOPvZSIuMoC6VfyoXsmH6pUqEOzriZ+XCS9zADnt3oTo+6m49CU8JjdjW9wEmvUcZvU6JaiFy+odF8byP07xzJdb+fGZdvjebi9QKePmA2HxMGcoJK+H+6be9J/ZrmD/qQye+3orh9Myea1PQ+5PiLD/TI78PNi3GLZMN/43oDo0HQrxg8Cvin1rwQjtPvHh9IkPR2tNyplL/JZ8lm3JZ9mfmsEve06RnJ5JnuXG63om9RJjfBbS6o/50OMxq88ikqAWLu21Pg25+51feXX+Dt5+IN46LxrVFkatgflPwMftofv/QVPr/8fnCBaL5rM1Sfzzpz20qlmJKQ83pYq/V8knWtOZJKPnvHWWMS2ufi8YPA+i2jnNghOlVMEKyd5xYQWP5+VbyMjO40JWHlm5+XiY3PAwuVGpgicept7GRUcbfE4kqIVL8/U08U7/eO7/aC0d6lSmT3y4lV64Mgz8GjZ+Cotehn0/G4tl/MNKPtdJpZzJZMyc7WxLOUti7xj6N4u032yOvGzYs9DoPR9cbiwqaf20cS3AhWbamNzdCPDxKHpzKBu1pwS1cHmNqwXyXNc6/PXb34kJq0itEF/rvLBSxpzZqLbw3ePwfgvjJroJjzhNz680tNZ8szmFv32/i/pV/fjpmfZUq3Rre3HctNS9Rjhv+wJyMqHhvfDoYmP2RBn4C8VeJKhFmTCqQ01+O3KGkZ9vYv7otrc/Xn2tkPrGFL71H8FP4+D3OcYtvyrfxoIbOzmdkc24b39nxR+pPN+tDsPa1Sh5zvDtyr0EO+cZAX1kLVSNg04vG3OOnWA5titynW6BEMVwc1P8+4F48i2aMXO2YfWFXO4m40LjE2vB3QM+bGWE9qUz1n0fK/ppxwm6T1rJ0TOXWPBUG0Z2qGnbkD6+HRa+AP+qCz+OMX7BjVgBI1dCs8ckpG+D9KhFmVHR28zHg5vS9/3VfLTiIKM62uAmAUHRMPg7+ONHY+x625fQ6a/QZAi4m63/frfgfFYuiQt2Mu+3ozzRsRZPd6ltu7ulZKYbf2H8NhNObIeI5nDn3415xR4VbPOe5ZAEtShT6ob6MbFfLE9/8RvVK/lwd6Oq1n8TpaBeD6h1B6z/GH55Dda8a2zuFNffoYG9ev9pxszZhqfZnTmPt6ZJ9ducX14YSz4cWAZbZxoXCL0DjZ/7vqkuMRzkiiSoRZnTMzaMI+mZ/OWrrVTx96RJdRvNKjB5QpunjZVza9+Hn8bCyonQ7nljqbPZftPeLuXk88+f9jB9bRKDW1Zn7F31Srybyk3RGo5vhR1zYce3kHHSWJTywOfGLyx3iRJbKrF1lVLTgJ7AKa11Q9uXJMTtG9WhJilnLjFs+ibmjmpNjcpWmglSGJ8g6PIKtHrSCOzF4+GXvxnDIc2G2XxK39bkszz31VYyc/KZPrQ57etUtt6Ln9x1OZznwplDEN7UWDrdqJ8xhVHYRWl+DX4GvAfMsG0pQliPUorXesdw/OwlHp62ga9HtiIswNu2b3olsNs+a0xHW/8xrH7H6HnG9Yfa3ay6JD0n38LERXv4aMVBesZW5bXeDanoc5vDLvl5kLzO2Npz7yJjp7rQRtDkEWPc2Ur7bYibU2JQa61XKqWi7FCLEFZlcnfj/UEJDPnvRgZOWceXI1rdeIsmW/D0NeZfN30MDi4zVuDNHQYmL4jpC/V6GqvwbmNoZOexczz99T7OZll4f2Bj7mx4G2PxZw7D4dXG/s77f4bsC8ZFwbgBRq0y7uxwpdqP+nJQ/1Dc0IdSagQwAqBr165Npk2bdksFpaenExTkOiuVbEna4nq32h6ZOfk8v+AgZ7PyebdvTYJ97X+xT+Vk4J20BJ8DC/E8vgGtTGSHtSArsi05VRqTG1j76l1HipGXr5m55RT/3XiCFuHejOtWg0DvmxgftuRiPrMf8+ldeJ7YjOfxjZgyjpHvFUh2WEuyIjuQFdkWi5cNLkLaUFn4byUiIqLIuZNWC+pryY0DrEPa4nq30x4XsnJ5ZNoGTl3IZuZjLYgKduDUsewLcHAF7FsEB5bDuSPg4QvhCcbikOA6xsb4wbWNGRWXV/DtOnael+Zu50h6Jq/1iSGhkoXIyMjC3yP3krGnRvpB4yttv7H5/smdkJ8DflUhsoWx6jKqrbGk24VXCpaR/1bkxgGifPPzMjNzWAuemLWF+z9aw2dDm9Mw3EELMDz9oH5P4wvgwglI3mDs1ndiB+z4Ds6nGM+ZvLH4hnAsz59j58286udPw/pV8E7+kYw9GeBphvxsYy+NzHS4mAqZp68uxDH7QFANY/533R7Q8a/GLwMH7E4nbp0EtSg3fDxMTHm4KS9+s51+H63l3w/E2Wae9c3yC4UGvY2vK7IvQNoBtu/5g4VrthGk0+le14MofzcjlLPO4paTBZ4BRm/cpxJUaWjctNenkvGaQTWMjfdduKcsDKWZnvcF0BEIVkqlABO01lNtXZgQtmB2d+PtB+KoFeLL6NlbGNWxJs91rWv7/S9u0sHzijeX5PHLbl+GtH6IQd3q3LB/SXpKCj6u/+e+KIXSzPoYYI9ChLAXpRRPdqpF/ap+PPPlVtYfTGfSg/FEBtlpR7linM3M4Z1f9vH52sM0jQpkwei2jhuiEU5DNmUS5VbnelX46S/tMbkr7nrnV2avP0J+IXfvsIfzWblM/mUf7d9axrI9p/hgUAJfDG8pIS0AGaMW5Vx4gDezh7Vk2upDvLFwF19uPEJi7xgSbvcejKWUlpHN5+sOM23VISp4mhjTvS4PNqtmu02UhEuSoBblnpubYli7GvSIrcrff9zDvR+soUu9EEZ3rnX7N80thNaa34+eY/qaw3y//RiVfT158c569GsagafJdjdyFa5LglqIy6pW9ObdAY0Z1jaad5fu554P1tCkeiD9mkTQI7Yqfl63t1DmQGoGP2w7zvfbj7H/VAZtawXz3oDGdKlfxekuZgrnIkEtxJ/ERQbw6SNN2XPiPF9uSOYfP+3h1QU7aR4VRPs6wcRHBlI31I+K3kUHd06ehf2nMth9/Dwbk9JZfeA0yemXqFPFlz5xYfSMCyPakYtuhEuRoBaiCPVC/UnsHcO4u+uxev9pVu49zVcbk3nzf3vQGoJ9PQn29SDQxwM3N8i3aDKy8zh1PpvTGdlYNFT28yQ+MoDH2kTTtnYwtUL8HP1jCRckQS1ECTxN7nSuV4XO9YzVfJk5eew9mUHS6YukX8zhTGYOWoObAh9PEyF+noT6e1En1I9gX+vtlifKLwlqIW6Sj4eJ+MgA4iMDHF2KKCdkDpAQQjg5CWohhHByEtRCCOHkJKiFEMLJSVALIYSTk6AWQggnJ0EthBBOToJaCCGcXKlubnvTL6rUp0DKLZ7eBNhsxXJcmbTF9aQ9riftcVVZaIskrfVnhT1hk6C+HUqpTVrrpo6uwxlIW1xP2uN60h5XlfW2kKEPIYRwchLUQgjh5JwxqD9xdAFORNrietIe15P2uKpMt4XTjVELIYS4njP2qIUQQlxDgloIIZyc0wa1UuoppdQepdROpdRbjq7H0ZRSzyultFIq2NG1OJJSauLlz8V2pdR3SqkAR9dkb0qpO5VSfyil9iulxjq6HkdSSkUqpZYppXZdzopnHF2TLThlUCulOgF9gDitdQzwLweX5FBKqUigG3DE0bU4gZ+BhlrrWGAvMM7B9diVUsodeB+4C2gADFBKNXBsVQ6VBzyvtW4AtASeLIvt4ZRBDYwC/qG1zgbQWp9ycD2ONgl4ESj3V3611ou11nmXv10HRDiyHgdoDuzXWh/UWucAX2J0asolrfVxrfWWy/++AOwGwh1blfU5a1DXAdoppdYrpVYopZo5uiBHUUr1AY5qrbc5uhYn9CjwP0cXYWfhQPI136dQBoPpViilooDGwHoHl2J1Dru5rVJqCRBayFMvY9QVhPGnTDPga6VUDV1G5xKW0BZ/xRj2KDeKaw+t9fzLx7yM8WfvLHvWJpyTUsoXmAv8RWt93tH1WJvDglprfUdRzymlRgHfXg7mDUopCxAMpNqrPnsqqi2UUo2AaGCbUgqMP/O3KKWaa61P2LFEuyruswGglBoC9AS6lNVf3sU4CkRe833E5cfKLaWUGSOkZ2mtv3V0PbbgrEMf84BOAEqpOoAHcNqRBTmC1vp3rXWI1jpKax2F8WduQlkO6ZIope7EGK/vrbXOdHQ9DrARqK2UilZKeQD9gQUOrslhlNGDmQrs1lq/7eh6bMVZg3oaUEMptQPjYskj5bDnJAr3HuAH/KyU2qqU+sjRBdnT5Qupo4FFGBfOvtZa73RsVQ7VBhgMdL78ediqlLrb0UVZmywhF0IIJ+esPWohhBCXSVALIYSTk6AWQggnJ0EthBBOToJaCCGcnAS1KPOUUlGXF8kI4ZJkep4o0y6vcn0a8AUOAP3L84Ih4ZokqEWZpZTywwjnO4FYYDmQBswE5mqtZyilRgLttdaDHFaoECVw2F4fQtiBBWNr2CAArXUSgFJqBLBaKXUIeB5j8y8hnJYEtSiztNYXlVLDgTeBUKVUQ+BVrfVJpdSrwDLgHq11ukMLFaIEcjFRlGla6wVAP+AtoDJGDxqgEcYwSJiDShOi1CSoRZmllPJVSlW//O2Vu3/4KaWaY9zKqjHwglIq2lE1ClEaMvQhyjIz8DFQCWM/8yPAwxjbgg7VWh9TSj0PTFNKdZYdGoWzklkfosy7fIumjlrrzxxcihC3RIY+RHlwFtjq4BqEuGXSoxZCCCcnPWohhHByEtRCCOHkJKiFEMLJSVALIYSTk6AWQggn9//VuLzQ/VORSAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "n = 6\n", "#create grid and get function for computing coefs\n", "xhat, θ, x, fcheb, S = cheby_dct_setup(xa, xb, n)\n", "#sample the example function at the grid points\n", "y = fun(x)\n", "#compute expansion coefficients super easily!\n", "a = fcheb(y)\n", "#see how well the interpolation does\n", "xx = linspace(xa, xb, 500)\n", "yfun = fun(xx)\n", "yche = cheby_sum(xx, a, xa, xb) #<-- evaluates cheby expansion\n", "plt.plot(xx, yfun, label='example function')\n", "plt.plot(xx, yche, label='cheby interpolation')\n", "plt.xlabel('$x')\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "stable-booking", "metadata": {}, "source": [ "Now let's see how much improvement we get if we add one more point." ] }, { "cell_type": "code", "execution_count": 4, "id": "fancy-bobby", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWoAAAEGCAYAAABM7t/CAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAyu0lEQVR4nO3dd3hUZd7G8e8zk0kvEJLQgoQOARJaIKFKR0ABBZcq0lEQUUGxINnVXV3ru1gXBBEFZQEFRKX33qSD1AAJJSGBEErqPO8fwQAqSQiZnJnJ73Ndc8lkzpy552G4PTlzznOU1hohhBD2y2R0ACGEELmTohZCCDsnRS2EEHZOiloIIeycFLUQQtg5mxT1jBkzNFCgW0JCQoGf62w3GQsZDxmPYjUWd5WvolZKlVBKzVNKHVZKHVJKReW2fExMTH5W+5fS0tIK/FxnI2NxJxmPO8l43OLsY+GSz+X+AyzRWvdUSrkCnjbMJIQQ4jZ5FrVSyg9oCTwJoLVOB9JtG0sIIcTv8rNFXQlIAL5USoUDO4FntdbXbl9IKTUcGA7Qvn17YmNjCxQoKSmpQM9zRjIWd5LxuJOMxy3OMBbBwcF3fSw/Re0CNACe0VpvVUr9B5gATLx9Ia31FGAKQHR0tM7tRe8ncHEjY3EnRxyPjIwMYmNjSU1NLdT1KqVISUkp1HU6KkcaC3d3d4KDg7FYLPl+Tn6KOhaI1VpvvXl/HtlFLYTIh9jYWHx8fAgJCUEpVWjrTU9Px9XVtdDW58gcZSy01iQmJhIbG0ulSpXy/bw8j/rQWp8Hziilatz8UVvgYMFiClH8pKamUqpUqUItaeGYlFKUKlXqnn+7yu9RH88As24e8XECGHSP+YQo1qSkxe8K8lnIV1FrrXcDje557fYsKxOuxMG1i3AtAVKTs3+uTODuB96B4FcBvAKMzSmEKPbyu0Xt2LSGhMNwYi3EboP4w5B4FLJuHmXo4gEeJW4ua4UblyHr5gH0PmWhTF2o1AqqdYCAaiBbR0Lcs5iYGLp27cr+/fvz/ZzJkyfz2Wef0aBBA2bNmlVoOTZt2kTfvn0B2LFjBzNnzmTy5MmFsn5bcJqizsyyci45lTOXrhObdIPr8Sfwj99M+UvbqHJ1FyWslzhnKst+l1BiTI0549GLSx4PkOURiKePH8ElPShfwoPKgV6ElvHFw3oVLp+C8/vh7C7Y/gUsexUCakCDJyC8D3iVMvptC+HUPv30U1asWFGoR/vExMQwe/bsnKJu1KgRjRrZ9w4DhyrqlNQMTiVe51TidWISr3E68TpnLmXf0i5foAn7aWY+QAvzQcpzgctmf054N2BDxadJDIzkhld5XEwKs1IEW62USMviWlomF6+mseHoReIu3+Bccipmk6JakDcNK5akdY12NGvfG4/OJkg8Bnv/B5s/hpX/gEaDoPlz4FPG6KERIlfffPMNkydPJj09nSZNmvDpp5+ya9cuhgwZwrZt28jKyqJx48bMmTOHkJAQunXrxqVLl8jIyODNN9+kW7duxMTE0KlTJyIjI9m0aRMREREMGjSISZMmER8fz6xZs2jcuDHR0dEcP36cY8eOcfHiRV588UWGDRt2R56srCwmTJjAmjVrSEtLY9SoUYwYMeKOZUaOHMmJEyd46KGHGDx4MMnJyXh7ezNu3DgA6tSpw+LFiwHo1KkTLVq0YNOmTZQvX56FCxfi4eHBsWPHGDlyJAkJCZjNZubOncuECRM4dOgQ9erVY+DAgdSvX5/33nuPxYsXk5SUxODBgzlx4gSenp5MmTKFsLAwoqOjOX36NCdOnOD06dOMHTuWMWPGFM1fHnZW1Bdjj/DbuUwOXbnAhSupXLiSStylG8QkXuNU4nUSr2XvqvB1dyGi5HVauB+jd9ZBKpn2UsL1GFZXXwhpjqnKC1CpFSUCa9BAKRrcQ4YrqRnsj01mb1wym48n8vTsXQC0qBpAz4bBtG35Mq6tXoJDi2Dtv2HnDIgaDS3HgcWj8AdFOJXff/MrDBkZGVgsmQCU9XPHxfzXB3EdOnSIOXPmsHHjRiwWC08//TSzZs3iiSee4JFHHuG1117jxo0b9O/fnzp16pCZmckPP/yAr68vFy9eJDIykkceeQSAY8eOMXfuXKZPn05ERASzZ89mw4YNLFq0iH/9618sWLAAgL1797JlyxauXbtG/fr16dKlyx2Zpk2bhp+fH9u3byctLY1mzZrRoUOHOw5Z+/zzz1myZAmrV68mICCA6Ojou47FsWPH+O6775g6dSqPP/448+fPp3///vTr148JEybQo0cPUlNTsVqtvP322znFDLBmzZqc9UyaNIn69euzYMECVq1axRNPPMHu3bsBOHz4MKtXryYlJYUaNWrw1FNP3dOx0PfDroo6dVo3Iq2J7KE6nq4PUMa9DJU9PenmZaFMieuUUlfwu34KS+JvcOli9pd+D0RB2AAIaYGpbDiY7+8t+bpbaFo1gKZVAxjZqgrX0zPZdCyRn/ad47n/7cbDYqZfk4o82awLAaHdYP98WDYR9s+DLu9D1XaFNBrCGZ1LTqXFO6sLfb3rX2xNBf+/noJn5cqV7Ny5k4iICABu3LhBUFAQAK+//joRERG4u7vn7KPVWvPKK6+wbt06TCYTcXFxXLhwAYBKlSpRt25dAGrXrk3btm1RSlG3bt07JmPr1q0bHh4eeHh40Lp1a7Zt20a9evVyHl+2bBl79+5l3rx5ACQnJ3P06NF7Orb4diEhITnrb9iwITExMaSkpBAXF0ePHj2A7BNN8rJhwwbmz58PQJs2bUhMTOTKlSsAdOnSBTc3N9zc3AgKCuLChQtFdgKWXRW1aehSLh1ZTlTaMdSlGLiyM/sLv2uAZ6nsW8UmEPEklK4DgTXBZNsptT1dXWgXWpp2oaX5e7faLNx9lmnrTzB1/Ql6R1TgmbbdCBjdEVb9E2b1gsbDof0/wMXNprmEYyrr5876F1sXyrqyt6gtOeu9G601AwcO5K233vrTY4mJiVy9epWMjAxSU1Px8vJi1qxZJCQksHPnTiwWCyEhITnH/bq53fpcm0ymnPsmk4nMzMycx/54CNof72ut+eijj+jYsWO+36+LiwtWqzXn/u3HIt+ey2w2c+PGjXyvN7/++Bq3v19bs6uiLlf+AWJ1R1TwEKOj/CVfdwsDIivSt/EDLNl/ng+W/8b8XXGMbFWZoe3fwr1mZ/h+OJzaBI9/Bf6VjY4s7IyL2XTXLd97ld+z8dq2bUu3bt147rnnCAoKIikpiZSUFCpWrMiIESN44403OHnyJC+99BIff/wxycnJBAUFYbFYWL16NadOnbrnbAsXLuTll1/m2rVrrFmzhrfffpv09FtzuXXs2JHPPvuMNm3aYLFYOHLkCOXLl8fLy+uu6wwJCcnZXbFr1y5OnjyZawYfHx+Cg4NZsGAB3bt3Jy0tjaysLHx8fO56unmLFi2YNWsWEydOZM2aNQQEBODr63vP77+w2VVROwqzSdElrCwda5dmzo4zfLDsCN//GsdbPerSZORG+H4oTG0LvWdDxVyn7hbC5kJDQ3nzzTfp0KEDVqsVi8XCJ598wtq1a7FYLPTt25esrCyaNm3KqlWr6NevHw8//DB169alUaNG1KxZ855fMywsjNatW3Px4kUmTpxIuXLl7tg1MnToUGJiYmjQoAFaawIDA3P2b9/NY489xsyZM6lduzZNmjShevXqeeb4+uuvGTFiBK+//joWi4W5c+cSFhaG2WwmPDycJ598kvr16+csHx0dzeDBgwkLC8PT05Ovvvrqnt+7LSitc72wQIFER0fr3Hb85yY2NtbhJt65fD2df/18iLk7YxkQWZFXOlXDffnL8OvX0O0TCHu8QOt1xLGwJUcdj0OHDlGrVq1CX6+9zm8RHR19x9EZRcFex+Ju7vKZuOsJGrJFXQhKeLryTs9wHgkvzwtzd7PtZBIf94mmaqmq8MOI7LMeGw/Le0VCCPEXpKgLUfNqAfzybEvGz93Dwx9v4v3He9D5ET9YNBoy06DpaKMjCmFzBf1tWtydXIW8kPl7ufLFwEaMblOV0bN38eHFCKw9psLy12HzJ0bHE0I4INmitgGlFKNaV6VakDdj5+zmVGgV3uv+OS4LRoCHP9TrY3REIYQDkS1qG+pQuwz/GxHF+qMXGf5rJTI6vA0LR8Hhn42OJoRwIFLUNlanvB9zR0Zx+NwV+u0NI63ZCzB/CJzbY3Q0IYSDkKIuApUDvZn3VFMuXk2j92+tyKzWGb7tAykXjI4mirEnn3wy5xTu/IiJiaFOnToFfr2hQ4dy8GDuF4dasGBBnssUhvy8l99n2fvdjh07inQipttJUReRciU8+G5YJEnXMxh86Qms3qVhTj/IKNwLngphr7744gtCQ0NzXaYgRW2rU7n/WNSNGjUybM5qKeoiFOTrzqyhTTialMkLppfQyXGweGz2hQ2EsKGZM2cSFhZGeHg4AwYMyPn5unXraNq0KZUrV75j6/rdd98lIiKCsLAwJk2alPPzzMxM+vXrR61atejZsyfXr19n1apVdO/ePWeZ5cuX50yEdLsHH3yQHTt2AODt7c2rr75KeHg4kZGRXLhwgU2bNrFo0SLGjx9PvXr1OH78OMePH6dTp040bNiQFi1acPjwYSD7t4GRI0fSpEkTXnzxRd544w0GDBhAVFQU1apVY+rUqUD2nCLjx4+nTp061K1blzlz5vwpV0xMDC1atKBBgwY0aNCATZs2ATBhwgTWr19PvXr1+PDDD1mzZg1du3YFICkpie7duxMWFkZkZCR79+4Fbp3Z+OCDD1K5cuXCK3atdaHfJk2apAvqzJkzBX6uozgen6IbvrFMvzd9trb+vZTWO7/6y+WKw1jcC0cdj4MHD966k5mhdVJModzSLhy9dT8z466vv3//fl2tWjWdkJCgtdY6MTFRa631wIEDdc+ePXVWVpY+cOCArlKlitZa66VLl+phw4Zpq9Wqs7KydJcuXfTatWv1yZMnNaA3bNigtdZ60KBB+t1339VWq1XXqFFDx8fHa6217tOnj160aNGfcrRq1Upv375da601kLPM+PHj9RtvvJGTae7cuTnPadOmjT5y5IjWWustW7bo1q1b5yzXpUsXnZmZqbXW+rXXXtNhYWH6+vXrOiEhQQcHB+u4uDg9b9483a5dO52ZmanPnz+vK1SooM+ePatPnjypa9eurbXW+tq1a/rGjRtaa62PHDmiGzZsqLXWevXq1bpLly45WW6/P3r0aB0dHa211nrlypU6PDxca631pEmTdFRUlE5NTdUJCQna399fp6en5/6ZuOWunSqH5xmgcqA3MwY1ptfnm2lQdQytf34RyjeC0rn/WiicwJU4+E9YoazqjhOmn90LJSv+5XKrVq2iV69eBARkX//T398/57Hu3btjMpkIDQ3Nmcp02bJlLFu2LGcOjKtXr3L06FEeeOABKlSoQLNmzQDo378/kydPZty4cQwYMIBvvvmGQYMGsXnzZmbOnJl7dlfXnK3Thg0bsnz58j8tc/XqVTZt2kSvXr1yfpaWlpbz5169emE2m3Pu/9XUqhs2bKBPnz6YzWZKly5Nq1at2L59O2Fht/4OMjIyGD16NLt378ZsNnPkyJFcs0PRT4cqRW2QOuX9mNynPkO+zmRDxaaUm/skDF8NrnefPUw4Ad/y2aVaCNIzMnD9feJ63/IFWsftU3fqm7vgtNa8/PLLf7riSkxMzF2nLx00aBAPP/ww7u7u9OrVCxeX3KvFYrHkPPduU4ZarVZKlCiRM3H/H/1xpr28pla9mw8//JDSpUuzZ88erFZrvuatzo0tpkOVfdQGah9amlc6h/LImT6k3bgKP483OpKwNbNL9pZvYdxKPHDrz7lcMKNNmzbMnTuXxMREIHv/am46duzI9OnTuXr1KgBxcXHEx8cDcPr0aTZv3gzA7Nmzad68OQDlypWjXLlyvPnmmwwaNKjAw3P7FKS+vr5UqlSJuXPnAtn/A9mz5+6HtS5cuJDU1FQSExNZs2YNERERtGjRgjlz5pCVlUVCQgLr1q2jcePGdzwvOTmZsmXLYjKZ+Prrr8nKyvpTlj/6fTpUoEimQ5WiNtiQ5pXo2KgmI1NHofd8C4d+NDqScDK1a9fm1VdfpVWrVoSHh/P888/nunyHDh3o27cvUVFR1K1bl549e+YUVo0aNfjkk0+oVasWly5d4qmnnsp5Xr9+/ahQocJ9zRTYu3dv3n33XerXr8/x48eZNWsW06ZNIzw8nNq1a7Nw4cK7Pvf3qVUjIyNzplbt0aNHzpeobdq04Z133qFMmTuvcfr000/z1VdfER4ezuHDh3O21G+fDvXDDz+84znR0dHs3LmTsLAwJkyYYPPpUGWaUzuQnmnl8f9u5onrX9FDr0SN2gpeAcVyLHLjqONRXKY5HT16NPXr12fIkKK/8MfEiRPx8/Mr0qlV78e9TnMqW9R2wNXFxKf9GvDvG92J1yXkkD3hcBo2bMjevXvp37+/0VGcknyZaCfKlfDg/T6NGfTlEBYffh3TvnngL1eHEY5h586dhr7+xIkT7eq3i8ImW9R2pHm1ADq1bc8n+lGsP43DdCPR6EiikNhiF6NwTAX5LEhR25lRrauyqfQAzmSVwHfLv42OIwqBu7s7iYmJUtYCrTWJiYn3fAig7PqwM2aT4t3eDXnp/4byzfGJcHwVVGljdCxxH4KDg4mNjSUhIaFQ15uZmZnn8crFhSONhbu7+z1/Ke4Y76yYCS7pyeM9HuWb+Wt5fMGzuD2zFVw9jY4lCshisVCpUqVCX6+jHgVjC84+Fvna9aGUilFK7VNK7VZK7bB1KAHd6pVnR8XhpFxNIXON7AIRoji7l33UrbXW9bTWjWyWRtxhxIPVec80BLX5I0jIe/4BIYRzki8T7ZivuwttHx3GlqyaXFk4To6tFqKYyu8+ag0sU0pp4L9a6yl/XEApNRwYDtC+fXtiY2MLFCiveQiKk6SkJGr5+zOl3GgiY0dxfuM3ZIa0NjqWYeSzcScZj1ucYSxy28ee36JurrWOU0oFAcuVUoe11utuX+BmeU+B7FPI72fHvjN/KXCvgoODGTugJ3PfW0THDf+mTOTj4OKW9xOdlHw27iTjcYszj0W+dn1oreNu/jce+AFonPszRGHy93LFr/Mk9I1LJK78j9FxhBBFLM+iVkp5KaV8fv8z0AHYb+tg4k6dGtVkof9gPLd8gL5yzug4QogilJ8t6tLABqXUHmAb8JPWeoltY4k/UkrRqvc4YqxBnJn3itFxhBBFKM+i1lqf0FqH37zV1lr/syiCiT+rXNqPA7XHU/70Aq6eKZyrhAgh7J8cnudguvboy05TGGfnvmR0FCFEEZGidjDuFjOq/d+pmryZEzuWGh1HCFEEpKgdUETUg2z1aUvW0oloq9XoOEIIG5OidlAVHv0nFdOPs3upba/VJoQwnhS1gwquXJNdZXoRuO3fpKWlGh1HCGFDUtQOLPRv/8BPJ7N1vpwEI4Qzk6J2YL7+QRyv+iTVf/uci5eTjY4jhLARKWoHV+fRl/A0ZbD5fx8YHUUIYSNS1A7OxbMEF8NG0CRuBsdi442OI4SwASlqJ1C583O4mTW//vCe0VGEEDYgRe0M3Ly51mg0rS9+y57jBZsHXAhhv6SonUS5dqOxuLhwYMF7aLkSjBBORYraWbh6kt70OR66MpdNB08anUYIUYikqJ1IYMvhKFdPTix+D6tVtqqFcBZS1M7E4o5qPpau1xey5NdjRqcRQhQSKWon49d0MBaLK6eWfkxGlkzYJIQzkKJ2NhYPTM3G0CvtB37YJlvVQjgDKWon5Nl0GF4WOLNyCumZslUthKOTonZGbt6oyKfom/k932+XI0CEcHRS1E7KvdlI/M2pnFg5TbaqhXBwUtTOyqMkOmIY/TLnM0+2qoVwaFLUTsy9xTOUM13it5Vfy1a1EA5MitqZeQVgbfAk/TPnMXfHKaPTCCEKSIraybm1HEsldZ59K74lLTPL6DhCiAKQonZ2vuWwhvWmT+b3/G/7GaPTCCEKQIq6GHBt8SxhHGPz6sVytqIQDkiKujgIqEZGtc70SpvPot1njU4jhLhHUtTFhGvL52itdvHTylUys54QDkaKurioEEF6+Ui6XJ3H0gPnjU4jhLgH+S5qpZRZKfWrUmqxLQMJ23Ft+RzdzBv5buUWuQqMEA7kXraonwUO2SqIKALVOqBLVqZ54jzWHb1odBohRD7lq6iVUsFAF+AL28YRNmUyYWkxlgGWVXy5crfRaYQQ+eSSz+X+D3gR8LnbAkqp4cBwgPbt2xMbW7CrYSclJRXoec7IJmPhH0Wgqxc1Y+fx07YKhJfzLvzXsBH5bNxJxuMWZxiL4ODguz6WZ1ErpboC8VrrnUqpB++2nNZ6CjAFIDo6Wuf2onm5n+c6G5uMRYsxPLXqQ8btH0CXxjULf/02JJ+NO8l43OLMY5GfXR/NgEeUUjHAd0AbpdQ3Nk0lbKvhk3ibMih5fAEHziYbnUYIkYc8i1pr/bLWOlhrHQL0BlZprfvbPJmwHXdfzI2HMNbjFz5ffdToNEKIPMhx1MVVk5GUsV7g+sElnEm6bnQaIUQu7qmotdZrtNZdbRVGFCHfsqi6jzHGYynTNsiFBYSwZ7JFXYypyFGEZ+5l7/Z1XLqWbnQcIcRdSFEXZ2XDsIa0ZJhlCV9vkQsLCGGvpKiLOVPT0XTUG/h5405SM+TCAkLYIynq4q5qe/CvRE/rEubvKthJSkII25KiLu5MJkxRo+hnXsnXaw+SJVOgCmF3pKgFhPfGzeJCVMpSlh+UKVCFsDdS1AIsHpgaD2WUx3Kmrj0qU6AKYWekqEW2xsMolRlPqbOr2XHqktFphBC3kaIW2byDUOGPM85nBf9de9zoNEKI20hRi1siR1E9dS8XDm/hWHyK0WmEEDdJUYtbSodClTa86LeCqevktHIh7IUUtbhT1Giapa1n8697iE9JNTqNEAIpavFHVdqgAqrztOdKvt4sp5ULYQ+kqMWdlEJFjeJRvYL5mw9zPT3T6ERCFHtS1OLP6vbC4upOd1Yzb6ecVi6E0aSoxZ9Z3FGNhzHCdSnT1x2T08qFMJgUtfhrEUPwzUqi7tUNLDsgp5ULYSQpavHXvAJQ4b0Z57OCKetPGJ1GiGJNilrcXeTTVLy+D85sZ+epJKPTCFFsSVGLuwusAdU68Ir/Kqask61qIYwiRS1yFzWaRtfXc+Dgfk5evGZ0GiGKJSlqkbtKLSEolPElVjNtg2xVC2EEKWqRO6VQUaPpnLGcX3YcIfFqmtGJhCh2pKhF3uo8hou7D4M81vHNltNGpxGi2JGiFnlzcUU1HsaT5iXM2nRcrlYuRBGTohb502gwXllXaGXdzA+/xhmdRohiRYpa5I+nP6peX573WsbUdcexymnlQhQZKWqRf02eosy1QwRe3s2qw/FGpxGi2JCiFvkXUBVV4yFeK7lSTisXogjlWdRKKXel1Dal1B6l1AGl1N+LIpiwU1GjqJOygfMxh9hz5rLRaYQoFvKzRZ0GtNFahwP1gE5KqUibphL2q2IzVNlwXiu1lqmyVS1EkcizqHW2qzfvWm7e5Juk4kopiBpN29RlbNh3lDNJ141OJITTy9c+aqWUWSm1G4gHlmutt9o0lbBvtbtj8ijBmBKbmL5RrlYuhK255GchrXUWUE8pVQL4QSlVR2u9//ZllFLDgeEA7du3Jza2YJdwSkqS6TR/Z89j4V2zN3/bM5Oobe3oVcsLX/d8fZTuiz2PhxFkPG5xhrEIDg6+62P39K9La31ZKbUa6ATs/8NjU4ApANHR0Tq3F83L/TzX2djtWJR6Fr37v/Rw286a2Fo8/WBIkbys3Y6HQWQ8bnHmscjPUR+BN7ekUUp5AO2BwzbOJeydRwlUgwGM9ljKjA0nSc+0Gp1ICKeVn33UZYHVSqm9wHay91Evtm0s4RCajCQg5TdCM/azaM9Zo9MI4bTy3PWhtd4L1C+CLMLR+FdC1erKK+dW8My6hjzWoDxKKaNTCeF05MxEcX+iRlP98gayEo6y7uhFo9MI4ZSkqMX9qdAEyjciOmgtU+W6ikLYhBS1uD9KQdQomqYsZf+xkxw4m2x0IiGcjhS1uH+1HsHkHcjLQZv4Yr2cACNEYZOiFvfP7AJNRtI9/WeW7jnFueQbRicSwqlIUYvC0eAJXHUaQ323MWNjjNFphHAqUtSicLj7oiKGMtS8mO+2xpCSmmF0IiGchhS1KDyRT+GTeo6OLjuYs/2M0WmEcBpS1KLweAeh6vfjBc+fmb7+BBlZclq5EIVBiloUrqbPEHT1MNVv7ObnfeeMTiOEU5CiFoXLvzIqtDsTSyxh6voTaC3XmBDifklRi8LXfCxVUrZjOrebzScSjU4jhMOTohaFr2w4VGnDJP8Vclq5EIVAilrYRvPnaHBtHSeP7OXIhRSj0wjh0KSohW2EtECVq8+rJVbwhVytXIj7IkUtbEMpaDaWtmkr2PjrAeJTUo1OJITDkqIWtlOzK6pkCKM8lzNz0ymj0wjhsKSohe2YTKhmY+ipl7Fg80Gup2canUgIhyRFLWwr7G9YPHzpo5Ywd0es0WmEcEhS1MK2XNxQLZ5niPlnZq0/QKacVi7EPZOiFrZXfwCubh50vPajXK1ciAKQoha2Z3HH1OI5Rrr+wrRV+8myymnlQtwLKWpRNBoMxN3VQovkH/llv0zWJMS9kKIWRcPVE3PzZxnt9jNTV+7HKlvVQuSbFLUoOo0G42GBiKQfWX7ogtFphHAYUtSi6Lh6YW76DGPcfuLzlQdkClQh8kmKWhStxsPwcrESFv8ja35LMDqNEA5BiloULTcfzFGjeM59MZ+uPChb1ULkgxS1KHpNhuNjzqT22e/ZdFwuLCBEXqSoRdFz98PcYiwvuC/i8+X7ZKtaiDzkWdRKqQpKqdVKqYNKqQNKqWeLIphwco2H4+FmoW7cd2w8JlvVQuQmP1vUmcALWutQIBIYpZQKtW0s4fRcPXFpNZ7Rrj/x2ZKdslUtRC7yLGqt9Tmt9a6bf04BDgHlbR1MFAMNB+Lq5Udk/HesOhxvdBoh7JbLvSyslAoB6gNb/+Kx4cBwgPbt2xMbW7ApLZOSkgr0PGdUHMbCs/5Ihq3/J/0XP0w17whMSt112eIwHvdCxuMWZxiL4ODguz6W76JWSnkD84GxWusrf3xcaz0FmAIQHR2tc3vRvNzPc52N049F2afI3DeDzonz2H85is51y+a6uNOPxz2S8bjFmcciX0d9KKUsZJf0LK3197aNJIoVswsubV/lCfNyZi7dJDPrCfEX8nPUhwKmAYe01h/YPpIodkK7owNr0OvKTBbujjM6jRB2Jz9b1M2AAUAbpdTum7fONs4lihOTCUunf/KoaS2Lli0nQ64CI8Qd8nPUxwattdJah2mt6928/VwU4UQxUrkVmZXbMSz1S77ddtroNELYFTkzUdgNS6c3iWIf25bPJSU1w+g4QtgNKWphP4JqYq0/gLH6a6asOWp0GiHshhS1sCsubV4lxBRP4qYZnE9ONTqOEHZBilrYF5/SmJqP5QXzXD7+5Vej0whhF6Sohd0xNRuDl6c75fd/wqFzfzq3SohiR4pa2B9XT9w7v8Uwl1+YvnCZTNjkzKxZN29WkL/nu7qnuT6EKDKh3UgvP4WHz/yH5Qea0aFO7qeWCzt2PQlid0D8QYg/BEnH4VoCXLsI6VdvLadM4BkA3kHgUwYCa0JQKJSpA6XrgMls3HswmBS1sE9K4dntfZp/2oxXF82gZY0XjU4k8suaBTEb4MhSiFkH5/eB2Q0Cq0NQbajWMbuMvQLA3Q+4ORFXVnp2eV+Lh+S47GLfNw+ungc3P6jYFCq1hJqdoWSIke+wyElRC/sVVIuMiGGM2j6d6au78khoSaMTibvRGuJ2wp5v4eBC9PUk0ss34UxAO3aXeY5dmZU4k5zB+ZhUUg5ncj09k9QMKyZTCq5mE24WMyU9LQT6VCDQuyoPlPKiWpg31Up7U8n9Km6xm+Hketj6OSx9GcrWg9BuEN4bfMsZ/e5tTopa2DW3tq9Qau9crBs+5HyF13De+dEcVPp12D8Ptn8B5/Zw3j+C1V79+eJGHY4f88DNxUTNMj48UEoRWtaXtjWD8PWw4Olqxt1ixqo16ZlW0jKtXLqWTsLVNC5cSWPz8Yt8tSmG5BsZuJgUtcuVpmHF4TRq8yJR3mcpefIX2DUTVr0B1TvhXrErlOvttLtHpKiFfXP3w+2hfzFiwSgmrWlHo1qVjU4kANJSYPs0sjZOJiMjnZ/Mbfko7QlSr1aiadVSPNUkgPBgPyoFeOFiLtgxC1prLl5N59C5K+w8dYmdpy4xZ/tprqVnERbcmta1evGw3wmqnJ5HqRVjYOvbEDUKGgwAV69CfsPGkqIWds8U/jeubp9NjzP/x6aj7WhaLcjoSMVX2lUyN3+KdePHXMsy8XFaZ7b6d6NNWCU+r1uGGqV9ULlc/OFeKKUI9HEj0CeQltUDAcjMsrInNpk1v8Wz8rd4/hOnCfDuQ8eKjzPSbyvB695FrX0bIoZB4+HgHVgoWYwmRS3sn1L49vyIsMkRfDLvfRqM+xfuFuf8FdduWbNI2TITtfoNrmdoplofJT28P49HVWNiOb8ii+FiNtGwYkkaVizJCx1qEH8llRWH4pm//SStDjeivFcjXir7K+1+nYPbpsmoBgOh2Rjwc+ydZlLUwjGUrMjlBs8wfOenTP+lC08/0tLoRMXGxX3Lyfz5ZXyvn+Y7t8fwaP8szzaqireb8fUR5OtO3yYP0LK8CTe/QJbsP8fMvf48ezGcft6/Mnr/QgJ3TEfV6wvNx4K/Y+46M36khcinjPABZJ5YQpUd/+Bgo+8JLedrdCSnlhgfx5lvnyMsaRnL3dpi6vwFAyPCMZsKZ9dGYQv0cWNAVAgDokI4l3yDxXtqM3j3g5S+vIbxe36k+q6vuVGjB15tX4SgmkbHvSdyZqJwHCYX/B7/jHamnfz47Sdy2S4bSU3PZOW3H2D+pDF+yYfY3mY27SfMo32TenZb0n9U1s+DYS0rs3hMS14e+zw/N/mG8e6vs+/QQayfRnLqs8e4fHyH0THzTbaohUNRZcO4Efk8Izd/xJxVrenbLtLoSE5Da83yzdvxX/48zfVhDlYfSd1er1HJ1d3oaPelapA3z3eogW5fnb2xg5i5/mdqHP0vFb9uy26PJiQ3GkPD5p3sYlfO3dhvMiHuwrv9BC4eWkKF9S8SE/4zIYHeRkdyeCcSrvLTtx8zMPE/JPtWJ6PPBuqXc6zdA3lRShFeoQThffuSZe3Dnm2rcdn4Aa3W92PzutpsrTCEkAYdaVYtkEAfN6Pj3kGKWjgeswX/fl/S+LPmfPnVvxj63JsFPla3uEvNyGLayr2U2TSJp03rSY58gQodX3baE0d+ZzYpwiPbQGQb0uL2EbLkLZqcGcfe2M+JTu9EbFArIqqWo94DJQgrX4IK/h75O+xQayikwxNvJ0UtHJIpqDqpD77OwNV/59ufWzPg4fZGR3I4u05fYuq383j5xnuU8jZj6r0E/weaGB2ryLmVr0vZIbPh4jHCNnzI5P1fkJUyjS0Ho5i/ox6vpFZHu5Wggr8nwSU9KFfCI+fMykyr5lpqOn6Jv1L7wo+UMl2j3vifCj2jFLVwWH4tnubCwSU02v48++oup25IGaMjOYSMLCsfrfiNjPX/4WOX/5EV2gPXRz64OUFSMRZQFVP3T6DzO5h/+4WW++bS4sQn4J5Kil914l0rcfZaEAmXXcnI0rhlXCE4K5ZqmUfwzUoi1q8hF0OftEk0KWrhuEwmSj/xFW4fNmbDrGeoMv5bPF3lI52bY/EpvDF7JU9ffpeG7icxP/wZ5vC/GR3Lvrh6Qd2eULcnKjMd4nbie2YLvhePUjU5Jvv0eQD3EuBfC8r2gWrtecC3HA/YKJJ8qoVj8yqFy+Nf0Wn2w3w740P6DxtXaKcwOxOrVfPV5hi2L/mGjy1TcC9TFZfHNzjsCSBFxsUVKkZl3wwk38AIh+ddvTnxES/xaNy7/LRqrdFx7M655BsM/WIdHkvH87H5A7ybDccybJmUtAORLWrhFMo99CJnYjZTZ90IDoUso1aVikZHsgsLd8fx9YLFfGj+iLK+mZh6/gghzY2OJe6RbFEL52AyETzka1zdPLg2qz+XU64ZnchQl6+n88zsXeyd9zZz1KuUr1Yfl6c3SUk7KClq4TSUuy8lhsynmj7J1s9Hkp5pNTqSIdYfTaD3hz/S/8Q4XnWfi7nr+5j+NhM8/Y2OJgpIilo4Fc/SVUjrMYO2137il2nRxeoK5jfSs4hedIAvvpzKPD2OiFIZmEasg4YDbXIShig6UtTC6QSFteNMi3foenYyK+d+anScIrHnzGW6T15FjT1vMcP1Hbwj+mIatiL7grLC4eX5ZaJSajrQFYjXWtexfSQh7l+ltkM5cDmeVnsnsnlZAFEdnPNY4YwsK5+sPsbS1av40vu/lHG/iuoxH6q2NTqaKET52aKeAXSycQ4hCl3tx17hQKUnqbdxNNvWLjY6TqE7mZRKz083kLHxExa7TaRcpVBMT2+WknZCeW5Ra63XKaVCiiCLEIWu3sD32TvlCrVXDWaP+oLwlo8YHem+Wa2a6RtP8tWSLXzuM41aLocxPfQ+1B8g+6KdVKEdR62UGg4MB2jfvj2xsbEFWk9SUlJhRXJ4MhZ3Kuh4lHzodX5dYKXRysFsuvxvHmjQsZCTFZ245DTeWnGaSomrWeY+DZNPZS48OI8sv4oQF2d0PMM4w7+V4OC7X9ex0Ipaaz0FmAIQHR2tc3vRvNzPc52NjMWdCjoewaO/YOPU52i8czwHLCbqPzSokJPZVpZVM3NzDDOWbuNdz5lEuGzjSvgI/Lr8g7JmOW8NnPvfivwNi2JBmUw0G/5/bJ7pS5Mtz7H70knq9fm7Q+wqOHj2Ci/P30PYxcUsd52FpWQVVLc1pGSWxE9KuliQv2VRbCilaDrwDTYsrEijXRPY/dFRag//Aou7l9HR/tKN9Cz+b+URVq3fwGS/2dR0OYhq/QpEjgKzCxRw96JwPHke9aGU+hbYDNRQSsUqpYbYPpYQttO821COdP4fZZK2cva9ZiSe2m90pDtYrZoffo3l4fd+4oHt/2Sp2wRqlfVFPbURmj2bXdKiWMnPUR99iiKIEEUprEkb4oLXcmHGIAK/bMuhRq9Qq8sYw3eF7Dx1iX/+uJfaFxax0P17PHy8MXX6Cmp2MTybMI78r1kUW+XLVyDwxSWs+/oNmm1/k+OHfsC/92eUrFCryLPsjb3M5OWH8T36A1O8FuHvcQVTs7HQdDRYPIo8j7AvUtSiWHO1uNBu8N/Zs6c71kVjqDCtOQcq9qHaY9G4+gbY9LW11mw+kchXaw/jd3wBb3guobTnRUxNRkLTMTKJksghRS0EEB5en4zaa1izYBrV971P2gfz+a1yX6p1fR53//KF+lpXUjNYvOccSzZsofGln3jPdTWeXhpzoyezC9o7sFBfTzg+KWohbrK4mGnfczjJHfuzdeFkqh6bgXnyl+wt2QafiD5UjHgYk8W1QOtOvpHB+qMJrN+5B48Ty+hm3khffiMjsCqWqIkQ3gfcvAv5HQlnIUUtxB/4+XjSrv8ErqU+z5YVc3DdN4saS4dzZZk7p3wakhrcDJ8qTShTJYwSfiUxmW59yae15sqNTOIu3+DIuSTiYw6SEvMrpS7toYV5H13VWdK8AnEJ6wnhn2MpV1++JBR5kqIW4i683F1p0XUAdB3AufPnObV5PvrEOkIOT6XsoX8BcF6X5Iry5YryJkubQGfirtMIUpfpqi7jgpUbLr6kVQjHO3Q4VG2LW+k6YJIZhkX+SVELkQ9ly5ShbI9RwCgArl6+SPyJPWQmHCf96kV80pIxYcXFYsHD0xuvgGDMARUgoDoevuXwkK1mcR+kqIUoAO8SAXg3aAvIlKLC9uT3LyGEsHNS1EIIYeekqIUQws5JUQshhJ2TohZCCDsnRS2EEHZOiloIIeycFLUQQtg5pbUu/JUq9QVQ0OsENQR2FmIcRyZjcScZjzvJeNziDGMRo7We8VcP2KSo74dSaofWupHROeyBjMWdZDzuJONxi7OPhez6EEIIOydFLYQQds4ei3qK0QHsiIzFnWQ87iTjcYtTj4Xd7aMWQghxJ3vcohZCCHEbKWohhLBzdlvUSqlnlFKHlVIHlFLvGJ3HaEqpF5RSWikVYHQWIyml3r35udirlPpBKVXC6ExFTSnVSSn1m1LqmFJqgtF5jKSUqqCUWq2UOnizK541OpMt2GVRK6VaA92AcK11beA9gyMZSilVAegAnDY6ix1YDtTRWocBR4CXDc5TpJRSZuAT4CEgFOijlAo1NpWhMoEXtNahQCQwyhnHwy6LGngKeFtrnQagtY43OI/RPgReBIr9N79a62Va68ybd7cAwUbmMUBj4JjW+oTWOh34juyNmmJJa31Oa73r5p9TgENAeWNTFT57LerqQAul1Fal1FqlVITRgYyilOoGxGmt9xidxQ4NBn4xOkQRKw+cue1+LE5YTAWhlAoB6gNbDY5S6Ay7uK1SagVQ5i8eepXsXP5k/yoTAfxPKVVZO+mxhHmMxStk7/YoNnIbD631wpvLvEr2r72zijKbsE9KKW9gPjBWa33F6DyFzbCi1lq3u9tjSqmngO9vFvM2pZQVCAASiipfUbrbWCil6gKVgD1KKcj+NX+XUqqx1vp8EUYsUrl9NgCUUk8CXYG2zvo/71zEARVuux9882fFllLKQnZJz9Jaf290Hluw110fC4DWAEqp6oArcNHIQEbQWu/TWgdprUO01iFk/5rbwJlLOi9KqU5k769/RGt93eg8BtgOVFNKVVJKuQK9gUUGZzKMyt6CmQYc0lp/YHQeW7HXop4OVFZK7Sf7y5KBxXDLSfy1jwEfYLlSardS6nOjAxWlm1+kjgaWkv3F2f+01geMTWWoZsAAoM3Nz8NupVRno0MVNjmFXAgh7Jy9blELIYS4SYpaCCHsnBS1EELYOSlqIYSwc1LUQghh56SohdNTSoXcPElGCIckh+cJp3bzLNcxgDdwHOhdnE8YEo5Jilo4LaWUD9nl3AkIA9YAicA3wHyt9Uyl1Aigpda6n2FBhciDYXN9CFEErGRPDesPoLWOAVBKDQc2KqVOAi+QPfmXEHZLilo4La31NaXUMOAtoIxSqg7wutb6glLqdWA10ENrnWRoUCHyIF8mCqemtV4E9ALeAQLJ3oIGqEv2bpByBkUTIt+kqIXTUkp5K6Uq3rz7+9U/fJRSjcm+lFV9YJxSqpJRGYXID9n1IZyZBfgvUIrs+cxPA0+QPS3oIK31WaXUC8B0pVQbmaFR2Cs56kM4vZuXaHpQaz3D4ChCFIjs+hDFwWVgt8EZhCgw2aIWQgg7J1vUQghh56SohRDCzklRCyGEnZOiFkIIOydFLYQQdu7/AexKv17wp6hbAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "n = 7\n", "#create grid and get function for computing coefs\n", "xhat, θ, x, fcheb, S = cheby_dct_setup(xa, xb, n)\n", "#sample the example function at the grid points\n", "y = fun(x)\n", "#compute expansion coefficients super easily!\n", "a = fcheb(y)\n", "#see how well the interpolation does\n", "yfun = fun(xx)\n", "yche = cheby_sum(xx, a, xa, xb) #<-- evaluates cheby expansion\n", "plt.plot(xx, yfun, label='example function')\n", "plt.plot(xx, yche, label='cheby interpolation')\n", "plt.xlabel('$x')\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "binary-target", "metadata": {}, "source": [ "Starting to look pretty good with 7 points. We can also directly test the convergence rate as the number of points increases." ] }, { "cell_type": "code", "execution_count": 5, "id": "published-illustration", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEGCAYAAACdJRn3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAro0lEQVR4nO3debxd0/3/8dc7kxDzVJpLE6RRVKtURJQYLjGmYox5DFVDW1UxVMO3KdVSXyWIr4ihhKqvBjFP8TXGLEFIU60bKjQ/SpBWfH5/rH3b47rDucO5+9y738/H4zzu2Wvvs/dnnZ3cz9177bWWIgIzM7PW6JF3AGZm1vU4eZiZWas5eZiZWas5eZiZWas5eZiZWasVInlMnjw5gADinXfe+ff7or2KXPei17/IdS96/dtZ9yYVInm8/vrr/36/aNGi/ALJWZHrDsWuf5HrDsWuf6XqXojkYWZmHcvJw8zMWs3Jw8zMWq1X3gG0haR+wATgn8CDEfG7nEMyMyuUqrnykDRJ0nxJMxuUj5A0W9IcSWOz4lHATRFxJLBbpwdrZlZwVZM8gMnAiNICST2Bi4EdgfWA0ZLWA2qAN7LNFndijGZmRhXdtoqI6ZIGNCjeFJgTEXMBJE0BRgJ1pATyHE0kQEljgDEAtbW11NXVAbBgwYIKRN81FLnuUOz6F7nuUOz6t6fuNTU1Ta6rmuTRhP785woDUtIYAlwIXCRpZ+DWxj4YEROBiQDjxo2L0i+huS/kCyZNglGjYPnlWxl6dWpV3buhIte/yHWHYte/EnWvpttWZYuIhRFxaER8r6KN5e+8A2eeCV/5CpxxBhT4rxczs1LVnjzmAWuULNdkZZ1jlVXgtdfg17+Ga66BAQPgtNPg3Xc7LQQzs2pU7cljBjBI0kBJfYB9gamdGkGfPnDkkfDqq3DBBTBlSkoiY8emKxMzswKqmuQh6XrgMWCwpDpJh0fEp8CxwF3Ay8CNETErlwB794bDDoPZs2HCBLj55pREfvxjePvtXEIyM8tL1SSPiBgdEatHRO+IqImIK7LyaRHx1YhYOyLG5x0nvXrBQQfBSy/BxIlw220wcCD88Ifw1lt5R2dm1imqJnl0Ob16wf77w6xZ6Ymsu+9OSeT442Fe5zXLmJnlwcmjvXr2hH33hRdfTI3qDzwAa60F3/8+/PWveUdnZlYRTh4dpUcP2GsveP55uP56eOQRWGcdOPpoKJlPxMysO3Dy6Gg9eqROhc8+C7//PcyYAYMGwRFHwNy5eUdnZtYhnDwqRYKRI+Gpp+CWW+CFF+CrX4VDD4U5c/KOzsysXZw8Kk2CnXeGJ55IT2a98goMHpye2Jo9O+/ozMzaxMmjs0gwYgQ8+ijceWe6hbXeeumJrZdfzjs6M7NWcfLobBLU1sLDD8M990BdHay/PhxyCMyfn3d0ZmZlcfLIiwTbbAMPPQT33QdPP51uZ116KSz2FCVmVt2cPKrB1lvDM8/AT38KJ50EQ4emZGJmVqWcPKpF797wox+l9o+vfAU23RSOPRbeey/vyMzMvsDJo9rU1KT+IbffnhrWBw+Ga6+FiLwjMzP7NyePajViBMycCccckzoYbr11GozRzKwKOHlUs7594Wc/S0mkb1/4xjfSPCILF+YdmZkVnJNHV7DOOnDHHWkiqmuvha99LfVa960sM8uJk0dXIcEee6QG9X32gT33hF139XhZZpYLJ4+uZpll4Fe/SgMvvv9+6mD485/DokV5R2ZmBdIlk4ek70q6XNINkrbPO55cfP3rMH166lT43/+dlu+5J++ozKwgOj15SJokab6kmQ3KR0iaLWmOpLHN7SMibomII4GjgX0qGW9Vk+Dgg9MAi9tuCzvskCamevPNvCMzs24ujyuPycCI0gJJPYGLgR2B9YDRktaT9HVJtzV4rVry0dOzzxXbiivCJZfA44/Da6/BuuvCBRfAp5/mHZmZdVO9OvuAETFd0oAGxZsCcyJiLoCkKcDIiDgb2KXhPiQJOAe4IyKeaew4ksYAYwBqa2upq6sDYMGCBR1Ukyr05S/DzTfT75prWO6MM/h04kTe+8Uv+OcmmwDdvO5lKHL9i1x3KHb921P3mpqaJtd1evJoQn/gjZLlOmBIM9sfB2wHLCdpnYi4tOEGETERmAgwbty4KP0SmvtCuoXTT4cjjqDPSSex6qhRcNhhcM45sOKK3b/uLShy/Ytcdyh2/StR9y7ZYB4RF0bExhFxdGOJw4DVVoNrrkkj9j7yCAwezFJTprhviJl1iGpJHvOANUqWa7Iya6+tt4bnn4eTTmKF00+HXXaBt97KOyoz6+KqJXnMAAZJGiipD7AvMDXnmLqPPn1g7FjevuMOePvt9FjvH/6Qd1Rm1oXl8aju9cBjwGBJdZIOj4hPgWOBu4CXgRsjYlZnx9bdfTpoEDz2WBpsce+90zzq77+fd1hm1gXl8bTV6CbKpwHTOjmc4undG846C3baCQ48EDbcEK66CoYPzzsyM+tCquW2lXW2zTaD555LSWSbbeDEE+GTT/KOysy6CCePIuvXL3UuvP12uO462GSTNGaWmVkLnDwMdtwRXnwxDfU+ZAicfTYsXpx3VGZWxZw8LFl5ZbjxRpg0KXUo3HJL+NOf8o7KzKqUk4f9hwQHHJCuQpZYIs1cePnl7lhoZl/g5GFftOaacO+98F//BccdB7vtlvqHmJllnDyscT16wA9/CE8/DfPmwQYbwP/+b95RmVmVcPKw5q2/fhrqfcyYNPXtoYfCP/6Rd1RmljMnD2tZnz4wfjw8/HCavXDDDdNPMyusZpOHpJ6SHuisYKzKbb55GmRxhx1Sj/Sf/MRzp5sVVLPJIyIWA59JWq6T4rFqt/TScNllMHUqXH01fPvb8MILeUdlZp2snNtWHwIvSrpC0oX1r0oHZlVul13SI73rrJN6pp97rjsWmhVIOQMj3py9zD5vlVXS0O5XX50e6b3ttjTI4sCBeUdmZhXW4pVHRFwFXA88nb2uy8rMUsfCgw9Ot6569EiN6dddl3dUZlZhLSYPScOB14CLgQnAq5K2rGxY1uUMGAD335/mTz/ooPRo78cf5x2VmVVIOW0e5wHbR8RWEbElsAPwm8qGZV1Sjx5w8snw4IMwbRoMHQqvvpp3VGZWAeUkj94RMbt+ISJeBXpXLiTr8rbYIg3tvtpqsPHGcMMNeUdkZh2snOTxtKT/kTQ8e10OPFXpwFoiqZ+kpyTtkncs1ohVVklXH6ecAvvtl6a+9WRTZt1GOcnjaOAl4Pjs9RLwvbYeUNIkSfMlzWxQPkLSbElzJI0tY1cnAze2NQ7rBD16wKmnpraQW26BYcM8zLtZN9Hso7qSegLPR8S6wPkddMzJwEXA1Q2OczFQC9QBMyRNBXoCZzf4/GHAN0hJrG8HxWSVtNVW6TbWAQfAt76V5gzZY4+8ozKzdmg2eUTE4uxqYM2I+GtHHDAipksa0KB4U2BORMwFkDQFGBkRZwNfuC2VPQHWD1gP+FjStIj4rCPiswr50pfgzjvTGFl77w3f/z786ldp3hAz63LK6SS4AjBL0pPAwvrCiNitA+PoD7xRslwHDGlq44g4DUDSIcC7jSUOSWOAMQC1tbXU1dUBsGDBgg4LuqupirofdhhLDB7Miscdx+Lp0/n7hAksXnPNTjl0VdQ/J0WuOxS7/u2pe01NTZPrykkeP23zkSssIiY3s24iMBFg3LhxUfolNPeFdHdVUfd99oGttqLn6NGsvvPOMHkyjBzZKYeuivrnpMh1h2LXvxJ1L6fN47KszaOS5gFrlCzXZGXWXa22Wpqt8MwzYdQoOOGENHd6nz55R2ZmZShnVN3Zkip9X2EGMEjSQEl9gH2BqRU+puWtZ0846yy44w649lrYckv4y1/yjsrMylDOo7r1bR73SZpa/2rrASVdDzwGDJZUJ+nwiPgUOBa4C3gZuDEiZrX1GNbFbL89PPdcuurYaKM0wKKZVbVOb/OIiNFNlE8DpnXksawL+fKXU3+QM85I7R8nnpiezOrtwQzMqlGTyUPSuhHxSkQ8JGmJiFhUsm6zzgnPCqVXL/jFL+A734EDD4RHHklDmxS4odOsWjV326p0XO3HGqybUIFYzJIdd0ydCgG++c3UJmJmVaW55KEm3je2bNax1lgjjc572GFp1sJTT4VPP807KjPLNJc8oon3jS2bdbzevdP0tn/8I1x6KWy7Lbz5Zt5RmRnNJ4+abL7y35a8r1/u30nxmaUrj2efhUWL0m2se+7JOyKzwmvuaauTSt43HII99yHZrWC+8hWYPj0N8T5iRHoq66c/TSP3mlmnazJ5eJ5yqzp9+sB556XJpg46CGbOhKuugqWWyjsys8Lxn23W9ey+e3qMd8aM9FhvNuilmXUeJw/rmjbcEJ58EpZcEr79bXjiibwjMisUJw/rulZdFe67L7WBbLUVXHddy58xsw7RYvKQ9NVsXKuZ2fKGkk6vfGhmZVhiiTQz4fjxqVf66afDZ54XzKzSyrnyuBw4BfgXQES8QBr11qw6SGksrKlT4cILYc894cMP847KrFsrJ3ksFRFPNihzV1+rPjvvDI89lkbo3WIL+GuHzJxsZo0oJ3m8K2ltsl7lkvYE3qpoVGZttf76qSF9ueVg001TMjGzDldO8vg+cBmwrqR5wA+AoysZlFm7rLxy6oW+yy4wfDhcc03eEZl1O+XM5/GXiNhOUj+gR0R8UOmgzNqtTx+4/HLYYAM45BCYNQuOOSbvqMy6jXKuPP4saSKwGeBWSOs6JPjBD9LMhJdcwkpHHAEf+G8fs45QTvJYF7iXdPvqz5IukrRFZcMy60A77giPP07vV1+FYcM8T7pZB2gxeUTERxFxY0SMAjYClgUeqnhkzZDUQ9J4Sb+VdHCesVgX8bWvMf/WW2GllVKP9EceyTsisy6trB7mkraSNAF4GugL7N3WA0qaJGl+fafDkvIRkmZLmiNpbAu7GQnUkPqeeGAjK8tnK6wAd98No0bB1lvD5Ml5h2TWZbXYYC7pdeBZ4EbgpIhY2M5jTgYuAq4uOUZP4GKglpQMZkiaCvQEzm7w+cOAwcCjEXGZpJuA+9oZkxVF795wySXpkd4jjkgN6eecAz175h2ZWZdSztNWG0bEPzrqgBExXdKABsWbAnMiYi6ApCnAyIg4G9il4T4k1QH/zBYXN3YcSWOAMQC1tbXUZSOvLliwoANq0TUVue7QoP67784SK67ISsccw6LnnmPBhRcSyyyTX3AV5nNf3Pq3p+41NTVNrmsyeUj6SUScC4yX9IVpZyPi+DZH9EX9gTdKluuAIc1sfzPwW0nfAaY3tkFETAQmAowbNy5Kv4TmvpDursh1hwb1339/2GQTltx1V/rvtVca3mSttfILrsJ87otb/0rUvbkrj5ezn1U3a2BEfAQcnncc1g0MHgyPPw577516pN98M2y5Zd5RmVW9JhvMI+LW7O1HEXFV6Qv4qIPjmAesUbJck5WZVd6KK8Idd8A++8B228EVV+QdkVnVK+dpq1PKLGuPGcAgSQMl9SGN2ju1g49h1rTeveHii+GCC+Coo+BHP4LFjTanmRnNt3nsCOwE9Jd0YcmqZWnHqLqSrgeGAytnDd8/i4grJB0L3EV6wmpSRMxq6zHM2uyYY9KtrL32gldegRtugG7ckG7WVs21ebxJau/YjdS/o94HwA/besCIGN1E+TRgWlv3a9Zhtt02TWu7006pP8jtt8OXvpR3VGZVpcnkERHPA89Lui4i/tWJMZnlb9Cg1At9l11g883hrrtgnXXyjsqsapTT5jFA0k2SXpI0t/5V8cjM8rbqqnD//ek21uabw4wZeUdkVjXKSR5XApeQ2jm2JvUMv7aSQZlVjaWXhj/+8T+3sO68M++IzKpCOcljyYi4D1BE/CUixgE7VzYssyrSuzdceSWccALsuitcfXXLnzHr5soZnmSRpB7Aa9kTUfOApSsbllmVkWD8eFh9dTj0UHjzTTj55FRuVkDlJI8TgKWA44H/ArYBPAy6FdOxx8Jqq8EBB8C8ealfiAdVtAJqMXlERH0r4YfAoZUNx6wL2HPP1Ji+227wt7+lOdL79s07KrNO1VwnwVuBLwyIWC8idqtIRGZdwZZbwsMPp1kKR4yAW26B5ZfPOyqzTtPclcevOy0Ks67o61+Hxx6DHXaA73wnPYnVv3/eUZl1iuY6Cf57qllJSwJrRsTsTonKrKtYYw34v/9Lt7CGDk0JZL318o7KrOJafFRX0q7Ac8Cd2fI3s1n+zAzSqLz33AMbbwxbbOH50a0QyunnMY400997ABHxHDCwYhGZdUVLLgk33fSfYd1vuSXviMwqqpzk8a+IeL9BWZMN6WaF1bMnTJgAp52Wnsi67LK8IzKrmHL6ecyStB/QU9IgUn+PRysbllkXJcHpp6fOhEcdlfqCnHmmOxNat1POlcdxwPrAIuB64H1Sx0Eza8rhh6cxsX79azjySPi0zVPgmFWlFpNHRHwUEadFxLcjYhPgGuCiyodm1sXtvHMalfeWW2D33eGjjp692Sw/TSYPSRtKulvSTEk/l7S6pD8A9wEvdV6IZl3YZpulp69mzkyTTL37bt4RmXWI5q48LgeuA/YA3iU9rvsnYJ2I+E3lQ2uapDUl3SJpkqSxecZi1qLBg+HRR+Hjj2HYMHj99bwjMmu35pLHEhExOSJmR8QFwMKI+ElEfNKeA2a/8OdLmtmgfISk2ZLmlJEQvg7cFBGHARu1Jx6zTrH66jB9eupUOHQoPPdc3hGZtUtzyaOvpI0kfUvSt0hDs5cut9VkYERpgaSewMXAjsB6wGhJ60n6uqTbGrxWBR4HDpd0P1nnRbOqt+yyMG1amlRqyy1Te4hZF9Xco7pvAeeXLP+tZDlIQ7O3WkRMlzSgQfGmwJyImAsgaQowMiLOBnZpuA9JPwZ+lu3rJtJshw23GQOMAaitraWurg6ABQsWtCXsbqHIdYcqqv8557Dc0kuz9IgRLPjNb/h45MiKH7Jq6p6TIte/PXWvqalpcl1zY1tt3eYjtl5/4I2S5TpgSDPb3wmMy/qfvN7YBhExEZgIMG7cuCj9Epr7Qrq7Itcdqqj+EyfC4MGsdPzxEJHmCamwqql7Topc/0rUvZxOglUnImYCe+Ydh1m7nHhimhfk0EPhgw/glFPyjsisbNWSPOYBa5Qs12RlZt3bgQdCv36w777w4Yfw85+7N7p1CdWSPGYAgyQNJCWNfYH98g3JrJOMGpV6o48aBQsXwm9+4wRiVa+s5CFpQ2BA6fYRcXNbDijpemA4sLKkOlLD9xWSjgXuAnoCkyJiVlv2b9Yl7bgj3HEH7LprugK57DLPjW5VrcXkIWkSsCEwC/gsKw6gTckjIkY3UT4NmNaWfZp1C8OHw733pmltFy6Eq6+G3r3zjsqsUeVceWwWEZ4azawzDBkCDz4ItbVpWPcbboC+ffOOyuwLyhlV9zFJTh5mneUb30i90Z9+Ot3GWrgw74jMvqCc5HE1KYHMlvSCpBclvVDpwMwKbd114eGH4U9/gh12gPcbzsdmlq9ybltdARwIvMh/2jzMrNIGDkwJZLvt0oi8d90FK62Ud1RmQHlXHu9ExNSI+HNE/KX+VfHIzAz694eHHkqTSQ0fDn/7W94RmQHlJY9nJV0nabSkUfWvikdmZsmqq8IDD6TOhN/5Dvz1r3lHZFZW8liSNAXt9sCu2esLgxWaWQWtsALccw/U1KQE8tpreUdkBddim0dEHNoZgZhZC5ZZJg3pvsceaUj3e++F9dfPOyorqHI6CV5J6hT4OdlETGbWmZZcMs2Jvt9+sNVWqRF9443zjsoKqJzbVrcBt2ev+4BlgQ8rGZSZNaNPH5gyBXbeGbbZJs2RbtbJyrlt9YfS5Wxsqv+rWERm1rJeveDKK2GppWD77dPAitttl3dUViBtGVV3ELBqRwdiZq3UowdMmABLL52uQm66KfVIN+sE5bR5fEBq81D282/AyRWOy8zKIcG556bG9FGj4NprYZ998o7KCqCc21bLdEYgZtZGEpxxRuoHst9+aSysw/w8i1VWp8/nYWYVcuKJ6RbWkUemBHLccXlHZN1Yp8/nYWYVdNRR6QrkkEPSpFKeF90qxPN5mHU3BxyQnsIqnRfdrINV/XwektaSdIWkm0rK+km6StLlkvbPKzazqjVqFEydCuefDz/4AXzmAbGtY1V0Pg9JkyTNlzSzQfmIbH9zJI1tbh8RMTciDm9QPAq4KSKOBHYrJxazwhkxAu68EyZNYvmxY51ArENVej6PycBFpAQEgKSewMVALVAHzJA0FegJnN3g84dFxPxG9luTxQOwuJUxmRXHVlvBvfeyVG0tjBkDEyem/iFm7VRO8ngnIqa2ZecRMV3SgAbFmwJzImIugKQpwMiIOJvyR+utIyWQ52ji6knSGGAMQG1tLXV1dQAsWLCgdZXoRopcdyhw/fv35+MJE1j7mGP4eOFC/t8vf1m4BFLYc0/76l5TU9PkunKSx7OSrgNuJQ3NDrTrUd3+wBsly3XAkKY2lrQSMB7YSNIpWZK5GbhI0s5ZXF8QEROBiQDjxo2L0i+huS+kuyty3aG49a8bPpwe99xDv+23p1+/foW8AinquYfK1L2c5FE6n0e9TntUNyL+DhzdoGwh4KHizVpjyBC4++40FpYEl11WuARiHSeP+TzmAWuULNdkZWZWaUOGpGHcd9ghLTuBWBs1mTwk/SQizpX0Wxqfz+P4Nh5zBjBI0kBS0tgX2K+N+zKz1tpss5RAts9uJjiBWBs0d+XxcvbzqbbuPBu+fTiwsqQ64GcRcYWkY4G7SE9YTYqIWW09hpm1wWabff4W1qWXOoFYqzSZPCKiviH6hoj4pHSdpJXL2XlEjG6ifBowrdwgzawCShMIOIFYq5TzL+VJSZvVL0jaA3i0ciGZWaepTyBTpsDRR7sjoZWtnKet9gcmSXoQ+DKwErBNJYMys05U3wayww7pFtYll/gKxFpUztNWL0oaD1wDfABsGRF1FY/MzDrP0KGffwrLCcRaUM6Q7FcAa5OGZf8qcJuk30bExZUOzsw6kROItUI5t61eBI6IiAD+LGkIcH5lwzKzXAwdmgZTHDEiLTuBWBPKuW11QYPl94GGo9yaWXex+eb/SSASTJjgBGJfUM5tq0Gk0W7XA/rWl0fEWhWMy8zyVJ9A6m9hOYFYA+X8a7gSuAT4FNiaNLz6tZUMysyqwOabpzaQ3/0OjjnGj/Ha55STPJaMiPsARcRfImIcsHNlwzKzqlCaQL7/fScQ+7dyGswXSeoBvJYNKzIPWLqyYZlZ1ahPIPW3sC6+2LewrKwrjxOApYDjgY1JswoeXMmgzKzK1LeBXHutr0AMKO9pqxnZ2w/xHBpmxTVs2OefwrroIl+BFFhzQ7I3O/VsROzW8eGYWVUrTSDgBFJgzV15DCVNF3s98ASgTonIzKrbsGFwxx2w445p+eKL05WIFUpzyWM1oBYYTZqs6Xbges+9YWZssUVKIPVXIE4ghdPk9WZELI6IOyPiYGAzYA7wYPbElZkV3RZbpFtYV1+dGtHjCxOOWjfWbIO5pCVIfTpGAwOAC4H/rXxYZtYl1CcQX4EUTnMN5lcDG5Bm/DszImZ2WlRfjGUt4DRguYjYMyv7LimxLQtcERF35xWfWaGVJpAlloDzz3cCKYDmHpM4ABhE6ufxqKR/ZK8PJP2j3ANImiRpvqSZDcpHSJotaY6ksc3tIyLmRsThDcpuiYgjgaOBfcqNx8wqYIstYOrUNArvGWfkHY11gubmMO+o5+8mAxeRxsQCQFJP4GJSg3wdMCN7NLgnaRDGUodFxPxm9n96ti8zy9M228DNN8N3vwv9+sHYZv8mtC6unOFJ2iUipksa0KB4U2BORMwFkDQFGBkRZwO7lLNfSQLOAe6IiGc6MGQza6uddoLrr4d99kkJ5Ljj8o7IKqTiyaMJ/Ul9SOrVAUOa2ljSSsB4YCNJp2RJ5jhgO2A5SetExKUNPjMGGANQW1tLXV2aOXfBggUdWY8upch1h2LXv1PrPmQIS513Hiv84Af8v0WL+GjffTvv2E3wuW+bmpqaJtfllTxaJSL+TmrbKC27kPT0V1OfmQhMBBg3blyUfgnNfSHdXZHrDsWuf6fW/YQToG9fVjzmGFbs3x9Gj+68YzfB575j5ZU85gFrlCzXZGVm1l0cdRR89BEceCAstRSMHJl3RNaB8koeM4BBkgaSksa+pF7sZtad/PCH8OGHsPfe6Wms+mHdrcur+Ihmkq4HHgMGS6qTdHhEfAocC9wFvAzc6GFPzLqp009PSWT33WH69LyjsQ7SGU9bNXqzMyKmkTogmll3JsHZZ8PChbDzznDffbDppnlHZe3UJRrMzayLk+C//zslkB12gAcfhG98I++orB08EL+ZdY4ePeDyy9MwJrW18PLLeUdk7eDkYWadp2fPNArv0KGw3XYwd27eEVkbOXmYWefq3RtuuAHWXx+23RbeeKPlz1jVcfIws87Xty/ccgussUa6Ann77bwjslZy8jCzfCy1FNx2Gyy3XEogf/973hFZKzh5mFl+ll02zQXSo0dqSH///bwjsjI5eZhZvlZcEe65Bz74IPUDWbgw74isDE4eZpa/VVdNnQfffDPNB/LJJ3lHZC1w8jCz6tC/f0ogL7+cxsL617/yjsia4eRhZtVj4MCUQJ54Ag44ABYvzjsia4KTh5lVl8GDUxvIPffAEUfAZ5/lHZE1wmNbmVn12XBDuOuu1IlwqaXgoovS+FhWNZw8zKw6ffvbcPvtaSDFfv3gl790AqkiTh5mVr2+8x344x9hl11g6aXhjDPyjsgyTh5mVt1qa+H3v4c99khXICeemHdEhhvMzawr2G03uPZa+MlP0qyE77yTd0SF5+RhZl3DPvukwRSnTYO11oKf/hTeey/vqAqr6pOHpLUkXSHppgbl/SQ9JWmXvGIzs062664waxZceCFcc03qFzJ+fBraxDpVRZOHpEmS5kua2aB8hKTZkuZIGtvcPiJibkQc3siqk4EbOzJeM+sCevWCQw+FV1+FX/wCJkxIVyLnnQcff5x3dIVR6SuPycCI0gJJPYGLgR2B9YDRktaT9HVJtzV4rdrYTiXVAi8B8ysbvplVrT594Hvfgzlz4NRT06O8a6+dksmiRXlH1+1V9GmriJguaUCD4k2BORExF0DSFGBkRJwNlHsLajjQj5R8PpY0LSI+1w1V0hhgDEBtbS11dXUALFiwoG2V6QaKXHcodv27fd332gvttBNLX3kly5x6Kp+dfTb/OOEEPtpzT+jVq/vXvxntqXtNTU2T6/J4VLc/UDrvZB0wpKmNJa0EjAc2knRKRJwdEadl6w4B3m2YOAAiYiIwEWDcuHFR+iU094V0d0WuOxS7/oWo+znnwNix9Dj/fFY880xWnDgRxo2DYcOKUf8mVKLuVd9gHhF/j4ijI2Lt7OqkdN3kiLgtr9jMrAotvzycdRb8+c9pePcjjuBL228PN98MEXlH123kkTzmAWuULNdkZWZmHWflleHcc+FPf+KTYcNg9GjYZJP0qK+TSLvlkTxmAIMkDZTUB9gXmJpDHGZWBKuvzvtnnQWvvQYbb5w6HA4bBvffn3dkXVqlH9W9HngMGCypTtLhEfEpcCxwF/AycGNEzKpkHGZmrLkmTJwIr7ySnsrabrs0au+jj+YdWZdU0eQREaMjYvWI6B0RNRFxRVY+LSK+mrVjjK9kDGZmn7POOqmD4cyZaf70YcPS3OnPPJN3ZF2KB0Y0s2Jab7004OKzz6bRejfeGEaNgjPPhA02qOyxP/sszZL4r3+l1z//+fmfjZU1t66Z7fsOHgxHHtnhVXDyMLNi22gjuPVWePzxNF7WhhvCllumToiffpp+ybf0Kne7+m3bolcv6N07xVX6s4Wynqus0rHfV304FdmrmVlXs9lmaerbBx+EO+6AHj2gZ8/Pv3r1+mJZa9Y33KalJFC6ro0TYS2sq2OFjv2mACcPM7PPGz48vaxZVd9J0MzMqo+Th5mZtZqTh5mZtZqTh5mZtZqTh5mZtZqTh5mZtZqTh5mZtZqTh5mZtZqiAOPaS/of0oyFABsDT+cYTp6KXHcodv2LXHcodv3bU/fXI2JyYysKkTxKSXoqIjbJO448FLnuUOz6F7nuUOz6V6ruvm1lZmat5uRhZmatVsTkMTHvAHJU5LpDsetf5LpDsetfkboXrs3DzMzar4hXHmZm1k5OHmZm1mqFSh6SXpf0oqTnJD2VdzyVJGmSpPmSZpaUrSjpHkmvZT8rMcFY7pqo+zhJ87Jz/5yknfKMsZIkrSHpAUkvSZol6YSsvNuf/2bq3u3Pv6S+kp6U9HxW9zOz8oGSnpA0R9INkvp0yPGK1OYh6XVgk4h4N+9YKk3SlsCHwNURsUFWdi6wICLOkTQWWCEiTs4zzkpoou7jgA8j4td5xtYZJK0OrB4Rz0hahtRB7LvAIXTz899M3femm59/SQL6RcSHknoD/wecAPwIuDkipki6FHg+Ii5p7/EKdeVRJBExHVjQoHgkcFX2/irSf6pup4m6F0ZEvBURz2TvPwBeBvpTgPPfTN27vUg+zBZ7Z68AtgFuyso77LwXLXkEcLekpyWNyTuYHHwpIt7K3v8N+FKeweTgWEkvZLe1ut0tm8ZIGgBsBDxBwc5/g7pDAc6/pJ6SngPmA/cAfwLei4hPs03q6KBkWrTksUVEfAvYEfh+dnujkCLdryzOPUu4BFgb+CbwFnBertF0AklLA38AfhAR/yhd193PfyN1L8T5j4jFEfFNoAbYFFi3UscqVPKIiHnZz/nA/5K+3CJ5O7snXH9veH7O8XSaiHg7+4/1GXA53fzcZ/e8/wD8LiJuzooLcf4bq3vRzn9EvAc8AAwFlpfUK1tVA8zriGMUJnlI6pc1oCGpH7A9MLP5T3U7U4GDs/cHA3/MMZZOVf9LM7M73fjcZw2nVwAvR8T5Jau6/flvqu5FOP+SVpG0fPZ+SaCW1ObzALBntlmHnffCPG0laS3S1QZAL+C6iBifY0gVJel6YDiwMvA28DPgFuBGYE3gL8DeEdHtGpabqPtw0i2LAF4Hjiq5/9+tSNoCeBh4EfgsKz6VdO+/W5//Zuo+mm5+/iVtSGoQ70m6MLgxIs7KfvdNAVYEngUOiIhF7T5eUZKHmZl1nMLctjIzs47j5GFmZq3m5GFmZq3m5GFmZq3m5GFmZq3m5GEdTlJIOq9k+cfZwIQdse/JkvZsect2H2cvSS9LeqBB+YDS0Xqb+fyplYvuc8cZJ+nHLWzzXUnrlSyfJWm7Djj28Oxc71pSdpuk4a3cx23tjcU6n5OHVcIiYJSklfMOpFRJL9tyHA4cGRFbt/FwrU4eknq28Vgt+S7w7+QREWdExL0dtO864LQO2pd1IU4eVgmfkuZN/mHDFQ2vHCR9mP0cLukhSX+UNFfSOZL2z+YneFHS2iW72U7SU5JelbRL9vmekn4laUY2+N1RJft9WNJU4KVG4hmd7X+mpF9mZWcAWwBXSPpVU5WUdIikmyXdqTRHxrlZ+TnAkkrzRvwuKzsgq8tzki6rTxSSPpR0nqTngaFKc86cm8X0pKR1su0GSLo/q9t9ktZsJJ4js/o/L+kPkpaStDmwG/Cr7Nhrl54DSdtKejY73iRJS2Tlr0s6U9Iz2bqmxkh6HnhfUm0j8TS17xGSXpH0DDCqZPt+2XZPZp8bmZWvX/LdvSBpUFPnxDpRRPjlV4e+SHNpLEvqybsc8GNgXLZuMrBn6bbZz+HAe8DqwBKk8XfOzNadAFxQ8vk7SX/4DCL95dsXGAOcnm2zBPAUMDDb70JgYCNxfhn4K7AKadSB+4HvZuseJM390vAzA4CZ2ftDgLlZHfuSem2vUVqv7P3XgFuB3tnyBOCg7H2QenrXb/s6cFr2/iDgtuz9rcDB2fvDgFuy9+OAH2fvVyrZz8+B45r4zieThqvoC7wBfDUrv5o0kGB9HPWfPwb4n0a+i+HAbcCWwENZ2W1ZeaP7LikfBIjU472+jr8g9X4GWB54FegH/BbYPyvvAyyZ979xv8JXHlYZkUYyvRo4vhUfmxFpPoZFpKGk787KXyT90q53Y0R8FhGvkX55r0saq+wgpeGonwBWIv2CAngyIv7cyPG+DTwYEe9EGrL6d6RfhK1xX0S8HxGfkK5svtLINtsCGwMzsvi2BdbK1i0mDeJX6vqSn0Oz90OB67L315CujBraILvKehHYH1i/hdgHA3+OiFez5av4fP3rB1R8ms9//58Taf6U+qFBWtr3uln5a5GywbUln9keGJt9Rw+SEs2awGPAqZJOBr4SER+3UC/rBK25B2zWWhcAzwBXlpR9Sna7VFIP0l+S9UrH2/msZPkzPv9vteGYOkH6K/a4iLirdEXWeLuwLcGXqTTmxTT+f0rAVRFxSiPrPomIxQ3Koon3LZlMunJ6XtIhpCuA9qivW1P1KjUeOJ10fttKwB4RMbtB+cuSngB2BqZJOioi7m/HcawD+MrDKibSoHs3khqf671O+isc0r343m3Y9V6SemTtIGsBs4G7gO8pDceNpK8qjZ7cnCeBrSStnLVBjAYeakM8jflXfSzAfcCeklbNYltRUmNXKPX2Kfn5WPb+UWDf7P3+pMH/GloGeCs77v4l5R9k6xqaDQyob1cBDqSN9Y+Iu4EVgA1b2PcrWXl9G9bokt3cBRwnSQCSNsp+rgXMjYgLSSPCbojlzsnDKu080ui29S4n/cJ+nnQrpi1XBX8l/eK/Azg6u2X0P6TbRs8oPUp7GS38tRxpVNWxpCGrnweejoiOGqZ8IvCCpN9FxEukv8rvlvQCaYa31Zv57ArZdifwn4cOjgMOzcoPzNY19FPSLbtHSL+k600BTsoaof/94EH2vR0K/D671fUZcGnrq/pv44E1mtt3Vj4GuD1rMC+dU+S/SH9MvCBpVrYMaf7xmdntrA1It0MtZx5V16yKSHqd1FD/bt6xmDXHVx5mZtZqvvIwM7NW85WHmZm1mpOHmZm1mpOHmZm1mpOHmZm1mpOHmZm12v8HCxybHHNSs7gAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "N = arange(6, 32, 2)\n", "err = zeros((len(N),))\n", "xx = linspace(xa, xb, 1000)\n", "for i,n in enumerate(N):\n", " xhat, θ, x, fcheb, S = cheby_dct_setup(xa, xb, n)\n", " y = fun(x)\n", " a = fcheb(y)\n", " yfun = fun(xx)\n", " yche = cheby_sum(xx, a, xa, xb)\n", " err[i] = max(abs(yche - yfun)/abs(yfun))\n", "plt.semilogy(N, err, 'r')\n", "plt.xlabel('Number of Interpolation Nodes')\n", "plt.ylabel('Maximum Relative Error')\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "arctic-orientation", "metadata": {}, "source": [ "After about 20 points, the relative error of the interpolating expansion is below about $10^{-10}$ everywhere along the curve. The classic \"super-exponential\" convergence or \"infinite order\" convergence is also on display. The error is bending downward as $N$ increases, even with the logarithmic y-axis. The error would decrease faster than any finite inverse power of $N$ for large $N$, but after about 25 nodes are used the error reaches machine precision." ] } ], "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.1" } }, "nbformat": 4, "nbformat_minor": 5 }