{
"cells": [
{
"cell_type": "markdown",
"id": "03e210dc",
"metadata": {
"editable": true
},
"source": [
"\n",
"\n",
"# Demo - The Tau method\n",
"**Mikael Mortensen** (email: `mikaem@math.uio.no`), Department of Mathematics, University of Oslo.\n",
"\n",
"Date: **June 15, 2021**\n",
"\n",
"**Summary.** Shenfun has primarily been developed for the spectral Galerkin method.\n",
"However, there are other methods out there that make use of global basis\n",
"functions and variational principles. One such method, which has a lot in\n",
"common with the spectral Galerkin method, is the Tau method. The\n",
"principle difference between a Tau method and a spectral Galerkin method\n",
"is in the choice of basis functions. The spectral Galerkin method is\n",
"usually defined through function spaces where the boundary conditions\n",
"of the problem are already built in. The tau-method, on the other hand,\n",
"usually considers\n",
"only the orthogonal basis, like pure Chebyshev or Legendre polynomials,\n",
"and derives differentiation matrices for these bases. The boundary conditions\n",
"are then usually fixed through manipulation of a couple of rows of the\n",
"differentiation matrix. In this demo we will show how the original\n",
"tau-method can be easily implemented using [shenfun](https://github.com/spectralDNS/shenfun)."
]
},
{
"cell_type": "markdown",
"id": "4032d1d0",
"metadata": {
"editable": true
},
"source": [
"## The tau method for Poisson's equation in 1D\n",
"\n",
"Poisson's equation is given on a domain $\\Omega = (-1, 1)$ as"
]
},
{
"cell_type": "markdown",
"id": "6860396a",
"metadata": {
"editable": true
},
"source": [
"\n",
"\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\nabla^2 u(x) = f(x) \\quad \\text{for }\\, x \\in \\Omega, \\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": "b722d3fb",
"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. To solve Eq. ([1](#eq:poisson)) with the tau method we choose either\n",
"Legendre of Chebyshev basis functions $\\phi_k$, and then look for\n",
"solutions"
]
},
{
"cell_type": "markdown",
"id": "c9f32970",
"metadata": {
"editable": true
},
"source": [
"\n",
"\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"u(x) = \\sum_{k=0}^{N-1} \\hat{u}_k \\phi_k(x), \\label{eq:u} \\tag{3}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "f10d071a",
"metadata": {
"editable": true
},
"source": [
"where $N$ is the size of the discretized problem and\n",
"$\\hat{\\mathbf{u}} = \\{\\hat{u}_k\\}_{k=0}^{N-1}$ are the unknown expansion\n",
"coefficients. For this function to satisfy the given boundary conditions, it is necessary\n",
"that"
]
},
{
"cell_type": "markdown",
"id": "fe13714e",
"metadata": {
"editable": true
},
"source": [
"\n",
"\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"u(-1) = \\sum_{k=0}^{N-1} \\hat{u}_k \\phi_k(-1) = \\sum_{k=0}^{N-1}\\hat{u}_{k}(-1)^k = a,\n",
"\\label{eq:dirichleta} \\tag{4} \n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "e0a0c7d7",
"metadata": {
"editable": true
},
"source": [
"\n",
"\n",
"\n",
"$$\n",
"\\begin{equation} \n",
"u(+1) = \\sum_{k=0}^{N-1} \\hat{u}_k \\phi_k(+1) = \\sum_{k=0}^{N-1} \\hat{u}_{k} = b,\n",
"\\label{eq:dirichletb} \\tag{5}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "c7ad046f",
"metadata": {
"editable": true
},
"source": [
"where we have use that $\\phi_k(1) = 1$ and $\\phi_k(-1)=(-1)^k$ for $k=0,1, \\ldots, N-1$,\n",
"for either Chebyshev or Legendre polynomials $\\phi_k$. This gives two equations\n",
"for the $N$ unknowns in $\\{\\hat{u}_k\\}_{k=0}^{N-1}$.\n",
"We will now use variational principles, like in the Galerkin method, in order to\n",
"derive equations that can be used to close the remaining $N-2$ unknowns.\n",
"To this end we multiply Poisson's equation by a\n",
"test function $v$, a weight $w$, and integrate over the domain"
]
},
{
"cell_type": "markdown",
"id": "389e6324",
"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. \\label{eq:varform} \\tag{6}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "29046d76",
"metadata": {
"editable": true
},
"source": [
"The weight function depends on the choice of basis functions. For Chebyshev\n",
"it will be $1/\\sqrt{1-x^2}$, whereas it is unity for Legendre.\n",
"\n",
"Finally, we insert the trial function $u=\\sum_{j=0}^{N-1}\\hat{u}_j \\phi_j$ and\n",
"the test function $v=\\phi_k$, to get"
]
},
{
"cell_type": "markdown",
"id": "c390a448",
"metadata": {
"editable": true
},
"source": [
"\n",
"\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\int_{-1}^1 \\nabla^2 \\sum_{j=0}^{N-1} \\phi_j \\, \\, \\phi_k \\, w\\, dx \\hat{u}_j = \\int_{-1}^1 f \\, \\phi_k\\, w\\, dx. \\label{eq:varform2} \\tag{7}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "9d4a27c7",
"metadata": {
"editable": true
},
"source": [
"This problem can be reformulated as a linear algebra problem,"
]
},
{
"cell_type": "markdown",
"id": "1dd2760c",
"metadata": {
"editable": true
},
"source": [
"\n",
"\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"a_{kj} \\hat{u}_j = \\tilde{f}_k, \n",
"\\label{_auto2} \\tag{8}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "f79f7562",
"metadata": {
"editable": true
},
"source": [
"\n",
"\n",
"\n",
"$$\n",
"\\begin{equation} \n",
"A \\hat{\\mathbf{{u}}} = \\tilde{\\mathbf{{f}}}.\n",
"\\label{_auto3} \\tag{9}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "3949a751",
"metadata": {
"editable": true
},
"source": [
"However, the matrix $A\\in \\mathbb{R}^{N \\times N}$ is singular because it\n",
"contains two zero rows. These two rows are used to implement the two boundary\n",
"conditions. Setting $a_{N-2, j}=1$ and\n",
"$a_{N-1, j}= (-1)^j$ for $j=0, 1, \\ldots, N-1$, and also fixing\n",
"the right hand sides $\\tilde{f}_{N-2}=a$ and $\\tilde{f}_{N-1}=b$, the\n",
"two boundary conditions will be satisfied."
]
},
{
"cell_type": "markdown",
"id": "fc08f897",
"metadata": {
"editable": true
},
"source": [
"## Implementation"
]
},
{
"cell_type": "markdown",
"id": "6f3913f0",
"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), and the [scipy.sparse](https://docs.scipy.org/doc/scipy/reference/sparse.html)\n",
"for handeling sparse matrices:"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "d2f50904",
"metadata": {
"collapsed": false,
"editable": true,
"execution": {
"iopub.execute_input": "2024-02-19T13:48:49.483400Z",
"iopub.status.busy": "2024-02-19T13:48:49.483112Z",
"iopub.status.idle": "2024-02-19T13:48:50.330161Z",
"shell.execute_reply": "2024-02-19T13:48:50.329453Z"
}
},
"outputs": [],
"source": [
"from shenfun import inner, div, grad, TestFunction, TrialFunction, Function, \\\n",
" project, Dx, Array, FunctionSpace, dx\n",
"import numpy as np\n",
"import scipy.sparse as scp\n",
"from sympy import symbols, cos, sin, exp, lambdify"
]
},
{
"cell_type": "markdown",
"id": "40fcaf01",
"metadata": {
"editable": true
},
"source": [
"We use `Sympy` for a manufactured solution and `Numpy` for testing.\n",
"The exact manufactured solution $u_e(x)$ and the right hand side\n",
"$f_e(x)$ are created using `Sympy` as follows"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "b7f66a59",
"metadata": {
"collapsed": false,
"editable": true,
"execution": {
"iopub.execute_input": "2024-02-19T13:48:50.333410Z",
"iopub.status.busy": "2024-02-19T13:48:50.332992Z",
"iopub.status.idle": "2024-02-19T13:48:50.343012Z",
"shell.execute_reply": "2024-02-19T13:48:50.342398Z"
}
},
"outputs": [],
"source": [
"x = symbols(\"x\")\n",
"ue = sin(4*np.pi*x)\n",
"fe = ue.diff(x, 2)"
]
},
{
"cell_type": "markdown",
"id": "739884be",
"metadata": {
"editable": true
},
"source": [
"Note that we compute the right hand side function `fe` that corresponds to\n",
"the manufactured solution `ue`."
]
},
{
"cell_type": "markdown",
"id": "377d957f",
"metadata": {
"editable": true
},
"source": [
"### Discretization\n",
"\n",
"We create a basis with a given number of basis functions,"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "12458aa3",
"metadata": {
"collapsed": false,
"editable": true,
"execution": {
"iopub.execute_input": "2024-02-19T13:48:50.345766Z",
"iopub.status.busy": "2024-02-19T13:48:50.345545Z",
"iopub.status.idle": "2024-02-19T13:48:50.355966Z",
"shell.execute_reply": "2024-02-19T13:48:50.354823Z"
}
},
"outputs": [],
"source": [
"N = 32\n",
"T = FunctionSpace(N, 'Chebyshev')\n",
"#T = FunctionSpace(N, 'Legendre')"
]
},
{
"cell_type": "markdown",
"id": "736218b2",
"metadata": {
"editable": true
},
"source": [
"Note that we can either choose a Legendre or a Chebyshev basis. The\n",
"remaining code will work either way."
]
},
{
"cell_type": "markdown",
"id": "37f9d1b5",
"metadata": {
"editable": true
},
"source": [
"### Variational formulation\n",
"\n",
"The variational problem ([6](#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": "fa34ca27",
"metadata": {
"collapsed": false,
"editable": true,
"execution": {
"iopub.execute_input": "2024-02-19T13:48:50.359251Z",
"iopub.status.busy": "2024-02-19T13:48:50.359054Z",
"iopub.status.idle": "2024-02-19T13:48:50.370371Z",
"shell.execute_reply": "2024-02-19T13:48:50.369768Z"
}
},
"outputs": [],
"source": [
"u = TrialFunction(T)\n",
"v = TestFunction(T)\n",
"# Assemble differentiation matrix\n",
"A = inner(v, div(grad(u)))\n",
"# Assemble right hand side\n",
"fj = Array(T, buffer=fe)\n",
"f_hat = Function(T)\n",
"f_hat = inner(v, fj, output_array=f_hat)"
]
},
{
"cell_type": "markdown",
"id": "a10787e4",
"metadata": {
"editable": true
},
"source": [
"Note that the `sympy` function `fe` is be used to initialize the [Array](https://shenfun.readthedocs.io/en/latest/shenfun.forms.html#shenfun.forms.arguments.Array)\n",
"`fj`, 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. An\n",
"[Array](https://shenfun.readthedocs.io/en/latest/shenfun.forms.html#shenfun.forms.arguments.Array) contains the solution evaluated on the\n",
"quadrature mesh. A [Function](https://shenfun.readthedocs.io/en/latest/shenfun.forms.html#shenfun.forms.arguments.Function) represents a global\n",
"expansion, like Eq. ([3](#eq:u)), and its values are the\n",
"expansion coefficients $\\{\\hat{u}_{k}\\}_{k=0}^{N-1}$."
]
},
{
"cell_type": "markdown",
"id": "4a47d951",
"metadata": {
"editable": true
},
"source": [
"### Fix boundary conditions\n",
"\n",
"We fix two rows of the differentiation matrix in order to satisfy\n",
"Eqs. ([4](#eq:dirichleta)) and ([5](#eq:dirichletb))."
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "7da80e99",
"metadata": {
"collapsed": false,
"editable": true,
"execution": {
"iopub.execute_input": "2024-02-19T13:48:50.373092Z",
"iopub.status.busy": "2024-02-19T13:48:50.372887Z",
"iopub.status.idle": "2024-02-19T13:48:50.382392Z",
"shell.execute_reply": "2024-02-19T13:48:50.381841Z"
}
},
"outputs": [],
"source": [
"A = A.diags('lil')\n",
"A[-2] = (-1)**np.arange(N)\n",
"A[-1] = np.ones(N)\n",
"A = A.tocsc()\n",
"f_hat[-2] = ue.subs(x, T.domain[0])\n",
"f_hat[-1] = ue.subs(x, T.domain[1])"
]
},
{
"cell_type": "markdown",
"id": "0e3b7fe8",
"metadata": {
"editable": true
},
"source": [
"Note that the last two lines uses evaluation of the sympy function\n",
"`ue` at the borders of the domain. Implemented like this it is\n",
"easy to change to a nonstandard domain size.\n",
"The sparsity pattern of the matrix A is now modified\n",
"with the typical tau-lines that we can visualize using [plotly](https://plotly.com/)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "454473fc",
"metadata": {
"collapsed": false,
"editable": true,
"execution": {
"iopub.execute_input": "2024-02-19T13:48:50.385039Z",
"iopub.status.busy": "2024-02-19T13:48:50.384822Z",
"iopub.status.idle": "2024-02-19T13:48:51.080843Z",
"shell.execute_reply": "2024-02-19T13:48:51.080389Z"
}
},
"outputs": [
{
"data": {
"text/html": [
" \n",
" "
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"application/vnd.plotly.v1+json": {
"config": {
"plotlyServerURL": "https://plot.ly"
},
"data": [
{
"hovertemplate": "x: %{x}
y: %{y}
color: %{z[0]}",
"name": "0",
"source": "data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAACAAAAAgCAAAAABWESUoAAAAbElEQVR4XpXMQRLAMAgCQNL//9ke0iSKGqZwwx1hBlEhYEqYEBPwGirEAry7CnEAX3aFeHBio6gHpQigEhEUgkAWDJJIgEUGJAoQRQWCKIEXNXCiAUd0YIsWLNGDT1zAFDcA9QFQHwAbvPzPCy6gNXA++w1YAAAAAElFTkSuQmCC",
"type": "image",
"xaxis": "x",
"yaxis": "y"
}
],
"layout": {
"margin": {
"t": 60
},
"template": {
"data": {
"bar": [
{
"error_x": {
"color": "#2a3f5f"
},
"error_y": {
"color": "#2a3f5f"
},
"marker": {
"line": {
"color": "#E5ECF6",
"width": 0.5
},
"pattern": {
"fillmode": "overlay",
"size": 10,
"solidity": 0.2
}
},
"type": "bar"
}
],
"barpolar": [
{
"marker": {
"line": {
"color": "#E5ECF6",
"width": 0.5
},
"pattern": {
"fillmode": "overlay",
"size": 10,
"solidity": 0.2
}
},
"type": "barpolar"
}
],
"carpet": [
{
"aaxis": {
"endlinecolor": "#2a3f5f",
"gridcolor": "white",
"linecolor": "white",
"minorgridcolor": "white",
"startlinecolor": "#2a3f5f"
},
"baxis": {
"endlinecolor": "#2a3f5f",
"gridcolor": "white",
"linecolor": "white",
"minorgridcolor": "white",
"startlinecolor": "#2a3f5f"
},
"type": "carpet"
}
],
"choropleth": [
{
"colorbar": {
"outlinewidth": 0,
"ticks": ""
},
"type": "choropleth"
}
],
"contour": [
{
"colorbar": {
"outlinewidth": 0,
"ticks": ""
},
"colorscale": [
[
0.0,
"#0d0887"
],
[
0.1111111111111111,
"#46039f"
],
[
0.2222222222222222,
"#7201a8"
],
[
0.3333333333333333,
"#9c179e"
],
[
0.4444444444444444,
"#bd3786"
],
[
0.5555555555555556,
"#d8576b"
],
[
0.6666666666666666,
"#ed7953"
],
[
0.7777777777777778,
"#fb9f3a"
],
[
0.8888888888888888,
"#fdca26"
],
[
1.0,
"#f0f921"
]
],
"type": "contour"
}
],
"contourcarpet": [
{
"colorbar": {
"outlinewidth": 0,
"ticks": ""
},
"type": "contourcarpet"
}
],
"heatmap": [
{
"colorbar": {
"outlinewidth": 0,
"ticks": ""
},
"colorscale": [
[
0.0,
"#0d0887"
],
[
0.1111111111111111,
"#46039f"
],
[
0.2222222222222222,
"#7201a8"
],
[
0.3333333333333333,
"#9c179e"
],
[
0.4444444444444444,
"#bd3786"
],
[
0.5555555555555556,
"#d8576b"
],
[
0.6666666666666666,
"#ed7953"
],
[
0.7777777777777778,
"#fb9f3a"
],
[
0.8888888888888888,
"#fdca26"
],
[
1.0,
"#f0f921"
]
],
"type": "heatmap"
}
],
"heatmapgl": [
{
"colorbar": {
"outlinewidth": 0,
"ticks": ""
},
"colorscale": [
[
0.0,
"#0d0887"
],
[
0.1111111111111111,
"#46039f"
],
[
0.2222222222222222,
"#7201a8"
],
[
0.3333333333333333,
"#9c179e"
],
[
0.4444444444444444,
"#bd3786"
],
[
0.5555555555555556,
"#d8576b"
],
[
0.6666666666666666,
"#ed7953"
],
[
0.7777777777777778,
"#fb9f3a"
],
[
0.8888888888888888,
"#fdca26"
],
[
1.0,
"#f0f921"
]
],
"type": "heatmapgl"
}
],
"histogram": [
{
"marker": {
"pattern": {
"fillmode": "overlay",
"size": 10,
"solidity": 0.2
}
},
"type": "histogram"
}
],
"histogram2d": [
{
"colorbar": {
"outlinewidth": 0,
"ticks": ""
},
"colorscale": [
[
0.0,
"#0d0887"
],
[
0.1111111111111111,
"#46039f"
],
[
0.2222222222222222,
"#7201a8"
],
[
0.3333333333333333,
"#9c179e"
],
[
0.4444444444444444,
"#bd3786"
],
[
0.5555555555555556,
"#d8576b"
],
[
0.6666666666666666,
"#ed7953"
],
[
0.7777777777777778,
"#fb9f3a"
],
[
0.8888888888888888,
"#fdca26"
],
[
1.0,
"#f0f921"
]
],
"type": "histogram2d"
}
],
"histogram2dcontour": [
{
"colorbar": {
"outlinewidth": 0,
"ticks": ""
},
"colorscale": [
[
0.0,
"#0d0887"
],
[
0.1111111111111111,
"#46039f"
],
[
0.2222222222222222,
"#7201a8"
],
[
0.3333333333333333,
"#9c179e"
],
[
0.4444444444444444,
"#bd3786"
],
[
0.5555555555555556,
"#d8576b"
],
[
0.6666666666666666,
"#ed7953"
],
[
0.7777777777777778,
"#fb9f3a"
],
[
0.8888888888888888,
"#fdca26"
],
[
1.0,
"#f0f921"
]
],
"type": "histogram2dcontour"
}
],
"mesh3d": [
{
"colorbar": {
"outlinewidth": 0,
"ticks": ""
},
"type": "mesh3d"
}
],
"parcoords": [
{
"line": {
"colorbar": {
"outlinewidth": 0,
"ticks": ""
}
},
"type": "parcoords"
}
],
"pie": [
{
"automargin": true,
"type": "pie"
}
],
"scatter": [
{
"fillpattern": {
"fillmode": "overlay",
"size": 10,
"solidity": 0.2
},
"type": "scatter"
}
],
"scatter3d": [
{
"line": {
"colorbar": {
"outlinewidth": 0,
"ticks": ""
}
},
"marker": {
"colorbar": {
"outlinewidth": 0,
"ticks": ""
}
},
"type": "scatter3d"
}
],
"scattercarpet": [
{
"marker": {
"colorbar": {
"outlinewidth": 0,
"ticks": ""
}
},
"type": "scattercarpet"
}
],
"scattergeo": [
{
"marker": {
"colorbar": {
"outlinewidth": 0,
"ticks": ""
}
},
"type": "scattergeo"
}
],
"scattergl": [
{
"marker": {
"colorbar": {
"outlinewidth": 0,
"ticks": ""
}
},
"type": "scattergl"
}
],
"scattermapbox": [
{
"marker": {
"colorbar": {
"outlinewidth": 0,
"ticks": ""
}
},
"type": "scattermapbox"
}
],
"scatterpolar": [
{
"marker": {
"colorbar": {
"outlinewidth": 0,
"ticks": ""
}
},
"type": "scatterpolar"
}
],
"scatterpolargl": [
{
"marker": {
"colorbar": {
"outlinewidth": 0,
"ticks": ""
}
},
"type": "scatterpolargl"
}
],
"scatterternary": [
{
"marker": {
"colorbar": {
"outlinewidth": 0,
"ticks": ""
}
},
"type": "scatterternary"
}
],
"surface": [
{
"colorbar": {
"outlinewidth": 0,
"ticks": ""
},
"colorscale": [
[
0.0,
"#0d0887"
],
[
0.1111111111111111,
"#46039f"
],
[
0.2222222222222222,
"#7201a8"
],
[
0.3333333333333333,
"#9c179e"
],
[
0.4444444444444444,
"#bd3786"
],
[
0.5555555555555556,
"#d8576b"
],
[
0.6666666666666666,
"#ed7953"
],
[
0.7777777777777778,
"#fb9f3a"
],
[
0.8888888888888888,
"#fdca26"
],
[
1.0,
"#f0f921"
]
],
"type": "surface"
}
],
"table": [
{
"cells": {
"fill": {
"color": "#EBF0F8"
},
"line": {
"color": "white"
}
},
"header": {
"fill": {
"color": "#C8D4E3"
},
"line": {
"color": "white"
}
},
"type": "table"
}
]
},
"layout": {
"annotationdefaults": {
"arrowcolor": "#2a3f5f",
"arrowhead": 0,
"arrowwidth": 1
},
"autotypenumbers": "strict",
"coloraxis": {
"colorbar": {
"outlinewidth": 0,
"ticks": ""
}
},
"colorscale": {
"diverging": [
[
0,
"#8e0152"
],
[
0.1,
"#c51b7d"
],
[
0.2,
"#de77ae"
],
[
0.3,
"#f1b6da"
],
[
0.4,
"#fde0ef"
],
[
0.5,
"#f7f7f7"
],
[
0.6,
"#e6f5d0"
],
[
0.7,
"#b8e186"
],
[
0.8,
"#7fbc41"
],
[
0.9,
"#4d9221"
],
[
1,
"#276419"
]
],
"sequential": [
[
0.0,
"#0d0887"
],
[
0.1111111111111111,
"#46039f"
],
[
0.2222222222222222,
"#7201a8"
],
[
0.3333333333333333,
"#9c179e"
],
[
0.4444444444444444,
"#bd3786"
],
[
0.5555555555555556,
"#d8576b"
],
[
0.6666666666666666,
"#ed7953"
],
[
0.7777777777777778,
"#fb9f3a"
],
[
0.8888888888888888,
"#fdca26"
],
[
1.0,
"#f0f921"
]
],
"sequentialminus": [
[
0.0,
"#0d0887"
],
[
0.1111111111111111,
"#46039f"
],
[
0.2222222222222222,
"#7201a8"
],
[
0.3333333333333333,
"#9c179e"
],
[
0.4444444444444444,
"#bd3786"
],
[
0.5555555555555556,
"#d8576b"
],
[
0.6666666666666666,
"#ed7953"
],
[
0.7777777777777778,
"#fb9f3a"
],
[
0.8888888888888888,
"#fdca26"
],
[
1.0,
"#f0f921"
]
]
},
"colorway": [
"#636efa",
"#EF553B",
"#00cc96",
"#ab63fa",
"#FFA15A",
"#19d3f3",
"#FF6692",
"#B6E880",
"#FF97FF",
"#FECB52"
],
"font": {
"color": "#2a3f5f"
},
"geo": {
"bgcolor": "white",
"lakecolor": "white",
"landcolor": "#E5ECF6",
"showlakes": true,
"showland": true,
"subunitcolor": "white"
},
"hoverlabel": {
"align": "left"
},
"hovermode": "closest",
"mapbox": {
"style": "light"
},
"paper_bgcolor": "white",
"plot_bgcolor": "#E5ECF6",
"polar": {
"angularaxis": {
"gridcolor": "white",
"linecolor": "white",
"ticks": ""
},
"bgcolor": "#E5ECF6",
"radialaxis": {
"gridcolor": "white",
"linecolor": "white",
"ticks": ""
}
},
"scene": {
"xaxis": {
"backgroundcolor": "#E5ECF6",
"gridcolor": "white",
"gridwidth": 2,
"linecolor": "white",
"showbackground": true,
"ticks": "",
"zerolinecolor": "white"
},
"yaxis": {
"backgroundcolor": "#E5ECF6",
"gridcolor": "white",
"gridwidth": 2,
"linecolor": "white",
"showbackground": true,
"ticks": "",
"zerolinecolor": "white"
},
"zaxis": {
"backgroundcolor": "#E5ECF6",
"gridcolor": "white",
"gridwidth": 2,
"linecolor": "white",
"showbackground": true,
"ticks": "",
"zerolinecolor": "white"
}
},
"shapedefaults": {
"line": {
"color": "#2a3f5f"
}
},
"ternary": {
"aaxis": {
"gridcolor": "white",
"linecolor": "white",
"ticks": ""
},
"baxis": {
"gridcolor": "white",
"linecolor": "white",
"ticks": ""
},
"bgcolor": "#E5ECF6",
"caxis": {
"gridcolor": "white",
"linecolor": "white",
"ticks": ""
}
},
"title": {
"x": 0.05
},
"xaxis": {
"automargin": true,
"gridcolor": "white",
"linecolor": "white",
"ticks": "",
"title": {
"standoff": 15
},
"zerolinecolor": "white",
"zerolinewidth": 2
},
"yaxis": {
"automargin": true,
"gridcolor": "white",
"linecolor": "white",
"ticks": "",
"title": {
"standoff": 15
},
"zerolinecolor": "white",
"zerolinewidth": 2
}
}
},
"xaxis": {
"anchor": "y",
"domain": [
0.0,
1.0
]
},
"yaxis": {
"anchor": "x",
"domain": [
0.0,
1.0
]
}
}
},
"text/html": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"import plotly.express as px\n",
"z = np.where(abs(A.toarray()) > 1e-6, 0, 1).astype(bool)\n",
"fig = px.imshow(z, binary_string=True)\n",
"fig.show()"
]
},
{
"cell_type": "markdown",
"id": "5a1ae665",
"metadata": {
"editable": true
},
"source": [
"### Solve linear equations\n",
"\n",
"Finally, solve the linear equation system and transform the solution from the spectral\n",
"$\\{\\hat{u}_k\\}_{k=0}^{N-1}$ 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": 7,
"id": "c55f7862",
"metadata": {
"collapsed": false,
"editable": true,
"execution": {
"iopub.execute_input": "2024-02-19T13:48:51.106078Z",
"iopub.status.busy": "2024-02-19T13:48:51.105386Z",
"iopub.status.idle": "2024-02-19T13:48:51.117236Z",
"shell.execute_reply": "2024-02-19T13:48:51.116694Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Error=6.1627552156954840e-11\n"
]
}
],
"source": [
"u_hat = Function(T)\n",
"u_hat[:] = scp.linalg.spsolve(A, f_hat)\n",
"uj = u_hat.backward()\n",
"ua = Array(T, buffer=ue)\n",
"\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": 8,
"id": "c5417c91",
"metadata": {
"collapsed": false,
"editable": true,
"execution": {
"iopub.execute_input": "2024-02-19T13:48:51.120542Z",
"iopub.status.busy": "2024-02-19T13:48:51.120304Z",
"iopub.status.idle": "2024-02-19T13:48:51.125610Z",
"shell.execute_reply": "2024-02-19T13:48:51.125085Z"
}
},
"outputs": [],
"source": [
"def main(N, family='Chebyshev'):\n",
" T = FunctionSpace(N, family=family)\n",
" u = TrialFunction(T)\n",
" v = TestFunction(T)\n",
"\n",
" # Get f on quad points\n",
" fj = Array(T, buffer=fe)\n",
"\n",
" # Compute right hand side of Poisson's equation\n",
" f_hat = Function(T)\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",
" A = A.diags('lil')\n",
" A[-2] = (-1)**np.arange(N)\n",
" A[-1] = np.ones(N)\n",
" A = A.tocsc()\n",
" f_hat[-2] = ue.subs(x, T.domain[0])\n",
" f_hat[-1] = ue.subs(x, T.domain[1])\n",
"\n",
" u_hat = Function(T)\n",
" u_hat[:] = scp.linalg.spsolve(A, f_hat)\n",
" uj = u_hat.backward()\n",
"\n",
" # Compare with analytical solution\n",
" ua = Array(T, 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": 9,
"id": "e2444fa7",
"metadata": {
"collapsed": false,
"editable": true,
"execution": {
"iopub.execute_input": "2024-02-19T13:48:51.128290Z",
"iopub.status.busy": "2024-02-19T13:48:51.128108Z",
"iopub.status.idle": "2024-02-19T13:48:51.140165Z",
"shell.execute_reply": "2024-02-19T13:48:51.139738Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"1.9741838920185597"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"main(12, 'Chebyshev')"
]
},
{
"cell_type": "markdown",
"id": "e42b9a11",
"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)."
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "1215522e",
"metadata": {
"collapsed": false,
"editable": true,
"execution": {
"iopub.execute_input": "2024-02-19T13:48:51.143099Z",
"iopub.status.busy": "2024-02-19T13:48:51.142906Z",
"iopub.status.idle": "2024-02-19T13:48:52.107805Z",
"shell.execute_reply": "2024-02-19T13:48:52.107155Z"
}
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAikAAAGHCAYAAAB1bcIdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAABmg0lEQVR4nO3dd3wUxf/H8delF0LogYRQRHqAQBAEpKkEgjQLKCJdFEHpKIiKnSIgKiICUiwIIkVUqj+qoFKjVBWIgAKGHgiQOr8/7stJSIAkJOwleT8fj33ozc7tfuY25D6ZnZ2xGWMMIiIiIk7GxeoARERERNKiJEVERESckpIUERERcUpKUkRERMQpKUkRERERp6QkRURERJySkhQRERFxSkpSRERExCkpSRERERGnpCRFstVvv/1G9+7dKVu2LF5eXuTLl49atWoxduxYTp8+bXV4ksXmzZtH1apV8fb2xmazERkZmapOmTJlsNlsN91mzZqV7fE2adIkxTm9vb2pUaMGEydOJDk5OUPH+uuvv25b3M4up34W58+f5/nnnyc8PJyiRYtis9l49dVX06x79c+Oi4sLfn5+3HnnnbRv356vv/46wz8/kjY3qwOQ3GvatGn06dOHihUrMnToUKpUqUJCQgJbt25lypQp/PTTTyxatMjqMCWLnDhxgs6dO9OiRQsmT56Mp6cnFSpUSFVv0aJFxMXFOV5Pnz6dTz75hOXLl+Pv7+8oL1eu3G2J+4477uCLL74AIDo6milTpjBw4ECOHTvGmDFj0n2cEiVK8NNPP922uCXrnTp1iqlTp1KjRg3atWvH9OnTb1j/6p+d2NhYoqKiWLx4Me3bt6dhw4Z8++23KX6mJROMSDbYtGmTcXV1NS1atDCXL19OtT8uLs588803FkSWdRITE9NsW171448/GsDMmzcvQ+8bOXKkAcyJEyeyKbLra9y4salatWqKsvj4eHPHHXcYHx8fEx8ff9tjyg2ioqIMYGbOnHlbz3vx4kWTnJyc6fcnJyc73n/ixAkDmJEjR6ZZN62fnStmzJhhANOhQ4dMxyJ2ut0j2eLtt9/GZrMxdepUPD09U+338PCgTZs2jtfJycmMHTuWSpUq4enpSbFixejSpQt///13ivc1adKEkJAQtmzZQsOGDfHx8eGOO+5g9OjRju7VEydO4OHhwcsvv5zqvPv27cNms/H+++87yo4fP87TTz9NyZIl8fDwoGzZsrz22mskJiY66lzpvh47dixvvvkmZcuWxdPTkzVr1gDwzTffUL16dTw9Pbnjjjt47733ePXVV7HZbCnOb4xh8uTJhIaG4u3tTcGCBXnkkUc4ePBghtt5xdmzZxk8eDB33HGH47Nr2bIl+/btc9SJj4/nzTffdHy+RYsWpXv37pw4cSLtC3iNJUuWUK9ePXx8fPDz86NZs2b89NNPjv3dunXjnnvuAeDRRx/FZrPRpEmTdB07LfPmzSM8PJwSJUrg7e1N5cqVGTZsGLGxsSnqNWnSJM3zdOvWjTJlymTq3O7u7oSFhXHx4kXH57Nr1y7atm1LwYIF8fLyIjQ0lNmzZ6d4X1q3OE6cOMFTTz1FcHCw43Nv0KABP/zwg6POjh07aNWqFcWKFcPT05PAwEAeeOCBFD/7ly9fZvjw4ZQtWxYPDw+CgoLo27cvZ8+eTRFDmTJlaNWqFcuXL6dWrVp4e3tTqVIlZsyYka62f/TRR9SoUYN8+fLh5+dHpUqVePHFF1PUSc9nca3Fixdjs9n4v//7vzTPabPZ+O233xxlW7dupU2bNhQqVAgvLy9q1qzJV199leJ9s2bNwmazsXLlSnr06EHRokXx8fEhLi4uXZ97Wq7cvrlV3bt3p2XLlsyfP59Dhw7d8vHyNKuzJMl9EhMTjY+Pj6lbt2663/PUU08ZwDz77LNm+fLlZsqUKaZo0aImODg4xV/YjRs3NoULFzbly5c3U6ZMMatWrTJ9+vQxgJk9e7aj3oMPPmiCg4NNUlJSivM8//zzxsPDw5w8edIYY8yxY8dMcHCwKV26tPn444/NDz/8YN544w3j6elpunXr5njflb8Mg4KCTNOmTc3XX39tVq5caaKiosyyZcuMi4uLadKkiVm0aJGZP3++qVu3rilTpoy59p9Yr169jLu7uxk8eLBZvny5mTNnjqlUqZIJCAgwx48fz3A7Y2JiTNWqVY2vr695/fXXzYoVK8yCBQtM//79zerVq40xxiQlJZkWLVoYX19f89prr5lVq1aZ6dOnm6CgIFOlShVz8eLFG16bL774wgAmPDzcLF682MybN8+EhYUZDw8Ps2HDBmOMMfv37zcffvihAczbb79tfvrpJ7N79+50Xfu0elLeeOMN8+6775rvv//erF271kyZMsWULVvWNG3aNMV7GzdubBo3bpzqmF27djWlS5e+6bmv99dwrVq1jJubm7l48aLZt2+f8fPzM+XKlTOffvqp+f77703Hjh0NYMaMGeN4T1q9B82bNzdFixY1U6dONWvXrjWLFy82r7zyipk7d64xxpgLFy6YwoULm9q1a5uvvvrKrFu3zsybN8/07t3b7Nmzxxhj/+u+efPmxs3Nzbz88stm5cqVZty4ccbX19fUrFkzRW9e6dKlTcmSJU2VKlXMp59+alasWGHat29vALNu3bobfhZffvmlAcxzzz1nVq5caX744QczZcoU069fP0edzH4WCQkJplixYqZTp06pzlunTh1Tq1Ytx+vVq1cbDw8P07BhQzNv3jyzfPly061bt1Sf7cyZMx3/Jp966imzbNky8/XXX5vExMSbfu7pcSs9KcYYM2XKFAOYzz77LN3nlNSUpEiWO378uAHMY489lq76e/fuNYDp06dPivJffvnFAObFF190lDVu3NgA5pdffklRt0qVKqZ58+aO10uWLDGAWblypaMsMTHRBAYGmocffthR9vTTT5t8+fKZQ4cOpTjeuHHjDOD4or3yS7dcuXKpbgHcddddJjg42MTFxTnKzp8/bwoXLpwiSfnpp58MYMaPH5/i/UeOHDHe3t7m+eefz3A7X3/9dQOYVatWmeu58uWzYMGCFOVbtmwxgJk8efJ135uUlGQCAwNNtWrVUiR858+fN8WKFTP169d3lK1Zs8YAZv78+dc9XlpudrsnOTnZJCQkmHXr1hnA/Prrr459WZWkJCQkmISEBHP06FEzbNgwA5j27dsbY4x57LHHjKenpzl8+HCK90ZERBgfHx9z9uxZY0zaSUq+fPnMgAEDrnv+rVu3GsAsXrz4unWWL19uADN27NgU5fPmzTOAmTp1qqOsdOnSxsvLK8XP86VLl0yhQoXM008/fcPP4tlnnzUFChS4YZ1b+SwGDRpkvL29HXWMMWbPnj0GMB988IGjrFKlSqZmzZomISEhxTlatWplSpQo4fg5vJKkdOnSJVWcN/vc0+NWk5Rly5alSt4k43S7Ryx35ZZJt27dUpTXqVOHypUrp+oiLl68OHXq1ElRVr169RTdqhERERQvXpyZM2c6ylasWMHRo0fp0aOHo+y7776jadOmBAYGkpiY6NgiIiIAWLduXYrztGnTBnd3d8fr2NhYtm7dSrt27fDw8HCU58uXj9atW6d473fffYfNZuOJJ55Ica7ixYtTo0YN1q5dm+F2Llu2jAoVKnD//fdzPd999x0FChSgdevWKc4bGhpK8eLFU533ar///jtHjx6lc+fOuLj89+siX758PPzww/z8889cvHjxuu/PrIMHD/L4449TvHhxXF1dcXd3p3HjxgDs3bs3S8+1e/du3N3dcXd3JzAwkPHjx9OpUyemTZsGwOrVq7nvvvsIDg5O8b5u3bpx8eLFFLe9rlWnTh1mzZrFm2++yc8//0xCQkKK/XfeeScFCxbkhRdeYMqUKezZsyfVMVavXu0439Xat2+Pr69vqn8foaGhlCpVyvHay8uLChUq3PS2Q506dTh79iwdO3bkm2++4eTJk2nGktnPokePHly6dIl58+Y5ymbOnImnpyePP/44APv372ffvn106tQJIMXPa8uWLTl27Bi///57iuM+/PDDabblRp/77WCMue3nzI2UpEiWK1KkCD4+PkRFRaWr/qlTpwD70xHXCgwMdOy/onDhwqnqeXp6cunSJcdrNzc3OnfuzKJFixz37WfNmkWJEiVo3ry5o96///7Lt99+6/iSurJVrVoVINUv6mtjPHPmDMYYAgICUsV0bdm///7rqHvt+X7++edU50pPO0+cOEHJkiVT1bv2vGfPnsXDwyPVeY8fP57ml9EVN7s2ycnJnDlz5obnz6gLFy7QsGFDfvnlF958803Wrl3Lli1bWLhwIUCK9meFcuXKsWXLFrZu3cquXbs4e/Ysn3/+ueOpjFOnTl23/Vf2X8+8efPo2rUr06dPp169ehQqVIguXbpw/PhxAPz9/Vm3bh2hoaG8+OKLVK1alcDAQEaOHOn4Yj116hRubm4ULVo0xbFtNhvFixfP1L+PtHTu3JkZM2Zw6NAhHn74YYoVK0bdunVZtWqVo86tfBZVq1blrrvucvzhkJSUxOeff07btm0pVKgQYP9ZBRgyZEiqn9U+ffoAN/83CTf/3G+HK0nhlc9GMkePIEuWc3V15b777mPZsmX8/fffN/0SvfJL9dixY6nqHj16lCJFimQqju7du/POO+8wd+5cHn30UZYsWcKAAQNwdXV11ClSpAjVq1fnrbfeSvMY1/6CuXZQXcGCBbHZbI5frle79hdikSJFsNlsbNiwIc3BxGmV3UzRokVTDS6+VpEiRShcuDDLly9Pc7+fn99133v1tbnW0aNHcXFxoWDBghmI+OZWr17N0aNHWbt2raP3BEg1SBTsvQTnzp1LVX6jxCutY9SuXfu6+wsXLnzd9gM3/PksUqQIEydOZOLEiRw+fJglS5YwbNgwoqOjHdejWrVqzJ07F2MMv/32G7NmzeL111/H29ubYcOGUbhwYRITEzlx4kSKRMUYw/Hjx7nrrrvS3dab6d69O927dyc2Npb169czcuRIWrVqxR9//EHp0qVv6bO4cvw+ffqwd+9eDh48yLFjx+jevbtj/5X3Dx8+nIceeijNY1SsWDHF67QGuqbnc89uS5YswWaz0ahRo9tyvtxKPSmSLYYPH44xhl69ehEfH59qf0JCAt9++y0A9957LwCff/55ijpbtmxh79693HfffZmKoXLlytStW5eZM2cyZ84c4uLiUvxCBGjVqhW7du2iXLly1K5dO9V2s7+CfH19qV27NosXL07RzgsXLvDdd9+lOpcxhn/++SfNc1WrVi3DbYyIiOCPP/5w3BJIS6tWrTh16hRJSUlpnvfaX/pXq1ixIkFBQcyZMydF93VsbCwLFixwPPGTla586VybtH388cep6pYpU4Y//vgjxbwrp06dYtOmTVkWz3333edInK726aef4uPjw913352u45QqVYpnn32WZs2asX379lT7bTYbNWrU4N1336VAgQKOOld+/q/997FgwQJiY2Mz/e/jRnx9fYmIiGDEiBHEx8eze/duRyy38ll07NgRLy8vZs2axaxZswgKCiI8PNyxv2LFipQvX55ff/01zZ/V2rVr3zCpTsvNPvfsMHPmTJYtW0bHjh1T3HqTjFNPimSLevXq8dFHH9GnTx/CwsJ45plnqFq1KgkJCezYsYOpU6cSEhJC69atqVixIk899RQffPABLi4uRERE8Ndff/Hyyy8THBzMwIEDMx1Hjx49ePrppzl69Cj169dP9YX8+uuvs2rVKurXr0+/fv2oWLEily9f5q+//mLp0qVMmTLlpj1Br7/+Og888ADNmzenf//+JCUl8c4775AvX74Us+o2aNCAp556iu7du7N161YaNWqEr68vx44d48cff6RatWo888wzGWrfgAEDmDdvHm3btmXYsGHUqVOHS5cusW7dOlq1akXTpk157LHH+OKLL2jZsiX9+/enTp06uLu78/fff7NmzRratm3Lgw8+mObxXVxcGDt2LJ06daJVq1Y8/fTTxMXF8c4773D27FlGjx6doXjTo379+hQsWJDevXszcuRI3N3d+eKLL/j1119T1e3cuTMff/wxTzzxBL169eLUqVOMHTuW/PnzZ1k8I0eOdIxdeuWVVyhUqBBffPEF33//PWPHjr3uZF3nzp2jadOmPP7441SqVAk/Pz+2bNnC8uXLHb0E3333HZMnT6Zdu3bccccdGGNYuHAhZ8+epVmzZgA0a9aM5s2b88ILLxATE0ODBg347bffGDlyJDVr1qRz585Z0s5evXrh7e1NgwYNKFGiBMePH2fUqFH4+/s7emsy+1lcUaBAAR588EFmzZrF2bNnGTJkSIqxTmBPRiMiImjevDndunUjKCiI06dPs3fvXrZv3878+fNveI70fO43smzZMmJjYzl//jwAe/bs4euvvwagZcuWKZLyS5cu8fPPPzv+/+DBgyxevJjvvvuOxo0bM2XKlJueT27CsiG7kidERkaarl27mlKlShkPDw/HY5OvvPKKiY6OdtRLSkoyY8aMMRUqVDDu7u6mSJEi5oknnjBHjhxJcbzrjai/3tMc586dM97e3gYw06ZNSzPGEydOmH79+pmyZcsad3d3U6hQIRMWFmZGjBhhLly4YIz572mFd955J81jLFq0yFSrVs14eHiYUqVKmdGjR5t+/fqZggULpqo7Y8YMU7duXePr62u8vb1NuXLlTJcuXczWrVsz1c4zZ86Y/v37m1KlShl3d3dTrFgx88ADD5h9+/Y56iQkJJhx48aZGjVqGC8vL5MvXz5TqVIl8/TTT5s///wzzTZdbfHixaZu3brGy8vL+Pr6mvvuu89s3LgxRZ2sfLpn06ZNpl69esbHx8cULVrUPPnkk2b79u1pThA2e/ZsU7lyZePl5WWqVKli5s2bd8uPIF9r586dpnXr1sbf3994eHiYGjVqpIrj2idaLl++bHr37m2qV69u8ufPb7y9vU3FihXNyJEjTWxsrDHG/khvx44dTbly5Yy3t7fx9/c3derUMbNmzUpx7EuXLpkXXnjBlC5d2ri7u5sSJUqYZ555xpw5cyZFvdKlS5sHHnggzXam9RTU1WbPnm2aNm1qAgICjIeHhwkMDDQdOnQwv/322y1/FldbuXKlAQxg/vjjjzRj+fXXX02HDh1MsWLFjLu7uylevLi59957zZQpUxx1rjzds2XLlhTvTc/nfiOlS5d2xHftFhUV5ah35Qm8K5uvr6+54447zCOPPGLmz5+favoDyRybMRqCLJLVEhISCA0NJSgoiJUrV1odjohIjqTbPSJZoGfPnjRr1szRTT5lyhT27t3Le++9Z3VoIiI5lpIUkSxw/vx5hgwZwokTJ3B3d6dWrVosXbr0hvOXiIjIjel2j4iIiDglPYIsIiIiTklJioiIiDglJSkiIiLilDRwNpOSk5M5evQofn5+aU7LLCIiImkzxnD+/HkCAwNTTeh3NSUpmXT06NFUK4GKiIhI+h05cuSGs3orScmkK+tHHDlyJEun4BYREcntYmJiCA4OvulaTEpSMunKLZ78+fMrSREREcmEmw2X0MBZERERcUpKUkRERMQpKUkRERERp6QxKSIiku2MMSQmJpKUlGR1KHIbuLq64ubmdstTdChJERGRbBUfH8+xY8e4ePGi1aHIbeTj40OJEiXw8PDI9DGUpIiISLZJTk4mKioKV1dXAgMD8fDw0ASYuZwxhvj4eE6cOEFUVBTly5e/4YRtN6IkRUREsk18fDzJyckEBwfj4+NjdThym3h7e+Pu7s6hQ4eIj4/Hy8srU8fRwFkREcl2mf1LWnKurLjm+qlxIvHxVkcgIiLiPJSkOIn4eKhaFZ58Ev74w+poRERErKckxUnMmQP798Mnn0ClStChA+zYYXVUIiKS03Tr1o127dpZHUaWyNNJynfffUfFihUpX74806dPtzSW2MOn8OcsAMbA/PlQqxa0aAHr1tnLRETk9slNX/Y5VZ5NUhITExk0aBCrV69m+/btjBkzhtOnT1sWT9/jL3OI0ozmBQI47ihfsQKaNIEGDeDbb5WsiIjIrYnPQQMg82ySsnnzZqpWrUpQUBB+fn60bNmSFStWWBfQiBH49+/OC17vE0VZJvMMZYhy7P7pJ2jTBqpXt98aSky0LlQRkbxuz549tGzZknz58hEQEEDnzp05efKkY//58+fp1KkTvr6+lChRgnfffZcmTZowYMAAR534+Hief/55goKC8PX1pW7duqxdu9axf9asWRQoUIAVK1ZQuXJl8uXLR4sWLTh27JijTlJSEoMGDaJAgQIULlyY559/HnPNX7NNmjTh2WefZdCgQRQpUoRmzZqlqw3OIMcmKevXr6d169YEBgZis9lYvHhxqjqTJ0+mbNmyeHl5ERYWxoYNGxz7jh49SlBQkON1yZIl+eeff25H6GkLCoKJEyEqCu8hz/KMz6f8SXk+4wmqsstRbdcu6NQJKlaEKVPg8mXrQhYRybTataFkydu/1a59y6EfO3aMxo0bExoaytatW1m+fDn//vsvHTp0cNQZNGgQGzduZMmSJaxatYoNGzawffv2FMfp3r07GzduZO7cufz222+0b9+eFi1a8OeffzrqXLx4kXHjxvHZZ5+xfv16Dh8+zJAhQxz7x48fz4wZM/jkk0/48ccfOX36NIsWLUoV8+zZs3Fzc2Pjxo18/PHH6WqDUzA51NKlS82IESPMggULDGAWLVqUYv/cuXONu7u7mTZtmtmzZ4/p37+/8fX1NYcOHTLGGPPVV1+Zvn37OuqPHTvWjBs3Lt3nP3funAHMuXPnsqQ9qURHGzNsmDH58pkkbOYbWpu72WTsN3z+24oXN2bMGGOyKwwRkVtx6dIls2fPHnPp0qWUO4KCTKpfaLdjCwpKd+xdu3Y1bdu2TVX+8ssvm/Dw8BRlR44cMYD5/fffTUxMjHF3dzfz58937D979qzx8fEx/fv3N8YYs3//fmOz2cw///yT4jj33XefGT58uDHGmJkzZxrA7N+/37H/ww8/NAEBAY7XJUqUMKNHj3a8TkhIMCVLlkwRd+PGjU1oaGiG2pAVrnvtTfq/Q3PsjLMRERFERERcd/+ECRPo2bMnTz75JAATJ05kxYoVfPTRR4waNYqgoKAUPSd///03devWve7x4uLiiIuLc7yOiYnJglbcQNGiMGoUDBmCy8SJtHn/fVrHfMs6GjOK4aykOQDHj8MLL9ir9u0L/fvb3yoi4tSKF8+x5922bRtr1qwhX758qfYdOHCAS5cukZCQQJ06dRzl/v7+VKxY0fF6+/btGGOoUKFCivfHxcVRuHBhx2sfHx/KlSvneF2iRAmio6MBOHfuHMeOHaNevXqO/W5ubtSuXTvVLZ/a1/Qg3awN18ZllRybpNxIfHw827ZtY9iwYSnKw8PD2bRpEwB16tRh165d/PPPP+TPn5+lS5fyyiuvXPeYo0aN4rXXXsvWuNNUuDC88QYMGoTt/fdpMnEiTc62YBu1GM0wFvAwBhfOnoW33oIJE6BXLxg8GEqVuv3hioiky9atVkeQacnJybRu3ZoxY8ak2leiRAnH7Zpr1yi6OnFITk7G1dWVbdu24erqmqLe1YmDu7t7in02my1VApIevr6+GWqDs8ixY1Ju5OTJkyQlJREQEJCiPCAggOPH7U/OuLm5MX78eJo2bUrNmjUZOnRoiuz1WsOHD+fcuXOO7ciRI9nahlQKFoSRI+Gvv+DNNwkr9Bfz6cBeKtOdGbiRAMClS/D++1CuHHTvDvv23d4wRURyu1q1arF7927KlCnDnXfemWLz9fWlXLlyuLu7s3nzZsd7YmJiUow1qVmzJklJSURHR6c6RvF09vb4+/tTokQJfv75Z0dZYmIi27Ztu+U2OItcmaRckVYWe3VZmzZt+OOPP9i/fz9PPfXUDY/l6elJ/vz5U2yW8PeHESPsycro0VQscpoZ9OQgd9CfifgQC9if/pk1C6pUgYcfztF/tIiIWObcuXNERkam2J5++mlOnz5Nx44d2bx5MwcPHmTlypX06NGDpKQk/Pz86Nq1K0OHDmXNmjXs3r2bHj164OLi4vgOqlChAp06daJLly4sXLiQqKgotmzZwpgxY1i6dGm64+vfvz+jR49m0aJF7Nu3jz59+nD27Nmbvq9v3743bIOzyJVJSpEiRXB1dXX0mlwRHR2dqnclx/Lzsw9G+esvGDeO4IAEJjKQQ5TmZV6nAGcA+0ixhQvhrrugWTNYvVpzrYiIpNfatWupWbNmiu2VV15h48aNJCUl0bx5c0JCQujfvz/+/v6ORfUmTJhAvXr1aNWqFffffz8NGjSgcuXKKVYDnjlzJl26dGHw4MFUrFiRNm3a8MsvvxAcHJzu+AYPHkyXLl3o1q0b9erVw8/PjwcffPCm7wsMDLxpG5yBzWTm5paTsdlsLFq0KMXMgHXr1iUsLIzJkyc7yqpUqULbtm0ZNWrULZ8zJiYGf39/zp07Z12vytUuXYJp02DMGDh6lPPk42OeZgKDOEZgiqp16sCLL0Lr1uBEP4sikgtdvnyZqKgox3QQeVVsbCxBQUGMHz+enj17Wh3ObXGja5/e79Ac+xV14cIFR9cbQFRUFJGRkRw+fBiwP6M+ffp0ZsyYwd69exk4cCCHDx+md+/eFkadjby9oV8/OHAAPvwQv5IFGMJ4DnIHH/MU5djvqLp5M7RrB9WqwWefQUKCdWGLiORGO3bs4Msvv+TAgQNs376dTp06AdC2bVuLI8thsuRhaAusWbPGAKm2rl27Oup8+OGHpnTp0sbDw8PUqlXLrFu3LsvOn+3zpNyqy5eN+fhjY0qXNgZMAq7mSx411YlMNW1A6dLGTJpkzMWLVgctIrnNjebKyM22b99uatWqZXx9fU3BggXN/fffb3777Terw7qtsmKelFxxu8cKTne753oSEuzdJW+9BQcPYoCltGQUw9nIPSmqFisGAwZAnz728bkiIrdKt3vyrjx9u0fSyd0devSA33+H2bOxlS/PAyzlRxqynoZE8N8o8uho+1iVUqVg+HD4918L4xYRkTxPSUpe4eYGXbrA3r3wxRdQuTIN+ZGlPMAOQnmUubhgf+wsJgZGj4YyZeyz2P71l6WRi4hIHqUkJa9xdYXHH7evVPjVVxASQii/MpeO7KMSvZiKh82+jPflyzB5Mtx5pz2/2b3b4thFRCRPUZKSV7m4QPv28Ouv9olUQkMpz36m8jQHTVkGMw5fl4sAJCXZh7WEhNifCvrlF2tDFxGRvEFJSl7n4gIPPgjbt8OSJVC7NkEcZRxDOZQczKuMpJDrWUf1b76Bu++Ge++FVas0MZyIiGQfJSliZ7PZZ3fbvBmWLoW6dSnMaUbyOoeSSjKBgQS5/TeD75o1EB5unxhuwQJITrYwdhERyZWUpEhKNhtERMBPP8HKlXDPPeQjloFM5EBiaabTk/LuUY7qW7fCI4/Y1wj69FP1rIhI3mKz2Vi8ePEtHaNJkyYMGDAgS+JJy6uvvkpoaGi2HT87KUmRtNls9sV+1q+3L/jTpAmexNOTGexNuJN5dCDU/b+RtL//Dl272u8cnTljYdwiIlno+PHjPPfcc9xxxx14enoSHBxM69at+b//+z+rQ8sTlKTIjdls0LSp/f7O+vXQrBmuJNOB+WxPCGEZLWjk8d8y4d98A2FhWnVZRHK+v/76i7CwMFavXs3YsWPZuXMny5cvp2nTpvTt29fq8PIEJSmSfg0b2m8BbdoEERHYgBasYF18Pb6lFYVc7F0oUVHQoIH98WXd/hGRnKpPnz7YbDY2b97MI488QoUKFahatSqDBg3i55//++Ps5MmTPPjgg/j4+FC+fHmWLFmS4jh79uyhZcuW5MuXj4CAADp37szJkydT1ElMTOTZZ5+lQIECFC5cmJdeeokrE8K//vrrVKtWLVV8YWFhvPLKK4B9teY6derg6+tLgQIFaNCgAYcOHUpR/7PPPqNMmTL4+/vz2GOPcf78ecc+Ywxjx47ljjvuwNvbmxo1avD1118DkJycTMmSJZkyZUqK423fvh2bzcbBgwcz+tGmX3bM158XOP3aPbfD5s3GtG7tWAToEMHmbjalWBfo0UeNiYmxOlARscr11m8JCzMmKOj2b2Fh6Yv71KlTxmazmbfffvuG9QBTsmRJM2fOHPPnn3+afv36mXz58plTp04ZY4w5evSoKVKkiBk+fLjZu3ev2b59u2nWrJlp2rSp4xiNGzc2+fLlM/379zf79u0zn3/+ufHx8TFTp041xhhz5MgR4+LiYjZv3ux4z6+//mpsNps5cOCASUhIMP7+/mbIkCFm//79Zs+ePWbWrFnm0KFDxhhjRo4cafLly2ceeughs3PnTrN+/XpTvHhx8+KLLzqO9+KLL5pKlSqZ5cuXmwMHDpiZM2caT09Ps3btWmOMMYMHDzb33HNPirYPHjzY1KtX77qfTVas3aMkJZOUpFxl+3ZjmjY1Bkwc7mYg41MkKhUqGPPrr1YHKSJWuN4XVVCQSbXY6e3YgoLSF/cvv/xiALNw4cIb1gPMSy+95Hh94cIFY7PZzLJly4wxxrz88ssmPDw8xXuOHDliAPP7778bY+xJSuXKlU1ycrKjzgsvvGAqV67seB0REWGeeeYZx+sBAwaYJk2aGGPsCRXgSCiuNXLkSOPj42NirvqLcejQoaZu3bqOmL28vMymTZtSvK9nz56mY8eOxhj7gok2m8389ddfxhhjkpKSTFBQkPnwww+v+9lkRZLiln19NJJn1KxpnzTljTfweP11JpjBNGQD3W2zOWfy88cfULcufPghdO9uH+YiInlb8eLOfV7zv1sttnT8wqpevbrj/319ffHz8yM6OhqAbdu2sWbNGvLly5fqfQcOHKBChQoA3H333SnOVa9ePcaPH09SUhKurq706tWLHj16MGHCBFxdXfniiy8YP348AIUKFaJbt240b96cZs2acf/999OhQwdKlCjhOF6ZMmXw8/NzvC5RooQjxj179nD58mWaNWuWIr74+Hhq1qwJQM2aNalUqRJffvklw4YNY926dURHR9OhQ4ebfj63QkmKZA1XV3j1VahfHzp14sGTi6luatKer9lBTS5fhp497WNvP/wQfH2tDlhErOTsg+vLly+PzWZj7969tGvX7oZ13d3dU7y22Wwk/2/yqOTkZFq3bs2YMWNSve/qJOJmWrdujaenJ4sWLcLT05O4uDgefvhhx/6ZM2fSr18/li9fzrx583jppZdYtWoVd999d7piBPj+++8JCgpKUc/T09Px/506dWLOnDkMGzaMOXPm0Lx5c4oUKZLuNmSGBs5K1goPhx07oH59ynGQTdTjGSY7ds+ebe9V2bvXwhhFRG6iUKFCNG/enA8//JDY2NhU+8+ePZuu49SqVYvdu3dTpkwZ7rzzzhSb71V/rV09EPfK6/Lly+Pq6gqAm5sbXbt2ZebMmcycOZPHHnsMHx+fFO+pWbMmw4cPZ9OmTYSEhDBnzpx0xVilShU8PT05fPhwqhiDg4Md9R5//HF27tzJtm3b+Prrr+nUqVO6jn8rlKRI1itZEtauhUGD8CKOyfRlDh3xtdn/oe/eDXfdBen89yMiYonJkyeTlJREnTp1WLBgAX/++Sd79+7l/fffp169euk6Rt++fTl9+jQdO3Zk8+bNHDx4kJUrV9KjRw+SkpIc9Y4cOcKgQYP4/fff+fLLL/nggw/o379/imM9+eSTrF69mmXLltGjRw9HeVRUFMOHD+enn37i0KFDrFy5kj/++IPKlSunK0Y/Pz+GDBnCwIEDmT17NgcOHGDHjh18+OGHzJ4921GvbNmy1K9fn549e5KYmEjbtm3Tdfxbods9kj3c3WH8ePuzyN270zFmLjXNDtrbvmaXCSE2Fjp1st/+mTgRvLysDlhEJKWyZcuyfft23nrrLQYPHsyxY8coWrQoYWFhfPTRR+k6RmBgIBs3buSFF16gefPmxMXFUbp0aVq0aIGLy3/9BF26dOHSpUvUqVMHV1dXnnvuOZ566qkUxypfvjz169fn1KlT1K1b11Hu4+PDvn37mD17NqdOnaJEiRI8++yzPP300+lu6xtvvEGxYsUYNWoUBw8epECBAtSqVYsXX3wxRb1OnTrRt29funTpgre3d7qPn1k2c2V0kGRITEwM/v7+nDt3jvz581sdjnPbv9++4nJkJBfxpi8fMovujt2hofD111CunHUhikj2uHz5MlFRUZQtWxYv/TVyS4wxVKpUiaeffppBgwZZHc5N3ejap/c7VLd7JPvdead9Argnn8SHS8ykBzPojrftEgCRkVCrFixcaG2YIiLOKjo6mgkTJvDPP//QvXv3m78hl1CSIreHtzdMm2YfOevtTXdm8YupQ0WXPwCIiYGHH4YBAyA+3tpQRUScTUBAAKNHj2bq1KkULFjQ6nBuGyUpcnt16QK//AIVK1KNXWxJDuMxvnTsfu89++z718zmLCKSpxljOHHiBI8//rjVodxWSlLk9qtWDbZsgUcfxY8LzOFxJvMMHjZ7F8rmzfb54b7/3uI4RUTEUkpSxBp+fvDllzBpEjZ3d55hCptMPcq6/AXAmTPQqhUMHw6JidaGKiK3Ts9o5D1Zcc2VpIh1bDbo2xd+/BFKlSKM7WxPDqUdixxVRo+Ge++Fo0ctjFNEMu3KTKcXL160OBK53a5c82tnu80IPYKcSXoEOYudPm0fr/L99xhgIgN43vYOif9bXqpoUfjiC7hmaQkRyQGOHTvG2bNnKVasGD4+PulaD0dyLmMMFy9eJDo6mgIFCqQ5/X96v0OVpGSSkpRskJwMY8fCiBGQnMxP3E0H1wX8nRQI2DteXnkFXn7ZvlSQiOQMxhiOHz+e7qnkJXcoUKAAxYsXTzMpVZKSzZSkZKO1a+Gxx+DffzlJYbrwGcuIcOy+/357r0qxYtaFKCIZl5SUREJCgtVhyG3g7u7uWHcoLUpSspmSlGx27Bh07Ajr1pGMjTG8wEu8STL2H/oSJWDuXGjUyOI4RUQkwzTj7E0cOXKEJk2aUKVKFapXr878+fOtDkmuVqIE/PADDB+OC4bhjGY191LcNRqw5zD33gtjxtjvEomISO6TZ3tSjh07xr///ktoaCjR0dHUqlWL33//PcXS2TeinpTb6PvvoXNnOHOGfynG47a5rDZNHbsfeMA+kW3hwhbGKCIi6aaelJsoUaIEoaGhABQrVoxChQpx+vRpa4OStD3wAOzYAXfdRQDRrDT38wqvYcPehfL99/a1f375xeI4RUQkSzltkrJ+/Xpat25NYGAgNpuNxYsXp6ozefJkx+qKYWFhbNiwIVPn2rp1K8nJyQQHB99i1JJtSpeGDRvg2WdxJZnXeJXltKCIqz2xPHzYPp3+e+9B3uwbFBHJfZw2SYmNjaVGjRpMmjQpzf3z5s1jwIABjBgxgh07dtCwYUMiIiI4fPiwo05YWBghISGptqNXzQx26tQpunTpwtSpU7O9TXKLPD3hgw/sI2bz5SOcVUQmVaOByyYAEhLsCxQ+8gicO2dtqCIicutyxJgUm83GokWLaNeunaOsbt261KpVi48++shRVrlyZdq1a8eoUaPSddy4uDiaNWtGr1696Ny5803rxsXFOV7HxMQQHBysMSlW+f13ezayaxcJuPESbzKWFxy7y5WD+fPtawCJiIhzydVjUuLj49m2bRvh4eEpysPDw9m0aVO6jmGMoVu3btx77703TVAARo0ahb+/v2PTrSGLVaxoH4TSpQvuJDKGYSyhNQVd7V0oBw5AvXrw8ce6/SMiklPlyCTl5MmTJCUlERAQkKI8ICCA48ePp+sYGzduZN68eSxevJjQ0FBCQ0PZuXPndesPHz6cc+fOObYjR47cUhskC/j4wKxZMG0aeHrSmu/YnlSDu1y3ARAXB717wxNPwIUL1oYqIiIZ52Z1ALfi2ql2jTHpXhPinnvuITkDE2x4enri6emZofjkNrDZ4MknoXZteOQRyhw4wIak+gzlHT6gHwBz5sD27fbbPyEhFscrIiLpliN7UooUKYKrq2uqXpPo6OhUvSuSR4SGwrZt8OCDeBLP+/TnK9rj5xoLwL59UKeOfT4VERHJGXJkkuLh4UFYWBirVq1KUb5q1Srq169vUVRiOX9/WLAAJkwANzfa8zXbkkKp4bYLgEuXoFs36NkTtGq8iIjzc9ok5cKFC0RGRhIZGQlAVFQUkZGRjkeMBw0axPTp05kxYwZ79+5l4MCBHD58mN69e1sYtVjOZoOBA2HdOggKojz7+SnxLnoxzVFlxgy4+277A0IiIuK8nPYR5LVr19K0adNU5V27dmXWrFmAfTK3sWPHcuzYMUJCQnj33XdpdJtWnNO0+DnAiRP2UbMrVwLwGU/Q23UaF5O8AMiXD6ZPh0cftTJIEZG8R6sgZzMlKTlEUhK89Ra8+ioYwx4q84jbN+xNLO+o0rcvjB9vnytORESyX66eJ0Uk3Vxd4ZVXYMUKKFKEKuxlS2IonW2fO6p8+CE0aGBfWVlERJyHkhTJG5o1g8hIaNAAXy4y23RmGk/i6RIP2B8MatQIDh2yNkwREfmPkhTJO4KCYM0aGDwYG/Akn/Bzch1Ku/0NwP79cM89GlArIuIslKRI3uLuDuPGwaJF4O9PKL+yIbEeFWx/APD33/bVlH/91eI4RURESYrkUe3a2e/x1KxJMH+zwdxDDeyZyYkT0KQJ/PSTpRGKiOR5SlIk7ypXDjZtgg4dKMYJ1tCEetgXqDx71j6MZfVqa0MUEcnLlKRI3ublZV/cp0cPCnKWlYRzHz8AEBsLLVvCt99aHKOISB6lJEXE1dW+knK/fuQjlu9oRWuWAPaVlB96CObOtThGEZE8SEmKCICLC0ycCCNG4EUcC3iYjswBIDERHn/cnseIiMjtoyRF5AqbDd58E8aMwZ1EPqMzvZgKgDHw1FP2tQtFROT2UJIicq3nn4fJk3ElmY95msGMc+waPNgxw76IiGQzJSkiaXnmGfj0U2wuLrzDUF7nZceu116zJytKVEREspeSFJHr6dwZ5s/H5u7Oy7zJuwxw7Hr3Xfvtn6Qk68ITEcntlKSI3MhDD9mfQfb2ZgDvMY0nsZEMwPTp0KkTJCRYHKOISC6lJEXkZpo3t6+i7OfHk3zCl3TEzZYIwLx59jzm0iWLYxQRyYWUpIikR8OG9ulnCxXiUb5isWmLpy0OgO++gwcegPPnLY5RRCSXUZIikl61a8O6dVC8OA+wlGWmBflsFwD74srNmsGZMxbHKCKSiyhJEcmIkBDYsAFKlaIpa/nB3EcB2zkAfvnFvjDhv/9aG6KISG6hJEUko+68E378ESpUoC6bWWcaUszlBAC//QaNGsHhwxbHKCKSCyhJEcmM4GBYvx6qV6c6O9mQ3IBg298A/PGHfQjLn39aHKOISA6nJEUkswIC7INR6talAn/yo6nPnbb9gL0npWFD2LnT4hhFRHIwJSkit6JQIVi1Cpo0oRRH2GDuoZrNnpn8+y80bgybN1sco4hIDqUkReRW+fnB0qXwwAMU51/WmsbUsdkzkzNn4L777A8FiYhIxihJEckK3t6wcCF06EAhzvCDuY/GrAXgwgVo0cKex4iISPopSRHJKh4eMGcO9OiBHxdYRgQt+R6Ay5ehbVuYP9/iGEVEchAlKSJZydUVpk2Dfv3w5jKLeJD2fAVAYiI89hjMmGFxjCIiOYSSFJGs5uICEyfCiBF4kMCXdKQHnwCQnAw9e8L771sboohITqAkRSQ72Gzw5pswZgyuJDONXvRnomN3//7w1ltgjHUhiog4OyUpItnp+edh8mRcMLzLQF7mdceul16CF15QoiIicj15Pkm5ePEipUuXZsiQIVaHIrnVM8/Ap59ic3HhdUYylqGOXe+8A3362G8DiYhISnk+SXnrrbeoW7eu1WFIbte5s/3RHnd3hjKOKTyNDXtmMmUKdOkCCQkWxygi4mTydJLy559/sm/fPlq2bGl1KJIXPPQQfPsteHvzNFP5nCdwJRGAL76A9u3tjyqLiIid0yYp69evp3Xr1gQGBmKz2Vi8eHGqOpMnT6Zs2bJ4eXkRFhbGhg0bMnSOIUOGMGrUqCyKWCQdmjeHFSvAz4/H+ZIFPIyHLR6Ab76B1q0hNtbiGEVEnITTJimxsbHUqFGDSZMmpbl/3rx5DBgwgBEjRrBjxw4aNmxIREQEhw8fdtQJCwsjJCQk1Xb06FG++eYbKlSoQIUKFdIVT1xcHDExMSk2kUxp2BBWr4ZChWjLEr43LfGxXQTghx8gPBzOnrU2RBERZ2AzxvmfLbDZbCxatIh27do5yurWrUutWrX46KOPHGWVK1emXbt26eodGT58OJ9//jmurq5cuHCBhIQEBg8ezCuvvJJm/VdffZXXXnstVfm5c+fInz9/xhslsmsXNGsGx4+ziXq0dFnOuWT7z1JoKKxcCUWLWhuiiEh2iImJwd/f/6bfoTkySYmPj8fHx4f58+fz4IMPOur179+fyMhI1mVwNbdZs2axa9cuxo0bd906cXFxxMXFOV7HxMQQHBysJEVuzf79cP/9cOgQOwgl3OUHTiYXBqBSJfsCyyVLWhyjiEgWS2+S4rS3e27k5MmTJCUlERAQkKI8ICCA48ePZ8s5PT09yZ8/f4pN5JbdeSds2AAVKlCTSDYkNyDI5SgA+/bZ7wwdOGBxjCIiFsmRScoVNpstxWtjTKqy9OjWrdsNe1FEslVwMKxfD9WrU4nf2ZDcgDtcogD46y97orJnj7UhiohYIUcmKUWKFMHV1TVVr0l0dHSq3hWRHCEgANasgbp1KctfbEhuQBWXvQAcOwaNGsG2bRbHKCJym+XIJMXDw4OwsDBWrVqVonzVqlXUr1/foqhEblGhQvZBKE2aEMgx1iU3JMxmz0xOnYJ777XfGRIRySucNkm5cOECkZGRREZGAhAVFUVkZKTjEeNBgwYxffp0ZsyYwd69exk4cCCHDx+md+/eFkYtcov8/GDpUnjgAYpwiv8z93KP7UcAYmL+m2ZFRCQvcNqne9auXUvTpk1TlXft2pVZs2YB9sncxo4dy7FjxwgJCeHdd9+lUaNGtyW+9I5MFsmU+Hj7VPpffcVFvHmIRaygOQDu7jB3rn0CWxGRnChXPYLsjJSkSLZLSoKnnoIZM4jDg058wQIeAcDFBWbOtK/5IyKS0+TqR5BF8gRXV5g2Dfr1w5N45vIYXZgN2FdN7toVJk+2OEYRkWykJEXEmbm4wMSJ8NJLuJHETLrTl/+WiujbF0aPti48EZHspCRFxNnZbPDGGzBmDC4YPuA5hvO2Y/fw4TB+vIXxiYhkEyUpIjnF88/D5MnYgLcZwSiGOXYNGQKzZ1sXmohIdlCSIpKTPPMMfPopuLgwjDG8xn8LYvbsCd9+a2FsIiJZTEmKSE7TuTPMnw/u7rzMGzzLB4D9YaAOHTThm4jkHpl6BPmff/5h48aNREdHk5ycnGJfv379siw4Z6ZHkMVyK1ZA27Ykx8XzBJ/zJY8D4O8P69ZBjRoWxycich3ZNk/KzJkz6d27Nx4eHhQuXDjFgn42m42DBw9mPuocREmKOIVvvoGHHiI+2ZU2LGEFLQAoXhw2boQ77rA4PhGRNGRbkhIcHEzv3r0ZPnw4Li55926RkhRxGjNnQo8exOLDffwfv3A3YE9QNm60JywiIs4k2yZzu3jxIo899lieTlBEnEr37jBmDL5c5HseoDJ7ADh4EFq0gLNnrQ1PRCSzMpxp9OzZk/nz52dHLCKSWc8/D0OGUJjTrCScUhwC4NdfoU0buHTJ4vhERDIhw7d7kpKSaNWqFZcuXaJatWq4u7un2D9hwoQsDdBZ6XaPOB1joEcPmDWL36nAPWzkJEUAe6KyYAG4uVkco4gI6f8OzfCvrLfffpsVK1ZQsWJFgFQDZ0XEIjabfa2fU6eo+O23LKMFTW1ruWDysWSJfa3CTz6xVxMRyQky3JNSsGBB3n33Xbp165ZNIeUM6kkRp3XpEjRvDhs28H/cS0uWEo8nAEOHwtixFscnInletg2c9fT0pEGDBrcUnIhkI29vWLIEatTgPlbzBZ2wYZ/P6J137JuISE6Q4SSlf//+fPDBB9kRi4hklQIFYPlyuOMOHmEBU+jt2PX88/anlkVEnF2Gx6Rs3ryZ1atX891331G1atVUA2cXLlyYZcGJyC0oXhxWroQGDXjq32mcoCgv8RYATz4JhQpB27YWxygicgMZTlIKFCjAQw89lB2xiEhWK1fO3qPSuDEvxrxNNMV4n/4kJ8Ojj9pn1m/c2OogRUTSlqEkJTExkSZNmtC8eXOKaxpLkZwhNBS+/RZbeDjvxg3kFIX5gieIi7M/mrx2LdSsaXWQIiKpZWhMipubG8888wxxcXHZFY+IZIdGjWDePFxcbMykOxEsBSAmxj4r7f79FscnIpKGDA+crVu3Ljt27MiOWEQkO7VtC9On404i82lPPTYBEB0N4eFw7JjF8YmIXCPDY1L69OnD4MGD+fvvvwkLC8PX1zfF/urVq2dZcCKSxbp3h5Mn8X3+eb6jFY1Yz25CiIqy96isW2d/MEhExBlkeDK3tBYWtNlsGGOw2WwkJSVlWXDOTJO5SY42dCiMG8c/BNKAjRyiDAD33GMfTOvjY214IpK7Zdu0+FFRUbcUmIg4gbFj4eRJgmbNYhXNaGDbxAlTlB9/tD/1s3AhXDO7gIjIbZfhnhSxU0+K5HiJifDww7BkCdupSRPbOs4bPwC6dLFP+JZGx6mIyC3LtmnxAQ4cOMBzzz3H/fffT7NmzejXrx8HDhzIdLAiYgE3N5g7Fxo2pBY7+Ma0wQP7k3uffmq/I6Q/YUTEShlOUlasWEGVKlXYvHkz1atXJyQkhF9++YWqVauyatWq7IhRRLLLVev8NGUtX9IRF+zjyiZMgDFjLI5PRPK0DN/uqVmzJs2bN2f06NEpyocNG8bKlSvZvn17lgborHS7R3KV48ehQQM4eJBpPMlTTHPsmjbNPo2+iEhWybbbPXv37qVnz56pynv06MGePXsyejhLRUVF0bRpU6pUqUK1atWIjY21OiQRa1xZ5ycggF5M522GO3Y9/TQsWmRhbCKSZ2U4SSlatCiRkZGpyiMjIylWrFhWxHTbdOvWjddff509e/awbt06PD09rQ5JxDpX1vnJn59hjGYgEwBIToaOHe3T54uI3E4ZfgS5V69ePPXUUxw8eJD69etjs9n48ccfGTNmDIMHD86OGLPF7t27cXd3p2HDhgAUKlTI4ohEnMBV6/yMixvCSYrwGV1SrPNTq5bVQYpIXpHhnpSXX36ZV155hQ8++IDGjRvTqFEjJk2axKuvvsqIESOyLLD169fTunVrAgMDsdlsLF68OFWdyZMnU7ZsWby8vAgLC2PDhg3pPv6ff/5Jvnz5aNOmDbVq1eLtt9/OsthFcrSr1vn5hJ48wHcAnD9vn5X2zz8tjk9E8owM96TYbDYGDhzIwIEDOX/+PAB+fn5ZHlhsbCw1atSge/fuPPzww6n2z5s3jwEDBjB58mQaNGjAxx9/TEREBHv27KFUqVIAhIWFpbkY4sqVK0lISGDDhg2O21QtWrTgrrvuolmzZmnGExcXl+JYMTExWdRSESd0ZZ2fHj34ig40ZwU/0pATJ+zr/GzcCIGBVgcpIrldjpjMzWazsWjRItq1a+coq1u3LrVq1eKjjz5ylFWuXJl27doxatSomx7zp59+4rXXXmP58uUAvPPOOwAMHTo0zfqvvvoqr732WqpyPd0judo778Dzz3MWfxqxnp3Y1+YKCYH166FgQYvjE5EcKdue7vn333/p3LkzgYGBuLm54erqmmK7HeLj49m2bRvh4eEpysPDw9m0aVO6jnHXXXfx77//cubMGZKTk1m/fj2VK1e+bv3hw4dz7tw5x3bkyJFbaoNIjjB0KAwZQgHOsYLmlMW+LMauXdCqFVy8aHF8IpKrZfh2T7du3Th8+DAvv/wyJUqUwGazZUdcN3Ty5EmSkpIICAhIUR4QEMDx48fTdQw3NzfefvttGjVqhDGG8PBwWrVqdd36np6eevpH8qb/rfNTYtYsVv5vnZ9oU4xNm6B9e1i8WOv8iEj2yHCS8uOPP7JhwwZCQ0OzIZyMuTZBurISc3pFREQQERGR1WGJ5C42m31Gt9OnuXPJEpab5jSxrSPG5GfpUujRA2bP1jo/IpL1MvxrJTg4GKuHsRQpUgRXV9dUvSbR0dGpeldEJAtctc5PTSJZYlrjyWUAPv8cBg/WOj8ikvUynKRMnDiRYcOG8ddff2VDOOnj4eFBWFhYqrWCVq1aRf369S2KSiSXu2qdn8asZy6POdb5mTgR0jFeXUQkQzJ8u+fRRx/l4sWLlCtXDh8fH9yvuRl9+vTpLAnswoUL7N+/3/E6KiqKyMhIChUqRKlSpRg0aBCdO3emdu3a1KtXj6lTp3L48GF69+6dJecXkTQUKGCflbZBA9od/IZp9KInMwAYMQKKFIGnnrI2RBHJPTKcpEycODEbwkht69atNG3a1PF60KBBAHTt2pVZs2bx6KOPcurUKV5//XWOHTtGSEgIS5cupXTp0rclPpE868o6Pw0a0OPfmZykCC8wFoBnnoHChSGNqY1ERDIsR8yT4oy0CrLkeZGR0LgxxMQwlLGMwz7HkIcHLFsG995rbXgi4ryybZ4UERHAsc4PXl6M5Xm6MROA+Hj7hLXbtlkbnojkfEpSRCTz/rfOj83VlWn0og3fAHDhAkREwB9/WByfiORoSlJE5Na0aQPTpuFGEnN5jEasA+DECWjWDP75x+L4RCTHUpIiIreue3cYOxZvLrOENtQgEoDDh+0LEmbRQ38iksdkKElJTEzEzc2NXbt2ZVc8IpJTDR0KQ4fiTwzLacEdHARgzx544AGIjbU4PhHJcTKUpLi5uVG6dGmSkpKyKx4RycnGjIFu3SjOv6ykGQG2fwH4+Wd45BH7oFoRkfTK8O2el156ieHDh2fZpG0ikotcWeenTRvKcZAVJhx/2znAPgdct26QnGxtiCKSc2R4npSaNWuyf/9+EhISKF26NL6+vin2b9++PUsDdFaaJ0XkBi5dgubNYcMGNnAP4bZVXDZeADz3HLz3nj2fEZG8Kb3foRmecbZdu3a3EpeI5AVX1vlp0oSGv/7IV6Y9D7KIJNz44AMoVgxeesnqIEXE2WnG2UxST4pIOhw/Dg0awMGDzKYL3Zjt2PXRR6CltkTypmzrSbli27Zt7N27F5vNRpUqVahZs2ZmDyUiudVV6/x0/fdTTlKEIYwHoE8f+zo/7dtbHKOIOK0MJynR0dE89thjrF27lgIFCmCM4dy5czRt2pS5c+dStGjR7IhTRHKqcuXso2YbN2ZwzAROUJQxDMMY6NTJnsc0bGh1kCLijDL8dM9zzz1HTEwMu3fv5vTp05w5c4Zdu3YRExNDv379siNGEcnprlrnZxTD6cEnACQkwEMPwV9/WRqdiDipDI9J8ff354cffuCuu+5KUb5582bCw8M5e/ZsVsbntDQmRSQTliyBhx4iMQke4HtW0hyAatVg40bw87M4PhG5LbJtFeTk5GTc3d1Tlbu7u5OsCRBE5EauWeenAr8DsHMndO6sOVREJKUMJyn33nsv/fv35+jRo46yf/75h4EDB3LfffdlaXAikgv9b52fgpzlW1pTgDMAfPMNvPyyxbGJiFPJcJIyadIkzp8/T5kyZShXrhx33nknZcuW5fz583zwwQfZEaOI5DZDhsDTT1OBP5nHo7hgX2rj7bdhzhyLYxMRp5HpeVJWrVrFvn37MMZQpUoV7r///qyOzalpTIrILUpIgGbNYN063uc5+vM+AJ6esH491KljcXwikm3S+x2aoSQlMTERLy8vIiMjCQkJyZJAcyolKSJZ4ORJqFMHExXF03zMNJ4CoEQJ2LIFgoIsjk9EskW2DJzVKsgikqWKFIFvv8WWLx+TeJZGrAPg2DFo1w4uXrQ2PBGxllZBFhFrVa0Kc+bgYUtkAQ9ThigAtm6FHj1AC3eI5F1aBTmTdLtHJIuNHQsvvMBOQqjPJi5gnzTljTe0GKFIbqNVkEUkZxk6FHbtotpnn/EFnWjHYgwuvPyyvbPlwQetDlBEbrcMJSmJiYkA9OjRg+Dg4GwJSETyKJsNpk6FP/6gzS/f8jYvMpzRADzxhH1G2tBQa0MUkdsrwwNnx40bp4GzIpI9vLxg8WIoWZIXGEMnPgfsA2jbtoXoaGvDE5HbK8MDZ++77z7Wrl2bDaGIiGBfFvmbb7B5ezOdJ6nDLwAcPmxfjDAuzuL4ROS2yfCYlIiICIYPH86uXbsICwtLNXC2TZs2WRaciORRtWrB7Nl4dejAYtpxF1v4h5Js3AjPPAOffGK/OyQiuVuGn+5xcbl+54vNZsszt4L0dI/IbfDqq/Daa2yjFg3ZwCV8ABg/HgYNsjY0Ecm8bF0F+XpbTktQ3n33XapWrUqVKlXo168fmVwhQESyyyuvwMMPE8Z2ZtHNUTx0KCxbZl1YInJ7ZDhJyS1OnDjBpEmT2LZtGzt37mTbtm38/PPPVoclIldzcYHZsyE0lA7M5xVeAyA5GR57DPbutTg+EclW6U5SWrZsyblz5xyv33rrLc6ePet4ferUKapUqZKlwWW3xMRELl++TEJCAgkJCRQrVszqkETkWr6+8M03UKwYI3mNh/kagJgYaN0aTp2yOD4RyTbpTlJWrFhB3FXD6seMGZNiavzExER+//33LAts/fr1tG7dmsDAQGw2G4sXL05VZ/LkyZQtWxYvLy/CwsLYsGFDuo9ftGhRhgwZQqlSpQgMDOT++++nXLlyWRa/iGShUqVg0SJcPNyZTVdqEAnAgQPQoYN9QWURyX3SnaRcO14ju8dvxMbGUqNGDSZNmpTm/nnz5jFgwABGjBjBjh07aNiwIRERERw+fNhRJywsjJCQkFTb0aNHOXPmDN999x1//fUX//zzD5s2bWL9+vXXjScuLo6YmJgUm4jcRvXrw9Sp+HKRJbShGP8CsHo1DBxocWwiki0y/Ajy7RIREUFERMR190+YMIGePXvy5JNPAjBx4kRWrFjBRx99xKhRowDYtm3bdd8/f/587rzzTgoVKgTAAw88wM8//0yjRo3SrD9q1Chee+21zDZHRLJC166waxelxo1jEQ/SlDXE48mHH9qnzn/mGasDFJGslO6eFJvNhu2aiQmufX27xMfHs23bNsLDw1OUh4eHs2nTpnQdIzg4mE2bNnH58mWSkpJYu3YtFStWvG794cOHc+7cOcd25MiRW2qDiGTS6NHQsiX1+YmPedpR/Nxz9l4VEck90t2TYoyhW7dueHp6AnD58mV69+7tmMwt7jZOA3ny5EmSkpIICAhIUR4QEMDx48fTdYy7776bli1bUrNmTVxcXLjvvvtuOBGdp6eno+0iYiFXV/jyS6hXj257ZrObqoxjKElJ0L49bN4MGl4mkjukO0np2rVritdPPPFEqjpdunS59Ygy4NqeHGNMhnp33nrrLd56662sDktEslv+/LBkCdSpw+jTw9hDFZbyAKdP25/4+flnexURydnSnaTMnDkzO+PIkCJFiuDq6pqq1yQ6OjpV74qI5FLlysHXX+MaHs6cxMepx0/spQp790LHjvYcxtXV6iBF5FbkyMncPDw8CAsLY9WqVSnKV61aRf369S2KSkRuu6ZN4f338SeGJbShEPZJU5YuhWHDLI5NRG6Z0yYpFy5cIDIyksjISACioqKIjIx0PGI8aNAgpk+fzowZM9i7dy8DBw7k8OHD9O7d28KoReS2e+YZ6NOHOznA1zyCG/ZJU8aNs09WKyI5V4YXGLxd1q5dS9OmTVOVd+3alVmzZgH2ydzGjh3LsWPHCAkJ4d13373uI8RZTQsMijiRhARo0QJWr+YjetOHjwDw8IA1a+xTrIiI80jvd6jTJinOTkmKiJM5fRrq1IEDB+jLJCbTF4BixWDLFvuktSLiHLJtFWQREadUqBB8+y3kz89EBnAv/wdAdDS0aQOxsRbHJyIZpiRFRHKPypVh7lzcXZKZT3vKsR+AX3+1T1abnGxxfCKSIUpSRCR3iYiAsWMpxBm+pTX5sa/evmABaGULkZxFSYqI5D6DBkG3blRmH3N5DBeSAHj9dfjqK4tjE5F0U5IiIrmPzQZTpkD9+kSwnLE879jVrRvcYO1REXEiSlJEJHfy9ISFCyE4mEFMoBv2WbMvXYK2beHYMYvjE5GbUpIiIrlXQAAsWYLNx4cp9KY+GwH45x9o186esIiI81KSIiK5W2gofPYZnsSzkIcoxSHAvlpyr16gmaJEnJeSFBHJ/R56CF5/nQCi+Ya2+GCfNOWLL2DsWItjE5HrUpIiInnDSy/Bo48Syq98RmdH8fDh9hWTRcT5KEkRkbzBZoMZMyAsjIdYxBu8BNhv93TqBDt3WhyfiKSiJEVE8g4fH1i8GIoXZwRv8ShzAbhwwT51/okT1oYnIikpSRGRvKVkSVi8GJunJzPoQRhbAfjrL3jkEYiPtzY8EfmPkhQRyXvq1oVPPsGHS3xDW0pwFID166FvXz3xI+IslKSISN7UqRMMG0YQR1lMOzy5DMD06fDBBxbHJiKAkhQRycveegtat6YOW5hBD0fxwIGwcqWFcYkIoCRFRPIyFxf7ZCkhITzOlwznbQCSk6FDB/j9d4vjE8njlKSISN7m52efKKVwYd7kJdqyGIBz5+xP/Jw5Y214InmZkhQRkbJlYcECXNxc+YzOhGCfNOWPP+DRRyEx0eL4RPIoJSkiIgCNG8PkyfhxgSW0oQj2SVNWrYIhQyyOTSSPUpIiInJFr17Qrx9l+YsFPIw79klT3nsPpk2zODaRPEhJiojI1caPh2bNaMQGJtPHUdynj30eFRG5fZSkiIhczc0N5s2D8uV5kk/oz0TAPi7l4YchKsra8ETyEiUpIiLXKlgQvv0W/P0ZxxDCWQHAyZP2J37On7c4PpE8QkmKiEhaKlaEefNwczHM41EqYJ80Zdcu+2S1SUkWxyeSByhJERG5nubNYcIECnCOb2lNAeyTpnz7Lbz0ksWxieQBSlJERG6kXz/o2ZMK/MlXdMAV+6Qpo0fDV19ZHJtILqckRUTkRmw2mDwZGjakGT8wgUGOXT16wN69FsYmksvliSTlwQcfpGDBgjzyyCOp9n333XdUrFiR8uXLM336dAuiExGn5+EBCxZA6dI8xwc8wWcAxMbCQw9pIK1IdskTSUq/fv349NNPU5UnJiYyaNAgVq9ezfbt2xkzZgynT5+2IEIRcXpFi8KSJdh8ffmYp6nGbwDs2wdPPgnGWByfSC6UJ5KUpk2b4ufnl6p88+bNVK1alaCgIPz8/GjZsiUrVqywIEIRyRGqV4cvvsCHSyzgYfJzDrCPTXnvPYtjE8mFLE9S1q9fT+vWrQkMDMRms7F48eJUdSZPnkzZsmXx8vIiLCyMDRs2ZMm5jx49SlBQkON1yZIl+eeff7Lk2CKSS7VtCy+9RHn28yldHMVDh8KPP1oYl0guZHmSEhsbS40aNZg0aVKa++fNm8eAAQMYMWIEO3bsoGHDhkRERHD48GFHnbCwMEJCQlJtR48eveG5TRr9szabLc26cXFxxMTEpNhEJI969VUID6ctSxjGKMA+I22HDnD8uLWhieQmblYHEBERQURExHX3T5gwgZ49e/Lkk08CMHHiRFasWMFHH33EqFH2Xw7btm3L1LmDgoJS9Jz8/fff1K1bN826o0aN4rXXXsvUeUQkl3F1hTlzICyMNw69zC/UZQ33cuwYPPoo/PADuLtbHaRIzmd5T8qNxMfHs23bNsLDw1OUh4eHs2nTpls+fp06ddi1axf//PMP58+fZ+nSpTRv3jzNusOHD+fcuXOO7ciRI7d8fhHJwQoXhq+/xs3Tjbk8RhB/A/ZFCIcPtzg2kVzCqZOUkydPkpSUREBAQIrygIAAjmegT7V58+a0b9+epUuXUrJkSbZs2QKAm5sb48ePp2nTptSsWZOhQ4dSuHDhNI/h6elJ/vz5U2wiksfVrg0ffkgxTjCf9rgTD9gXUv76a4tjE8kFLL/dkx7XjhMxxlx37EhabvTETps2bWjTpk2mYxORPK5nT/j5Z+pNn84EBvEc9vF13btDtWr2JYBEJHOcuielSJEiuLq6puo1iY6OTtW7IiJimQ8+gLAw+vIhj/MFABcu2Cd6u3DB4thEcjCnTlI8PDwICwtj1apVKcpXrVpF/fr1LYpKROQaXl6wYAG2QoWYylOEsBOAPXugVy9N9CaSWZYnKRcuXCAyMpLIyEgAoqKiiIyMdDxiPGjQIKZPn86MGTPYu3cvAwcO5PDhw/Tu3dvCqEVErlG6NHz5Jb42+0RvftinKZg7197RIiIZZzNpTRZyG61du5amTZumKu/atSuzZs0C7JO5jR07lmPHjhESEsK7775Lo0aNbnOkKcXExODv78+5c+c0iFZE/vPWW/DSSyzkQR5mIQBubrB2LTRoYG1oIs4ivd+hlicpOZWSFBFJU3IytGsH337L84zhHZ4HIDAQtm8HDacTSf93qOW3e0REchUXF/j0UyhXjrd5kSasAeDoUXjsMfvMtCKSPkpSRESyWoECsHAhbt4ezOUxArHPbL12LYwYYWlkIjmKkhQRkexQvTpMnUoA0cynPW4kADB2LCxaZHFsIjmEkhQRkezyxBPw7LPU5yfGM9hR3LUr/PGHhXGJ5BBKUkREstP48VCvHs/xAY/xJQDnz8PDD0NsrMWxiTg5JSkiItnJwwPmz8dWrBjT6EUVdgOwaxc89ZQmehO5ESUpIiLZLSgI5s0jn+tlFvAw+TgPwJw58OGHFscm4sSUpIiI3A5NmsCYMVTid2bS3VE8aBD89JN1YYk4MyUpIiK3y6BB8MgjPMICBjMOgIQEaN8eoqMtjk3ECSlJERG5XWw2mDEDKlViNMNoxDoA/vkHOnbURG8i11KSIiJyO/n52Sd6y+fNPB6lBEcBWL0aXn7Z4thEnIySFBGR261yZZgxg+L8y1d0cEz0Nno0fPONxbGJOBElKSIiVmjfHgYP5h428g5DHcVdusCff1oYl4gTUZIiImKV0aOhcWP68x4dmAdATIwmehO5QkmKiIhV3Nxg3jxsgYFM50kqsReAnTuhd29N9CaiJEVExEoBATB/Pn5ul1nIQ46J3j7/HKZMsTg2EYspSRERsVr9+vDuu1RmHzPo4Sju3x9++cXCuEQspiRFRMQZ9O0LnTrRnq8ZyATAPtHbI4/AiRMWxyZiESUpIiLOwGaDjz+GatUYwwvcwwYA/v4bHn8ckpIsjk/EAkpSREScha8vLFyIu78vX9GBAI4D8MMP8MorFscmYgElKSIizuTOO+HTTynBcb6iA67Y58p/+21YssTi2ERuMyUpIiLOpk0bGDGCRmxgDC84irt0gQMHLIxL5DZTkiIi4oxeew2aNWMQE3iE+QCcO2ef6O3iRYtjE7lNlKSIiDgjV1eYMwdbqVLMoAcV2QfAr79Cnz6a6E3yBiUpIiLOqkgRWLAAP494FvIQvlwAYPZsmDrV4thEbgMlKSIizqx2bfjwQ6qwl0/o6Sju1w+2bLEwLpHbQEmKiIize/JJ6NmTR/mK/kwEID7ePj7l5ElrQxPJTkpSRERygkmTICyMsTxPfTYCcOSIJnqT3C1PJCkPPvggBQsW5JFHHklRfuTIEZo0aUKVKlWoXr068+fPtyhCEZGb8PKCr7/Go5AfX9GBYvwLwKpV9geBRHKjPJGk9OvXj08//TRVuZubGxMnTmTPnj388MMPDBw4kNjYWAsiFBFJhzJlYM4cgmzHmMejjone3ngDvv/e2tBEskOeSFKaNm2Kn59fqvISJUoQGhoKQLFixShUqBCnT5++zdGJiGRA8+bw+us0YR2jGO4ofuIJOHjQwrhEsoHlScr69etp3bo1gYGB2Gw2Fi9enKrO5MmTKVu2LF5eXoSFhbFhw4Ysj2Pr1q0kJycTHByc5ccWEclSL74IrVoxhHE8xAIAzp61D6S9dMna0ESykuVJSmxsLDVq1GDSpElp7p83bx4DBgxgxIgR7Nixg4YNGxIREcHhw4cddcLCwggJCUm1HT16NF0xnDp1ii5dujBVEw+ISE7g4gKffYatXDlm0p0K/A5AZCT07auJ3iT3sBnjPD/ONpuNRYsW0a5dO0dZ3bp1qVWrFh999JGjrHLlyrRr145Ro0al+9hr165l0qRJfP311ynK4+LiaNasGb169aJz587XfX9cXBxxcXGO1zExMQQHB3Pu3Dny58+f7jhERLLMr79CvXrsunQHdfmFi/gC9oneevWyODaRG4iJicHf3/+m36GW96TcSHx8PNu2bSM8PDxFeXh4OJs2bbrl4xtj6NatG/fee+8NExSAUaNG4e/v79h0W0hELFejBkydSgi7mcZ/Wcmzz8LWrRbGJZJFnDpJOXnyJElJSQQEBKQoDwgI4Pjx4+k+TvPmzWnfvj1Lly6lZMmSbPnfNI0bN25k3rx5LF68mNDQUEJDQ9m5c2eaxxg+fDjnzp1zbEeOHMl8w0REssoTT0DfvjzOlzzLB4B9ordHHoFTpyyOTeQWuVkdQHrYbLYUr40xqcpuZMWKFWmW33PPPSQnJ6frGJ6ennh6eqb7nCIit82ECbBtG+N/Hsw2wviJ+hw6ZM9fvvvOvlahSE7k1D0pRYoUwdXVNVWvSXR0dKreFRGRPMvDA+bPx6NYQb6iA0WJBmD5cvscKiI5lVMnKR4eHoSFhbFq1aoU5atWraJ+/foWRSUi4oRKloR58yjpepy5PIYL9rnyX38dli2zODaRTLI8Sblw4QKRkZFERkYCEBUVRWRkpOMR40GDBjF9+nRmzJjB3r17GThwIIcPH6Z3794WRi0i4oSaNIHRo7mXNbzNi4D9ceROnSAqytrQRDLD8keQ165dS9OmTVOVd+3alVmzZgH2ydzGjh3LsWPHCAkJ4d1336VRo0a3OdKU0vv4lIjIbWUMtG+PWbCAh1jIYh4EoFYt2LjRvgSQiNXS+x1qeZKSUylJERGndf481KnDuX1Hqc1W9lMegCefhGnTLI5NhFwyT4qIiGSCnx8sXIh/vmQW8DDeXARg+nSYMcPi2EQyQEmKiEhuVLkyzJhBdXYylaccxX36wPbtFsYlkgFKUkREcqv27WHwYJ7gC/rwIQBxcfbxtffeC889B1OmwIYNoAXgxRlpTEomaUyKiOQIiYlw//3ErfuJxqzjF+6+btXixSEkBKpWTbn5+9/GeCVP0MDZbKYkRURyjH//hVq1OH40iQFMZLVnS07Epf/3VlBQ6sSlShXQrz7JLCUp2UxJiojkKJs2QePG9p4V4ARF2E3VlJsthFOmcLoPGRycMnEJCbEPhcmXL7saIbmFkpRspiRFRHKcyZOhb9/r7jZANMVSJy9U5QyF0n2aMoXOUbX0BaqWT6BqdReq1vam8t0F8PF3z4JGSG6gJCWbKUkRkRxp+3bYtw9OnLBv0dGp/3vuXIq3GOAYJdJMXmJI34AVG8mUdTlEVe8oqhY8StXip6ha+gKV7kzEO6gQFC0KxYr999/ChcEtR6yBmzsZYx9lffHif1uZMvZ1orKAkpRspiRFRHKt+Hg4eTLtBOaq5MZEn+Cff93YfaFUisRlD1U4T/p+L7qQRDkOpEp/KvIHnoV8UyYuRYva54Dx9My6zc0NbLZs/kCzWFJSyuQho1tsbPrqXZse/Pkn3HlnljRBSUo2U5IiIvI/cXEpkhrzbzRH/rzM7n2u7D7ky+5jhdl9ugR7LpYm1vim65CuJHIn+1MlLyU4hheX8eIyriTfeuw2mz1Z8fDI2uQnrS0xMf0Jwo0SjPj4dDcvCRdi8eUC+Rz/vfr/b/bfq/9/2Uo3SjarfOufOUpSsp2SFBGRjElOhsOHYfeOeHZvucju35LY9bsbew/7cik+47d23EjAkzhH0uLF5dv+OksSJSABt3QlChlJKmLx5TLeWRIfwK5lR6jaIjhLjpXe71Dd8BMRkdvCxcU+rKFMGQ8eePC/sQ3JyfZVmnfvTrnt2weXL1//eIm4k4g7sVj3OJEriSmSluslNp7EkYD7dZOKeDwta0N6+PjA5aJZk6BkhHpSMkk9KSIi2SspCQ4e/C9p2bMHzp61Jy6XL9vvMl35/2tf/+9J6zzH19f+CPiV/179/5nd5+NjTzCzknpSREQkR3N1hfLl7Vu7dhl7b1LSjZOY9CQ6mXl9ddn1EiWbLf3JQkYSCm/vrE8mrKYkRUREch1XV3sPgI+PdTEkJtoTlitJi7u7PaHw8sp5DxRZRUmKiIhINnBzs2++6XugSdKQyzqGREREJLdQkiIiIiJOSUmKiIiIOCUlKSIiIuKUlKSIiIiIU1KSIiIiIk5JSYqIiIg4JSUpIiIi4pSUpIiIiIhTUpIiIiIiTknT4mfSlcWjY2JiLI5EREQkZ7ny3Xnlu/R6lKRk0vnz5wEIDg62OBIREZGc6fz58/j7+193v83cLI2RNCUnJ3P06FH8/Pyw5cLlLGNiYggODubIkSPkz5/f6nCyldqa++SVdoLamlvl9rYaYzh//jyBgYG4uFx/5Il6UjLJxcWFkiVLWh1GtsufP3+u/AeSFrU198kr7QS1NbfKzW29UQ/KFRo4KyIiIk5JSYqIiIg4JSUpkiZPT09GjhyJp6en1aFkO7U198kr7QS1NbfKS229EQ2cFREREaeknhQRERFxSkpSRERExCkpSRERERGnpCRFREREnJKSlDxu/fr1tG7dmsDAQGw2G4sXL06xv1u3bthsthTb3XffbU2wt2DUqFHcdddd+Pn5UaxYMdq1a8fvv/+eoo4xhldffZXAwEC8vb1p0qQJu3fvtijizEtPW3PLdf3oo4+oXr26Y8KrevXqsWzZMsf+3HJNb9bO3HI90zJq1ChsNhsDBgxwlOWW63q1tNqZm69reilJyeNiY2OpUaMGkyZNum6dFi1acOzYMce2dOnS2xhh1li3bh19+/bl559/ZtWqVSQmJhIeHk5sbKyjztixY5kwYQKTJk1iy5YtFC9enGbNmjnWacop0tNWyB3XtWTJkowePZqtW7eydetW7r33Xtq2bev4wsot1/Rm7YTccT2vtWXLFqZOnUr16tVTlOeW63rF9doJufO6ZogR+R/ALFq0KEVZ165dTdu2bS2JJztFR0cbwKxbt84YY0xycrIpXry4GT16tKPO5cuXjb+/v5kyZYpVYWaJa9tqTO69rsYYU7BgQTN9+vRcfU2N+a+dxuTO63n+/HlTvnx5s2rVKtO4cWPTv39/Y0zu+7d6vXYakzuva0apJ0Vuau3atRQrVowKFSrQq1cvoqOjrQ7plp07dw6AQoUKARAVFcXx48cJDw931PH09KRx48Zs2rTJkhizyrVtvSK3XdekpCTmzp1LbGws9erVy7XX9Np2XpHbrmffvn154IEHuP/++1OU57brer12XpHbrmtGaYFBuaGIiAjat29P6dKliYqK4uWXX+bee+9l27ZtOXYmRGMMgwYN4p577iEkJASA48ePAxAQEJCibkBAAIcOHbrtMWaVtNoKueu67ty5k3r16nH58mXy5cvHokWLqFKliuMLK7dc0+u1E3LX9QSYO3cu27dvZ8uWLan25aZ/qzdqJ+S+65oZSlLkhh599FHH/4eEhFC7dm1Kly7N999/z0MPPWRhZJn37LPP8ttvv/Hjjz+m2mez2VK8NsakKstJrtfW3HRdK1asSGRkJGfPnmXBggV07dqVdevWOfbnlmt6vXZWqVIlV13PI0eO0L9/f1auXImXl9d16+X065qeduam65pZut0jGVKiRAlKly7Nn3/+aXUomfLcc8+xZMkS1qxZQ8mSJR3lxYsXB/77K+2K6OjoVH+x5RTXa2tacvJ19fDw4M4776R27dqMGjWKGjVq8N577+W6a3q9dqYlJ1/Pbdu2ER0dTVhYGG5ubri5ubFu3Tref/993NzcHNcup1/Xm7UzKSkp1Xty8nXNLCUpkiGnTp3iyJEjlChRwupQMsQYw7PPPsvChQtZvXo1ZcuWTbG/bNmyFC9enFWrVjnK4uPjWbduHfXr17/d4d6Sm7U1LTn1uqbFGENcXFyuuqZpudLOtOTk63nfffexc+dOIiMjHVvt2rXp1KkTkZGR3HHHHbniut6sna6urqnek5Ova6ZZNmRXnML58+fNjh07zI4dOwxgJkyYYHbs2GEOHTpkzp8/bwYPHmw2bdpkoqKizJo1a0y9evVMUFCQiYmJsTr0DHnmmWeMv7+/Wbt2rTl27Jhju3jxoqPO6NGjjb+/v1m4cKHZuXOn6dixoylRokSua2tuuq7Dhw8369evN1FRUea3334zL774onFxcTErV640xuSea3qjduam63k91z71kluu67WubmdeuK7poSQlj1uzZo0BUm1du3Y1Fy9eNOHh4aZo0aLG3d3dlCpVynTt2tUcPnzY6rAzLK02AmbmzJmOOsnJyWbkyJGmePHixtPT0zRq1Mjs3LnTuqAz6WZtzU3XtUePHqZ06dLGw8PDFC1a1Nx3332OBMWY3HNNb9TO3HQ9r+faJCW3XNdrXd3OvHBd08NmjDG3u/dGRERE5GY0JkVERESckpIUERERcUpKUkRERMQpKUkRERERp6QkRURERJySkhQRERFxSkpSRERExCkpSRERERGnpCRFREREnJKSFBHJNbp164bNZmP06NEpyhcvXozNZrMoKhHJLCUpIpKreHl5MWbMGM6cOWN1KCJyi5SkiEiucv/991O8eHFGjRpldSgicouUpIhIruLq6srbb7/NBx98wN9//211OCJyC5SkiEiu8+CDDxIaGsrIkSOtDkVEboGSFBHJlcaMGcPs2bPZs2eP1aGISCYpSRGRXKlRo0Y0b96cF1980epQRCST3KwOQEQku4wePZrQ0FAqVKhgdSgikgnqSRGRXKtatWp06tSJDz74wOpQRCQTlKSISK72xhtvYIyxOgwRyQSb0b9eERERcULqSRERERGnpCRFREREnJKSFBEREXFKSlJERETEKSlJEREREaekJEVERESckpIUERERcUpKUkRERMQpKUkRERERp6QkRURERJySkhQRERFxSv8PJNL4t962OWsAAAAASUVORK5CYII=",
"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 Tau 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": "914025e4",
"metadata": {
"editable": true
},
"source": [
"## Complete solver\n",
"\n",
"\n",
"A complete solver, that can use either Legendre or Chebyshev bases, chosen as a\n",
"command-line argument, can also be found [here](https://github.com/spectralDNS/shenfun/blob/master/demo/poisson1D_tau.py).\n",
"\n",
""
]
}
],
"metadata": {
"kernelspec": {
"display_name": "shenfun24",
"language": "python3",
"name": "shenfun24"
},
"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.10.0"
}
},
"nbformat": 4,
"nbformat_minor": 5
}