{ "cells": [ { "cell_type": "markdown", "id": "191f070c", "metadata": { "editable": true }, "source": [ "\n", "\n", "# Demo - 1D Poisson's equation\n", "**Mikael Mortensen** (email: `mikaem@math.uio.no`), Department of Mathematics, University of Oslo.\n", "\n", "Date: **April 13, 2018**\n", "\n", "**Summary.** This is a demonstration of how the Python module [shenfun](https://github.com/spectralDNS/shenfun) can be used to solve Poisson's\n", "equation with Dirichlet boundary conditions in one dimension. Spectral convergence, as\n", "shown in the figure below, is demonstrated.\n", "The demo is implemented in slightly more generic terms (more boundary conditions)\n", "in [poisson1D.py](https://github.com/spectralDNS/shenfun/blob/master/demo/poisson1D.py), and\n", "the numerical method is is described in more detail by J. Shen [shen1](#shen1) and [shen95](#shen95).\n", "\n", "\n", "\n", "\n", "\n", "

Figure 1: Convergence of 1D Poisson solvers for both Legendre and Chebyshev modified basis function.

\n", "" ] }, { "cell_type": "markdown", "id": "5c9e928b", "metadata": { "editable": true }, "source": [ "## Poisson's equation\n", "\n", "Poisson's equation is given as" ] }, { "cell_type": "markdown", "id": "88a950ff", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "\\nabla^2 u(x) = f(x) \\quad \\text{for }\\, x \\in (-1, 1), \\label{eq:poisson} \\tag{1}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "04929b21", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation} \n", "u(-1)=a, u(1)=b, \\notag\n", "\\label{_auto1} \\tag{2}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "6dff101a", "metadata": { "editable": true }, "source": [ "where $u(x)$ is the solution, $f(x)$ is a function and $a, b$ are two possibly\n", "non-zero constants.\n", "\n", "To solve Eq. ([1](#eq:poisson)) with the Galerkin method we need smooth continuously\n", "differentiable basis functions, $v_k$, that satisfy the given boundary conditions.\n", "And then we look for solutions like" ] }, { "cell_type": "markdown", "id": "40cbe4bb", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "u(x) = \\sum_{k=0}^{N-1} \\hat{u}_k v_k(x), \\label{eq:u} \\tag{3}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "75b794d4", "metadata": { "editable": true }, "source": [ "where $N$ is the size of the discretized problem,\n", "$\\{\\hat{u}_k\\}_{k=0}^{N-1}$ are the unknown expansion\n", "coefficients, and the function space is $\\text{span}\\{v_k\\}_{k=0}^{N-1}$.\n", "\n", "The basis functions of the function space can, for example, be constructed from\n", "[Chebyshev](https://en.wikipedia.org/wiki/Chebyshev_polynomials), $T_k(x)$, or\n", "[Legendre](https://en.wikipedia.org/wiki/Legendre_polynomials), $L_k(x)$, polynomials\n", "and we use the common notation $\\phi_k(x)$ to represent either one of them. It turns out that\n", "it is easiest to use basis functions with homogeneous Dirichlet boundary conditions" ] }, { "cell_type": "markdown", "id": "daad7929", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "v_k(x) = \\phi_k(x) - \\phi_{k+2}(x),\n", "\\label{_auto2} \\tag{4}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "49939654", "metadata": { "editable": true }, "source": [ "for $k=0, 1, \\ldots N-3$. This gives the function space\n", "$V^N_0 = \\text{span}\\{v_k(x)\\}_{k=0}^{N-3}$.\n", "We can then add two more linear basis functions (that belong to the kernel of Poisson's equation)" ] }, { "cell_type": "markdown", "id": "c2eb52a2", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "v_{N-2} = \\frac{1}{2}(\\phi_0 - \\phi_1), \n", "\\label{_auto3} \\tag{5}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "19b3c0c6", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation} \n", "v_{N-1} = \\frac{1}{2}(\\phi_0 + \\phi_1).\n", "\\label{_auto4} \\tag{6}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "f9753564", "metadata": { "editable": true }, "source": [ "which gives the inhomogeneous space $V^N = \\text{span}\\{v_k\\}_{k=0}^{N-1}$.\n", "With the two linear basis functions it is easy to see that the last two degrees\n", "of freedom, $\\hat{u}_{N-2}$ and $\\hat{u}_{N-1}$, now are given as" ] }, { "cell_type": "markdown", "id": "7ee320eb", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "u(-1) = \\sum_{k=0}^{N-1} \\hat{u}_k v_k(-1) = \\hat{u}_{N-2} = a,\n", "\\label{eq:dirichleta} \\tag{7} \n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "5bfbba11", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation} \n", "u(+1) = \\sum_{k=0}^{N-1} \\hat{u}_k v_k(+1) = \\hat{u}_{N-1} = b,\n", "\\label{eq:dirichletb} \\tag{8}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "a377dc2d", "metadata": { "editable": true }, "source": [ "and, as such, we only have to solve for $\\{\\hat{u}_k\\}_{k=0}^{N-3}$, just like\n", "for a problem with homogeneous boundary conditions (for homogeneous boundary condition\n", "we simply have $\\hat{u}_{N-2} = \\hat{u}_{N-1} = 0$).\n", "We now formulate a variational problem using the Galerkin method: Find $u \\in V^N$ such that" ] }, { "cell_type": "markdown", "id": "bc16986e", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "\\int_{-1}^1 \\nabla^2 u \\, v \\, w\\, dx = \\int_{-1}^1 f \\, v\\, w\\, dx \\quad \\forall v \\, \\in \\, V^N_0. \\label{eq:varform} \\tag{9}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "c56593b5", "metadata": { "editable": true }, "source": [ "Note that since we only have $N-3$ unknowns we are only using the homogeneous test\n", "functions from $V^N_0$.\n", "\n", "The weighted integrals, weighted by $w(x)$, are called inner products, and a\n", "common notation is" ] }, { "cell_type": "markdown", "id": "e7f282db", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "\\int_{-1}^1 u \\, v \\, w\\, dx = \\left( u, v\\right)_w.\n", "\\label{_auto5} \\tag{10}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "2aae170f", "metadata": { "editable": true }, "source": [ "The integral can either be computed exactly, or with quadrature. The advantage\n", "of the latter is that it is generally faster, and that non-linear terms may be\n", "computed just as quickly as linear. For a linear problem, it does not make much\n", "of a difference, if any at all. Approximating the integral with quadrature, we\n", "obtain" ] }, { "cell_type": "markdown", "id": "0597485b", "metadata": { "editable": true }, "source": [ "$$\n", "\\begin{align*}\n", "\\int_{-1}^1 u \\, v \\, w\\, dx &\\approx \\left( u, v \\right)_w^N, \\\\ \n", "&\\approx \\sum_{j=0}^{N-1} u(x_j) v(x_j) w_j,\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "id": "1820492f", "metadata": { "editable": true }, "source": [ "where $\\{w_j\\}_{j=0}^{N-1}$ are quadrature weights.\n", "The quadrature points $\\{x_j\\}_{j=0}^{N-1}$\n", "are specific to the chosen basis, and even within basis there are two different\n", "choices based on which quadrature rule is selected, either Gauss or Gauss-Lobatto.\n", "\n", "Inserting for test and trialfunctions, we get the following bilinear form and\n", "matrix $A=(a_{jk})\\in\\mathbb{R}^{N-2\\times N-2}$ for the Laplacian" ] }, { "cell_type": "markdown", "id": "0f03aa15", "metadata": { "editable": true }, "source": [ "$$\n", "\\begin{align*}\n", "\\left( \\nabla^2u, v \\right)_w^N &= \\left( \\nabla^2\\sum_{k=0}^{N-3}\\hat{u}_k v_{k}, v_j \\right)_w^N, \\quad j=0,1,\\ldots, N-3\\\\ \n", " &= \\sum_{k=0}^{N-3}\\left(\\nabla^2 v_{k}, v_j \\right)_w^N \\hat{u}_k, \\\\ \n", " &= \\sum_{k=0}^{N-3}a_{jk} \\hat{u}_k.\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "id": "08f3c800", "metadata": { "editable": true }, "source": [ "Note that the sum runs over $k=0, 1, \\ldots, N-3$ since\n", "the second derivatives of $v_{N-2}$ and $v_{N-1}$ are zero.\n", "The right hand side linear form and vector is computed as $\\tilde{f}_j = (f,\n", "v_j)_w^N$, for $j=0,1,\\ldots, N-3$, where a tilde is used because this is not\n", "a complete transform of the function $f$, but only an inner product.\n", "\n", "By defining the column vectors $\\boldsymbol{\\hat{u}}=(\\hat{u}_0, \\hat{u}_1, \\ldots, \\hat{u}_{N-3})^T$\n", "and $\\boldsymbol{\\tilde{f}}=(\\tilde{f}_0, \\tilde{f}_1, \\ldots, \\tilde{f}_{N-3})^T$\n", "we get the linear system of equations" ] }, { "cell_type": "markdown", "id": "65e0bcc0", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "A \\hat{\\boldsymbol{u}} = \\tilde{\\boldsymbol{f}}.\n", "\\label{_auto6} \\tag{11}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "998ba60a", "metadata": { "editable": true }, "source": [ "Now, when the expansion coefficients $\\boldsymbol{\\hat{u}}$ are found by\n", "solving this linear system, they may be\n", "transformed to real space $u(x)$ using ([3](#eq:u)), and here the contributions\n", "from $\\hat{u}_{N-2}$ and $\\hat{u}_{N-1}$ must be accounted for. Note that the matrix\n", "$A$ (different for Legendre or Chebyshev) has a very special structure that\n", "allows for a solution to be found very efficiently in order of $\\mathcal{O}(N)$\n", "operations, see [[shen1]](#shen1) and [[shen95]](#shen95). These solvers are implemented in\n", "shenfun for both bases." ] }, { "cell_type": "markdown", "id": "81e844aa", "metadata": { "editable": true }, "source": [ "## Method of manufactured solutions\n", "\n", "In this demo we will use the method of manufactured\n", "solutions to demonstrate spectral accuracy of the `shenfun` Dirichlet bases. To\n", "this end we choose an analytical function that satisfies the given boundary\n", "conditions:" ] }, { "cell_type": "markdown", "id": "a7b27de2", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "u_e(x) = \\sin(k\\pi x)(1-x^2) + a(1-x)/2 + b(1+x)/2, \\label{eq:u_e} \\tag{12}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "14a5e799", "metadata": { "editable": true }, "source": [ "where $k$ is an integer and $a$ and $b$ are constants. Now, feeding $u_e$ through\n", "the Laplace operator, we see that the last two linear terms disappear, whereas the\n", "first term results in" ] }, { "cell_type": "markdown", "id": "d45e3287", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", " \\nabla^2 u_e(x) = \\frac{d^2 u_e}{dx^2}, \n", "\\label{_auto7} \\tag{13}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "d0f28dd7", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation} \n", " = -4k \\pi x \\cos(k\\pi x) - 2\\sin(k\\pi x) - k^2 \\pi^2 (1 -\n", "x^2) \\sin(k \\pi x). \\label{eq:solution} \\tag{14}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "91f8eb49", "metadata": { "editable": true }, "source": [ "Now, setting $f_e(x) = \\nabla^2 u_e(x)$ and solving for $\\nabla^2 u(x) = f_e(x)$,\n", "we can compare the numerical solution $u(x)$ with the analytical solution $u_e(x)$\n", "and compute error norms." ] }, { "cell_type": "markdown", "id": "fc08f897", "metadata": { "editable": true }, "source": [ "## Implementation" ] }, { "cell_type": "markdown", "id": "dc126f5a", "metadata": { "editable": true }, "source": [ "### Preamble\n", "\n", "We will solve Poisson's equation using the [shenfun](https://github.com/spectralDNS/shenfun) Python module. The first thing needed\n", "is then to import some of this module's functionality\n", "plus some other helper modules, like [Numpy](https://numpy.org) and [Sympy](https://sympy.org):" ] }, { "cell_type": "code", "execution_count": 1, "id": "659b5c0d", "metadata": { "editable": true, "execution": { "iopub.execute_input": "2026-06-03T08:49:42.681553Z", "iopub.status.busy": "2026-06-03T08:49:42.680941Z", "iopub.status.idle": "2026-06-03T08:49:43.086984Z", "shell.execute_reply": "2026-06-03T08:49:43.086733Z" }, "tags": [ "thebe-init" ] }, "outputs": [], "source": [ "from shenfun import inner, div, grad, TestFunction, TrialFunction, Function, \\\n", " project, Dx, Array, FunctionSpace, dx\n", "import numpy as np\n", "from sympy import symbols, cos, sin, exp, lambdify" ] }, { "cell_type": "markdown", "id": "cfc9c744", "metadata": { "editable": true }, "source": [ "We use `Sympy` for the manufactured solution and `Numpy` for testing." ] }, { "cell_type": "markdown", "id": "40d4aa04", "metadata": { "editable": true }, "source": [ "### Manufactured solution\n", "\n", "The exact solution $u_e(x)$ and the right hand side $f_e(x)$ are created using\n", "`Sympy` as follows" ] }, { "cell_type": "code", "execution_count": 2, "id": "38c3edf0", "metadata": { "editable": true, "execution": { "iopub.execute_input": "2026-06-03T08:49:43.088919Z", "iopub.status.busy": "2026-06-03T08:49:43.088663Z", "iopub.status.idle": "2026-06-03T08:49:43.100902Z", "shell.execute_reply": "2026-06-03T08:49:43.100680Z" }, "tags": [ "thebe-init" ] }, "outputs": [], "source": [ "a = -1\n", "b = 1\n", "k = 4\n", "x = symbols(\"x\")\n", "ue = sin(k*np.pi*x)*(1-x**2) + a*(1 - x)/2. + b*(1 + x)/2.\n", "fe = ue.diff(x, 2)" ] }, { "cell_type": "markdown", "id": "a66d4067", "metadata": { "editable": true }, "source": [ "These solutions are now valid for a continuous domain. The next step is thus to\n", "discretize, using a discrete mesh $\\{x_j\\}_{j=0}^{N-1}$ and a finite number of\n", "basis functions.\n", "\n", "Note that it is not mandatory to use `Sympy` for the manufactured solution. Since the\n", "solution is known ([14](#eq:solution)), we could just as well simply use `Numpy`\n", "to compute $f_e$ at $\\{x_j\\}_{j=0}^{N-1}$. However, with `Sympy` it is much\n", "easier to experiment and quickly change the solution." ] }, { "cell_type": "markdown", "id": "e2cfe44f", "metadata": { "editable": true }, "source": [ "### Discretization\n", "\n", "We create a basis with a given number of basis functions, and extract the computational\n", "mesh from the basis itself" ] }, { "cell_type": "code", "execution_count": 3, "id": "9328ec31", "metadata": { "editable": true, "execution": { "iopub.execute_input": "2026-06-03T08:49:43.102247Z", "iopub.status.busy": "2026-06-03T08:49:43.102161Z", "iopub.status.idle": "2026-06-03T08:49:43.113237Z", "shell.execute_reply": "2026-06-03T08:49:43.113007Z" }, "tags": [ "thebe-init" ] }, "outputs": [], "source": [ "N = 32\n", "SD = FunctionSpace(N, 'Chebyshev', bc=(a, b))\n", "#SD = FunctionSpace(N, 'Legendre', bc=(a, b))" ] }, { "cell_type": "markdown", "id": "db6ef1ac", "metadata": { "editable": true }, "source": [ "Note that we can either choose a Legendre or a Chebyshev basis." ] }, { "cell_type": "markdown", "id": "37f9d1b5", "metadata": { "editable": true }, "source": [ "### Variational formulation\n", "\n", "The variational problem ([9](#eq:varform)) can be assembled using `shenfun`'s\n", "[TrialFunction](https://shenfun.readthedocs.io/en/latest/shenfun.forms.html#shenfun.forms.arguments.TrialFunction), [TestFunction](https://shenfun.readthedocs.io/en/latest/shenfun.forms.html#shenfun.forms.arguments.TestFunction) and [inner()](https://shenfun.readthedocs.io/en/latest/shenfun.forms.html#shenfun.forms.inner.inner) functions." ] }, { "cell_type": "code", "execution_count": 4, "id": "d3f0a739", "metadata": { "editable": true, "execution": { "iopub.execute_input": "2026-06-03T08:49:43.114498Z", "iopub.status.busy": "2026-06-03T08:49:43.114425Z", "iopub.status.idle": "2026-06-03T08:49:43.123259Z", "shell.execute_reply": "2026-06-03T08:49:43.123028Z" }, "tags": [ "thebe-init" ] }, "outputs": [], "source": [ "u = TrialFunction(SD)\n", "v = TestFunction(SD)\n", "# Assemble left hand side matrix\n", "A = inner(v, div(grad(u)))\n", "# Assemble right hand side\n", "fj = Array(SD, buffer=fe)\n", "f_hat = Function(SD)\n", "f_hat = inner(v, fj, output_array=f_hat)" ] }, { "cell_type": "markdown", "id": "5ecdef9b", "metadata": { "editable": true }, "source": [ "Note that the `sympy` function `fe` can be used to initialize the [Array](https://shenfun.readthedocs.io/en/latest/shenfun.forms.html#shenfun.forms.arguments.Array)\n", "`fj`. We wrap this Numpy array in an [Array](https://shenfun.readthedocs.io/en/latest/shenfun.forms.html#shenfun.forms.arguments.Array) class\n", "(`fj = Array(SD, buffer=fe)`), because an Array\n", "is required as input to the [inner()](https://shenfun.readthedocs.io/en/latest/shenfun.forms.html#shenfun.forms.inner.inner) function." ] }, { "cell_type": "markdown", "id": "ee8d601e", "metadata": { "editable": true }, "source": [ "### Solve linear equations\n", "\n", "Finally, solve linear equation system and transform solution from spectral\n", "$\\boldsymbol{\\hat{u}}$ vector to the real space $\\{u(x_j)\\}_{j=0}^{N-1}$\n", "and then check how the solution corresponds with the exact solution $u_e$.\n", "To this end we compute the $L_2$-errornorm using the `shenfun` function\n", "[dx()](https://shenfun.readthedocs.io/en/latest/shenfun.utilities.html#shenfun.utilities.dx)" ] }, { "cell_type": "code", "execution_count": 5, "id": "3efea157", "metadata": { "editable": true, "execution": { "iopub.execute_input": "2026-06-03T08:49:43.124585Z", "iopub.status.busy": "2026-06-03T08:49:43.124506Z", "iopub.status.idle": "2026-06-03T08:49:43.137796Z", "shell.execute_reply": "2026-06-03T08:49:43.137606Z" }, "tags": [ "thebe-init" ] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Error=2.3565371191241489e-10\n" ] } ], "source": [ "u_hat = A.solve(f_hat)\n", "uj = SD.backward(u_hat)\n", "ua = Array(SD, buffer=ue)\n", "print(\"Error=%2.16e\" %(np.sqrt(dx((uj-ua)**2))))" ] }, { "cell_type": "markdown", "id": "3c8908a6", "metadata": { "editable": true }, "source": [ "### Convergence test\n", "\n", "To do a convergence test we will now create a function `main`, that takes the\n", "number of quadrature points as parameter, and prints out\n", "the error." ] }, { "cell_type": "code", "execution_count": 6, "id": "c69c9d1e", "metadata": { "editable": true, "execution": { "iopub.execute_input": "2026-06-03T08:49:43.139045Z", "iopub.status.busy": "2026-06-03T08:49:43.138971Z", "iopub.status.idle": "2026-06-03T08:49:43.141093Z", "shell.execute_reply": "2026-06-03T08:49:43.140875Z" }, "tags": [ "thebe-init" ] }, "outputs": [], "source": [ "def main(N, family='Chebyshev'):\n", " SD = FunctionSpace(N, family=family, bc=(a, b))\n", " u = TrialFunction(SD)\n", " v = TestFunction(SD)\n", "\n", " # Get f on quad points\n", " fj = Array(SD, buffer=fe)\n", "\n", " # Compute right hand side of Poisson's equation\n", " f_hat = Function(SD)\n", " f_hat = inner(v, fj, output_array=f_hat)\n", "\n", " # Get left hand side of Poisson's equation\n", " A = inner(v, div(grad(u)))\n", "\n", " f_hat = A.solve(f_hat)\n", " uj = SD.backward(f_hat)\n", "\n", " # Compare with analytical solution\n", " ua = Array(SD, buffer=ue)\n", " l2_error = np.linalg.norm(uj-ua)\n", " return l2_error" ] }, { "cell_type": "markdown", "id": "8199b9c0", "metadata": { "editable": true }, "source": [ "For example, we find the error of a Chebyshev discretization\n", "using 12 quadrature points as" ] }, { "cell_type": "code", "execution_count": 7, "id": "e2444fa7", "metadata": { "editable": true, "execution": { "iopub.execute_input": "2026-06-03T08:49:43.142236Z", "iopub.status.busy": "2026-06-03T08:49:43.142166Z", "iopub.status.idle": "2026-06-03T08:49:43.149893Z", "shell.execute_reply": "2026-06-03T08:49:43.149669Z" }, "tags": [ "thebe-init" ] }, "outputs": [ { "data": { "text/plain": [ "1.1654916318123927" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "main(12, 'Chebyshev')" ] }, { "cell_type": "markdown", "id": "41161308", "metadata": { "editable": true }, "source": [ "To get the convergence we call `main` for a list\n", "of $N=[12, 16, \\ldots, 48]$, and collect the errornorms in\n", "arrays to be plotted. The error can be plotted using\n", "[matplotlib](https://matplotlib.org), and the generated\n", "figure is also shown in this demos summary." ] }, { "cell_type": "code", "execution_count": 8, "id": "be4c522a", "metadata": { "editable": true, "execution": { "iopub.execute_input": "2026-06-03T08:49:43.151092Z", "iopub.status.busy": "2026-06-03T08:49:43.151020Z", "iopub.status.idle": "2026-06-03T08:49:43.870264Z", "shell.execute_reply": "2026-06-03T08:49:43.870018Z" }, "tags": [ "thebe-init" ] }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAikAAAGHCAYAAAB1bcIdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABjhUlEQVR4nO3dd1yV5f/H8ddhunFgKIqznLhCJSy3gjhS09Qyd5maqamZVra/OSpHmZoNR2na0sxMpVS03IOsNFvkSFypgAsErt8f5+fJIw5A4D7A+/l43A89132dc973udHz4b6v+7ptxhiDiIiIiItxszqAiIiIyLWoSBERERGXpCJFREREXJKKFBEREXFJKlJERETEJalIEREREZekIkVERERckooUERERcUkqUkRERMQlqUgRl7Znzx769etHxYoVyZcvH4UKFeLOO+9k8uTJnDp1yup4ksmWLFlCzZo1yZ8/PzabjaioqGv2W79+PTabzbG4u7vj5+fH/fffz759+9L9vn379qVChQq3Fj6XyKmfxYoVK+jduze1atXC09MTm812zX5X/+x4eXlRsmRJ7r77bp555hkOHDiQzcnlRlSkiMt69913CQoKYvv27Tz55JOsWrWKpUuXcv/99zN79mwGDBhgdUTJRCdOnKBXr15UrlyZVatWsXnzZqpUqXLD57z66qts3ryZdevW8dRTTxEREcHdd9/NP//8k673Hj9+PEuXLr2V+GKxpUuXsmXLFmrUqEGdOnVu2v/Kn53333+fZs2a8cEHH1C9enUWLlyYDYklTYyIC9q0aZNxd3c3bdq0MRcvXky1PiEhwXz55ZcWJMs8SUlJ19y2vOr77783gFmyZMlN+65bt84A5tNPP3Vqf//99w1gXnnllayKmev16dPHlC9fPtvf9/z587f0/OTkZMffH3vsMXO9r7fr/ewYY8y///5r6tWrZzw8PMyePXtuKY9kDh1JEZf06quvYrPZmDNnDt7e3qnWe3l5ce+99zoep6SkMHnyZKpVq4a3tze33XYbvXv35vDhw07Pa9asGYGBgWzfvp3GjRtToEABKlWqxMSJE0lJSQHsv9F7eXkxfvz4VO/766+/YrPZePPNNx1tR48e5dFHH6Vs2bJ4eXlRsWJFXnzxRZKSkhx9/v77b2w2G5MnT+aVV16hYsWKeHt7s27dOgC+/PJLateujbe3N5UqVWL69Om88MILqQ5ZG2OYOXMmdevWJX/+/BQrVoyuXbvy119/pXs7Lztz5gyjRo2iUqVKjs+ubdu2/Prrr44+iYmJvPLKK47Pt2TJkvTr148TJ05cewdeZfny5YSEhFCgQAEKFy5M69at2bx5s2N93759ueeeewDo3r07NpuNZs2apem1r3TXXXcBOA7Zp/Xn4lqnOD799FOCg4Px8fFxfH79+/d3rE9JSeGVV16hatWq5M+fn6JFi1K7dm2mT5/u9Drff/89LVu2pHDhwhQoUIBGjRrx9ddfO/WZN28eNpuNdevWMXjwYHx9fSlRogT33XcfR44cuel2//XXX/To0QN/f3+8vb3x8/OjZcuWTqfL0vpZXK1evXo0btw4VXtycjJlypThvvvuc7Sl9eekQoUKtG/fni+++IJ69eqRL18+XnzxReDmn/v1uLnd+tdZ8eLFeeedd0hKSmLq1Km3/HqSCayukkSulpSUZAoUKGCCg4PT/JyBAwcawAwdOtSsWrXKzJ4925QsWdIEBASYEydOOPo1bdrUlChRwtxxxx1m9uzZJiIiwgwZMsQAZv78+Y5+nTt3NgEBAU6/nRljzJgxY4yXl5c5efKkMcaYmJgYExAQYMqXL2/eeecd8+2335qXX37ZeHt7m759+zqeFx0dbQBTpkwZ07x5c/PZZ5+ZNWvWmOjoaPPNN98YNzc306xZM7N06VLz6aefmuDgYFOhQoVUvw0+8sgjxtPT04waNcqsWrXKLFq0yFSrVs34+fmZo0ePpns74+LiTM2aNU3BggXNSy+9ZFavXm0+//xzM3z4cLN27VpjjP031DZt2piCBQuaF1980URERJj33nvPlClTxtSoUeOmvwEvXLjQACY0NNQsW7bMLFmyxAQFBRkvLy+zceNGY4wxf/zxh3n77bcNYF599VWzefNm88svv1z3Na/32/CXX35pAPP0008bY9L+c3H10YNNmzYZm81mevToYVauXGnWrl1r5s6da3r16uXoM2HCBOPu7m6ef/55891335lVq1aZadOmmRdeeMHRZ/369cbT09MEBQWZJUuWmGXLlpnQ0FBjs9nM4sWLHf3mzp1rAFOpUiXz+OOPm9WrV5v33nvPFCtWzDRv3vyGn68xxlStWtXcfvvt5sMPPzSRkZHm888/N6NGjTLr1q1z9MnoZzF9+nQDmN9++83pPVeuXGkAs3z5cmNM+n5Oypcvb0qXLm0qVapkPvjgA7Nu3Tqzbdu2NH3uaZHRIymXlS5d2lSuXDld7ylZQ0WKuJyjR48awPTo0SNN/fft22cAM2TIEKf2rVu3On1hGWP/8gbM1q1bnfrWqFHDhIWFOR4vX77cAGbNmjWOtqSkJOPv72+6dOniaHv00UdNoUKFzIEDB5xe7/XXXzeA44v2cpFSuXJlk5iY6NS3QYMGJiAgwCQkJDja4uPjTYkSJZz+o928ebMBzBtvvOH0/EOHDpn8+fObMWPGpHs7X3rpJQOYiIgIcz0ff/yxAcznn3/u1L59+3YDmJkzZ173ucnJycbf39/UqlXLqeCLj483t912m2nUqJGjLS1fHlf3XbJkibl06ZI5f/682bBhg7n99tuNu7u7+fHHH9P1c3H1F/Pl/XfmzJnrZmjfvr2pW7fuDXPedddd5rbbbjPx8fGOtqSkJBMYGGjKli1rUlJSjDH/FSlXZ508ebIBTExMzHXf4+TJkwYw06ZNu26fW/ksTp48aby8vJz6GGNMt27djJ+fn7l06ZIxJn0/J+XLlzfu7u5m//79Tn3T8rmnxa0WKcHBwSZ//vy3lEEyh073SI53+ZRJ3759ndobNmxI9erV+e6775zaS5UqRcOGDZ3aateu7TSqPzw8nFKlSjF37lxH2+rVqzly5IjToecVK1bQvHlz/P39SUpKcizh4eEAREZGOr3Pvffei6enp+PxuXPn2LFjB506dcLLy8vRXqhQITp06OD03BUrVmCz2XjooYec3qtUqVLUqVOH9evXp3s7v/nmG6pUqUKrVq24nhUrVlC0aFE6dOjg9L5169alVKlSqd73Svv37+fIkSP06tXL6XB8oUKF6NKlC1u2bOH8+fPXff7NdO/eHU9PTwoUKECTJk1ITk7ms88+o3bt2un+ubhSgwYNAOjWrRuffPLJNQfiNmzYkB9//JEhQ4awevVq4uLinNafO3eOrVu30rVrVwoVKuRod3d3p1evXhw+fJj9+/c7PefKU5hg31/ADa84KV68OJUrV+a1115jypQp7N69O9UpvVv5LEqUKEGHDh2YP3++43VPnz7Nl19+Se/evfHw8ADS/3NSu3btVAOj0/K5ZwdjjCXvK6mpSBGX4+vrS4ECBYiOjk5T/3///ReA0qVLp1rn7+/vWH9ZiRIlUvXz9vbmwoULjsceHh706tWLpUuXcubMGcA+bqB06dKEhYU5+h07doyvvvoKT09Pp6VmzZoAnDx50ul9rs54+vRpjDH4+fmlynR127Fjxxx9r36/LVu2pHqvtGzniRMnKFu2bKp+V7/vmTNn8PLySvW+R48eTfW+V7rZvklJSeH06dM3fP8bmTRpEtu3b2fXrl0cPHiQv/76i06dOqXpva/+ubhSkyZNWLZsGUlJSfTu3ZuyZcsSGBjIxx9/7Ogzbtw4Xn/9dbZs2UJ4eDglSpSgZcuW7NixA/hv317v/a/MeNnV++zyeKwr99nVbDYb3333HWFhYUyePJk777yTkiVLMmzYMOLj42/5swDo378///zzDxEREQB8/PHHJCQkOBU96f05uVaWtHzu2eHgwYOOfSTW8rA6gMjV3N3dadmyJd988w2HDx++6Zfo5f/YY2JiUvU9cuQIvr6+GcrRr18/XnvtNRYvXkz37t1Zvnw5I0aMwN3d3dHH19eX2rVr87///e+ar3H1f3RXD4QtVqwYNpuNY8eOpXru0aNHnR77+vpis9nYuHHjNQcTX6vtZkqWLHnTgZOXB3GuWrXqmusLFy583edeuW+uduTIEdzc3ChWrFg6EjurVKkS9evXv+l7Z+TnomPHjnTs2JGEhAS2bNnChAkTePDBB6lQoQIhISF4eHgwcuRIRo4cyZkzZ/j22295+umnCQsL49ChQxQrVgw3N7frbjuQ4Z/Nq5UvX573338fgN9++41PPvmEF154gcTERGbPnn3Ln0VYWBj+/v7MnTuXsLAw5s6dS3BwMDVq1HD0Se/PyfXmMbnZ557Vtm3bxtGjRzXFgYvQkRRxSePGjcMYwyOPPEJiYmKq9ZcuXeKrr74CoEWLFgB89NFHTn22b9/Ovn37aNmyZYYyVK9eneDgYObOncuiRYtISEigX79+Tn3at2/Pzz//TOXKlalfv36q5Wa/jRUsWJD69euzbNkyp+08e/YsK1asSPVexhj++eefa75XrVq10r2N4eHh/Pbbb6xdu/a6fdq3b8+///5LcnLyNd+3atWq131u1apVKVOmDIsWLXI6hH7u3Dk+//xzxxU/WSGzfi68vb1p2rQpkyZNAmD37t2p+hQtWpSuXbvy2GOPcerUKf7++28KFixIcHAwX3zxhdORkJSUFD766CPKli1703lgMqJKlSo8++yz1KpVi127dgG3/llcPkW1bNkyNm7cyI4dO1JdcXMrPyfXkpbPPbOdOnWKQYMG4enpyRNPPJHl7yc3pyMp4pJCQkKYNWsWQ4YMISgoiMGDB1OzZk0uXbrE7t27mTNnDoGBgXTo0IGqVasycOBA3nrrLdzc3AgPD+fvv/9m/PjxBAQE3NJ/Nv379+fRRx/lyJEjNGrUKNV/tC+99BIRERE0atSIYcOGUbVqVS5evMjff//NypUrmT179k2PBL300ku0a9eOsLAwhg8fTnJyMq+99hqFChVymlX37rvvZuDAgfTr148dO3bQpEkTChYsSExMDN9//z21atVi8ODB6dq+ESNGsGTJEjp27MjYsWNp2LAhFy5cIDIykvbt29O8eXN69OjBwoULadu2LcOHD6dhw4Z4enpy+PBh1q1bR8eOHencufM1X9/NzY3JkyfTs2dP2rdvz6OPPkpCQgKvvfYaZ86cYeLEienKmx638nPx3HPPcfjwYVq2bEnZsmU5c+YM06dPx9PTk6ZNmwLQoUMHAgMDqV+/PiVLluTAgQNMmzaN8uXLc8cddwAwYcIEWrduTfPmzRk9ejReXl7MnDmTn3/+mY8//vi6RxPSY8+ePQwdOpT777+fO+64Ay8vL9auXcuePXsYO3bsLX8Wl/Xv359Jkybx4IMPkj9/frp37+60/lZ+Ti5Ly+d+PQcOHGD79u0A/PnnnwB89tlngP2S56uPuP3+++9s2bKFlJQU/v33X7Zu3cr7779PXFwcCxYscJyyFYtZOGhX5KaioqJMnz59TLly5YyXl5cpWLCgqVevnnnuuefM8ePHHf2Sk5PNpEmTTJUqVYynp6fx9fU1Dz30kDl06JDT6zVt2tTUrFkz1ftcbwKr2NhYkz9/fgOYd99995oZT5w4YYYNG2YqVqxoPD09TfHixU1QUJB55plnzNmzZ40x/13d89prr13zNZYuXWpq1aplvLy8TLly5czEiRPNsGHDTLFixVL1/eCDD0xwcLApWLCgyZ8/v6lcubLp3bu32bFjR4a28/Tp02b48OGmXLlyxtPT09x2222mXbt25tdff3X0uXTpknn99ddNnTp1TL58+UyhQoVMtWrVzKOPPmp+//33a27TlZYtW2aCg4NNvnz5TMGCBU3Lli3NDz/84NQnI1f33KxvWn8urv5cVqxYYcLDw02ZMmWMl5eXue2220zbtm0dl0wbY8wbb7xhGjVqZHx9fR37bcCAAebvv/92eu2NGzeaFi1aOPbXXXfdZb766iunPpev7tm+ffs1t/PKS4mvduzYMdO3b19TrVo1U7BgQVOoUCFTu3ZtM3XqVJOUlHTLn8WVGjVqZADTs2fPa65P689J+fLlTbt27VI9Py2f+/Vc/gyvtfTp08fR7/Jnennx8PAwJUqUMCEhIebpp59Otf/EWjZjNIxZxNVcunSJunXrUqZMGdasWWN1HBERS+h0j4gLGDBgAK1bt6Z06dIcPXqU2bNns2/fvlSzl4qI5CUqUkRcQHx8PKNHj+bEiRN4enpy5513snLlyhvOXyIiktvpdI+IiIi4JF2CLCIiIi5JRYqIiIi4JBUpIiIi4pI0cDaDUlJSOHLkCIULF86UCZlERETyCmMM8fHx+Pv7O9189GoqUjLoyJEjBAQEWB1DREQkxzp06NANZ+VWkZJBl2+WdejQIYoUKWJxGhERkZwjLi6OgICAG96gFFSkZNjlUzxFihRRkSIiIpIBNxsuoYGzIiIi4pJUpIiIiIhLUpEiIiIiLkljUkREJMsZY0hKSiI5OdnqKJIN3N3d8fDwuOUpOlSkiIhIlkpMTCQmJobz589bHUWyUYECBShdujReXl4Zfg0VKSIikmVSUlKIjo7G3d0df39/vLy8NAFmLmeMITExkRMnThAdHc0dd9xxwwnbbkRFioiIZJnExERSUlIICAigQIECVseRbJI/f348PT05cOAAiYmJ5MuXL0Ovo4GzIiKS5TL6m7TkXJmxz/VT40JSUqxOICIi4jpUpLiIg/svEFD4NKO6HSRqVwrGWJ1IRETEWipSXMTHz//KkfPFmPJpOeoFuVHL7zgTnzrNoUNWJxMRkZykb9++dOrUyeoYmUJFios4ufNvPEl0PP7lxG2Mm1yM8uVSaF7jKB/MSiA21sKAIiJ5TG76ss+pVKS4iNd+Cufoe18zu8ab3MNGR7vBjfX7SjFgiDelSiTSvdVJvlpuuHTJwrAiIpJjJSYm3ryTi1CR4iry5aP4gM48+sswNh6qyJ+jZvJy8alUYb+jy8VkLz75zpd7O9rwL3aeof3OsXUrGr8iIpLN9u7dS9u2bSlUqBB+fn706tWLkydPOtbHx8fTs2dPChYsSOnSpZk6dSrNmjVjxIgRjj6JiYmMGTOGMmXKULBgQYKDg1m/fr1j/bx58yhatCirV6+mevXqFCpUiDZt2hATE+Pok5yczMiRIylatCglSpRgzJgxmKu+FJo1a8bQoUMZOXIkvr6+tG7dOk3b4ApUpLiismWp9PoQnj05gl83nGBbh5d53HMWJTnu6HLyXAHenleQu+6CqmXieem5JP7808LMIiLpUb8+lC2bvUv9+pkSPSYmhqZNm1K3bl127NjBqlWrOHbsGN26dXP0GTlyJD/88APLly8nIiKCjRs3smvXLqfX6devHz/88AOLFy9mz5493H///bRp04bff//d0ef8+fO8/vrrfPjhh2zYsIGDBw8yevRox/o33niDDz74gPfff5/vv/+eU6dOsXTp0lSZ58+fj4eHBz/88APvvPNOmrbBJRjJkNjYWAOY2NjY7HnD+HiT+N58syLwKdODRSYf5439GIrzElI73sx8O8WcPJk9sUREbuTChQtm79695sKFC84rypRJ/R9YVi9lyqQre58+fUzHjh1TtY8fP96EhoY6tR06dMgAZv/+/SYuLs54enqaTz/91LH+zJkzpkCBAmb48OHGGGP++OMPY7PZzD///OP0Oi1btjTjxo0zxhgzd+5cA5g//vjDsf7tt982fn5+jselS5c2EydOdDy+dOmSKVu2rFPupk2bmrp166ZrGzLDdfe9Sft3qGaczSkKFcJzQG/aDYB2f/5J3JwpfPHeKT481ZZ1NMf8/0GxzXsKsfkxGD4smbatk3hogDft20MGJ/sTEckapUrl2PfcuXMn69ato1ChQqnW/fnnn1y4cIFLly7RsGFDR7uPjw9Vq1Z1PN61axfGGKpUqeL0/ISEBEqUKOF4XKBAASpXrux4XLp0aY4ftx9Vj42NJSYmhpCQEMd6Dw8P6tevn+qUT/2rjiLdbBuuzmUVFSk5UeXKFJn0DH0npNB37VoOzxzBoq8K82FSD36mFgCXkt35cpU7X64Cn4KXuL+7O736uHHPPaCJH0XEcjt2WJ0gw1JSUujQoQOTJk1Kta506dKO0zVX36PoysIhJSUFd3d3du7cibu7u1O/KwsHT09Pp3U2my1VAZIWBQsWTNc2uAoVKTmZmxu0akXZVq0YExvLmCVL+HHGZD76qTYL6UkM/gDEnvPkvQ/gvQ+gvH8iPft68dBDUL26xflFRHKgO++8k88//5wKFSrg4ZH6a7Ry5cp4enqybds2AgICAIiLi+P333+nadOmANSrV4/k5GSOHz9O48aNM5TDx8eH0qVLs2XLFpo0aQJAUlISO3fu5M4777ylbXAV+p06t/DxgYEDqbPnQ17b255Do98kolg3ejOfgpx1dDtwxItXX4UaNaB+vSSmTYNjx6yLLSLiymJjY4mKinJaHn30UU6dOsUDDzzAtm3b+Ouvv1izZg39+/cnOTmZwoUL06dPH5588knWrVvHL7/8Qv/+/XFzc3McXalSpQo9e/akd+/efPHFF0RHR7N9+3YmTZrEypUr05xv+PDhTJw4kaVLl/Lrr78yZMgQzpw5c9PnPfbYYzfcBlehIiU3ql4d99cm0ur4Iuav8OVYx0dZ6N6bcFbiTpKj284oD554Asr4pxDexrBoEZw7Z2FuEREXs379eurVq+e0PPfcc/zwww8kJycTFhZGYGAgw4cPx8fHx3FTvSlTphASEkL79u1p1aoVd999N9WrV3e6G/DcuXPp3bs3o0aNomrVqtx7771s3brVcfQlLUaNGkXv3r3p27cvISEhFC5cmM6dO9/0ef7+/jfdBldgMxk5uSXExcXh4+NDbGwsRYoUsTrOzZ08CYsWcWzOlyz+JZAP6cVOUl+OV6hgCvd1ceOhh6BFC7jqVKmISLpcvHiR6OhoKlas6PQFndecO3eOMmXK8MYbbzBgwACr42SLG+37tH6Huk65JFnL1xeGDcPv5+8YvrsfO4Z9yF6fEJ7mf5TjgKPb2XNuLFgAoaEQUDaF0aMhKkoTxomIpMfu3bv5+OOP+fPPP9m1axc9e/YEoGPHjhYny1lUpORFdevC9OlUPx7J/z6rRnS7x4l0a84jzMGHM45uMUfdeOMNqFcPatc2TJoEhw9bllpEJEd5/fXXqVOnDq1ateLcuXNs3LgRX19fq2PlKDrdk0E57nTPzcTEwEcfcfH9hXy9vzIf0ouVtOUSXk7dbDZDs2Y2evWCLl0gN2y6iGQdne7Ju3S6RzJP6dLw5JPk27ebLlvGsOzRVcQUrspMBtOIHxzdjLGxbh307w9+foYePWDFCnTDQxERyXQqUsSZzQbBwTB7NiWO7WXwwsb80OoF/uB2XuQ5bue/e0pcvGhjyRLo0AH8/Q1jxkBcnIXZRUQkV8nTRUrnzp0pVqwYXbt2tTqKa8qfHx58ECIiqHxgLc+95MlvFduwhWAeYwYl+O9umSdP2njtNahRw7B8uYWZRUQk18jTRcqwYcNYsGCB1TFyhnLlYPx4bH/8TnDka8zou5OYAreznA50YwneXATgn39sdOwI3brB0aMWZxYRkRwtTxcpzZs3p3DhwlbHyFnc3KBJE5g7F8+jh+jwwX0safw2+6hOKKsd3T79FKpXN7z3ni5fFhGRjHHZImXDhg106NABf39/bDYby5YtS9Vn5syZjlHDQUFBbNy4MfuD5mWFC0O/frBhAxW3LmFVzdF8yEP4cgKAM2dsPPIING8Ov/1mcVYREclxXLZIOXfuHHXq1GHGjBnXXL9kyRJGjBjBM888w+7du2ncuDHh4eEcPHjQ0ScoKIjAwMBUy5EjR9KdJyEhgbi4OKdFrtCwIbZdO3nopars86hNL/47jRYZaZ9n5X//g8RECzOKiGSy6/0SnR7NmjVjxIgRmZLnWl544QXq1q2bZa+flVy2SAkPD+eVV17hvvvuu+b6KVOmMGDAAB5++GGqV6/OtGnTCAgIYNasWY4+O3fu5Oeff061+Pv7pzvPhAkT8PHxcSzpubdCnuHlBePH4/vjdywImc1qQqnIXwAkJNh49lkICoKtWy3OKSKSRkePHuXxxx+nUqVKeHt7ExAQQIcOHfjuu++sjpYnuGyRciOJiYns3LmT0NBQp/bQ0FA2bdqUJe85btw4YmNjHcuhQ4ey5H1yhRo1YONGQqe356cCdzGa13DDflfNn3+GkBDDsGEQH29xThGRG/j7778JCgpi7dq1TJ48mZ9++olVq1bRvHlzHnvsMavj5Qk5skg5efIkycnJ+Pn5ObX7+flxNB2XlISFhXH//fezcuVKypYty/bt26/b19vbmyJFijgtcgPu7jBsGAV/2cZrod+ynQbUYxdgnxDurbegZk3DihUW5xQRuY4hQ4Zgs9nYtm0bXbt2pUqVKtSsWZORI0eyZcsWR7+TJ0/SuXNnChQowB133MHyq+Zh2Lt3L23btqVQoUL4+fnRq1cvTp486dQnKSmJoUOHUrRoUUqUKMGzzz7L5QnhX3rpJWrVqpUqX1BQEM899xxgv1tzw4YNKViwIEWLFuXuu+/mwIEDTv0//PBDKlSogI+PDz169CD+it8UjTFMnjyZSpUqkT9/furUqcNnn30GQEpKCmXLlmX27NlOr7dr1y5sNht//fVXej/aNMuRRcplNpvN6bExJlXbjaxevZoTJ05w/vx5Dh8+TIMGDTI7olSoAKtWcef8EWwr1obXGE1+zgNw6JCNDh2gRw84dszamCKSverXh7Jls3epn/rG79d16tQpVq1axWOPPUbBggVTrS9atKjj7y+++CLdunVjz549tG3blp49e3Lq1CkAYmJiaNq0KXXr1mXHjh2sWrWKY8eO0a1bN6fXmz9/Ph4eHmzdupU333yTqVOn8t577wHQv39/9u7d6/SL9J49e9i9ezd9+/YlKSmJTp060bRpU/bs2cPmzZsZOHCg0/fhn3/+ybJly1ixYgUrVqwgMjKSiRMnOtY/++yzzJ07l1mzZvHLL7/wxBNP8NBDDxEZGYmbmxs9evRg4cKFTpkXLVpESEgIlSpVSvsHm14mBwDM0qVLHY8TEhKMu7u7+eKLL5z6DRs2zDRp0iRbMsXGxhrAxMbGZsv75QpHjxrTvbv5k4qmFWuM/eJk+1KsWIp5/31jUlKsDikimenChQtm79695sKFC07tZcoYp/8DsmMpUybtubdu3WqAVN8zVwPMs88+63h89uxZY7PZzDfffGOMMWb8+PEmNDTU6TmHDh0ygNm/f78xxpimTZua6tWrm5Qr/gN86qmnTPXq1R2Pw8PDzeDBgx2PR4wYYZo1a2aMMebff/81gFm/fv01Mz7//POmQIECJi4uztH25JNPmuDgYEfmfPnymU2bNjk9b8CAAeaBBx4wxhiza9cuY7PZzN9//22MMSY5OdmUKVPGvP3229f9bK63741J+3dojjyS4uXlRVBQEBEREU7tERERNGrUyKJUclN+frB4MZW+nMaa0n2ZT2/HrLWnT9sYMABatoTff7/J64hIjleqFJQpk71LqVJpz2f+/1RLWo7O165d2/H3ggULUrhwYY4fPw7YL+BYt24dhQoVcizVqlUD7Ec3Lrvrrruc3iskJITff/+d5GT7eL5HHnmEjz/+mIsXL3Lp0iUWLlxI//79AShevDh9+/YlLCyMDh06MH36dGJiYpwyVqhQwWlesNKlSzsy7t27l4sXL9K6dWunnAsWLHBkrFevHtWqVePjjz8GIDIykuPHj6c6IpTZPLL01W/B2bNn+eOPPxyPo6OjiYqKonjx4pQrV46RI0fSq1cv6tevT0hICHPmzOHgwYMMGjTIwtSSJvfei61pU3qPGUP4nOo8wVQW8hAA69ZBrVqG55+3MXo0eHpanFVEssSOHVYnuLE77rgDm83Gvn376NSp0w37el71H5XNZiMlJQWwj+fo0KEDkyZNSvW80qVLpzlPhw4d8Pb2ZunSpXh7e5OQkECXLl0c6+fOncuwYcNYtWoVS5Ys4dlnnyUiIoK77rorTRkBvv76a8qUKePUz9vb2/H3nj17smjRIsaOHcuiRYsICwvD19c3zduQITc8zmKhdevWGSDV0qdPH0eft99+25QvX954eXmZO++800RGRmZbPp3uySTr1hlz++1mFaGmAn85HZqtXduYrVutDigit+JGh/xdXZs2bUyZMmXM2bNnU607ffq0MSb1cARjjPHx8TFz5841xhjz9NNPm6pVq5pLly5d930un+650tixY1O1jRkzxrRu3dq0b9/eDBw48IbZ77rrLvP4448bY+yne+rUqeO0furUqaZ8+fLGGGPi4uKMt7e3WbBgwQ1f86+//jKA2bFjhylatKj5+OOPb9g/V5/uadasGcaYVMu8efMcfYYMGcLff/9NQkICO3fupEmTJtYFloxp1gz27CFsTF1+dqvDSN5wXK68Zw/cdZdhxAg4e9bSlCKSB82cOZPk5GQaNmzI559/zu+//86+fft48803CQkJSdNrPPbYY5w6dYoHHniAbdu28ddff7FmzRr69+/vOJUDcOjQIUaOHMn+/fv5+OOPeeuttxg+fLjTaz388MOsXbuWb775xnGqB+xnGsaNG8fmzZs5cOAAa9as4bfffqN69eppyli4cGFGjx7NE088wfz58/nzzz/ZvXs3b7/9NvPnz3f0q1ixIo0aNWLAgAEkJSXRsWPHNL3+rXDZIkXykPz5YdIkCm5bxxt1PmQrwdRlN2C/XHn6dKhZE1autDiniOQpFStWZNeuXTRv3pxRo0YRGBhI69at+e6775wmDr0Rf39/fvjhB5KTkwkLCyMwMJDhw4fj4+ODm9t/X8G9e/fmwoULNGzYkMcee4zHH3+cgQMHOr3WHXfcQaNGjahatSrBwcGO9gIFCvDrr7/SpUsXqlSpwsCBAxk6dCiPPvpomrf15Zdf5rnnnmPChAlUr16dsLAwvvrqKypWrOjUr2fPnvz444/cd9995M+fP82vn1E2Y3T7t4yIi4vDx8eH2NhYzZmSmS5dgtdf59IL/2Nq4hCe50Uu8t8/hB49YPp0uO02CzOKSJpdvHiR6Ohox33WJOOMMVSrVo1HH32UkSNHWh3npm6079P6HaojKeJaPD1h3Dg89+xkzD2b+ZlAWvKtY/XixVCtGsydq7sri0jecfz4caZMmcI///xDv379rI6TbVSkiGuqWhUiI6n89igiCnZmLn0pzr8AnD4N/ftD69ZwxQVgIiK5lp+fHxMnTmTOnDkUK1bM6jjZRkWKuC43NxgyBNveX+jb9gT7qM4DLHKs/u47++XKkybZzxKJiORWxhhOnDjBgw8+aHWUbKUiRVxfuXKwYgW3LZzGIt/hrCScctjvSXHxoo2xY6FBA9efd0FERNJHRYrkDDYbPPgg7N1L+IPF+YWajGCq43LlH3+E4GDDyJFw7pzFWUUkFV2jkfdkxj5XkSI5S8mSsHAhhVYsYWrZKWzhLmrzIwApKTamTrVfrrxqlcU5RQT4b6bT8+fPW5xEstvlfX71bLfp4bLT4ovcULt28MsvNBg3jh0z6/MGo3iR57lIfg4cgPBw+4GXadPsdY2IWMPd3Z2iRYs67hNToECBdN2tXnIeYwznz5/n+PHjFC1aFHd39wy/luZJySDNk+JCNm6Ehx/m999SeJR3WEcLx6rixWHKFOjd237GSESynzGGo0ePcubMGaujSDYqWrQopUqVumZRmtbvUBUpGaQixcVcvAgvv4yZOIl5Kb0YxRucprhjdatWMHs2VK5sYUaRPC45OZlLuhQvT/D09LzhERQVKVlMRYqLioqCAQM4tuswI5jGYh5wrMqfH154AUaOBA+d6BQRsYxmnJW8qW5d2LoVv8mj+Thff1bQjgAOAnDhAjz1FDRsCLt2WRtTRERuTkWK5D4eHvDkk7BnD+2anuMXajKM6dhIAWD3bvu8KqNH63JlERFXpiJFcq877oC1ayn8zhtML/IcmwmhFnsASEmBN96AWrVgzRqLc4qIyDWpSJHczc0NBg6EvXsJvrcUOwnifzyNNxcBiI6GsDDo21dHVUREXI2KFMkbypSBZcvwXLKQp297nz3UphnrHKvnz4cmTeDwYQszioiIExUpknfYbNCtG+zdS5XeIaylBe8xgMLEAfbBtA0bwrZtFucUERFARYrkRSVKwPz52FatYkC5b9lMCBX5C4CYGGja1LB4scUZRURERYrkYWFh8Msv1OxRm60E05gNgP3Oyg88AM89Zx9gKyIi1lCRInlboUKwaBElXxrGt7SiP+87Vr38sv3skAbUiohYQ0WKiM0G48fjteQj3vMeyhuMdMyp8vnnGlArImIVFSkil3Xrhm3jBkaWXswK2jsNqG3QQANqRUSym4oUkSs1aADbt9P2zmNOA2qPHtWAWhGR7KYiReRqZcrAhg3U7FKdbTRMNaB2/HgNqBURyQ4qUkSupWBB+OQTfJ8dzLe0YgDvOVa98ooG1IqIZIc8W6TEx8fToEED6tatS61atXj33XetjiSuxs0NXn4Zr4XzeNdrKFN4AjeSAfuA2saNNaBWRCQr2YwxxuoQVkhOTiYhIYECBQpw/vx5AgMD2b59OyVKlEjT8+Pi4vDx8SE2NpYiRYpkcVqx3ObN0KkTK48H0YPFxGPf56VKwbJlEBxsbTwRkZwkrd+hefZIiru7OwUKFADg4sWLJCcnk0frNUmLkBDYto22tf+5xoBa+Phji/OJiORCLlukbNiwgQ4dOuDv74/NZmPZsmWp+sycOZOKFSuSL18+goKC2LhxY7re48yZM9SpU4eyZcsyZswYfH19Mym95Erly8MPP1Dz3tvZRkOaEAlAQgI8+KAG1IqIZDaXLVLOnTtHnTp1mDFjxjXXL1myhBEjRvDMM8+we/duGjduTHh4OAcPHnT0CQoKIjAwMNVy5MgRAIoWLcqPP/5IdHQ0ixYt4tixY9fNk5CQQFxcnNMieVChQvDFF/iOGUAErTWgVkQkC+WIMSk2m42lS5fSqVMnR1twcDB33nkns2bNcrRVr16dTp06MWHChHS/x+DBg2nRogX333//Nde/8MILvPjii6naNSYlD5s3D/PIQKYlPcZoXicFdwDq1YMvv4SAAIvziYi4qFw9JiUxMZGdO3cSGhrq1B4aGsqmTZvS9BrHjh1zHA2Ji4tjw4YNVK1a9br9x40bR2xsrGM5dOhQxjdAcoe+fbGtW8sTvh85zVC7ezc0bAhbt1qcT0Qkh8uRRcrJkydJTk7Gz8/Pqd3Pz4+jR4+m6TUOHz5MkyZNqFOnDvfccw9Dhw6ldu3a1+3v7e1NkSJFnBYR7rkHtm0jvOahaw6oXbTI4nwiIjmYh9UBboXNZnN6bIxJ1XY9QUFBREVFZUEqyXMqVoRNm6jZowfbvmlIFz5nA01JSICePWHvXnjpJfu0KyIiknY58r9NX19f3N3dUx01OX78eKqjKyLZokgR+OorfEf0IoLWPMx/kwP+739w//0aUCsikl45skjx8vIiKCiIiIgIp/aIiAgaNWpkUSrJ89zdYepUvN6ZwRz3IUxlhGOG2i++sM9Qq6FMIiJp57Kne86ePcsff/zheBwdHU1UVBTFixenXLlyjBw5kl69elG/fn1CQkKYM2cOBw8eZNCgQRamFgEGDsR2xx2M6NKFqqf304PFxOHD7t32myx/+aVmqBURSQuXvQR5/fr1NG/ePFV7nz59mDdvHmCfzG3y5MnExMQQGBjI1KlTadKkSbbk07T4clO//w4dOrB3vxsd+Iq/qAyAtzd88IF9AjgRkbword+hLlukuDoVKZImZ85At26cjNhFVz4jkmaOVc88owG1IpI35ep5UkRyjKJFYeVKfId0Zw2hqQbUdu2qAbUiItejIkUkq3l4wNtv4zVjKnPcBjsNqF261D7VigbUioikpiJFJLs89hi2b1YywmceK2hPEWIBiIqyD6jVDLUiIs5UpIhkp9BQ2LKF8Mq/s4W7qMSfABw7phlqRUSupiJFJLtVqwZbt1K9WSm20ZCmrAdwzFD7zDOQkmJtRBERV6AiRcQKJUrA6tWUeKQLawjlEeY4Vr36qn1A7dmzFuYTEXEBKlJErOLlBe+8g9fUybxjG8w0hjsNqNUMtSKS16lIEbGSzQYjRmD7ajnDC8/la9qlGlC7ZYu1EUVErKIiRcQVtGsHmzbRpsJ+tnAXlbHfEuLYMWjWDBYutDaeiIgVVKSIuIrAQNi2jer3+LKVYJqxDrAPqH3oIQ2oFZG8R0WKiCspWRK+/ZYSfTqwmjANqBWRPE1Fioir8faGuXPxmvQK7zCI6QxLNUPtwYMWZxQRyQYqUkRckc0GY8ZgW7aUYQU/cBpQ++OP0LAhbN5scUYRkSymIkXElXXsCD/8QJuAvakG1DZvrgG1IpK7qUgRcXV16tgH1Ab7XHdArTEWZxQRyQIqUkRyglKlYP16SjzYhjWEMpB3HKtefRUef1yFiojkPipSRHKKfPngo4/wfOUFZv//gFob9muS334bBg/WJcoikruoSBHJSWw2eOYZbJ9+yrD877GA3o4rf955BwYOVKEiIrmHihSRnKhrV9i4kYdKruEjHsKdJADefx/694fkZIvziYhkAhUpIjlVUBCsW8cDt63lYx5wFCrz50OfPpCUZHE+EZFbpCJFJCerWRPWr+f+Ut/zCd3w4BJgvzT5oYdUqIhIzqYiRSSnq14d1q/nPv+tfE4XPEkEYMkS6NEDLl2yOJ+ISAapSBHJDapWhchI7i27m6V0xosEAD7/HLp1g8REi/OJiGSAihSR3OL22yEyknblfmY59+LNRQCWLbOPs01IsDaeiEh6qUgRyU0qVYLISMIq/MYK2pOf8wB89RV07gwXL1qcT0QkHVSkiOQ2FSrA+vW0qhTN17SjAOcA+OYb+62ALlywNp6ISFrl6SLFw8ODunXrUrduXR5++GGr44hknvLlYf16mt9+mG8IpyBnAVizBtq3h/PnLc4nIpIGNmPy7h0/fH19OXnyZIaeGxcXh4+PD7GxsRQpUiSTk4lkkn/+gRYt+OE3X8L5hnjsP6tNm8KKFVCokMX5RCRPSut3aJ4+kiKS65UpA+vXc3e1U6whlCLEAhAZCeHhEB9vcT4RkRtw2SJlw4YNdOjQAX9/f2w2G8uWLUvVZ+bMmVSsWJF8+fIRFBTExo0b0/UecXFxBAUFcc899xAZGZlJyUVcTOnSsH49d9WI5ztaUpTTAHz/PYSFQWysxflERK7DZYuUc+fOUadOHWbMmHHN9UuWLGHEiBE888wz7N69m8aNGxMeHs7BgwcdfYKCgggMDEy1HDlyBIC///6bnTt3Mnv2bHr37k1cXFy2bJtItvPzg3XrqF8rke9oSXH+BWDzZggNhTNnrI0nInItOWJMis1mY+nSpXTq1MnRFhwczJ133smsWbMcbdWrV6dTp05MmDAh3e8RHh7Oyy+/TP369a+5PiEhgYQrJpqIi4sjICBAY1IkZzl5Elq14scfDa34lpOUBOy3AVqzBooXtzifiOQJuXpMSmJiIjt37iQ0NNSpPTQ0lE2bNqXpNU6fPu0oOg4fPszevXupVKnSdftPmDABHx8fxxIQEJDxDRCxiq8vrF1LnTs9WEdzSnIcgJ07oWVLew0jIuIqcmSRcvLkSZKTk/Hz83Nq9/Pz4+jRo2l6jX379lG/fn3q1KlD+/btmT59OsVv8GvkuHHjiI2NdSyHDh26pW0QsUzx4vDttwTWz896muGH/d9MVBS0aAEnTlgbT0TkMg+rA9wKm83m9NgYk6rteho1asRPP/2U5vfy9vbG29s7XflEXFaxYhARQY02bYjc2pTmrCMGf376CZo3h+++sw9jERGxUoaKlH/++YcffviB48ePk5KS4rRu2LBhmRLsRnx9fXF3d0911OT48eOpjq6IyHUULQpr1lA1PJzITU1pwVoOE8Avv0CzZrB2rf3CIBERq6S7SJk7dy6DBg3Cy8uLEiVKOB25sNls2VKkeHl5ERQUREREBJ07d3a0R0RE0LFjxyx/f5Fco0gRWLWKO9q2JfJ7+xGVg5Tn11//K1TKlLE6pIjkVekuUp577jmee+45xo0bh5tb1g1pOXv2LH/88YfjcXR0NFFRURQvXpxy5coxcuRIevXqRf369QkJCWHOnDkcPHiQQYMGZVkmkVypcGH45hsqtW9PZKS9UPmbivz2m31m2nXrQOPERcQK6b4EuUSJEmzbto3KlStnVSYA1q9fT/PmzVO19+nTh3nz5gH2ydwmT55MTEwMgYGBTJ06lSZNmmRprss0Lb7kOufOwb33cnDt7zRnHX9h/zdesaL9iEqFCtbGE5HcI63foekuUsaMGUPx4sUZO3bsLYfMyVSkSK50/jx06sThiL20YC2/UwWAcuXsR1RucJW+iEiaZVmRkpycTPv27blw4QK1atXC09PTaf2UKVMyljiHUZEiudaFC3DffRxZ9SMtWMt+qgFQtqy9ULn9dovziUiOl9bv0HSPSXn11VdZvXo1VatWBUg1cFZEcrj8+WHpUvy7dmX9181oyXfspSaHD/83RqVKFatDikhekO4jKcWKFWPq1Kn07ds3iyLlDDqSIrleQgJ068bx5ZtpyXf8TC0ASpWyj1GpXt3ifCKSY2XZtPje3t7cfffdtxRORHIAb2/49FNu63wP62hOHaIAOHrUfnnyzz9bmk5E8oB0FynDhw/nrbfeyoosIuJqvLxgyRJ8uzZnLS24k50AHD9un5l2zx6L84lIrpbu0z2dO3dm7dq1lChRgpo1a6YaOPvFF19kakBXpdM9kqckJUGvXpxevIowVrOdhoDjNkDUq2dxPhHJUbJs4GzRokW57777bimciOQwHh7w4YcUc+9LxMLWtGEVWwjh1Cn7TQkjIqB+fatDikhuk64iJSkpiWbNmhEWFkapUqWyKpOIuCIPD5g/Hx/3/qxeEEZbVvID93DmDLRqBatXQ3Cw1SFFJDdJ15gUDw8PBg8eTEJCQlblERFX5u4OH3xAkf73s4o2NCESgNhYaN0aNm2yOJ+I5CrpHjgbHBzM7t27syKLiOQE7u7w7rsUeuRBVtKW5qwFID4ewsJg40aL84lIrpHuMSlDhgxh1KhRHD58mKCgIAoWLOi0vnbt2pkWTkRclJsbzJ5NQQ8PVsxqTyeWEUEoZ89Cmzbw9df2y5RFRG5Fuq/uudadj202G8YYbDYbycnJmRbOlenqHhHAGBg+nAtvvct9fMEqwgH7pLVffQUtW1qcT0RcUpZd3RMdHX1LwUQkF7HZYPp08ru7s3RaZ7ryGV/TngsXoH17+PJLCA21OqSI5FTpLlLKly+fFTlEJKey2WDKFPJ5ePD5613ozhK+pBMXL8K998IXX0DbtlaHFJGcKN0DZwH+/PNPHn/8cVq1akXr1q0ZNmwYf/75Z2ZnE5GcwmaDyZPxHjuST+jGfXwO2G//07mz/dSPiEh6pbtIWb16NTVq1GDbtm3Url2bwMBAtm7dSs2aNYmIiMiKjCKSE9hs8OqreD37FIvpQTeWAJCYCF26wNKlFucTkRwn3QNn69WrR1hYGBMnTnRqHzt2LGvWrGHXrl2ZGtBVaeCsyA28+CJJL7xMH+aziJ6AfS64jz+Grl0tziYilsuyuyDv27ePAQMGpGrv378/e/fuTe/LiUhu9PzzeLz8AgvoTS8WAPbb//ToAYsXW5xNRHKMdBcpJUuWJCoqKlV7VFQUt912W2ZkEpHc4NlncZ/wP+bSj358AEByMvTsafjoI4uziUiOkO6rex555BEGDhzIX3/9RaNGjbDZbHz//fdMmjSJUaNGZUVGEcmpxo7F3cOD9558GE8uMYdHSUmx0bu3ITnZRp8+VgcUEVeW7iJl/PjxFC5cmDfeeINx48YB4O/vzwsvvMCwYcMyPaCI5HCjR+Pm4cGsJwbjQRIzeQxjbPTrZ0hKsnGNs8ciIkAGBs5eKT4+HoDChQtnWqCcQgNnRdJpxgzM448zgmm8yXAAbDbDokU2evSwOJuIZKssm3H2SnmxOBGRDBo6FJu7O9OGDMGDJKYwCmPsp358fW20amV1QBFxNekeOHvs2DF69eqFv78/Hh4euLu7Oy0iItc1eDC2OXN4ndE8zLsAXLpko3NnyCOzF4hIOqT7SErfvn05ePAg48ePp3Tp0thstqzIJSK51SOPYHN3Z9aAgRznNpbTkbNnITwcfvgBbr/d6oAi4irSXaR8//33bNy4kbp162ZBHBHJE/r3x+P0aRaP7kEoa/iexhw/DmFh9kKlVCmrA4qIK0j36Z6AgABuYayty9i/fz9169Z1LPnz52fZsmVWxxLJO0aNIv/ooSznXmryMwB//WW/GWFcnMXZRMQlpPvqnjVr1vDGG2/wzjvvUKFChSyKlb3Onj1LhQoVOHDgAAULFkzTc3R1j0gmSEmBPn3456O1NGITB7HfZb1FC1i5Ery9Lc4nIlkiy6bF7969O+vXr6dy5coULlyY4sWLOy050fLly2nZsmWaCxQRySRubvDBB5RpU5vVhFGcfwFYuxZ69bLPUCsieVe6x6RMmzYtC2KktmHDBl577TV27txJTEwMS5cupVOnTk59Zs6cyWuvvUZMTAw1a9Zk2rRpNG7cON3v9cknn9C7d+9MSi4i6eLpCZ9+SrUWLfh6eztasJYLFODTT8HPD958036DZRHJe9JdpPTJpnmsz507R506dejXrx9dunRJtX7JkiWMGDGCmTNncvfdd/POO+8QHh7O3r17KVeuHABBQUEkJCSkeu6aNWvw9/cH7IecfvjhBxbrrmci1ilUCL7+mrvuuYfPfuvKvSwnGQ9mzIDSpeHpp60OKCJWuKUZZ7OLzWZLdSQlODiYO++8k1mzZjnaqlevTqdOnZgwYUKaX/vDDz9k9erVfHSTO54lJCQ4FTxxcXEEBARoTIpIZvr7b2jUiPkxrenLfEfzu+/Cww9bF0tEMleWjUlxBYmJiezcuZPQ0FCn9tDQUDZt2pSu1/rkk0/o3r37TftNmDABHx8fxxIQEJCu9xGRNKhQAVatok+RZUxijKP50UcNy5dbF0tErJEji5STJ0+SnJyMn5+fU7ufnx9Hjx5N8+vExsaybds2wsLCbtp33LhxxMbGOpZDhw6lO7eIpEHt2rB8OU96vckTTAEgJcVG9+72OVREJO/IkUXKZVfPdmuMSdcMuD4+Phw7dgwvL6+b9vX29qZIkSJOi4hkkaZNsX28iNd5kgdZCMDFi9C+Pfzyi8XZRCTbpKtISUpKwsPDg59//jmr8qSJr68v7u7uqY6aHD9+PNXRFRHJoe67D7dZbzOXfoSyGoAzZ+yz0h48aG00Ecke6SpSPDw8KF++PMkWT17g5eVFUFAQERERTu0RERE0atTIolQikukGDcLr+af5jK7UZzsA//xjL1T+/dfibCKS5dJ9uufZZ59l3LhxnDp1KivyOJw9e5aoqCiioqIAiI6OJioqioP//yvUyJEjee+99/jggw/Yt28fTzzxBAcPHmTQoEFZmktEstnzz1N44IN8TTvu4DcAfv3Vfurn3DmLs4lIlkr3Jcj16tXjjz/+4NKlS5QvXz7VLK27Mul+6+vXr6d58+ap2vv06cO8efMA+2RukydPJiYmhsDAQKZOnUqTJk0y5f1vRtPii2Sj5GS4/36il+6mEZs4SmkA2rWDpUvt88GJSM6R1u/QdBcpL7744g3XP//88+l5uRxLRYpINrt4EUJD+XFjLE3YQBw+APTpA3PnalZakZwky4oUsVORImKBM2egSRPW/1ScMFaTiP0OhE89BRMnWhtNRNIuyydz27lzJx999BELFy5k9+7dGX0ZEZG0K1oUVq2iWbloFvEgNlIAmDQJsum2YiKSjdJ9757jx4/To0cP1q9fT9GiRTHGEBsbS/PmzVm8eDElS5bMipwiInb+/rB6NV3uuYeZ/w5hMLMBeOIJ+w0JH3jA4nwikmnSfSTl8ccfJy4ujl9++YVTp05x+vRpfv75Z+Li4hg2bFhWZBQRcVatGnz9NYMKfMhz/DdOrk8fw5o1FuYSkUyV7jEpPj4+fPvttzRo0MCpfdu2bYSGhnLmzJnMzOeyNCZFxAV88w2mfQcGpbzNHB4FoGBBWLcOrvovSkRcSJaNSUlJScHzGtf7eXp6kpKSkt6XExHJuPBwbHM/YCZD6MwXgH3ulLZt4fffLc4mIrcs3UVKixYtGD58OEeOHHG0/fPPPzzxxBO0bNkyU8OJiNxU7964T57IIh6kMRsAOHnSPittTIzF2UTklqS7SJkxYwbx8fFUqFCBypUrc/vtt1OxYkXi4+N56623siKjiMiNjR5NvieGsJx7qcUeAKKjITwcYmMtziYiGZbheVIiIiL49ddfMcZQo0YNWrVqldnZXJrGpIi4mJQU6NWLI4vW0YhNHKACAM2awTffQL58lqYTkStkyWRuSUlJ5MuXj6ioKAIDAzMlaE6lIkXEBSUmQocO7F/zN3fzA//iC0CXLrBkCbi7W5xPRIAsGjjrKndBFhG5Ji8v+OwzqtYvwte0owD2OxB+/jkMGwaaX1skZ3HZuyCLiGRI4cLw9dcE336Kz+mCB5cAmDkTXnnF4mwiki4uexdkV6fTPSIuLjoaGjXiw6Ot6M2HjuZ33oGBAy3MJSJp/g5N97T4nTp1upVcIiLZo2JF+OYbejVpwvH4UYzmDQAGDzbcdpsN/Vcm4vrSVaQkJSUB0L9/fwICArIkkIhIpqlbF778klFt2hCTWJo3GE1Kio0ePWDNGmjSxOqAInIj6R44+/rrr2vgrIjkHM2bw8KFTOYpHvr/0z4JCXDvvfDTTxZnE5EbSvfA2ZYtW7J+/fosiCIikkW6dsXt7bf4gP604RvAPslbmzZw4IDF2UTkutI9JiU8PJxx48bx888/ExQUlGrg7L333ptp4UREMs2QIXjGxPDpK/fTgrVspyFHjtinz//+e/D1tTqgiFwt3Vf3uLld/+CLzWbLM6eCdHWPSA5kDAwcyIn3lnEP3/MbVQFo2BDWrrXfQVlEsl6W3gX5ekteKVBEJIey2WDWLEre24jVhFEa+41St22Drl3h0iWL84mIk3QXKSIiOZqHByxeTIV7AlhFG3w4A8CqVTBggP0WQCLiGtJcpLRt25bYK24n+r///Y8zZ844Hv/777/UqFEjU8OJiGSJ/Plh+XJq10xhOffizUUAPvwQxo61OJuIOKS5SFm9ejUJCQmOx5MmTXKaGj8pKYn9+/dnbjoRkaxSrBisWkWTgL9ZxIO4YT9d/dpr8MYbFmcTESAdRcrV42vTOd5WRMT1lC0Lq1dzX/FIZjLE0Tx6NHz0kYW5RATQmBQRyeuqV4cVK3g0/4e8yHOO5n79DKtWWZhLRNJepNhsNmw2W6o2EZEcLyQEPv2U8W6vMohZACQl2eja1X7lj4hYI82TuRlj6Nu3L97e3gBcvHiRQYMGOSZzu3K8Sk7x+uuvM3fuXGw2G2PHjuWhhx6yOpKIWKVdO2zvv8eMfgM4zm18QRfOnYN27eyTvVWtanVAkbwnzUVKnz59nB5f6wu9d+/et54om/z0008sWrSInTt3Avbp/tu3b0/RokWtDSYi1unbF/ejR1k4ridtKEEkzTh50j4r7aZN4O9vdUCRvCXNRcrcuXOzMke227dvH40aNSJfvnwA1K1bl1WrVtGjRw+Lk4mIpZ56inxHj/Ll9I40JZIfqcuBA/b7/GzYAPo9RiT7uOzA2Q0bNtChQwf8/f2x2WwsW7YsVZ+ZM2dSsWJF8uXLR1BQEBs3bkzz6wcGBrJu3TrOnDnDmTNnWLt2Lf/8808mboGI5Eg2G0yZgk+PtnxDOBWIBux3TO7YES5etDifSB7iskXKuXPnqFOnDjNmzLjm+iVLljBixAieeeYZdu/eTePGjQkPD+fgwYOOPkFBQQQGBqZajhw5Qo0aNRg2bBgtWrSgc+fONGjQAA+PdN9vUURyIzc3mDeP0q0CWU0YvpwA7EdSHnwQdAcQkeyR7hsMWsFms7F06VI6derkaAsODubOO+9k1qxZjrbq1avTqVMnJkyYkO73ePjhh+ncuTPt2rW75vqEhASnwcFxcXEEBAToBoMiuVl8PDRrxvZdbjRnHecoBMDAgTB7tv2gi4ikX5bdYNAVJCYmsnPnTkJDQ53aQ0ND2bRpU5pf5/jx4wDs37+fbdu2ERYWdt2+EyZMwMfHx7EEBARkLLyI5ByFC8PKlTSofJovuA8P7HcgnDMHXnzR4mwieUCOLFJOnjxJcnIyfn5+Tu1+fn4cPXo0za/TqVMnatSowUMPPcTcuXNveLpn3LhxxMbGOpZDhw5lOL+I5CB+frB6NaG3/cg8+jqaX3xRs9KKZLUcPQjj6snkjDHpmmAuPUddvL29HXPEiEgeU7kyrFpFz6ZNORY/klFMAeDhhw133GEjONjifCK5VI48kuLr64u7u3uqoybHjx9PdXRFRCRT1KsHy5bxhMcMBvIOAAkJNjp1gsOHrY0mklvlyCLFy8uLoKAgIiIinNojIiJo1KiRRalEJNdr0QLbRx/yFo/TlPUAHD1qvzT5/Hlro4nkRi57uufs2bP88ccfjsfR0dFERUVRvHhxypUrx8iRI+nVqxf169cnJCSEOXPmcPDgQQYNGmRhahHJ9bp3x+vXX/nsha40ZBvRVGLXLujbF5Ys0RU/IpnJZS9BXr9+Pc2bN0/V3qdPH+bNmwfYJ3ObPHkyMTExBAYGMnXqVJo0aZIt+dJ6+ZSI5EIpKdCjBz9/upcQNnOWwoB9MO1zz93kuSKS5u9Qly1SXJ2KFJE87vx5uOcevtpdho58ifn/s+effgpdu1qcTcTF5ep5UkRELFegAHz5JR38tjORsY7m3r0Nu3dbmEskF1GRIiKSUQEBsGwZT3q9SS8WAHDhgo1777UPqBWRW6MiRUTkVtx1F7b33mUOA7mLzYD9kuTOnXUzQpFbpSJFRORW9epFvjHDWUpnymKfjXrLFnj0UdCoP5GMU5EiIpIZXn2VUu0bsJx7KcA5ABYsgNdftziXSA6mIkVEJDO4u8PChdSreYn59HE0P/WUYcUKC3OJ5GAqUkREMkuRIrB8OV1LrOcFngfAGBsPPAA//2xtNJGcSEWKiEhmqlQJPvuM8e4TuJ9PADh7Fu69F06etDibSA6jIkVEJLM1a4bbzBnMoy93shOA6Gj7JG+JiRZnE8lBVKSIiGSFgQMp8PjDfElHShEDQGQkDB2qK35E0kpFiohIVpkyhbKtqrOMTnhjnzTl3XdhxgyLc4nkECpSRESyiocHfPIJwXec5j0edjQ/8YQhIsLCXCI5hIoUEZGsVKwYLF/OQz4rGMsEAJKTbXTrBr/9ZnE2ERenIkVEJKtVqwZLlvA/23g6sByAM2egQwc4fdraaCKuTEWKiEh2CAvDbcrrLKQngfwE2I+k9OgBSUkWZxNxUSpSRESyy/DhFB7QneXciy8nAFizBkaPtjiXiItSkSIikl1sNpg5k4qNA/icLnhwCYDp0+1X/YiIMxUpIiLZycsLPv+cJuUPMovBjuYhQwyRkRbmEnFBKlJERLJbyZKwfDkPF1zMcKYBkJRko0sX+8y0ImKnIkVExAq1a8PChbzOk4SyGoB//7Xf4yc+3uJsIi5CRYqIiFU6dsTj1ZdYQneqsB+w3y25Z09ITrY4m4gLUJEiImKlsWMp+mA7vqIDRbFPmvLVV/DssxbnEnEBKlJERKxks8F771GlQVE+oRvu2CdNmTgRPvrI4mwiFlORIiJitfz5YdkyWvvvZRojHM0PP2zYssW6WCJWU5EiIuIK/P3hyy95zPt9HmU2AAkJNjp1gkOHrI0mYhUVKSIirqJ+fWzz5/EWj9OU9QAcOwYdO8K5c9ZGE7GCihQREVfSvTuez47lM7pSkb8A2L0b+vaFlBRro4lktzxRpHTu3JlixYrRtWvXdK0TEbHEiy/i27kJX9GBQtgnTfnsM3j5ZYtziWSzPFGkDBs2jAULFqR7nYiIJdzcYMECatbx5GMewIb9EMoLL9iLFZG8Ik8UKc2bN6dw4cLpXiciYplCheDLL2lfchuTeMrR3Lu3/fSPSF5geZGyYcMGOnTogL+/PzabjWXLlqXqM3PmTCpWrEi+fPkICgpi48aN2R9URCS7lS8PS5cy2mM6vZkPwIUL9qnzjx61OJtINrC8SDl37hx16tRhxowZ11y/ZMkSRowYwTPPPMPu3btp3Lgx4eHhHDx40NEnKCiIwMDAVMuRI0cyLWdCQgJxcXFOi4hIlrv7bmzvzuEdHiWETQAcPgydO8PFixZnE8liHlYHCA8PJzw8/Lrrp0yZwoABA3j44YcBmDZtGqtXr2bWrFlMmDABgJ07d2Z5zgkTJvDiiy9m+fuIiKTSty/5fvqJL6bcR0O2cYhybNkCAwfC/Pn2SWtFciPLj6TcSGJiIjt37iQ0NNSpPTQ0lE2bNmVrlnHjxhEbG+tYDml2JRHJTpMnU6pNPb6kIwWwT5ry4Yfw2msW5xLJQi5dpJw8eZLk5GT8/Pyc2v38/DiajhOyYWFh3H///axcuZKyZcuyffv2NK27kre3N0WKFHFaRESyjbs7LF5MvWoXWUBvR/PYsYavvrIwl0gWsvx0T1rYrjqWaYxJ1XYjq1evztA6ERGX4uMDX31Fl4YNefH0czzPSxhj48EHYfNmCAy0OqBI5nLpIym+vr64u7unOmpy/PjxVEdXRETyhNtvh88+Y7zbq3RjCQBnz9qv+Dl50uJsIpnMpYsULy8vgoKCiIiIcGqPiIigUaNGFqUSEbFYixbY3nqTufTjTuwXDkRHQ9eukJhocTaRTGR5kXL27FmioqKIiooCIDo6mqioKMclxiNHjuS9997jgw8+YN++fTzxxBMcPHiQQYMGWZhaRMRiQ4ZQYHBfvqQjpYgBIDIShg4FYyzOJpJJbMZY++O8fv16mjdvnqq9T58+zJs3D7BP5jZ58mRiYmIIDAxk6tSpNGnSJJuTOouLi8PHx4fY2FgNohURa1y6BGFhbF13jqZEkkA+AN58Ex5/3OJsIjeQ1u9Qy4uUnEpFioi4hH//heBgFv4ZzEMsBMDNzfDNNzaumr1BxGWk9TvU8tM9IiJyC0qUgOXL6Vn4K8bxKgApKTa6dYP9+y3OJnKLVKSIiOR0NWrA4sW8wng6sgyA2Fj7FT+nT1sbTeRWqEgREckN2rbF7bVJfEgvarEHgN9+g+7dISnJ4mwiGaQiRUQktxg1isJ9urCce/HlBAARETB6tMW5RDJIRYqISG5hs8E771AhxJ/P6YIn9klTpk+Hd9+1OJtIBqhIERHJTby9YelSmgT8zSwGO5qHDDFERlqYSyQDVKSIiOQ2fn6wfDkDCixmBFMBSEqy0aWLfWZakZxCRYqISG5Uty4sWMBrPEkYqwD7lCodOkBcnLXRRNJKRYqISG7VpQseLz3PYnpQlV8B+OUX6NIFYmIsziaSBipSRERys2efpWj3NnxFB4pxCoBvv4WqVWHaNF2eLK5NRYqISG5ms8EHH3BHkA9L6Uxx/gUgPh6eeALuvBM2brQ4o8h1qEgREcntChSAZctoWuo3fqMKA3kHGykA/PQTNGkCffrAsWMW5xS5iooUEZG8oGxZWLGCEqW9eYdBbOEugtjhWL1ggf0U0IwZOgUkrkNFiohIXhEUBHv2QMeONGQ7WwlmFoMo5h4L2O/38/jj0KABbNpkcVYRVKSIiOQtvr6wdCnMmoV7fm8G8Q77k2+nv9s8R5eoKLj7bujfH06csCypiIoUEZE8x2aDQYNgxw6oU4eSnOT9lH5sIoS6Rf50dJs7F6pUgVmzIDnZwrySZ6lIERHJq2rUgK1b7Zf5ACFsYUdcFd4qNA6fgpcAOHMGhgyB4GDYts3CrJInqUgREcnLvL1hyhT45hvw88OdFIaencj+c2XpU/2/qmTnTrjrLhg40D5zrUh2UJEiIiLQpo19UG3btgD4cZx5+4LZULkvtapcBMAY+92Uq1Sx/5mSYmVgyQtUpIiIiN1tt8GKFfDmm/YjLEDjP+ez64Av07p+T+HCBoBTp+xHVEJC7EdYRLKKihQREfmPzWa/DnnbNvuYFcAj4RzDP2vM/pB+9Ox60dF12zb75cpDhtgLF5HMpiJFRERSq13bfvXPY485mkqvmc9HP1Ri/es7qFnT3maM/eqfqlXhgw90Ckgyl4oUERG5tvz57VPQLl9un18FICaGpqMbsDtsLK9PTKJQIXvzyZMwYADcc499nhWRzKAiRUREbqxDB/ug2tatHU2eUyYx6pNgfv36T7p3/6/r5s32iW2HDbNfvixyK1SkiIjIzZUuDatWwRtvgKenvW3XLsqE12Zx6/f5NsJQrZq9OSUF3nrLfgpowQL7KSGRjFCRIiIiaePmBiNH2ieAq1rV3nb+PDz8MC3f6caP608zcaL9pssAx4/b767cpIn9bssi6ZUnipTOnTtTrFgxunbt6tQeHx9PgwYNqFu3LrVq1eLdd9+1KKGISA5Sr5792uOBA/9r++wzvOrX5qm7Ivn1V7jyv9vvv7c/5YknIC4u++NKzpUnipRhw4axYMGCVO0FChQgMjKSqKgotm7dyoQJE/hXUymKiNxcwYLwzjvwxRdQvLi97fBhaN6cgNnP8OmiS6xaBXfcYV+VnAzTptkPwCxapFNAkjZ5okhp3rw5hQsXTtXu7u5Ogf8/Lnnx4kWSk5Mx+pcjIpJ2nTvDjz9C8+b2x8bAq6/CPfcQdvuf/PQT/O9/9guFAI4ehZ497d1/+cW62JIzWF6kbNiwgQ4dOuDv74/NZmPZsmWp+sycOZOKFSuSL18+goKC2LhxY6a9/5kzZ6hTpw5ly5ZlzJgx+F6+zE5ERNKmbFmIiICJE8HDw962bRvUrYv3kgU8Pc6wdy906vTfUyIjoW5dePJJiI+3IrTkBJYXKefOnaNOnTrMmDHjmuuXLFnCiBEjeOaZZ9i9ezeNGzcmPDycgwcPOvoEBQURGBiYajly5MhN379o0aL8+OOPREdHs2jRIo4dO3bNfgkJCcTFxTktIiLy/9zd4amnYNMmuP12e9vZs/aRsz17UqFYLEuXwtdfQ+XK9tVJSfD661CtGixZolNAcg3GhQBm6dKlTm0NGzY0gwYNcmqrVq2aGTt2bLpee926daZLly437DNo0CDzySefXHPd888/b4BUS2xsbLpyiIjkevHxxvTrZ4y97rAv5csb8/33xhhjLlww5sUXjcmXz7lLy5bG7NtnbXTJHrGxsWn6DrX8SMqNJCYmsnPnTkJDQ53aQ0ND2bRp0y2//rFjxxxHROLi4tiwYQNVL19Wd5Vx48YRGxvrWA4dOnTL7y8ikisVKmSfI3/xYvDxsbcdOGC/FvmFF8jnkcRzz9nHpLRr99/TvvvOPhv/2LFw7pw10cW1uHSRcvLkSZKTk/Hz83Nq9/Pz4+jRo2l+nbCwMO6//35WrlxJ2bJl2b59OwCHDx+mSZMm1KlTh3vuuYehQ4dSu3bta76Gt7c3RYoUcVpEROQGune3D6q95x7745QUePFFaNoU/v6bSpXsN11evhwqVLB3uXQJJk2C6tXh8891Ciivc+ki5TKbzeb02BiTqu1GVq9ezYkTJzh//jyHDx+mQYMGgH0sS1RUFD/++CN79uxh8ODBmZpbRCTPK18e1q2Dl16yj1sB+7iVOnXg448B+6z7v/wC48eDl5e9y6FD9rlW2rSB33+3KLtYzqWLFF9fX9zd3VMdNTl+/HiqoysiIuKiPDzsFcjGjf8dMomLgwcftA+sjY+nQAF7HfPLL/bC5LI1ayAwEJ591j65rWS95GT7fZcOHLDfsmnjRli50posHta8bdp4eXkRFBREREQEnTt3drRHRETQsWNHC5OJiEi6hYTYb5H82GOwcKG9bcEC+5S0ixZBcDC3327/Qly2DEaMgIMHITHRPtfKRx/Z51gpVMg+l9zlP6+3FCoE+fJBOg6853hJSfb6LzY2Y0tc3PUvCb906b8rzLOL5UXK2bNn+eOPPxyPo6OjiYqKonjx4pQrV46RI0fSq1cv6tevT0hICHPmzOHgwYMMGjTIwtQiIpIhPj72aqNNGxgyxP6N+NdfcPfd9vEqY8dic3enc2cIDbXPC/faa/YvyAMH7I/Tw2a7cRGTkXVZVQAlJjoXCxkpMrLyaFNc3H+TC2cXmzHWDktav349zS/PVHiFPn36MG/ePMA+mdvkyZOJiYkhMDCQqVOn0qRJk2xO6iwuLg4fHx9iY2M1iFZEJCP++st+aGTLlv/amjSxFzEBAY6m/fth6FD49lsLMt5AegqgAgXgwoUbFxgXL1q3LQUK2OvHGy3Dh0Nmfd2l9TvU8iIlp1KRIiKSCS5dgpdftp/PSUmxtxUtCu++63SXQmNg716IibFfnnz2rP3Pay03WnfunL1YyE0KFbp5gXGtpUiR//709MzezCpSspiKFBGRTLRxo/2oypVzUA0YYL8rYaFCmfpWycn20yIZKXBuVhyl52iIzWYvEC4XCxlZChf+76KpnERFShZTkSIikslOn4ZBg+CTT/5rq1LFPqg2KMi6XOlwowIof37nAqNQIXBz6Wtss46KlCymIkVEJAsYA/PmweOP/zftrKen/XTQqFF591s9l0nrd6j2toiIuA6bDfr1g927oX59e9ulSzBmjP1yn82b7QNu4+I0HW0eoCMpGaQjKSIiWSwxEZ5/3j5P/rW+qjw9oUQJ8PX9788r/36ttiJF8tbEKS5Kp3uymIoUEZFssnYt9OoFR47c+mt5eNgLlhsVNVev8/HRaaZMltbvUMsncxMREbmhFi3s87PPnQvR0XDyJPz7r/Ofab2sJikJjh2zL2nl7m6fxSwtR2ou/1m0aPYUNsbYtykpyT5q9/Lf07Kkp39yMvTvn+1TzupISgbpSIqIiAs5f95erFyrgLn859VtWTk9q5vbf4XN5SM3bm6ZX0hcnlsmO5w7Z5/1LRPoSIqIiOQdBQpAuXL2Ja0uXLh+AXO9IufyFUc3k5Ly33Nyi6SkbH9LFSkiIpI35c8PZcval7S6eNFetNzsKM2Vf17vjn1XcnOzn0q5cnF3T912vSU7+np7Z/yzziAVKSIiImmVLx+UKWNf0ioxEU6dsv/9WgWCu7sG5l6HihQREZGs5OUFpUpZnSJHUukmIiIiLklFioiIiLgkFSkiIiLiklSkiIiIiEtSkSIiIiIuSUWKiIiIuCQVKSIiIuKSVKSIiIiIS1KRIiIiIi5JRYqIiIi4JE2Ln0HGGMB+u2kRERFJu8vfnZe/S69HRUoGxf//XS0DAgIsTiIiIpIzxcfH4+Pjc931NnOzMkauKSUlhSNHjlC4cGFsNpvVcTJdXFwcAQEBHDp0iCJFilgdJ0tpW3OfvLKdoG3NjfLCdhpjiI+Px9/fH7cb3AFaR1IyyM3NjbJly1odI8sVKVIk1/4juZq2NffJK9sJ2tbcKLdv542OoFymgbMiIiLiklSkiIiIiEtSkSLX5O3tzfPPP4+3t7fVUbKctjX3ySvbCdrW3CivbGdaaOCsiIiIuCQdSRERERGXpCJFREREXJKKFBEREXFJKlJERETEJalIyeM2bNhAhw4d8Pf3x2azsWzZMqf1ffv2xWazOS133XWXNWFvwYQJE2jQoAGFCxfmtttuo1OnTuzfv9+pjzGGF154AX9/f/Lnz0+zZs345ZdfLEqccWnZ1tywX2fNmkXt2rUdE16FhITwzTffONbnlv0JN9/W3LA/r2fChAnYbDZGjBjhaMtN+/aya21nbt6vaaUiJY87d+4cderUYcaMGdft06ZNG2JiYhzLypUrszFh5oiMjOSxxx5jy5YtREREkJSURGhoKOfOnXP0mTx5MlOmTGHGjBls376dUqVK0bp1a8d9mnKKtGwr5Pz9WrZsWSZOnMiOHTvYsWMHLVq0oGPHjo4vq9yyP+Hm2wo5f39ey/bt25kzZw61a9d2as9N+xauv52QO/druhiR/weYpUuXOrX16dPHdOzY0ZI8Wen48eMGMJGRkcYYY1JSUkypUqXMxIkTHX0uXrxofHx8zOzZs62KmSmu3lZjcu9+LVasmHnvvfdy9f687PK2GpM792d8fLy54447TEREhGnatKkZPny4MSb3/Vu93nYakzv3a3rpSIrc1Pr167ntttuoUqUKjzzyCMePH7c60i2LjY0FoHjx4gBER0dz9OhRQkNDHX28vb1p2rQpmzZtsiRjZrl6Wy/LTfs1OTmZxYsXc+7cOUJCQnL1/rx6Wy/LTfsT4LHHHqNdu3a0atXKqT237dvrbedluW2/ppduMCg3FB4ezv3330/58uWJjo5m/PjxtGjRgp07d+bY2RCNMYwcOZJ77rmHwMBAAI4ePQqAn5+fU18/Pz8OHDiQ7Rkzy7W2FXLPfv3pp58ICQnh4sWLFCpUiKVLl1KjRg3Hl1Vu2p/X21bIPfvzssWLF7Nr1y62b9+eal1u+rd6o+2E3LdfM0JFitxQ9+7dHX8PDAykfv36lC9fnq+//pr77rvPwmQZN3ToUPbs2cP333+fap3NZnN6bIxJ1ZaTXG9bc8t+rVq1KlFRUZw5c4bPP/+cPn36EBkZ6Vifm/bn9ba1Ro0auWZ/Ahw6dIjhw4ezZs0a8uXLd91+OX3fpmU7c9N+zSid7pF0KV26NOXLl+f333+3OkqGPP744yxfvpx169ZRtmxZR3upUqWA/35Lu+z48eOpfmPLKa63rdeSU/erl5cXt99+O/Xr12fChAnUqVOH6dOn58r9eb1tvZacuj8Bdu7cyfHjxwkKCsLDwwMPDw8iIyN588038fDwcOy/nL5vb7adycnJqZ6Tk/drRqlIkXT5999/OXToEKVLl7Y6SroYYxg6dChffPEFa9eupWLFik7rK1asSKlSpYiIiHC0JSYmEhkZSaNGjbI77i252bZeS07dr1czxpCQkJCr9uf1XN7Wa8nJ+7Nly5b89NNPREVFOZb69evTs2dPoqKiqFSpUq7YtzfbTnd391TPycn7NcMsG7IrLiE+Pt7s3r3b7N692wBmypQpZvfu3ebAgQMmPj7ejBo1ymzatMlER0ebdevWmZCQEFOmTBkTFxdndfR0GTx4sPHx8THr1683MTExjuX8+fOOPhMnTjQ+Pj7miy++MD/99JN54IEHTOnSpXPdtuaW/Tpu3DizYcMGEx0dbfbs2WOefvpp4+bmZtasWWOMyT3705gbb2tu2Z83cvVVL7lp317pyu3MC/s1LVSk5HHr1q0zQKqlT58+5vz58yY0NNSULFnSeHp6mnLlypk+ffqYgwcPWh073a61jYCZO3euo09KSop5/vnnTalSpYy3t7dp0qSJ+emnn6wLnUE329bcsl/79+9vypcvb7y8vEzJkiVNy5YtHQWKMblnfxpz423NLfvzRq4uUnLTvr3SlduZF/ZrWtiMMSa7j96IiIiI3IzGpIiIiIhLUpEiIiIiLklFioiIiLgkFSkiIiLiklSkiIiIiEtSkSIiIiIuSUWKiIiIuCQVKSIiIuKSVKSIiIiIS1KRIiK5Qt++fbHZbEycONGpfdmyZdhsNotSicitUJEiIrlGvnz5mDRpEqdPn7Y6iohkAhUpIpJrtGrVilKlSjFhwgSro4hIJlCRIiK5hru7O6+++ipvvfUWhw8ftjqOiNwiFSkikqt07tyZunXr8vzzz1sdRURukYoUEcl1Jk2axPz589m7d6/VUUTkFqhIEZFcp0mTJoSFhfH0009bHUVEboGH1QFERLLCxIkTqVu3LlWqVLE6iohkkI6kiEiuVKtWLXr27Mlbb71ldRQRySAVKSKSa7388ssYY6yOISIZZDP6FywiIiIuSEdSRERExCWpSBERERGXpCJFREREXJKKFBEREXFJKlJERETEJalIEREREZekIkVERERckooUERERcUkqUkRERMQlqUgRERERl6QiRURERFzS/wHWRWvMCowT1wAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "%matplotlib inline\n", "\n", "import matplotlib.pyplot as plt\n", "\n", "N = range(12, 50, 4)\n", "error = {}\n", "for basis in ('legendre', 'chebyshev'):\n", " error[basis] = []\n", " for i in range(len(N)):\n", " errN = main(N[i], basis)\n", " error[basis].append(errN)\n", "\n", "plt.figure(figsize=(6, 4))\n", "for basis, col in zip(('legendre', 'chebyshev'), ('r', 'b')):\n", " plt.semilogy(N, error[basis], col, linewidth=2)\n", "plt.title('Convergence of Poisson solvers 1D')\n", "plt.xlabel('N')\n", "plt.ylabel('Error norm')\n", "plt.legend(('Legendre', 'Chebyshev'))\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "09bd3e41", "metadata": { "editable": true }, "source": [ "The spectral convergence is evident and we can see that\n", "after $N=40$ roundoff errors dominate as the errornorm trails off around $10^{-14}$." ] }, { "cell_type": "markdown", "id": "d778fc61", "metadata": { "editable": true }, "source": [ "## Complete solver\n", "\n", "\n", "A complete solver, that can use any family of bases (Chebyshev, Legendre, Jacobi, Chebyshev second kind),\n", "and any kind of boundary condition, can be found [here](https://github.com/spectralDNS/shenfun/blob/master/demo/poisson1D.py).\n", "\n", "\n", "\n", "1. **J. Shen**. Efficient Spectral-Galerkin Method I. Direct Solvers of Second- and Fourth-Order Equations Using Legendre Polynomials, *SIAM Journal on Scientific Computing*, 15(6), pp. 1489-1505, [doi: 10.1137/0915089](https://dx.doi.org/10.1137/0915089), 1994.\n", "\n", "2. **J. Shen**. Efficient Spectral-Galerkin Method II. Direct Solvers of Second- and Fourth-Order Equations Using Chebyshev Polynomials, *SIAM Journal on Scientific Computing*, 16(1), pp. 74-87, [doi: 10.1137/0916006](https://dx.doi.org/10.1137/0916006), 1995." ] } ], "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.12.3" } }, "nbformat": 4, "nbformat_minor": 5 }