{ "cells": [ { "cell_type": "markdown", "id": "veterinary-massage", "metadata": {}, "source": [ "# Chebyshev Boundary Value Problem Example" ] }, { "cell_type": "markdown", "id": "average-signature", "metadata": {}, "source": [ "This notebook shows how to use chebyshev expansions to solve a second order boundary value problem (BVP). I define a simple test function that represents the second derivative of another function that we want to find. Because the test function is simple, we can analytically integrate it twice to get the exact solution and compute the error of the chebyshev approximation." ] }, { "cell_type": "code", "execution_count": 17, "id": "powered-librarian", "metadata": {}, "outputs": [], "source": [ "from orthopoly.chebyshev import *\n", "from numpy import *\n", "from scipy.linalg import solve\n", "import matplotlib.pyplot as plt\n", "%matplotlib inline" ] }, { "cell_type": "markdown", "id": "opening-daisy", "metadata": {}, "source": [ "Define the functions, where $S = d^2F/dx^2$" ] }, { "cell_type": "code", "execution_count": 66, "id": "short-broad", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEKCAYAAAAVaT4rAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAnIklEQVR4nO3deXxU5d3+8c83CSFEIAlbkATCJmGTVfZNUYQq4lKtSuGxVMW6UKlV2+pTm6d9uv20VutOtcUW1KqAuACKApZFoiD7vkMCSYCwB7Levz8y+KQIZIDMnJnM9X698sqcM8mc6xi5cnLPOfcx5xwiIhI5orwOICIiwaXiFxGJMCp+EZEIo+IXEYkwKn4RkQij4hcRiTABK34z+5uZ5ZnZ6grr6pnZbDPb5PucFKjti4jI6QXyiH8iMOyUdT8HPnPOXQJ85lsWEZEgskBewGVmzYEPnXMdfcsbgMudc3vM7GJgnnMuPWABRETkW4I9xp/snNvje5wDJPvzTRMnTnRAWH3s3bvX8wzaZ+2z9jni9/m0Ys70RKA555yZnTGYmY0FxgIMGTKEq666KmjZqkJ+fj6FhYVexwgq7XNk0D6Hj9TU1NOuD3bx55rZxRWGevLO9IXOuQnABICMjAx3ph0IZeGY+UJpnyOD9jm8BXuo533gDt/jO4DpQd6+iEjEC+TpnG8CXwDpZpZlZncCfwCGmNkm4CrfsoiIBFHAhnqcc7ef4akrA7VNERGpnK7cFRGJMCp+EZEIo+IXEYkwKn4RkRB0+EQxby/ZxYni0ip/bRW/iEgIytyazxPTV2NW9a+t4hcRCUGLt+6ne1oSNWOiq/y1VfwiIiFo8db99G5RPyCvreIXEQkxBwuKWLvnMH1aqfhFRCLCl9vyqRkTRafUxIC8vopfRCTELN6az2Vp9YiNCUxFq/hFRELM4q376d2yXsBeX8UvIhJCDhYUsS7nML1bBmZ8H1T8IiIhJXNbPnEx0QEb3wcVv4hISFm8dT+XNU8K2Pg+qPhFRELKws37AjrMAyp+EZGQkXv4BBtzjzLgkgYB3Y6KX0QkRCzcvI/E+Bp0aJIQ0O2o+EVEQsSCTfvo16oB0VEBmJmtAhW/iEgIcM6xYPM++rUO7DAPqPhFRELCxtyj5B0pDPj4Pqj4RURCwvxNe0mrH0/TevEB35aKX0QkBCwM0jAPqPhFRDxXVFJG5rZ8Bqj4RUQiw9c7D3CiuJS+rVT8IiIRYcGmfVyamkhCfI2gbE/FLyLisQWb99G/dWCnaahIxS8i4qFDBcWszDpI/9YNg7ZNFb+IiIcWbdlHzZhouqUlBm2bKn4REQ/N3ZBHv9b1qRkTHbRtqvhFRDxSVuaYu2EvV7RtFNTtqvhFRDyyZvdh9h4p5Ip0Fb+ISESYsz6Pto3r0CSxVlC360nxm9lPzGyNma02szfNLM6LHCIiXpqzIY/BQR7mAQ+K38xSgB8DlznnOgLRwG3BziEi4qV9RwtZmXUwMorfJwaoZWYxQDyw26McIiKe+HzDXhJq1aBrs6Sgbzvoxe+cywaeAnYCe4BDzrlPgp1DRMRLczbkMahNw4Dfbet0YoK9QTNLAq4HWgAHgXfMbJRzbtIpXzcWGAswZMgQsrKygh31guTn53sdIei0z5FB+3zhSkodn6/P5aFBqQHtttTU1NOuD3rxA1cB25xzewHMbCrQF/iP4nfOTQAmAGRkZLgz7UAoC8fMF0r7HBm0zxdm8db9FBSXcVOfdBLjY6vsdf3lxRj/TqC3mcWbmQFXAus8yCEi4om56/Po1izJk9IHb8b4M4F3ga+BVb4ME4KdQ0TEK3PW5wX9at2KvBjqwTn3K+BXXmxbRMRL2/cdY1PeUa5s513x68pdEZEg+nhNDs3rx5OeXMezDCp+EZEgmrUmh6EdG1P+Fqc3VPwiIkGSc+gEy3YeZFiHxp7mUPGLiATJJ2tzaFw3js6piZ7mUPGLiATJrNU5DO2QTJQHV+tWpOIXEQmC/GNFZG7LZ2hHb4d5QMUvIhIUn67LpW5cDD2b1/M6iopfRCQYPl6dw5D2ycREe1+73icQEanmjhaWMH/TPoaFwDAPqPhFRAJu7vo8YmOi6NuqgddRABW/iEjAzVy9h8vTGxJXI9rrKICKX0QkoI4WlvDZujyGd2ridZRvVDpJm+9G6MOBAUAT4DiwGvjIObcmsPFERMLbp2tziY2O4vL0hl5H+cZZi9/M/ofy0p8HZAJ5QBzQBviD75fCT51zKwOcU0QkLL2/YjdXd2gcMsM8UPkR/5e+KZRP52kzawQ0q+JMIiLVwsGCIv69cS+v/aCH11H+w1mL3zn30anrfPfMPejK5VH+V4CIiJxi5uoc6taqQb9W9b2O8h/O+uaumT1hZm19j2ua2VxgC5BrZlcFI6CISLj6YMVurrm0cUhctFVRZWluBTb4Ht/h+9wQGAT8LlChRETCXe7hE3yxdT8jOqd4HeVbKiv+Iuec8z0eCrzlnCt1zq3Do9s2ioiEg/eWZdMkoRaXpSV5HeVbKiv+QjPraGYNgSuATyo8Fx+4WCIi4cs5x5Svs/hutxTPp2A+ncqKfzzwLrAe+LNzbhuAmV0DLAtsNBGR8LRm92E25h7lxm6pXkc5rcrO6lkMtD3N+hnAjECFEhEJZ+8uzaJ7WhItGlzkdZTTquysnlF2ljsCm1krM+tf9bFERMJTcWkZ76/YzU3dQu9N3ZMqe4O2PrDczJYCS4G9lF+525ryM3v2AT8PaEIRkTAyb8NejhaWMPzS0Jmb51SVDfU8a2bPA4OBfkAnyufqWQeMds7tDHxEEZHw8c6SXQxpn0xCfA2vo5xRpadkOudKgdm+DxEROYO8wyf4bH0eE8eE1hQNpwqty8lERMLYO0uzaJIYR78QueHKmaj4RUSqQFmZ462vdnLrZU1D8tz9ilT8IiJVYNGW/ew+eIJbLmvqdZRK+TXtgpnVBL4LNK/4Pc65XwcmlohIeHnzq51ckd6I5LpxXkeplL9H/NOB64ES4FiFDxGRiLf3SCGfrMnh9p6hf7QP/k+0luqcGxbQJCIiYeqtL3fSOCGOy9MbeR3FL/4e8S8ys0sDmkREJAwVl5YxOXMno3unER3ib+qe5G/x9weWmtkGM1tpZqvM7Lzvs2tmiWb2rpmtN7N1ZtbnfF9LRMRLs9fmcvB4Ed8Lgzd1T/J3qOc7VbzdZ4FZzrmbzSwWTfEsImFq4qLt3NAlhcT4WK+j+M2vI37n3A4gEbjO95HoW3fOzCwBGAi85nvtIufcwfN5LRERL63bc5gvt+Uzuk+a11HOiV/Fb2YPApOBRr6PSWY27jy32YLyyd7+bmbLzOxVMwvNuUtFRM7i7wu30aN5Eh2aJHgd5Zz4O9RzJ9DLOXcMwMz+CHwBPHee2+wGjHPOZZrZs5TP8PnLil9kZmOBsQBDhgwhKyvrPDblnfz8fK8jBJ32OTJon8vtO1bMtK+z+fWwtJDtp9TU098Ixt/iN6C0wnKpb935yAKynHOZvuV3Oc3Uzs65CcAEgIyMDHemHQhl4Zj5QmmfI4P2Gd78eD2p9eK5tX/7kJ+i4VT+Fv/fgUwzm+ZbvgHfGP25cs7lmNkuM0t3zm0ArgTWns9riYh44VhhCZMW7+TRYelhV/rgZ/E75542s3mUn9YJMMY5dyH33B0HTPad0bMVGHMBryUiElRvL9lFTJTx3RC9p25lzlr8ZlbXOXfYzOoB230fJ5+r55w7r8E+59xy4LLz+V4RES8Vl5bx2oJtjO6TRlyNaK/jnJfKjvjfAIZTfttFV2G9+ZZbBiiXiEhIem9ZNgeOFXFHn+ZeRzlvld16cbjvc4vgxBERCV2lZY4X521hVJ80ki4Knwu2TuXvefyf+bNORKQ6+3DlbvYcOs7dA8J7sKOyMf44yqdTaGBmSfzfKZx1gZQAZxMRCRllZY7n52xmZM80GtSu6XWcC1LZGP89wHigCeXj/CeL/zDwfOBiiYiElhmr97Ajv4BJg8L7aB8qH+N/FnjWzMY5587nKl0RkbBXUuZ4evZGRvVKC4s7bFXG32mZy8ws8eSCmSWZ2X2BiSQiElpmrc8n99AJ7r+ilddRqoS/xX93xRk0nXMHgLsDkkhEJIScKC7lb1/mcmf/FtQP87H9k/wt/mgz++a6ZDOLBsL3XCYRET9NztzJieIy7hoY/mP7J/k7V88s4F9m9opv+R7fOhGRauvQ8WKen7OJUd0bUTeuhtdxqoy/xf8zysv+Xt/ybODVgCQSEQkRz322iTpxNbi5cwOvo1QpfydpKwNe8n2IiFR72/Yd4/UvtvPc7V2JjS6t/BvCiL9X7vYzs9lmttHMtprZNjPbGuhwIiJe+d2MdXRrlsTQDo29jlLl/B3qeQ34CeUXcVWvX30iIqeYv2kvn67L5YMH+lPhvJZqw9/iP+ScmxnQJCIiIeBEcSm/fG81o3ql0TElvO6l6y9/i3+umT0JTAUKT650zn0dkFQiIh55ad4WjhaW8vDQdK+jBIy/xd/L97nizVMcMLhq44iIeGfbvmO8NG8LT97SiYRa1ef0zVP5e1bPFYEOIiLipbIyx2NTV9GjRRIjOjfxOk5A+VX8ZvbE6dY7535dtXFERLwxOXMHy3cd5OPxA6vlG7oV+TvUc6zC4zjKb8e4rurjiIgE3879Bfx+5np+/p22NKsf73WcgPN3qOdPFZfN7Cng44AkEhEJorIyx6NTVtApNYHRvdO8jhMU/h7xnyoeSK3KICIiXvjr/K2syjrEzAcHEhVVvYd4TvJ3jH8V5WfxAEQDDQGN74tIWFux6yBPfryB/3dzp4gY4jnJ3yP+4RUelwC5zrmSAOQREQmKIyeK+fFby7iucxNu6hZZAxhnnavHzD4BcM7tAEY653Y457JV+iISzpxz/GLqKgB+fX0Hj9MEX2WTtDWs8PiWQAYREQmW1xZsY/baXF4Y2Y061WiefX9VNtTjKnleRCSsLNqyj9/PXM9Tt3SqtnPxVKay4m9pZu8DVuHxN5xzIwKWTESkiu3cX8C4N5YxuncaN3aNrHH9iior/usrPH4qkEFERALpUEExYyZ+SbuL6/L4te28juOpsxa/c+7zYAUREQmUwpJS7pm0hJioKF4c1Y0a0X7dg6raOt8LuEREwkJpmeOht1ewZe8x3ru/X7W6afr5iuxfeyJSrZWVOX4xdSXzN+7lHz/sSUpiLa8jhQTPit/Mos1smZl96FUGEam+nHP85qO1fLhyD6//sCftLq7rdaSQ4e+UDW2AR4C0it/jnLuQG7E8SPkMn/ppiEiVKitzZHywhreX7OJvP+hB12ZJXkcKKf6O8b8DvAz8lSq42bqZpQLXAr8FHrrQ1xMROam0zPH4tFV8sGI3r4/pSa+W9b2OFHL8Lf4S59xLVbjdZ4BHgTpV+JoiEuGOF5Xyk38tZ+GWffzzrl5005H+aflb/B+Y2X3ANP7zZuv557pBMxsO5DnnlprZ5Wf5urHAWIAhQ4aQlZV1rpvyVH7+Of+nCXva58gQqvucX1DMzz7aTn5BMS/c2JJGUcfIyjpW+Tf689ohus+VSU09/UVq/hb/Hb7Pj1RY54CW55GlHzDCzK6h/G5edc1sknNuVMUvcs5NACYAZGRkuDPtQCgLx8wXSvscGUJtnzfnHeG+yV+RFB/LB+P60KhuXJVvI9T2+UL4eweuFlW1QefcL4BfAPiO+B8+tfRFRPz18ZocHn5nBb1a1Ocvt3chPlaXJ1XG37N6agD3AgN9q+YBrzjnigOUS0TkrIpKyvjjrPX8feE2HriiNQ9e1YboCLmD1oXy91fjS0AN4EXf8mjfursuZOPOuXmU/xIREfFb1oECHnhjGbvyC5g4picD2zSs/JvkG/4Wfw/nXOcKy3PMbEUgAomInIlzjje/3MXvZqyjfZO6zHhwAMkBGM+v7vwt/lIza+Wc2wJgZi2pgvP5RUT8tSu/gJ9NWcnXOw/w8NXpjOnXQkM758nf4n8EmGtmWymfmz8NGBOwVCIiPoUlpUxcuJ1nP9tEx5QEZj04kOYNLvI6Vljz96yez8zsEiDdt2qDc67wbN8jInIhnHN8ti6P//1oLYeOF/P4te24vUczonSUf8HOWvxmNtg5N8fMbjrlqdZmhnNuagCziUiEWp19iD/OWs+iLfsZ3TuNn1zVhoR4TadcVSo74h8EzAGuO81zDlDxi0iV2ZBzhD/P3sisNTkMbtuImQ8OoE2yZnapapXdgetXvoe/ds5tq/icmVXZRV0iEtlWZR3i5X9vYcaqPfRr1YCp9/XVPDsB5O+bu1OAbqesexfoXrVxRCRSOOf496Z9vPL5FhZt2c+ASxrwxl296dNKs2kGWmVj/G2BDkDCKeP8dSmfZ0dE5JwcLyrlg5W7+duCbWzKO8p1nS5mxo8H0L6Jbs0RLJUd8acDw4FE/nOc/whwd4AyiUg1tCHnCG9k7mDqsmwAbunelFfvuIzUpHiPk0Weysb4pwPTzayPc+6LIGUSkWriRHEpH63cwxtf7mTpjgN0bZbIE8PbM7xTE2rFRnsdL2L5O8Y/1sy+dYTvnPthFecRkWpgc94R3sjcxZSvsygtc9zYNYX/vaGj7nsbIvwt/oo3RI8DbgR2V30cEQlXhSWlzFqdw+TMnXy5LZ/OqQk8dk1bruvcRFMlhxh/r9ydUnHZzN4EFgQkkYiEla17j/Lmlzt5d2kWRSVljOiSwhPD29MxJcHraHIG5/tr+BKgUVUGEZHwUVxaxqdrc5mUuYOFm/fToUldHh6azvVdUqhdU0f3oc7fG7EcofxKXfN9zgF+FsBcIhKC9hw6zmuZOXy0fj2HjhdzXecmvHd/PzqnJmCmOXTChb9DPbpmWiRClZU5Fm7Zx6TFO/h0XR5N6sYydmBLbu6eSmJ8rNfx5DxUdgHXqVfr/gfn3NdVG0dEQsWJ4lKmL8/mr/O3sW3fMa5un8w/ftiTZjWP07RpU6/jyQWo7Ij/T2d5zgGDqzCLiISAA8eKmJy5g4mLdlBYXMrI3s34Qd/mXJxQC4CsrCyPE8qFquwCriuCFUREvLUrv4BX52/l7SVZ1Lsolh8NasmtPZpSJ07TIVc3/r65WwO4FxjoWzUPeMU5VxygXFXi0PFinHPUio0mNjpKbz6JnMau/AJenLeZd5Zkkd64Dn+8uRPXdGxMTHSU19EkQPw97+oloAbwom95tG/dXYEIVVWe/Hg9kxbvBCAmyoiPjSa5bhwpSbVISaxFWv14OjZJoENKAgm1dFQjkaVi4be7uC6vjO7O4LaNdIAUAfwt/h7Ouc4VlueY2YpABKpKD17ZhpE90zheXEJBUSlHT5SQe/gEWQeOk33wOEt3HOCPszZQWuZo0eAiujVLYlB6Qwa0bkDSRTpbQaqn/UcL+ctnm5icuVOFH6H8Lf5SM2vlnNsCYGYtgdLAxaoaDevUpGGdmmf9mhPFpazdc5iVuw6SuS2fx6eu4mhRCZ1TE7nm0saM6JxC4wTNQC3h73hRKX9buI2X5m0huW5NXvh+N65un6zCj0D+Fv8jwFwz20r5RVxpwJiApQqiuBrRdGuWRLdmSfygXwuKS8tYvusgc9bn8fqiHfx+5nr6tqrPTV1TubbTxcTV0IyCEl7KyhxTvs7iT59spKSsjJ9/py239WiqMfwI5u8FXJ+Z2SWUz88PsME5Vxi4WN6pER1Fj+b16NG8Ho9cnc6SHQeYtiybjPfX8LsZ6xjZqxmje6fRqK7+CpDQtyrrEP89fTUbcg4zdkBLxg5qpSkVpNILuHoAu5xzOc65QjPrAnwX2GFmGc65/GCE9EpUlNGzRT16tqjH49e2Y8rSLCYu2s7Ln2/h+i4pjBvcmrT6F3kdU+RbDhYU8dQnG5icuZNhHRrz0ve70SSxltexJERU9qv/FeAqADMbCPwBGAd0ASYANwcyXCipXTOGO/o2Z3TvNOZtzOO5OZsZ/KfPublbKg8Mbk3TerqLkHjPOcc7S7P4w8z1JNSqwcQxPRnUpqHXsSTEVFb80RWO6m8FJvimaJ5iZssDmixERUUZg9smc0V6I+Zt3MszszdyxVPzGNU7jfFXXaK5S8QzO/cX8POpK1m64wDjBrfm7oEtqRmj96Tk2yotfjOLcc6VAFcCY8/he6s1M+OK9EZc3qYhn6zN5fcz1vHe8mweGtKGkT2beR1PIkhpmeP1Rdt58uMNdEpN4OPxA2neQEOQcmaVlfebwOdmtg84DswHMLPWwKEAZwsLZsbQDo25PL1h+T++WRuYtHgH4/s3JjXV63RS3W3OO8rPpqxk/Z7DPH5tO0b2bEZUlE7PlLM76/lczrnfAj8FJgL9nXOuwveNC2y08FIzJpqxA1sx95HL6ZiSwP1TNvP4tFUcPhHSs1pImHLOMXHhNq75y3xq14zhk4cGMap3mkpf/FLpcI1zbvFp1m083w2aWVPgH0Ay5TN8TnDOPXu+rxdqGtSuydPf68KApjX58/wcZq/9nN/c0JGhHRp7HU2qibwjJ3jknZV8uS2fjOs6cHvPproIS86JF1dwlAA/dc61B3oD95tZew9yBFSPpnX4ePxAbuyWwr2TlvLIOys4oqN/uUCfrs1l2DPzyT9WxEc/7s/IXs1U+nLOgl78zrk9J2/g4pw7AqwDUoKdIxhqxUbzi++04+17+pC5LZ/vPDufr7ZX60sfJEAKikp4bNoqxv5zCbf1aMqUe/vSsmFtr2NJmPL0mm0zaw50BTK9zBFolzWvx4wHB9C/dQNufeUL/jhrPUUlZV7HkjCxbs9hhj+3gHnr83jj7t48OqwtsTGabkHOn2enZJpZbWAKMN45d/g0z4/Fd/rokCFDwu6uP/n53z6yf6BXPbo0iuYPc3Ywd+1ufj0sjSZ1zz6JXDg53T5Xd4He5w/X7ufpz7Pp27wuPxvcgjqxxz3/t6Cfc/hIPcOphZ4Uv+/GLlOAyc65qaf7GufcBMqvDiYjI8OdaQdC2eky35aayuAurXjwzeXc+fZmnrqlc7V64zccf04XKhD7fLyolF9OX8305dk8dk07ftC3eUiN5evnHN6C/veilf/f+xqwzjn3dLC3Hwoa1Ylj0l29GNOvBfdOWspvPlyroR/5xua8o9zwwkK+2LKft+/pw5h+LUKq9CX8eTFQ2I/yO3gNNrPlvo9rPMjhqego46EhbZg4pifvLcvme698QfbB417HEo9NX57NiOcX0CQxjg/H9adrsySvI0k15MVZPQucc+ac6+Sc6+L7mBHsHKFiYJuGfPTjAdSINkY8t4DFW/d7HUk8UFxaRsb7a3jo7RU8MLg1r93RQ3eBk4DRqQEhoHFCHG/c3ZvhnS7m+69mMnHhNv7vImmp7vYdLWTUq5lMX57NP3/Yk/sub60rcCWgInqitVBSIzqK/7m+Ix1SEvjvaatZlX2Y397YUXf8quZWZR3inn8uIemiWD4Y15/UJE3vLYGnI/4Q873LmvKve3qzYPNevvfKF+zWuH+1NWVpFt99eRE9W9Tj3R/1VelL0Kj4Q1DXZkl8MK4/sdFRjHh+AV9uC89ziOX0To7nPzplJY8OTefPt3ahVqz+spPgUfGHqEZ1ysf9h3ZozMi/LuaNzJ1eR5IqcOp4/l0DWupUTQk6jfGHsNiYKH5746W0u7guT0xfzbo9h3niuvbUiNbv63Ck8XwJFWqQMDCqdxqT7+rFR6v2MPq1TPKPFXkdSc7RtGVZ3KzxfAkRKv4w0atlfabf34+DBcWMeH4B63O+Nb2RhKCS0jJ++9Fafvr2Ch6+WuP5EhpU/GGkab14pt7Xl0tTErjpxUXMWp3jdSQ5i4MFRYyZ+BVvL8li4pie3D1Q4/kSGlT8YSY+NoYXRnbjR4Nacd/kpTzz6UbKynSxV6jZmHuE619YSO7hE7z/QD8GtmnodSSRb+jN3TAUFWX8+MpLaJNch4feXs6GnCM8dUtnLqqpH2co+HhNDg/9azn9Wjfg6Vu7UFs/FwkxOuIPY8M6NmbqfX1ZvfsQ331pEbvyC7yOFNHKyhzPfLqReyct5e6BLXl5VHeVvoQkFX+Ya9u4LtPv709SfCzXv7BQk7x55GhhCT+atJS//nsrL43qzvir2mi+HQlZKv5qoN5Fsfzjzp5c1+liRr2ayaTFO7yOFFF27D/GTS8uZEPuEabd369a3VhHqif9HVpNnJzkrW2Fi71+dV0H3Zs1wL7aeYSM2WvplJrAc7d3JTFeUylL6FPxVzO392xG60a1+dE/l7Ip7ygvfb8b9WtXn/v6hoqyMseL8zbz9Oyt3DWgJY8OTSdGV1RLmND/qdVQj+b1eH9cf44VljDi+YWs3a2LvarSoYJi7vrHEl6at4WMoWk8dk07lb6EFf3fWk2lJNbi3R/1pUuzRL770iJmrtrjdaRqYXX2IYY/P59d+QVMf6A/g1sneh1J5Jyp+KuxWrHRPH97Vx4Y3Jr73/iap2frYq8L8a+vdnLTS4vo2jSJ9+7vR+tGtb2OJHJeNMZfzZkZ91/RmjbJdRj/1jI25BzmqVs6UyeuhtfRwkZBUQkZ769h6tfZ/HJ4e/6rT5qmXpCwpiP+CDGkfTLT7u/HxtyjXPuXBazMOuh1pLCwdvdhrntuAQs27eNf9/Thjr7NVfoS9lT8EaRNch0+GNefbr5x/9cW6KbuZ+KcY+LCbdzw4kLaJNdh5oMD6Z6W5HUskSqhoZ4IU7tmDH++tQt9Wzfgiemr+WLLfp68uRNJF+n885PyjxXx6LsrWLB5H/8zogO39Wiqo3ypVnTEH4HMjO9d1pQPHujPrvwCrvnLfP69ca/XsULC7LW5DH3m32QdOM4HD/Tn9p7NVPpS7aj4I9glyXV4zzfFwB1//5LHpq3iaGGJ17E8cbCgiPFvLeOefy7hpm4pvHd/Py5JruN1LJGA0FBPhKsVG03GiA5c3SGZR99dyb837uXJmzvTp1V9r6MFzSdrcnhs2moS42sw5d6+dG2msXyp3nTELwD0bdWAWeMHMuCShox8dTEZ76/hyIlir2MFVO7hE4x7cxk/mrSUm7un8uG4/ip9iQg64pdv1K4Zw+9vupRhHRvz+LRVzFi1h18Ob8/wThdXq3Hu4tIyXl+0nWc+3URa/Xim3tePLk0TvY4lEjQqfvmWQW0aMvsng3hh7mYeens5//pqFxkj2tO6UfiPeS/YtI/ffLiW3YeO8+iwdL7fK41ozZsvEUZDPXJatWKjeXhoOrPGD8QMhj4zn19MXUXu4RNeRzsvq7MPMfq1TP7rb5l0bprA3Icv57/6NFfpS0TSEb+cVauGtfnHD3vy+ca9/GHmegY9OZc7+7fg7gEtw2Lu+Y25R3h+zmY+WLmbwemNmPngQNIbh/9fLiIXQsUvlTIzLk9vxMBLGvLe8myenr2Rvy/czu09m3Fn/xY0SazldcRv+XrnAV6cu4VP1+XSs3k93rq7N71aRs6ZSiJn40nxm9kw4FkgGnjVOfcHL3LIuYmKMm7qlsp1nZvw4crdvDxvK68v2s6Izk0Y2asZydHeTv9woriUj1bu4c0vd7JkxwGuateIKff2oXtaPU9ziYSaoBe/mUUDLwBDgCzgKzN73zm3NthZ5PzUiI7ixq6p3NAlhbkb8vj7wu3c8soXNE2oyff7FHFD1xSS68YFJUtZmWPZroO8vzybacuycQ6u79qE3954qYZ0RM7AiyP+nsBm59xWADN7C7geUPGHGTNjcNtkBrdNJutAAa/NWcs/vtjB72eup1NqAkPaJXNlu2TaNq5DVBW+iXq0sIQl2/OZvTaX2WtzyTtSSPe0JP7bd+ppfKxGMEXOxot/ISnArgrLWUAvD3JIFUpNiufOXo355Y3dWZl9iE/X5vLRqj38afZG6sbF0LlpIl2bJXFpSgJp9eNJTapVaUE759h3tIgd+4+xbd8xVmUfYsn2A6zPOUx0lNG7ZX1+fOUlDGmfHLS/MESqg5A9NDKzscBYgCFDhpCVleVxonOTn5/vdYSgO7nPDQxu63ARt3VoSc6RItbkHGNNTgGfrc7m5XmbKSotfy8gqVYMibViqBljxMVEERttnChxFBSXcqyojEPHSygoLgOg4UU1aNUgjn5pF3Fv74a0axRPXI3ys5GLD+8jy6PbCkfyzzmShOs+p6amnna9F8WfDTStsJzqW/cfnHMTgAkAGRkZ7kw7EMrCMfOFOnWfU4HL2v3fclmZY9/RQnYdOE7WgQIOFhRzvLiUE8WlnCguo1aNaGrHxVC7ZjQJtWJp3iCeZvXiQ3r4Rj/nyFCd9tmLf01fAZeYWQvKC/82YKQHOcQDUVFGo7pxNKobpxubiHgk6MXvnCsxsweAjyk/nfNvzrk1wc4hIhKpPPn72Tk3A5jhxbZFRCKd5uoREYkwKn4RkQij4hcRiTAqfhGRCKPiFxGJMCp+EZEIY855O5WuP8zsVcrn9Akn3YGlXocIMu1zZNA+h4/tzrmJp64Mi+IPR2a2xDl3mdc5gkn7HBm0z+FPQz0iIhFGxS8iEmFU/IEzwesAHtA+Rwbtc5jTGL+ISITREb+ISIRR8QeBmf3UzJyZNfA6S6CZ2ZNmtt7MVprZNDNL9DpToJjZMDPbYGabzeznXucJNDNramZzzWytma0xswe9zhQMZhZtZsvM7EOvs1QVFX+AmVlT4Gpgp9dZgmQ20NE51wnYCPzC4zwBYWbRwAvAd4D2wO1m1t7bVAFXAvzUOdce6A3cHwH7DPAgsM7rEFVJxR94fwYeBSLizRTn3CfOuRLf4mLK775YHfUENjvntjrnioC3gOs9zhRQzrk9zrmvfY+PUF6GKd6mCiwzSwWuBV71OktVUvEHkJldD2Q751Z4ncUjPwRmeh0iQFKAXRWWs6jmJViRmTUHugKZHkcJtGcoP3Ar8zhHlQrdO1iHCTP7FGh8mqceBx6jfJinWjnbPjvnpvu+5nHKhwYmBzObBJ6Z1QamAOOdc4e9zhMoZjYcyHPOLTWzyz2OU6VU/BfIOXfV6dab2aVAC2CFmUH5kMfXZtbTOZcTxIhV7kz7fJKZ/QAYDlzpqu/5wtlA0wrLqb511ZqZ1aC89Cc756Z6nSfA+gEjzOwaIA6oa2aTnHOjPM51wXQef5CY2XbgMufcPq+zBJKZDQOeBgY55/Z6nSdQzCyG8jevr6S88L8CRjrn1ngaLICs/AjmdSDfOTfe4zhB5Tvif9g5N9zjKFVCY/xS1Z4H6gCzzWy5mb3sdaBA8L2B/QDwMeVvcr5dnUvfpx8wGhjs+9ku9x0NS5jREb+ISITREb+ISIRR8YuIRBgVv4hIhFHxi4hEGBW/iEiEUfGLiEQYFb+ISIRR8YucB9+89EN8j//XzJ7zOpOIvzRXj8j5+RXwazNrRPkslSM8ziPiN125K3KezOxzoDZwuW9+epGwoKEekfPgm331YqBIpS/hRsUvco7M7GLK7zNwPXDUNyOpSNhQ8YucAzOLB6ZSfu/ZdcBvKB/vFwkbGuMXEYkwOuIXEYkwKn4RkQij4hcRiTAqfhGRCKPiFxGJMCp+EZEIo+IXEYkwKn4RkQjz/wGu7pAuQ7vqlgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "#the solution function\n", "F = lambda x: sin(x) + exp(x/2) - x/3\n", "#the source function (second derivative of F)\n", "S = lambda x: -sin(x) + exp(x/2)/4\n", "#domain to look at\n", "xa = -5\n", "xb = 5\n", "\n", "x = linspace(xa, xb, 250)\n", "plt.plot(x, F(x))\n", "plt.xlabel('$x$')\n", "plt.ylabel('Solution Function (S)')\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 81, "id": "orange-sympathy", "metadata": {}, "outputs": [], "source": [ "#number of points to use\n", "n = 16\n", "#setup elements for solving BVP\n", "xhat, θ, x, A = cheby_bvp_setup(xa, xb, n, 0, 2, 0)" ] }, { "cell_type": "code", "execution_count": 82, "id": "signal-teacher", "metadata": {}, "outputs": [], "source": [ "#source array\n", "b = zeros((n,))\n", "b[0] = F(xa)\n", "b[1:-1] = S(x[1:-1])\n", "b[-1] = F(xb)" ] }, { "cell_type": "code", "execution_count": 83, "id": "square-receipt", "metadata": {}, "outputs": [], "source": [ "#solve for the coefficients of the solution cheby expansion\n", "a = solve(A, b)\n", "#evaluate the cheby solution and the exact solution\n", "xx = linspace(xa, xb, 10000)\n", "yexact = F(xx)\n", "ycheby = cheby_sum(xx, a, xa, xb)" ] }, { "cell_type": "code", "execution_count": 86, "id": "social-constant", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Maximum error = 6.2541e-08\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEMCAYAAAA8vjqRAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABLy0lEQVR4nO2dd3wcd5n/P8/2lbTqzZZkW+4ljmOnOL0Rg0MSAoQcDiGUBAKhHNxxB+G4A8Md5Do/OEIJJISQUEJCSSA9IT1Oc+LeLcmSrL6SdrW9fH9/7M5qJW2Z2Z3Z2V0979dLL1ur3ZlnZuf7fb5P/ZIQAgzDMAyjBIPeAjAMwzClBysPhmEYRjGsPBiGYRjFsPJgGIZhFMPKg2EYhlEMKw+GYRhGMaw8GIZhGMWw8mAYhmEUY9JbgFwgogsAXI+Y/GuFEOfqLBLDMMy8ouCWBxHdRUTDRLR31utbiegQER0lolszHUMI8YIQ4tMA/gzgF1rKyzAMw8yFCt2ehIguBDAF4B4hxCnx14wADgPYAqAPwOsArgNgBHDbrEPcKIQYjn/ufgA3CSHcBRKfYRiGgQ6WhxDieQDOWS+fBeCoEOK4ECII4DcArhZC7BFCXDnrR1IciwBMylEcd999twBQcj8jIyO6y8DXzNfM1zyvrzctxRLzaAPQm/R7H4DNWT5zE4Cfp/sjEd0M4GYA2LJlCy677LJ8ZSw4TqcTgUBAbzEKCl/z/GC+XXOpXm97e3vavxWL8lCMEOIbWf5+B4A7AGD79u0i000oZkpV7nzga54fzLdrLrfrLZZU3X4AHUm/t8dfywsiuoqI7jh06FC+h2IYhmGSKBbl8TqAFUTUSUQWANsAPJTvQYUQDwshbl61alXeAjIMwzDT6JGq+2sArwBYRUR9RHSTECIM4HMAHgdwAMD9Qoh9KpyLLQ+GYRgN0CPb6johxAIhhFkI0S6EuDP++iNCiJVCiGVCiG+rdC62PBiGYTSgWNxWDMMwTAlR1sqD3VYMwzDaUNbKo5zdVpPeEB7fN4gn9g3C5Q/pLQ7DzDuEENjbP4k/vd2P3X0TKHS3Dr0p2TqP+czvd/bhG3/ah4gQEAIwGQjb37MO15xeXnnkDFOsDLv9+NL9u/DCkVE0VFow5gni/OWN+N6209BQZdVbvIJQ1sqDiK4CcNW2bdv0FkU1HnizD19+YBe+snU1bjy/E0IAP3+pC//4wC74QhF8+OzFeovIMGXNuCeI63/6KqxmA578uwuxosWBYyNT+Ntfv4Xrf/YqHrjlXFRZy3pqBcBuq5Li2MgU/vmPe/C1K9biUxctg9logMVkwKcuWobvvG89tj+0D2+dGNdbTIYpW4QQ+PKDu2Egwn2fOBsrWhwAgGVNVbj3ps3wBiP414f36yxlYSh/9VhG3PbIAZy+uA43nrdkzt+2nbUILx8bw9f+sBcPf/58GA2U83miUYEn9g/iucMjCEcEzlxSj/ecthA2szEP6RmmsOw/6cKfdvVjcNKPzsZKXHtGB9pq7Xkd8y97BvDcoRE89PnzUGM3z/hbXaUF/3HNqfjQz3bgg2d1YNOiurzOVeyUteVRTrzZM46nDw7jq5evAVFqxfDPV6xBz5gHf3gr984uTk8QH7zjFXzhN2/D6QnCH47i248cwDu/+zz29k/mfFw5THiDcHpDiETnV+BxPjEVCOPkhA+hSFSzc0SiAv/25/244v9ewOtdTliMBjyyZwCX/PezuP/13uwHSEM4EsX/PHEYHz13MVa3Vqd8zznLGrBlTQtuf+ZozucpFcra8iinmMcvX+nGRSubcEpbTdr3NFfbcP3Zi/GT547h/RvbYFBofXgCYVx3xw6YjISnv3QR2usqAABufwi3/n4PrvvpDvz6k2dnlEEp/lAEd77Yhd+8fgK9Th8AoMZ+BO9a14JbLl6OzsZK1c7F6EMoEsUfdvbjnh3d2NvvAgAYDcDFK4fw+XeswGkdtaqdSwiBf3xgF57aP4Sff+xMXLyqOfH6Pa/04Nbf7wYRcO0ZHVmONJc/7x7AkMuPWy5envF9t1y8DO/74cs4MuROuLXKkbK2PMol5uH0BPHInkF8eHP2YPiN53Wia9SD54+MKD7PP/9xL7yhMO77xOaE4gAAh82M72/biItXNeNTv3wTkz51UoN7xjy44vsv4OcvdeP6zYvx8OfOxz3XrcI3rlqLYyMevOv/PY+7X+qadymQ5USv04trf/wKvvXn/ThveSN+9+lz8MyXLsJt7+6EyUh43w9fwvefPqLad3zH88fx6J5B3PeJsxOKAwCICB89dwm+fuVafO0Pe3FkSPn+cffu6MH7N7WhvtKS8X0bF9XhlLZqPLgz796uRU1ZK49iY8IbxDf+tBef/dVO7O6bkP25x/cNotpuwsWrmrK+t7XGhktWN+P3Ch/cl4+N4o9v9+P72zaitmLu4DAaCP9xzXpYzQZ8S4WAYK/Ti2t+9AoW1trx9JcuwqcvWob17TVY2mDD+ze144FPn4NvvmcdvvPoQXz7LwdYgZQgR4en8P4fvQyryYCn/v4ifPXyNThzST2WNlXh3CXV+MkNZ+BH15+OHz17DP/+2MG8z3dsZAr/8+RhfOvqdVjfnto6/ui5S3Du8gb80x/2KHqmDg+58UbPOD50lrxsxvdtbMef3u5HVKYLNhCO4HtPHcEnfvEGHts7KFsuPWHlkQfHRqZkxwFCkSg+ctdreOHoKLyBMD74kx3Yf9Il67NP7h/CZWtaYDLK+7qu2dSGx/cNwi2zeFAIgX9/9CCu2dSOjRmCfBUWE/79/afiwZ192NU7IevYqfAGw7jpF69jVWsV7vzomXMCj0BspXjdWYvw84+diXtf7cEPnz2W8/mYwjPk8uNDP92BDe21uOems9BaY0v5vq2ntOInN5yOO1/owoNv9uV1zm89vB+bO+vxgQz1TkSE7Vetw84TE3jhyKjsYz+86yTWLazG2oWpYx2zuerUBRiY9GOXjEWiEAJffmA37nmlGw6bCZ+57008umdA1nmmAmG8cGQEgXBE1vvVpKyVh5btSV45NoZ3fvd5XPl/L8p66H/xcjd6xrz47c3n4K6PnYkLVzbiq7/fnXX14wmE8eLRUWxZ2yJbtktWN8NkIDx/WN7geK3LiT39k/j8pZl9uQBwVmc93rWuBf/9RO739LtPHobLF8btH9oEiynzI3je8kZ8929Ow/88cQjPHBzK+ZxM4QhFovjMfTuxsNaO26/fCKspc5behSub8A/vWoXtD+1D/4Qvp3PuPDGO5w6P4CtbV6dNKJFY0liJ957Whh8oCGo/sW8IW9e1yn5/c7UNaxdUyxqDzx4ewcO7TuLuj5+F737wNHzm4uX4xkP74AtmVgjeYBjX/PBl3HDna7jp7jdkWzlqUdbKQ6uYhxAC33nkAN6/sQ1f3roK//aX/Rm/6HAkirte7MLNFy5Fk8MKIsI/X7EWe0+68NzhzLGJN3rGIYTAucsaZctnNRlx/opG/PXQsKz33/1yNy5b04LFDfKC05+7ZAVeODKKg4PyLKdkDg66cOeLXfjX956S0j2WisvXL8AnLliKrzy4BxPeoOJzMoXlJ88dw7GRKfzow5uyKg6JT16wFEubq/Afj+bmvvrhX49iy9oW2ckcHz9vCV7rduLo8FTW93aPenBoyI13KlAeAHDRqiY8dzj7GPzeU0fwwTM7Eq62Wy5eFksyyJI1ee+OHox5gvjTZ8/Da11OPHNQ3nhXi7JWHlpxcNCNPf2T+Mwly/Gxc5cgKmL53+l48egoRqYC2HbmdIZHR30Ftq5rxf1vZE4dfKPbifVtNbBblNVYXLKqGc8eGs66GnH5Q3j6wPAM2bKxvr0Gmzvr8fMXuxXJBMQGyrnLGhVZUgDw91tWosZuxm2P5O8bZ7Tj+MgUvv/MUfzzFWuxoEZ+TYXRQPiny1fjoV0nFaeED0z68MzBYdx4Xqfsz5zSVoO1C6rxuyzjDwCeOzyC9jo7VrZUKZLrguWNeLt3At5gOO17Dg668HbvBD567pLEa5VWE64+rQ2/ezO9bFL22I3nL8GGjlpcvj77XKI2rDxy4JmDw1jZUoXOxkpUWEy4bE0LntiXPsj19IFhnL20YU7Pm/dubMNTB4YzNjZ8rcuJM5fUK5bxwpVNGJ0K4lCWrJIn9w3BZjbgghXZg/HJfGjzIvxlzwD8Ifm+1sNDbjy6dxBfvGyFonMBgM1sxPar1uH+N3tzsniYwvCfjx3CpkW1uGZTm+LPbl7agAtWNOKO548r+tzv3ujDovoKnL1U2Th5/6Y2PLJ3IKvr+NWuMZy9tCGrO2w2p3bUQgDY05deGf7p7ZM4tb1mTt3Ieze24a0TExh2+1N+7sCAG33jPlx16kIAwLvWteKlo6MIhrWrn5kNK48cePnYKC5MmmwvW9OMF4+OIpyi8EkIgWcODuPS1c1z/nbhykYYifDKsbGU5wlHBN7uncDpi5VXqi6staOt1o43ezK3K3l07yDeua41a+xhNlvWtiAqBJ4+IN9Uvm9HD05fXIczclCGAHD+ikZcuKIJ//kYt9gvRvb0TeKxfYP4soy4QzpuPK8Tj+wZwLAr9aQ5GyEE/vBWP649o0PxOd+xpgW9Th+OZHBdCSHwWpcTZ3Uqf2arrCasbHbg7QzJJc8dGsFla+Za4evbalBbYcZLR1PHTP56aBirWhzoqI+l1J+3vBHeUERWgF4tWHkoJNaG2YVTkwqbTl9SB28wgoODc1f5/RM+9E/4cMGKuTELq8mIszrr0z4gJyb8CISj2JBjEdWmxXXYmUF5hCNRvHp8DBetVGZ1ALHMq3eubcHDu07Ker8/FMEf3upX5B5Lxd9vWYlnDg7jUIp7zejL954+gnesbs6rLcdFK5vQVmfH72RmXh0bmULXqAeXn6IsHgEAnY2VWNZUiacOpE/EODbiwehUEGd3Nig+PgBsXFSLt05MpPzbsMuP/QOulOPPaCCct6wRLx5JvbDc2TOOMzun73ON3YyljZUZrRy1KWvloUW2Vd+4D5O+EE5JStlrdtjQXmdP2ZRwV+8kHFYTljam9peet7whrfI4OupHbYUZzY7cWjyfsbgOb2RQHnv6J+EOhHHOstwGxmVrW/CiTFP5qQNDiEQFrjh1QU7nktjQUYvNnfX46QvKXBuMtpwY8+Lpg0P49MXL8jqOwUC4esNC/Hm3vFTVx/cNYVlTJZY2KYtHSFywogk7jjvT/n1X7wQaqyzoqM+tJ9YpbTU4kMbN+mqXEw6bKW2Qf+Oi2pT1YEII7OqbwGkdM5X0KW012HuSlYcqaJFtte/kJCotRiyZlZm0oaMWu1No/V19Ezi1oyZtq5BNi+pwbMSTsibj+Jgfq1sdObsANnTU4oTTi0lv6pjKy8fGsLrVgcYc9x84f3kjvMEwdsro5PvEviFcvLoZFZb8O+J88oKl+NPb/RibCuR9LEYd7n21B6tbq3FGDi7W2Vy5YSEODLhwbCR7JtTTB4awZa1yq0PizCX12Nkznraf2v4BF9YurMl5DK5orsIJpxeBFAusPf2TOLW9Jm0T0/VtNTg6MgVPYGbAvX/Ch9Gp4Jy2LusWVsuuHVODslYeWtA95kVnU+UcZbCiuSrlw76nbxLr22rTHm/twmoYCNiX4ks/PuZL24BNDitbqkCEtEHzt06M5xSMl6itsGDjojo8eyhzunEoEsVfDw1jSwrfbi5csroZdRUWPCTTZaaEcCRa1tXsoUg0ZWwuHwLhCH77ei8+cs7inCfZZFa2OLCsqRJP7Mtc1+MJhLG7bxLnL5efxj6bM5bUYSoQTpuEcWDAhbULch+DK1ocEAI4MT53obOrdyLj3LAubpHMnhuODE/BYjTM6fvW2RhTVIV6fll5KKRnzIvF9XPrIZY3V+Ho8NScL+7YyBRWNKc3qSssJixrqkqZntjlDGBlHo3VKiwmLK6vSDsw9vRPYn2eTQ7PXdaA17vTm/0A8HqXE95gBJesmps0kAtGA+F9G9vw4M78KpIl3joxji/85i2c8W9PYvnXHsWqf3kMV9/+Em7/61GMlrh1I4TAs4eG8dlf7cSZ334KK772KFb/y2O4/Hsv4MfPHZuzqs2FFw6PwheK4Mo8XZLJXLiyCS8ezbwokWIJmxbX5nyeluqYyzlVbFAIEbc8clce9ZUWNFRa0OWcmQAQjca2sN2Qpo0KEAu4t9fZ59SidI14sLihYo7FsrihAt5gBCMFembLuquuFvQ6vSl9lMubq+DyhzEyFUCzI9aKYSoQxrA7gM6mzMV3axdW48DATOsgFIliaCqIJQ0VaT4lj9Wt1SkD+cNuP4ZcAaxry31gADG320+eO45AOJK2IOyV42M4tb0GNRVz25DkyjWnt+Mnzx/Pq3NpMBzFv/55P365owfvWN2MWy9fg446OzzBMN46MYFfvXoCt//1KD5/6Qp84oJOmGW2hykW9p904et/2ou3eyewZW0LvrJ1NRbVV8AXimBnzzjufqkbv3i5Gz/40KacMvokHt59EpeuaobDpt73e/7yRtz36gn4Q5G0+8i81u3EKW01ebtC1y2sxoEUY2TQ5ceEN4S1C/LrjLusuQrd4zOVR/+ED55gBKuzWDWdjVXoGp2lPEY9WJpiTlkUz7zqGfMm5iAtYeWhkB6nB+9eP3eFJcVAukenv7juUQ8AYGmWtuKdjZVz+uycnPAhKpBIxcuVVa0OvJgiIL/vpAsWowErmvMbGBsX1SIYiWLfSVfaLJvXu504Kw/3WCoSro39Qzkpj2A4ik/e8wb2nZzErz65eU4F/6WrW/DFy1bit6/34j8fP4hnDg7hBx/ahJZq7Qdlvggh8KvXTuCbD+/HRSub8MyXLsaiWYuQi1Y24VMXLcW3Ht6P6366Az/9yBk5Zd35ghE8uX8I//WBDWqJDyBW8xGNCrzRPY7zU2QqAjGLNpcU2tmsbq1OOUaOj3hgNJDszgvpWFRfgQHXXAVgMhA66jIH4jsbKtA16p0p1+hUSneXzWxEa7UN3aOevNzRcimtpZTORKICJyf8KTMvbGYjGqusGJic7s1zfNSD+kpL1jYcnY2V6IorGolepw9GAhakaSgnlyWNFTjh9M55/cCACytaqhTXd8ymtsKCZU2VaVOCg+Eo3joxkXNtRyYuW9uSMc0yE//yx73Yd3ISD95ybtrWL0YDxYoh//YChCIC7739JRyXEcTVE6nJ5Tcf2o+vX7kWd9xw+hzFIVFhMeG296/Hx89dgs/c+2ZO6c8vHxtFOCpS1jHlQ5XVhDULqtNuqyyEwN6TsYBzvqxZ4MChQfecbgw9Y16019nztjjb6+wYdM9sq9M16sGihoqszU5jc8PMZ6571IvOxtTfaUuNDcPuwritWHkoYNwbRCQq0pqEC2ttMxq79Tq9siyHpY1VcHqCM7KiTji9aHFYZHfSTcei+kqMuANzWiR0jXhyTm+czfq2GuwfSB1X2XdyEoFwVJUsnNm8c20L3u5NX4Wbjkf3DOCBnX244yNnyFpVttXa8Zubz8aaBdX4m5/skNUPSQ+EEPjmw/txzys9uPNjZ+DDZ2cPYBMRvrJ1Nc5Z1oC/++3binf4e/7wCDZ31itunyOHU9trsCtN3ULfuA9ufzivYLbEqtZqTAXCc5oy9ox5Eq6gfGirtWPIPTPjsWvUg04Zz97ihkr0jvsSsdRoVGDY7UdrmtYvzQ6r7ALLfClr5aF2nYcUPG2sSm1JLKyxY2Bi+osbcvmxQIabY0l8FdE1Nm199I57saBaXuPATEgPv7RLn0T3mAedecZTJFa2OtKuWvcPuNBRb0ddlg10cuG0jjrU2M14UUFrbW8wjK8/tA+3XLRMUTGbzWzEjz98Ok7rqMHHfv6aYoVVCH703DH86rUTuPvjZypqN2MwEG57/6k4OenD3S91Kzrn80dGc3J3yWFDey129U2kzB46MOCC3WzM26UEAB11dpgMhJ6xmRZ695hnTkp+LrTV2THqCc1QzF2jHlm7ZLZU2xAMRzEeX1g6vUGEIgKtaeaVZoeVLQ81ULvOY9QdhNFAqEvjhlpQa5vhthqc9KOlOnsNhcNmRo3djJMTsz5blX8AsrHKggqLET1jM91iXaNeLFFpi9fVrQ4cGZ5KmSt/YMCVV7pxJowGwubOeryaochrNj+PT46fuUR5MZvFZMD3r9uIhkoLPvGLNxT19dKav+wewH8/fgjf/ZvTsHmp8qLPJocVn790BW5/9mjGXmvJnBjzomvUgwu1Uh4dtRhxBzDkmjsZ7h9wYfUCR9oaCSWYjAa01dnR45w5RnrGvFiswgKro64CUREb0xJ94/K8EtI+KNJnh+JWRbp5pdnBbquiZGTKj4ZKS9qCv4U1dpxMtjzcAbTIjFm0VttmPFwj7gAaKvNXHkSERfUz4x5ufwijUwHVlMfKFgeC4Si6ZykoADg44MaaVu32cT57aQN2dKVu4TAbfyiCO54/js9fujznDJ0Kiwk/++iZGJz04zuPHMjpGGrTPerBlx/Yhb+7bGVeFfzXb14Em8mI+3ackPX+V7vG0FhlzZiKng9LmyphMhAOp6hTOjToVnVRMnuMCCFwQqbbORtSksVA0vgedgVkJV/UVZhhMRkSSmPI5YfVZEi5gRoANFdbC2YVs/JQwKg7mLEau6HKAqdnOjA27PKjRWbKXEuNDYOumcqjrkKdZLiFtfYZD65knsvxucqhrdaOKqsJh2e5roQQODjozpqOmA+bOxvQM+adYfGl4+FdJxEVIuNOc3JocljxvW0bce+OHtk7vmlFIBzB5369E5sW1+Gzl2TfzCsTNrMRN5yzGPfu6ElbcZ3MzhMT2LSoVpXCwFSYjQYsbqhIGWPqHvNmzWJUQkd9BXqTlMdUIAxvMJJ3wgoQs1gdVmPC7e0JhOEOhGV5JYgILdXWxNwwFFc66e55s8OKEbY8io/RqQAaM/SZqq+MKQ8hBCJRgWF3IO32m7NprbbOtDymAmhQSXnM9oMOTvphNxtRq1LdBRFhccPcrK6BST+mAmGs0tDyWN3qgMNmyto9GADuffUEPnB6uyotUs5Z1oDPX7oC//SHPboWEv7f00cxOOnH//7NaWktYiX8zRkdGHb78XyWTcqAWHHlJg0SIZJZ3lyFo7My3IQQsWC2SjE7IGZ5JMc8JFeZWvUS9RWmRDsdyYqQPzfYZritMimd2goz/KFoQVyqrDwU4PQEUZ9hwq2vtCAYicITjCQys5pkNjVsrbEnVhehSBTj3sznUsLsDIwhd+wBVHPF2FE3V3n0jHlBFPubVhgMhFMW1mBPlg2Eep1e7OqdwDWb8rM6kvnsJcvRUm3DNx/er9oxlXBw0IUfP3cM33zPKbKfs2w0Oay4aGVT1saELn8Ih4bceXXQlcPy5iocm2V5jE4F4Q1GVIlHSCystc9YvA27/SBKnxyjlDq7CSNTMa/EoEs6trzvrL7SkthBc8IbShtzBZBwZ7l88uJW+cDKQwFufxjVaXyNQOxLBgDnVBAT8ewIuav7lmprYkUSs15iqxU1aKqeGUQbcgXQrHKx26KGCvSOz3Qd9Tq9WFhjz7uWJBvr22uy7j732N5BdNTbsS6PVhOzsZgM+I9rTsVfdp8s+P7qkajArQ/uwcWrmvDu9bk3BkzFu9cvwJP7BzN2S97TNwkDUd7tbbKxrKkKx0ZmxtJOxAPbaqTRSjRVWeOZTLFrHnYF0FhlzTtVXqIuyfKQji23fqTWbklkW036QmnjHQAS89MkK4/iwuUPoTpDC4aE8vAGMemLrRQyfdEzPlthwXg8XiL5LOvsKrqtki2PSb/qldIddfYZPmMA8YBjbq2slXBKWw329rsyNoR7fN8gtq5rVd0/v6GjFh85Zwn+9c8HCrqL2wNv9uLwkBvfvPoU1a/pHWta4A1G8FpX+iy2g4NuLG2s1KS+I5n2ugqMTgVmuGFiXRysqrgfJZqrrRACGItbB0Muf85bIaSizm5KuDdHpwKKOlnXVpoxEVcGE95gxgVpDSuP4sTtD8NhS//A2s1GWE0GjHuCmPSFYDMb0vZ7mk1thQUufxjhSBQT3hCMBoLDqs7AbHZY4QlGEk3whtx+tMoI1imho74C/eO+GVW6J5xeVVeH6VjfVoNJXwh946mD5lOBMN7qnVCtMeNsvnjZCjg9Qdy7o0eT48/GGwzjf544jE9duAxtteor5xq7Gad11OLlY+nrZw4NurBSw1iWRFu8fUdywkf/hA/tWdp6KEVy+0kLt2G3vGwoudRXmDE6lex6ku+SrrVPu62yWR5WkxE2s4GVRzqIaBER/ZGI7iKiWwt1Xpc/lLH5GxGhtsKMCV9MedTa5ftL6+JpuS5/OH4ek2orSslFJQ2MIZlpgkpor7MjGInO6OhZKOWxuL4CVpMBR4ZTFyq+2TMOAwEbNfLP11ZY8LfvWIHvPX0kMci15GcvdAEAPnlhp2bnOGdZA15Osz0yEE+VzaPjs1xaHFYYCDNqoIZcftnBZrk4rCbYzIZEmuvoVAANKha2OqzGRBxiwheU7ZEAYq5vyQ0+4QuhJku7o5iyKUPlEZ/wh4lo76zXtxLRISI6KkMhrAfwgBDiRgAbNRN2FrGYR2ZTudJqwlQggglv5hXCbKQg2Lg3CJcvs3tMKbVxOSTTV6nZLIcmx0wFBcQG/EINVsazMRgIS5uqcGx4bp0JAOw4PoYN7bWaulhuOHsxaivMuON5bXc4HJ0K4MfPHcPfb1mpqttmNucsbcDuvomUm5RFowKHh6Y0zaKTMBkNaK2e2fZn2B1QvWssEc0osJv0hVTtilBpMSSKLye8IUWZjnUVZozHFyWuLJYHELMcy9XyuBvA1uQXiMgI4HYAlwNYC+A6IlpLROuJ6M+zfpoB7ABwExE9A+CxQggthIDLl9nyAGIrmCl/OKt5ORvpvRPeoCwlpYQKixFmI2HSF4IQApMKH145VNtMsJgMCeURiQqMeYIFaQ0NxArK0u0890a3Ot1XM2ExGfDZS5bjFy93a2p93PliF5od1rxrVbJx2qJaAMDe/rk9y3rHvfCFIpp1DphNW50d/UkuyWGXH80qu12BWGbVaJLyUDJ+s1FlNcLtDycdW75iqrFb4I67tGMejcxyVdtNiXNpScGVhxDieQCzI3FnATgqhDguhAgC+A2Aq4UQe4QQV876GQbwcQDfEEJcCuCKQsjtD0URjgpUZ4h5ADHLwxOIKw8FE7TNbITdbMS4JxRzW1nVe3CJCDV2Mya8QfhCEQQjUVUHhnSOpqrpAiWlqcr5EsvKmas8olGB/SddqnRfzcb7NrahvsqCu17s0uT4k94QfvlKD265eJlqWUDpkDYp25diT+zuMS/MRkrEI7Smpdo2o2p6yKW+5QHEFnCSdTDpDWXMrFRKlcUIbzAyrQAUzA1SnHV0KtbXKptcdosJ3pD2yqNY9vNoA9Cb9HsfgM0Z3v8YgO1E9CEA3aneQEQ3A7gZALZs2YK+vvx2nRv1xB4q7+QY+vpSu0cAwChCGBybwFQgAhNB0XkdVgO6+odwcsQLM8JwOuX3bMpGhYnQMzCCg7bY5O53OTNeRy7UWglH+ofRt4BwdDS2Ugy7x9AXzpxGm0yu11xnDODwoGvO/e6bCMATjKCOvHk/A3L40IYG/ODF47h8mU12woPca7779SFUmAlnNSt7rnJlaZ0Zrx4ZwNbOmavk3cdG0VJlxsDJ/pyPreR7tiGEvlEP+vr6EBUCI24/jAG36vfALEI4OTqJvr4+OD0BRH1zn6dcifhjC5tDx09g1OVD1D8l+9juidiY3Xss1jZmanwUfZS+szOFgxh2Tqoie3t7egu3WJSHIoQQewF8IMt77gBwBwBs375dZLoJcgjEV7XLFrVnDNY11YwCBgJFw2hyWDPe/NnUVByDtaoawhlFSx2hvr5e0ecz0VjdA4O1Evaa2N4Vqzs7VO9029YwiCDFrvm4dwRmI2HNssWKK59zuebTUYXJJ06gurFlRrxo9/gAHFYTzlyzVJUK7Gzc1LoQd785gpdORvCJCxbL/ly2a/YFI3hwz3588bKV6Fy8KF8xZXHWihDue7Vnjmzu3W4sa6nJ+9mU+/lFLT50TY6gvb0dI+4AIgJYu7Qd7SrHXFobJjDk8qOtrQ1Tgd3obGtFe7s6TR8HXUEAQ3DUN2MqeAidC5vR3i6vD5m52g/gIIyVtQCAJR0L0Z6h8LahZhREpNrckY5iybbqB9CR9Ht7/LW8ULMlu5Rnbk+zJaZEZTzm4Q1Gsr53NhVWE6YCsWwrNU1mAHG3VSgRSFP7+EAs3VHKthpxx4LyhZiwASQC8ydn7clwYMCFNQuqCyaHxWTAhzcvxi9e6ZbVH0ouD+86iXBU4NoztJ0QklnV4kDPmHdO/Ypa+1zIpaHKkiiwkwLH9Rq0+JcCzZ5gBOGoUDUuWBm3Ql3+EKYCYUVb9lbEEz2kGpRsiRJ2ixG+AritikV5vA5gBRF1EpEFwDYAD+V7UDVbskvKw2bJfMuqbCZ4gmH4QhHF2T1VViO8gQhcvsz1JLlQW2HBpC+ECW8sDViNVtZzzzGdUjiiQUZXJhoqLbCaDHOUx7GRKSzTqOtrOq7bvAhDkwE8neMuh7MRQuCeHd249vQOTTOsZrO0qRKRqEDv+OziT19hlUelFWOe6ToHQH7xrRJiyiOsyTkqzLF5Y2QqgKiAorlB+s6lIsOKLJ+ttJrgDZZhbysi+jWAVwCsIqI+IrpJCBEG8DkAjwM4AOB+IcS+QsuWCX8oCiLAkiVQWWWNZTr4QzlYHpaY4vGGwqhUeZJw2Exw+UNwKQzWKSE5RXDcG9RkA6h0EBHaamdm5QCx/lpLVOyBJIfGKiuu2rAQd7/crcrx3uqdwN5+F244R74bTA1aq22wmQ3omtUepK9AnQMkGqpi2UaBcAST3hAqLEZNWt5U281w+UKJHT3VVB5GA8FuNiZZD/LnBqOBYDMbMDYVABFgzXLtdrOxPJWHEOI6IcQCIYRZCNEuhLgz/vojQoiVQohlQohvq3Eutd1WNpMxa+GezWyEPxSJua0UWh6VFiM8gTD8oShsZnW/Gns820PtTK5kkpWHyxfOmpmmNm11dvQn7acS676qzoY+SrnhnMV4+dgYTszanS4X7n2lBxesaJS185yaGAyEzsYqHE/aQ1tqJ55uG1QtkIr1nPHODVpYHcD08ytt2ay2lWe3GBMtiHJZWI5NBWE3Z5+DKizGOdtOa0GxuK00QU23lVw3lNVkQDAchS/HmIcnGIEvGIFV4WezHttsgi8YU2pKVj1KmKE8NIjbZCO2Gde05eH0BDEVCGNRfWEnXQDY0F6D5c1VeHBnfhkvnkAYj+4dxHVnFSZIPpvOxgp0jU4rQKmITs2+T9mojRfQTvpCsQprjZ6rSosJ3mAsXmkykOrWjc1kSLjflI7BCosRo56grDmlwlKmlkchUdfyiMIm42GymY0IhKM5xjxM8AZi5rlSxZONSmvsgcpFLrlU281w+8OIRIXqVfJyaK2ZuRtjT7xRo5r7PsiFiPCB09vx4M6+Gf2+lPLY3kGYjIRLV2vTlysbbbV2DE7OLNAD5LcTV4Mqa8wCyKX4Vgl2ixFREVNSWowRW7LlodgrEevKK+dzFZbYQlFrylp5qG152GRaHoFwJDfLw2KEJxD7rE1l5RHLwMhNLrlIg9rtD8GlcpW8HBodVox6ptuj9I37UF9pSUw+heZ9G9twcsKHVzN0p83GH9/uxxXrF6j+PMiltWbmLpQjU4HE1qiFwmIywGoywO0Py2rPkSvTWU0BTcaIzWSE0ysvY2rOZy1GONnyKE0C8ZhHNqwmA3w5rvArLEZ4gmH4w+rHPCri8RRvMKyp2wqIrdzcOlgejZWWREASiLexKKB7ZTYt1Tact7wRD+3KLet8yOXHS0dH8b6NbSpLJp8Fs7ZHHtaoujsbDlus+ltTyyM+MY95gpqMESnmYTEZFGc7Wo2xTrly5LLEXedaU9bKQ/WAuYwJ3Wo2whOUVxMy57Om2IohEhWqr3wkUzYWyNdmJS7FOFw+bWpVstFQZcWkL5QYOCPuQMHao6Tj8lMW4Mn9QznVfDy6ZwAt1TacuUTbvlyZaK2xYcIbSrhBhnW6p9W2eA2UT7vnKmF5eIKajBGb2SDbepiNxWSAV6ZHwmIyIBhh5ZEXegTMk+MiSk17i8mQaNustpuiwmKEN+620srykNKLvcGwLtlWDVXTnYkBbbqvKuWd61ow5gnK2mN9Nk8eGMKWtS0FK3BMxYJ4NwXJ+hhxB3Sx5qpssRR4TzCMSo2eX/sMt5X6U6PdHHNb5TL+zMbYMyBnTrEYDYhEhapFqqkoa+WhJrGAuQy3VdKkn60mZDYW43TbZi3cVpGowIRM0zcXjPEMlXFvEMFIVFEVrRo0VsYmtdGkKnctuq8qobHKijOX1OOxvYOKPjfpC+HV405sWduikWTyaKqyggiJRIRxb1CT6u5sOGyxzg2+UFQzy1mKQzg9QU2KMa1mIya8uQXjJaUhZ+ta6b1au65Yecgk5raSF/OQUNr51GIyIBQR8eOoHDA3xwbD2FRA0+Cr3WxMpCNWqrQTolyq7SaYDJSIewy79Y15SGxd14rH9w1m3CZ3Ns8eGobdbMTmzgYNJcuOyWhIdGQG1G9VLheH1Qy3PwRfMKyJVQBML35ibisNYh7xcZfL2LbEP2OSYYVKCkZr11VZKw81Yx6hSDRhOmYiWXnIeX8yySap2g+vZMloaXkAMQvHGZ+85VhqakJEqKu0wOmZdlvpHfMAgItXNaF/wofjo/K7GD91YBgXrWoqaFZTOuoqLIksoUkNYw6ZiCWTxBJRtGzRIlWBa5JtFR+DFoXzQuwzcctDxvNgZcsjf9SMeYQjQpYlYZmhPJTd3mTFk60FgVIkWeRmbOSK3TKdjqjlzn3pcNhMcPtDiEYFJn0hXVwss+lsrERbrR0vHkm/J3gy0ajAy0dHcbFGe64rpS6pZ5leloc13rnBF5SXMp8rdrMxVuehgfIwGWJjMJe9WBJuKxmWR8JtxZaH/oQjUYSiQpYlkaww5JiYyeSjeLIhKSMh1HeJJVMRz0cHCm95AEC1zQyXP9ZCQwgUPF04FUSEC1Y04gWZyuPQkBtjniDOWaavy0qirmLamtNLedjMBvhDsc4NFRq6Xc2m2Ji1auAak8a30nkBmB6/Shaw/lBEkatUKaw8srCrdwLLv/YoRtz+xMohE8n523JMzGSsGiqP5OOpfexkKsymaeWRpQOxFjjiWTkuDbuv5sL5Kxqx4/gYQjJWgy8fG8OShgq0FWD/dznUVVow7gnCH4ogGI7q4raKdW6IwKthhwRgemzIGetKkZRGLq5IaeGqZAH77b8cwKX/85zic8mlrJWHGjEPKcXy0KAbJjlfXNJDZ1b4AFqM04NC7exMy4xAvnapn3aLEePeoKwOxFpQbTPHK9y127ckF85b1oipQBh7+rPvqvjKsVGcs6yxAFLJo77SgnFvUNN26NmQ3ElCYTtzpSRiCxqMEZMxd8sjl2yrZw4Oo0tBnE0pZa081Ih5BOJBp3FvSNYXZ0x66JRO0skTfLbOmUpJll3LSV0KmMvp/qkFkuUx6QuBCHDo1JpkNnWVFixtqsTOLPUe0ajAq13OonFZAfGNxHwhXZWHVGAHKC++VYLZKH+SVnxsg2Q95BDzMErZVvLqPApBWSsPNfAltTaW01IgeVWRT8xDbcx5KDUlWE2xNgp69WKqtpvj+5aEUWU16VpgN5tNi+rw1omJjO85PuqB2x/GpkW1BZFJDtJWAZIrsNA9y4CY22o8x460SkjEJTS0PHJRHiYFbqvZykPayE5tWHlkwZPUYExOpkOywlC68tZSeRBRkt9USyVlgEfD5ovZcMQ343L59QnsZmLTojrsPJHZ8tjVO4HGKkvRxDuA2M50nkAEnmAEZiNpmnCRDptpuu2PlgsTLceIkrjFbAzxuUSOUjMYaMY85Alos7cHK48sJG+qIifTIZ/tXXPxhSpBUmZa+HMlJAWoRbaKHKRKZD1awmdj0+JaDEz6MTDpS/ueXX0T2NBeq4vLLx1V1vjWyhp0e5ZL8vOkpVvGrGHMIxGMz8XyUOjySp6HtErZZeWRBU9g2vKQo/XzGfQGjScM6ehaWx6Atn7pTEit511+9feBz5cVzQ5UWIzY3Zc+aL6rdwIbOmoLJ5QMYpZHblsrq0Wy0lKaxagEi4YxD1MeVo1BofJIztDVqliwrJWHGtlWvlCy20rb26W95SGdR7vrkNKN9aqMtpljysMf0q4BZK4YDYSVLQ4cGnSn/HsoEsWBATdOba8psGSZqbSaEIoIzQtMMzGj7Y+G4yQf6yDrsQ25WzXSNcu99mRrg5VHDqiRbRVO+hK0DDQDMzO1tICgvdsqYfZrrGjTYTMb4Q9GZPciKzSrWx04OOhK+bfuUQ+CkShWt1YXWKrMSJtpjbi17YuWiULVKSmp5FaKGpZHLkotwMpDH8JJbY21WI0kY9TYbRWN27KFcFtprWjTYTcb4c9xJ8dCEFMeqS2PI8NTcNhMaNG5E/BspAaXozK3QdWC5BV3PnHFbGiaqpvH2JiOeSj/LCsPnYgmOQ+1WI0ko+WgAADpSrSc2C0K2ihogd1iRCgiMBUIa9oDKVdWtVaje9STMn3y8JAbK5qriipYDkzv0zI6FdDNbVWoxYjFJD+rSSmJPTlyGBvSwjKXOYLdVjoRjhTQ8tBYeUTjVpS22SrKfLNqI/XTGvcGdemtlY3VrQ5EBXB0eGrO344MTWFli0MHqTIjuarGvdo0DJSDlnG6ZKSkFS3GiFFqjJjDtUhuq1ySajjbSieSd+PSekIslNtKSyWYT/M3NbDH+2nFNt0pvse7rtKCGrsZ3WNz20YcGXZjeXOVDlJlRvpOJ7xB3WIekiWgtVFmzCO2kA1pSORyaGk85XL5UY12FCy+0VVkRESy5aHtk6t1NbT0DGla56Ghz1gOUgGb01OclgcALGmoQM+Yd8Zr0ahA95gXy5qKT3kYDbEC00lfWDe3VaGeJ4OGtVDSsXNxS+ZjeWi1HS0rjyxEoiIxIWptGRSKQgTMtXbBpUMK6Oa63WchWNRQiZ5ZlsfIVADBcBQd9cVTWZ6MxWjAhFebTZLkUChLNmF5aOAmk6aPXKaRhOWh8LNmI81YAKtJWSsPNeo8IlGRMNuLqU9SPmhZjChZZ3plW0lulWAkWpSpugCwuH6u5dE3Hvu9rbZCD5GyYjUbEY4KWHVO1dX6qZKUhxYhFmMe1kOuVovVZGS3VS6oUeeRrDzKQ3Wo3+49GWmA6FXnkVxMVrTKI4XbqtfpQ2OVtWitJavOsaxCLUYS1oEGo11SALncwmnFo+xzFpOBLQ+9iCTtIKh1+5BCoeV1JFIKdbI8kic3m079tbLRUV+BQZd/Rgpl37gX7XXF6bIClO1kpwWFzrbSQkdKwy6X8WfK0WoxEHHMQy9muq10FiZPpAGhpQ5M9ODRaYWa3D1Yj82o5NBabQMQi3NI9I37ilx5xCwii06LAi2TPJKRzqKFizqfRZshx5iH0TCzVk1NinN0FRERIZL8raVteeST7SEXyfLQa4UKTK9S9cr4ykZLXHkMufyJ1/onfGgrYuWhd/GndF6tCyinrQP1jy25nnKZy6cXfsoEMxJBozIPVh7ZCEems61K3WtVCPmns1X0u1lmnYP22bBbjKi2mTA0Oa08RtwBNFUVV1uSZHSPeSQmXm1W0RLSAlELJZXPrUvIpfScBmLLQy8iUZEYOKUe8yiE5TTdwE1P5aFvurAcWqptMyyP0akgmhzFqzyU7KGtBQnlofF5DAnLQ/1nJx+FlGu8xEDE2VZ6McNtVbxzkSwKYnlIbisdA0TSdRar2wqIKY9BVyzmEYkKOD3FbXmY8mjqpwYFWwhoGDCXJn6RgwrM1Z1mNHCdh27MCJiXuPYohPiSztDTbSUZ9/rKkJmWahuG45bHpD+MqAAai9jyULqTndoUqllkImCuwfmkhVUuc/m0O03Z5wzE7UlmQERrieh+IvoREX1Ay3PFUnXLw/IohPIzUm4PuZoUYtOrfKmvNMPpDQIAnN7YVseNRWx5JOp3dI4jaRzyyKsKXO6xC/lZo6GMUnWJ6C4iGiaivbNe30pEh4joKBHdmuUwlwP4PyHELQA+opmwiGnt6aZkpa09CiG9NMlo9LwqolgD5gBQW2HBuDcEIKY8TAZCrb249lxPZnonu+JVyGowHZgurlRdSvE/uefUaizqscnz3QB+AOAe6QUiMgK4HcAWAH0AXieihwAYAdw26/M3AvglgG8Q0XsANGgtMGnoBy03DAnlob/20HuVnIm6Cgsm45bHuC+M+kpLUbe/KRbLQ2sSsQUNdKR0zFxGBuVh0Ws1EguuPIQQzxPRklkvnwXgqBDiOAAQ0W8AXC2EuA3AlWkO9dm40vm9ZsLGmX6gSnvgFGI6z8evqxbSt2Qs4lVybYU5YXm4/WHUVhSv1QHMH8tDem61cPEa8hgbuU49RKRZerMelkcq2gD0Jv3eB2BzujfHlc8/AagE8F9p3nMzgJsBYMuWLejr68tJsKiIIuCPBTZHRkbQZ/Jm+cQ0uZ5T+qzT6cz586mQAmf5yJWN0ZHY/Zl0uXI6jxrXHI3GdukbHR6C2W/J+3haEPZMweULoedEL4YnpmA1RDX9XvIl4PcBACbHx9DXF8r7ePl8z1rep0nXJABgeGgIFaFJ1Y7rdDrhpniCxOSk4msYGYl1YR4bG0NfX1j258KhICYmlJ9Por29Pe3fikV5KEII0Y24YsjwnjsA3AEA27dvF5luQiYMtA8VdjuASbQ0N6O9vU7Gp3YByHzj5X42V7lTQnvUP+YsXAYXgCNwOBw5nydf+QyGgwDC6GhbiNYaW17H0gq30QWBY6huaEHUOIDGaqum30u+OKrGAIxjQUsz2tsbVTmm8uvNZ1zJw3HAC2AICxa0ol3lvVWqrLUADqGmpkbxNQxFxgEcRWNDA9rbW2V8InavLBZLTueTQ7HYoP0AOpJ+b4+/lhdqtGSPHSf2b4l7rQqCsQhiHlKAsLgD5jE31bg3iKlgBNVFHCwHpgPIeradKQTSU6up2yqPOo9cPpfL+eRQLE/C6wBWEFEnEVkAbAPwUL4HVaMle+w4sX8LlWuuFYWYz6W5Rc+Yh+TjLeY6j5q4spj0hTAViMBhKw0nQDFX7atC/NnRtEgwpzqP+L8K5yACaTYW9UjV/TWAVwCsIqI+IrpJCBEG8DkAjwM4AOB+IcS+QsuWiuT7XurjphDzuTRA9EzVLcRe7fliNxthIMAbjMQsD1txWx5SnLzclYeWlkc+NSS5yqPlelePbKvr0rz+CIBH1DwXEV0F4Kpt27apcrxSrzAvhPaY7hxaBG6rIp7oiAiVFhM8gTCmAiXgtpon6erTXgbtzpFLDUk+8mg1Eot3aaYCqrmt4re/9HWH9hN6YmDoeK8ky6PYV8kVViM8wXBJuK0SLcFLvFBWLlq4qPOxPHLtqkvQzoWcVXkQkYGIztXm9NqiVsA8cbwSHziFMAaKQcEm6jyKQZgMVFpN8AQimApG4Shyt1UxtJ0pBNICS5udBHNTAEAeCz8i/QLmQogoYtXfJYfaAfMyr49SlWJQtMU+0Uluq0A4iooi3W9dIp8K51JCyyLBfI6YqztNV8sjztNEdA2VerpRjmgZRCskhYhCFMMtmp7oikCYDFRYjAnlYSty5TG9t3dx39N8kcaIlo0RC3kPNd1yWub7PgXgdwCCROQiIjcRubQTSx3Ud1uVNoUIYhfDCrVU5rcqqwkTvhCiArBbSsOsLXvloaHlIVHIWxizPHSs8xBCOIQQBiGEWQhRHf+9WhOJVEQNt1XyjS/2lWw2CmJ5FOAc2SiVCa7CasLYVKw5otVU3JaHRInc2pyZjnlo4bbK/Zi5WkSx3lY5nzYjslM84h1sL4z/+qwQ4s/aiFR8TK9G9JUjXwoZMNfzVpXK12Q3G9DrjPWMKna3lUSpj4GsaDjW1dBHShUQQedUXSL6dwBfALA//vMFIprdKr3sKXXLoxBwoFw+VpMRLn+syaDNXBpuq3IfA9MrfO0C5rkcO1fXUzHEPN4NYIsQ4i4hxF0AtgK4Qjux1EG9mId26XvlRnHco6IQIisWkwGTvpjysJeI5VEadzZ3hIbtSfQqgdI72woAapP+X6OyHJqgeqpuma+6VIFvkWysJgNcPsnyKA3lUe5joFz62EkQtKvzkBvz+A6At4jor4hNDxcCyLZVLDMPKQa3VXFYP9mxmAxw+WN7M5SK8iiTOTUr2hgeUlfdAkLaWR5ZlQcRGQBEAZwN4Mz4y18RQgxqI1LxUuo7CRaCfFowqC1DsSNlWJkNVPStVCTK3fLQtKeVDrdO14B5vML8y0KIASHEQ/GfklAcasQ8km98eQ8bdSiGe1QM1o8cLKbY8LOaSkNeoHQUc65o6a6SjpxL8DvxCcWputDM9JAb83iKiP6BiDqIqF760UQiFVEr5iFR7qsuNZju36PfvSqVr0lSHtK/pYDesQCtT6/l4fW4d7GYhzbIjXl8MP7vZ5NeEwCWqitOcVIqk1ExwLdKPta40ijm1vGz0VtUzU+vaSv2wkNFEPO4VQjxW21EKAVKZ3BnQssHKfkcjDwk5VEq8Q5Af+u7UKv3YnuOE1lgCj+n6za08ZjHP2py9hKjEPthlDp6uzUAfbfAVUJCeeh/y2Sjt6ha61kt3a3S0Cjk86nl9ZR1zEMtimA+VIVCXEa53KtCIMU6TCWkPfReHGgdSyvEDoKFXoTq3duqJGMeam9DWyor2nRQAfxW0y0YND1NRkrFQjQaJMujlJSH3gKU8OETFebzqCW7EKIzxU9RKw5AvWyrRIpd/iLpSmEsD71nl9JR8lKgvLQC5npbHhofv3S+Ct3JqDyI6MtJ/7921t++o5VQxUTyRMTPVXZo1r96UCK6IxEoL6Xi0xISNS+Kr1Yo9lQXw+JMIpvlkezv+eqsv21VWZaipYi+r7wo6HWUy03TkITy0FkOJRTfpKouhbi+chka2Z5bSvP/VL+XLaXiBskGD4ziQlIepXTPqJQ0XQ6U0nehN9keBZHm/6l+L3tK/sEqoPy6uq1K5MkspVhHNNGqvHRkzoVivbpc6zy0JFu21Yb4XuUEwJ60bzkBsGkqGaM+hdhJUI/OoXMoDe2RsDx0lkMJpSRrPmihI8tN72ZUHkKI0ugTzRQNxTBASsXySFSWF8E9y4aYJ5ZHUTzAJUJZezDV2kmwROYipsQwJSyP4p+woolNkvSVY74yvT2urmLMoKyVhxp1HskFZ6UwyJnSUfZSkWApIBKpojoLwhQNpfP06sh0kWCpTEupKaj8OvqOctkvQQ+kyvJSmJDnzVbMJfLsFAOsPBhVKfe5RU2MxtIJmCfcJrpKMX+ZzrYqnm+AlYcCiumLy4VSl18upbJ2LKVU3XljeTCyYeUhg1KZjIqBYlBQpeJ5kLKtoiUhL8c8mJmw8mA0Qc/5sNRiHqUiL1BcvZU0ocivr5jEk9uSnVHAf33g1KLcHa4QAfNieriLHSnmURqWxzxBQ0Wez6FP66jFx85dgjOW1KknUJ6w8sjCjK66MifGa8/o0EaYEqAYdEeprI4lyyNaQpbHfKHYHiGLyYDt71mntxgzYLeVDIrsOSoJeD7MTikFn/n7zJ8S+rplUfTKg4iWEtGdRPRA0muVRPQLIvopEV2vtQw8buRTKqv+YoBv1fyi3BSwpsqDiO4iomEi2jvr9a1EdIiIjhLRrZmOIYQ4LoS4adbL7wfwgBDikwDeo7LYTIlTKpNyqcjJqEu5LLC0jnncDeAHAO6RXiAiI4DbAWwB0AfgdSJ6CIARwG2zPn+jEGI4xXHbAeyJ/z+issxMHpTHsCgMxZDWLJdyWzUz+aOp8hBCPE9ES2a9fBaAo0KI4wBARL8BcLUQ4jYAV8o8dB9iCuRtFND1VjpDXX9KvZVLISiTBWhZooliL7MhoUe2VRuA3qTf+wBsTvdmImoA8G0AG4noq3El83sAPyCiKwA8nOZzNwO4GQC2bNmCvr6+nIQVQsDn8wEATg4MIOgy53QcpfT19cHpdKp6TGn1mOu9kHeO2ElcLndO51HjmkU0CkDb61SDSDxHNxyOFL2sHq8HgHr3NNfvWQih6b1yud0AgP7+PpiN6q1LnU4nXP4wAGByYqJg37ff74fbbcj5fO3t7Wn/VvSpukKIMQCfnvWaB8DHs3zuDgB3AMD27dtFppuQEdoDu90OwIWFCxaguVrrPbB2AZj+0nKWOwVEu1U/Zmp2o7rakfN58pXPaDwAIFKA68yPaFQA2A2TyVj0slZWOAE4VZVT+bF2gYg0vVfV+z0AhtDW1g6LSV2nhqOhBcA+1NTWFuD7js0jNpsNDkfuYzETemRb9QNILoRoj7+mOmrt51Eu3oVC+q3ZR56dUnJbzRc3pKZXWWa3UA/l8TqAFUTUSUQWANsAPKTFidTYzwMou++cKRLKJeumHNHyqymXb13rVN1fA3gFwCoi6iOim4QQYQCfA/A4gAMA7hdC7NPo/KpYHtMHVOcwejFf5qp5cpkMoyuaKg8hxHVCiAVCCLMQol0IcWf89UeEECuFEMuEEN/W8PyqWB4MoyVs2c4Pys31V/QV5ox6FDTmUbhTMYxqsNUqn6LPtsoHIroKwFXbtm3L6zgc/FWOnoOw1NxzJSaubvzre0/BmlaHpufQcqiXUlGoHMra8lDFbZXcVbfMvnyt+NbV63DLxct0lKC0vqdSWJsUwwLqhrMX44wl9QU5lxZPULm5rcra8lCLUlvJ6s1HzlmitwgMU7SUy3xS1spDLbfVVRsWosZuRn2lRR3BGIZhNOSnHzkDgy4/ntg3qNk52G0lgxq7Gf997Yai3B2QmUu5rOyKifJyuOhDIV1/W9a24IazF2t6jrJWHgzDMIw2lLXyUL1IkCkJ2PBg8kWL6v9ys4jLWnlwkSDDqMMFKxr1FqHkKYaMNTUp64A5M5Mye3bLhlKYVK4+rQ1Xn9amtxhlQbkYIGVteahBueVmzwfKzT3AMMUIWx4y4LmotPjh9afj1a4xvcWQDSu74kFLK9BojH3RDlthNpTTmrJWHmrVeRSST124FMubqzQ59nyZo05fXIfTF9fpLYZs5sv3Ukpo8Z1U28z49SfPxubOwlTJa01Zu61KMWD+1XevwbVndGR/I1M28L4e84dzljXAUCb1YmWtPJiZcPSmOCmPqaQ8YD0uH1YeDMMwcUoh861YYOXBMDrDq93ig7+T7JS18lCjwpxXIozW8DzFlCJlrTzUCpjzKoTRFH6+mBKkrJUHw5QCvMkYU4qw8mAYnWHLlilFWHkwjM6w7ig+uPYmO6w8GIZhGMWw8phHCE4dYxhGJVh5ZIGnW0Zr+BljSpGyVh5q7STI2TAMMz/gLRjkU9bKoxQbI2oJBwEZhlGLslYeDMMwjDaw8phHcMCcYTLDLmr5sPJgGIZhFMPKg2EYJg4HzOXDyiML7OphtIYfMaYUYeUhg3JJUlq3sEZvERiGKRNMegvAFI5ffXIz3P6w3mIwDFMGsPKYRzhsZjhsZr3FYJiixWhgZ4xcil55ENFSAF8DUCOE+EC61xiGYfLlUxcuRUu1VW8xSgJN1SwR3UVEw0S0d9brW4noEBEdJaJbMx1DCHFcCHFTttcYhmHypdJqwvWbF+stRkmgteVxN4AfALhHeoGIjABuB7AFQB+A14noIQBGALfN+vyNQohhjWVkGF3h9FCmFNFUeQghnieiJbNePgvAUSHEcQAgot8AuFoIcRuAK7WUJxd4WDMMw8xFj5hHG4DepN/7AGxO92YiagDwbQAbieirQojbUr2W4nM3A7gZALZs2YK+vr7cpBXA6OgI+mz+3D6fB06ns+Dn1Jv5eM3hcDj357NEmW/fs17X6/f74XYbcn6+2tvb0/6t6APmQogxAJ/O9lqKz90B4A4A2L59u8h0EzJCu9DU2IT29sbcPp8nOctdwsyva94Fs8k8z645xny7Zj2u12brh8Ph0OTceuSl9QPoSPq9Pf6a6qi1nwfDaAnHPJhSRA/l8TqAFUTUSUQWANsAPKTFiXg/D4ZhGG3QOlX31wBeAbCKiPqI6CYhRBjA5wA8DuAAgPuFEPs0Oj9bHgzDMBqgqfIQQlwnhFgghDALIdqFEHfGX39ECLFSCLFMCPFtDc/PlgdT/LDXiilBuBafYRiGUUxZKw813FbcLpthGGYuZa08VHNblUlLdoZhGLUoa+XBMAzDaENZKw/OtmIYhtGGslYenG3FMAyjDWWtPBimFOCcDKYUKWvlwW4rhmEYbShr5cFuK6YU4GQ+phQpa+WhFsTDm9EQdlsxpQgrD4ZhGEYxrDwYhmEYxZS18uCAOcMwjDaUtfLggDnDMIw2lLXyYBiGYbSBlUcGBLfUZRiGSQkrDxkQZ+oyGsJrFKYUKWvlwQFzhmEYbShr5cEBc4ZhGG0oa+XBMKWAxcR+Uab0MOktAMPMZ35yw+loMfn0FoNhFMOWB8PoyLvWtaKx0qy3GAyjGFYeGeAsGIZhmNSw8pABe6QZhmFmwsqDYRiGUUxZKw+u82AYhtGGslYeXOfBMAyjDWWtPBiGYRhtYOXBMAzDKIaVB8MwDKMYmg9tx4noZwD69JYjB04H8KbeQhQYvub5wXy75lK93m4hxN2p/jAvlEepQkRvCCHO0FuOQsLXPD+Yb9dcjtfLbiuGYRhGMaw8GIZhGMWw8ihu7tBbAB3ga54fzLdrLrvr5ZgHwzAMoxi2PBiGYRjFsPIoEYjoS0QkiKhRb1m0hoj+i4gOEtFuIvoDEdXqLZMWENFWIjpEREeJ6Fa95dEaIuogor8S0X4i2kdEX9BbpkJBREYieouI/qy3LGrByqMEIKIOAO8EcEJvWQrEkwBOEUKcCuAwgK/qLI/qEJERwO0ALgewFsB1RLRWX6k0JwzgS0KItQDOBvDZeXDNEl8AcEBvIdSElUdp8F0AXwYwLwJUQognhBDh+K87ALTrKY9GnAXgqBDiuBAiCOA3AK7WWSZNEUIMCCF2xv/vRmwybdNXKu0honYAVwD4md6yqAkrjyKHiK4G0C+E2KW3LDpxI4BH9RZCA9oA9Cb93od5MJFKENESABsBvKqzKIXg/yG2+IvqLIeqmPQWgAGI6CkArSn+9DUA/4SYy6qsyHTNQog/xd/zNcRcHfcVUjZGW4ioCsCDAL4ohHDpLY+WENGVAIaFEG8S0cU6i6MqrDyKACHEZaleJ6L1ADoB7CIiIOa+2UlEZwkhBgsoouqku2YJIvoYgCsBvEOUZz55P4COpN/b46+VNURkRkxx3CeE+L3e8hSA8wC8h4jeDcAGoJqI7hVCfFhnufKG6zxKCCLqBnCGEGJUb1m0hIi2AvhfABcJIUb0lkcLiMiEWDLAOxBTGq8D+JAQYp+ugmkIxVZAvwDgFEJ8UWdxCk7c8vgHIcSVOouiChzzYIqRHwBwAHiSiN4moh/rLZDaxBMCPgfgccQCx/eXs+KIcx6AGwBcGv9e346vyJkShC0PhmEYRjFseTAMwzCKYeXBMAzDKIaVB8MwDKMYVh4MwzCMYlh5MAzDMIph5cEwBYaIlsSLIBmmZOFUXYYpIER0C4C/BVAF4BiAbaXeLYCZn7DyYJgCQUQOxBTGVgCnAngWwBiAewE8KIS4h4g+BeBCIcT1ugnKMDLg3lYMUziiiLXVrwcAIUQ3ABDRzQBeIqIuAF9CbK8LhilqWHkwTIEQQniI6JMAbgPQSkSnAPi6EGKIiL4O4K8A3ieEcOoqKMPIgAPmDFNAhBAPAbgWwH8CaELM0gCA9Yi5sBbqJBrDKIKVB8MUCCKqIqLF8V+lnfQcRHQWYtvRbgTwD0TUqZeMDCMXdlsxTOEwA/gJgAYAjYjtSf8RAA8B+LgQ4iQRfQnAXUR0aZnuY8KUCZxtxTAFJr4F68VCiLt1FoVhcobdVgxTeCYAvK2zDAyTF2x5MAzDMIphy4NhGIZRDCsPhmEYRjGsPBiGYRjFsPJgGIZhFMPKg2EYhlEMKw+GYRhGMf8fWy0gveZXsVYAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "#compute and plot the relative error\n", "err = abs(ycheby - yexact)\n", "print('Maximum error = %g' % max(err))\n", "plt.semilogy(xx, err)\n", "plt.xlabel(\"$x\")\n", "plt.ylabel('Error')\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "attractive-suite", "metadata": {}, "source": [ "With 16 points, the boundary value problem is solved with a maximum error of about a millionth of a percent!" ] } ], "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 }