num0_ex02_DanielSchiller_ThomasBlankenburg_DennisJakob.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "| Übungszettel 2 | Einführung in die Numerik |\n",
    "|---|---|\n",
    "| Dozent      | Prof. Kanschat |\n",
    "| Tutoren     | S. Meggendorfer und J. Witte |\n",
    "| Abgabe      | bis 29.04.21 23.15 Uhr |\n",
    "| Studierende | Daniel Schiller, Thomas Blankenburg, Dennis Jakob |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "---\n",
    "\n",
    "### Aufgabe 1: Skalarprodukt\n",
    "\n",
    "Definieren Sie eine Funktion `l2prod(coeff_p, coeff_q)`, die das $L^2([-1,1])$-Skalarprodukt\n",
    "$$\\int_{-1}^1 p(x) q(x) \\,dx$$\n",
    "berechnet, wobei die Polynome $p$ und $q$ als Potenzreihe (aufsteigend) durch die Koeffizienten `coeff_p` und `coeff_q` bestimmt sind. Das heißt, bestehe der Vektor `coeff_p` aus den Skalaren $\\alpha_0,\\ldots,\\alpha_n$, so ist das Polynom $p$ bestimmt durch\n",
    "$$ p(x)=\\sum_{k=0}^n \\alpha_k x^k .$$\n",
    "\n",
    "**Verwenden Sie** entweder die importierte Potenzreihenklasse `PowSer` aus numpy, um das $L^2$-Skalarprodukt zu berechnen. **Oder** überlegen Sie sich, wie das $L^2$-Skalarprodukt endlicher Potenzreihen explizit berechnet wird."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "---\n",
    "\n",
    "Eine kurze Einführung zu den Polynomklassen in `numpy.polynomial` (insbesondere `PowSer`) finden sie [hier](https://numpy.org/doc/stable/reference/routines.polynomials.classes.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numpy.polynomial import Polynomial as PowSer\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Anbei finden Sie eine Funktion `polyplot(coeffs)`, die eine Liste von Koeffizientenvektoren `coeffs` entgegen nimmt und entsprechende Polynome plottet."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def polyplot(coeffs):\n",
    "    n = 1\n",
    "    for coeff_p in coeffs:\n",
    "        p = PowSer(coeff_p)\n",
    "        x, y = p.linspace(100)\n",
    "        plt.plot(x, y, lw=2, label=\"$p_{%d}$\"%n)\n",
    "        n += 1\n",
    "    plt.legend(loc=\"upper left\")\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Als Beispiel plotten wir alle Monome bis zum Grad 3:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD4CAYAAADhNOGaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABNGklEQVR4nO3dd3xT1fvA8c9Jmu5FaWlLSylbQUABAUURBGTIcONeIA4QRUBAUEBEQEQEAZWh4kb9oqIgKEtABGXKkl269x5pM87vjxv4VWihI2k6zvv16qs5ueshDXlyz73nOUJKiaIoilJ36ZwdgKIoiuJcKhEoiqLUcSoRKIqi1HEqESiKotRxKhEoiqLUcS7ODqAiAgMDZWRkpLPDUBRFqVH27t2bKqUMuvj5GpkIIiMj2bNnj7PDUBRFqVGEEOdKel51DSmKotRxKhEoiqLUcSoRKIqi1HE18hpBSUwmE7GxsRiNRmeHUm7u7u6Eh4djMBicHYqiKHVQrUkEsbGx+Pj4EBkZiRDC2eGUmZSStLQ0YmNjadKkibPDURSlDrJL15AQ4iMhRLIQ4nApy4UQYqEQ4pQQ4h8hRIdiy/oJIY7blk2saAxGo5H69evXqCQAIISgfv36NfJMRlGU2sFe1wg+AfpdZnl/oIXtZwTwPoAQQg8sti1vDTwghGhd0SDKkgSs0lrR3TtMTUteiqI4R2FepkP2a5euISnlNiFE5GVWGQJ8KrWa17uEEP5CiFAgEjglpTwDIIT42rbuUXvEdVGMHE6KQeiNWAuDgOr14ZuUUUD/iWudHYaiKNXUnXmruXv3Tta3bccbK1bZdd9VdY0gDIgp1o61PVfS811K2oEQYgTa2QQRERHlDsAqrQh9AUKY0RkysZrqlXsfiqIoVc2fHMa6fYLvvnPUz4ZmBf/a/RhVlQhK+votL/P8pU9KuRRYCtCpU6dyz6aj1+lpXq8xZ7POYtXnE+7rR4BHQHl34zDHcjyImn27s8NQFKW6kBKO/oB13RRmpwvujHbB6CYYunC93Q9VVeMIYoFGxdrhQPxlnncIdxd3Gno3BCAxP5F8U77dj3H//fczdOhQunTpQuPGjVm7VnX3KIpSTjmJsOph+PZxPnQppM1u7Tu737An8AwMtfvhquqMYA0wynYNoAuQJaVMEEKkAC2EEE2AOOB+4MHKHiyyzH3tZ8u137J8Yz948CB33HEHq1atYseOHbz00kvcfrv6pq8oShlICfs/h18ngzGLHT71+CvNi5fjrVj9fYkc9pxDDmuXRCCE+AroAQQKIWKBqYABQEr5AbAOGACcAvKBJ2zLzEKIUcAGQA98JKU8Yo+YnKGgoIDU1FSmTp0KQOvWrcnIyODMmTPMnDmTrKwsvvvuOydHqShKtZR+Fn56Ac7+DkBM855MIJHXfsgCIPS5Uei8vBxyaHvdNfTAFZZLYGQpy9ahJQq7Kcs3d5PFxJmsM5itZup71CfEK6TSxz18+DAtWrTA3d0dgH379tG+fXuaNm3KihUruOeeeyp9DEVRahmrBXZ/CJtngCkfPALIv20GL0Svpt2OHCJSwCU0BP/7hzoshDpba8igNxDuE45AkFaQRmZhZqX3efDgQaKjozEajeTl5TF16lTGjBlT+WAVRamdkv+Fj/rChklaErjmHuTIv5iW/Q9nU0/w8A7tIzpo1Ch0rq4OC6PWlJioCC+DF8FewSTmJRKfG4+b3g0PF48K7+/gwYM89NBD9OjRg+zsbF555RW6detmx4gVRakVzEXwx7uwbS5YisCnIQx8B1r159MjK/kl6hfuOOBCQEYhbi2a4zdkiEPDqbNnBOcFuAfg7+aPlJKY7BjMVnOF93Xw4EGefPJJdu3axdGjR3n44YcBSEtL45lnnmH//v3MmjXLXqErilITxe2FpT1gy0wtCXR8HEbuglb92Rm3k3f2voNXgeT+P/UABI0di3Bx7Hf2On1GAFp5h1DvUAothRSYC4jJiaGxb2N0ovw58vTp07Ro0eKS5+vXr88HH3xgj3AVRampivJh6yz4cxFIK9SLhEELoektAERnRzNu2zis0srU0+3Q5e7H8/rr8b7lFoeHVucTAYBO6Gjk04gzWWfIN+WTlJdEqHf579WNi4tzQHSKotR4UTtgzfOQfgaEDm4YBT0ng6snAHmmPEZvHk1OUQ4DPbvSZMNfSKDB+HFVUotMJQIbg95AI59GRGVHkW5Mx93FnXruqgyFoiiVYMyGjVNhz0daO+hqGLIYwjteWMUqrUzaPonTWadp6teUZ3b4kl9UhO+A/ni0a1clYapEUIynwZNQr1Dic+NJyEvAVe+Kl8Ex9+0qilLLndgAP4+B7DjQGaD7OLjpJXD5790/i/YvYkvMFnxcfZgfPJL8n19AGAwEvfhilYWqEsFF6rnXo9BSSFpBGjE5MTT1a4qr3nG3bSmKUsvkpcH6iXDoG60d1hEGL4LgSyvsrz2zlmWHlqEXeuZ1fxvduEUABDz2KK4VKK5ZUSoRlCDYM5hCSyG5RblE50TTxLcJep3e2WEpilKdSQmH/we/vAz5aeDiAbdOga7PQgmfH4dSDvHaH68BMKHzBNoczCRu/3709etT/5lnqjR0lQhKIIQg3DucM1lnKDQXEpsbS4RPhJpARlGUkmXHw88vwYlftHbkzTB4IQQ0LXH1hNwEnt/8PEXWIu5reR/3Nb6Ds88NBCDohdHovb2rKnJAjSMolV6nJ8I3Ar1OT25RLon5ic4OSVGU6sZqhT0fw+IuWhJw89VuCX3sp1KTQJ4pj5GbR5JmTKNLSBcmdplIxspPMcXH49aqFf53313F/wh1RnBZbno3Gvk04lz2OdIL0nHTuVWrOQwURXGitNNakbio7Vq71QC4fR74Nix1E4vVwsvbXuZkxkkifSOZ12MeJKWS+uGHAARPmojQV303tEoEV+Bl8KKhV0PicuMu3Enk7Vq1p22KolQjVgvsWgKbZ4K5ADwDYcBb0OYuuEL38dt73mZb7Db83PxY3Gsxfm5+xL41DVlQgE+/fnh17VpF/4j/UomgDPzd/Sm0FJJakEpMTgxN/Jrg7uJe4rr3338/UkqioqJITExkyZIlaj4CRaktko7CmlFamQiAdkOh7yzwqn/FTb849gWfH/scF50L7/Z4lwjfCPJ27SLnl/UIDw+CXx7v4OBLVzsTwTQ/u+8yGCgae5jswmzOZZ+jqV9TDHrDJeupiWkUpRYyF8H2edqP1QS+YTBwPrTsW6bNt0RvYc5fcwB4/cbX6RTSCWkykfjGGwAEPv00hoaldyk5Wu1MBA4S5h2G2Wom35RPdE40kb6R/7mttLSJaX744QfWrl1LcnIyI0eO5LbbbnPWP0FRlPKK3audBSQf1dodn4A+r4O7b5k2P5J6hAnbJyCRjLx2JIOaDQIg/fMvKDp1GkPjCAKefMJR0ZeJvWYo6wcsQJtlbLmUcvZFy8cDDxU75tVAkJQyXQgRBeQAFsAspexU6YCmZVV6FyXRAY18GnE26yxGs5HY3Fga+TS6UKCutIlp7rjjDu644w4yMjIYN26cSgSKUhMU5WsVQnct0YrEBTSFwe9B5E1l3kVsTiwjN42kwFzAkGZDeLrd0wCYEhNJfe89AIInTXLoXANlUelEIITQA4uBPmiT0f8thFgjpTx6fh0p5Vxgrm39QcAYKWV6sd30lFKmVjaWquCicyHCN4KzWWfJLcolIS+Bhl4NEUL8Z2Iai8XC1KlTeeutty5s+8YbbzByZIkTtSmKUp2c3aYVicuI0orE3Tgaer4ChrLPV5JpzOTZjc9qt4mGdmHqDVMvjEVKmjkTa34+Pn1649Ojh2P+DeVgjzOCzsApKeUZANsE9UOAo6Ws/wDwlR2O6zRuejcifCKIyo4i05iJQWeggWeDUiemkVIyceJE+vfvT4cOHZwdvqIopTFmwa+vwr6VWjv4Gu0sIKx8/2+NZiOjNo8iKjuKlvVaMr/H/AvXFHM2byHnt43oPD0JnjzZ3v+CCrFHIggDYoq1Y4EuJa0ohPAE+gGjij0tgV+FEBL4UEq5tJRtRwAjACKqsAZHaTwNnjTyaUR0djQp+SkYdAYOHjzIsmXLmDNnzn/Wfe+999i4cSNZWVmcOnWKZ6p4+LiiKGXw7zpY+xLkJIDeFbqPh24vXlIk7kosVgsTt0/kYMpBQrxCeL/3+/i4+gBgzc8n8Y0ZgDaC2BBS+bnS7cEeiaCkG2dlKesOAv64qFuom5QyXgjRAPhNCPGvlHLbJTvUEsRSgE6dOpW2/yrl4+pDqHcoCbkJxOfGc+r0qRInphk9ejSjR492QoSKolxRbopWH+jIaq0dfr1WJK7BVeXelZSSN3a/waboTfi4+vBB7w9o4NngwvKURYsxxyfg1vpq6j300GX2VLXskQhigUbF2uFAfCnr3s9F3UJSynjb72QhxPdoXU2XJILqKsA9ALPVTEp+Chv/2UiBpQAvnSpdrSjVnpRw6Fv4ZQIUpIPBE3q9Bp1HlFgkriyWHFzCdye+w03vxqJbF9HMv9mFZQWHDpP+yScgBKHTpzt8+snysEetob+BFkKIJkIIV7QP+zUXrySE8ANuAX4s9pyXEMLn/GPgNuCwHWKqUkEeQdRzr4eUkujsaIxmo7NDUhTlcrJi4cuhsPopLQk07QHP/VlqpdCy+Prfr/ng4AfohI653efSIfj/rytIk4mEKVPAaiXg0UfxaNvWTv8Q+6h0SpJSmoUQo4ANaLePfiSlPCKEeMa2/PxkvXcCv0op84ptHgx8b7uS7gJ8KaVcX9mYqpoQglCvUCxWC9lF2oCzJn5N1DwGilLdWK2w92P4bSoU5YC7H/R9E6596IrlIS5n/dn1vLn7TQCm3jCVnhE9/7M8bflyCo8fx9CoEUEvVL9uYrucm0gp1wHrLnrug4vanwCfXPTcGaC9PWJwNiEEYT5hWLIt5JnyiMqOoolvkxJHHyuK4gRpp2HNaDi3Q2tfNVArEudTuQu222O3M2n7JCSSFzq8wF0t7vrP8sJTp0hd8j4AoTNeR+fpWanjOUL16aSqBXRCd2HeY6PZyLmcc0T6RuKiUy+zojiNxQy7FsOWN8FsBK8gGDAXWt9RqbMAgH1J+3hp60uYpZnH2zzOsGuG/We5tFhImDwFaTLhf+89TisqdyXqE8rO9Do9jX0bczbrLIXmQqJzomns01jNcKYozpB4CH4cBQkHtHb7B6HvTPCsfDn5Y2nHGLVpFEaLkbta3MVLHV+6ZPKq9I8/puDgQVwaNKDBeOcVlbsSlQgcwEXnQqRvJGezzlJgKiAmJ4YI34gLpSgURXEwcyFsmws75oPVDH6NYNC70Ly3XXZ/OvM0T//2NDmmHPo07sNrXV+7JAkYj58gZcFCAEJnvoHet2y1iZxBJQIHMegNNPZrTFRWFHmmPGJyYv5Tl0hRFAeJ+Us7C0g9rrU7j9BuC3Xzsc/us2N46tenyCjM4Kawm5hz85xLzvhlURHxkyZqXUJDh+J98812ObajqETgQG56Nxr7NiYqO4rcolzicuMI9w5Xcx8riiMU5sLmN2D3B4CE+s21gWGNb7DbIRLzEnnqt6dIKUjh+pDr/1M6orjUDz6k8OgxDOHh1bpL6DyVCOyspIlpbu17K1FZUWQXZhNHHGHeYSoZKIo9nd4CP42GzGgQeug2Gm6ZCIaSJ5CqiKS8JJ7c8CRxuXG0C2zHe7e+V+IEVfn792tTTwpB6Jsz0XtX/wGmtTIRtF3pmMEahx47dMV1SpuYprFvY85lnyOrMAuBoKF3Q5UMFKWyCjLh18mw/3OtHdJWOwtoeK1dD5NakMrwX4cTkxPD1QFX836f9/EyXPoBb8nNJX78y2CxEPDkk3h17mzXOBylViYCZyltYppjx46xYMECEpMTaX9De+574j4QXChfrShKBRz7GdaOhdxE0LtBjwlauWg7j91JK0hj+IbhRGVH0apeK5bdtgxf15Iv/CbNeANTbCxura8m6MUX7BqHI9XKRFCWb+6OUNrENFdffTUffPABVquVJ4Y9gRCCTGMmAm1EskoGilIOucmwbjwc/UFrN+qinQUEtbT7odIK0hi2YRins07T3L85S29bip9byVPhZq1dS9aPPyLc3QmbO9fpk82Uh7qFxY6KT0yTl5fH1KlTGTNmDABr1qzhpptuom+fvkT4RCCEIMOYQUJeAlJWi2KqilK9SQkHvoJF12tJwOAF/efCE+sdkgRSC1IvJIFmfs1YfttyAtxLHn9QFBtL4rTpAARPnIBbs2YlrlddqURgR8Unprn++ut59tln6datGwCDBw9m586dfPHFF3i7ev8nGcTnxSNLrdytKAqZMfDFPfDDM2DMhGa9tCJxXUaAzv4fY6kFqQzfMPzCmcCKviuo71G/xHVlURFxY17CmpODd+9e+A8davd4HK1Wdg05S2kT02zdupXVq1dTWFjIgAEDAC4kg+icaDKNmWQZszBbzaochaIUZ7XCnhWwcRoU5YK7P/SbBe0fqHR5iNIk5SUx/FftmkBz/+Ysv215qUkAIHnePIyHDmFo2JCGM2fWyK5e9aljR6dPny5xYpoePXrQo4R5Sb1dvWns25jo7GgKzAVM3D6RWTfPwqBTheoUhdST2rzB0X9q7dZDtK4gn2CHHTI+N55hG4YRmxtLy3otWdpn6WWTQM6mTaSv/BRcXAib/w56v5KvH1R3KhHYUVxcXLm38TJ40di3MYkikQ1RGyiyFPH2LW+rEtZK3WUxwc73YOtssBSCVwOtSmjrwQ49bEx2DMN+HUZCXgKt67dmaZ/SLwwDFMXEED/pFQAajB2LR/uaW0hZXSOoBjwNngS4B+Dj6sOWmC08v/l58k35zg5LUapewj+w7FbYNF1LAtc+DCN3OzwJnMo4xWPrHyMhL4F2Qe1YdtuyyyYBa0EBsc+PxpqdjfettxLw+GMOjc/RVCKoJlz1rnzc92MC3APYGb+TZzY+Q05RjrPDUpSqYTLCptdhaQ9I/Af8IuDh1XDHYrtUCr2cI6lHeHzD4xfKRizts7TUcQKgzUucOG0ahf/+i6FxBA1nz6qR1wWKs0siEEL0E0IcF0KcEkJMLGF5DyFElhDigO3ntbJuW5e0CmjFJ/0+IdgzmP3J+xm2YRjpxnRnh6UojhW9Cz64CbbPA2mFLs9odwQ17+XwQ+9J3MOwX4eRVZhF9/DuLOm1pMQRw8VlfPElWT+uQXh4EP7ee9W6qmhZVToRCCH0wGKgP9AaeEAI0bqEVbdLKa+1/bxezm3rjCZ+Tfi0/6c08mnEsfRjPPbLY8Tnxjs7LEWxv8JcbWDYR/0g7SQEtoQn10P/OeDm7fDDb4newjMbnyHPlEe/yH682/PdEmsHFZf/998kzZ4NQMOZb+De0v7jF5zBHmcEnYFTUsozUsoi4GtgSBVsW2s19G7Ip/0/pVW9VkRlR/HIL49wOvO0s8NSFPs5tQmW3AB/LQWhg5vHwdPbIaJqZvD64dQPjNk6hkJLIfe1vI/ZN8++4t16RbGxxI5+AcxmAh5/HF/breC1gT0SQRgQU6wda3vuYjcIIQ4KIX4RQrQp57YIIUYIIfYIIfakpKTYIezqLdAjkI/6fUSHBh1Izk/msfWPcSD5gLPDUpTKyU+HH56Dz++CrGgIbQ8jtkKvV+1aKbQ0Uko+Pvwxr/7xKhZp4el2TzOl65QrziBoyc0j9tnnsGRk4HXzzTQYP87hsVYleySCkq6SXDxMdh/QWErZHngP+KEc22pPSrlUStlJStkpKCioorHWKL6uvnzY50N6NOpBVmEWw38dzubozc4OS1Eq5ugaWNwFDnyhFYnrPQ2Gb4bQdlVyeIvVwpy/5/DO3ncAmNh5IqOuG3XFC73SaiX+5ZcpPHkS16ZNCXtnHkJfu6aetUciiAUaFWuHA//p1JZSZkspc22P1wEGIURgWbat69xd3JnfYz53t7ibQkshY7aO4Zvj3zg7LEUpu5wkWPUIfPMI5CVDxI3w7E64aQzoq2YoU6GlkPHbxvPFsS8w6AzM7T6Xh65+qEzbJr81l9zNm9H5+dFoyWL0PvaZ6aw6scdf4W+ghRCiCRAH3A88WHwFIUQIkCSllEKIzmgJKA3IvNK2NU1JE9Pcfvvtldqni86FqTdMJdgrmCUHljBj1wzic+MZ3WG0mvpSqb6khANfwoZJYMwCV2/oMx06PumQ+kClyTRm8sKWF9iXvA9vgzcLei6gc2jZ5glI//wL0j/5BAwGwhe8i2tkpENjdZZKJwIppVkIMQrYAOiBj6SUR4QQz9iWfwDcAzwrhDADBcD9Uiu5WeK2lY3p2FVXV3YXJbr632NXXKe0iWkqSwjBs+2fJdgzmNf/fJ0Vh1cQlxvHGze9gZverdL7VxS7yjgHP70AZ7Zo7eZ9YOB88G90+e3sLDo7muc2Pce57HM08GzAkl5LaBXQqkzb5mzeTNKbbwIQOuN1vLpWzYVsZ7DLeZmtu2fdRc99UOzxImBRWbetqUqbmAYgLy+P7t27M336dAYOHFjhY9zV4i6CPYMZ+/tY1ketJyk/iQU9F1DPvZ5d/g2KUilWC/y9HDZOB1MeeARot4O2vddhReJKcyD5AKM3jyajMINW9VqxuNdigr3KVqeo4J9/iBs7DqxWAp8fhf8ddzg2WCerlbWGyvLN3RFKm5gGYM6cOdx33312OU63sG6s7LeSkZtGsj95Pw+ufZDFvRbT1L+pXfavKBWSclwrEhezW2u3uQv6vwXeVX9zx89nfua1P17DZDXRLawb826Zd8WBYucVnj5NzFMjkAUF+N11F4HPPefgaJ1PdTDbUWkT02zcuJHWrVsTHGy/qomtAlrxxYAvaF2/NbG5sTy87mF2xu202/4VpcwsJtg2VxsdHLMbvEPg/i/h3o+rPAlYpZWF+xYyafskTFYTQ1sNZdGti8qcBEzx8UQPG44lKwvvHj0InT6txpePKAuVCOyotIlptmzZwq5du/jyyy9ZtmwZVqvVLscL9grmk36f0KdxH3JMOTy36Tm+OPaFmvFMqTrxB2BpT9j8BliK4LpHtCJxV1X+ulh55ZvyGff7OJYdWoZO6JjUeRJTuk4p8xwf5vR0oocNx5yYiEfHjoS9Ox9hqBsl4Wtl15CzlDYxzcyZMwH45JNPCAwMRGfHOyY8XDx4+5a3WbR/EcsOLWP2X7M5nn6cKV2nqFLWiuOYCrQy0TvfA2kB/8YweCE07eGUcGJzYnlhywucyDiBt8Gbt295m25h3cq8vSUzk+gnh1F09ixuLVvS6P0l6NwdP8CtulCJwI5Km5jmvMcff9whx9UJHaM7jKZFvRa89sdrfH/qe85knWF+j/kEedaNwXdKFYr6Q7sWkH4aEND1Obh1CriWrfvF3v5K+Iuxv48lszCTSN9IFt66kCZ+Tcq8vSUnh+jhT1H477+4Nm5Mo+XLakUhufJQXUN2FBcXZ9dv++XVv0l/VvZfSYhXCAdTDjL056GqLIViP8Zs+Pkl+GSAlgSCroJhv2lTRzohCUgpWXlkJSN+G0FmYSY3h93MF7d/Ua4kYM3LI2bE0xgPH8YQHk7Eyk8wNGjgwKirJ5UIapnW9Vvz9e1f0ym4EykFKTyx/gm+PPalum6gVM7J37QicXtWgM4FbpkAT2+DRtc7JZx8Uz4vb3uZt/e8jUVaGHbNMN679b3LziNwMUtuLtFPjaBg/35cGoYS8cknGEJCHBh19aUSQS1U36M+S29byqOtH8Uszcz6axaTdkxSs54p5ZefDqufhi/ugexYaHgdjPgder4CLs4ZyHgm6wwPrXuI9VHr8TJ48W6Pd3mx44tXLBxXnCUri+gnh1Gwbx8uoaE0/uQTXMNLrHdZJ9SqRFBTv/U6Im6DzsD468czt/tcPFw8WHtmLQ+sfUCVs1bKRko48j0s7gz/fA0u7tBnBgzbCCHXOC2sdWfWcf/P93Mq8xRN/Jrw5e1f0qtx+SawMWdkEP3Ekxj/+QdDWBiNP/sM14gIB0VcM9SaRODu7k5aWlqNSwZSStLS0i4MQrO3fk368dXtX9HMrxlnss7wwNoH+On0Tw45llJLZCfAqofh28chLwUa36QVies2usqKxF2s0FLIG7veYML2CRSYCxjQZABf3/41Tf3KN4jSlJTEuUcewXj0KIaICBp//lmdPhM4T9S0D06ATp06yT179vznOZPJRGxsLEaj0UlRVZy7uzvh4eEYHHjPcr4pnxm7ZvDzmZ8BGNxsMJO7TMbT4OmwYyo1jJSw/zPYMAUKs8DVB26bAR0eq9IicRc7k3mG8dvGcyLjBAadgYmdJ3Jvy3vLPdCr8OxZoocNwxyfgGvzZkSsWIHBjoM8awIhxF4pZadLnq8tiUC5Mikl35/6nlm7Z2G0GIn0jWTuLXO5KuAqZ4emOFv6Wa1I3NnftXaLvlqROD/nfVuWUvLDqR+Y9dcsCswFRPhEMPeWubSuX/7ZbAsOHyFmxAgs6el4tG9Pow8/QO/vb/+gq7nSEkGt6RpSrkwIwV0t7uLrgV/T3L85UdlRPLj2QVYeWYlV2me0s1LDWC3w5xJ4/0YtCXgEwF3L4cFVTk0CWYVZjPt9HK/tfI0CcwEDmw7km0HfVCgJ5GzewrlHHsGSno7XzTcT8fFHdTIJXI46I6ijjGYjb+95m1XHVwHQJbQLM7vNLHN1RqUWSD4GP46CONv/pWvu1orEeQU6Nay/Ev5i0o5JJOcn4+niyeSukxncbHCF9pX++RdaKWmrFb8hQwid8TrCte6OuFddQ0qJtsZsZerOqaQb0/F19eWVLq8woMmAOlFoq84yF8GO+VqhOKsJfBrCwHegVX+nhmU0G1mwbwGfH/scgPZB7Zl10ywa+ZZ/DgNpNpM05y0yPvsMgMDnRxH43HN1/n3t0EQghOgHLECbXGa5lHL2RcsfAibYmrnAs1LKg7ZlUUAOYAHMJQV5MZUI7Cu1IJWpO6eyLXYbAH0a92FK1ykEuAc4OTLF7uL2wo/PQ7Jt/qeOT2izhrn7OTWsQymHeGXHK0RlR6EXep5u9zRPtXuqzAXjijNnZBD30kvk/7kLDAYavjEDvyFDHBB1zeOwRCCE0AMngD5ocxD/DTwgpTxabJ0bgWNSygwhRH9gmpSyi21ZFNBJSpla1mOqRGB/UkpWn1zNW3+/Rb45nwD3AF7p8gp9I/s6OzTFHoryYeub8OdikFao10QrEteku1PDKrQUsuTAElYeWYlFWmjq15Q3b3qTNoFtKrQ/44kTxI4chSkmBn39+oS/txDPDh3sHHXN5ciLxZ2BU1LKM1LKIuBr4D/pV0q5U0qZYWvuQpukXqlGhBDc3fJuVg9ZzfUh15NuTGfc7+MYs2UMqQVlztFKdXR2u3YxeOd7WvuGUdq4ACcngQPJB7hnzT18dPgjrNLKo60fZdXAVRVOAllr1hA19H5MMTG4t2lDk+++VUmgjOwxOiQMiCnWjgW6XGb9YcAvxdoS+FUIIYEPpZRLS9pICDECGAEQUcdHATpSmHcYy29bzncnvuOdve+wMXojuxN3M67TOO5sfmed72OtUYxZ8NtU2Pux1m7QGgYvgvCOTg0rtyiXBfsWsOr4KiSSpn5NmX7jdK5tcG2F9mctLCTpzVlkrtJufPAbMpiQ6dPrVBnpyrJH19C9QF8p5XBb+xGgs5Ty+RLW7QksAW6SUqbZnmsopYwXQjQAfgOel1Juu9wxVddQ1UjITWD6run8EfcHAB2DO/LaDa+VezSn4gQnNsDPYyA7DnQG6D4ebhoDLs67Y0ZKyaboTczaPYvkgmRchAtPXPMET7d/Gjd9xeoWFZ45Q9zYcRQeO4ZwdSV4ymT87y3/YLO6orSuIXucEcQCxS/rhwPxJQTQDlgO9D+fBACklPG238lCiO/RupoumwiUqhHqHcr7vd7nl7O/MOfvOexN2svda+7miTZPMLztcDUquTrKS4X1E+HQt1o7rKN2FhBc/vvv7SkmO4bZf8++cENCu6B2TL1hKi3rtazQ/qSUZK76hqTZs5FGI4ZGjQh7dz4ebSrWrVTX2eOMwAXtYnEvIA7tYvGDUsojxdaJADYDj0opdxZ73gvQSSlzbI9/A16XUq6/3DHVGUHVyyrMYv7e+fzv5P8ACPUKZcL1E7g14lb17as6kBIO/w9+eRny08DFA3q9Cl2egXJU5bQ3o9nIR4c/YsWhFRRZi/A2eDO6w2jua3lfuaqFFmdOSSFh2nRyN20CwG/IEIJfnYLe29ueoddKjr59dADwLtrtox9JKWcKIZ4BkFJ+IIRYDtwNnLNtYpZSdhJCNAW+tz3nAnwppZx5peOpROA8B5IP8ObuNzmWfgzQBqJNuH4CLeqVPjOb4mBZcbB2LJywXXpr0h0GLYSAsk/QYm9SSn499yvv7HmH+Dytg2BQ00G81OklAj0qNmBNSkn22nUkzZiBJSsLnbc3IdOm4Tew6udHrqnUgDLFbixWC9+c+IZF+xeRXZSNTui4t+W9PHftc2rsQVWyWmHfSvjtNSjMBjdfuO0N6PAoOPEs7WjaUd76+y32Ju0FoGW9lkzqPIlOIVccIlQqU1ISiTNmkLtROwvw6taN0DdmYAgNtUvMdYVKBIrdZRozWXxgMd+c+AartOJt8GZY22E8fPXDuLuoOzYcKu20ViQuarvWbjUAbp8Hvg2dFlJ8bjwL9y9k7Zm1APi7+fP8dc9zd4u7K9wNJC0WMr78ipR338Wal4fOy4sGEyfgf889qkuyAlQiUBzmZMZJ5u2dd+HuomDPYJ5t/yxDmg+p0MhQ5TIsZtj9PmyeCeYC8AyEAXOhzZ1OOwvIMGaw4tAKvvr3K4qsRRh0Bh686kGeavcUfm4VH7Gcv38/STPfxHj4MADevXsRMnmyOguoBJUIFIf7M/5P3tn7Dv+m/wtApG8kI68dyW2Rt6ETqtBtpSUd0YrExe/T2m3vg36zwau+U8LJLcrls6OfsfLoSvJMeQD0b9Kf0deNJtyn4mNGTUlJJM+bR/YabQIll5AQQl6dgk+v8s1EplxKJQKlSlillfVn17PowCJicrRxhs39m/Ns+2fp3bi3SggVYS6E7e/A9nlakTjfMBj4LrS8zSnh5Bbl8uW/X7LyyEqyi7IB6BbWjeeve5429St++6YlN5e0FStI/2QlsqAA4epKwJNPEPjUU+i8vOwVfp2mEoFSpUxWE9+f/J6l/ywlKT8J0BLC8LbD6RvZV3UZlVXsHu0sIEW7S4tOT0Lv6eDuW+WhZBVm8eW/X/L50c8vJIAODTrw/HXPV+pCsNVoJHPVKlLf/wBLZiYAPn360GDCy7iGq2o09qQSgeIURZYifjj1A8sOLSMxLxGAcO9wnrjmCQY3G6wuKpemKE+7DrBrCSAhoCkMfg8ib6ryUJLykvj06Kd8e+JbCswFgJYAnrv2OTqHdK7wRVur0UjmN9+Qtmw55pQUADw6dqTB2LF4drjObvEr/08lAsWpiixFrDm9ho8Pf0x0TjQAAe4BDG01lKGthlLfwzn93NXSma2wZjRkngOhgxufhx6TwOBRpWEcSzvGp0c/Zf3Z9ZilGYAbQm9geNvhXB9yfYUTgCU7m4yvV5H+2adYUrSChm6tryZo1PN49+yh7gZyIJUIlGrBYrXwW/RvfHToowuD0lx1rvRv0p8HrnqgwpUna4WCTPjtVdj3qdYOvkY7CwirugqaJquJLdFb+Orfr9iTpP0f0wkdvSJ6MaztsEpdAyiKjSXjs8/J/PZbrPn5ALi3bk3gqJF49+ypEkAVUIlAqVaklOxJ2sOnRz/l95jfkWjvw3aB7bi31b30jeyLh0vVfgN2qn/XwdqXICcB9K7Q/WW46UXQG6rk8Il5iXx/6nu+O/4dyQXJAHgZvLirxV08eNWDFb4LSFqt5P3xBxlffEnu779rpTAAzxu6Un/YcLy63agSQBVSiUCptqKzo1l1fBXfn/qenKIcALwN3tze9HbubH4nreu3rr0fFrkpWn2gI6u1dvj1WpG4Blc5/NAmi4ltcdtYfXI1O+J2YJVWQLvt94GrHmBws8F4u1asfo8pLo7M738ga/VqTPFaiQlhMOA7oD/1Hn1UFYdzEpUIlGqvwFzAL2d/4X8n/sc/qf9ceL6ZXzMGNRvE7U1vJ8QrxIkR2pGUWoXQXyZAQToYPKHXa9B5hEOLxEkpOZR6iJ9O/8T6qPVkFmYC4KJzoVdEL+5ucTddQ7tWKPFaMjPJ3vAr2WvXkv/33xe+/RsaNsR/6FD8770HlwBVgsSZVCJQapQTGSdYfXI1686sI6Mw48Lz1zW4jr6Rfbmt8W0EeQY5McJKyIrV5go4+avWbtoDBi2AepEOOZyUkuMZx1l/dj3ro9YTlxt3YVlz/+bc0fwOBjUbVKE6UebUVHI2bybnt43k7doFJhMAwtUVnz598L/7Ljy7dkXo1PiR6kAlAqVGMllN/BH3Bz+d/onfY3+n0FIIgEDQLqgdt0bcSs9GPYn0jaz+3UdWK+z9CH6bBkU52oTxfd+Eax+ye3kIs9XMgeQDbI7ZzObozf/58A/yCKJvZF+GNB9Cq3qtyvW6SasV49Fj5G77nbxt2yk4ePDCN390OrxuuAHfgQPx6d0LvY+PXf9NSuWpRKDUePmmfLbGbGV91Hp2xu+8kBRAG5twc/jN3BR2Ex2DO+JlqGYjUdNOw5rn4ZxWj4mrBmpF4nzs19WVmJfIroRd7Ijbwc74nReutwDUd69P78a96RvZlw4NOpS5CJy0Wik6fZr8PXvI27Wb/N27Lwz6Aq3f3+vGG/Hp0xvvW29VXT/VnEoESq2Sb8pnZ/xONkVvYkfcjgt93QAuwoW2QW3pHNKZDsEduDboWufNpmYxw5+LYOssMBvBKwgGvA2th1T6LCApL4n9yfvZk7SH3Qm7icqO+s/ySN9IejTqQa+IXrQNbFumD39zejrGI0coOHSIgoMHKThwEGtW1n/WcQkNxbt7d7y734xnl67ovatZ0lVKpRKBUmtZrBYOpx1mW+w2dsXv4nDa4Qt3wADohZ6W9VrSNrAt1wReQ5vANjTxa4JB5+BbMxMPaeUhEg5o7fYPQt+Z4Fn+b815pjyOpx/ncOphDqce5p/Uf/7T3QPg6eJJp5BOdGvYjZvDbqaRb6NS9gbWvDwKo6IoOnOWwpMnKTx+HOOJE5gTEi5Z1yUkBM8OHfDs2gWvLl0wRERU/244pUSOnqGsH7AAbYay5VLK2RctF7blA4B84HEp5b6ybFsSlQiUy8kpymFP4h72JO1hf/J+jqUduzAy9jyDzkBz/+a0rNeSpv5NaerXlEjfSMK8wzBU9t59cyFsmws75oPVDH6NtCJxLXpfcdPcolxicmI4k3WGs1lnOZ15muMZxy8U8CvO2+BN+wbt6dCgA51DOtMmsM2F5GbJzcOckow5KQlTfAKmxARMMbGYYmIoionBnJRU4vGFpyfura/Go00b3K9pi2fHDhgaOm+OA8W+HJYIhBB6tDmL+6BNZP838ICU8mixdQYAz6Mlgi7AAilll7JsWxKVCJTyyDflcyTtCEdSj3Ao9RBH0o5c8m36PJ3Q0cCzAWHeYQR7BhPsGUwDzwYEuAfg7+6Pv5s/Pq4+eBm88DJ44apz/e+345i/tLOA1OMAyOufouCW8eTpdOSZ8sguyiazMJMMYwapBakk5yeTnJ9MYm4CCVmx5BZk4mIBFwu4msFgAVcTeFr0RLqG0sI1jEhDCOG6AAJMbsicXCzZ2VgyMrBkZGDOyMCcmoq0jdwtlcGAW2RjXJs0xbVZU9xbXYVbq5a4RkQg9M6b41hxrNISgT1KQHYGTkkpz9gO9DUwBCj+YT4E+FRqWWeXEMJfCBEKRJZhW0UpM2kyYc7IwJKRiSUzE2tONpasbJrl5hCZl0//vBCs+b4U5eWQnZ1Kfm46xoJcTAV5mAuNWM0m9JZY9NZY9FbQSdDZfiMhC8iWIIp/fxIgpNBWuMD2LfqHtQjW/v+qErwk+EhoLLX9nj/O5VnQpvzWpv2WQNpl1hbu7rgEBmIICcElNBRDSAiG8HBcG4VjaNQIQ2gowkVVgFU09ngnhAHFz1tj0b71X2mdsDJuC4AQYgQwAiAiIqJyESs1krWoCFNcHKbYWExx8ZgSEjAnJmBKSsacmoI5JfWSC5uX4wL42n4qzw7X2lz0CBcDwmBAuLoi3FzRGVwRnp7o3N3Rebij8/JC5+ml/fb1Qe/rh97XB329etqPfz1cggLReXurfnylzOyRCEp6t138v6K0dcqyrfaklEuBpaB1DZUnQKVmseTmUnjiJIUnTlB45jRFZ85SdOYMpoSE/79nvTQ6ne1D0R+9v7/2Qenjg87XF52Xp+2D1BOdhyc6D3eEmxvCzQ2dq6v24WswgIsLwsVF6yLRuyB0AvR6tK/+aIOjhABjFpYtb2I+spoiATKoJfSaCg1aoxd6XPUGXHQGhBD//6Es/n9fQie0Y+n1oNOpD27FaeyRCGKB4rcnhAPxZVzHtQzbKrWYNS+PgkOHMB4+TMHhIxgPH8YUG1vyyjodhoYNMUQ00n6HNsQQGopLcDAuQUG4BAWi9/evmlGsx36GtWNxyU3EzdMNrx4T4MbRVVYkTlHsyR6J4G+ghRCiCRAH3A88eNE6a4BRtmsAXYAsKWWCECKlDNsqtYglM5O8v/8m/6+/Kdi7F+Px42Cx/GcdYTDg2qI57i1a4taiOa5Nm+Ia2QTXRuHaN3Znyk2GdePh6A9au1FXrVR0UEunhqUolVHpRCClNAshRgEb0G4B/UhKeUQI8Yxt+QfAOrQ7hk6h3T76xOW2rWxMSvUhzWYKDhwg9/dt5P3xB8Zjx/7bvaPX496mDR7t2+He5hrcr2mDW9Omzv/Av5iU8M8qWD8RCjLA1Rt6T4NOw0DV0VFqODWgTLE7a14eudt3kLNxI7nbtmHNzr6wTBgMeFx3HZ5dOuPZsRMe7dqi83TSqN+yyozWisSd2qi1m/WCQe+Cv7ppQalZHHn7qKJgLSwkd+vvZK9dS+7vvyML/78OkGuTJnh3vxmvm7vj2akjOvcaMk+x1Qp7VsDGaVCUC+7+0G82tL/f7kXiFMWZVCJQKkxKScGBA2T+73/krN+ANTf3wjKPa6/Fp09vfHr3xrVxYydGWUGpJ7UicdF/au3WQ7QaQd4NnBuXojiASgRKuVmyssj8/nsyv/2OotOnLzzv3ro1vgMH4jugP4aQGjqBjMUEOxfC1jlgKQTvYK1K6NWDnB2ZojiMSgRKmRmPnyDj88/J+uknpNEIgD4wEL8hg/G/6y7cmjVzcoSVlHAQfhypFYsDuPZh6PsGeNRzblyK4mAqESiXJaUkf/du0lZ8RN727Ree97rxBvwfeACfHj2q3x0+5WUywu9z4I8FIC3aReBBC6DZrc6OTFGqhEoESomklORu3UrqkvcxHtK+IQsPD/zvvJN6Dz+EW9OmTo7QTs79CWtGQdopQECXZ+DWV8GtYpO2K0pNpBKB8h9SSnK3bCV18WKMR7QhHfqAAAIeeRj/++/HpV4t6SYpzIGN0+Hv5YCEwJYweBFElFjqSlFqNZUIlAvy9+0j+e15FOzbB4A+KJDA4cPxv+8+dB4eTo7Ojk5uhJ9fhKwY0LlAtxeh+3gw1JDbWhXFzlQiUCiKiiJp7tvkbtoEgL5ePeo/PYJ6999fc+75L4v8dNjwChz8SmuHXgtDFkFIW6eGpSjOphJBHWbJzSXtgw9IW/kpmEwIT0/qP/4YAU8+id67lvWRH/kB1o2DvBTQu0HPV+CGUaBX/wUURf0vqIOklOT88guJs2ZhSUkFwO+uuwh68QUMDWrZgKmcRC0BHPtJa0fcqBWJC2zu3LgUpRpRiaCOKYqJIXH66+Tt2AGAR/v2BE+ZjEfbWtY9IiUc+ELrCjJmaUXi+kyHjk+qInGKchGVCOoIabGQ8fnnJM9/F2k0ovP1pcH4cfjffXfV1O+vShnn4KcX4MwWrd28j1Ykzi/cqWEpSnWlEkEdUHTuHPGvTKZg714AfAcOJHjiBFwCA50cmZ1ZLfDXMtg0HUz52ojg/m9B23tVkThFuQyVCGoxKSWZq1aRNHsO0mhEHxRI6PTX8bm1p7NDs7+U4/DjKIj9S2u3uUtLAt5Bzo1LUWoAlQhqKXNGBglTXr1wS6jvwIGETJmM3t/fuYHZm8UEf7wLv78FliLwDoGB78BVtzs7MkWpMSqVCIQQAcAqIBKIAu6TUmZctE4j4FMgBLACS6WUC2zLpgFPASm21V+RUq6rTEwK5P31F/Fjx2FOSUHn40PItKn43V4LPxjj98OPz0OSrUhch8egz+vg4e/UsBSlpqnsGcFEYJOUcrYQYqKtPeGidczAWCnlPiGED7BXCPGblPKobfl8KeXblYxDAaTVStqy5aQsWABWKx4dOxL21hwMYWHODs2+TAWwdRbsfA+kFfwba7eENr3F2ZEpSo1U2UQwBOhhe7wS2MpFiUBKmQAk2B7nCCGOAWHAURS7sWRlEf/yBHJ//x2A+k8/TdDzoxAutaz3L+oPbcKY9NMgdNqgsJ6vgKuXsyNTlBqrsp8SwbYPeqSUCUKIy45GEkJEAtcBu4s9PUoI8SiwB+3MIaOUbUcAIwAiItRcscUVnjpFzMiRmM5Fo/Pzo+Gc2fj06OHssOzLmK1NGblnhdYOugqGLIbwS6ZfVRSlnK6YCIQQG9H69y82uTwHEkJ4A/8DXpRSnp/N/H1gBiBtv+cBT5a0vZRyKbAUtMnry3Ps2ixn8xbix4/HmpeH29VXE/7ee7iG17KuoBO/apPHZ8dqReJuHqv9uLg5OzJFqRWumAiklL1LWyaESBJChNrOBkKB5FLWM6AlgS+klKuL7Tup2DrLgJ/LE3xdJqUkbflyUt6ZD1LiO6A/oTNn1q4qoXlpsGES/LNKaze8TjsLCG7j3LgUpZapbNfQGuAxYLbt948XryCEEMAK4JiU8p2LloWe71oC7gQOVzKeOkGaTCS+PoPMb78FIQgaM4b6I55C1JZBU1LCke9h3XjITwUXd7h1CnR5VhWJUxQHqOz/qtnAN0KIYUA0cC+AEKIhsFxKOQDoBjwCHBJCHLBtd/420beEENeidQ1FAU9XMp5az5KbS9yLY8jbsQPh5kbDt97Ct+9tzg7LfrITYO1YOL5Waze+CQYvhPo1fD5kRanGKpUIpJRpQK8Sno8HBtge7wBK/KoqpXykMseva8xpaUQ/9RSFR4+hDwig0ZLFeFx7rbPDsg8pYf9nsGEKFGaBm682JqDDY6pInKI4mDrPriGKYuOIGTaMonPnMDSOIGLZMlxry91T6Wfhp9FwdpvWbtkPbn8H/GrZRW9FqaZUIqgBCk+eJHrYcMzJybhdfTURy5bWjoJxVgvs/gA2v6EVifOsr9UHuuZuVSROUaqQSgTVnPHoUaKfHIYlMxPP668nfMli9D4+zg6r8pKPaUXi4vZo7Wvugf5zwKsWJDhFqWFUIqjGCg4dInrYcKzZ2Xjd0p3wBQtq/hzC5iLYMR+2zQWrCXwawsD50KqfsyNTlDpLJYJqKn//fmKeGoE1Nxfv3r0If+cdhKurs8OqnLi92llAsq26SMcntFnD3P2cG5ei1HEqEVRDBf/8Q8zwp7Dm5eHTvx9hb72FMBicHVbFFeXD1jfhz8Vakbh6TbQicU1udnZkiqKgEkG1Yzx6lGhbEvAd0J+Gb71VswvHnd2uFYnLOFusSNxkcPV0dmSKotjU4E+Y2sd4/ATRTw7Dmp2NT5/eNJwzp+YmAWMW/PYa7P1EazdoDYMXQXhHp4alKMqlauinTO1TFB1N9HDt7iDvW24hbN68mtsddPwX+PklyIkHnQG6j4ebxoBLDb/GoSi1lEoE1YApOZnoYcOxpKTi2bUrYQsX1MwLw3mp8MsEOPyd1g7rBEMWQYOrnRuXoiiXpRKBk1mys4l5agSmmBjcr7mG8EWL0LnVsPLKUsLh/8EvL0N+Grh4QK9XocszoNM7OzpFUa5AJQInshYWEvvcSAqPH8e1SRMaLf0QvXcNm2krK04rEnfiF63dpDsMWggBTZwbl6IoZaYSgZNIq5WESa+Qv2cPLsHBRKxYjktAgLPDKjurFfat1C4IF2aDmx/cNgM6PKrKQyhKDaMSgZOkzH+X7HXr0Hl50WjphxgaNnR2SGWXdhp+egGitmvtVrfD7fPAN9S5cSmKUiEqEThBxqpvSFu2DPR6whYswL1VK2eHVDYWM+xaAltmgtkInoEwYC60uVOdBShKDVapRCCECABWAZFoE8vcV9Lk80KIKCAHsABmKWWn8mxfm+T9+SeJr78OQOj0aXjf1M3JEZVR0hGtPET8Pq3dbij0mw2eNag7S1GUElV2xo+JwCYpZQtgk61dmp5SymvPJ4EKbF/jFUVFEfviGLBYqP/UcPzvucfZIV2ZuRC2vAkfdteSgG8YPPgt3LVUJQFFqSUq2zU0BOhhe7wS2ApMqMLtawxLdjYxzz6HNSsL7549CRozxtkhXVnsHu0sIOWY1r5+OPSaCu6+zo1LURS7qmwiCD4/+byUMkEI0aCU9STwqxBCAh9KKZeWc3uEECOAEQARNWxmLmmxEDd2HEVnz+LWogUN585FVOfpF4vyYPNM7XoAEgKaaUXiImtIN5aiKOVyxUQghNgIhJSwaHI5jtNNShlv+6D/TQjxr5RyWzm2x5Y8lgJ06tRJlmdbZ0tZ+B5527ej9/cn/P0l1XuswJnftWkjM6JA6OHG56HHRDB4ODsyRVEc5IqJQErZu7RlQogkIUSo7dt8KJBcyj7ibb+ThRDfA52BbUCZtq/JcjZtIu3DD0GnI+zd+biGhzs7pJIVZMJvr8K+T7V28DVaeYiG1zk1LEVRHK+y/RNrgMdsjx8Dfrx4BSGElxDC5/xj4DbgcFm3r8kKz5wl/mXtkkeDsS/h1bWrkyMqxb9rYXEXLQnoXeHWKTBiq0oCilJHVPYawWzgGyHEMCAauBdACNEQWC6lHAAEA98L7T5zF+BLKeX6y21fG1jz8ogd/bw2uUzfvgQ8+aSzQ7pUbopWH+jIaq0d3lk7CwiqIeMaFEWxi0olAillGtCrhOfjgQG2x2eA9uXZvjZIfH0GRadO49qsGaEzZyKq04ArKeHQt1ql0IJ0MHhB76naXUGqSJyi1DlqZLEDZK7+nqwff0S4uxP+7vzqdXE4KxZ+HgMnf9XaTXvCoAVQr7Fz41IUxWlUIrCzwpMnL4wcDnntNdxatHByRDZWK+z9CH6bBkU52oTxfWfBtQ+q8hCKUsepRGBH1oICYseMQRqN+A0Zgv9ddzo7JE3qKW3e4OidWvvqQTBgHvgEOzcuRVGqBZUI7Chp9hztukDTpoRMfc3Z4WhF4v5cBFtnaUXivBrA7W9D6yHOjkxRlGpEJQI7ydm4kcxVqxAGA2Hz30Hn6encgBIPwY8jIeGg1m7/IPSdqeoDKYpyCZUI7MCUlEzClFcBaDBurHPLSpuMsG0u/PEuWM3gFwGD5kPzUscFKopSx6lEUEnaTGMTsWRm4tWtG/UeecR5wUTvhjWjIPUEIKDzCK1InJu382JSFKXaU4mgkjI+/5y8nX+ir1eP0FlvOqeYXGEubJ4Buz8EJNRvoQ0Mi6imI5kVRalWVCKohMIzZ0me9w4AoTNex9Cg1OKpjnNqE/z0ImRFa0Xiur0At0wAg3vVx6IoSo2kEkEFSbOZ+EkTkYWF+A0Zgk/vKu6DL8iADVPgwOdaO6QtDFkMoSUO4lYURSmVSgQVlLbiI4wH/8ElOJjgya9U7cGProF14yA3CfRu0GMC3Dga9IaqjUNRlFpBJYIKMB4/QcqiRQCEzpyJ3reKZuzKSdISwLE1WrtRV23CmKCWVXN8RVFqJZUIykmazSRMngwmE/5Dh1bN5PNSwsGvYP0kMGbaisRNsxWJq8YznSmKUiOoRFBO6Ss/xXj4MC6hoTQYP87xB8yM1i4Gn96ktZv1gkHvgn/Nmq5TUZTqSyWCciiKiiJl4UIAQqdPQ+/twPvzrVb4ezlsnAamPHD3h36zoP0Dqkicoih2pRJBGUmrlYRXX7PdJTQY7+7dHXew1JPw4yiI2aW1W98BA+aCtxNuT1UUpdarVAezECJACPGbEOKk7Xe9EtZpJYQ4UOwnWwjxom3ZNCFEXLFlAyoTjyNlfvsd+X//jb5+fRpMnOiYg1hMsH0evN9NSwLewTD0c7hvpUoCiqI4TGXPCCYCm6SUs4UQE23tCcVXkFIeB64FEELogTjg+2KrzJdSvl3JOBzKnJpK8rx5AIRMfgWXepfku8pLOKgViUs8pLWvexhuewM8HHAsRVGUYiqbCIYAPWyPVwJbuSgRXKQXcFpKea6Sx61SSbPnYM3Oxuvmm/Hp39++OzcZ4ffZ8MdCkBbtIvCghdCsp32PoyiKUorKJoJgKWUCgJQyQQhxpf6L+4GvLnpulBDiUWAPMFZKmVHShkKIEcAIgIiIqrtjJvePP8j++WeEmxshr71q37mHz/2pTRiTdhIQ0PU5uHUKuFajqS0VRan1rniNQAixUQhxuISfcs1uIoRwBQYD3xZ7+n2gGVrXUQIwr7TtpZRLpZSdpJSdgoKCynPoCrMWFl6YdjLwuedwbdTIPjsuzIG14+DjfloSCGwFw37V7gpSSUBRlCp2xTMCKWWpRXSEEElCiFDb2UAokHyZXfUH9kkpk4rt+8JjIcQy4OeyhV010pYuw3QuGtfmzaj/xOP22enJjfDzi5AVAzoXuGkMdB8PLm722b+iKEo5VbZraA3wGDDb9vvHy6z7ABd1C51PIrbmncDhSsZjN0UxMaQtWwZA6NSpCFfXyu0wPx02vKKNEAatONyQxVqxOEVRFCeqbCKYDXwjhBgGRAP3AgghGgLLpZQDbG1PoA/w9EXbvyWEuBaQQFQJy50mafYcZFERvoMG4Xn99RXfkZRw9EetRlBeCri4Q49JcMMo0KthHIqiOF+lPomklGlodwJd/Hw8MKBYOx+oX8J6TpzOq3S527aRu2kTOk9PGoyrRBmJnERYOxb+tfV4RdyoFYkLbG6fQBVFUexAfSW9iLWoiKSZbwIQOGoUhuAKDOSSEg58oXUFGbPA1Qf6TIeOT6gicYqiVDsqEVwkfeVKis6dw7VZMwIeebj8O8iIgp9egDNbtXbzPlqROL9wO0apKIpiPyoRFGNOSSHt/Q8AbQSxMJRjoherBf5aBpte14rEeQRA/znQ9l5VJE5RlGpNJYJiUhYuxJqfj/ett+J1443l2PC4ViQu9i+t3eYu6P8WeFfNeAdFUZTKUInAxnjsGJnf/Q9cXMo+z4DFBDvehW1vgaUIfELh9nfgqmpbO09RFOUSKhEAUkqS3pwFUhLw0EO4NWly5Y3i92tnAUm2oQ8dHoU+M8DD36GxKoqi2JtKBEDOxo1aiWl/fwKfe/byK5sKYOss2LlIKxJXL1IrEtf0liqJVVEUxd7qfCKQJhPJb2tVsANHjULv51f6ylE7YM1oSD8NQqcNCus5GVw9qyhaRVEU+6vziSDjm2+0ekKRkdQbel/JKxmztSkj96zQ2kFXw5BFEN6pyuJUFEVxlDqdCCy5uaQuXgJA0NiXSr5d9MSvWpG47DitSNzNY7UfVSROUZRaok4ngvSPPsKSno7Hddfh0/uiIqt5abB+Ihz6Rms3vE4rEhfcpuoDVRRFcaA6mwhMScmkffwJAA3Gj/v/CWekhCOrYd3LkJ+qFYnrOVmbNEYViVMUpRaqs59sqYsWIQsK8OnTG88OHbQnsxNg7UtwfJ3WjrwZBi2A+s2cF6iiKIqD1clEUBQVRebq1aDXEzRmjHYWsO9T+PVVKLQVibttBnR4TBWJUxSl1quTiSDlvUVgseB3z924+Qv4dDCc3aYtbNlPGx3sF+bcIBVFUapInUsExn//JXvtWoTBQFBXb1hyA5gLwLO+Vh/omrtVkThFUeqUSvV7CCHuFUIcEUJYhRCl3lQvhOgnhDguhDglhJhY7PkAIcRvQoiTtt/1KhNPWaQsWAiAf1tXDHtmaUmg7b0w8i9oe49KAoqi1DmV7QA/DNwFbCttBSGEHliMNnl9a+ABIURr2+KJwCYpZQtgk63tMPl7/iZ3yxaEiyQw4jT4NIQHVsHdy8Er0JGHVhRFqbYqlQiklMeklMevsFpn4JSU8oyUsgj4GhhiWzYEWGl7vBK4ozLxXI4syCJlwpMABLTMxaXbYzByF7Tq56hDKoqi1AhVcY0gDIgp1o4FutgeB0spEwCklAlCiFLnhRRCjABGAERERJQ7CCkN6Hx80LllUH/6Cmjbt9z7UBRFqY2umAiEEBuBkBIWTZZS/liGY5TU6S7LsN1/N5ByKbAUoFOnTuXeXufpSaOv1mHOzEUfqqaNVBRFOe+KiUBK2ftK61xBLNCoWDsciLc9ThJChNrOBkKB5Eoe6/I8/HFR8wUoiqL8R1WMlvobaCGEaCKEcAXuB9bYlq0BHrM9fgwoyxmGoiiKYkeVvX30TiFELHADsFYIscH2fEMhxDoAKaUZGAVsAI4B30gpj9h2MRvoI4Q4CfSxtRVFUZQqJKQsd3e703Xq1Enu2bPH2WEoiqLUKEKIvVLKS8Z8qUI6iqIodZxKBIqiKHWcSgSKoih1nEoEiqIodVyNvFgshEgBzlVw80Ag1Y7h2IuKq3xUXOWj4iqf6hoXVC62xlLKoIufrJGJoDKEEHtKumrubCqu8lFxlY+Kq3yqa1zgmNhU15CiKEodpxKBoihKHVcXE8FSZwdQChVX+ai4ykfFVT7VNS5wQGx17hqBoiiK8l918YxAURRFKUYlAkVRlDquViYCIcS9QogjQgirEKLU26yEEP2EEMeFEKeEEBOLPR8ghPhNCHHS9rueneK64n6FEK2EEAeK/WQLIV60LZsmhIgrtmxAVcVlWy9KCHHIduw95d3eEXEJIRoJIbYIIY7Z/uYvFFtm19ertPdLseVCCLHQtvwfIUSHsm7r4LgessXzjxBipxCifbFlJf5NqyiuHkKIrGJ/n9fKuq2D4xpfLKbDQgiLECLAtswhr5cQ4iMhRLIQ4nApyx373pJS1rof4GqgFbAV6FTKOnrgNNAUcAUOAq1ty94CJtoeTwTm2Cmucu3XFmMi2iAQgGnAOAe8XmWKC4gCAiv777JnXEAo0MH22Ac4UezvaLfX63Lvl2LrDAB+QZuVryuwu6zbOjiuG4F6tsf9z8d1ub9pFcXVA/i5Its6Mq6L1h8EbK6C16s70AE4XMpyh763auUZgZTymJTy+BVW6wycklKekVIWAV8DQ2zLhgArbY9XAnfYKbTy7rcXcFpKWdFR1GVV2X+v014vKWWClHKf7XEO2pwXYXY6fnGXe78Uj/dTqdkF+Att5r2ybOuwuKSUO6WUGbbmLrRZAh2tMv9mp75eF3kA+MpOxy6VlHIbkH6ZVRz63qqViaCMwoCYYu1Y/v8DJFhKmQDaBw3QwE7HLO9+7+fSN+Eo26nhR/bqgilHXBL4VQixVwgxogLbOyouAIQQkcB1wO5iT9vr9brc++VK65RlW0fGVdwwtG+W55X2N62quG4QQhwUQvwihGhTzm0dGRdCCE+gH/C/Yk876vW6Eoe+t644Z3F1JYTYCISUsGiylLIsU16KEp6r9L20l4urnPtxBQYDk4o9/T4wAy3OGcA84MkqjKublDJeCNEA+E0I8a/tm0yF2fH18kb7D/uilDLb9nSFX6+SDlHCcxe/X0pbxyHvtSsc89IVheiJlghuKva03f+m5YhrH1q3Z67t+s0PQIsybuvIuM4bBPwhpSz+Td1Rr9eVOPS9VWMTgZSydyV3EQs0KtYOB+Jtj5OEEKFSygTb6VeyPeISQpRnv/2BfVLKpGL7vvBYCLEM+Lkq45JSxtt+Jwshvkc7Ld2Gk18vIYQBLQl8IaVcXWzfFX69SnC598uV1nEtw7aOjAshRDtgOdBfSpl2/vnL/E0dHlexhI2Ucp0QYokQIrAs2zoyrmIuOSN34Ot1JQ59b9XlrqG/gRZCiCa2b9/3A2tsy9YAj9kePwaU5QyjLMqz30v6Jm0fhufdCZR4h4Ej4hJCeAkhfM4/Bm4rdnynvV5CCAGsAI5JKd+5aJk9X6/LvV+Kx/uo7Q6PrkCWrUurLNs6LC4hRASwGnhESnmi2POX+5tWRVwhtr8fQojOaJ9HaWXZ1pFx2eLxA26h2HvOwa/XlTj2vWXvq9/V4QftP30sUAgkARtszzcE1hVbbwDaXSan0bqUzj9fH9gEnLT9DrBTXCXut4S4PNH+Q/hdtP1nwCHgH9sfO7Sq4kK7K+Gg7edIdXm90Lo5pO01OWD7GeCI16uk9wvwDPCM7bEAFtuWH6LYHWulvdfs9DpdKa7lQEax12fPlf6mVRTXKNtxD6JdxL6xOrxetvbjwNcXbeew1wvtS18CYEL77BpWle8tVWJCURSljqvLXUOKoigKKhEoiqLUeSoRKIqi1HEqESiKotRxKhEoiqLUcSoRKIqi1HEqESiKotRx/wcLmQW6WjRmIAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "monomials = ([[1.],[0.,1.],[0.,0.,1.],[0.,0.,0.,1.]])\n",
    "polyplot(monomials)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dementsprechend ist $p_{k+1}(x) = x^k$.\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def l2prod(coeff_p, coeff_q):\n",
    "    p = (PowSer(coeff_p) * PowSer(coeff_q)).integ()\n",
    "    return p(1) - p(-1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "---\n",
    "\n",
    "### Aufgabe 2: Gram-Schmidt\n",
    "\n",
    "Schreiben Sie den aus der Vorlesung bekannten Gram-Schmidt Algorithmus als Funktion `gram_schmidt(Coeffs)`. Dabei sei `Coeffs` eine Matrix, sodass jede Spalte ein Polynom in Potenzreihenform bestimmt: bestehe `Coeffs` aus den Elementen $\\alpha_{ij}$, so sind durch\n",
    "$$ p_j(x) = \\sum_{i=0}^{m-1} \\alpha_{ij} x^i, \\qquad j = 0,\\ldots,n-1, $$\n",
    "$n$ Polynome mit maximalem Grad $(m-1)$ bestimmt. Für den Algorithmus verwenden Sie `l2prod()` aus Aufgabe 1.\n",
    "\n",
    "- **Überprüfen** Sie den maximalen Orthogonalisierungs-Fehler nach Anwenden des Gram-Schmidt Algorithmus auf die Monombasis $x^0,\\ldots,x^{n-1}$ für $n=5,10$ und $20$.\n",
    "\n",
    "- **Lohnt sich** eine Re-Orthogonalisierung? Falls ja, wie oft?\n",
    "\n",
    "- *Bonus:* Überprüfen Sie auch das Verhalten für den modifizierten Gram-Schmidt Algorithmus aus der Vorlesung. Vergleichen Sie Ihre (re-)orthogonalisierten Basen mit der entsprechenden Legendre-Polynombasis bis zum Grad $(n-1)$.\n",
    "\n",
    "---\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gram_schmidt(Coeffs):\n",
    "    n, m = Coeffs.shape\n",
    "    for j in range(n):\n",
    "        delta = np.zeros(m)\n",
    "        for i in range(j-1):\n",
    "            r = l2prod(Coeffs[:, j], Coeffs[:, i])\n",
    "            delta += r*Coeffs[:, i]\n",
    "        Coeffs[:, j] -= delta\n",
    "        norm = np.sqrt(l2prod(Coeffs[:, j], Coeffs[:, j]))\n",
    "        Coeffs[:, j] /= norm\n",
    "\n",
    "def modified_gram_schmidt(Coeffs):\n",
    "    n, m = Coeffs.shape\n",
    "    for j in range(m):\n",
    "        for i in range(j-1):\n",
    "            r = l2prod(Coeffs[:, j], Coeffs[:, i])\n",
    "            Coeffs[:, j] -= r*Coeffs[:, i]\n",
    "        norm = np.sqrt(l2prod(Coeffs[:, j], Coeffs[:, j]))\n",
    "        Coeffs[:, j] /= norm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  gs: n =  5, maxg = 1.110e-15\n",
      "  gs: n = 10, maxg = 5.400e-13\n",
      "  gs: n = 20, maxg = 1.039e-01\n",
      "\n",
      " dgs: n =  5, maxg = 1.998e-15\n",
      " dgs: n = 10, maxg = 1.023e-12\n",
      " dgs: n = 20, maxg = 5.988e-05\n",
      "\n",
      " tgs: n =  5, maxg = 2.220e-16\n",
      " tgs: n = 10, maxg = 6.690e-12\n",
      " tgs: n = 20, maxg = 8.338e-05\n",
      "\n",
      " mgs: n =  5, maxg = 4.441e-16\n",
      " mgs: n = 10, maxg = 2.881e-12\n",
      " mgs: n = 20, maxg = 4.266e-05\n",
      "\n",
      "dmgs: n =  5, maxg = 8.882e-16\n",
      "dmgs: n = 10, maxg = 1.819e-12\n",
      "dmgs: n = 20, maxg = 3.620e-05\n"
     ]
    }
   ],
   "source": [
    "for n in [5,10,20]:\n",
    "    v = np.identity(n)\n",
    "    gram_schmidt(v)\n",
    "    g = np.zeros((n, n))\n",
    "    for i in range(n):\n",
    "        for j in range(n):\n",
    "            g[i, j] = l2prod(v[:,i], v[:,j]) if i != j else 0\n",
    "    print(f\"  gs: n = {n:2d}, maxg = {abs(max(np.min(g), np.max(g), key=abs)):.3e}\")\n",
    "\n",
    "print()\n",
    "\n",
    "for n in [5,10,20]:\n",
    "    v = np.identity(n)\n",
    "    gram_schmidt(v)\n",
    "    gram_schmidt(v)\n",
    "    g = np.zeros((n, n))\n",
    "    for i in range(n):\n",
    "        for j in range(n):\n",
    "            g[i, j] = l2prod(v[:,i], v[:,j]) if i != j else 0\n",
    "    print(f\" dgs: n = {n:2d}, maxg = {abs(max(np.min(g), np.max(g), key=abs)):.3e}\")\n",
    "\n",
    "print()\n",
    "\n",
    "for n in [5,10,20]:\n",
    "    v = np.identity(n)\n",
    "    gram_schmidt(v)\n",
    "    gram_schmidt(v)\n",
    "    gram_schmidt(v)\n",
    "    g = np.zeros((n, n))\n",
    "    for i in range(n):\n",
    "        for j in range(n):\n",
    "            g[i, j] = l2prod(v[:,i], v[:,j]) if i != j else 0\n",
    "    print(f\" tgs: n = {n:2d}, maxg = {abs(max(np.min(g), np.max(g), key=abs)):.3e}\")\n",
    "\n",
    "print()\n",
    "\n",
    "for n in [5,10,20]:\n",
    "    v = np.identity(n)\n",
    "    modified_gram_schmidt(v)\n",
    "    g = np.zeros((n, n))\n",
    "    for i in range(n):\n",
    "        for j in range(n):\n",
    "            g[i, j] = l2prod(v[:,i], v[:,j]) if i != j else 0\n",
    "    print(f\" mgs: n = {n:2d}, maxg = {abs(max(np.min(g), np.max(g), key=abs)):.3e}\")\n",
    "\n",
    "print()\n",
    "\n",
    "for n in [5,10,20]:\n",
    "    v = np.identity(n)\n",
    "    modified_gram_schmidt(v)\n",
    "    modified_gram_schmidt(v)\n",
    "    g = np.zeros((n, n))\n",
    "    for i in range(n):\n",
    "        for j in range(n):\n",
    "            g[i, j] = l2prod(v[:,i], v[:,j]) if i != j else 0\n",
    "    print(f\"dmgs: n = {n:2d}, maxg = {abs(max(np.min(g), np.max(g), key=abs)):.3e}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Damit lohnt sich eine Re-Orthogonalisierung erst bei hohem n ($n \\approx 20$). Dabei macht eine weitere Re-Orthogonalisierung den höchsten Unterschied, somit reichen insgesamt 2 Re-Orthogonalisierungen aus."
   ]
  }
 ],
 "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}