{ "cells": [ { "cell_type": "code", "execution_count": 38, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "from math import sin, pi\n", "\n", "DEBUG = 1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Exercise 1: Interpolation with Hierarchical Basis" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Solve interpolation problem:\n", "* given regular grid poitns $x_n$ and corresponding data points $u_n = \\{ u_i = f(x_i) \\}$\n", "* using **hierarchical basis** $\\psi_{k}(x)$\n", "* requiring: $\\displaystyle f_n = \\sum_{k=1}^{N} c_k \\psi_k (x_n)$\n", "\n", "Steps:\n", "1. set up linear system of equations to obtain $c_k$ for n=3\n", "2. derive formula for $c_k$ (depending on $f_n$)\n", "3. repeat steps 1 and 2 for $n=7$\n", "4. Generalize this to $n=2^L-1$ " ] }, { "cell_type": "code", "execution_count": 39, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 39, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD/CAYAAAD/qh1PAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAADvNJREFUeJzt3HGQFvV9x/H3g5xG5LhINI0lmVCJZohGq6hYQfPIlKRTZGL4Iy02jlJRTNM6hT8csaNckkn+SeMMZhqTUmqiGFrTahXpgG25p3LigFQ6YmsSodrJTM9iK+UIJ+GAp3/sHjw+PHfH7R7Ps/fl/ZrZYXd/t7tfvnCf2/vt8zwgSZIkSZIkSZIkSZIkSZKkUTQT6GqwfynwWjrWBVzczKIkSaPjXuBVYEuDsceBK5pbjiSpkXE5jt0FLABKDcZmAPcDm4H7clxDktRiU4GXGux/AJgMtAHPAfOaWJMkqUaeO/qhrATeBfqB9TiNI0ktM/4UnLODZO7+U0AfMAdYXf9F06ZNq+7evfsUXF6SQtsNfGIkB4zGHX01/XMhcCewj2Revgt4geTVNxvqD9q9ezfVarXwy4oVK1peg3Vao3Va58ACTBtpSOe9o38LuC5dX1uzf23dtiSpRU7VHL0kqSAM+mGUy+VWl3BSrHP0jIUawTpH21ipM4tGr4Fvlmo63yRJOkmlUglGmN3e0UtScAa9JAVn0EtScAa9JAVn0EtScAa9JAVn0EtScAa9JAVn0EtScAa9JAVn0EtScAa9JAVn0EtScAa9JAVn0EtScAa9JAVn0EtScAa9JAVn0EtScAa9JAVn0EtScAa9JAVn0EtScAa9JAVn0EtScAa9JAWXN+hnAl0N9s8HtgFbgMU5ryFJyqGU49h7gS8BvwCuq9nfBvw7cBXQB7wI3ATsqTu+Wq1Wc1xekk4/pVIJRpjdee7odwELGlxwejq2D+gHuoEbclynMA4dOsRdd93D5MkfZcqUT/L440+0uiRJGtb4HMc+BUxtsH8SScgP2A905LhOYSxbdj9r1vyM997bzN69Pdx99xeZMuUC5syZ0+rSJGlQp+Jh7D6gvWa7Hdh7Cq7TdE89tY733vtT4NeA6+jr+0OeeebvW12WJA0pzx39YH4CXAScCxwgmbb5VqMv7OzsPLZeLpcpl8unoJzRM2lSBz09bwKXAtDW9iaTJ3+0tUVJCq1SqVCpVHKdI8/DWEimbn5E8jB2ITARWEXy8PVBkt8YVgOPNDh2zD2M3bhxIwsW3MrBg79PW1sPH/zgZnbu3Mr555/f6tIknSayPIzNG/R5jLmgB3jllVdYt+45Jk48h9tuu43zzjuv1SVJOo0Y9JIUXLNfXilJGgMMekkKzqCXpOAMekkKzqCXpOAMekkKzqCXpOAMekkKzqCXpOAMekkKzqCXpOAMekkKzqCXpOAMekkKzqCXpOAMekkKzqCXpOAMekkKzqCXpOAMekkKzqCXpOAMekkKzqCXpOAMekkKzqCXpOAMekkKzqCXpOAMekkKzqCXpOCyBv044HvAFqALmFY3vhR4LR3rAi7OWqAkKZ/xGY+7GTgTuA6YCXw73TfgSuBWYEeu6iRJuWW9o58FbEjXtwJX1Y3PAO4HNgP3ZbyGJGkUZA36SUBvzfaRunOtBZYAc4DZwLyM15Ek5ZR16qYXaK/ZHgccrdleyfEfBOuBK9I/36ezs/PYerlcplwuZyxHkmKqVCpUKpVc5yhlPG4BMB9YBFwLPMDxu/YO4FXgU0Af8CSwmuNTPQOq1Wo14+Ul6fRUKpVghNmdNehLwHeBy9LtRSTz8hOBVcBCklfe/BL4R+CrDc5h0EvSCDUz6EeDQS9JI5Ql6H3DlCQFZ9BLUnAGvSQFZ9BLUnAGvSQFZ9BLUnAGvSQFZ9BLUnAGvSQFZ9BLUnAGvSQFZ9BLUnAGvSQFZ9BLUnAGvSQFZ9BLUnAGvSQFZ9BLUnAGvSQFZ9BLUnAGvSQFZ9BLUnAGvSQFZ9BLUnAGvSQFZ9BLUnAGvSQFZ9BLUnAGvSQFlzXoxwHfA7YAXcC0uvH5wLZ0fHHm6iRJuY3PeNzNwJnAdcBM4NvpPoA24CHgKqAPeBF4FtiTq1KpYA4ePMijjz5KT8/bXH/9bObOndvqkqSGsgb9LGBDur6VJNQHTAd2AfvS7W7gBuBvMl5LKpxDhw4xe/bneP31c+jru4oJE+7i619fyrJl97S6NOkEWaduJgG9NdtHas41ieMhD7Af6Mh4HamQ1q9fz09/epi+vueAr9HXt4nly5dz9OjRVpcmnSDrHX0v0F6zPQ4Y+B++r26sHdjb6CSdnZ3H1svlMuVyOWM5UnP19vZSrX6c4/c3H+PIkcP09/dz1llntbI0BVOpVKhUKrnOUcp43AKSB66LgGuBB4B56Vgb8G8kc/cHSB7Izgd66s5RrVarGS8vtdZbb73FpZdezYEDfw5cTVvbN7nmmt10d29sdWkKrlQqwQizO2vQl4DvApel24uAGcBEYBVwE/Agye3OauCRBucw6DWmdXd3s2jRPezZ08OsWbNZs+b7TJ48udVlKbhmBv1oMOglaYSyBL1vmJKk4Ax6SQrOoJek4Ax6SQrOoJek4Ax6SQrOoJek4Ax6SQrOoJek4Ax6SQrOoJek4Ax6SQrOoJek4Ax6SQrOoJek4Ax6SQrOoJek4Ax6SQrOoJek4Ax6SQrOoJek4Ax6SQrOoJek4Ax6SQrOoJek4Ax6SQrOoJek4Ax6SQrOoJek4MZnOOZsYA1wPrAfuA34n7qvWQnMSserwM1Ab/YyJUlZlTIcswyYCHwN+B3gN4A/rvuazcDngXeHOE+1Wq1muLwknb5KpRKMMLuzTN3MAjak6xuA32xwzouAVUA3sCjDNSRJo2S4qZs7OPFu/b85Pg2zH+ioG58APAw8lJ6/C9gO7MxVqSQpk+GCfnW61PpboD1dbwf+r268jyToD6bbm4DLaRD0nZ2dx9bL5TLlcvkkSpak00elUqFSqeQ6R9Y5+nbgq8DvAtcDX6kZnw6sBa4EzgAqwGLg9brzOEcvSSOUZY4+y6tuHgF+SPLA9ZfALen+pcAuYB3wGPAS0A/8gBNDXpLUJFnu6EeLd/SSNELNetWNJGkMMeglKTiDXpKCM+glKTiDXpKCM+glKTiDXpKCM+glKTiDXpKCM+glKTiDXpKCM+glKTiDXpKCM+glKTiDXpKCM+glKTiDXpKCM+glKTiDXpKCM+glKTiDXpKCM+glKTiDXpKCM+glKTiDXpKCM+glKTiDXpKCM+glKTiDXpKCyxP0XwCeGGTsTuBl4CVgXo5rSJJyGp/xuJXAZ4EdDcY+AvwRMAM4G+gG/gE4lPFaklR41WqVt99+m4kTJ9Le3t7qct4n6x39i8CXgVKDsWvS8X6gF9gFXJbxOpJUeD09PVxyyTVceOFlfOhDH2H58hWtLul9hgv6O4CddcsM4MkhjmkH9tVs7wc6ctQoSYW2cOGdvPHGXA4e3EN//5t85ztP8uyzz7a6rGOGC/rVwKfrln8Z5phekrAf0A7szVqgJBXdjh3bOXz4KySTHB/mwIEv8vLL21td1jFZ5+iHsg34BnAW8AFgOvBaoy/s7Ow8tl4ulymXy6egHEk6taZM+Ti9vZuAW4F+Jkx4galTvzQq565UKlQqlVznaDTHfrI+AywBbkm3l5LMx68DFgN3kfzG8A3g6QbHV6vVao7LS1Ix7Nixgxtv/G2q1cs5evTnXH31hTz//NOMHz/699KlUglGmN15gj4vg15SGO+88w5bt26lo6ODWbNmMW7cqXmbkkEvScFlCXrfGStJwRn0khScQS9JwRn0khScQS9JwRn0khScQS9JwRn0khScQS9JwRn0khScQS9JwRn0khScQS9JwRn0khScQS9JwRn0khScQS9JwRn0khScQS9JwRn0khScQS9JwRn0khScQS9JwRn0khScQS9JwRn0khScQS9JwRn0khRcnqD/AvDEIGMrge1AF7AJmJTjOpKkHLIG/Urgm0BpkPErgc8CNwJzgN6M12m5SqXS6hJOinWOnrFQI1jnaBsrdWaRNehfBL5M46AfB1wErAK6gUUZr1EIY+Uf3zpHz1ioEaxztI2VOrMYLujvAHbWLTOAJ4c4ZgLwMPB7wG8BfwB8OnelkqRMxg8zvjpdRqKPJOgPptubgMtJfkhIksaQMrC2wf7pwL+S/LbQRjLNM73B1+0Cqi4uLi4uI1p2MULD3dEPZeCiA5amBawDHgNeAvqBHwCvNzj+EzmuLUmSJEmSpFEzk+TNU/XmA9uALcDiplbU2GB1LgVeS8e6gIubWVSNNuBx4AVgK0n/ahWln8PVWZR+ngH8JcnLgDcDl9SNF6Wfw9VZlH4O+DDw8wZ1FKWfMHiNRerlKzV11L8opki9BOBe4FWSgmq1AW8AHen6NpLmt8pgdUISWlc0t5yGbgceStfPBf6zZqxI/bydweuE4vTz88BfpOufAf6uZqxI/RyqTihOPyHp1dPAT3h/SBapn4PVCMXp5QdIgr6REfeyGZ91swtYwIlvrpqeju0jeWjbDdzQhHoGM1idkLx34H6Su6n7mllUnR8DD6br44DDNWNF6udQdUJx+vkMsCRdnwrsrRkrUj+HqhOK00+AbwGPAD11+4vUz8FqhOL08nKS9yRtBP6JZLZhwIh72Yygf4oTv9Eh+fybfTXb+0l+QrXKYHVC8jLSJSQf5zAbmNesouocAH4BtJOE6Z/UjBWpn0PVCcXpJ8ARkleGPQz8qGZ/kfoJg9cJxenn7cA7wPPpdu1NU1H6eTuD1wjF6eUBkh9InwPuJvlcsYG8HnEvW/nplftIgmBAOyfeqRTFSuBdkp+e62ntr3YfI3kT2mPAX9XsL1o/B6sTitVPSL75Lyb52I6z031F6yc0rhOK089FwFySOeVfB37I8SmFovRzqBqhOL38Gcc/NPIN4H+BC9LtovTyBFNJXldfq43kL3MucCbJp11eQGtN5cQ6O0jmmM8h+en/Y5KPdmiFXyF5T8KNDcaK1M+h6ixSP28Flqfrk4D/IJkbhWL1c6g6i9TPWvUPMovUzwH1NRapl0uAP0vXf5Xk++mMdLuIvQSSAB14yLkQuDNdv4nkQcJ2kg9Ja7WpNK5zIUmdm4EVzS/rmJXAf3H8SXwXcAvF6+dwdRaln2cDfw38M8m/+3yK+f9zuDqL0s9aXcAnKWY/BzSqsSi9HM/xV669AFxLsXspSZIkSZIkSZIkSZIkSZIkSZIkSc3x/5i5VfYgedh3AAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Generate some data points\n", "# Let's use any simple function, say, f=sin(x) over interval [0, 2*pi]\n", "n = 3\n", "dmin = 0\n", "dmax = 2*pi\n", "ndata = 21\n", "\n", "x = np.linspace(dmin, dmax, n+2)[1:-1]\n", "xExtend = np.linspace(dmin, dmax, n+2)\n", "y = np.sin(x)\n", "\n", "plt.scatter(x,y)\n" ] }, { "cell_type": "code", "execution_count": 40, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 40, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD/CAYAAADoiI2GAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAE1VJREFUeJzt3X+MHGd9x/H3Gts0cewEOSTlR1ubo0kpbewoAZuzcx4rRq0KFiFA1VQCOzS0BdFEpVKSVkDO4o9WrYxASBBsoIa2UEpFqGibQIDbC/aZBIxIimhJY0TUqhWggGJ0BcmJt3/MPvGyufPs3s7MMz/eL2l1Mzd7e9/V+T43/uyzuyBJkiRJkiRJkiRJkiRJknRO24C5JT6/F3gAWABuKnUiSVKmW4GHSEN60BrgP4EL+9sPAJeUO5oktdeqEa7zCHA90Bn6/Iv6xx4HTgNHgZlcp5MkLWuUAP808MQSn99AGt7Bj0nPxiVJJRglwJfzOLB+YH898KPJxpEkjWr1BF/7H8AvA88CFknrk78avtLU1FTv5MmTE3wbSWqlk8ALz3WFcc7Ae/2PNwBvIu293wZ8jvQBzg8D//u0CU6epNfrNfZyxx13RJ/B+zf+5fOf77F2bY+pqWbev6b//Jp+33q9HsBUViiPegb+XWC6v/2Jgc//c/8i1Uq3C298I3zkI3D6NKxZE3siaXyTdOBSbXW78JrXwEUXwde/HnsaaWUM8AklSRJ7hEI18f4tLsI3vgEvexns3JnQ7caeqDhN/PkFTb5voxpe212EXr/PkSrh3nvhwAE4ehTuugsOHYK77449lfSzOp0OZGS0Z+BqnW4XwsnbzAwcO5b24FLdGOBqncEA37gRNm+2B1c9GeBqlcH+O0gSGt2Dq7kMcLXKwgJceSWsW3f2cwa46soAV6sM1ieBPbjqygBXqywV4PbgqisDXK2xVP8dJAnMLfWWJVKFGeBqjYUF2Lr1Z/vvIElgfr70kaSJGOBqjW4Xdu9e+pg9uOrIAFdrLNV/B/bgqiMDXK1wrv47cDmh6sYAVysstf57mAGuujHA1Qrnqk8Ce3DVjQGuVpifX/4BzMAeXHVjgKvxRum/A2sU1YkBrsY7fjztv88/P/u6BrjqxABX443Sfwf24KoTA1yNN06A24OrTgxwNdo4/XdgjaK6MMDVaOP034EvbKW6MMDVaOPUJ8HMTPrEH3twVZ0BrkbrdmHXrvG+xh5cdWGAq7FC/z09Pf7X2oOrDgxwNdZK+u/AAFcdGOBqrJX034HrwVUHBrgaa5IAtwdXHRjgaqSVrP8eZo2iqjPA1UjHj6fvf7mS/jswwFV1BrgaaZL6JLAHV9UZ4GqkPALcHlxVZ4CrcfLovwNrFFWZAa7GCf33ud7/clQGuKrMAFfj5FGfBPbgqjIDXI2zktc/WY49uKosK8BXAXcCC8AcMDV0/NXAV4EHgD/MfTppTJO8/slyrFFUVVkBfh2wFpgGbgcODh1/N/ByYAfwJ8CFeQ8ojSPP/jswwFVVWQG+A7inv30/cPXQ8dPARcB5QAfo5TqdNKY8++/AHlxVlRXgG4BTA/tPDn3NQeAE8E3gs0PXlUpXRIDbg6uqVmccPwWsH9hfBZzpb/8i8Fbgl4D/A/4WeC3wj8M3Mjs7+9R2kiQkef+GSeS7/ntYqFG2bcv/tiWAbrdLd8yurpNx/HpgL3AjsB14B/CK/rHLgH8AXkJapbyH9Ez8Q0O30ev1bFZUvC98AWZn4ejR/G/7rrvg0CG4++78b1taSqfTgYyMzjoDv4v0Qcpj/f0bgRuAC4DDwEdJV6j8FHgEOLLiaaUJFVGfBDMzsG9f2oOvWVPM95DGlXUGngfPwFWKnTvhwAG49tpibn/LlvQs3BpFZRjlDNwn8qgRiuy/A5cTqmoMcDXCJO9/OSoDXFVjgKsRiuy/A9eDq2oMcDVCGQHuenBVjQGu2iuj/w6sUVQlBrhqL4/3vxyVAa4qMcBVe90u7N5dzveyB1eVGOCqvTL678AeXFVigKvWyuy/A2sUVYUBrlorY/33MANcVWGAq9bKrE8Ce3BVhQGuWosR4PbgqgoDXLUVo/8OrFFUBQa4aitG/x0Y4KoCA1y1FaM+CezBVQUGuGorZoDbg6sKDHDVUsz+O7BGUWwGuGopZv8dGOCKzQBXLXW7sGtX3BlCD/7EE3HnUHsZ4KqlMl/Aajn24IrNAFftVKH/DqxRFJMBrtqpQv8dGOCKyQBX7cRcPjjMHlwxGeCqnSoF+MaNsGmTPbjiMMBVK1XqvwNrFMVigKtWqtR/Bwa4YjHAVStVqk8Ce3DFYoCrVqoY4PbgisUAV21Usf8OrFEUgwGu2qhi/x0Y4IrBAFdtVLE+CezBFYMBrtqocoCHHvzEidiTqE0McNVClfvvIElgfj72FGoTA1y1UOX+O7AHV9kMcNVCleuTwB5cZTPAVQt1CHDXg6tsWQG+CrgTWADmgKmh4y8B7gO+DPw9sDbvAaU69N+BNYrKlBXg15GG8jRwO3Bw4FgHOATsB64Bvghszn9EtV0d+u/AAFeZsgJ8B3BPf/t+4OqBY5cBjwFvA7rARcC3c55PqkV9EtiDq0xZAb4BODWw/+TA11xMemb+PmAPcC0Q+V0K1UR1CnB7cJUpK8BPAeuHrn+mv/0Y8AjpWfcTpGfqVyPlqE79d2CNorKszjh+DNgLfArYDjw0cOw7wAWkD2yeJO3BP7TUjczOzj61nSQJSV1OpxRdnfrvIEng8GG49dbYk6hOut0u3TH/8ndGOP5+4Ir+/o3AVaTBfZi0MvmL/vWOAX+8xG30er3eWENJwdvfDp0OvOtdsScZ3WOPwQtekH5cnXWKJC2j0+lARkZnBXgeDHCt2M6dcOAAXHtt7EnGs2VLehb+0pfGnkR1NUqA+0QeVVYd++/AHlxlMMBVWcePw9at9eq/gySBubnYU6jpDHBVVrcLu2u6MHVmBhYWXA+uYhngqqw6rf8e5npwlcEAVyXVuf8O7MFVNANclVTH9d/DDHAVzQBXJdW5Pgl8XRQVzQBXJTUhwO3BVTQDXJXThP47sEZRkQxwVU4T+u/AAFeRDHBVThPqk8AeXEUywFU5TQpwe3AVyQBXpTSp/w6sUVQUA1yV0qT+OzDAVRQDXJXS7cKuXbGnyJc9uIpigKtS6vwCVsuxB1dRDHBVRhP778AaRUUwwFUZTey/AwNcRTDAVRlNWj44zB5cRTDAVRlNDvDQg584EXsSNYkBrkpocv8dWKMobwa4KqHJ/XdggCtvBrgqocn1SRB68NOnY0+ipjDAVQltCPCNG2HzZteDKz8GuKJrQ/8dWKMoTwa4omtD/x0Y4MqTAa7o2lCfBPbgypMBrujaFOD24MqTAa6o2tR/B0kCc3Oxp1ATGOCK6vhx2Lq1Hf13kCQwPx97CjWBAa6o2lSfBPbgyosBrqjaGOD24MqLAa5o2th/By4nVB4McEUT+u9162JPUj4DXHkwwBVNG+uTwB5ceTDAFU2bA9weXHkwwBVFm/vvwBpFk8oK8FXAncACMAdMLXO9Q8Cf5ziXGq7N/XdggGtSWQF+HbAWmAZuBw4ucZ0/AH4N6OU7mpqszfVJYA+uSWUF+A7gnv72/cDVQ8engZcCHwQ6+Y6mJjPA7cE1uawA3wCcGth/cuBrngO8E3grhrfGYP99ljWKJpEV4KeA9UPXP9Pffi1wMfCvwG3A7wJvyHtANY/991m7dvnCVlq51RnHjwF7gU8B24GHBo69r38B2Af8CvCxpW5kdnb2qe0kSUja/n/nlrM+OWtmBvbvT3vwNWtiT6OYut0u3TH/O5ZVfXSA9wNX9PdvBK4CLgAOD1xvH3A58GdL3Eav1/PxTZ21cyfMzsKePbEnqYYrroDDh2HbttiTqEo6nQ5kZHQZ3bUBrqcsLsKll8L3vmeFEtx8MzzveXDbbbEnUZWMEuA+kUelsv9+Oh/I1EoZ4CqV/ffTuR5cK2WAq1QG+NNdfDFs2uR6cI3PAFdpXP+9PGsUrYQBrtIsLNh/L8cA10oY4CqN9cny7MG1Ega4SmOAL88eXCthgKsU9t/ZrFE0LgNcpVhYgCuvtP8+FwNc4zLAVQrrk2z24BqXAa5SdLvpK+9pefbgGpcBrsKF/nt6OvYk1WeNonEY4Cqc/ffoDHCNwwBX4ey/R2cPrnEY4CqcAT46e3CNwwBXoVz/PT5rFI3KAFeh7L/HZ4BrVAa4CmV9Mj57cI3KAFehDPDx2YNrVAa4CmP/vXLWKBqFAa7C2H+vnAGuURjgKoz1ycrZg2sUBrgKY4CvnD24RmGAqxCLi/Dgg/bfk0gSmJuLPYWqzABXIXz/y8klCczPx55CVWaAqxDWJ5OzB1cWA1yFMMAnZw+uLAa4cuf67/zs3u1yQi3PAFfuXP+dH9eD61wMcOXO+iQ/9uA6FwNcuTPA87NxI2zebA+upRngypX9d/6sUbQcA1y5sv/OnwGu5RjgypX1Sf7swbUcA1y5MsDzZw+u5Rjgyo39d3GsUbQUA1y58fVPiuMLW2kpWQG+CrgTWADmgKmh4zcAXwGOAh8AOnkPqProdtNnDip/MzPpH0h7cA3KCvDrgLXANHA7cHDg2HnAu4AE2AlcCLwy/xFVF/bfxbEH11KyAnwHcE9/+37g6oFjPwVe1v8IsBr4Sa7TqTbsv4tnD65hWQG+ATg1sP/kwNf0gB/0t/8IWAd8IdfpVBuu/y6eAa5hqzOOnwLWD+yvAs4M7f8l8ELgNcvdyOzs7FPbSZKQ+P/sxrE+Kd7MDOzbl/bga9bEnkZ563a7dMf8C531oOP1wF7gRmA78A7gFQPHD5NWKDeTnpEvpdfrLXdITbFjBxw4AHv2xJ6k2bZsgUOHYNu22JOoaJ1OBzIyOivAO8D7gSv6+zcCVwEXAF/rX+4buP57gc8M3YYB3nCLi3DJJfD971uhFO2WW+C5z4Xbbos9iYo2SoBnVSg94M1Dn3t4YPsZ44+lprH/Lk+SpGfgBrjAJ/IoB/bf5fF1UTTIANfEDPDyuB5cgwxwTcT13+VzOaECA1wTsf8unwGuwADXRKxPymcPrsAA10S6Xdi1K/YU7WIPrsAA14qF/nt6OvYk7WONIjDANQH773gMcIEBrgnYf8djDy4wwDUBAzwee3CBAa4Vcv13fNYoMsC1Ivbf8RngMsC1Ir7/ZXz24DLAtSL23/HZg8sA19gWF+HBB+2/q8Aapd0McI0t9N/nnx97Ehng7WaAa2zWJ9VhD95uBrjGZoBXhz14uxngGov9d/VYo7SXAa6xHD9u/101SQLz87GnUAwGuMZifVI9oQd/4onYk6hsBrjGYoBXz8aNsGmTPXgbGeAama9/Ul324O1kgGtk9t/VZYC3kwGukVmfVJc9eDsZ4BqZAV5d9uDtZIBrJPbf1WeN0j4GuEZi/119Bnj7GOAaifVJ9dmDt48BrpEY4NVnD94+Brgy2X/XhzVKuxjgynTsGGzdav9dB0kCX/pS7ClUFgNcmT75SXjVq2JPoVHs2QNf+Qr84AexJ1EZOiV8j16v1yvh26gIi4vw/OfDt74Fz3lO7Gk0ije8Aa66Cm65JfYkmkSn04GMjPYMXOf06U/Djh2Gd53s3w9HjsSeQmUwwHVOR46kgaD6SBL44Q/TB57VbFkBvgq4E1gA5oCpoeN7gQf6x2/KfTpF9eij6bvv7N0bexKNY9Uq2LcPPvrR2JOoaFkBfh2wFpgGbgcODhxbA7wbeDmwC/h94JICZqy0boPXbH3sY3DNNV2e+czYkxSnqT+/ffvg4x+He+/txh6lME392Y0jK8B3APf0t+8Hrh449iLgEeBx4DRwFJjJe8Cqa+o/ojNn0vrk2c/uxh6lUE39+U1NweWXw4c/3I09SmGa+rMbR1aAbwBODew/OfA1G0jDO/gxcGF+oymmo0fhvPN88LLO9u+3B2+61RnHTwHrB/ZXAWf6248PHVsP/GipG2lyh/rtb8OJE7GnyN/DD8NNN6XLCFVPr3sdvOUtzf39a+rvXp6uB/66v70d+JeBY2uAh4FnkfbkXwOWOl97BOh58eLFi5exLo8woQ7wAeBY/3IZcAPwpv7xV5KuQvka8OZJv5kkSZIkSVJ7rQP+CZgH7gWeG3ec3F0IfBbokj6haXvUaYrzauDvYg+Rk6wnqDXFNtL71zRrgL8B7iNd3ty0h2mfAXyEdGn2l4EXxxzmFuDt/e19wHsizlKEWeDm/vZlQBMfF38v8O/Ax2MPkpPrSX9BIA25z0ScpSi3Ag+R/pFqmv2kTyKEdBHFo/FGKcSrgA/1t3dRgX+fYe34O4E7Yg5SgAuBn+tvv5j0r2bT/DaQAJ+IPEdeDpLep+C/Yw1SoOuBFwLHYw9SgHXABf3tjcDJiLMU5Rn9j/s4uxKwcL8H/NvQ5ar+sS8C3wO2lDVMAc51/34e+DpwTZzRcnGu+5fQnAA/DPzmwP6jNPOF3TbRzAAP1gNfAn4n9iAFOUL6fJuXR57jKZeTw9rGCvp14JvAb8QepEAJzQnwg8DrBvb/K9YgBdtEcwP8F4CvktYpTXYp8F3gvKUOlnHW8afA6/vbi0DT3jP7V4FPka6P/1zkWTSaY8Bv9be3k3bFqo9Lgc+T9vxH4o5SiNeT5ibAT0if/X5m+asX6xLgbtJHw+8DmvbWuJ8BvkN6/+aAu+KOU5hdNOdBzKWeoNZEm2jmg5jvBf6Hs79zc5x9HKoJzgM+Sbpyb4HmrbKRJEmSJEmSJEmSJEmSJEmSJEmS6uX/ASXBCi3o/H8zAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Define hat function (basis function of our choice)\n", "\n", "# Recall: the mother of all hat functions: h(x) = max{1-|x|, 0} \n", "\n", "def motherhat(x): \n", " return max(1-abs(x), 0);\n", " \n", "# TEST\n", "xx = np.linspace(-3, 3, 101)\n", "yy = list(map(lambda t: motherhat(t), xx))\n", "plt.plot(xx, yy)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Step 1. Solve for n=3 (3 hierarchical basis functions)" ] }, { "cell_type": "code", "execution_count": 43, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD/CAYAAADoiI2GAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAEgVJREFUeJzt3X+snfVdwPH3aUvDpC1bRLrRTm7XiUwZYii2tsLOgjAy1gRxmvAXrW4qcU7HHwZHgGsIDDWdIyuDDaZuGlnECHHg0CA9Iq2DQeYYcQNaetoaZ2iXQbc6OkaPfzznoWeHe+9zfzw/vt/v834lJ33Oj3vu94x73zv3ez7nXpAkSZIkSZIkSZIkSZIkSZrRemDHFJdvBh4HdgEfqHVFkqRCfwg8RRbpUScAzwEnD48fB06td2mS1F6LZnGb3cDlQGfs8ncMr3sJeAV4FLig1NVJkqY1m4D/A/DDKS5fQRbv3HfJno1Lkmowm4BP5yVg+cj55cB3FrYcSdJsLVnAx34T+CngTcARsu2TPxu/0dq1awd79uxZwKeRpFbaA7x9phvM5Rn4YPjvFcAHyfa9rwb+mewFzs8C33rdCvbsYTAYJHu64YYbGl+Dj2/up5dfHgADLrwwzceX+n+/1B/bYDAAWFsU5dk+A+8DG4fHd49cfv/wJEVl//7s3xdfbHYd0kIsZA9cila/D4sXG3DFzYAvULfbbXoJlUr18fX7cO658IMfdJteSqVS/e8HaT+22Rqf7a7CYLifIwXj2mvhlVdg+3Y4cgQ6dXwnSHPQyb4oZ/zK9Bm4Wqnfh7POgpNOghdeaHo10vwYcLVSvw8TE9mp3292LdJ8GXC1kgFXCgy4WufoUTh0CFatMuCKmwFX6+zfD6tXZ2OEBlwxM+BqnXz7BAy44mbA1ToGXKkw4Gqd0YCffjrs2we+VUExMuBqndGAL1vmLLjiZcDVOqMBB7dRFC8DrtYx4EqFAVer5DPgp512/DIDrlgZcLXK6Ax4zoArVgZcrTK+fQIGXPEy4GoVA66UGHC1Sr+fzX6PchZcsTLgapV+H9as+dHLnAVXrAy4WmWqLRRwG0VxMuBqFQOulBhwtcZUM+A5A64YGXC1xlQz4DkDrhgZcLXGdNsnYMAVJwOu1jDgSo0BV2vMFHBnwRUjA67WmCngzoIrRgZcrTFTwMFtFMXHgKs1DLhSY8DVCjPNgOcMuGJjwNUKM82A5wy4YmPA1QpF2ydgwBUfA65WMOBKkQFXK8wm4M6CKzYGXK0wm4A7C67YGHC1wmwCDm6jKC5FAV8E3AHsAnYAa8eu/xXgK8DjwO+UvjqpJAZcKVpScP1lwFJgI7Ae2Da8LPdx4OeBI8B/AXcDL5W/TGn+ZjMDnjPgiknRM/BNwIPD48eAdWPXvwK8EXgD0AF8+UfBmc0MeM6AKyZFAV8BHB45/+rYx2wDngSeBr44dlspCLPdPgEDrrgUbaEcBpaPnF8EHBse/yTwIeB04P+AvwHeD/z9+J1MTk6+dtztdul2u/NdrzRn/X42IjgbExOwd2+Vq5Gm1uv16PV6c/qYTsH1lwObga3ABuA64NLhdWcAfwecR7aV8gmyZ+J3jd3HYOBgrRp07bVw4olw3XXFt/3e9+DUU+HIEegUfXdIFepkX4AzfhUWPQO/F7gI2Dk8vxW4AlgG3Al8jmxC5WVgN/BX816tVJF+Hy65ZHa3HZ0FX7my0mVJC1YU8AFw1dhlz44c//nwJAVrLnvgcHwf3IArdL6RR8mbb8Cl0BlwJW0uM+A5A65YGHAlbS4z4DkDrlgYcCVtrtsnYMAVDwOupBlwpcyAK2nzCbi/F1yxMOBK2nwC7u8FVywMuJI2n4CD2yiKgwFX0gy4UmbAlaz5zIDnDLhiYMCVrPnMgOcMuGJgwJWs+W6fgAFXHAy4kmXAlToDrmQtJODOgisGBlzJWkjAnQVXDAy4krWQgIPbKAqfAVeyDLhSZ8CVpIXMgOcMuEJnwJWk/fth1ar5zYDn/Av1Cp0BV5L6fVizZmH3MTGRTaJIoTLgStJC97/BLRSFz4ArSWUE3Flwhc6AK0llBNxZcIXOgCtJZQQc3EZR2Ay4kmTA1QYGXMkpYwY8Z8AVMgOu5Czk94CPM+AKmQFXcsraPgEDrrAZcCXHgKstDLiS0+9nM9xlcBZcITPgSs6+feU9A3cWXCEz4EpOmVso4DaKwmXAlRwDrrYw4ErK0aNw8GA5M+A5A65QGXAl5cCB7PeAL1lS3n0acIXKgCspZW+fgAFXuIoCvgi4A9gF7ADWjl1/HvAI8O/AF4ClZS9QmgsDrjYpCvhlZFHeCFwDbBu5rgN8BtgCnA/8K7DAv4EiLUwVAXcWXKEqCvgm4MHh8WPAupHrzgC+DVwN9IA3As+UvD5pTqoIuLPgClVRwFcAh0fOvzryMaeQPTP/JPDLwIXAu8teoDQXVQQc3EZRmIoCfhhYPnb7Y8PjbwO7yZ51/5Dsmfo6pAaV+Tb6Uf6FeoWoaNhqJ7AZuAfYADw1ct3zwDKyFzb3kO2D3zXVnUxOTr523O126Xa7812vNK18BnzVqvLv279Qr6r1ej16vd6cPqYzi+s/BZw9PL8VOJcs3HeSbZncMrzdTuAjU9zHYOCrP6rB7t1w8cXw/PPl3/dtt8HTT8Ptt5d/39JUOp0OFDS66Bn4ALhq7LJnR453AOvnvDKpAlXtf0N2v/ffX819S/PlG3mUjL17qw24L2IqNAZcyej3YU1F70RwFlwhMuBKRpVbKM6CK0QGXMmoMuDgNorCY8CVDAOutjHgSsLRo3DoULm/B3ycAVdoDLiSsH8/rF4NixdX9zkMuEJjwJWEqrdPwIArPAZcSTDgaiMDriTUEXBnwRUaA64k1BFwZ8EVGgOuJNQRcHAbRWEx4EqCAVcbGXBFr44Z8JwBV0gMuKJXxwx4zoArJAZc0atr+wQMuMJiwBU9A662MuCKXp0BdxZcITHgil6dAXcWXCEx4Ipev589M67LxET259ukphlwRa/KP6U2lYmJbBtFapoBV9TqnAHP+UKmQmHAFbU6Z8BzBlyhMOCKWp0vYOYMuEJhwBU1A642M+CKWhMBdxZcoTDgiloTAXcWXKEw4IpaEwEHt1EUBgOuqBlwtZkBV7SamAHPGXCFwIArWk3MgOcMuEJgwBWtprZPwIArDAZc0TLgajsDrmg1GXBnwRUCA65oNRlwZ8EVAgOuaDUZcHAbRc0z4IqWAVfbFQV8EXAHsAvYAayd5nafAT5W4rqkGTU5A54z4GpaUcAvA5YCG4FrgG1T3Oa3gbMAX85RbZqcAc8ZcDWtKOCbgAeHx48B68au3wj8AvBpoFPu0qTpNb19AgZczSsK+Arg8Mj5V0c+5i3A9cCHMN6qmQGXYEnB9YeB5SPnFwHHhsfvB04B/gl4M/BjwDeAz5e8Rul16v5L9FM5/fRsHYMBdHwKowYUBXwnsBm4B9gAPDVy3SeHJ4ArgTOZJt6Tk5OvHXe7Xbrd7rwWK+X6fbjkkmbXsGxZdnrhBVi5stm1KH69Xo9erzenjyl63tABPgWcPTy/FTgXWAbcOXK7K4GfBj46xX0MBr5dTSXbtAluuQXOP7/ZdZx3HmzfDuvXN7sOpaeT/Vg3Y6OLnoEPgKvGLnt2itt9bvbLkhYuhD1wOL4PbsDVBN/Io+iEMAOe84VMNcmAKzohzIDnDLiaZMAVnVC2T8CAq1kGXNEx4FLGgCs6IQXc3wuuJhlwRSekgPt7wdUkA67ohBRwcBtFzTHgio4BlzIGXFEJaQY8Z8DVFAOuqIQ0A54z4GqKAVdUQts+AQOu5hhwRcWAS8cZcEUlxIA7C66mGHBFJcSAOwuuphhwRSXEgIPbKGqGAVdUDLh0nAFXNEKcAc8ZcDXBgCsaIc6A5wy4mmDAFY0Q/hL9dCYmYO/eplehtjHgika/D2vWNL2KqU1MZKOEUp0MuKIR6guY4Cy4mmHAFY2QA+4suJpgwBWNkAMOvpCp+hlwRcOASz/KgCsKIc+A5wy46mbAFYWQZ8BzBlx1M+CKQujbJ2DAVT8DrigYcOn1DLiiEEPAnQVX3Qy4ohBDwJ0FV90MuKKwb1+4vwdllNsoqpMBVxRieAYOBlz1MuAK3tGjcPBg2DPgOQOuOhlwBe/AAVi1CpYsaXolxQy46mTAFbxYtk/AgKteBlzBM+DS1Ay4ghdTwJ0FV52KAr4IuAPYBewA1o5dfwXwZeBR4HagU/YCpZgC7iy46lQU8MuApcBG4Bpg28h1bwBuBLrALwEnA+8rf4lqu5gCDm6jqD5FAd8EPDg8fgxYN3Ldy8AvDv8FWAJ8v9TVSRhwaTpFAV8BHB45/+rIxwyAg8Pj3wNOAh4qdXVqvZhmwHMGXHUpmqw9DCwfOb8IODZ2/k+BtwO/Ot2dTE5Ovnbc7XbpdrtzXKba6sCBLN4xzIDnJibg619vehWKTa/Xo9frzeljil50vBzYDGwFNgDXAZeOXH8n2RbKh8mekU9lMPAlec3TQw/BzTfDww83vZLZe+AB2L4dvvSlpleimHU6HShodNHzmnuBi4Cdw/NbySZPlgFPAL8BPALk3163AvfNb7nS68W2/w1uoag+RQEfAFeNXfbsyHHAf+BKKdi7N76Aj86CdxysVYV8I4+C1u/DmjVNr2JunAVXXQy4ghbjFgq4jaJ6GHAFzYBL0zPgCtbRo3DoUFwz4DkDrjoYcAVr/35YvRoWR/hSuQFXHQy4ghXr9gkYcNXDgCtYBlyamQFXsGIOuL8XXHUw4ApWzAF3Flx1MOAKVswBB7dRVD0DrmAZcGlmBlxBinkGPGfAVTUDriDFPAOeM+CqmgFXkGLfPgEDruoZcAXJgEvFDLiClELAnQVX1Qy4gpRCwJ0FV9UMuIKUQsDBbRRVy4ArSP1+tgURu4mJ7M/CSVUw4ApOPgO+alXTK1m4iYlsH1yqggFXcFKYAc+5haIqGXAFJ5X9bzDgqpYBV3AMuDQ7BlzBSSngzoKrSgZcwUkp4M6Cq0oGXMFJKeDgNoqqY8AVHAMuzY4BV1BS+D3g4wy4qmLAFZSUZsBzBlxVMeAKSmrbJ2DAVR0DrqAYcGn2DLiCkmLAnQVXVQy4gpJiwJ0FV1UMuIKSYsDBbRRVw4ArKAZcmj0DrmCkOAOeM+CqQlHAFwF3ALuAHcDases3A48Pr/9A6atTq6Q4A54z4KpCUcAvA5YCG4FrgG0j150AfBy4CHgX8FvAqRWsMWi9Xq/pJVSqzsfXxPZJXY+vqYCn/PWZ8mObraKAbwIeHB4/Bqwbue4dwG7gJeAV4FHggrIXGLrUv4gMeDkMePlSfmyzVRTwFcDhkfOvjnzMCrJ4574LnFze0tQ2qb6ACc6CqxpLCq4/DCwfOb8IODY8fmnsuuXAd8pb2tx99atw/fX1fs5nnoEnn6z3c9apzsf3ta/BTTfV87nqtmxZdnrve2FJ0XddiVL++qz7sZ1zDtx4Y32frwyXA385PN4APDBy3QnAs8CbyPbJnwDeMsV97AYGnjx58uRpTqfdLFAHuB3YOTydAVwBfHB4/fvIplCeAK5a6CeTJEmSJElqr8XArWRjho8DlzS7nMqcCbxI9npASk4Gvgj0yN6wtaHR1ZSj6A1qsTsB+GvgEbLx383NLqcypwIHyLZ2U/NHZF+fXwGubHIhW4DbhsenAR9pbimVWUH2Au//kl7AJ4EPD4/PAFKYabgc+Ivh8XrgvgbXUoUtZG+yg2zIYF9zS6nMCcC9wDdJL+Bd4B+HxycBf9zcUuBvgY8C95NF7m1NLqYCHeBu4GxgL+kF/GTgxOHxz5L9JBW7bcCvj5z/76YWUpGTgGXD4x8H9jS4lqp8AriY7Ceo1AJ+M3AT2ROLh4Fzp7th2ROpvwn8wdhlB4Hvk02sXEA2lviukj9vXaZ6fPuALwBPDc93al1RuaZ6fFvInnW/mezH8t+veU1VmO4Nasemvnl0jgz/XQ7cA1zb4FqqsIWsK/9CttUQ8/fcVH4CeCtZM99G9mz8zKYWczfZj6y5bzW1kIo8R/YsYAfZ/1H1Gl1NNd4JPA28p+mFlGQb8Gsj5w80tZAKvZVs/3RLw+uowr+RfZ/tIHvz4JeBlU0uqGQfA64eOf+fwCkNrYXfBe4aHv8c2f/YqUpxC+VnyPYZ39n0Qko00xvUUrAS+Abw7qYXUoMUt1AuJfvpArLXDZ+jwZ8ylgKfBf5jeDqnqYXU4HnSC/h9ZI8r/ynj3maXU4qp3qCWkluB/+H4f7MdHH8dIzUpBhzgTzj+JsmLGl6LJEmSJEmSJEmSJEmSJEmSJEmS1G7/D4sxqimnJbuiAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Now, let's define any hat function : hat_n,i(x) = motherhat( (x-x_n,i)/h )\n", "\n", "def hat(a,b,x):\n", " h = (b - a)/2\n", " m = (a+b)/2\n", " if (h > 0):\n", " return motherhat((x-m)/h)\n", " else:\n", " return 0\n", " \n", "# test\n", "xx = np.linspace(-5, 5, 101)\n", "yy = list(map(lambda t: hat(-1,1,t), xx))\n", "plt.plot(xx, yy)\n", "\n", "\n", "dmid = (dmin+dmax)/2\n", "\n", "# 1. construct basis functions\n", "h0 = lambda t: hat(dmin, dmid, t)\n", "h1 = lambda t: hat(dmin, dmax, t)\n", "h2 = lambda t: hat(dmid, dmax, t)\n", "\n", "basis = [h0, h1, h2]" ] }, { "cell_type": "code", "execution_count": 44, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[1. 0.5 0. ]\n", " [0. 1. 0. ]\n", " [0. 0.5 1. ]]\n" ] } ], "source": [ "# 2. construct System matrix A:\n", "# | h0(x0) h1(x0) h2(x0)|\n", "# | h0(x1) h1(x1) h2(x1)|\n", "# | ... |\n", "# | h0(xm) h1(xm) h2(xm)|\n", "\n", "A = np.array([[h(i) for h in basis] for i in x])\n", "\n", "print(A)\n" ] }, { "cell_type": "code", "execution_count": 45, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[ 1.00000000e+00 1.11022302e-16 -1.00000000e+00]\n" ] } ], "source": [ "c = np.linalg.solve(A, rhs)\n", "\n", "print(c)\n" ] }, { "cell_type": "code", "execution_count": 46, "metadata": {}, "outputs": [], "source": [ "# 4. Define the resulting approx. function with 3 function supports (basis functions)\n", "\n", "def func1(x):\n", " result = 0.0\n", " result += c[0] * h0(x)\n", " result += c[1] * h1(x)\n", " result += c[2] * h2(x)\n", " return result\n", " " ] }, { "cell_type": "code", "execution_count": 47, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 47, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD/CAYAAAD/qh1PAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAGblJREFUeJzt3XuUVNWZ9/FvNc1NoDsg4z0TvHIHBREExQKJxgsqZqKiMkoU0XFMRhEVs5RmHDUEcYKOmnmNb5JJ1GgSLwuNokZKURRQMIBAOzjeJROXKPcWpGv+2IV22m5aqrr7nDr9/axVy6o+VNUvO/r04Tn77A2SJEmSJEmSJEmSJEmSJElqRIOBuXX8/HJgee7YXOCQ5gwlSWocVwFLgfl1HPs1cFjzxpEk1aWkgPeuBk4HUnUcGwhcC8wDringOyRJEesGvFTHz68DugCtgceAk5oxkySphkLO6HdmFrAW2AY8jm0cSYpMaRN8Zjmhd98L2AyMBO6p/YcOPPDA7JtvvtkEXy9JifYmcNCuvKExzuizuX+OBSYA6wh9+bnA84TZN0/WftObb75JNpuN/WPq1KmRZzCnGc1pzh0P4MBdLdKFntG/DQzNPb+/xs/vr/VakhSRpurRS5JiwkLfgHQ6HXWEr8WcjacYMoI5G1ux5MxHXXPgm0s212+SJH1NqVQKdrF2e0YvSQlnoZekhLPQS1LCWeglKeEs9JKUcBZ6SUo4C70kJZyFXpISzkIvSQlnoZekhLPQS1LCWeglKeEs9JKUcBZ6SUo4C70kJZyFXpISzkIvSQlnoZekhLPQS1LCWeglKeEs9JKUcBZ6SUo4C70kJZyFXpISzkIvSQlnoZekhCu00A8G5tbx89HAQmA+cGGB3yFJKkCqgPdeBZwLbASG1vh5a2AFcDiwGXgROBn4a633Z7PZbAFfL0ktTyqVgl2s3YWc0a8GTq/jC3vmjq0DtgEvAMML+J7Y2Lp1Kxdd9AO6dNmPffftzq9/fW/UkSSpQaUFvPchoFsdPy8jFPkdNgDlBXxPbFxxxbX85jdvsGXLPD75ZA0XX3wG++67NyNHjow6miTVqykuxq4DOtV43Qn4pAm+p9k99NBstmy5BdgfGMrmzZfx6KN/jDqWJO1UIWf09VkFHAx0BjYR2jYz6vqDFRUVXzxPp9Ok0+kmiNN4ysrKWbPmLaBP7ifnsnZtJsJEkpIuk8mQyWQK+oxCLsZCaN3cR7gYOxboCNxNuPh6PeFvDPcAd9Xx3qK7GDtnzhxOP30cVVXfp7R0De3bt6VLlzvp0aOUW26BXr2iTigp6fK5GFtooS9E0RV6gMWLFzN79mN07NiB8847j7KyrtxxB9x0E5xxBkybBl27Rp1SUlJZ6CP08cdQUQG//S1ccw1cdhm0aRN1KklJ09zTK1XD7rvD7bfDvHkwd25o4zz8MCTod5mkIuUZfRN5+mm44orwC+Df/x0OOyzqRJKSwDP6GPn2t2HJEhg7Fk44AcaPhw8/jDqVpJbIQt+ESkth4kSorIQ99oC+feGGG2Dz5qiTSWpJLPTNoLwcpk+HV16BZcugRw+4916oro46maSWwB59BF54AS6/HEpKQv9+6NCG3yNJ4PTKolJdHc7qr702FPrp06Fbt6hTSYo7L8YWkZISGDcOVq0KUzEHDoQpU2D9+qiTSUoaC33EOnSAqVNh6VJYswa6d4e774bt26NOJikpbN3EzKuvhv79unVw661w7LFRJ5IUJ/boEyKbhYcegquugt69YcaMcKYvSfboEyKVgu9+F1asgKOPhmHD4Ic/hLVro04mqRhZ6GOsbVuYPDkU/K1bw/z7WbNg27aok0kqJhb6IrDHHnDXXWGxtCeegD59YPZsF0yT9PXYoy9CTzwBkybBPvuEC7b9+kWdSFJzsUffQpxwAvz5zzBmTFg8bcIE+Mtfok4lKa4s9EWqdWu49NJww1VZWWjn3HwzVFVFnUxS3Fjoi1znzjBzJrz8MixaFC7YPvCA/XtJX7JHnzCZTNjwpH37sGDaEUdEnUhSY7JHL9LpcGZ/4YWhh3/OOfDee1GnkhQlC30CtWoVdrSqrIQDDoBDD4XrroONG6NOJikKFvoE69gx7Gi1ZAm89VZYRuEXv3DDE6mlsUffgixYEBZMq6oK/ftjjok6kaRd5aJmalA2Cw8+CFdfDQMGwE9+AgcdFHUqSV+XF2PVoFQKzjwTVq6EQYNgyJBwl+2nn0adTFJTsdC3UO3bhx2tli8Pu1p17w533AGffx51MkmNzdaNgLCkwhVXhF2uZs4MyyxIih979CpINguPPQZXXgn77x8Kfu/eUaeSVJM9ehUklYLRo2HZsnBGP2IEXHIJfPRR1MkkFSLfQl8C/AyYD8wFDqx1/HJgee7YXOCQfAOq+bVpE3a0WrUqPO/ZM2xn+NlnUSeTlI98C/1pQBtgKHANMLPW8QHAOGBE7vFGvgEVnS5dwo5WL74I8+ZBr17whz+4YJpUbPLt0c8EFgAP5l6/D+xX4/gK4HVgL+Bx4Md1fIY9+iLzzDPhgm3nzmHDk4EDo04ktTzN2aMvA9bXeL291mfdD0wERgJHASfl+T2KkVGjwnIK554LJ58M550HH3wQdSpJDSnN833rgU41XpcANVdQmcWXvwgeBw7L/fNvVFRUfPE8nU6TTqfzjKPm0qpV2NHqzDPDRif9+sEPfhBm6nToEHU6KXkymQyZTKagz8i3dXM6MBoYDwwBruPLs/ZyYCnQC9hMaO/cAzxZ6zNs3STA22+H5RTmz4ebbgrLIpc4l0tqMs05jz4F3Ans2JZ6PDAQ6AjcDYwlzLz5DHgGmFbHZ1joE2T+/LBgWnV1WDDtqKOiTiQlkzdMKVLV1XD//WFphSFDYPr0cOOVpMbjDVOKVElJaN2sWgV9+8Lhh4e2zvr1Db9XUtOx0KvR7bZb2NFq2TL461/hkEPgP//TBdOkqNi6UZNbvDj079euDevnHHdc1Imk4mWPXrGVzcIjj8DkyWFJ5JkzoUePqFNJxccevWIrlYIxY+D112HkSDj6aLjsMvj446iTSclnoVezats27Gi1YkWYpdOjR5iOuXVr1Mmk5LLQKxJ/93dhR6vnnoOnn4Y+feDRR10wTWoK9ugVC3PmhAXT9twzLJh26KFRJ5LiyR69itbxx4ftDL/3vfD8ggvCtoaSCmehV2yUloYdrSorw1r4ffrAjTfCli1RJ5OKm4VesfONb4QdrRYuDHPwe/QISyvY6ZPyY49esff88+GGqzZtwgydIUOiTiRFxx69Emn4cFi0CC6+GP7hH+Dss+Gdd6JOJRUPC72KQklJ2NGqshIOPhgGDIAf/Qg2bIg6mRR/FnoVlQ4dYNo0eO01ePfdsJzCPffA9u1RJ5Piyx69itrChWH+/caNYf79yJFRJ5KalouaqUXKZuF3vwtr3/frB7fcEto7UhJ5MVYtUioFZ5wBK1fC0KFw5JFhls4nn0SdTIoHC70So127cFb/+uuweXPo399+O2zbFnUyKVq2bpRYy5aF/v3774d2zoknhrN/qZjZo5dqyWbhj38MSyP//d+HDU/69o06lZQ/e/RSLakUnHRSOLsfPRqOPRYmTgx72UothYVeLULr1mFHq1WrwublvXrB9OlQVRV1MqnpWejVonTpEtbLmT8/PHr1ClMz7SIqyezRq0V79tkwFbOsLNxwNWhQ1ImknbNHL+2ikSPDUsjnnw+nngr/+I9hlo6UJBZ6tXitWoUdrSor4ZvfhP79YepU2LQp6mRS47DQSzmdOoUdrRYvhjfeCDdc/epXUF0ddTKpMPbopXq89FLo32/bFi7gDh8edSKpeXv0JcDPgPnAXODAWsdHAwtzxy/M8zukSB15ZCj2V14J48bBd78Lb74ZdSpp1+Vb6E8D2gBDgWuAmTWOtQZuBb4NHANcBOxRQEYpMqkUjB0b5t8PGABHHAGTJ8O6dVBVVcVdd93F9ddP5emnn446qlSvfAv9MODJ3PMFwOE1jvUEVgPrgG3AC4B/6VVRa98+7Gi1fDmsXQvdu2fp3v1OJk16nBtuyHLaaRdx6623RR1TqlO+hb4MWF/j9fYan1VGKPI7bADK8/weKVb23jvsaDV58lw++CDNli2zgX9l8+ZnmTJlCtVeuVUMleb5vvVApxqvS4Ad/4avq3WsE1DnyuAVFRVfPE+n06TT6TzjSM2ra9f3aNv252zefF/uJ9/i888PYtu2bbRt2zbSbEqWTCZDJpMp6DPynXVzOuGC63hgCHAdcFLuWGvgdWAwsIlwQXY0sKbWZzjrRkXr7bffpk+fQWza9P+AQZSUvExJySguuugbTJsGXbtGnVBJ1Zyzbh4GqoAXCRdiLwfGAhMIffkrgDmEIn8PXy3yUlHr1q0bTz75MAcddANlZYM4/vgHWLEiS0kJ9OwZlkP+7LOoU0qB8+ilRrZyZZiSWVkJM2bAaae54YkajxuPSDHy1FNhh6uuXcOCaQMGRJ1ISeCiZlKMHHccvPZamId/4okwfjx8+GHUqdQSWeilJlRaGna0qqyEPfYI2xjecEPYvFxqLhZ6qRmUl4cdrRYtgqVLoUcPuPdeF0xT87BHL0XghRfCgmklJWHBtKFDo06kYuHFWKmIVFeHs/prrw2Ffvp06NYt6lSKOy/GSkWkpCSsirlqVdi7duBAmDIF1q9v+L3SrrDQSxHr0CHsaLV0aZiV07073H03bN8edTIlha0bKWZeeSXMv//00zD/ftSoqBMpTuzRSwmRzcJDD8FVV4W2zi23hDN9yR69lBCpVNjRasWKsIXhsGHwwx+GtfClXWWhl2Ksbduwo9WKFbB1a5h/P2tW2MdW+rps3UhFZPlymDQJ3n47tHNOPtkF01oae/RSC5DNwpNPhoK/zz5hSeT+/aNOpeZij15qAVIpOOEE+POfYcyYsHjahAnwl79EnUxxZaGXilTr1nDppeGGq7Iy6NMHbr4ZqqqiTqa4sdBLRa5z59C+efnlsGhajx7wwAOhxSOBPXopcTKZcMNVu3ZhwbTBg6NOpMZkj14S6XQ4s58wAU4/Hc45B959N+pUipKFXkqgVq3CjlaVlXDAAXDYYXDddbBxY9TJFAULvZRgHTuGHa2WLIG33grLKPziF2540tLYo5dakAULwoYnVVVhwbR0OupE2lXeMCWpQdksPPggXH11aOnMmAEHHRR1Kn1dXoyV1KBUCs48M8y/P+IIGDIk3GX76adRJ1NTsdBLLVS7dmFHq+XLw65W3bvDHXfA559HnUyNzdaNJCAsqTBpUtjlaubMsMyC4scevaSCZLPw2GNw5ZWw//5hhcw+faJOpZrs0UsqSCoFo0fDsmXhjH7kSLjkEvjoo6iTqRAWeklf0aZN2NFq1aqw+UnPnmF2zmefRZ1M+cin0LcH/gA8DzwOdK3jz8wCXgHmAs8CZfkGlBSdLl3gpz+FF1+EefNCwf/9710wrdjk06O/AugI/CtwJnAk8C+1/sw84FRgZztc2qOXiswzz4QLtuXlYcG0gQOjTtTyNFePfhjwZO75k8CoOj7zYOBu4AVgfB7fISmGRo2CxYth3LiwjeF558EHH0SdSg1pqNBfACyr9SgH1ueOb8i9rmk34DbgHOA7wD8BfRspr6SItWoVVsasrIR994V+/WDaNNi0Kepkqk9pA8fvyT1q+gPQKfe8E1D7frrNhEK/Y5+bZ4H+hF8Sf6OiouKL5+l0mrQLb0hFo6wMbroJLroIrrkmbHhy441w7rlQ4jSPRpPJZMhkMgV9Rr49+k7ANOAs4Gjg0hrHewL3AwOAVkAGuBBYWetz7NFLCTJ/flgwrbo69O+POirqRMnUXDdMtQd+BewNfAacDfwVuBxYDcwm/DI4E9iW+7N31/E5FnopYaqr4f77w9IKgwfDT34SbrxS4/HOWEmxsHlzWEbhpz+FCy+Ea68NM3VUOO+MlRQLu+0WdrRatizcVdu9O/zsZy6YFhXP6CU1ucWLw4blH38czvSPOy7qRMXL1o2k2Mpm4ZFHYPLkcIZ/yy3hTlvtGls3kmIrlYIxY2DFCjj2WBg+HC67LJzlq2lZ6CU1qzZtQhtn5cpwlt+jR5iOuXVr1MmSy0IvKRJdu8J//Ac89xw8/TT07h1aO3Z0G589ekmxMGdOONPfc0+49VY49NCoE8WTPXpJRev448N2hmecAd/5DlxwAaxZE3WqZLDQS4qN0lK4+OKwYNruu0PfvmH9nC1bok5W3Cz0kmKnvDwsn7BgASxZEi7Y3nef/ft82aOXFHvPPx/696WlYYbOkUdGnSg69uglJdLw4bBwYdio/Hvfg7Fj4Z13ok5VPCz0kopCSUnY0aqyMtxZO2BAWCxtw4aok8WfhV5SUenQASoqwgyd998PRf/nP4ft26NOFl/26CUVtUWLwoYnGzeG+fcjR0adqGm5qJmkFimbhd//Hq6+OkzJnDEDDjkk6lRNw4uxklqkVCpcpF2xAoYNg6FDw1n+2rVRJ4sHC72kxGjXDq66KhT8LVvC/PvbboNt26JOFi1bN5ISa9kymDQJ3n03rH9/0knh7L+Y2aOXpFqyWXjiiVDw99svXLDt2zfqVPmzRy9JtaRScOKJsHQpnHoqjBoFEyfC//5v1Mmaj4VeUovQujX88z/DqlVhLn7v3vDjH0NVVdTJmp6FXlKL0rlzaN+89BK8/HLYt/bBB5O9YJo9ekkt2ty5YSpmx45hwbRBg6JOtHP26CVpF40YAa++Ct//fujhjxsXllZIEgu9pBavVatQ6Csr4Vvfgv794frrw7IKSWChl6ScTp3g3/4tbHayenW44eqXv4Tq6qiTFcYevSTV4+WXQ/9+69bQvx8+POpE3jAlSY0um4UHHggLph1+eNji8MADo8vT3BdjxwD31nNsArAIeAk4qYDvkKRIpVJw1llh/v3AgTB4MEyeDJ9+GnWyry/fQj8LuIm6f6vsBVwGDAWOB24G2uT5PZIUC+3bhx2tli+HTz4JG57ceSd8/nk4ns1mWbNmDRtiuOVVvoX+ReAS6i70R+SObwPWA6uBfnl+jyTFyl57hR2t5swJa+D37w/33beW3r2P4IAD+rH77nsxZcrUqGP+jdIGjl8A/Eutn50PPAik63lPJ2BdjdcbgPI8sklSbB16KPzpTzB7Npx11haqqn5DNnsI8BG3334MRx45kFNOOSXqmEDDZ/T3AH1rPV5t4D3rCcV+h07AJ/kGlKS4SqXglFOgtHQI2eyehCbHHmzadAaLFr0SdbwvNHRGn4+FwI1AW6Ad0BNYXtcfrKio+OJ5Op0mnU43QRxJalr77bcPK1fOBsYB29htt+fp1u3cRvnsTCZDJpMp6DMKmV55DDARODv3+nJCP342cCFwEeFvDDcCD9fxfqdXSkqEJUuWMGLEiWSz/amufo9Bgw7gqaceprS08c+lnUcvSRH56KOPWLBgAeXl5QwbNoySkqZZeMBCL0kJ5+qVkqSvsNBLUsJZ6CUp4Sz0kpRwFnpJSjgLvSQlnIVekhLOQi9JCWehl6SEs9BLUsJZ6CUp4Sz0kpRwFnpJSjgLvSQlnIVekhLOQi9JCWehl6SEs9BLUsJZ6CUp4Sz0kpRwFnpJSjgLvSQlnIVekhLOQi9JCWehl6SEs9BLUsJZ6CUp4Sz0kpRwhRT6McC99RybBbwCzAWeBcoK+B5JUgHyLfSzgJuAVD3HBwDHASOAkcD6PL8ncplMJuoIX4s5G08xZARzNrZiyZmPfAv9i8Al1F3oS4CDgbuBF4DxeX5HLBTL//nmbDzFkBHM2diKJWc+Gir0FwDLaj0GAg/u5D27AbcB5wDfAf4J6FtwUklSXkobOH5P7rErNhMKfVXu9bNAf8IvCUlSEUkD99fx857Aa4S/LbQmtHl61vHnVgNZHz58+PCxS4/V7KKGzuh3ZseX7nB5LsBs4L+Al4BtwC+BlXW8/6ACvluSJEmSJEmNZjDh5qnaRgMLgfnAhc2aqG715bwcWJ47Nhc4pDlD1dAa+DXwPLCAMH41xWU8G8oZl/FsBfx/wjTgeUDvWsfjMp4N5YzLeO6wB/BeHTniMp5Qf8Y4jeXiGjlqT4qJ01gCcBWwlBCoptbAfwPluecLCYMflfpyQihahzVvnDqdD9yae94ZeKfGsTiN5/nUnxPiM56nAj/PPT8GeKTGsTiN585yQnzGE8JYPQys4m+LZJzGs76MEJ+xbEco9HXZ5bFsjrVuVgOn89Wbq3rmjq0jXLR9ARjeDHnqU19OCPcOXEs4m7qmOUPV8jvg+tzzEuDzGsfiNJ47ywnxGc9HgYm5592AT2oci9N47iwnxGc8AWYAdwFrav08TuNZX0aIz1j2J9yTNAf4E6HbsMMuj2VzFPqH+Op/6BDWv1lX4/UGwm+oqNSXE8I00omE5RyOAk5qrlC1bAI2Ap0IxfRHNY7FaTx3lhPiM54A2wkzw24D7qvx8ziNJ9SfE+IznucDHwFP5V7XPGmKy3ieT/0ZIT5juYnwC+l44GLCumI76vUuj2WUq1euIxSCHTrx1TOVuJgFrCX89nycaP9q903CTWj/Bfy2xs/jNp715YR4jSeE//gPISzb0T73s7iNJ9SdE+IznuOBbxN6yocCv+LLlkJcxnNnGSE+Y/kGXy4a+d/Ax8DeuddxGcuv6EaYV19Ta8L/mM5AG8Jql3sTrW58NWc5ocfcgfDb/3eEpR2isCfhnoQRdRyL03juLGecxnMcMCX3vAz4H0JvFOI1njvLGafxrKn2hcw4jecOtTPGaSwnAnfknu9D+O+pVe51HMcSCAV0x0XOscCE3POTCRcSXiEskha1btSdcywh5zxgavPH+sIs4EO+vBI/Fzib+I1nQznjMp7tgQeA5wj/v48mnv9+NpQzLuNZ01ygO/Eczx3qyhiXsSzly5lrzwNDiPdYSpIkSZIkSZIkSZIkSZIkSZIkSVLz+D8hMXCye/TurwAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# 5. Plot result\n", "\n", "# compute y1 with func1, that's our approximation\n", "x1 = x\n", "y1 = list(map(lambda t: func1(t), x1))\n", "\n", "plt.plot(x,y) # original data seires\n", "plt.scatter(x,y1) # compare with our approximation with n=3\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Step 2. Generalize it to arbitrary number of basis functions: $n = 2^l - 1$" ] }, { "cell_type": "code", "execution_count": 48, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[(0.0, 0.25), (0.0, 0.5), (0.25, 0.5), (0.0, 1.0), (0.5, 0.75), (0.5, 1.0), (0.75, 1.0)]\n" ] } ], "source": [ "# 1. Gererate hierarchical intervals (for hat functions)\n", "def gen_interval_recursive(a, b, max_level, l, intervals):\n", " m = (a+b)/2.0\n", " if (l > max_level):\n", " return intervals\n", " else:\n", " gen_interval_recursive(a, m, max_level, l+1, intervals)\n", " intervals.append((a,b))\n", " gen_interval_recursive(m, b, max_level, l+1, intervals)\n", " \n", "def gen_intervals(a,b,max_level):\n", " intervals = []\n", " gen_interval_recursive(a, b, max_level, 1, intervals)\n", " return intervals\n", " \n", "# TEST\n", "intervals = gen_intervals(0.0, 1.0, 3)\n", "print(intervals)" ] }, { "cell_type": "code", "execution_count": 54, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[1.00000000e+00 5.00000000e-01 0.00000000e+00 2.50000000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.25000000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 6.25000000e-02\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", " [0.00000000e+00 1.00000000e+00 0.00000000e+00 5.00000000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 2.50000000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.25000000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", " [0.00000000e+00 5.00000000e-01 1.00000000e+00 7.50000000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 3.75000000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.87500000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", " [0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 5.00000000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 2.50000000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", " [0.00000000e+00 0.00000000e+00 0.00000000e+00 7.50000000e-01\n", " 1.00000000e+00 5.00000000e-01 0.00000000e+00 6.25000000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 3.12500000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", " [0.00000000e+00 0.00000000e+00 0.00000000e+00 5.00000000e-01\n", " 0.00000000e+00 1.00000000e+00 0.00000000e+00 7.50000000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 3.75000000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", " [0.00000000e+00 0.00000000e+00 0.00000000e+00 2.50000000e-01\n", " 0.00000000e+00 5.00000000e-01 1.00000000e+00 8.75000000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 4.37500000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", " [0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 5.00000000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", " [0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 8.75000000e-01\n", " 1.00000000e+00 5.00000000e-01 0.00000000e+00 2.50000000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 5.62500000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", " [0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 7.50000000e-01\n", " 0.00000000e+00 1.00000000e+00 1.11022302e-15 5.00000000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 6.25000000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", " [0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 6.25000000e-01\n", " 0.00000000e+00 5.00000000e-01 1.00000000e+00 7.50000000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 6.87500000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", " [0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 5.00000000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 7.50000000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", " [0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 3.75000000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 7.50000000e-01\n", " 1.00000000e+00 5.00000000e-01 0.00000000e+00 8.12500000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", " [0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 2.50000000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 5.00000000e-01\n", " 1.11022302e-15 1.00000000e+00 1.11022302e-15 8.75000000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", " [0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.25000000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 2.50000000e-01\n", " 0.00000000e+00 5.00000000e-01 1.00000000e+00 9.37500000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", " [0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", " [0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 9.37500000e-01\n", " 1.00000000e+00 5.00000000e-01 0.00000000e+00 2.50000000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.25000000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", " [0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 8.75000000e-01\n", " 1.11022302e-15 1.00000000e+00 1.11022302e-15 5.00000000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 2.50000000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", " [0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 8.12500000e-01\n", " 0.00000000e+00 5.00000000e-01 1.00000000e+00 7.50000000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 3.75000000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", " [0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 7.50000000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00\n", " 0.00000000e+00 1.11022302e-15 0.00000000e+00 5.00000000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", " [0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 6.87500000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 7.50000000e-01\n", " 1.00000000e+00 5.00000000e-01 0.00000000e+00 6.25000000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", " [0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 6.25000000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 5.00000000e-01\n", " 2.22044605e-15 1.00000000e+00 0.00000000e+00 7.50000000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", " [0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 5.62500000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 2.50000000e-01\n", " 0.00000000e+00 5.00000000e-01 1.00000000e+00 8.75000000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", " [0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 5.00000000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", " [0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 4.37500000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 8.75000000e-01\n", " 1.00000000e+00 5.00000000e-01 0.00000000e+00 2.50000000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", " [0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 3.75000000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 7.50000000e-01\n", " 0.00000000e+00 1.00000000e+00 2.22044605e-15 5.00000000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", " [0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 3.12500000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 6.25000000e-01\n", " 0.00000000e+00 5.00000000e-01 1.00000000e+00 7.50000000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", " [0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 2.50000000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 5.00000000e-01\n", " 0.00000000e+00 1.11022302e-15 0.00000000e+00 1.00000000e+00\n", " 0.00000000e+00 1.11022302e-15 0.00000000e+00]\n", " [0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.87500000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 3.75000000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 7.50000000e-01\n", " 1.00000000e+00 5.00000000e-01 0.00000000e+00]\n", " [0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.25000000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 2.50000000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 5.00000000e-01\n", " 2.22044605e-15 1.00000000e+00 0.00000000e+00]\n", " [0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 6.25000000e-02\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.25000000e-01\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 2.50000000e-01\n", " 0.00000000e+00 5.00000000e-01 1.00000000e+00]]\n" ] } ], "source": [ "# 2. construct System matrix A:\n", "# | h0(x0) h1(x0) h2(x0)|\n", "# | h0(x1) h1(x1) h2(x1)|\n", "# | ... |\n", "# | h0(xm) h1(xm) h2(xm)|\n", "\n", "\n", "maxlevel = 5 # number of basis functions n = 2**level -1\n", "n=2**maxlevel - 1\n", "x = np.linspace(dmin, dmax, n+2)[1:-1]\n", "xExtend = np.linspace(dmin, dmax, n+2)\n", "y = np.sin(x)\n", "# generate list of hierarchical intervals\n", "ints = gen_intervals(dmin, dmax, maxlevel)\n", "\n", "A = np.array([[hat(l[0], l[1], i) for l in ints] for i in x])\n", "\n", "print(A)\n" ] }, { "cell_type": "code", "execution_count": 55, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[ 3.74860583e-03 2.91300418e-02 1.06751262e-02 2.07106781e-01\n", " 1.59764555e-02 7.03261419e-02 1.88455141e-02 1.00000000e+00\n", " 1.88455141e-02 7.03261419e-02 1.59764555e-02 2.07106781e-01\n", " 1.06751262e-02 2.91300418e-02 3.74860583e-03 1.22464680e-16\n", " -3.74860583e-03 -2.91300418e-02 -1.06751262e-02 -2.07106781e-01\n", " -1.59764555e-02 -7.03261419e-02 -1.88455141e-02 -1.00000000e+00\n", " -1.88455141e-02 -7.03261419e-02 -1.59764555e-02 -2.07106781e-01\n", " -1.06751262e-02 -2.91300418e-02 -3.74860583e-03]\n" ] } ], "source": [ "# 3. Solve A*c = y\n", "# But this is NOT a square system, it's an overdetermined system\n", "# We will solve it as a least square problem: x,res,rank,s = np.linalg.lstsq(A,y)\n", "\n", "# OR\n", "\n", "# 3. Solve A^T * A * c = A^T * y\n", "# Transform this into a square system and use x = np.linalg.solve(M,rhs)\n", "\n", "cc= np.linalg.solve(A,y)\n", "\n", "print(cc)" ] }, { "cell_type": "code", "execution_count": 56, "metadata": {}, "outputs": [], "source": [ "# 4. Define the resulting interpolation function of f=sin(x)\n", "def func2(x):\n", " result = 0.0\n", " for i in range(len(cc)):\n", " result += cc[i]*hat(ints[i][0], ints[i][1], x)\n", " return result\n" ] }, { "cell_type": "code", "execution_count": 57, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 57, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD/CAYAAADllv3BAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XuczHX7x/HX7LLssuRQOXUiuUVJKce0VHIskZ90IukuOnJXNyWJSoW6VXepjQ4qd6UbiSSHJVYokWNFd1FIzuwBuzu/P67d7K49zB5mvjPfeT8fj3k0h+/MXLuba75zfT6f6wMiIiIiIiIiIiIiIiIiIiIiIiIiIq7WHFiUx/2DgfWZjy0CzgtkUCIiUjyPAN8DiXk8NgVoGthwREQEIKIEz90C9AA8eTx2CfAo8BUwtATvISIiAXY2sDyP+x8HqgJlgc+ALgGMSUQkrJXkjL0gE4B9wHFgNirLiIgETBk/vGZlrPZ+PpAMtAcm5T6oXr163q1bt/rh7UVEXG0rcG5BB5TGGbs38799gDuBg1hdfRGwBJsdM/ekyLZuxev1Bv3liSeecDwGxakYFafizLoA9QpLyiU9Y/8FaJV5fWq2+6fmui0iIgHirxq7iIg4RIm9EHFxcU6H4BPFWXpCIUZQnKUtVOL0RV5z0APFm1kvEhERH3k8Higkd+uMXUTEZZTYRURcRoldRMRllNhFRFxGiV1ExGWU2EVEXEaJXUTEZZTYRURcRoldRMRllNhFRFxGiV1ExGWU2EVEXEaJXUTEZZTYRURcRoldRMRllNhFRFxGiV1ExGWU2EVEXEaJXUTEZZTYRURcRoldRMRllNhFRFxGiV1ExGWU2EVEXEaJXUTEZZTYRURcpqSJvTmwKI/7uwErgURgQAnfQ0REisBTguc+AtwCHAFaZbu/LLARaAYkA8uArsDuXM/3er3eEry9iEj48Xg8UEjuLskZ+xagRx5v0DDzsYPAcWAp0LYE7yMO2fPDD3SsfjaVPRWoH1WVZa+95nRIIuKDMiV47n+Bs/O4vxKW1LMcBiqX4H3EXzIy4OOPYcsWuOgi6NIlx8NdL7mC1UnXcZzHOXR8JR0H3cqGiy7izJYtTxy0aRPMmgXR0XDTTVCtWoB/CBHJrSSJPT8Hgdhst2OB/X54HykJrxd694bPP4eUFEvMd98N48ZBWhopK1dyVtKfdOJUFrGFVmziEPXYMHYsZ777LlSsCEuXwjXXwLFjUKYMPP00fP89nHaa0z+dSFjzR2LfDNQHqgBJWBlmbF4Hjhw58q/rcXFxxMXF+SEcydPatTBnDiQn2+2kJHjxRZg7F37+mfI1a3ITXr7nECmU51XuYQxv03L5ckvcNWrA3r0nnp+WZpcXX4QxY5z7uURcJiEhgYSEhCI9pySDp2ClmA+wwdM+QEUgHhssHYHV8CcBeRVnNXjqpAULoFs3O1vPEhUFU6faWXiFCoztch0j56wgmQGUoT5RERfy+64GnFK1HGzdCq1bw549OV/31lvh3XcD+7OIhBFfBk9LmthLQond37xe2LAB9u+HJk2gUiW7b9YsePRR2LwZ0tPt2IgIqF3bEnbZsn+9xBdPP82yL+ZRq04dEiPeYt+BKGbMsMoLDzwA8fEnPhwiI6FCBRg2DO67z65nZMB338HRo9C0qZV8RKTYlNjDmdcLN98MM2daFo6MtBLJpEmQmgqjR8NZZ8GNN8K2bdCwoQ2k1q2b70sePw5du9rTXn8dPMePwaBB8OGHUK4cjBoF7drByJGweDH84x8wfbrV3SMj4ZRTIDHRPkBEpFiU2MPZRx9B//5WO88SFWVlkl697Ay9GA4fhrZt4YYb4LHHCjhw7Vro08dmzWQpUwY6drRvDCJSLP6exy7B7KefTgxsZomIsJkwxUzqALGxNub65pvwzjsFHNikiV2yS0uDjRuL/d4i4hsldrc65xzwZPtQj4iwckspqFnTkvsjj8C8eQUceNllEBOT876zziqVGEQkfyrFuNFPP1kxPCbGzpDLlbNT7cWL4dxzS+1tli6FHj0suV90UR4HpKXZAfPn2wdL1ar2LeLll61MIyJFphp7OFq40JLmqFFw112wcyccOAD16lmNvZRNmwYPPgjLluVzMu712uDs0aMWw4YNcO210LcvPPFEicpCIuFIid3tDh6EL76w6x062IDp44/bXPT27QMWxr/+BW+8Ycm9ShUfnvDHH3D99VCnDkyebDNldu2C5s2hQQO/xysSypTY3WzHDrjkEjhyxM6KMzKs+D13LtSvH/BwhgyBb76xskz58j48ITUVBgyAzz6zlgSRkTan/r33rHwjInnSrBg3GzbMVn0eOWJTGlNSrNDtQFIHazFTsybcdpt9xhSqfHmbQ5+SYpcjR+y//frZB5WIFJsSe6j69VcbnMxu505nYsFK5e+8YxWVhx/28Um7dmUuYc0mKclWQolIsSmxh6r27a18kSUmBq6+2rl4sJPwGTOsYeS//uXDEy677OSz87p1/TLIKxJOlNhDkddrZ+ennmqnypGRcN11MHy405FRtaol9nHjbMZMgS680EZdy5e3n+HUU+3+ffv8HqeIm2nwNBQ9+6zNfFmyxJpqeb02Vz2IrFljE3U++QQuv7yQgzMyrL5eoYL1l1m5Er780sdRWJHwolkxbjRlik1pTEyEWrWcjqZA8+ZZF9+EhCIses3IsJ2Y0tKsuVj2cpOIaFaM63z5JTz0kK3nD/KkDnbG/vzz0LlzEcZ1s0Zh9+61OZT68BcpMiX2YJaRAb//bguR1qyxNrzTpsH55zsdmc/69oU77rDtVA8f9vFJ5cpZu9+FC2H8eJvnvm2bzX0XkUKpFBOsfv/dZr5s327T/8qXh7fesn65Icbrte4G27ZZx95s+3gUbPt2aNbMPhE8Hrv85z/WB0ckTKkUE8puvNF2M0pJsXrzsWNBN0DqK48HXn3VpqzfdVcRqitVqpxYuJScbHPce/eG3bv9Gq9IqFNiD1Zr1pzYtg4ssX/zjXPxlFCZMjYWum6dbbDkk//97+TB07Jl4YcfSjs8EVdRYg9WZ5yR83aFCtZjPYRVqGCtYd57zzbqKFStWievQk1NhTPP9Et8Im6hxB6sHnrIahgVK1pGbN3a5g6GuNNPtwVMw4fb5J4CVasGEybYXP3YWJsx07y5NusQKYQGT4PRvn3QtKltPn3KKXZp0cJVvcuXL7e27J9/buOjBdqyxfZOrV7d9mudNAmuuSYgcYoEGy1QCkVer818OeMMHxuuhK4ZM2DQINuJqW5dH5+0aJFN+1y9GmrU8Gt8IsFIs2JC0euvw88/w3PPOR2J33XvDo89Bp06WQdin7RrB/372wR5n/oDi4QfnbEHk/XrLXEtXRpWOwkNHWptbxYssHJ6odLS4Ior7JPB5x7BIu6gUkwoSU62NrYPPWSbTYSRjAwbF05JgY8/9rE9zK+/wqWXwuzZ9l+RMKHEHkruvhsOHYL337fZMGHm6FEryTRqBC+95OOvYNo0O91fvRoqVfJ7jCLBQIk9mO3ZY7Xib7+FypVtheX69WGdoA4csBa/ffvaFxefDBgAixfbN57q1eHf/4Y2bfwap4iTfEnsZQp6UPwkI8P6wGzebAtwduyw5B7OH3TYrM7PP4dWraB2bejTx4cnHT5srRe8Xvs9XnONncGH0RiFSG6aFeOE336zudnZV1V6vbBqlXMxBYk6daxs/sADNrOxUJ9+mvMDMT3dPh1EwlhxE3sEMBFIBBYB9XI9PhhYn/nYIuC84gboSjExOfvAgJ3Fx8Q4E0+QueACa+LYu7dVpwqUuzFaRIR+jxL2ipvYuwNRQCtgKDA+1+MXA7cC7TIvPxY3QFeqXt32KM0SHW37fzZv7lxMQaZ9e1uf1bmzdTDO1+jRJxK5x2PtjXv3DkiMIsGquIOn44EVwEeZt38D6mR7fCOwAagBzAaezeM1wnfw1Ou1nuLlylmS/9vf4J57QrYtrz899xx88IHNc69cOZ+DPvvMLmCbrK5bp1Wp4lr+HDytBBzKdjsdO/vPWgo4Ffg3cBiYDnTBErwAfPSR7Trx7bcQFeV0NEHtkUdsv42ePa1pWJ6/rq5dT2y+UaWKFeg//DCgcYoEk+Im9kNAbLbb2ZM6wAROJP7ZQFPySOwjszXmjouLIy4urpjhhJB9++DBB23rNyX1Qnk81uCxZ0+bHTplSiFz3EeMsLLWrFnQrVvA4hTxl4SEBBISEor0nOKWYnoA3YDbgRbA49hZOUBl4HvgfCAZK9dMAubmeo3wLMX0728taCdMcDqSkJKcDFdeaR0XnnmmkIMXLbLJ8GG+LkDcyZ8LlDzAq8CFmbdvBy4BKgLxQB9sZsxRYD7wZB6vEX6JfcECS+zr11tylyLZs8fmuA8eDAMHFnLwgAE2kPrKKwGJTSRQtPI0mCQnW4lgwgTo0qXw4yVPW7fa6tSJE62fe77274fGja35TKtWAYtPxN/UtjeYPPmkNatSUi+RevVg5kw7IV+xooADq1Sx+ZJ33mmNaETCiM7Y/c3rtY2pO3a0aXinneZ0RK4we7Yl9yVLoH79fA7yeq217yWXwOOP231h2GBN3EVn7E6aP982+IyMtP1K//lPJfVS1KWLfQnq1Al2787nII/HWkWOGWNrBGJibP5kOJxQSFhTYveHX36xlaW7d1sSSUmBN95wOirX+fvfrVFY166QlJTPQW+/be0bjh+H1FTr/vjaa4EMUyTglNj94euvT94t4uef4eBBZ+JxsVGj4Pzz4cYbbWOlk8yalbPZWnKyFelFXEyJ3R9OPfXkr/seD1So4Ew8LubxQHy8jY/ee28eVZaaNXPW1SMirCewiIspsftDu3bQsKEllKgoq+2++CKUUft7fyhb1jZTWrHCyuk5jBtni5Sio+1bVFSUFedFXEyZxh/S021HpEcesTPG5s2hRQuno3K1SpVspkyrVtbT/bbbMh9o0AA2brTyS2qqLVvduxfOOMPReEX8SdMd/eGFF2DePNvwQdPrAmrjRvvC9N57cPXVeRzwxhvWcGbJEv1tJCRp5akTdu2ynSKWLtX2bA5ZsgRuuAG+/BKaNMn1YHo6XHYZDBkCN9/sSHwiJaHE7oTbb7fB0+efdzqSsPbxx5a7ly2DM8/M9eDy5dCrF2zapJ49EnKU2APt66+tv+ymTeoqGARefBHefNO+PFWpkuvBvn1t/OPZvPaAEQleSuyBlJFhg6T33w+33up0NJJp8GBYvdqGPHJsULVzp5XMEhPhPG3JK6FDLQUC6a23bN7dLbc4HYlkM368dXLo29c+e/9SsyYMHWqbnrjpBEMEnbGXjgMHbN/SOXPg4oudjkZySU21GTItWsDYsdkeOHbMWimPHavdliRkqBTjbxkZNsvioYcse7z+utMRST727bNebAMHWrXsL/PmwaBB1nmzTBn71iUSxPy5mbWMHg1PPWV9SMqUgR9/dDoiKUDVqrasoE0b6yjQs2fmAx062E5LsbFWkmnb1vajPeUUR+MVKQnV2Ivjk09sNsWxY5YMvF4YPtzpqKQQZ59tPcHuvtumQQIwd65ty5Sebt/AEhOhXz8HoxQpOSX24liwwLoEZklLg4ULnYtHfNa0qa1K7dkTNm8GEhKsjJbl2DFb4SQSwpTYi+OMM3LNnQNq1HAmFimya66xL1ydO8Ou2PrWICw7bYgiIU6Dp8Vx+DA0awY7dlgZxuOBRYvsPgkZo0fDjP9msPhYSypu22hn62lpdhZ/+eVOhyeSJw2e+ktsrO1jOmuWlWSuvFLdAkPQ8OGwbVsEvbYt59PHp1P2yH54+WXYvt3p0ERKRGfsEtbS0mwXwxo1rP2AZ+lX1hxs82broy8SZLTyVKQQZcrAhx/C2rW2zR6XX26tIV54wenQRIpNZ+wiwB9/QMuWVp7pH/czXHqpLVqqVcvp0ERy0MpTkSL48Udbn/T229Bx0T/hzz9h8mSnwxLJQYldpIgSE6F7d5j78WEu7l1f/X8k6CixixTD9Olw772wdOD7nDM/3qayahs9CRIaPBUphuuvh2HDoNN7N7F313GYMcPpkESKpLiJPQKYCCQCi4B6uR7vBqzMfHxAsaMTcci998J113m4NvIzUv4xHI4edTokEZ8V9/tlD6Ar0B9oDgwDumc+VhbYCDQDkoFlmcfuzvUaKsVIUMvIsH1TUhcu5aE6/+Dz6HJc2vZyrn36aadDkzDmzxr7eGAF8FHm7d+AOpnXLwSeAzpl3n4BO3Oflus1lNgl6KUkpVG38gq6pa/mE37hKB9zX4tzGLN8sdOhSZjyZ429EnAo2+30bK9VCTiY7bHDQOVivo+Io1ZNfo3D6XfzGT1pS2uS+IbxXy/n8I4dTocmkq/i9oo5BMRmux0BZO0oeTDXY7HA/rxeZOTIkX9dj4uLIy4urpjhiPjHoT//JJKKJBPNSppTi6PsoxxHdu8mVouXJAASEhJISEgo0nNKUmPvBtwOtAAeB7pkPlYW2IDV3pOwMkw3YGeu11ApRoLen5s2Uf/8ZhzkVW7hCLPoTeUy/fjl6Kd4IjSpTALPnzV2D/AqVk8HS/CXABWBeGywdAR2Jj8JeC2P11Bil5Dw3dSp3NLvPv44lswrntbcX3Emi5bH0KiR05FJONICJZHSNnMm79+zjEcjn2P5co9ayUjAaYGSSGm79lpuPu8bBjZdQefOcOhQ4U8RCTSdsYsU1Zo1eDtcwz3dtvHTtnLMng1RUU4HJeFCpRgRf7njDtJPqUaPLc9TuTK8847ayUhgKLGL+MvOndC4McmLV9F+QF2uugqeesrpoCQcqMYu4i81a8KQIcQ88TCzZtkuTK+/7nRQIkZn7CLFlZICf/sbTJnC1tptadMG4uOha1enAxM30xm7iD9FR8Ozz8KQIdQ7J4OZM6F/f1i50unAJNwpsYuUxI03QtmyMGUKl13qZfJkuO462LLF6cCkREK8mqDELlISHg88/DAMGACRkXTtfxojb9xMp062ZaqEmO3boWlTKFMGqlWzrRFDkGrsIiXh9UK9evC//524r0IFhvf/nfkrK7NwIcTEOBeeFFGjRvDDD5CebrdjYmDtWjj3XGfjykY1dhF/+/NPyN3CNzKS0XELadAA+vQ5kSMkyB05Aj/+mPMPFhEBX3/tXEzFpMQuUhKVK59cj01Px3PaqcTHQ3Iy3HdfyJdsw0NMDERGnnz/qacGPpYSUmIXKYly5WDcOEsK5ctbzf3ii6F1a6Ki4JNPIDERnnvO6UClUBER9rcE6xFRsSK0aQNXX+1sXMWgGrtIaVixAlatssG3OXPgu+9sAA6r1LRqZStTb7nF4TilYA88ANu2wVVXQa1aNsUpyPruq6WASKB5vZYUune3GkymDRugfXv44AO48koH45P8rV8P7drBxo1BXX5RYhdxwoYNEBd3UoJYvBh69YL58+HCC/N/ujjA67VP3B494N57nY6mQJoVI+KERo2s5vLooznuvuIKePll6NLFKjYSRKZNgz174O67nY6kVOiMXcQfDhyAhg1h1ixo1izHQ+PHw1tvwdKlcMopDsUnJyQn29/q3Xft0zfIqRQj4qTJk60r2LJlOQbgvF4YPNjWvcydaxNrxEEjRtj89f/8x+lIfKLELuKkjAxo2RIGDYK+fXM8lJ4OvXtbm5n33w+6iRfh4+ef4bLLYM0aqFPH6Wh8osQu4rSVK22GzObNUKlSjodSUmyKdOvWmufumO7dLbHnGg8JZkrsIsGgf3+oWhWeftrmtmdb3bh3ryX2e+8N+skY7uH1QmoqLFkC99xj0xzLl3c6Kp9pVoxIMBg2DF56yVanRkdbgs9UrRp8/jmMGQPTpzsYY7hYtMh+6RUq2I4o998fUkndVzpjF/G3nj1h5sycHQOnToVrr/3rkG+/hU6dYMYMW6UqfrB7t3XiPHLkxH3VqtnS4Kgo5+IqIp2xiwSDJUtydgxMToaEhByHXHKJzbbr0cO6xoofrFt3cpOvo0fh11+dicePlNhF/K1GjZy3y5WDM8446bCOHeGZZ6BzZ/jjjwDFFk5q1IBjx3Led+wYVK/uTDx+pMQu4m+TJlmnwIoVbfC0QgW46648D+3fH267zVanZq8YSClo1Mh6wXg89jeIibHBjSpVnI6s1KnGLhII27db+cXjgaFDrRtY27Z5Hur12k57u3ZZaT6zSaSU1IEDltyHDbPEfuGFVgMLMZruKBKMZsyw5L5mTb4zMo4ft7HVOnXgjTfs80BK6K67rMb+6qtOR1Ii/krs0cB7wKnAYaAvsCfXMROA1pmPe4HuwKFcxyixS/jq2RPOPx9Gj873kCNHrHVJ9+7w+OMBjM2NFi+2xmzr19uuVyHMX4l9CFARGAX0BloCD+Y65ivgOmBfAa+jxC7ha8cOuOgiWLAALrgg38N27bLpjyNGQL9+gQvPVVJToUkTeP552zgjxPlrumNrYG7m9bnAVXm8Zn0gHlgK3F6M9xBxt1q1bEulO+8scLfrGjVsAdPQofDFFwGMz02eeso+PF2Q1H1VWGK/A1iX61KZE2WVw5m3s4sBXgJuBjoCg4D8T0lEwtWAATb1sZCab4MGtnfqrbfC6tUBis0t1q2zQYqXX3Y6koAqbLx9UuYlu0+A2MzrscCBXI8nY4k9NfP2QqAJ9qGQw8iRI/+6HhcXR1xcnA8hi7hERIQlnTZt7GzyzDPzPbR1a5g4Ebp1sy7AZ58duDBDVnq6fXg+/TTUrOl0NMWWkJBAQq4FbYUpbo09FngSuBG4HLgn2+MNganAxUAkkAAMADbleh3V2EXAEk9iop1Vpqfbsvd8+vi+/LKd4C9bZn3FJA87dtjUxrlz4dNPYeFCV/VF9uesmHeAmsBR4CZgNzAY2ALMwpJ/b+B45rHxebyOErsIwOHDdkZ57Jg1aG/QwJpV5TN74+GH4euv4csvXdm/qvi8XmvqFR9v0xpTUy25X32105GVKs1jFwkFTz5pKyCPHrXb5cpBnz62f14eMjLgppvs5P7DD111MloyM2fCzTdDUpLd9nhs5pHLBibUBEwkFKxceSKpg13/9tt8D4+IgHfegT//hH/8IwDxhYp162z3kixer21wEoaU2EWcdsEFOWsqkZG29L0A5cpZ//Yvv4QXX/RzfKGifv2cG8h6PHDOOc7F4yCVYkSclpQEcXF2dpmRYbXhZcugRYtCn7ptm82YGT8e/u///B9qUEtNtfUBR45Yg6/ISFtx2rix05GVKtXYRUJFerr1jklLg1WrYPJkWL485xloPtautfHBadPy7SsWHgYPtt7qo0bBwYP2TSjXPrNuoMQuEoq8XrjhBpvX7mOdZf58GzdctMha0ISd2bNh0CD47jvXzwNVYhcJVfv2QdOmNmm9SxefnjJlijULS0y0ikTY2LHD2u9+9BFcfrnT0fidErtIKFu61M7cV6/2OVM/84zltyVLXFmFOFl6OnToYG0wR4xwOpqAUGIXCXWjRtkA4Lx5J+/XmQevFwYOhJ9/tupE2bIBiNFJY8bYIqSFC336/biBErtIqEtPhyuvtLPSRx/16SlpabYpdtWqtsbJtZt0LF8O118P33xjO5KECSV2ETf47Teb4REdbQtwrrrKZs3Exub7lKQkaN8errnGTvpdweuF556zAeX0dLu89ZbtRBJGtPJUxA0OH4bkZNi505pbzZplLQcKUKGCHfbBB9Y6xRXi423Hqd27Ye9em6++c6fTUQUlJXaRYDd/fs6GMEeP+rTrxmmnWfl5xAirt4e8qVPtAy5LWpp9cslJlNhFgl2lSicPDMbE+PTUc8+1vbP79bN1TyEt9/x0j8f1c9aLS4ldJNj16gW1a5/oJ+PxQMeOPj+9eXOYNMn28ti61U8xBkKvXvbfiAj7oKtQwXrZy0k0eCoSCo4csQHT3but98nDDxe5QczEifDCC7aAqXp1P8bqD1u22OKjUaPgjz9sILVPH/tKEmY0K0bErbIaxHz8sS3O8dGjj1rbgQULfK7mOG/3bmjVCh55BP7+d6ejcZwSu4ibLVhgO24sWOBzB0OvF267zb4ATJsWAmt6kpKgXTubtzl6tNPRBAUldhG3e/99GDbM6is+LtI5dgw6d7Yd+F55JYgXMKWl2cDAaadZGSpoAw0sJXaRcPD889YBrGNH6+nerBkMHVpgy9+DB63F7803W4UjaHzwgZWXqlWzrxUHD9qG1K7vjeA7JXaRcJCebmfru3fbRh3R0TYVZuHCAs9yf//dStdjxlhFx3EvvgjDh9tcdY/HZr9s2mQ7I8lftPJUJBxs2ACHDllSB2s7sHIl/PhjgU+rXdsWLg0ebAOqjnvmmRMLkLxeS+4zZzobU4hSYhcJdceP51yZCpYU09IKfWrjxvDhh9C7t+0F7ajjx3Pe9nptQECKTIldJNRdcIH1a8+qQ0dEWJL0ccpLXBxMmGD7efz2m//CLFBSkg2SZv+AKl/e+tFLkanGLuIGe/bAfffZafdFF8HFF9ug6tSpNl3QB2PHwrvv2v4elSv7Od7sfvvNZr80agR168L06VClCowbZwPBkoMGT0XC2cKFtjpz9GifFvZ4vfDAA7B+vTUPi4oKQIyrVllP9fvvt9W0mtJYKCV2kXD344/QrZtt1BEVZV0ha9Wy2kvDhicdnp5uLVliYuzsPXfpvkQOHYIHH4QVK2wSfYcOtklrfHzY9VQvCSV2EYH9+y2J//mnzZzxeKxj5KZNULPmSYenpNimTW3bwrPPllIMXq/NrfzuO2s7nDWd8auvoGXLUnqT8KDpjiJiBfM9e05Mh/R6bXB1zpw8D4+OtjVB06fDv/9dSjH8/jusWWNJPSuG6GhITS2lN5DsyjgdgIj4mcdjM2TS00/cd/Qo7NiR71OqV4fPP4c2bWy+e4krJQcOnDydEaCMUpA/lOSM/Xrg/XweuxNYBSwHupTgPUSkpDweGxXNaueY1cv81VdtOuGGDXb/2rXWd2bkSNi2jbp17cz9zjtt3+hCpabCSy/BkCG2uwfYFnZDh1pd56yz7CwdrN1B3brQokVp/7RC8WvsE4AOwHdA7sXINYB5wCVANLAUaAbkXmmgGrtIoHi98OabttS0Th0btIyNteT21lXxAAAGHUlEQVQ+diw0aWL17qNHrfZdsSKsXg116zJnDvTvD0uWwOnlt5Fy4ACnN26MJ/vI6vHjVivfsMESfEyMTVXcsME+PB57zOr548bZp0TDhtY+oGJF534nIcqfg6f/B+wG7gJy76p7LdAJGJh5+7/AM8A3uY5TYhcJBocPw/nn51yd5PHYFMmJEwGIfyODh+7dQfLxlkSyn0bRMcxbu4xqWX1c5syxTT+SknK+xqZNNgNGSo0vib2wAtcdwIO57usHfATE5fOcWOBgttuHgUAudxCRooiNtVky2Xm9Nt/x11+hcWMarlpFpePXc5itpOMhNeVm3rzscv754ECb+J6YmDOpg5V8atcO3M8hfykssU/KvBTFISy5Z4kF9hfxNUQkkG66KWcTruhom+teowasW0fUd2uZzQ+M4xTKcYxHSeSnA3utdNOjBwwcaKtHjxyx50dF2epXlVoc4Y8h6ZXA00A5oDzQEFif14EjR47863pcXBxxcXF+CEdECjVsmNXG4+MtKT/xBNxxhz3WrRtLEpfz+OxU0uhNFGlMIokWFZ+nwzPPnHiNxYuhXz/YudPq7W+/7cRP4joJCQkkJCQU6TklWaB0BVZjzxo8HQxsAWYBA4C/Y7Nungam5/F81dhFQkTqgQNcUacBG5NiiaAaZT0bWPrZR/ytc2enQws7WnkqIqUmLTWVpRMnknzwIC1uuYWq9eo5HVJYUmIXEXEZtRQQEQlDSuwiIi6jxC4i4jJK7CIiLqPELiLiMkrsIiIuo8QuIuIySuwiIi6jxC4i4jJK7CIiLqPELiLiMkrsIiIuo8QuIuIySuwiIi6jxC4i4jJK7CIiLqPELiLiMkrsIiIuo8QuIuIySuwiIi6jxC4i4jJK7CIiLqPELiLiMkrsIiIuo8QuIuIySuwiIi6jxC4i4jJK7CIiLlOSxH498H4+j00AvgEWAQuBSiV4HxERKYLiJvYJwDOAJ5/HLwY6AO2A9sChYr6P4xISEpwOwSeKs/SEQoygOEtbqMTpi+Im9mXAQPJO7BFAfSAeWArcXsz3CAqh8sdWnKUnFGIExVnaQiVOXxSW2O8A1uW6XAJ8VMBzYoCXgJuBjsAg4IISRyoiIj4pU8jjkzIvRZGMJfbUzNsLgSbYh4KIiASxOGBqHvc3BNZg3wbKYmWbhnkctwXw6qKLLrroUqTLFgpR2Bl7QbLeJMvgzDecBbwLLAeOA28Dm/J4/rkleG8RERERERERESm2ghY7OSUCmAgkYgut6jkbToGaYzEGq7LAFGAJsALo5mw4+YoEJmPTdL8CGjkbTqFOA7YD5zkdSAFWY/9vLqLoEzECaRj2b30V0NfhWPLTlxO/y6+BFIJ48ecErAb/gdOB5NID+0cOljhnOBhLQR4Bvsf+pwxW/YAXMq9XAX51LpQCXQe8mXn9CoL3bw72YTkd2EzwJvbyWGIPdnHAp5nXKwBPOheKz14BBuT3YDD0iilosZOTWgNzM6+vAJo5GEtBtmAfQsH2+8vuY2BE5vUIIM3BWAoyE7gr8/rZwH7nQinUWOA1YKfTgRSgCbau5QtgAXaCFIw6YNOxZ2CTPz4t+HDHNcO+Tb5Z2IGBkN9iJ8h/6qST4rEFVll+JTg+CPNyNjYLKdjFYusabnQ6kEK8DRwErnY4jvz0Ax7LvL4IaOBcKAVqjP27B1uNvoXg/DcUj53ElcG+/Wx2NpxC/Rf7Rhn04gi+xD4e6JXt9nanAvHB2QR/Yj8Dq1/2czgOX50O/AJEOxxHXhYDCVhS34/VW093MqB8RGHlmCwrgNoOxVKQMcCQbLfXANUdiqUwpwDrCzsoGD89g8UyoHPm9RZYHVuK53RgHjYe8LazoRToVmwQDWxgKiPzEmyuwE6G2mFJ6DbgDycDysft2AkSQC1soC8YS0dLOfHtvBZWZ9/rXDgFaouVtQpUkgVKpSn3YqdgMB37Kr4s83awNzMLtt9fdo8ClbE6e1atvRMn2k4Ei2nYB89ibHDyAeCokwGFuEnAW9hsKLB/Q8H4QTkbS5grsZPdQQTvv6fzgK1OByEiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiLvb/ndumabUQ3wUAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# 5. Plot result\n", "y2 = list(map(lambda t: func2(t), x))\n", "\n", "#plt.plot(x,y) # original data seires\n", "plt.plot(x,y2, color='red') # compare with our 2nd approximation with n = 2^maxlevel basis functions\n", "plt.scatter(x,y2, color='red')\n", "plt.plot(x1,y1) # compare with our 1st approximation with n = 3 basis functions\n", "plt.scatter(x1,y1) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Lecture Exercise 2: Interpolation with Partial-Hierarchical Basis" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Consider new (partly hierarchical) basis (slide 20): see figure in slide\n", "\n", "Solve the intepolation problem:\n", "1. set up linear system of equations to obtain $c_k$ for n=7\n", "2. Generalize this to $n=2^L-1$ \n", "2. how does this extend to solving the interpolation problem for a hierarchical basis\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Step 1. Solve for system of $n=7$" ] }, { "cell_type": "code", "execution_count": 69, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[1. 0.5 0. 0. 0. 0. 0. ]\n", " [0. 1. 0. 0. 0. 0. 0. ]\n", " [0. 0.5 1. 0.5 0. 0. 0. ]\n", " [0. 0. 0. 1. 0. 0. 0. ]\n", " [0. 0. 0. 0.5 1. 0.5 0. ]\n", " [0. 0. 0. 0. 0. 1. 0. ]\n", " [0. 0. 0. 0. 0. 0.5 1. ]]\n", "[[1. 0.5 0. 0.25 0. 0. 0. ]\n", " [0. 1. 0. 0.5 0. 0. 0. ]\n", " [0. 0.5 1. 0.75 0. 0. 0. ]\n", " [0. 0. 0. 1. 0. 0. 0. ]\n", " [0. 0. 0. 0.75 1. 0.5 0. ]\n", " [0. 0. 0. 0.5 0. 1. 0. ]\n", " [0. 0. 0. 0.25 0. 0.5 1. ]]\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 69, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD/CAYAAADllv3BAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAGIRJREFUeJzt3XmUVPWZh/GnEWSRBjdUTJxDBgTBbVSUVS2NO6CIGgX3fZmYETUuELU10UzGaMTJJGaMMZs6atzHjBq1C1kUREVxIdBtNM5JJs64YUSE0MwfbzU0bUt3V3X1rbr1fM6pYy1ddd/S9svlve/9XZAkSZIkSZIkSZIkSZIkSUq1EUBtC89PBV7NvVYLDO7MoiRJ+bkEeAWY28JrvwJ269xyJEkAXQp4bx0wCahq4bU9gGnALOCyArYhSepkA4BnW3j+CmBzoBvwn8C4TqxJkipaIXvsGzIDeB9YBTyKbRlJ6jRdi/CZfYne+zBgObA/cFvzHxo4cOCa+vr6ImxeklKtHhi0oR/oiD32Nbl/TgbOBD4i+uq1wDPEdMxjn6usvp41a9ak9nbVVVclXoPfz+/m90vfDRjYWigXusf+FjA6d/+uJs/f1eyxJKmTFKvHLklKiMFeJJlMJukSiirN3y/N3w38fpWgpRn0zrIm1y+SJLVRVVUVtJLd7rFLUsoY7JKUMga7JKWMwS5JKWOwS1LKGOySlDIGuySljMEuSSljsEtSyhjskpQyBrskpYzBLkkpY7BLUsoY7JKUMga7JKWMwS5JKWOwS1LKGOySlDIGuySljMEuSSljsEtSyhjskpQyBrskpYzBLkkpY7BLUsoY7JKUMoUG+wigtoXnJwDzgbnAGQVuQ5LUDlUFvPcS4ATgr8DoJs93A14HhgPLgTnAeODdZu9fs2bNmgI2L0mVp6qqClrJ7kL22OuASS1sYGjutY+AVcBsYJ8CtqO2WrYMJk2Cvn3hK1+B3/0u6YokJaCQYL8f+FsLz/chQr3Rx0DfArajtjr2WPjtbyPg33oLJk6EN95IuipJnawYB08/AqqbPK4GPijCdtTUmjWxh/7ZZ+ueW70annoquZokJaJrET5zMbA9sBnwCdGGub6lH6ypqVl7P5PJkMlkilBOhXj9dahq1hVbuRIWLoRVq6Bbt2TqklSQbDZLNptt13sKOXgKMAC4kzh4OhnoDdxKHCy9kvgbwW3Aj1t4rwdPO8Krr8K3vw3ZLIwdC//1X7BiBfToAVttFb32t9+G6dPhpJMMeKnMteXgaaHBXgiDvRCLFsE118Azz8BFF8F550Hv3hHwtbXQrx+cdhr06gWzZsHVV0N9PUybBiefDBtvnPQ3kJQHgz2NXnklAn32bLj4Yjj3XNhkk7a9d86cCPglSyLgTznFgJfKTLHHHdWZFi6MUcaDD4ZRo2Lv++KL2x7qAGPGwBNPwJ13wn33wfbbw09+Er14SalhsJe6l16CI4+Eww6DvfeOQL/oovYFenOjR8Pjj8N//Ac8+CAMGgQ//vH6EzWSypbBXqpefDHm0MeNg333hbo6mDo1euYdZdSoONh6773wyCMR8D/6kQEvlTmDvdS88AIcfjhMmAD77x976Bdc0LGB3tyIEXFi0333waOPRsD/8IcxXSOp7BjspWLBggjzI46AAw+MPfRvfAN69uy8GvbaK4L9/vujVTNoEPzrvxrwUpkx2JM2f360WyZOjAOjdXVw/vmdG+jN7blntGYefDDOZh04EG6+GT79NLmaJLWZwZ6UefPigOhRR0Ww19XB178eJxaViuHD4eGHI+SffjoC/qabDHipxBnsne255+DQQ+GYY6L1UlcXJxeVUqA3t/vusff+6KMwc2YE/A9+AMuXJ12ZpBYY7J1l7txotRx7bLRdli6Nk4u6d0+6srbbbTd44IE40Dp7dgT8DTfAJ58kXZmkJgz2YpszBw46CKZMibbL0qVw9tnlFejN/cM/xATN44/Ds89GwH//+wa8VCIM9mKZNQsOOABOOAG+9rU4jf+ss9J1Cv8uu8BvfhMHWOfNi4C//noDXkqYwd7RnnkGvvrVWGhr8uQI9DPOSFegN7fzznGS05NPxtjm3/89fO978Ne/Jl2ZVJEM9o4yc2acUHTqqXD88fD738Ppp1fWMrk77QR33x0TNC+9FHvw//zP8PHHSVcmVRSDvVDZLGQyEeInnQSLF8dyuZUU6M3tuGOsQ1NbG6tRDhwI110Xl+yTVHQGe2vWrIE//xnef3/9555+OtZwOfPM2EtfvDiWwa3kQG9u2LBYSXLmTHjttTiT9dpr1w/4lSvhj3/07FapAxnsG/LhhzByZPSM+/eHE0+MA4X77APnnBO98zfeiH5612JcZTAlhg6FO+6IA8qLF0fAf+c7sQDZVlvF61tsEbPykgrmhTY2ZMqUGOtrXK+8SxfYcsuY3T7uOMM8X0uWQE1NtGua/g706hUnbPXvn1hpUqnzQhuFeu659S9C0dAA++0XI4yGev4GD4Yrr/z8ejjdusXfgCQVxGDfkIEDYy+9UY8e0TdW4fr3jz8om/rsM/i7v0umHilFbMVsyB/+ED32FSsihLbfPk6lL+ba6JXkZz+Lhc+6dYt1Z7bdNto05XxWrlRkXsy6IyxbFqfN9+gRl5Rz6qVj1ddH+2XAgOi79+4Nt98OVUn+akqly2BXefnkk7iu6+TJ8M1vJl2NVJLaEuweAVTp2GQTeOihaH/tsEMsayyp3dxjV+mZNw/Gj4+TwHbeOelqpJLiuKPK04gRMGNGXNT73XeTrkYqO+6xq3R961uxFs9TTzkpI+V48FTlraEhLiHYp0+MRjopI9mKUZnr0gV++Ut4+eW4QpOkNnEqRqXNSRmp3fLdY+8C3ALMBWqBgc1enwq8mnutFhicb4ES220H998f69wvWpR0NVLJyzfYJwIbA6OBy4Abmr2+O3AisF/utiTfAiXASRmpHfIN9jHAY7n784DhzV7fA5gGzCKCXyrclClx2cFJk2LBMEktyjfY+wBNr3O2utln3QWcDewPjAXG5bkdaX3XXBMX5zjnnPXXcpe0Vr4HT5cB1U0edwGarsE6g3XB/yiwW+6f66mpqVl7P5PJkMlk8ixHFaNLF/jVr2Ds2LjgycUXJ12RVFTZbJZsNtuu9+Q7GDwJmACcCowErmDdXnlf4BVgGLAcuAe4jXWtm0bOsSt/77wTkzK33OKkjCpKMU9QqgJ+BOySe3wq0VfvDdwKTCYmYz4DngSubuEzDHYVZt68CPWnnnJNGVUMzzxV+t15J0yfDvPnQ79+SVcjFZ3BrsowfTo88ww8+aRryij1DHZVhoYGOPpo6NvXNWWUeq4Vo8rQOCmzcGFMykgVzrVilA6bbAIPP7xuTZnx45OuSEqMrRily3PPxbIDTz8NO+2UdDVSh7MVo8ozciT84AcxBvm//5t0NVIi3GNXOjkpo5RyKkaVq3FSZtNN4bbbnJRRatiKUeVqvPrSiy86KaOK41SM0qt375iUGTXKSRlVFFsxSj8nZZQitmIkWDcpc/jhTsqoIrjHrsoxbRrMmuWkjMqaUzFSUw0NcNRRsNlmTsqobNmKkZpqXFPmxRfhxhuTrkYqGqdiVFkaJ2VGjoQhQ5yUUSrZilFlevZZOOIIJ2VUdmzFSF9k1CgnZZRa7rGrsjkpozLjVIzUGidlVGZsxUitcVJGKeRUjNR0UmaHHWDcuKQrkgpiK0Zq5KSMyoCtGKk9Ro2KdoyTMipz7rFLzU2bBrNnx6TMxhsnXY20HqdipHw0Tspsvjn89KdOyqik2IqR8tE4KfPCC3ESk1RmnIqRWtJ8TRknZVRGbMVIG9I4KVNbCzvumHQ1UlFbMV2AW4C5QC0wsNnrE4D5udfPyHMbUvIaJ2UmTHBSRmUj3z32ScB44DRgBHA5MDH3WjfgdWA4sByYk/vZd5t9hnvsKh+XXx6TMlOmwJ/+FC0a2zNKQDGnYm4A5gH35B7/N/Dl3P1dgO8Bh+Ye30jsuf+m2WcY7Cofq1ZB//7w4YewejX06gWXXAJXXZV0ZaowxWzF9AGWNXm8usln9QE+avLax0DfPLcjlYaZM2HFigh1gOXL4dpr4zmpxOQ7FbMMqG7yuAvQkLv/UbPXqoEPWvqQmpqatfczmQyZTCbPcqQiW7YMNtpo/eeqquCTT6BHj2RqUkXIZrNks9l2vaeQHvsE4FRgJHAF0Nhw7Aa8RvTePyHaMBOAPzf7DFsxKh//8z8weDB8/PG65wYPhsWLPYFJnaqYrZgHgBXEgdEbgKnAZOBMYBVwIfA4Eeq38flQl8rLNtvAU0/F6o99+sDQodGGee+9pCuTPsc5dilfl18Oc+fC737nmjLqNK4VIxVTQwNMmgRbbOGaMuo0rhUjFVOXLvDrX8OCBa4po5LiWjFSIRrXlBk1Kvrvhx2WdEWSrRipQ8ydCxMnuqaMis5WjNRZRo+GG26Iqy/93/8lXY0qnHvsUke67LJYEdJJGRWJUzFSZ2togCOPhH794NZbnZRRh7MVI3W2xkmZ55+Hm25KuhpVKKdipI5WXb1uUmbIECdl1OlsxUjFMnduXH0pm3VSRh3GVoyUJCdllBD32KVic1JGHcipGKkUOCmjDmQrRioFjZMy8+c7KaNO4VSM1Bmqq+GRR+Ii2DvsAIce2vp7pDzZipE605w50ZZxTRnlyVaMVGrGjIHvf99JGRWVe+xSEi69FJ57zkkZtZtTMVKpWr06WjJbbw3//u9OyqjNbMVIpWqjjeCOO2DePJgxI+lqlDJOxUhJaTopM2SIkzLqMLZipKQ1TspkszBsWNLVqMTZipHKwZgxcP31MGGCkzLqEO6xS6Xi0kuj5/7EE07K6As5FSOVEydl1Aa2YqRy4qSMOohTMVIpaX71JSdllAdbMVIpmj0bJk1yUkafYytGKldjxzopo7zls8feE/g10A/4GDgZaP6bNwMYk3t9DTARWNbsZ9xjl1pzySWxjruTMsop1lTMhUBv4BrgWGAUcEGzn5kFHAG8v4HPMdil1qxeDRMnQv/+8JOfOCmjorVixgCP5e4/BhzQwmduD9wKzAZOzWMbkiAmZe68M66ZevPNSVejMtHaVMzpfH5v/C+sa6t8DPRt9nov4Gbgxtzn1wILgEUFVSpVqsY1ZUaNgsGDnZRRq1oL9ttyt6buA6pz96uBD5u9vpwI9hW5x08Du9JCsNfU1Ky9n8lkyGQybShZqkADBsC998akzMyZMHRo0hWpk2SzWbLZbLvek2+PvRq4GjgO2Bv4xyavDwXuAnYHNgKywBnAG80+xx671F4//zl85ztxEtMWWyRdjRJQrIOnPYFfAP2Bz4ApwLvAVKAOeIQI/2OBVbmfvbWFzzHYpXw4KVPRXCtGSiMnZSqaJyhJaeSkjFrhWjFSOWo6KTNkCBxySNIVqYTYipHKWeOaMk7KVAxbMVLajR0L//IvsabMe+8lXY1KhHvsUhp885uwYAE8/riTMinnVIxUKZyUqRi2YqRK4aSMmnAqRkoLJ2WUYytGShsnZVLNVoxUiRonZfbfH770Jdh8czjlFFixotW3qo1WroRzz431er785WiDlRBbMVIaDRsWl9T729/i8d13w5o18ItfJFtXWlx0Ufy7/PTTeHzGGXHger/9kq0rxz12KY1++9t1oQ6xt/7QQ8nVkzb3378u1CHuP/hgcvU0Y7BLadS3L3Tvvv5zq1bBm28mU0+avPMOLF++/nNdu8JmmyVTTwsMdimNTj0V+vWLcK+qgh494spLe+4Jp50G9fVJV1h+/vjH6KvvuisccED8O+3SJU4I23xzOO+8pCtcyx67lEabbgqvvAK33w4ffACHHRZjkO+/DzfdBHvtBYcfDtOnw6BBSVdb2t5+G777XbjnHjjrLPj97+MPzRdegIcfht694+B0v35JV7qW445SJfrgA5gxA374Qxg/PgJ+++2Trqq0vPUWXHcd3HdfBPpFF8GWWyZdleOOkr7AZptBTQ3U1cFXvhJ78yedBEuWJF1Z8v7wBzjzTNhjj9gLX7Ik9thLINTbymCXKtmmm8JVV0XADxoEo0fDiSdGu6HSvPkmnH46DB8OW28dgX7ttWV5bVmDXVIE/JVXxkHVIUPiJKfjj4fFi5OurPjq6+OA8l57xQldS5fGBcPLMNAbGeyS1unbF771rQi7YcNg771hyhR4442kK+t4dXUxPTRiBGy3XQT6NdfEhEuZM9glfV6fPnFA9c03YeedYd99YfJkeP31pCsr3NKlcPLJMHIkDBgQAX/11SU1h14og13SF6uuhssvjz34XXeNU+aPPRZefTXpytpvyZI4QDx6NAwcGIF+1VXRhkoZg11S66qr4bLLIuD32AO++lX42tdg0aKkK2vd4sVwwgkwZkyMdNbVxfGEFAZ6I4NdUtv17g2XXBIBv+eecOCBcPTRcTJUqVm8OA4A7713LF9cXw9XXBHHEVLOYJfUfr17x3VW6+tjBv7gg+Goo+Dll5OuLA70TpkC++wDO+4YNU6fHscNKoTBLil/m2wSZ2TW10er45BD4iIfCxd2fi2vvQbHHQeZDOyyS9Q0bVpFBXojg11S4Xr1ggsvjDDdZ59Ym2biRHjppeJv+9VX44Du/vvDbrtFDZddFscFKpTBLqnj9OoFF1wQ4brffrEOzRFHxIJZHW3RIjjmmFhpcfjw2Oall0abqMIZ7JI6Xs+e8E//FBMoBxwQ4T5hAixYUPhnv/JKHLA98MA4uai+Pvr9BvpahQT7kcAdX/DamcDzwLPAuAK2Iamc9ewJ558fAX/wwdGeGT8enn++/Z/18stxgPbgg2MW/c034eKLo8+v9eQb7DOA62h56chtgPOB0cDBwHeBjfPcjqQ06NEDvv71CPjDDosDrOPGwfz56//ce+/FtVqbLun90ktw5JFxoZCxY2MP/cILo+2jFuUb7HOAc2k52PfKvb4KWAbUAbvkuR1JadKjR1xpqK4u9tyPPjoCe/bs2JvfdttYiOuQQ2Du3HV7+JlMBPrUqQZ6G7QW7KcDi5rd9gDu2cB7qoGPmjz+GEj/GQGS2q5797jM3NKl0X8fPx4eeQRWrozbk0/CQQfFpEtdXfTre/ZMuuqy0Vqw3wbs3OzW2uHtZUS4N6oGPsi3QEkp1r07nHNOjCk2NKx7vqEBdtgBvvENAz0Pxbjm6XzgWqA70AMYCrS4YlBNTc3a+5lMhkwmU4RyJJW8oUOj9bJyZTzu2jXWhRfZbJZsNtuu9xRyzdN9gbOBKbnHU4l++iPAGcBZxN8IrgUeaOH9XvNUUnj//Rhd/MtfoKoq1nOZPx+22SbpykpOW6556sWsJZWGTz+Ng6gNDTH94hhjiwx2SUqZtgS7Z55KUsoY7JKUMga7JKWMwS5JKWOwS1LKGOySlDIGuySljMEuSSljsEtSyhjskpQyBrskpYzBLkkpY7BLUsoY7JKUMga7JKWMwS5JKWOwS1LKGOySlDIGuySljMEuSSljsEtSyhjskpQyBrskpYzBLkkpY7BLUsoY7JKUMga7JKWMwS5JKVNIsB8J3PEFr80AFgC1wNNAnwK2I0lqh3yDfQZwHVD1Ba/vDhwE7AfsDyzLcztlK5vNJl1CUaX5+6X5u4HfrxLkG+xzgHNpOdi7ANsDtwKzgVPz3EZZS/svV5q/X5q/G/j9KkFrwX46sKjZbQ/gng28pxdwM3A8cAhwHrBzwZVKktqkayuv35a7tcdyIthX5B4/DexK/KEgSSphGeCuFp4fCiwk/jbQjWjbDG3h5+qANd68efPmrV23OlrR2h77hjRupNHU3AYfAX4JPAusAn4OvNHC+wcVsG1JkiRJkiS1SRfgFmAucQLTwGTLKZoRxPdLk27Ar4BngHnAhGTL6XAbAT8jxnRnATsmW07RbAW8AwxOupAieJH4/66W9g9+lLrLidx8Hjg54Vo+ZxLxPw9E+D2YYC3FcgnwCvEfIU1OAW7M3d8MeDu5UoriCOCnufv7ks7fzW7AA8Bi0hfsPYhgT6MM8HDu/ibA1Rv64STWihkDPJa7Pw8YnkANxVZH/AH2RWfmlqt7gStz97sAf0uwlmJ4CDg7d38A8EFypRTN9cCPgT8nXUgR7EqcR/M48BSx45gWBxEj4w8SAyoPb/jHO9+txIlLjd4mnYuRDSAmg9Komjg/4bikCymSnwMfAQcmXEdHOwWYnrtfCwxJrpSi2Ik4qRLi7Pc60pMttxI7xF2Jv2ktTracz7sBOKbJ43eSKqTIBpDOYN+O6PGdknAdxbY18BbQM+E6OtJMIEuE+gfAc8T3TIuNiXZMo3nAlxKqpaN9F7iwyeOFwJYJ1dKiScDtufsjgUcTrKWYBpC+YN+aOCdhv6QLKZITiQNUECuSvgl0T66coqolfT32s4F/y93flvhdTcse+zjgidz9bYGllFirt4ro8c3J3dL2y9VoAOk7eDoD+BPrpg5qWX8Pqdz1BO4m9mznkr6pn6bSGOxdWTe19Qyx45gm3wPmE0uip61NKEmSJEmSJEmSJEmSJEmSJEmSJKlU/D/XNuTXv8dBWQAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Partly hierarchical basis: The level 1 basis function becomes level 2\n", "\n", "maxlevel = 3\n", "n=2**maxlevel - 1\n", "x = np.linspace(dmin, dmax, n+2)[1:-1]\n", "xExtend = np.linspace(dmin, dmax, n+2)\n", "y = np.sin(x)\n", "# 1. Gererate hierarchical intervals (for hat functions)\n", "\n", "# 1.1 Manually change interval3 (for h3)\n", "def gen_intervals_partial(a,b,max_level):\n", " intervals = []\n", " points = np.linspace(a,b,2**(max_level-1)+1)\n", " for i in range(2**(max_level-1)):\n", " intervals.append((points[i],points[i+1]))\n", " for i in range(2**(max_level-1)-1):\n", " intervals.append((points[i],points[i+2]))\n", " return intervals\n", "\n", "# 2. Construct A\n", "\n", "ints = sorted(gen_intervals_partial(dmin, dmax, maxlevel))\n", "\n", "A = np.array([[hat(l[0], l[1], i) for l in ints] for i in x])\n", "\n", "print(A)\n", "#multiply by matrix of next hierarchization step to get matrix from exercise 1\n", "B = np.identity(7)\n", "A = np.matrix(A)\n", "B[3,1]=0.5\n", "B[3,5]=0.5\n", "print(A*B.T)\n", "\n", "# 3. Solve for c\n", "\n", "cc= np.linalg.solve(A,y)\n", "\n", "# 4. Define interpolating function func3\n", "\n", "def func3(x):\n", " result = 0.0\n", " for i in range(len(cc)):\n", " result += cc[i]*hat(ints[i][0], ints[i][1], x)\n", " return result\n", "\n", "# 5. Plot\n", "y3 = list(map(lambda t: func2(t), x))\n", "\n", "#plt.plot(x,y) # original data seires\n", "plt.plot(x,y3, color='red') # compare with our 2nd approximation with n = 2^maxlevel basis functions\n", "plt.scatter(x,y3, color='red')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 2", "language": "python", "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.6" } }, "nbformat": 4, "nbformat_minor": 1 }