From 063f4a063f838ffd952585f0319359f49cb11146 Mon Sep 17 00:00:00 2001 From: Kees van Kempen Date: Tue, 15 Mar 2022 10:11:42 +0100 Subject: [PATCH] 10: Ugh --- Week 6/10 Hyperbolic PDEs.ipynb | 79 ++++++++++++++++++++++++++------- 1 file changed, 64 insertions(+), 15 deletions(-) diff --git a/Week 6/10 Hyperbolic PDEs.ipynb b/Week 6/10 Hyperbolic PDEs.ipynb index f029c25..18c5b09 100644 --- a/Week 6/10 Hyperbolic PDEs.ipynb +++ b/Week 6/10 Hyperbolic PDEs.ipynb @@ -66,7 +66,7 @@ }, { "cell_type": "code", - "execution_count": 1, + "execution_count": 15, "metadata": { "deletable": false, "nbgrader": { @@ -83,7 +83,8 @@ }, "outputs": [], "source": [ - "import numpy as np" + "import numpy as np\n", + "from matplotlib import pyplot as plt" ] }, { @@ -162,7 +163,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 35, "metadata": { "deletable": false, "nbgrader": { @@ -185,7 +186,7 @@ " n = len(t)\n", " m = len(x)\n", " \n", - " # TODO: The stepsize is defined by fixed h (k), but the input x (t)\n", + " # TODO: The stepsize is defined by fixed h and, but the input x and t\n", " # could allow variable step size. What should it be?\n", " #h = x[1] - x[0]\n", " l = x[-1]\n", @@ -197,14 +198,14 @@ " \n", " w = np.zeros((n, m))\n", " \n", - " # Set initial values w[i, 0] = f[x[i]] and w[i, 1] = ...\n", - " # TODO: Fix previous comment.\n", - " w[:, 0] = f(x)\n", - " # TODO: What does the length of x[1:]?\n", - " w[:, 1] = (1 - λ**2)*f[x[1:m - 1]] + λ**2/2*f(x[2:m]) - λ**2/2*f(x[2:m]) + k*g(x[1:m - 1])\n", + " # Set initial values for w[i, 0] and w[i, 1]\n", + " w[0] = f(x)\n", + " w[1, 1:m - 1] = (1 - λ**2)*f(x[1:m - 1]) + λ**2/2*f(x[2:m]) - λ**2/2*f(x[0:m - 2]) + k*g(x[1:m - 1])\n", " \n", " for j in range(2, n):\n", - " w[:, j] = " + " w[j, 1:m - 1] = np.dot(A, w[j, 1:m - 1]) - w[j - 1, 1:m - 1]\n", + " \n", + " return w" ] }, { @@ -251,7 +252,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 39, "metadata": { "deletable": false, "nbgrader": { @@ -266,12 +267,60 @@ "task": false } }, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0.9999862004036565\n" + ] + }, + { + "data": { + "text/plain": [ + "[]" + ] + }, + "execution_count": 39, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABEyElEQVR4nO2dd3hcR7n/P7Pqvax6ry5yt2VbbrEdx4kTBydAAglJCNwUAoQWegtccuFHuXAvFxIgBC4ESEIghPimOcV2XCVb7rFsWbJky5LVe5d2d35/zK4tHMtW2d2zq53P8+g52t2z57xnz8z3vPPOOzNCSolGo9Fopj4mow3QaDQajXvQgq/RaDQ+ghZ8jUaj8RG04Gs0Go2PoAVfo9FofAR/ow0Yjbi4OJmVlWW0GRqNRuNVHDhwoEVKGX+5zzxW8LOysigtLTXaDI1Go/EqhBBnR/tMh3Q0Go3GR9CCr9FoND6CFnyNRqPxEbTgazQajY+gBV+j0Wh8BC34Go1G4yM4RfCFEL8XQjQJId4d5fO7hBBHhRDHhBB7hBDznHFejUaj0YwdZ+Xh/wH4JfD0KJ9XA6ullO1CiBuBJ4GlTjq3z9DaM8ju063UtPbiZzIxLTGcZblmQgM9djiFxsuobOpmX3U7rT2DRAT7Mzc9mvlp0ZhMwmjTNE7AKUohpdwhhMi6wud7RrwsBtKccV5fobFrgB+/Xs5Lh+uw2P51/YKIIH/uWZbJp9bmER6khV8zMfZVt/HD105wsKbjPZ9lmkN5ZP00Ns1LQQgt/N6MEQpxH/Da5T4QQjwIPAiQkZHhTps8lrfKGvnC84cZtNi4uyiT9y9IZXpSBBab5Oi5Dp7ZV8MT20/zyrF6nrhrIbNSoow2WeNFWG2SH79+kt/sqCIlKphHby7gupmJJEcH0947xK7KFn63q5rPPXeYV47W89MPzSMiOMBoszUTRDhrxSu7h/+ylHL2FfZZCzwBrJRStl7peIWFhdLXp1Z4pqSGb/7zGLNSIvnFnQvJjgu77H77qtv47LOH6B4Y5rf3FrI8N87Nlmq8kSGLjc88e5Atxxv5yNIMvr2xgJBAv/fsZ7VJfr+rmh++fpJpiRH8+b4lmMODDLBYMxaEEAeklIWX+8xtWTpCiLnAU8AtVxN7DfzjYC3fePEYa6bF8/eHlo8q9gBLsmN56eEVpMaEcN8fSjlU0+5GSzXeiJSSR54/zJbjjXz75gJ+8P45lxV7AD+T4IFrcvj9xxZT1dzD3b/bR9+Qxc0Wa5yBWwRfCJEB/AO4R0p5yh3n9GYO1bTzlb8fZUWemV/fs4jggMtXxJEkRgbz5/uXkhAZxANPH6Cpe8ANlmq8lV+9c5qXj9bztRtncN/K7DF9Z/W0eH59zyLKG7r48t+PotfD9j6clZb5LLAXmC6EqBVC3CeEeEgI8ZB9l0cBM/CEEOKwEMK3YzVXoLNvmIefOURSVDBP3LWIIP+ri72DhIhgnrynkJ7BYb7w18NYbbpCat7Lzopm/nNLOe+bl8InrskZ13fXTk/gKxtm8MrRen6zo8pFFmpchbOydO68yuf3A/c741xTGSklX33hKI1dA/z9k8uJChl/59j0pAi+t2k2X3nhKI9vq+Sz6/JdYKnGW2nvHeLzzx0mPyGCH31wzoSybj5xTQ7H6jr5yZZyVubFMTtVJwp4C3qkrQfx6rEGXj/ewJdvmM789OgJH+f2wjRumZ/CL7ZWUNnU4zwDNV7PD149QWf/MP99x/wJj98QQvCDW+cQExrI1/9xDIvV5mQrNa5CC76H0Dto4T9eKaMgOZL7V42vmX0pQgi+fXMBIQF+PPrSuzrWqgFg7+lW/naglvtX5TAzOXJSx4oKDeC7mwo4VtfJH/accY6BGpejBd9D+MXWSuo7B3js1ln4OWFUY1x4EF/eMIM9p1vZfOS8EyzUeDPDVhvf/Ocx0mND+JyTwnwb5ySzbkYCP33jFE1dOknAG9CC7wGcbe3ld7uquG1RGosyY5123I8syWBuWhT/79WTDAxbnXZcjffxt9Jaqpp7efTmWaOmX44XR0ty2GrjF1srnXJMjWvRgu8B/M/blZiE4Cs3THfqcf1Mgq9tmEFD1wDP7qtx6rE13sPAsJWfv32KhRnRXDczwanHzooL444l6Ty7r4YzLb1OPbbG+WjBN5iq5h5ePFTL3UWZJEQGO/34y/PiKMqJ5fFtp+kf0l6+L/KnvWdp7BrkyzfMcMlcOJ+9Np8APxM/e1MPsfF0tOAbzC+2VhLob+Kh1bkuO8cj66fT0jPIn4rPuOwcGs+kb8jCE9srWZUfx7Jcs0vOkRAZzMdXZLH5yHkqm7pdcg6Nc9CCbyDVLb28dLiOe5dlER/hurlJlmTHsio/jid3VOlYvo/x/P5ztPcN8/nrXDse476V2QQHmPjtjmqXnkczObTgG8jvd1XjbzJx36qxDW2fDJ9cnUtLzxAvHa5z+bk0noHFauN3u6tZlBnj1GSAy2EOD+JDhem8eKiORp2x47FowTeIjr4h/n6glk3zU0iIcH7s/lKW5ZqZmRzJUzurdV6+j/D68QbOtfXzwCTHdYyV+1fmYLHZ+N/dZ9xyPs340YJvEM/uO0f/sJV/W+F67x5UCt39K7OpaOphR0WLW86pMQ4pJb/dUUV2XBjrCxLdcs4Mcyg3zknmL8Vn6RnUs2l6IlrwDWDYauOPe86wPNdMQcrkRjyOh/fNSyEhIoindupJr6Y6h891cKS2k39bkeWUgXxj5b6V2XQPWnTo0EPRgm8AbxxvpKFrwG3evYNAfxN3F2Wys6KFs606Z3oq85eSGsIC/Xj/QveuJrogPZqC5Ej+XFyjQ4ceiBZ8A3hufw0pUcGsneHcQTBj4UOF6fiZBM/tP+f2c2vcQ2f/MC8fPc+m+aluX+dYCMFdRRmcqO/i0LkOt55bc3W04LuZc2197Kps4Xa78LqbpKhg1k5P4G+ltQzrWQ6nJC8erGVg2MZdS41ZF/qW+amEBfrxl2I9utvT0ILvZp4vVZ71hxanG2bDnUvSaekZ5O0TjYbZoHENUkqe2VfDvLQow+apDw/y59YFqbx89DydfcOG2KC5PFrw3YjFauNvpbVckx9PanSIYXasnhZPclQwz+7TYZ2pxpHaTk419nDnEmO8ewd3Lslg0GLj5WN6plZPQgu+G9lZ0UJD1wB3LjHOuwfw9zNx+6I0dlQ060EyU4x/HKwlyN/ExrnJhtoxKyWSaYnh/OOgztbxJLTgu5EXD9URHRrAtTPckxd9JW5dkIqUsPmw9sCmCkMWG5uPnOf6WUlEBI9/eUxnIoTgAwvTOHC2Xc+i6UE4axHz3wshmoQQ747yuRBC/I8QolIIcVQIsdAZ5/UmegctvFnWyE1zkgn0N/45mxMfzrz0aF48pD2wqcL28iY6+ob5wIJUo00B4Nb5qQihWh0az8BZyvMHYMMVPr8RyLf/PQj8yknn9RreKGugf9jK+z2kMgJ8YEEqZfVdlDfoGQ6nAv84WEdceCCr8uOMNgVQGWEr8+L4x6E6bDadk+8JOEXwpZQ7gLYr7HIL8LRUFAPRQghjg4xu5sVD50mNDmFRRozRplzg5rnJ+JkE/9SjIr2ezr5htp5s4pb5qfj7Gd+CdPCBhanUtvdzoKbdaFM0uC+GnwqMTAmptb/3LwghHhRClAohSpubm91kmutp7h5kV0Uzty5IwWRA7v1omMODWD0tnpe0B+b1bClrYMhqY9O8FKNN+RfWFyQR5G/ilaP1RpuiwcM6baWUT0opC6WUhfHx8Uab4zRee7cem1QDUjyN981L5nznAIdrO4w2RTMJXjlaT3psCHPTjMm9H43wIH/WTk/glWP1WLVTYTjuEvw6YGQuYpr9PZ/gtWMN5MaHMS0xwmhT3sO6mYkE+AleO6Y9MG+lvXeI3ZUtbJyT4pIlDCfLxrnJNHcPsv/MlaK+GnfgLsHfDHzUnq1TBHRKKX1CYVp7BimpbuXG2Z7ZZREZHMDKvDhee7dBT3blpbxR1oDFJrnZ4Nz70bh2RgLBATqs4wk4Ky3zWWAvMF0IUSuEuE8I8ZAQ4iH7Lq8CVUAl8FvgU844rzfwZlkjNgkbZicZbcqo3Dg7mdr2fo6f7zLaFM0EePloPZnmUGa5cart8RAW5M+6GYm89m49Fj1/k6E4ZSo9KeWdV/lcAp92xrm8jdePN5AeG+KxlRFgfUEifi8KXnu33rD5VzQTo613iD2nW3nwmhyPDOc42Dg3mVeO1bPvTBvLcz0jbdQX8ahO26lGZ/8wuytbuHF2skdXxpiwQIpyYnntmA7reBtbTzZhtUlu9OAWJKj5mwL9TbxZpifsMxIt+C5k68lGhq3So8M5Dm6cnUxVSy+nGnuMNkUzDt4sayApMpg5Ht4yCwvyZ2VeHG+WNWqnwkC04LuQ1441kBgZxPy0aKNNuSrXz0pECJVCqvEOBoat7DjVwnUFCR7dgnSwviCR2vZ+TuqR3YahBd9F9A5aeOdUMxtmJXnUYKvRSIgIZnFmLK+/22C0KZoxsruyhf5hK+sLPL8FCbBuZgJCoMM6BqIF30W8c6qZQYuNDR6ajnk5NsxO4mRDt17v1kt460Qj4UH+FOXEGm3KmEiICGZ+erQWfAPRgu8i3j7RRFRIAIuzPGfunKuxbqZaY3frySaDLdFcDZtN8taJJlZPiyfI389oc8bM+oJEjtV1Ut/Zb7QpPokWfBdgs0m2l6vK6EkTWV2NTHMYufFhWvC9gMO1HTR3D7K+wPi1FcbD9XZ739JeviF4jxp5EUdqO2jtHbrgMXsT62YmUlLVRs+gxWhTNFfgzbJG/EyCtdO9q4zlxoeTHRfGG1rwDUELvgvYdrIJk1C5x97G2ukJDFlt7KpoMdoUzRV4s6yRJVmxRIUau7LVeBFCsL4gkeKqVroG9ALn7kYLvgt4+2QTCzNiiA4NNNqUcVOYFUNEsD/bdFjHYznT0ktlU4/XhXMcrC9IZNgq2XFq6kyB7i1owXcyjV0DHD/fxbVeGM4BCPAzsXpaPFvLm/Qc+R7KtnL1MPbGkCHAgvRoIoP9eadcC7670YLvZBye8bUzvLMygrK9uXtQT6bmobxzqpnsuDAyzWFGmzIh/P1MrJoWzzunmvWoWzejBd/JvH2yiZSoYKZ74Nz3Y2XNdDVA5u2TumPN0xgYtlJc1eqV/UMjWT0tnqbuQcrqtVPhTrTgO5FBi5XdlS2sneEdQ91HIzYskAXp0To90wMpqW5jYNjG6uneLfhr7A+sd3Qc361owXciJVVt9A1ZvTa2OpJ1MxM5WttJU/eA0aZoRvBOeTOB/iaKss1GmzIpEiKDKUiOZLuO47sVLfhOZLu9Mi7L8f75vh0hA52e6VlsP9VEUY6ZkEDvGV07Gmumx3PgbLtOz3QjWvCdyM6KZpZmx06JyliQHIk5LJCdWvA9hnNtfVQ193p9/N7B6mnxWG2SPZW6jLkLLfhOor6zn4qmHlble793D2AyCVbmx7GzokWnZ3oIjnj3Gi+P3ztYmBlDRJC/Duu4ES34TsLhCa/KnxqVEdS1tPQMcqJBZ1J4AtvLm0mPDSEnzjvTMS8lwM/Eyvw4tpfr9Ex34axFzDcIIcqFEJVCiK9d5vMMIcQ2IcQhIcRRIcRNzjivJ7GzooX4iCBmJHlvOualXGNvreiwjvEMWWzsOd3C6mnxXp0Bdimrp8XT0DWgV1pzE5MWfCGEH/A4cCNQANwphCi4ZLdvAc9LKRcAdwBPTPa8noTNJtlV0cyq/LgpVRkTIoOZkRTBzgrd5Daa0rMqA2z1NO/PABuJI71UT7PgHpzh4S8BKqWUVVLKIeA54JZL9pFApP3/KOC8E87rMRw/30V73zDXTKFwjoNV+XHsr26nf8hqtCk+ze7KFvxNgmW53p2OeSnJUSHkxoex+7RuRboDZwh+KnBuxOta+3sj+S5wtxCiFngV+MzlDiSEeFAIUSqEKG1u9p4n/g67B7wib2p02I5kVX48Q1YbxdWtRpvi0+yqbGV+ejThQf5Gm+J0VubFUVLVxpDFZrQpUx53ddreCfxBSpkG3AT8SQjxnnNLKZ+UUhZKKQvj473HW95Z0UxBciTxEUFGm+J0lmTHEuRvYucp7YEZRWf/MMdqO1g+BR0KgOV5cfQPWzlU0260KVMeZwh+HZA+4nWa/b2R3Ac8DyCl3AsEA1Oi9PYOWjhwtp1V06bE5byH4AA/lmTH6ji+gRRXtWKTyhOeihTlmDEJ2H1atyJdjTMEfz+QL4TIFkIEojplN1+yTw2wDkAIMRMl+FNCQUqqWxm2SlZPwfi9g2vy46lo6tHrkBrE7soWQgP9mJ8ebbQpLiEqJIC5adHs1gOwXM6kBV9KaQEeBrYAJ1DZOMeFEN8TQmyy7/ZF4AEhxBHgWeBjcook3u6saCE4wMQiL1qsfLw4+ib2VGoPzAh2V7awJDuWQP+pO2xmZV4ch8910K2nWXApTilBUspXpZTTpJS5Usrv2997VEq52f5/mZRyhZRynpRyvpTyDWec1xMormpjUWYMQf7eP53CaMxIiiA2LFBnUhhAfWc/p5t7p2w4x8HyPDNWm6Skqs1oU6Y0U9dlcAMdfUOcbOjy+pkLr4bJJFiWY2bv6VY9ItLN7La3qpbnTm3BX5gRQ3CASTsVLkYL/iQoqW5DSqZcbvTlWJZrpr5zgDOtfUab4lPsqWzBHBY4pUZwX47gAD8WZ8XqOL6L0YI/CYqrWgkOMDE3LdpoU1zOhTi+9sDchpSSXZUtLM+Lw2SaOiO4R2NlXhynGnto6tJrMLgKLfiToLiqjcLMqd2Z5iDLHEpyVLDuuHUjlU09NHUPssIHWpBw0anQYR3XMfWVykW09w5xor6LopxYo01xC0KoYf17q1r1dMluwhHemIojuC9HQXIk0aEBF/otNM5HC/4EKalW2QRFOb7hfYHqOGzrHaK8sdtoU3yCXZWtZMSGkh4barQpbsFkEqzIjWN3ZYtODnARWvAnSHFVKyEBfj4Rv3ew3B5a2KNHRLoci9VGSVWrz3j3DlbkxVHfOUBVS6/RpkxJtOBPkOKqVgqzYnwifu8gJTqE7LgwvSSdGzha10n3oIUVeb7TggQuXK/O1nENvqNWTqS9d4iTDd0+Fc5xsCzXTEl1GxarntnQley1t6KW+VgZy4gNJTU6hOIq3Yp0BVrwJ8DF+L1vdNiOZHmumZ5BC8fqOo02ZUpTXNXK9MQIzOFTbwbWKyGEYGlOLMVVbTqO7wK04E8AR/x+Tmq00aa4HYfHqeP4rmPYaqP0TLtPOhSgEiHaeoeoaNLLHjobLfgTwBfj9w7M4WrdXj0Ay3Ucre2kf9jKUh8L5zhwOBU6rON8fE+xJkmbD8fvHSzPjaP0TDsDw3rZQ1dQYl9dbEm2b3r4aTEhOo7vIrTgj5N99sroy4K/Is/MoMXGoZoOo02ZkhRXtTEtMZw4H4vfOxBCUJRj1nF8F6AFf5wUV7XZ8++jjDbFMAqzYjEJ3eR2BcNWGwfOtLF0is/AejWKcmJp6x3iVKOO4zsTLfjjxBG/D/Dz3Z8uKiSAgpTIC6EHjfN4t66T3iGrT7cg4WILWjsVzsV3VWsCtPYM+nz83kFRtplDNR06ju9kiu0LgCz10QwdB+k6H98laMEfB/t8cP6c0Viao+L4R851GG3KlKKkupW8BN+N34+kKEcN8tOT9TkPLfjj4OL8Ob4bv3ewJCsWIS4OQtNMHovVxv7qNp/Nv78URxxf5+M7D6cIvhBigxCiXAhRKYT42ij7fEgIUSaEOC6EeMYZ53U3xVVtPh+/dxAVGsDMpEjd5HYix8930Ttk9fkOWwc6ju98Jq1cQgg/4HHgRqAAuFMIUXDJPvnA14EVUspZwOcne15309ozSHljt08sZzhWlubEcrCmnSGLnlfHGTiEzdfj9w7SY0NJi9FxfGfiDFd1CVAppaySUg4BzwG3XLLPA8DjUsp2ACllkxPO61Z0/P69LM02MzBs42hth9GmTAmKq1rJjQ8jISLYaFM8BpWPrxfdcRbOEPxU4NyI17X290YyDZgmhNgthCgWQmy43IGEEA8KIUqFEKXNzc1OMM157K1qJTTQjzmpOn7vYKl9JKj2wCaPxT5/jq9OpzAaRTlm2vuGOdWkF91xBu4KRvsD+cAa4E7gt0KI6Et3klI+KaUslFIWxsfHu8m0saHy72N1/H4EMWGBzEiK0B23TqCsvovuQYtuQV6CowO7WE/W5xScoV51QPqI12n290ZSC2yWUg5LKauBU6gHgFfQ0jPIqcYenT1xGYpyzJSeaWdYz48/KUrs+fdFPjp/zmikxYSSHhtyYXyCZnI4Q/D3A/lCiGwhRCBwB7D5kn3+ifLuEULEoUI8VU44t1vQ8fvRWZodS/+wlaO1en78yVBc1UpOXBgJkTp+fylF2WZKqnUc3xlMWvCllBbgYWALcAJ4Xkp5XAjxPSHEJvtuW4BWIUQZsA34spTSa9poxVWthOn4/WVxzOiop1mYOFabZF91m47fj4KO4zsPf2ccREr5KvDqJe89OuJ/CTxi//M6dPx+dMzhQUxLDKekqo1PrTHaGu/kxIX4vQ7nXA5Hmure063MSIo02BrvRivYVbgYv9fe12gszTZTekavcztRLuTf6wFXl8URxy/RcfxJowX/KlzoTNPe16gszYmld8jKu+e7jDbFKymuaiU7LoykKB2/Hw0dx3cOWvCvgiN+P1vH70fF4ZnqfPzxcyF+r7NzrshSHcd3Clrwr4KO31+d+IggcuPDKNGCP25O1HfRNaDz76/GhUF+Oh9/UmgVuwItPYNUNPXo+XPGgCMfX8fxx4dj0JqeP+fKOObV0YP8JocW/CtwMX6vBf9qLM0x0z1ooaxex/HHQ3FVK5nmUJKjQow2xeNZmq3nx58sWvCvwN6qFhW/T9GpYFfDMUJUZ1KMHZs9fl+ks3PGhJ4ff/Jowb8CxVVtLM6OxV/H769KQmQwOXFhegDWODjZ0E1n/7AO54wRR0tbl7GJo5VsFJq7B6ls0vn342FpTiwl1W1YdZN7TFyc/16XsbGg17mdPFrwR8HhRWjBHztLs810D1g4oeP4Y6KkupUMu4hpxsbSnFhKqtpQg/c140UL/igUV7USHuSv4/fjwBGa0B7Y1bHZJCU6/37cFOWYae0dolLH8SeEFvxRKK5qY3FWjI7fj4PkqBAyzaE6dW4MlDd209E3rFuQ46RID/KbFFrNLkNT94CO30+Qomwz+3Tq3FUp0evXToj02BBSooL1/PgTRAv+ZdD59xNnaU4snf3DnGzQQ+CvRHFVG2kxIaTFhBptilchhKAoR82ro+P440cL/mVwxO9n6fj9uHFknOgm9+jYbJJ9Z9q0QzFBlubE0tIzxOlmHccfL1rwL0NxVauO30+Q1OgQMmJDda70Faho6qGtd0h32E4Qx4Nyrw7rjButaJfQ1D3A6eZe7X1NgiJ7Pr6O418eR+tHl7GJkREbSnJUsJ6sbwJowb8ER/xeT5g2cYpyzHT06Tj+aBRXtZIaHUJ6rI7fTwQhBEuzYynW+fjjRgv+JeytaiUiyJ+CZB2/nyhFOo4/Ko78e+3dT46iHDMtPYOcbu412hSvwimCL4TYIIQoF0JUCiG+doX9PiiEkEKIQmec1xUUV7Xq+XMmSUq0ysffqwX/PTji93oFtcmhnYqJMWlVE0L4AY8DNwIFwJ1CiILL7BcBfA4omew5XUVT1wBVzb26MjoBnY9/eXT83jlkmkNJigzWg/zGiTPc2CVApZSySko5BDwH3HKZ/R4DfgQMOOGcLqG4WuffO4uiXJWPf6JBz6szkr2nW0mL0fH7ySKEYGlOLMVVOh9/PDhD8FOBcyNe19rfu4AQYiGQLqV85UoHEkI8KIQoFUKUNjc3O8G08VGs4/dO42KTW3tgDlT8vlU7FE6iKMdMc/cgVS06jj9WXB6oFkKYgJ8BX7zavlLKJ6WUhVLKwvj4eFeb9h6Kq1pZouP3TiE5KoQscyh79RqkFzjV1E173zDLtOA7hQvz42unYsw4Q9nqgPQRr9Ps7zmIAGYD24UQZ4AiYLOnddw2Xojf68roLIpyzOyrbtXz49txPPz0/DnOIcscSkJEkO64HQfOEPz9QL4QIlsIEQjcAWx2fCil7JRSxkkps6SUWUAxsElKWeqEczsN3ZnmfIpyzHTp+fEvUFzVSnqsnj/HWTjm1dFx/LEzacGXUlqAh4EtwAngeSnlcSHE94QQmyZ7fHdRXNWm4vd6/hynoVPnLnIh/16vX+tUinLMNHUPcqa1z2hTvAJ/ZxxESvkq8Ool7z06yr5rnHFOZ+OI3/uZhNGmTBmSooLJjgujuKqV+1flGG2OoZxsUPPf6xHczmXkojvZcWEGW+P56N5JoL6zn+qWXl0ZXUCRXucW0OvXuoqcuDDidRx/zGjB52JnmhZ851OUo9a5LTvv23H84iq9fq0ruDA/vp5XZ0xowUcJfnRoADOTdPze2eg4/sj5c3R2jitYmh1LQ9cAZ3Uc/6powQf2nG6lKNuMScfvnU5iZDA59ji+r3KioYvOfh2/dxXaqRg7Pi/459r6qOvoZ3meroyuYmmOmlfHYrUZbYohOEYbL9UZOi4hNz6MuPAgPa/OGPB5wd9zugVAj350IctyzXQPWijz0Xz84qpWMs2hpOj4vUvQ8+qMHS34p1uJCw8iLyHcaFOmLEXZF1PnfA2rTVJS1aodChdTlGOmvnOAmjYdx78SPi34Ukr2nm5lea4ZIXT83lUkRAaTEx/mk/PqnKjvomvAokdwu5hlOb7rVIwHnxb80829NHUP6s40N7Asx8z+M+0+F8e/mH+vM3RcSW58OHHhgXp21qvg04LvWJFpuRZ8l1OUY6Zn0MJxH8vH33u6lSxzKMlROn7vSlQc38ye0y06jn8FfFvwT7eQEhVMhl6MwuU4PFxfWvZw2GqjpLqNFXlxRpviE6zMi6Oxa5DTzT1Gm+Kx+Kzg22yS4qo2luXG6fi9G0iICCY/IZzdlS1Gm+I2jtZ20DNoYaUWfLfg+J13VfhOGRsvPiv45Y3dtPUO6fi9G1mRF8f+M20MDFuNNsUt7KpoRQg9ZYe7SI8NJSM2lN0+mBwwVnxW8Pfo+XPczqr8OAaGbRw82260KW5hd2ULc1KjiA4NNNoUn2FFXhzFp1t9LjlgrPis4Ds60/RkVu5jaY4Zf5Nglw+EdXoHLRysadfxezezMi+O7kELR+s6jTbFI/FJwbdYbZRUt2rv3s2EB/mzICPaJwR/X3UbFpvU8Xs3syzXjBCwW8fxL4tPCv6R2k66ByyszHP/Qum+zoq8OI7VddLRN2S0KS5lV2ULQf4mFmXGGG2KTxEbFsislEifcComgk8K/s6KZoTQ+fdGsCo/Dikv9qFMVXZXtrA4K5bgAD+jTfE5VuTFcbCmnb4hi9GmeBw+Kfi7KlRnWkyY7kxzN3PTogkP8p/SHlhT9wAnG7p1/N4gVubFMWyV7NOzZ74Hpwi+EGKDEKJcCFEphPjaZT5/RAhRJoQ4KoR4WwiR6YzzToTugWEOnetgVb6ujEYQ4GeiKMc8pXOl91Sq1ouO3xvD4qxYAv1NPjXmY6xMWvCFEH7A48CNQAFwpxCi4JLdDgGFUsq5wN+BH0/2vBOluEqtr6rj98axMs9MTVsfNVN0haJdlS1EhwZQkKJXUDOC4AA/FmXEsLtyaocNJ4IzPPwlQKWUskpKOQQ8B9wycgcp5TYppaN2FwNpTjjvhNhZ0UxIgB8LM6ONMsHnWZmvHrZTMawjpWR3ZQvLc8346RXUDGNlfhxl9V209gwabYpH4QzBTwXOjXhda39vNO4DXrvcB0KIB4UQpUKI0ubmZieY9l52VbSwNCeWIH/dmWYUufFhJEUGT8kmd1VLL/WdAzp+bzCO33+qJweMF7d22goh7gYKgZ9c7nMp5ZNSykIpZWF8vPNDLnUd/VS19LIqX4dzjEQIwcr8OHafbsFqm1ozGzoeYjp+byxzUqOICPafkk7FZHCG4NcB6SNep9nf+xeEENcB3wQ2SSkNaWftqlCtBt1hazwr8+Lo6BumbIpNl7yrooW0mBA9A6vB+JkEy3PN7KzQ0yWPxBmCvx/IF0JkCyECgTuAzSN3EEIsAH6DEvsmJ5xzQuyoaCExMoh8vZyh4Tia3DsrXRO6M4Ihi409p1tZlR+vZ2D1AFZPS6Cuo5/KJj1dsoNJC76U0gI8DGwBTgDPSymPCyG+J4TYZN/tJ0A48DchxGEhxOZRDucyrDbJnsoWVuTp6ZA9gfiIIGYkRbDj1NQR/NKzbfQMWlg7XYcMPYE19vuwvXzqlLHJ4u+Mg0gpXwVeveS9R0f8f50zzjMZjtR20N43zJrpCUaborGzZnoCT+2sontgmIjgAKPNmTTvlDcT4CdYruP3HkFKdAjTEsPZfqqJB67JMdocj8BnRtpuP9mEScA1On7vMaydHo/FJqdMx9r28mYWZ8USHuQUP0rjBNZMT2B/dTu9g3qaBfAhwd9W3szCjBg9N7kHsTAzhohgf7ad9P4m9/mOfsobuy+EETSewZpp8QxZbTo9045PCH5T9wDH6jpZO0OHczyJAD8T1+THs/1Uk9dnUrxj74vQIUPPYlFWDKGBfmwvNyxXxKPwCcF/p9xRGbX35Wmsnh5PY9cgJ+q7jTZlUmw72URKVLDOAPMwgvz9WJ4bx/byZq93KpyBTwj+9vJmEiKCKEjWc5t4GmumqYfwNi/2wIYsNnZXtrBmRoLOAPNA1kyPp66jn9PNOj1zygv+sNXGjopm1k7XldETSYgMZnZqpFc3uUvPttE7ZL3w8NJ4Fjo98yJTXvAPnm2ne8Ci4/cezNrpCRw4205n37DRpkwInY7p2aTFhJKXEK4FHx8Q/G32yrgiT69u5amsmZ6ATcKOCu+skFtPNul0TA9n7fR4Sqpb6R7wTqfCWUx5wd9erirjVBjYM1WZnx5NdGgA2056X1jnTEsvFU09XDcz0WhTNFdgfUESw1bp817+lBb8c219nGzoZq1OlfNo/EyCtdMT2FrehMVqM9qccfFmWSMA6wu04HsyizJjiA0LvHC/fJUpLfhbjjcAcP0sXRk9nRtmJdLRN8y+M961DukbZQ3MTI4kXc+O6dH4mQTrZiSwrbyJIYt3ORXOZEoL/hvHG5mRFEGmOcxoUzRX4Zpp8QT5m3jjuPd4YC09gxw428712rv3CtYXJNI9YKGk2ndH3U5ZwW/pGWT/2Taun5VktCmaMRAa6M+q/HjeON7gNQNktp5owiZ1OMdbWJUfT3CAyafDOlNW8N8qa0RKFSrQeAc3zErkfOcA79Z5x6Iob5Q1kBodwiy9WLlXEBLox6r8eN4sa/Qap8LZTFnB33K8gbSYED261otYNzMRk7jY9+LJ9A1Z2FnRwvqCRD2gz4tYX5BIvRc5Fc5mSgp+98AwuytbuWFWkq6MXkRsWCBLsmO9QvB3nGph0GLT8XsvY92MBEwC3izz/DLmCqak4G8vb2bIauMGHb/3Om6YlURFUw9VHj7vyZtljUSFBLA4O9ZoUzTjwBweRGFmLG/4aBx/Sg4N3HK8AXNYIIsyY/71A8sQmPxhoAMCQqGnEYIioLsBgsKhpxn8A2GwG2xWsFnU/wGh0HkOIlOh8RgkzoEzOyB7DRx/EebcBvufgqUPwc6fwuqvwFvfhfWPwWtfgY0/hZc/D+/7OWz+HGz6H/g/+/aVL8HG/4Q3vg3XPwbbfwirvwrFv1LHO/oczP0wVL4Neeug/jAkL4DOGojJhsEuCDWDtIF/CPgHgRAQGAY2m7o+aYOAYHffhgmxviCRf/+/Mt4oa+Sh1Z458+SQxcZbJxq5dkYCAX4jfCZHXHi4T/1vHQLLICBhoEvdg846iEiCxuMQPwPO7ITMFVD2EhTcAqW/h8J/g10/g5VfgLf+HdY9qsrRhh/Cy19Q5eX/Pgcbf6bK1cafwStfVJ+/+W21//YfwapHVDlafD8c/SvMej+c3grZ10DDUUicDZ21EJUGQ72qDpgClL0hMeq9yFToa4XYbOhpgthc6K6H2BzobVbfHeyGkGj334gJcsPsJB57uYyq5h5y4j2zjLmKKSf4g+11fLz8kxRl3oXfX56AZQ/DW99R220/gMX3wf7fQcEmOPkKZCyH2n1gzof2M6qgD/eBbRgCI+yFOxsajkH6UlVhpt8Ex/8BvS1w5FlAwsmXISgSqrZBeCKcK4HS3ymB3vek+n7Jb9QDo+TXqsIV/xrqStX7Z3ZCyZNw6nUIjYNjzyvxPvQnJRoH/whdtXD4GZh9G5S/BrlrobYUkuZAR40SkuF+8PNXtvS3gzkPWiogayWcK4bZH4TKt2DBPVD+Kiz6GJx8FRZ+FCrfhDm3w5ndkHcdNJ9Qx+5vV7+LG0iLCWV2aiSvvdvAQ6tz3XLOC0ipHpZd59U9rD8CcdOgegekzIdTWyB1IbW7X2DBYCBfsDTBgYVQU6KEr/mketAO9cDwAASGqvsSmwN1ByBjmbq/M98HR55Tv/n+p2DxA7D/t9DToATfMqAE2i8QTr2mHuhndsLex1VZ3fsEnD+kyk39EVW+6g/Dgf+F2v1w6M9wdpeyqfJN9f3j/1DOzrsvQH8bVLylysT5g5AwU9kZFq/sBrsD1ATmXGg6ocp+9Q5Vb959QT1E9j0Jq78GO36sHjZvfAtu/RW8+R3Y9AvY9h9w03+qh9e676prW/4ZOLFZleHzByG1EHqb1O/tpvDrTXOU4L98tJ7Prst3yzn/hYFO8A9W2hIYDh1nISBM1W/hp+q1i3CK4AshNgA/B/yAp6SUP7zk8yDgaWAR0Ap8WEp5xhnnvpR2WzjzxSky5C4lzqYAJbZHnlU/7KktaltTDG1VqoK2nlaVvbNWeSrD/corDgxTXgxSbZvK1MOg7oD6vGavOumZnfbtLrU9u9v+eve/vn9hu/uS/Rzfc3y+45Lj2rdn96pWR/1hGOxUAtNVBwEhqsIOdcNgD5j81Hv9HeqvtVJ9r+4ACBNUvQNWixKTgS4lBl3n4fCfoekk7PsNXPMV2PETuOkn8PrX4LbfKy/y9j/C1sfg5v9SArTm66qVs+AudfzMFarlFJWubBamcVfkjXNS+NHrJznX1je5AU0DXapCNRyFhAIo+ydM2wC7fw6FH4ct34RVX4SXPgXX/Tv87WNw44+UF73ma/D2Y7DsU7Dnl7DwHiWkBbeSefwl/i1wLunVJ2CoEuqPQmwWdNUrobRalGgHR6n7YxmC9rOqtdXfDnUHlffvKAfV76htlX1bfen9d7y+pHxdWn4u3V66X00xSKt6WAz3qvLT2wxtdvEZ6lXOhZSq7Pe1AvZ6ERShhPncPnUNp7cq4TqxWb0+9Gd1rOJfQdNx2Pmfap+3v6d+d5tV1UHHQ625XD0I1j+mWsM3/5favu/n6hirvwIn/g/m3anKe9ZK5WDF5qhrCI5Soin8lIMzDpKjQlicFcPLR89PTvAtg8opaz6lHKvKN1X5P/QnmLFR1Y85t6tW/7w7YMd/qtcH/wh569XDP3meKp+Rqare+AXBJ3dN3KarICabniSE8ANOAeuBWmA/cKeUsmzEPp8C5kopHxJC3AG8X0r54Ssdt7CwUJaWlk7MqF8uVgJoGVBektUeyrFZVAGRVkAAcsR2qmG/LmFSDyfHdZsCVOvFLwisgypcNdynWjND3RAcrUJeYfGqAkemKtGKyYb2aiWcTWXKM6srhZw1ULVdhQuOvwiF98HBp2Hl51XFvubLqiW14B7ljeasUWGouHxlV6hZ2eAXoCoxakqMVT/expdvmM6n1+apy7Fa1H20DauK7xeoBCsyRYW7kuepFtWs98PW/4Aln1Ae5/LPKPFZ+FFl1/SbVMsmc6V6wCbPV4ISNx1ayiE6Q5Wd8CQlTqFmJXzBUTDQiQwMRwz1MCSCCJSDF3/HC+XL/nvDiP+nYhkbY/nyDwFLv3rwDvVc+B0JjYO+FohIVg+b6EzliMXPUPc1ZYF6MGWtUg++/OtVOcu/Xj24clYroU2ep84TkaTulclfHWO4X20HOlVZG+yCyDRVhgKC+eOeM3xn83He+MI1TEuMuHhZ1mH7/U9UD87URepBNW2DarEs+QS88ogq1//8FFz7TXjz0YuttIJbVHgu+xp7y3ChasnEz1Qt5pgsFUmISIHu8xfLV2CEegib/OEb51V9mOidEeKAlLLwsp85QfCXAd+VUt5gf/11ACnl/xuxzxb7PnuFEP5AAxAvr3DySQn+8x9VP7rGOBxC4BAAx4PF8UAJjVOVMCRGeUoBIcqrtA5BRAq15+toJ5I5sTZVCYZ7VWslyB5mC0tQTeCwBOV5BkWpVo9/sHrQOwR4Soqt5sJ9ddxnk//FkJzDyQuOUv0L4UmqjJjz1UMleR7W9hpK2sPIjbCRGIoqey0VEJOpIgLmfGituCjMjnLrKMeOcj3yAe8kTkcsptUvniWff3ZC37+S4DsjSycVODfida39vcvuI6W0AJ3Ae+YrFkI8KIQoFUKUNjdPYla7+JkT/65m8gg/VQn8Au19IeGqkoTEqkoTnqQ8r7AEe+gjUgk1KE/POoTFPJ2+YSu9IUmq2RwUobzBwS61T0+j2vY22d/vVJXSMqAqpc1ysQNSTMlkNN9FmFBiH6Dus1+Q3bnwt5c5ixL7oV7l+fe1Qvx0lXiRvhS6G/BLmUdsMJwdCEFGpsBQn+qnaamAtMVK7BMKlNhHpqlyGxKrynFAmL0VE3ixdeNEcrv3Y7K5Zhpnj+q0lVI+CTwJysOf8IESZqitozl5odnteCrbm5/vaXJ7szd4mWt4TzPb7vk4vOCAMOU5B0Ve9Lb72y96zZFpyouOzVH9HQmzVHw2bbHqHMxZqzqpHeGcxffbwzhfUh3Ta74B7/5dNYNr9sD0jSpMkjRXnT8iSdkUEKJsESYwKXEO7R7gzh+8zcNpeTxy/XR1PY7YuLSqzKrAMBUmismG8lcgvQj2/lLFfd/4Nqz4LLz+dRVe2vGTEWGdjWp/R7jAET5whBPe0+xW4QcZHI0Y6GDAFEqwre+it3dpWAfeG+qYEuGd0cI49ut+T/kaa7gwS/3ejrBH6iLVH+QIi0y/SfW9zXyfCutM26DCJBlF0NemkioCwpRTEJut+gvMuUrwzXnKy49KUy3JwIt9QgdLavjGi8fYvHYFc9Oi1ZvD/aq+tFWpcGHVdvWQOPQ0zHifPQHk07D5s6r/6qVPwdpvqGyqxfepTviZm1TfhiPc6Qh/OsKhjvCoo3w5HCH/EGyWIUxymKCFd7jkDjrD9akD0ke8TrO/d9l97CGdKFTnrWtwePj569XW0eudvUptM5apbeoiVXgTZ6tKG5OpQgNh8coLDQhVXqQwqYIDKmsDLr6OzVHbmGz7NkttozPHtnXs7/i+43iO45vtnUpxdtFLKLBvZ118bfJXBTwoUhXsiBQlUuZcVdmS56rryyhS15KzRn132g1qO3OT2s65XW0X3QsIJZimAJXmF2qGm36sxH7TL5So3vI4rPsO3PzfcMezcP334dP7VFrplypg6Sfg/rdg7u2qUy7/Olhwt7InfYmKl0ckKW/Mz/+C2AMkRASzLNfMPw+fvzgM3s9fdYoGRymPLSpNdY4lFqiYavYq+MhfYdat8IVjsOQB+NIp9eC5/2246adwyxPw/l/Dis+rjJKZ71P2J86G6/9D3f9VX1LOwuL71G877w4QfpyJX8uw9KM/Y7UStYxlqvWSPE9VWnO++u0j09R9DbLbaQpQxwe1L6iHHlx831FmHffZcd8d5eFCOZto+coa5XjZqkxEpqnyHhKryn9AqPqOyf+ibQ7bUxepbfpStc1crrZZ9vqVu05tL5Svm9V29gfVdsE9arv0IXXuVV9SZffabypxXPcdVU7WfRdu+IF6fdfzcMP34YGtqmP9wXdUOuodf1H3be3XYfnDqjwUbILEWaqMhcaqem3y+xexB9g4N5lAfxMvHKi9+GZAiCpn8dNUWZt5M4THq879+Glw57OqA/mzB1W5/uoZlT77yT1w449VPbjll8rum/9LOTvX/4dyMlY+opyLBXep33T6BqUnmcvVNmkOdQHptBFJwcpbcAXO8PD3A/lCiGyUsN8BfOSSfTYD9wJ7gduArVeK30+auHyYf7eqsJYB5Xn2NKpUsI5zkHet8mDTFqunf6JdPCNTlOfhF6higQOdEJ6gYoPmPNU0TJylMi7SFqvsl8wVyhvIWqme2pkrlbeStRIOn736duT+7dVq6zhea6USsdYKyFqhtqmL1PkTZyl74vKVdxWbqyprcKSqpNYhFeroqlPfCYlVBSwsXlW8qDSY9xElSPPugLRC5annXQe516qHQNIc9ZuFmdVnoAQclGiCyvUGmHHTxd8eVAWbJB9cmMYjzx9hX3UbS3MmuGJZqH1gVJo9pLngLrVd/+9q++E/q+0n7ZksX61WtuetUy2d3GvVvc9ayU/eHiI7OIUvrb8Nqharh23CTOVV1h9VD8WOmovn7mtRxzD5KWHvb1cPu7YqJZaNx9V9bnxXebPNJ9T9bikfcd9Xqf3/pZydGVF+Vqi4tGM7WvnKXKG2qYVqLEDiLHvabq7yiKPTlYccEKzKjyMcIkyqxTzUo0Ie3Q32Mlqt7nlrBcy9A1oq1QO+t0kJoHUI1n5LHfPaR9WDZfln1H2Y+2HljGUsU2NYIpJVFhSoBzCoTlmANPvDxdFqdyJRIQGsL0hk85HzfHNjAYH+E/B/A+0z8SbaHTFHPVj2abW96cdqm2l3MufanarFD4xIzQyDrjo6hwTf/u3r3Fhg5sP+rlmwadKCL6W0CCEeBrag0jJ/L6U8LoT4HlAqpdwM/A74kxCiEmhDPRRch8kPbn1c/X/X39TW4e3P/bAqxIX/pry4ZZ9W3kxfq4oV93eozy0DqgkoTBezNGbcrLyFlIWq8EckK4H0C1CF3jKg0v3629VDprteVYKOGij6lKpwyx5Wgr3sYZXutvxhVUmWfVo9fJZ+UjWXlzyovMd5d6imb/71yitMK1T/m/OUJx6eqCpXYJhqYpv8lZeCVN+XNnVdDhbfr7YzNqqto0IVflxtHV5Zst0DDTNuacgbZyfz6EvHeb60duKCP14cD6rIFLVNXQjAufjVvFqzjS+u/zgiNf/C+0y/8fLHcaT2WodVuRAmVS4Cw2D+Xars5K5VD4yEAiXwEUmqdeIfrJwV6zCs+JwKt638gnJGVnxeCfbyz6iY9LJPq7Db0k+o8EDhfWq/hR9V5Wneh5Udsz+gykPeOuW5pyxQIbDYbFXmQ2Ls/S4BFzvPw+LVMRyDr6IzVTmOSleCFZ4IC+9VTobjQeooP3f/XW03/UJt13xVbefbfUFHq8DxOxvEbQvTeOVoPdvKm9w7Mt8xUC3W3rIPT2Dz3jNst87hK2tXuey0k87ScRWTytLRTBm+9sJRXjp8nv3fus7QNWN//lYF//XWKXZ9dS1pMXqxk6mCxWqj6P9tZWFGNE9+9LKJLW7jlsd3Mzhs5fXPXzOp47g6S0ejcRm3F6bRP2zl1aP1htlgs0leOFjL8lyzFvsphr+fiVvnp7CtvImWnkHD7Dh+vpMj5zr4UGH61XeeBFrwNR7NwowYcuLDeL703NV3dhG7T7dQ09bn8sqoMYY7lqQzbJWGlrFnSmoI8jfxwYVpLj2PFnyNRyOE4I7F6ZSebedEvTFzmP9p71liwwK5cY6efXUqkpcQQVFOLM+U1GC1uT/E3Tto4aXD59k4N5moUNd01jrQgq/xeD5UmE6Qv4mn9551+7nPd/Tz1olGPrw4nSB/5w6w0XgOdxdlUtvez45TkxjwOUH+78h5egYt3LU0w+Xn0oKv8XiiQwO5dX4q/zxUR2efa0YgjsZz+2qQwEeWuL4yaozj+oIk4sKD+HOxe50KKSV/LjnLtMRwFma4fkZaLfgar+CeZZn0D1v52wH3xVkHLVae3X+OtdMTJjdrp8bjCfQ3ccfidLaWN3G2tddt5y2uauPdui4+tjzbLavzacHXeAWzU6MozIzh6b1nsVidO1nVaPzzUB3N3YN8bHmWW86nMZZ7lmUSYDLx1M5qt53ztzurMIcF8oGFl04/5hq04Gu8hvtXZVPT1scrx1yfomm1SX6zo4pZKZGsyo9z+fk0xpMYGcz7F6TyfOk5t6RoVjZ1s/VkEx9dlkVwgHv6h7Tga7yG6wuSyE8I54ltp7G5OJvizbIGqpp7eWh1rlua2hrP4MHVOQxZbTy954zLz/XUzmqC/E3cXeS+/iEt+BqvwWQSfGptLuWN3bx1wnWLUEsp+dU7VWSaQ7lxtk7F9CVy48O5viCRP+49S8+gxWXnOdfWxwsHa/lQYTrm8KCrf8FJaMHXeBXvm5tCRmwov9xWiaumBXn7RBNHznXwiWty8ffTVcTX+PTaPDr7h3lqZ5XLzvHztyswCXFxRTc3oUuzxqvw9zPx8No8jtZ2uiSWb7Ha+NHrJ8mJC+P2QteOetR4JnPTotk4J5nf7qiiudv5sfzKpm7+cbCWjy7LJCkq2OnHvxJa8DVexwcXpTEjKYIfvX6SgWGrU4/9wsFaKpp6+PIN0wnQ3r3P8sXrpzFgsfHLrRVOP/bP3jxFSIAfn1zjXu8etOBrvBA/k+BbGws419bPH53YudY7aOFnb55iQUY0G3Ts3qfJiQ/nw4vT+UtJDZVNPU477q6KFl491sCD1+QSGxbotOOOFS34Gq9kZX4c185I4JdbKznf0e+UY/5kSzlN3YN8a+NMnZmj4QvXTSMsyJ+v/+OoU7LCBoatfPuld8mOC+MTq3OcYOH40YKv8Vq+874CLDbJV184OukO3NIzbfxx7xnuXZbFosxYJ1mo8WbiI4L41saZ7D/Tzl9KJj/lwhPbT1Pd0stjt8x2W979pWjB13gtmeYwvrFxJjsrWvhLSc3VvzAKfUMWvvLCUVKjQ/jyDdOdaKHG27ltURqr8uP44WsnqWntm/BxSs+08cS2Sm6dn8JKAwfyacHXeDV3L81gVX4c33/lBO/WdY77+zab5Et/O8KZll5+9MG5hBm4qpbG8xBC8IP3z8Hfz8QDT5fSO4Hc/ObuQT79zEHSYkL43q2zXWDl2NGCr/FqhBD89PZ5xIYF8vE/7Ke2fXxe2P9sreDVYw18/caZrMjTUyho3kt6bCiPf2QhFU3dPPL84XHF8wctVh5+5iCd/cP86u5FRAa7dr77qzEpwRdCxAoh3hRCVNi375nfUwgxXwixVwhxXAhxVAjx4cmcU6O5lITIYP7344sZGLby8f/dT1P3wJi+9/TeM/z3WxXctiiN+1dlu9hKjTezMj+Ob9w0ky3HG/nmP4+NSfQHhq18+i+HKKlu44cfmMvM5Eg3WHplJuvhfw14W0qZD7xtf30pfcBHpZSzgA3Afwshoid5Xo3mX5iWGMGT9xRS297PB57Yw7Ha0cM7QxYb33+ljEdfOs51MxP5/vtn66wczVW5b2U2D6/N49l953jwT6V09A2Num99Zz93P1XCWycaeeyWWdy6wD2zYV4NMZnsBiFEObBGSlkvhEgGtkspr9jrJYQ4AtwmpbziiIbCwkJZWlo6Yds0vsmRcx188s8HaOoe5K6lGdy7PIvsuDCEEAwMW3nrRCM/f6uCiqYeProsk2/fXKAHWGnGjJSSp/ee5bGXy4gODeQz1+Zx6/zUC0sTNncP8nzpOX79zmmsNsmPb5vLzXNT3GqjEOKAlLLwsp9NUvA7pJTR9v8F0O54Pcr+S4A/ArOklO+Z1FwI8SDwIEBGRsais2fdv6Sdxvvp6Bvix1vK+ev+c1htkrjwQEID/WnoHGDIaiM7LoxvbZzJupmJRpuq8VKOn+/kOy8dp/RsOyYByVEh2KSkvlOFE6+dkcC3by4gOy7M7bZNSvCFEG8Blxt2+E3gjyMFXgjRLqW87DpdjhYAcK+UsvhqRmsPXzNZGrsGeP3dBsrOdzFgsZIUFcyK3DhW5MXhZ9IhHM3kkFJyrK6TrSebONvahxCQlxDO+pmJ5CdGGGbXlQT/qjloUsrrrnDgRiFE8oiQTtMo+0UCrwDfHIvYazTOIDEymHv1alUaFyGEYG5aNHPToo02ZcxMNni5GbjX/v+9wEuX7iCECAReBJ6WUv59kufTaDQazQSZrOD/EFgvhKgArrO/RghRKIR4yr7Ph4BrgI8JIQ7b/+ZP8rwajUajGSeT6rR1JTqGr9FoNOPnSjF8nY+m0Wg0PoIWfI1Go/ERtOBrNBqNj6AFX6PRaHwELfgajUbjI3hslo4QohmY6NwKcUCLE83xBvQ1+wb6mn2DyVxzppQy/nIfeKzgTwYhROloaUlTFX3NvoG+Zt/AVdesQzoajUbjI2jB12g0Gh9hqgr+k0YbYAD6mn0Dfc2+gUuueUrG8DUajUbzXqaqh6/RaDSaS9CCr9FoND6C1wq+EGKDEKJcCFEphHjP4ulCiCAhxF/tn5cIIbIMMNOpjOGaHxFClAkhjgoh3hZCZBphp7O52nWP2O+DQggphPD6FL6xXLMQ4kP2+31cCPGMu210NmMo3xlCiG1CiEP2Mn6TEXY6CyHE74UQTUKId0f5XAgh/sf+exwVQiyc9EmllF73B/gBp4EcIBA4AhRcss+ngF/b/78D+KvRdrvhmtcCofb/P+nt1zzW67bvFwHsAIqBQqPtdsO9zgcOATH21wlG2+2Ga34S+KT9/wLgjNF2T/KarwEWAu+O8vlNwGuAAIqAksme01s9/CVApZSySko5BDwH3HLJPregFkwH+Duwzr7Qurdy1WuWUm6TUvbZXxYDaW620RWM5V4DPAb8CBhwp3EuYizX/ADwuJSyHUBKednlRb2IsVyzBCLt/0cB591on9ORUu4A2q6wyy2olQKlVEvDRtuXkp0w3ir4qcC5Ea9r7e9ddh8ppQXoBMxusc41jOWaR3Ifyjvwdq563fambrqU8hV3GuZCxnKvpwHThBC7hRDFQogNbrPONYzlmr8L3C2EqAVeBT7jHtMMY7x1/qpcdRFzjfchhLgbKARWG22LqxFCmICfAR8z2BR3448K66xBteR2CCHmSCk7jDTKxdwJ/EFK+VMhxDLgT0KI2VJKm9GGeQve6uHXAekjXqfZ37vsPkIIf1QTsNUt1rmGsVwzQojrgG8Cm6SUg26yzZVc7bojgNnAdiHEGVSsc7OXd9yO5V7XApullMNSymrgFOoB4K2M5ZrvA54HkFLuBYJRk4xNVcZU58eDtwr+fiBfCJEthAhEdcpuvmSfzcC99v9vA7ZKe0+Il3LVaxZCLAB+gxJ7b4/pOrjidUspO6WUcVLKLCllFqrvYpOU0psXRB5L+f4nyrtHCBGHCvFUudFGZzOWa64B1gEIIWaiBL/ZrVa6l83AR+3ZOkVAp5SyfjIH9MqQjpTSIoR4GNiC6t3/vZTyuBDie0CplHIz8DtUk68S1TFyh3EWT54xXvNPgHDgb/b+6Rop5SbDjHYCY7zuKcUYr3kLcL0QogywAl+WUnptC3aM1/xF4LdCiC+gOnA/5s1OnBDiWdRDO87eL/EdIABASvlrVD/FTUAl0Ad8fNLn9OLfS6PRaDTjwFtDOhqNRqMZJ1rwNRqNxkfQgq/RaDQ+ghZ8jUaj8RG04Gs0Go2PoAVfo9FofAQt+BqNRuMj/H/NJHuK1DrVTwAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], "source": [ - "# Do your own testing here ...\n", + "# There is an 𝑙 in the boundary conditions we assume should be a 1. (-1 pt)\n", + "a = 1\n", + "x = np.linspace(0, 1, 300)\n", + "t = np.linspace(0, 1, 300)\n", + "f = lambda x: np.sin(2*np.pi*x)\n", + "g = lambda x: 2*np.pi*np.sin(2*np.pi*x)\n", "\n", - "# YOUR CODE HERE\n", - "raise NotImplementedError()" + "w = pdeHyperbolic(a, x, t, f, g)\n", + "u = lambda x, t: np.sin(2*np.pi*x)*(np.cos(2*np.pi*t) + np.sin(2*np.pi*t))\n", + "\n", + "#plt.plot(t[1:], w[1:, ])\n", + "sub = np.arange(1, 300)\n", + "#plt.plot(t[sub], w[sub, sub])\n", + "#plt.plot(t, u(x, t))\n", + "print(np.max(w))\n", + "\n", + "# Now we hope that w[j, i] \\approx u(x[i], t[j]).\n", + "jlist = np.arange(1, 300)\n", + "ilist = np.arange(1, 300)\n", + "\n", + "plt.plot(t[jlist], u(x[ilist], t[jlist]))\n", + "plt.plot(t[jlist], w[jlist, ilist])" ] }, {