{
"cells": [
{
"cell_type": "markdown",
"id": "81b1d2f0",
"metadata": {},
"source": [
"---\n",
"title: \"FEM: Error Estimates\"\n",
"subtitle: \"Lecture 17\"\n",
"date: 2026-04-13\n",
"abstract-title: Overview\n",
"format:\n",
" html:\n",
" other-links:\n",
" - text: This notebook\n",
" href: L17.ipynb\n",
"---"
]
},
{
"cell_type": "markdown",
"id": "63cf0867",
"metadata": {},
"source": [
"::: {.callout-note}\n",
"\n",
"I encourage you to play around with the juptyer notebook for this lecture - you can copy the code with the ```this notebook``` button on the side of this page.\n",
"\n",
":::\n",
"\n",
"::: {.callout-warning}\n",
"\n",
"These notes are mainly a record of what we discussed and are not a substitute for attending the lectures and reading books! If anything is unclear/wrong, let me know and I will update the notes.\n",
"\n",
":::\n",
"\n",
"::: {.callout-tip}\n",
"\n",
"This chapter is mostly based on \n",
"\n",
"* @Burden Chapter 11, \n",
"* @FNC Chapter 10, \n",
"* Lecture 16: @FNC [Galerkin method](https://fncbook.com/galerkin/), \n",
"* Today: Chapter 14.4 of @Suli: Theorem 14.6. \n",
"\n",
":::"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "f8408d86",
"metadata": {},
"outputs": [],
"source": [
"#| echo: false\n",
"\n",
"using Plots, LaTeXStrings, LinearAlgebra, Interpolations, SparseArrays\n",
"\n",
"# IF YOU ARE READING THIS, THANK YOU\n",
"# FOR DOWNLOADING THE LECTURE NOTES.\n",
"# THE FIRST PERSON TO LET ME KNOW,\n",
"# WINS SOME CHOCOLATES OF YOUR CHOICE"
]
},
{
"cell_type": "markdown",
"id": "81473756",
"metadata": {},
"source": [
"## Recap\n",
"\n",
"We are considering linear boundary value problems of the following form: \n",
"\n",
"\\begin{align}\n",
" \\mathcal L[u] &:= -u'' + p u' + q u = f &&\\text{on } (a,b) \\nonumber\\\\\n",
" &u(a)=0, \\qquad u(b) = 0 .\n",
" \\tag{BVP}\n",
"\\end{align}\n",
"\n",
"where $p, q$ may be functions of $x$. By approximating the continuous problem with finite differences, we obtained a discrete problem: find $U \\in \\mathbb R^{n+1}$ (with index starting at $0$) such that \n",
"\n",
"\\begin{align}\n",
" \\begin{split}\n",
" &\\mathcal L_h [U] = F \\\\\n",
" &U_0 = U_n = 0 \n",
" \\end{split}\\tag{$\\text{BVP}_h$}\n",
"\\end{align}\n",
"\n",
"where $\\mathcal L_h$ depends on the choice of finite difference scheme. Here is the implementation we saw last time:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "31b2938a",
"metadata": {},
"outputs": [],
"source": [
"\"\"\"\n",
" fdm(n, p, q, f; a=0,b=1, A=0,B=0)\n",
"\n",
"Approximate the solution to \n",
" -u′′ + p u′ + q u = 0 on (a,b)\n",
" u(a) = A, u(b) = B\n",
"using central finite differences\n",
"\"\"\"\n",
"function fdm(n, p, q, f; a=0.,b=1., A=0.,B=0.)\n",
" # mesh size\n",
" h = (b-a)/n \n",
" # interior grid points \n",
" x = [ a + j *(b-a)/n for j ∈ 1:n-1 ] \n",
"\n",
" # Finite Difference Coefficients\n",
" # (using a second order central difference) \n",
" d = (2/h^2) * ones(n-1) + q.(x)\n",
" d₋ = (-1/h^2)*ones(n-2) - p.(x[2:end])/(2h)\n",
" d₊ = (-1/h^2)*ones(n-2) + p.(x[1:end-1])/(2h)\n",
" Lh = Tridiagonal(d₋, d, d₊)\n",
" F = f.(x)\n",
" F[1] += (1/h^2 + p(x[1])/(2h))*A\n",
" F[end] += (1/h^2 - p(x[end])/(2h))*B\n",
" U = Lh \\ F\n",
"\n",
" # add the boundary conditions\n",
" x = [a, x..., b]\n",
" U = [A, U..., B]\n",
"\n",
" return x, U\n",
"end;"
]
},
{
"cell_type": "markdown",
"id": "c75dde21",
"metadata": {},
"source": [
"Then, we reformulated the problem as: find $u$ such that\n",
"\n",
"\\begin{align}\n",
" &- (\\alpha u')' + \\beta u = f \\qquad \\text{on } (a,b) \\nonumber\\\\\n",
" &u(a) = 0, \\quad u(b) = 0. \\tag{BVP}\n",
"\\end{align}\n",
"\n",
"We assume that $\\alpha$ has a continuous derivative, $\\alpha(x)\\geq c>0$, $\\beta$ is continuous with $\\beta(x) \\geq 0$, and $f$ is square integrable: $\\int_a^b |f|^2 < \\infty$. \n",
"\n",
"We then saw that we can write down a weak formulation of the problem: find $u\\in V$ such that \n",
"\n",
"\\begin{align}\n",
" \\mathcal A(u,v) &:= \\int_a^b \\left[ \n",
" \\alpha(x) u'(x) v'(x) + \\beta(x) u(x) v(x)\n",
" \\right] \\mathrm{d}x \\nonumber\\\\\n",
" %\n",
" &= \\int_a^b f v =: (f,v)\n",
"\\end{align}\n",
"\n",
"for all $v \\in V$. Here, $V$ is a space of functions $v$ satisfying the boundary conditions $v(a) = v(b) = 0$ and for which we can define $\\mathcal A(u,v)$ for all $u,v \\in V$ (it takes a bit more work to define exactly what we mean by $V$ - see a course including Sobolev spaces - it is some space that includes smooth funcions satisfying the boundary conditions).\n",
"\n",
"This is called the weak formulation because, if $u$ solves (BVP), then $u$ solves the weak formulation.\n",
"\n",
"The idea behind the Galerkin method is to solve the problem projected to a finite dimensional space: Find $u \\in V_h$ such that $\\mathcal A(u,v) = (f,v)$ for all $v \\in V_h := \\mathrm{span}\\{\\phi_1,\\dots,\\phi_m\\}$, which leads to the following matrix equation: \n",
"\n",
"\\begin{align}\n",
" A u = F \\quad \\text{where} \\quad \n",
" &u(x) = \\sum_{j=1}^{m} u_j \\phi_j(x) \\nonumber\\\\\n",
" &A_{ij} := \\mathcal A(\\phi_i,\\phi_j) \\nonumber\\\\\n",
" %\n",
" &F_i := (f, \\phi_i)\n",
"\\end{align} \n",
"\n",
"As a result, we have taken the continuous problem (find $u\\in V$ such that $\\mathcal A(u,v) = (f,v)$ for all $v\\in V$) and approximated it by a discrete problem (find $u\\in V_h$ such that $\\mathcal A(u,v) = (f,v)$ for all $v\\in V_h$) that can be written as a linear system of equations $Au = F$. Since we are now in \"Linear Algebra World\", we can use results from Chapter 2: Iterative Methods (and last semester Chapter 5: Direct Methods) to solve this system of equations. \n",
"\n",
"Last lecture, we saw an example of *(i)* a subspace $V_h$ (piecewise linear functions), *(ii)* building the matrix $A$ (assembly) and *(iii)* solving the linear system of equations:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "fbf57a22",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAfIAAAHyCAIAAACf89uHAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nOzdd1gTSRsA8EmFQELvIEhHig0EFBRRVECxgO2s2FHs/eyi2BWs6NmwYjkrglRp0kQUG0VAeu8hCZC23x97Xy4HiKiQEJjfc889yezszpu4vk5mZ2cxCIIACIIgqLfACjsACIIgqCvBtA5BENSr4IUdAPRjDAYjNja23U16enp6enoAgJqampSUlO8dwdLSUk5ODgDQ0NCAVhMXF7e1tW1bMycnJz8/HwCgqalpYGDQFeF35Pjx48ePH7948aKrq2uXHHDHjh1Xrly5f/++vb19lxyw1+ByuVgstuMSqJdAoB4vKyvre398e/bsQeuEhYV18Kf86tUrtFpSUhJagsFgvn371rYtGxsbtML69esF8NH27NkDALhz587P7rhr1y4HB4f8/PxW5atXrwYAvHz5sosC/LHz5887ODgkJiYKrMVfcObMGTKZnJyczCuZPn26hoYGg8EQYlQCw+Vyd+/evX379vXr16enpws7nG4He+siQ0ZGZt++fa0Kra2t+d9qa2uvW7eu7b76+vr8b0kkUlNT0+3bt3fv3s1fnp2dnZCQgG7tmqC7TVpaWkREBJ1Ob1U+ffp0PT29AQMGCCySr1+/RkRErF27VmAt/oLS0lIajVZVVcUrKS4uLi4u5nA4QoxKYC5cuJCUlHT+/PmBAweSSKRDhw4JO6LuBdO6yCCTye2mbH6qqqo/rAMAGDp0aFlZ2fXr13ft2oXBYHjl/v7+CIK4urreuXPnd8MVEjs7Ozs7O2FH0eP069eP939eSVZWFplM/qnj5Ofn9+/fv2tjEwAfHx93d3cCgeDi4rJw4UJhh9PtYFrvi7BY7Lx587y8vF6/fj1y5Ei0kMvl3rp1y8DAwMrKqjNpnUaj/fXXX6mpqYWFhZKSkkpKSra2tjNmzJCVleXVKSgouHTpUmpqKp1O19LSmj59+tSpU/n/IWlrz5494uLiO3bs4C9MS0u7d+/e2LFjx40b19TUtH///oyMDADAqVOnFBQU0E+EdsGCgoISExMXLlzI/wOlqanp2rVr4eHh1dXVsrKy9vb2y5cv589oOTk5V65csbS0nDx58uXLl4ODg+vr6/X19desWTNkyJAOovXy8oqPjwcA3LhxA30BANi+fbuMjAz6OjY29tatW1+/fsVisWZmZkuWLBk0aNAPv9suhyZ0DQ0N/hL+t510586dnTt3dmVk3a+0tDQ3N9fa2rp///4PHjwQdjgCIexRIOjH0LF1DQ2NDuqgY+sjRozo+FDo2PrIkSO/ffuGwWCWLFnC2xQSEgIAOHLkyJkzZ8CPxtZramrQC6oyMjJDhgwxMTGhUCgAgMDAQF6dFy9ekEgkNPKBAwcSCAQAwJQpU1paWnh12o6tS0pKqqmptWoO/Wdm7969CILU1dXJysqiR5OSkpKVlZWVlVVUVERrth1bLy4uRsdkKBTK0KFD0YSrpaX19etXXp3w8HAAwMKFC9ELrRoaGtLS0gAAMTGx2NjYDr4HHR0dMTExAACZTJb9v8LCQnTrpk2bAABYLNbY2Bi9so3FYs+cOdPBAbvJ+/fvJSQk+Et8fHycnZ1/9jje3t5dF5SAPH36FABQWVkp7EAEB14HFxkcDif/v0pKSlrVYTKZJW2UlZW1PZq2tvbIkSPv379Po9HQEn9/fywWO3fu3M4E4+fn9/Xr1/Xr11dWVr579+7z5891dXWhoaG8PnJRUdGcOXNYLNatW7eKioo+fPiQl5dnYWHx7NkzLy+v3/gagIyMTG1t7YQJEwAAiYmJtbW1tbW1lZWV36s/d+7cjIyMRYsWVVRUpKamVlRUrF27tqCgYMaMGa1GlgMCAlpaWvLy8oqKiurq6vbu3dvS0rJly5YOgsnNzfXw8AAA3L17t/b/0K7x3bt3T548qamp+eHDhy9fvmRnZ7969YpCoaxfv57Xr+8S5eXlfn5+vIvhqMrKyvr6et7btn3zX+ut/4KUlJQLFy60Oldzc3MFMKxfUlKSmpr68uVLOTm5wsLC1NRUNpvd3Y32BDCti4yysjLt/2o7iPz27VuNNoyMjNo94MKFC2k02pMnTwAADQ0Nz549mzBhQif/qufk5AAA5s+fj/aaAQA4HG78+PGGhobo2wsXLlCp1AULFsybNw8tUVdX9/f3x2AwZ86caXups5skJibGxMSoqan5+fmhPx2IRKKPj4+RkdGHDx9evnzJX5lIJD58+BAdO8ZgMLt379bQ0EhJSWEymb/Q9JEjRwAAPj4+pqamaIm9vf2uXbu4XO7Ro0d/83PxlJSUnD9/fvbs2a6uru/fv+eVjx8//sSJE7y38vLyraar9uvXT11dvavC+J7AwMDs7GxTU9ORI0dyuVy0sKCgwMDA4O3bt93dem5ubmpqakxMjLa2dmpq6rt37/rIJWI4ti4yJCUl58+fz1+ipKTUqo6ysvKkSZNaFUpISLR7wJkzZ65bt+7GjRvz588PCAhoampyd3fvZDBor3zv3r0nT55sd3p7ZGQkAGDp0qX8hSYmJiNGjIiPj3/z5o1g5pWjYSxYsAAdKkFhsdjFixdv3bo1IiKC/+saPny4mpoa7y0OhzM1NS0uLi4qKtLV1f2pdisrKz99+qSgoDB16lT+8qVLl27ZsiUqKorD4eBwuF/8VHzOnj27Z88eJpNZWVn57ds39EpARUXFx48fd+3axV9z1apV/G/19fV5ebabNDc3x8XFHTt27M6dO8XFxc3NzeipGBERgX637e4VGBiYnp7+w4OLiYmtWbOm4+9w1KhRo0aN8vLymj59+vLly3/tU4gimNZFhqysrJ+fX8d1dHV1r1y50skDkslkNze3mzdv5uXl+fv7S0tLu7i4dHJfT0/Pe/fuvXjx4sWLFzo6Ovb29i4uLhMnTsTj/zmjCgoKAADogDI/PT29+Ph4dKsAFBYWAgDaJmU0MHQrj6amZqtq6Ah7Q0PDz7aLfkAdHZ1W9/vIyMgoKChUV1fX1tYqKir+7GFbYbFYEhISEhISgYGBAIDRo0ej5TExMQCAUaNG8Vd2cnLifysrK9tqdmyXi42NRc+o4OBga2trXvciJibGwsJCUlKy3b1cXFw6fx7+UFVVVUlJyeDBg7vqgCIBpvU+zd3d/caNG7t27UpOTl65ciU6TNEZ0tLSb9++DQgICAwMjI6Ovnr16tWrVwcNGhQYGIiOLKO/dsXFxVvtiDYhsCHOjsNgsVj8hV3SfeZvl/8nQqum2/0GamtrR44c2dzc/MPjz5w58/DhwwQCAb3mHBAQYGdnJy8vj26NiYkxNTVt+2Pup1RUVPj7+7cqjI2Nbfst2dvbW1patj3C+PHjAQAMBuP58+cHDhzglUdHR/OG5robOjDV8XSm3gem9T7Nzs5OR0fn7t27AIDOj8CgiETiwoULFy5cyOVy3759u2fPntDQ0K1btwYEBAAAlJSUampqCgsLTUxM+PdCu7HKysrfO6yYmFjbvFZRUfFTsfGgqa1Vr7wzYfwmtN2ioqJW5SwWq7S0FI/H81IwPzk5uS9fvvxsWw0NDSEhIb6+vryS6OhoBweHnz1OK8rKytu2bWtVyOFw2hZ2LCgoiE6nu7m5oW+/fftWVFQ0ZsyY3wyvk9LS0qSkpLS1tQXTXA8B03qfhsFgVq5c6efnp6Wl1W6HqzOwWKylpeWdO3cUFBR418Gsra0zMjKeP3/On9Zramri4uJwOFwHbWlqaqalpTU0NKADIKi4uLhW1dD5lD+8GxYdZ3j+/HmrifDopLfhw4d35gN2DJ3/3ioSbW1tZWXlgoKCDx8+8E9UDwwM5HA4lpaWRCLx95tGpaamtrS0jBs3Dn1bWVmZkZHh7e0NAPj48ePAgQO7qqFfk5CQMGDAAN6dUNHR0UQi0cbGhkajVVZW6ujotKofFBT0+fPnHx5WXFx89erVP/yBlZaWNmTIkI5vleh9YFrv6zZv3rx58+af3evFixdWVlb8o8MfPnwAfP3f5cuX+/v7+/j4zJ49G+0rIQjy559/0mi06dOnd9BNNjU1TUtL8/X13bt3L1py//59NAvzQ9NEcnKyubl5B3E6OTlpamomJyffuXOHN3czMDAwNDRUTk5u5syZP/vB20JH5JOTk/mPhsFgVqxY4eXltWXLlsDAQHQ0pq6uDl2tAZ0T2VUKCgpIJBKvQxodHY0gyIgRI2g0WnR0tNDTen5+Pv9aDlFRUUOGDCGRSI8fP271Sw41ceLEiRMndlXr79+/R+fC9ilwgmOvkpycLNueixcvdm1Dp0+f1tDQmDBhwurVq/fs2TN37lwXFxc8Hs+7BdHa2nr79u1VVVVDhgxZtmzZli1bhg0bdvnyZS0tLf7hgrY2bNhAIBD27ds3bty4VatW2dnZzZ07t+1QrJubGw6HW716tbGxsYWFhZWVVbtHIxKJ/v7+JBJp/vz5Li4u27dvd3NzmzZtGh6Pv3r1Ktrl/03Ozs5kMvnUqVP6+voWFhYWFhbojQLbt2+3srIKDw83MzPbsGHDmjVrTExM0tPTp02b9rPjXR0zNjbmcDgtLS0AgIaGhsuXL1MoFCUlpb///rsnZDRjY2PeT5nExMTXr1+jV7DfvHnDmw7bTRgMRnZ2dl8bWAewty4SSCSSlZVVxwPBMjIy30ttAADeDf2SkpLm5uYd/3VSUlIyNzfveAL70qVLpaSkkpKSIiIiuFyuoqLiyJEjd+7cyVuKAABw6NAhExOTU6dOXb16FUEQeXn5FStW7N+/n/+D9OvXz8rKCl0AADV06NCgoKAtW7ZERERERkZaWFiEhoayWKz09HRVVVVeNUtLy+joaPReJ/57kbS1ta2srHj37gMA7O3tExIS9u7dGx4e/uLFCyKR6ODgsG/fPv55IFJSUubm5m1nwmhra5ubm39vhihKQ0MjMTHxr7/+ysvLKy8vR/7/uDESiRQZGXn48OGbN2+i/5IZGBhs2rRp3bp1XTsmYGVltXPnzuXLlxsbG7NYrICAgEWLFm3ZskVTU7O782ZnbNmyZdmyZdu2bSOTyerq6oGBgZ6enjt27HB0dOzupj99+sThcPjPyT4Cg8CH3kG/Ae0ndpz42Gx2c3Pzz64q1dzcjMViu3AMGgBQX1/Pn/EFpqmpCYvFtjsxRrR4e3uL0JowFy9e9PHx6WBd694KDsJAvwWHw3Wc0wEAeDz+Z3M6AEBcXLxrczoAQCg5HQBAIpF6QU4HAEyePFnYIfxYTEwMOu8zMTFRYDMpexTYW4cgqFexsrLKz8/PyckZNmxYbGzsb87fF0UwrUMQ1Kv4+fl9+/atubnZzc2Nd+dtnwLTOgRBvQ2LxcLj8X1tujoPTOsQBEG9CrxkCkEQ1KvAtA5BENSrwLQOQRDUq8C0DkEQ1KvAtA5BENSrwLQOQRDUq8C0DkEQ1KvAtA5BENSrwLQOQRDUq8C0DkEQ1KvAtA5BENSrwLQOQRDUq/SJtF5TU8NgMIQdBQRBkCD0ibS+f//+O3fudFyHTqcLJhgRxWKxmEymsKPo0eAp1DEOh9Pc3CzsKHq0rjqF+kRa53A4HA6n4zpcLlcwwYguuIZzx+Ap9EPwFOpYV51CvSSte3p6YjCY5ORkYQcCQb0ZlUp9+vRpQkKCsAOBOoIXdgBdICEhISsri0KhCDsQCOrNqqqqBpkPxquKNVXSp09y9TvrJ+yIoPaJfG+dwWCsXLny4sWLffYBVxAkGI8eP8Lrk7TXmhrutvC/fgMOlPdYIt9b371797x58/T09Dqo09LSUlhYmJqair4lEolmZmbt1lywYEFZWVnXRyn60FHRrv23E4PBXL16tV+/fl14TKjLtXCYqeVp8cVv7n+510xtAghgM1hshL029E9LpSFTB03kNLKlpKRIJJKwI4X+IdppPTk5OTo6OikpqeNqubm5YWFhISEh6Fsikfj8+XMJCQn+OnQ6HYPBPH/+/NatW/AEFYx169bl5OTIysoKO5CugZ5Cwo6ia9BotMyCrAYK/VbsvRpSg5GS/vtnKYwP9f1wKm/WRnCa2V77vULuh3/TywksCa8trUZyWZe2np01YYaRkdGVK1d8fX0tLCxGjx7Nf0wOh8NkMn84eaEv68wpJCYmRiAQOq4j2ml98+bN1tbW169fBwAwmcynT59KSEi07YkbGxvPmjXLw8Ojg0MhCEImkwEAo0ePhsP0giEjI0MikdCvvRfgnUKi6+HDh4HxwWPcx9+IusOVxtjp2ZhJGQ2SMZnqNIVqQ5WUlMRisdnZ2fLy8vLy8n+C7QAALsL9Up0ZVRC/582xwd62A6WNC1qKSRKkjx8/2tjYGBgYzJ49+8SJE6mpqYMGDSISibDP1IGuOoVEO607ODiUlpaioytcLjczM7OiouJ7AywQBLWSn5+vpqZ29PjRK8+vLz+yOoH5BjMMU9Vc4z1zzzDVoQTsv/lBSkoKfWFgYMB/BCwGa6ZobKZovNZiWV5DYXRB/Ik3F5gGTBsNq2J2We63bw319VVVVVu3bh02bNjKlSuTk5OnTp0Ke07dSrTT+t69e3mv7927t337disrKyHGA0ECExMTs2b50ura2qXLV3h5H/qpfd++fauqqnrtzvXbUfcmeU7L0Mseu2uSmrzqqUEHDeU6ukzVMW1pTe2BmosG/pHXUJhYknLp/Y0iaqml2lB7LZvo2BgMAjIyMp49e8bhcExMTEpLSydNmoTD4X65Oeh7RDut83N3d1dSUhJ2FBAkCAiCzJ7uesJWW1dObfnNqyOsLMePH48Vl+h4r/T0dC6Xm5KVeu75pYEuQyv6V7tumzVCw3Kn+kY5Ulde4dCW1tSW1pxj7FZOr3xdlHz3y6PDiaetVM0tlYfcvndbHC+ekJBw7tw5RUVFLBarqamppqbWha1DvSetnz59WtghQJCA0Gg0BqNplJYCBgALBdK7C14Dk25zGTQMDo8RE8fgiRiCGIYohsETMOISbA6XzsXkiTMflX6p0ZVGiHg3F0srCc2hkmOIBHEMVQxD+9wsLoHB4jAkSQwWhxEnYXAEDFEcQyBiCMTfiVNFUmm6kct0I5f6lobEopSIghjfd5cGKZmM1rR5GvxMkiBx8uRJV1fX7OxsNpstLS3d7kGYTGZ5ebm6ujrs2ndS70nrENR3EMtzTRUpW15l6EqTggvqY+6FqhkaopsQFhNhMRFWC8JmFpbl57YUH3l0Hmsoo0VRtuw/2lJcQxcjBdgshMVEqPUsVgvCZiGslv/vxURYTKSpEWGx0CMgLc0Ih43mdwyeiJUgYwhiAE/8p+Sf/8QweAKGIIYhELEk8r+FBAKvsiSBOEbKaJTpAJY4LrkmPaEk5VzqVX05nRETLb94povhxQYMGDBr1qzjx4+3+qSJiYmukycRMUBMkhwRE6epqSnwL1v0wLQOtcbhcDrfLeJyuRgMptdM7BMBXC415Bb9TfjjJ0+uBL2qrqwI8nU3/H9OBwBgCMTylrq3VWkXg69QJRnmmoNXLtg0WstGgST3aw3+k+uZzQiHhTQzEC4XaaIhXC7SzEA4bKSl6Z9/GJgtCIvJphUjXC63iQa4CLeJBricfyowWxBWC5fVApgtJhy2mQS5BYtNl/uWmJt1PeWqBoe4b9eEwU0tmcfX5xQWmw+zlJCVBxgsliS5bf3ufZaaE/WVT6fkHfU+cP7S5S76HnszmNZ7m5CQEEtLSzm5X/w7XFVVpaamxmKxvlehqKgoLy9v1KhR6Ntp06a5ubktWLDg15qDfgqHWlt78wgGT1DefK6cyiipKC2vrKitrQUAcBEkuy43oTglPCuquL5kiOLAJaMWjDOyp4j/7oS5f4ZiJH73OOi8dXSCI7eZAbjcfk2NE7jcZkbDl/pvCVUf99d8omDFFPHiRMDWa2LJyUhzm+g0Ol2VIgsAUJUkpjU0/GYMfQRM673NunXrbt682X0zghISEq5du8ZL69u2bdPQ0OimtiB+LTkfa28fk7R2lJowF2Aw42xtGP05YqoSbnOmL/XxzJMoYdOYmoja7AFTZFQoo0aOEna8HUGv7mIlyAAAMlC30jK2ApPW/X8K/KWihMbqhoa4mit7/GYt95zlvZ/ZxCKK4x8/PCbswEUDTOtClpiY+PfffzMYDDc3NwcHBw6Hc+TIkSlTppiamgIAHj9+zGQyZ8+enZube/fu3dzcXCUlpaVLl/LmDpeVlf311195eXmamporV66Mi4urrq4+f/78kydPXF1dLS0teQ3V19dfuHAhIyNDWlp66tSpDg4OAICsrKzr16/X1dWNHTt25syZrWI7ePCgh4eHgoICACAwMFBCQsLc3DwgICA7O3v79u1kMnnXrl0FBQWSkpLoiOfbt28DAgKampomT57s6OgIACgoKHj06NHIkSMvX75MJpPXr18Px0Z/BYI0vnpIi3kiO2ezuJE5AKCxsTHvW/6wLWMBAPQ8am5clvd2rxP7j3tuXWJsbCzscH8Rbwr8avOlGTVZQQZh+2OOx94PVZ+ip+rQv+DW59DH95ycnIQdpggQ+aW+RNrjx48XLFgwdOhQBweHtWvXPnz4EIfD6evru7q6UqnU1NTU5cuXDx06FADw/v17BQWFOXPmqKmp2dralpeXAwAKCwsHDx7MZrNnz56tpKRUWFiora0tJiZmZGRkbm6OpmOeJUuWFBQUuLu729vbFxcXAwCys7NHjBihqKjo6Oh44MCBAwcOtArv3Llz6A98AEBkZOTr16+JRKK2tjaFQjE3Nx84cCAA4N69ex8+fAAAxMfHOzk5GRoajho1ysPDA731t7i4+ODBg0ePHnV2dmaxWM7OznDF7Z/FpTVUX9rVnPFWafN5NKcDAIqZZVgJXMXr4sa8+tp35bScelNtE39/f9HN6fywGIyJgtFWu7X+jmdZFc0yJgpYApY8QD7rfSKA509nIH3AqlWr/Pz8Oq5DpVIRBJGWlkZfIAgSHR2to6Nz7ty5yMhIHR2dCxcuhIWF6ejoXLp0KSQkREdH5/Lly0FBQTo6OteuXQsMDNTR0blx48bTp091dHRu3br16NEjHR2dgICABw8e6Ojo3L9/PyAgwM3Njb/RAQMGREREoK/DwsJGjBiBvl66dOmMGTN0dXUfPHjAX7+qqio3N3fKlClXrlxBEGT16tUeHh6tPoiBgUFSUlLbD2hqanr//n3+Ek9PT09PT/T127dvJSUl2Wx2ZWUlHo9HC5WVlbOystDX69at27dvH4Ig9+7dGz9+PO8gkydPvnHjBoIgU6ZMOXToEFr49OlTPT09BEFev35NJpNpNBqCIGw2m0wmFxcX8/YdMWJEYmLid/5ARA/vzOlCzdkfS/fMqX9+FeFw0BIWh33ny99T/p634+JeipxUP13Nkz6nurzd7sBmsxkMxs/utddrn5KBivZMI3FFiRPu4xgfE7ojth6iq04hOAjzXTY2NuHh4crKykQikf+FiooKHo/nf6GqqorD4dAXWCzWzMwMfTF48GD0Pgtzc3P0xdixY3nHZ7FYX79+XbJkCTrthMVi8S5U+vj4qKmpjR07dsaMGWhJWFjYqlWrVFRUSCRSVlaWhYUFACAjI2P+/Pmd/DgHDx5cuXLlrl27nJycNm3apKmpmZ2dzRt4GTx4cHNzc0lJyS8v2ZGdnc1bdcfc3Pzbt2/ox9HQ0JCUlAQA4HA4WVnZ2tpadXX1X2uib0EQWuyzxsgHsnM28TrpefUFhxJ98UzcbKKLoaW+/X1bdDCtF9u7a4+JkXFcXFz+rJJ3dlIlgX/pmVgCLJzA3hGY1r8Lj8fr6Oigr7vqhbi4OO/4BAKBTCaHhIQYGRm1avrkyZMDBgx4/fr1ly9fTExMAABr1669fPmyvb09AOCPP/5AH44lKytbU1PTyY8zZcoUFxeXd+/e+fn5OTk5ffnyRUZGpr6+Ht3a0NDA4XBkZWX5F9EmEAi8f2moVCq61CIGg0Ha+yEsIyPT8P+JCnV1dRQKBV1nDouFA30/jUun1t4+hrBalDadxUnLAwA4COd+xtMHGU/nGk7f5bbV4dDwIUOGCDtMQcBgMDNmzED7Nxse7zgs/n53yIP+zn8IO64eDf6VE6bp06cfOHAAffQzm81OT08HAMTGxl66dOnJkydHjhyZNWsWg8EAAKD/BwBkZWW9ePECfe3m5nbx4sWKigoAAJPJrKurAwAoKCgUFBS0bevjx49YLNbCwsLT0xPdxdnZ+fr161QqFQDg6+s7cuTIVgsw6evrR0VFAQAKCwt5jSoqKpaUlLDZ7FbHd3Z2vnjxYktLC4Igp0+fdnZ27qIvqc9pyf1UcXwVXlFdYeVhNKfnNxStCt36IjXk7Y5op/5js7Ky3NzchB2mEByfeoCrpvE853Hi6zhhx9KjwbQuTCdPnsThcP379x88eLCWltbff/9dX1/v6el5/fp1NTW1JUuWDBs2bN++fQAAb29vNze3wYMHL1q0aNq0aehQyezZs+fPnz9w4MBBgwbp6+t/+vQJAPDnn3/u3r1bV1f33r17/G2tXLlSU1PTwsJi2rRp586dAwDMnz9/3LhxBgYG+vr6ISEhV69eBQBgMBje0qDodVQjI6MFCxZMnToVbdTW1tbY2FhPTw+dZkMikdBe+caNG/v166erq6ujo5Ofn3/q1CkAAB6P5638BwCQlpaG9393BEFoMU9rbxyWnbVOxnUlBodHABKYE+r5cqsxVm+F1rzXoXGivvbv78BjcaennXotC/xv7YKLhXSg/R/UvYynp6eZmVnH6603NjZSKBQZGZmioiIBrxrKZrNramrQZY86qMZkMuvr69suZ8blcqurq+Xk5PD4HwypMRgMKpWqrKzMf1Moi8ViMBjfW44DPXjn11BraWlhMpmd/AJtbGxOnjxpbW3dyYP3cOgp9Mu7c+nU2jvHkeYmuYV/op30MlrFkcTTNAYt8fCryycv2djYdF2wQsB/O9LvyMx+vTXx+JZB6zGNBOiMrHQAACAASURBVAsLC/6BTVH3m6cQDxxbFz48Hq+srPzDakQisd30isViO5l2JSQkWj0TCgBAIBC+l9N/6uAoMTExMTGxzteHUMzCrNobh8VNrKSnLEM76S9yws4lXaHGV8VfiMJOgWsz/MtI33bFm8Az6X/JJkmdPHnyyZMnwo6ox4FpHYKEijfjZfYGceNhAIAKeuW+qGNMhOXMsXNYbY+F6+20McFp49dbHjmO6vs1V3/58kVRUREuys0PpvXf0tjYGBMTk5WVJSsra2lpid4aCkGdxG2m1wX4cOoqFdf74OWUAQDRhfEnE8/nPE0/vfjY+IXjhR1gD4WXU16oMe5Qw/uopgSJ99hTp059/vwZXrbhgWn91wUEBGzatKmsrAx9i8FgZs+e7evrK+COw6tXr+rr65WVlTMzM//444+2wyxQz8Qs+lrrf1jcxFJuwXYMDl/bVLc38sjX0pwzrkeVxsv3mid3dxOZCXNXHovxIr2fbj/trsXdmpoaLperoqIi7Lh6BDgT5hcFBATMnTuXl9MBAAiCBAQETJgwAZ2wKDBv3rzx8/M7efKkvr4+zOmiAUFoMU+rL+2WnrwEnfESXRi/OHhdWXrJbHEXQwU9mNN/CEuWVrR23krvd/XDbYw6ISkpafjw4by1Lvo42Fv/FRwOZ/PmzQiCUMTw+0cZje6vUMVgHkvIjsyrSktLu3HjxrJlywQWjJqaWnh4uMCag34Tt5lRF+DDqa1Q2nAaL69S19xwOM4nKTPFmTR6+9Ytwo5OlJDHTJc9uGT7H4sPvD7h53ji6tWrBAKBTqejdzX3ZbC3/ivS0tJKS0sBAJus9WYYqytKiBkrUC5NHCxPIgIAgoODO3mcly9f7ty5MyUlJSoq6siRIykpKb8QDIfDCQ0NDQ4Ojo2N/YXdIUFiFmVXnvDEScsprj+Fl1eJKng976mHCklpufwf25ZsFnZ0IgYrLkkZM10/OXnGgCk7og+OGGXz+vVrS0vLhj6/LDvsrX9Xy9f3Nf7e7W5KzylHXwxQ+HeSqRgOqysnWVPCzH8TW7pjers7EjUNFTz+OSaHw2lubqbRaJmZmfPnz6fT6cHBwcOGDeNVzsnJyc/Pb3UEDAZja2vLP4lQX1/f3NycRCLNmjVLQkICXS4G6oHoCcENwTdkp68mDR5JY9Ivvr0U+jmSHlS17uEKnC283PcryCMnl8c9dx0zo0Cu6HCi7z6nraWlpc3NzR3M2e0LYFr/LjH9wSq7/NvdpPsmBTwcAwDIq2eM0PjnOURsLlJQzwAAqA0c9r0dseL/3ouBw+GmTp16/Pjxo0ePAgCSkpL4czoAQE9PT09P74dx6uvro7d4mJqa3r17F6b1HojbzKi778upLlPa4IuXV31T+u5AzAktrPoVRx/1hfDJy78OQyBKTZjT8OLqhpXea8N33Et/smTJkujo6MWLFz958oRI/K3na4sumNa/D4PBfudBX8NsbOXk5Gpra08m5mhJk0ZoyNU3sw7Ff62gtwAAxjs6fW/HVurq6iQkJNDb5EJDQ3fu3Jmfn9+/f390a7u9dQDAyJEjeb312tpac3PzoqIiDAaDxWL77Hnck7GKc2r8vcUMhsit92Fwmb7J59+Uviu5nb13/0YtdS1hRyfyJC3HN0Y95uZ89rbb6RGySUu636hRox4/flxYWNiZXlHv1CXL+/Zwv7beesfQZVVQYrh/L1Ho6urS6fROBvbkyRN0EfOqqipXV9e4uLjCwsJO7ovicrn+/v4IgrBYLEdHx+zs7J/aXbj6wnrr9Dfhpbtm099FIwjytixt2sOFdrscC0p/7k+5d/i19dY7g/E+tuLEaoTL/VyVMfnhvG/1BQiCJCUl+fj4dEdz3Qeuty5knp6eTU1N+/bto9PpLRwuWjhq1Khr1651fpahlpYWuryqgoLCokWLZGRk+vXr91NhYDCYiRMn3r9/n06nnz17tu92T3oepKWp7r4vq7JEcd1JjqzCpfc3QvOirNmDFLSlNVV/7k8Z6hhpkG1j1KOmD69NBo9cNXTR7tjDFyec0NPTc3d3d3Fx0dXVFXaAggbT+q/bvHnzH3/8ERgYmJGRoaCgYGlpOX78eMzP3OrNv2T2pEmTfi0MBQWFWbNm/dq+UDdhFefW+HsTtQcorTv1pf7b4SCvhm+1ku+wW29uEHZovREGI+W8oP7v86SBIybojMmqzd33+tgx+70ZGRk5OTmfPn0yMzMTdogCBdP6b1FXV+94YUioD2KkRDQ8vyo9bQV+0Ii/Pga8zInwMHMvry2aebH1Q8ChriJuOBQnq0h/Ey5p7bjafMm2aK/LabdWDFlYXl4+a9aszMxMAS/LKlxw3joEdRmkpan25pHG6CeK604WaKktCV7/pSgjeUukHJXi7u4O7wHuVtIui6kvbyHMFiwGu9tmU0xRQkR+jK2t7bNnzwgEAofDEXaAggN767+lvLz8xYsXWVlZ8vLyw4YN439UKdTXcCqLKx/6Evrpya05/uDby78znptS9Ta7rdsQtUJbW1vY0fV+xH4GRC0j2utAypjpUkTKYbtd6yJ2alDULCwsjhw58vnz59u3bws7RgER+d761q1bx4wZY2VlNW/evM+fPwuyaV9fXwMDg2XLlp04ceLPP/90cHAYO3ZsXl6eIGMAADCZzICAgLdv36Jv37179+TJE39//8rKSgFH0pe99Duxf9GMKLxGnfOMla92fKj4Qo7AFMV+k5KSgjldYKQnLWp89ZDLoAEAtKT7bbby3Bl7qLqpdv369f3790f6wCODUCKf1tXU1Pbt2+fn56elpWVvb9/Y2CiYdv38/DZs2IA2h8H9c5n01atX48aN4z13VACoVOrdu3efPXtWXFwMAKiurj59+vS0adPmzp27cePGvnMeCxHCYl7f4D5n776nypxVF89MXuuqy+w3oEDzL5+L169fF3Z0fQteSYNkYk2LfoS+tdWwctEbvzfuCI6AO3jw4MGDB1NTU4UboWCIfFpfv379qFGjhg4devDgQRqNlpOTI4BGWSzW7t27AQBEaTGzrda21yZanx+vPLIfACA3N/fKlSsCiAElJSXl7u7OW4/0+fPnJiYmAAACgcBkMj9+/CiwSPomVkVh5am1D5JTVWcZaU7V115k3PKu8d7em6ampr3pYWwiRMp5AS0+iEOtQ98uNJstJy7r+/YSAMDU1BT9a9vriXxaBwCUlZVlZGQcPnxYX18fTWrd7d27dzU1NQAALTdD2YGKGByGKC1msHQQUVYcABAWFtbJ4zx9+nTt2rWRkZFRUVH79u2LiYn5zcAyMjJ4jzCmUCgZGRm/eUCoA4yUiKozm8l2U5WGmFO/1iIcpD69RllR+cOHD3AJB2HBSctLDnNoDA9A32IAZseIDenVWU++Bk2bNg1dEa/X/4rtDZdMvby8oqKiysvLL1682O7d858+fbp//z669AoAQExMLDY2ttW0BDqd3mrKeWr5h71xR9ttsfJdKfqCpPLvEqAYHIakJMGsa47PSpr0cE67OxrJ658Ysx99zeFwiEQim81uaGhwdXWlUqkxMTF2dna8yp1ZPKAVGo1GIBB4H5NKpbZbrYdAEKSpqYlGowk7kJ+GsJhNQdfYBRmSi/fFcysqhzFk31MSlryUV5Y/F3hLFD+RAKCPqBbAjBTsCBeq71rssAlYuX8eEbx72MZNcXuViYoDFYy3b98+Z84cd3f37g7jF7TNQm2JiYnx/o5/T29I635+fgCA9PT04cOH6+npte0oDRgwYPz48XPm/JNq233WM4IgvH4uylxl0IsZd9ttMVUn1eKIBQCAXkSVMVZAC7ksLqOMBgCwM7F9+p0d+eFwOGdn5/37958+fRoAEBcX5+DgwF+hk0t98ZOVleXllMbGRhkZmZ/aXcAwGAyJRGr1tfd87IqiGn9vgoau4pbzj76F38944irj6I8pqaut61OTo38WmtbRZem6F5mM2E1lxz2Wm7vl/wXk/SO37X197Pz4ow8ePKisrOyZZ13bLPRrekNaRxkbG5uZmb1//75tWsfj8QoKCjo6Ol3V1uDBg9XU1EpLSwsefSWQibKDlFj1LXkPMllUJgDA2dm5k8epqKiQk5MjEAgIgkRFRR06dCgrK8vQ0BDdmpub2+68mg5662ZmZl+/fkVfU6nUvnZznQAwUiLrn16SnuhOtB5/JPlcTl3+6THej27+ff78eWGHBv2LMtq13HsJqySXoP7PygEDlUzmmcz4M/rAhQnH2Wz26NGjw8LCeuvSeKKd1isqKurr69E8GBsbm5aWxr8CV/fB4XDHjx+fO3cum8HK9HvPv2nQoEGd/30XGxs7atQoAEBtba2enl5sbOzAgQN5W3V1dX+4nAWCIKWlpXl5eUpKSra2ti4uLsuXL+dyuY2NjZKSkgMGDPi5DwZ9H8JiNgRebc5MVVx9rElOfvOrvdJESvmlnItJfl5eXgAAgc3Cgn4II0aijJ3ZEHxLYdk+XqGb4aS8+oLDiacPjNquq6v78ePH3noJRLQvmdbW1jo6OsrKyiooKMybN+/cuXODBw8WTNNz5sy5deuWsrIyf+GMGTNCQkI63wUwMzNbsmQJAEBeXt7T01NbW/tnH2/N4XAyMjLWrFljaWmZm5tLJpNPnjwZHh4eGRn5119//dShoA6wq0oqfdZxGY3Km85WkImeYVs1xdQ2D/GcOMF5z549wo4Oaoek7SR2eUFLzif+wvXDPKgt1Juf7l+9ehWPx6elpQkrvG4l2r31AQMG5OXl0el0Lpcr+GHNefPmTZ48OTIyMjMzE13qa9CgQT91BCMjI95rtNv+s/B4fKsReVVVVVVV1V84FPQ9TR9e1/99nuIwi2w39W1Z2sGEU06yo/dM/9Px1ehVq1YJOzqofRgcXmrCHGqwv+Lak7xCPBa3f+S2FSGbtKT7saubVqxYkZOT81PL84kE0U7rKCE+kVZKSmratGnCah3qbv8MvGS8VVh5iKCmHZgTeu3j3TlKU621LKzDzAUzmxb6ZRLDHBpjnjSnp4gb//vcMVlxGW+7nVte7T019uDdu3d7X04Hoj4IA0Hdh11VUum7gVNfrbTpDE5V69L7G4+yXuh+UvbddkJTU/Nnf5lBQoDBSDnObwi8Cv47UV1fVmf9MI8dMQeNBg8YN25cXFycsALsJjCtQ1A7mj4mVJ3ZLGk5Tn7x7hY8bmfMoS8VmQukXZfOWhwZGSns6KDOIpkNx4hLMN5FtyofrWkzWtNm/+vjGzZtiI+PF0Zo3QimdQj6D4TNqn/s1xB4VcHjINluahWjenXYdjlxmeC1D79+zjIyMhLEzGuo60i7LKYG30A47FblywcvJOII3xRLV65c+fjxY6HE1k1gWoegf3FqK6vObuHUVyltPE1Q102vzvII2azNVJ+m4BgVGQUfmSKKxHRM8Uoa9MSQVuVYDGa3zeaUsvehBVGrV69OT08XSnjdAaZ1CPpH06eESt91EkNHyy/egyWRXxXE7Yg5OEna4cr68xgM5mcfMwv1HNIuSxrDAxBmc6tySYLEgVE77nx95PfkL01NTaHE1h16w0yYrpWfny/EqTV9SnNz679mwoJw2NTgm4z30fJL9hK1jBCA+H+8F/rt1cBCvVlOrnOyZnzvtl5IJBDUtMV0TWkxTynjZrfapCmlvtNmw9Gks/EhcboqOitWrBBKhF0LpvX/GDRo0NSpU4UdRQ+FIEjXzgbDYrGKiopdeMBfw6mrrLlxCEeRVd58AStBZnKYx5LOldHLdT4pZ376Ir5OHIuFP2pFntTERZU+6yRtJmIlWt/gYqk61NVgYgiILHocA9N6L/T7S+P2ViwWi8vl9r5Oa9PnpLp7vhR7V8qYGQCDqWmq2xnjTUEk+mcobtm9RdjRQV0GL69CGmTbGPlQ2mVx261zTNy+1RdgV2KCgoImTpwo+PC6FuyGQH0Vl9MQeK3+8QWFpfsoY2cCDOZbfcGq0C1WauZppxLVlNWEHR/UxaQmzKUnhXDqq9tuwgDMNus1BQ3FO+979YK+HeytQ30Rp76qxv8QjiyNDrwAAJJLUw8n+g5mGEllE8PDw4UdINT1cFJyktaO1LC7sjPXtt1KxBEPjd61pHFdgzRd8LF1Ldhbh/qc5i9JFSfWkMyGyy/Zi+b0vzMDjyef/3PwukDfxwJbLQ4SPIrDzKaP8eyKona3KpDkvO13nv98/bT/WQEH1rVgbx3qS7gcalgA/U2YwpI9RG1jAAAH4Zx5e/ljxRepKKyMMaWPPMK4z8KSyJTRrtSQ23IL/2y3gqmi0SztKZeir5fvLsMimBXLV4jixEfYW4f6Ck59deXZrcyibOXN59GcTmU2bn61t7apzoVrT2jBdeGDVqAei2w3rSUvnVmU/b0Kc8ynZ/ul3Xn34O7nR8OGWzY0NAgyvC4B0zrUJ7R8fV/ps07ccIjC0n3oFLciasmq0K3yHJn3xxNcJ7teunTph4+IhHoBDIFIGTeLGuT/vQpfvnwRI4trzzHWnGEgpkpKSUkRYHRdA6Z1qLfjcqght2sDfOTdd0g5zgMYDADgU1X6+shdc02mcxLp7gsXCjtESKAkrR3Z1WUt2e0/Q0NDQ4NeRWMUNzZV0KmF9VpaWgIO7/fBsXWoN+M01NTeOATwRKWNp3EUWbQwKDf8yofb5nVGYWdf+Pj4CDdCSPAwOLyU84KGwOtKG3xBmzvslJWVz/qe3bp9C5XROG2+m76+vlCC/B2wtw71NtHR0ZPGjXWd5Jz4JKDy1FoxgyGKKw+hOR0ByPWPAXe/PDrtcOhbfPaaNWuEHSwkHBJD7BAuu+lzYrtbF85fUFFSccD/cJFylYAD6xKwtw71KiUlJbNdp+4b3r+JxZ22YMmXxGgpU0t0UxO7+WD8STqTwX5Y/5oWc+3aNeGGCgkTBiPt7F7/9BLJxApgce1WWTdpVSwzpYxWoUpWbrdCjwV761Cv8unTp4EqMpP0VWYMUFOVlylg/tNxqWbUrA3/U0Zcer3h8v5qWtOnTxdunJDQiRsPw1FkGW9ffa8CAUcYrT7CedNUBoMhyMB+H0zrUK8yQIHyoagiJLficWZpOa0FHRhNr/66InTzQPKAKwvO4jDYo0ePEolEYUcKCZ/0lKUNIbcQNut7FWYPdJUfoZKRlSHIqH4fTOtQ78FIiSA8PXvt5JGHNEo4ovw0KJhCoUQXxu+IObDJcqUxU9fLywsumw7xEDUNiWo69Pig71VQI6sM1DB9U/m+sbFRkIH9Jji2DvUG6JPqWnI/Ka45NlFZc+LydQAABCB30x89/Ro8jm2zauLytLT2J7RBfZnURPfqC39KWI3Hiku0W2GyvuPxoNPFsXne3t4Cju2Xwd46JPI4DTVV57ZwG+uUNpwmKGue9D0lqyinpKY0Y8+8uKKks2MOExqw9+7dE3aYUE9EUO0vZjSUFv3dZ5mO7GdNVCKZjhwkyKh+E0zrkGhryf1UeWotyWyE/OI9WHGJjIyMg0cOGuwy11pnEnTmSd3NIu9dB5cvX25kZCTsSKEeSsppAS3uOZdW3+5WHAY3SW98Ojf7xIkTAg7sl8G0DokwekJw7Y3DcnM3owumAwCKiookNaTEFUgSGhRpJRl1FfWjR48KO0yoR8PLKUuY21MjHnyvwiS98e9pX06ePtVzHtPYMTi2DokkpKWp7p4vq6pEcb0PXu7facU4TbGa/Cr8wywWjQkY3MOHD8NJL9APSY2fU354GXnkZLy8StutypKKZkoDVkQs7NqHPnYf2FuHRA+7urTSdwPA45XWneLP6X9nBl7I8A94EjBV22mkovXHdx9hToc6A0uWJttOagy7+70KLvqOL/MjDQ0N8/PzBRjXLxL53npNTU1ERERJSYm2traLiwseL/KfCOpYc3pKXcApisMsst2/DxNHl03/VJk+W3zSiunL0tPT5eXlhRgkJHIoY2aUey9hlRcSVNpZYN1azeJ0yiWPXZ45OTn9+/cXeHQ/R+R764aGhg8fPiwrKzt06JCNjY2oDH5BvwJBGiMf1D04I79kL39OpzIbN0XuKadVbjZcqaPU/9WrVzCnQz8LI0Yi27tRg2+0uxWLwUzUG0ccSKZQKNXV7TwNtUcR+b5tZmamgoICAKClpUVbWzs8PNzFxUXYQUFdj9vMqLtzgkunKm08g5OS5ZUXN5b+GX1gmOrQ/IDM9Z/WhYWFCTFISKSRbV3KY58x8zOI/Qe03TpJd/zCF6trM0pryqr9/PwEH17niXxvHc3pAAAikYjD4XC49lftgUQau6Ko0mcdliKj4HmEP6e/LUtbHbbdVWfSZIVxE52dnzx5IsQgIVGHIRClJsxpCGx/DTg5kuxQlYEDp5lPmTJFwIH9LJHvrfOcOXOGTCaPGTOm7aaSkpLs7Ozc3Fz0raSk5KZNm1pdTGtpaYGX1zrAYrG4XK5Qmm5Oi6U9+4syabGY+RgmmwPYHLQ8OC/ixpf7W81XL5u4aPHixatXrwYAtLS0CCVIAE+hH+FwOEwmE4vt0V1J/CA79qtHjZ+SiAZD2m516j/20ocbKxTmnjhxojtWde7MKYTH43/Yee0laf3Zs2dHjhyJiIgQFxdvuxWPx0tISMjK/tPLIxKJeDy+1emFxWJ7+AknXOiXI+iviMuhhdxu/hgvt/IwXuXfh9RwEa5fmv/7yk9rNRfriGueO3fO1tZWoIG1B55CHUMQRAS+IiyW4jiX9vKGuOHQtk/YGKoykPme1SzFOnbsmKenZ5dP0OjM99OZSZa9Ia2/fPlyxYoVQUFBJiYm7VZQVlY2MzPz8PDo4CAEAgE+yrJjXC5XkF8Rl06tuXEIgycobzqLlSDzyhmsJq+E4xwu10N13lTnKeHh4fb29gKLqgPwFOoYFotFEKTnf0WEIXaM2Kfs9GTS4JFtt07Wd0yqf5ednU0ikbq+6S46hXr2v5ydEBER4e7u/uzZM3Nzc2HHAnUZZlF25ck1xH4GCsu8+HN6Ka3cI3SzAknertnCSMcwMTFx4MCBQowT6oUwGOlJSxqC/AGX03ajk87Y18XJbDxHS0uroaFB8NF1hmindTabPWXKFAUFhZMnT86cOXPmzJlBQd9dYxMSFYyUiOpLu6WnrZB2Wcz/Q/hTVfrqsG3TDJwxr5v37NojKSmpo6MjxDih3kpMzwwnq0RPbmdWlZQYZbj6sFdFr+fNm5eR0UPXYRftQRgsFvvs2TP+EgMDA2EFA/0+hM2qf3Sh5dtnpTXH8cr/WRg9KDf8ctqtzUNWNWbVzZ0718PDA16fhLqPtMuimiv7JczHYIhirTZN1nc8lnT2pvf5d6nvWCxWDxxWEvm07uDgIOwooK7Bqa+uuX4QL6+svOkshvjvpW8uwr2cdiu2KPHMuEMLp8y3sbFxcnISYpxQX0DsZ0DsP4D2OpAypvXzEc0UBxBxxA8Vn48cOTJ+/Phly5YJJcIOiHZah3qNltxPtTePkEdNoYydyV/exG4+EH+imd2yRG5W0aeCBw8eKCoqCitIqE+RnrSo8vRGSWtH/qs7qEl6455nh3h5eRUXFwslto6J9tg61Dv8f33dLa1yehWjek3Ydjlx2U3GK+dM/wOLxcKcDgkMXlGdZGpNi3rUdtME7TEp5Wkq2qocDic5OVnwsXUMpnVImJCWppobh2gJwUrrfcUMBvNv+lKd6RGy2V7LFvO6GYOA/Pz84cOHCytOqG+SclpASwjiUOtalUsQSKP6DQ/OjaDRaLt37xZKbB2AaR0SGnZVSaXvBgyeoLTuFE5OiX9TZEHczhjvTVae3HdNQUFB0tLS7d5oBkHdCictL2k5rjE8oO2maQYTn2W/nOY67cyZM4IPrGMwrUPC0Zz+purMZvKoKXJzt2AI/85pQQBy/WPAlbRb203X/n0yYP78+a9evZKUlBRiqFBfRnGYxXgXza4ua1WuJ6stJy77tjyNSqV6enoKJbbvgWkdEjh0fd2H5+SX7pUc/p85Lc3s5t2xR1LL0y5MOP7XUb+ev7A11OthJaXIo6ZQQ2613TRZ3/F5doiJiUlYWBiLxRJ8bN8DZ8JAAsVtZtTdOc5l0JQ2nsZRZPk3VTfV7og5qCvT36hE8+Yl/2vX2l9ID4IEjDLatdx7Caskl6Cuy18+tv+oi+/9aYCenZ1dX18vIyMjrAhbgb11SHBYJd8qT3hiKbIKqw63yuk5dd9WhW4ZrjZs8zDPK5cuw9sRoJ4DI0aiOMxqCL7ZqlwMR3TobxecG0mlUg0NDel0ulDCawumdUhAGO+iq/x2SE9cJDtzLQb3n5+J0YXxm1/t9RjoHucbHh8fHx0d/b1V2yBIKCRtJrLLC1tyPrYqn6Lv+CInVIIs6eHhUVlZKZTY2oJpHep+XE5D4DXqy1uKnkdIQ0bxb0EAcjf90fnUq8fs90nXSVAoFCsrK2GFCUHfg8HhpRznNgS1fiSelnQ/dYpqYsmb/fv3Z2ZmIggilPBagWkd6l5cWkOV305WWb7SxtME1f78m1gc1qEEn7iipLX9Fs8c42ZiYnL69GkxsdZLcEBQTyBhMRZhNjWnv2lVPlnf8Xl2KABgz5494eHhwgitNZjWoW7ELMquPLWWqGmgsGw/lvSfO7BrmurWhP/J5nJ8HbxzP+d4eXnBmelQj4bBSDktaAi8Bv7bJbfTtMmpzytuLL148aKWltb39hYkmNah7kJPCK6+uLPt+roAgNz6/FWhW6zUzJkh9QvnLnB3d588ebKw4oSgTiKZWmPEJRjvovgLCVi8o/aYFzlh5ubmwcHBhYWFwgqPB6Z1qOshbFbdPV/a60ClDadJZiNabU0qfbspcvfCAbP+MJgmLyffwx/iDkH8pF0WU4NuIBw2f6GL/oSQb5FMDrOqqury5cvCio0HzluHuhinvrrm2gG8gqrSeh/+9XVRf2cG3st4sm3gGvdJ8w8fPrxx40ahBAlBv0ZMxxSvrElPw/vMTwAAIABJREFUfEm2deEVqpFV9GV1Y4sS9+zZQ6VShRgeCvbWoa7UkvOp0mcdaZCN3ILtrXI6i8M6knj65bfItVqL9WS0T5w44erqKqw4IeiXSbssbgy7i7Q08RdO1p/wPDtEXFz82rVrDx8+FFZsKJjWoS6CILSYp7U3D8vN+3d93ZaWFt/TvstWLn/0/PHmV3ub2M2WZcYL3eZLSEi4uLh0fDwI6pkIatpieoNosf95LtsIDctSWkVefcGwYcOePn0qrNhQMK1DXQBdX5fxNlJpva+Y/r/r6y5fteL4bZ/opsT5yxa0fKaNZFrYjxydkpICl+6CRJrURPfGmCdc+r/jLTgMbqLuuMCcsLFjx167dq25uVmI4cG0Dv0udH1drARFcb1Pq/V1g4JeaLobqo3TVnfWeXL8fnJSkoGBgYKCgrBChaAugZdXIQ2ybXz1n8GWSXrjwvOjm9ktt27dEu6ajjCtQ7+l+Uty1ZnNZLupbZcEoDHpEiqU6uQyZn0z7WPdWs+1O3fuFFacENS1pCbMpSeFcuqreCWKEgqDlExeFcRNmzZNuL9HYVqHfhW6vu7f5+WX7pW0dmy18V35x8XBa//YuYD4GaRujh6hZ7l161ahhAlB3QEnJSc53Ika9p8nbLjoTXiW/VJeXn7Pnj3p6enCig2mdehXcBmN1Zd2NWekKG08TdQy4t/E4rAuvb/hnXBqvvb0kme575NSmU3MR/cfEQgEYUULQd2BMnZG08d4VsW/9x9Zqg1tZNKyanMyMzNnz54trMBgWod+GqvkW+WptXglDYWVrdfXzWso9AjdUkorvz7x7OfgNHNzcyKR+L3jQJBIw5LIlNGu1JDbvBIMwEzUHfc8O8TGxmbVqlVCC0xYDUMiipEaVeW3Q3rSYhnXlfyD6QhA/s4MXB++c6qBk1gc58Aurz///HPDhg1CDBWCuhvZbhozL51ZmMUrmaQ3PrYokcFumjx58vXr14USFUzrUKeh6+uG3Fb0PEIaPJJ/SyWjekPErlcFcafHeDtrO5SVlW3atElYYUKQwGAIRKlxf1D5nrAhLSY1TGVwWF4UHo/ftGmTUGY6wrQOdQqX1lB/ZS+rvKDt+rrRhfHLX24cpGTq2W+h3ZCRKSkpZ8+eVVFREVKkECRQksMd2bUVLdlpvJLJ+o5PvwYrKSnFxMQIZaFpmNahH2PmfanxWUvoZ6CwdB//+rp0FsM7wef6x7tH7HYPYGgT8ITr169bW1sLMVQIEjQsTspxHv+CvYOVzTAA86kqQ1VV1d3dXQgRCb5JSLTQE4Jrrh+UmrFW0mkB//q678o/LgpaI4YjXnI8dXTLIW9vbwMDgzFjxggxVAgSCokhdgiX2/QpkVcyUW/88+wQOTm5xMTEkpISAccj8mnd19d3yJAhRCJx27Ztwo6lt0FYzLp7PrT4F4rrTokZmfPK2VzO9Y8B3gmnNlquNGcYpyS+2b59u9DXwYAgocFgpCcubHhxDXA5aIGTztjEkhQqq/HTp0/q6uoCDkfk07qhoaGPj4+bm5uwA+ltOPVVVWc3Iyym0jofvLwqrzyvoXBFyKb8hsLrE8+qsZUWLlxIIpEMDQ1xOJwQo4Ug4RIfMAwnLU9PiUTfkomSthpWod+iuFyurq4uk8kUZDAin9adnJxGjx4Nl47qWi05Hyt91pMG2crN34Yh/nPNBwFIYE7oxsjdUw2cFvWb5bFohbS0dGFhoYWFhXCjhaCeQNplMTX0NsJmoW8n6zs9/RosJi5ua2tbUFAgyEj6xGM0mExmYWFhamoq+pZMJhsaGgo3pJ4LQRpfPaTFPJGbv41/LcbaprrjKecZ7Kbz44+qkVXOnj1raWkpJSUlxEghqEchahoS1XXp8S/IdtMAAMYKBmSi5PuKjzdu3BDwk/D6RFrPzs4ODQ0NCQlB3+Lx+BcvXkhISPDXodPpmP8+b7MvYjbTH51H6qsklx9iySiyaDS0+HVp8oUP1x01x/xhOM1jiYeUlJSPjw8AgE6nCzXcngWeQh3jcDhMJpPD4Qg7kG6Et59JvbYfMbXFiJEAAOP7jX6c8UJfQnvo0KHJycnKysod796ZU0hMTOyH63D0ibRuYmIye/ZsDw+PDuogCEImkzuo0Ouxq0pqrnoRdUxkFu/k3T5KZzH83l3/VJV+xG53SwkDj8Pb29u7u7uLi7d+mh0ET6GOoWmdRCIJO5DuRB7AGWCBpIRSHOcBACYZTbiZ8YBJYHt5eXG53B+eHl11Con82DrUJZq/JFWd2Uwe7cq/vu7nqsylwesBAJccT32IfDdx4sTq6moPDw+Y0yHoe6SdFtDinnMa6wAAJLy4vZZtcG74qlWrMjMzBRaDyKf1z58/P3z4MC8vLzMz8+HDh1+/fhV2RKLmn/V1L8gv2ydpPQEtQ6cw7o07sm7Ycgcxm2t/XR0zZkxqaqq2trZwg4WgHg4npyRhbt8Y8QB9O8XAKTA7lItwt2zZkpKSIpgYRD6t5+bmRkRE6OnpqaioRERECPjShKjj0qnVl3Y1Z7xV3nSGqPnPZeS8hkKP/09htFazWLVqlZKSkpycnJycnHCjhSCRIDV+DuNtJLumDACgK9NfUUIhqTT18uXLqqqqP9y3S4j82PqUKVOmTJki7ChEEqskt+baQXETS5mpywEWBwBAAPIiJ+zax7uLB84xxuhNdHA+duxYQkICAIDFYgk7XggSDViyNNnWpTH0ruycTQCAyfqOgdkhh0fv9vb23rp1qwAePCDyvXXo1zBSX1X57ZSevFjGdSWa02ub6rZFeb3MjTwz9pAhV7uysnLFihXDhw8XdqQQJHooY6Y3Z6ayygsAAPZatunVX8vplaGhoVFRUQJoXeR761DnHTty5NG9O5qaWnumjFIsz1JcfZSgooVuii6M9025NEXfaaHZrDH2YzQ0NG7fvg1zOgT9GowYiTxmOjX4hvziPWI44njt0S9ywu7evSuYWz1gb72vePz48Z3zPruMyHq12UsPnVbaeBrN6XQW40Ty+esf7x4Ysb008FtmRubly5dv3779wwNCENQBss0kZlEOMz8DADDVwDk4N1xFTdXDw6OhoaG7m4ZpvW9AkJTIUCctmSEq0kuHaH0qqUHX1+WfwpgVn56dna2mpqavry/scCFI5GEIRCnHuQ1B/gAAdYqqlnS/18VJCIJER0d3d9NwEKaX4zTUMN5G0pNCBzeV7fxaoS8rGV/aMHK4FZvLufX5wYuc0C3Wq28evDb9gNuLFy+mT58u7HghqPeQHDaOFvWoOTNV3Mh8sp7j8+yQmzdvCmBRPNhb750QNqspLa768t6KYyvZNeWyf2ycdS143qr1h7/U5lHU9188jk5hnI2fpEfsP2LEiICAAGGHDEG9DhYr5bSgIfAqQJCR/YYXUIvLmyqtrKyKioq6t9luPTokeMyi7PrHfmX75tMSgiUGj1Tde0t25loxHZM3b974XjiNWEtm0PNmLJs9xcBJ6SP5+MFjLS0t8+bNo1Aowg4cgnoh0iBbDIHYlBaHx+IcdcYE5oQ6Ozt/+fKlWxuFgzC9BJfxP/buO6CJ830A+Ht3WSQhEPbeSwFBESfWAbgAt4KD0WUdrVZrv63Wqm2/1tHaat1WreIWR924quJmWVREZYvI3hmQXO5+f1y/+VGcKEkYz+evy+XGEzge3rz33vPWyf++Krl+ilbU8/0Gms1dwzL6V12hEydPiAdY2IQ4UwrVrennLq47s3r16pkzZ+oqYAA6CIPQD6sOrOZ16RPmMuSTM1/ELdpWWVap0TNCWm/jKKohK01y43RDZhqvs7/hyI+4rr7o30XgKJq6/OT6lbrbtQ8rySBl1f0yQyPDb7/9VlchA9ChcF28CbGZ7PY5iz7DPYxdLz+5/tnQacePH9dceXBI622VsviJLOmCLOk8YWwl8A80mjiXqQX6r20o8q+8hN3pcQZc0Y9zvp+R8EniZxccnB3PHD9jbGysk7AB6IAMwt6v2Podv3vgCNeh+x4cnjt3rkaHOUJab2MouUR+J0GadFFVVcL3G2Q66xeWyQsKTciU8tPZF/ZnHHUVO870+nDD4t9Suycm3Uhks9k4DjdUANAqjq0bx6Gz5Orx3oPGrkneEjjhk+L0Z5o7HaT1NoKmGzL/liZdqE+/zXXrajBsyvOdLYyq+po/H58+lnna37Lbkp5fPr6VYedm5evr++GHH3K5XO0HDgBACBmExpSumSvoPTTEOfhE1rnNH63et29f165dNXEuSOutHVlSILtzRXr7HCE0EPQeJh736fOdLYwiSUncw+MX8q4EOfTfPPQXc4HpwIEDPTw8xowZM3fuXC2HDQBojGVqrefVu+7S4bBBI6NPfrry1580N/wM0norRdVL6+/dlCZdJEue6PkEmEz9nm3p8LKNs6pyD2T8mVR0Z6TrsN1hG+8l3w3w7fPHH39op64QAOBNiIZFlqycbh4Q1tXCG7PknD9/3tnZWRMTJUJab2VouiH3gSz5ovzvqxzHzsI+w/W69GEqLL7QvbIHe9MPZ1fnT/AY+UWP6aePn/7z76P+/v47d+7s06ePNgMHALwaYWAs6BFcd37/iL5DN975I2XzFV9fX00U1IO03lqoqstkKZekN05jbA7fP8jim2244KXF3iiavvUsadf9uJqG2jFuod8FfPX44WNpjXTz5s2LFi3q1KmTNiMHALwh/eCJxT9+5DNgdINKsfT35e7OGhnjCGldx2iloj79tuTGaWVhtp5PgPEH37KtnV+xPTNmcc+DwzyCO9Y9NNhxII5h48aNKygouHnz5tmzZ7UWOQCguXC+UNhvRN2Z3SF+wfdrH5/77syvv/7a4oPTIK3rjKIgU3rzjDz1CtvOTdhnOM+7t3pu6BeSk/Wnss7vzzhqITCd1jW6j3WPa9eu+YzssnHjxi1btsCMdAC0CfoDxxT/94PBAcOi0w8VJT2+c+eOn59fy54C0rq2/VNS8WY8xmLx/YMsFm7DhYav3qW6oeboo9NHH5/yNHFf+t4CJ5H9xo0bH+mn+/v7r1+/PiAgQDuRAwDeHcbh6QeF15+P6+XZ3ea3UB8fnxY/BaR1LWE6W6RJFxR5GXo+AeJJX3CdPF+7V7G09GDGsXO5l3pb+28YspKsUNw4e0PwHu/evXuzZ8/28vLSQuQAgJYl6BsiSfhzmH/31U/j936zI/5MfMv2w0Ba1zhFQaYs6YIs9TLbykngH2gc8w3G5rx2r5zq/H0PjtwqTA52HLAzdP3Tx08s+eaBYYFBQUGTJk36/ffftRA5AEATMIIlGjIFTzhNuBASYX1WVpabm1sLHh/SuqaoaivldxKkiedopYLfbYD53LWEkdmb7MiMWXxYmTXCZeiWQauM9Y3mzZt3/vz5hIQELcyrAgDQAn73QXWXDw8TeDrNtW/ZnI4grbc4WkU2PEyWJl38X0nFj1/2lH/THRF9szBpT/qhSnn1GLfQhb3mnjl5pnP/TocOHVq5ciWXy9XEYwsAAN3AMNGwqB5n/oi1r/Xo2vlByv0W7IeBtN5ilMX5sqSLssTzhImVwD/QaNIXGIf3so2lUmlsbGxdXd3kyZPNLS0u5l3Z++AwG2eP9xjhRNkOGzosIyrqiy++yM/P185U5QAALdPz6iW4GBcgMKaH1ZeVlZmbm79+nzcDaf1dUTKJ/O8EyY3TlKSG7zfQ9PNfWMYvKKnYGE3TAQP7VXCrCRF76c8/Dvg5xMbc+gOvSYuiF6zVzzp8+PCRI0e6dOminfgBALpiEPbBoLgf7/YwIVWqFjwspPW3pZ6/4mEK18PPcMSHb9jZUi6ruJF+Oys/q9svAxBCDVXyop1Z6fm312et/Omnn7p27cpmsyGnA9ARcJ08ncVO3Pr8wOjBGefutdRhIa03m7LkiSzxgizpPGFgIug97IXzV6jVk/V5NQVZVbk51fk51XnZVXkETpjSYqVcKS+WskUcxbP6L1bMHT16NEKoR48eWvwcAADdMwj7cPDeL6tHejU0NLTUMSGtvymqXipPvSJNuqiqLOZ3DzSdtYplYvX8ZuXyyseVWXk1BbnVTx5XZj2TlFgJzR0M7GwEVva1FjE+4ScOHNu2bdviBYt+/GEZSZIzZs4IDw/X/scBALQGbEuHgWbd9hGpSenJvm4t82gSRtN0ixxIV4qKihYuXJiRkeHr6/vDDz+8cC63mTNnent7T5s27RXHqaure3H54/+fvyKR6+Yr8A/kdfJXl1SUKmU51XnqJJ5Vlctn6zkY2DkY2LkbOTsY2nElrA3rNvTr108ul+/YsWPp0qWNC+erVCqCeGl1xlZFqVRSFAUTcbzCSy8hgBBCSKVSKRQKPb2XfrXtsMiK4uV7pp3OKDu7Mr5FLqE231ofPXp0t27d1q9f/+uvv0ZGRp4+fbqljkyWFcpSLskSz+MCEb97oOGoT2i+oERalleU8qgim2mS1zTU2uhb2RvYuhu59LXp4WzocDcxjc/n83i8ccPHRUREjB8/3tDQ0Nvb28HBYcKECU1O0VZyOgBAc1jGFqMset0yuK2iW+bGadturScmJg4ePLi0tJTD4dTW1pqbm9+7d8/FxaXJZq9urZMkOW3atGs3r77Xt/+mTZuQor7+3g1m/gpVl55lnTpn4/K8mid5NU+yqnLFPEMHA1s3IxdHQzt7A9uKrDIVSbq5uc2YMcPGxmbGjBnvv/9+dHT05MmTc3JyXF1d203ihtb6a0Fr/dWgtf4KqtrK/vMGF9ws9XPtvmPHjncc1ty2W+tpaWldu3blcDgIIZFI5OHhcffu3efT+quNHDXi6r3rFkEO+84cyBlw6+MQlyJz42c2wsfmhFJ10/FpoRluYsY2nuobdX7fWRYiPvnkkyFDhhgbG+/YsSN6zpSgoKBu3bpFRER07drVzs5O/SCoh4dHi39eAEC79P2qNSmHH9iMcElIvtazd4+M9IfvcrS2ndZLS0vFYrH6pVgsLikpeX6z+/fvHzhwYMWKFcxLgUDw119/qVsNNxJvOE71MuxkrGcuSD36eJPEvCG3avmXSz8a876YZ7jv/OaQkBB7e/uIdWOqy6ucnZ2VSuWKFSusrKwwDDt16hRCqL6+ftCgQQihuro6jX9mHdm3b19NTc2r7090cGFhYQcPHoTW6MucO3cuNTX166+/1nUgrdHBQwdsR7paBtob+1mk/ufyKzIJj8djs9mvPlrbTusikUgmk6lfSiQSAwOD5zdzd3cPDg6eNGkS85LNZpuZ/X95FgtT88rEIr6VsCLxmRXbbP+nu/X19fl8/uO0R8wGV69eZRZ++OEHZqFXr16a+DitWUVFRVVVFXQyvML9+/dZLBb8iF6mqqrq2bNn8PN5IQ83j6spN427mZffLBQKBe/4U2rbad3BwSEzM5NZVqlUubm5Dg4Oz2/GZrNNTEycnJxeeJBjR0/0D3wvKeGiuZnpsUsXW/ARXgAAeBO7du328++WNO8voVBw+MCRdzxaC0+2pGWDBw+uq6tjekL27dsnEoneoh3t5uZWVFBMK+nCJ0XOzq+acA4AADRBKBQ+ynjsbO+ceCMpMDDwHY/WtlvrXC53586d77//vkAgUCgUe/fufWEVtIaGhidPnqSkpLz6aK/doCMrLCysra2FH9ErkCSZlpYGpdleJj8/v6KiAi6hV6ivr09PT3/1LTorKytLy9dUnWrbAxwZJEmWl5ebmpq+bDRhZGRkYmLiq7urcnJyXtZLAxBC1dXVKpXqhU97AQbTBwj1k19GIpHIZLLGt7VAE/n5+dbW1izWq1rbISEh33333auP0x7SOgAAALW23bcOAACgCUjrAADQrkBaBwCAdgXSOgAAtCtte4Dj29m3b19cXJxIJJo1a1a3bt2e3yA5OXnt2rUSiSQ8PPz5sovtnkQiWb58+d27dzt16rRgwYLnH9yNjY198OABs2xoaNjRHge/d+/e9evX8/LyQkJC+vXr98JtTpw4ERsby+Fwpk2b9rJt2iuSJP/666+UlJSampr58+e/8MHvlStXVlZWMsudOnWKjo7Wbow6duXKlaNHj+bl5VlZWU2fPt3b2/v5bR49erRq1aqSkpKQkJCPP/64WSOsOlxr/cCBA/PmzYuJienWrVtgYODTp0+bbJCfnx8UFOTv7x8dHf35558fPnxYJ3HqUHR09L1792bNmpWfn//C/2pHjx4tLi52cnJycnKys7PTfoS6tWLFioSEhMOHDycmJr5wgwsXLsTExEyYMCEwMDA0NDQ9PV3LEepWWVnZt99+m5OTs2LFitra2hdus2nTJhaLxVxCFhYWWo5Q51auXGlhYfHRRx+ZmJj07t374cOmhb1qa2v79+9vY2Mzbdq01atX//bbb807Ad3B9OzZc9u2bcxyRETEokWLmmywYMGCSZMmMcu///573759tRqfruXk5HC53MrKSpqmZTKZQCC4f/9+k21GjRq1fft2XUTXioSGhv78888ve2vFihXM8meffTZjxgwtxtVaSKVShNCTJ09e+K6jo2NKSoqWQ2qdBg0atGrVqiYrN27c2K9fP2b5zJkz9vb2FEW9+TE7VmudoqjU1NSAgADmZd++fZOSkppsk5yc3LdvX/UGycnJdEca2p+amurh4cHUxdTT0+vWrVtycvLzmx04cODDDz/86aefampqtB5ja9fkEnr+GgMIoVWrVk2dOnXLli1KpVLXseiMSqV68uSJjY1Nk/XJycnqNBUQEJCfn19aWvrmh+1Yab2yslKpVKpr+RobGxcXFzfZpqSkxMjISL1BQ0NDVVWVVqPUqcYfH73kR9S/f/+RI0f269fv0qVL3bt3b8fliN8CRVFlZWXqn6GRkdHzP0AQHh4eGBjYvXv3TZs2hYaGdqiWU2Pff/+9gYHBmDFjmqwvLi5WX0JCoZDL5TbrKupYt0wFAgFCqL6+nnkpl8ufryjA5/Plcrl6AwzDmL06CIFAoP75IIRkMplQKGyyzeeff84sREZGdunS5fDhwzExMVqLsJXDcVxPT6/xJQSlaJ+3bNkyZmH8+PE2Njapqal+fn66DUn7NmzYsHv37itXrjxfLUAgEKgvIZIkFQpFs66ijtVa19PTMzExyc3NZV7m5uba2to22cbOzi4vL49ZzsnJMTc371AzvTEfX916ys3NfcVNUYIgnJycmvX1sCOwtbVVX2N5eXnPX2NATSwWm5iYlJWV6ToQbdu2bdvKlSsvXLjwfA8M+ncWysvLIwjiteW9GutYaR0hFBER8fvvvyOE6urq9u/fHx4ejhCSSCRr165lbvKEh4fv27dPIpEghLZu3cps0HH069cPw7CTJ08ihK5cuVJRUREUFIQQunPnDjMoiCRJ9QWXlpZ2+fLlPn366C7e1iIvL2/r1q3Mcnh4+LZt2yiKamhoiI2N7WiX0MtcuHDhypUrCKGqqqqKigpm5ZEjR8rLy319fXUamrbt3Lnz22+/jY+Pd3R0bLx+7dq1zPxu4eHhx48fZ5a3bt06cuTI5s261TJ3c9uOoqIib29vLy8vKyur6OholUpF03R+fj5CqLCwkKZpkiQnT55sbW3t5eXl6+tbXFys65C17dixYyYmJj169DA2Nt6/fz+z8qeffmIGBUkkEpFI5Orq6uXlpa+vv2zZMp0GqwOzZ88Wi8VsNltPT08sFsfFxdE0feLECRMTE2aDqqqq3r17u7u729vbjxgxoqGhQafx6oC9vb2hoSFCyMDAQCwWK5VKmqYjIyOZQUFJSUlCobBz585ubm6mpqYHDx7UdbzaZm5uzuPxxP/zzTff0DRNURRC6NatW8w28+bNMzMz69atm7Oz8+PHj5t1/I5YwZGiqEePHunr66u//tA0LZfL9fT01GP+CwoKpFKpu7t7x6yzKpVKs7OzHR0d1T16JEmqVCqmP4okyZycHJIknZyceDyeTiPVAZlM1tDQoH4pEAg4HA7TNm/cpHr8+DGHw3nhdF3tXnV1dePEwgxSUCgUGIYx83DK5fLc3FwOh2Nvb//amTnbn5qaGiaJM3g8HnPlyGQyHo+nnjSipKSkrKysU6dOLys5/jIdMa0DAEA71uH61gEAoH2DtA4AAO0KpHUAAGhXIK0DAEC7AmkdAADaFUjrAADQrkBaBwCAdgXSOgAAtCuQ1gEAoF2BtA4AAO0KpHUAAGhXIK0DAEC7AmkdAADaFUjrAADQrkBaBwCAdgXSOgAAtCuQ1gEAoF2BtA4AAO0KpHUAAGhXIK0DAEC7AmkdAADaFUjrAADQrkBaBwCAdoWl6wAAeDGapi9cuHDp0qWKigqRSBQcHBwcHIxhmHqD+/fvZ2Rk9OzZ09zc/PDhw4mJiQqFYsWKFeXl5cnJyV5eXu7u7qdOnbp8+bJUKl2wYIGdnR1CqLKyMi4u7v79+yRJuri4jBs3zt7evvF5r127VlRUNGTIEJVKtW/fvoyMDJFItHTp0leEWlFR8eeff969e1ehUFhYWPTp02fgwIEsFgshlJ2dnZqa6uPj4+bm1niXCxcuVFVVjR49mtmsqqrqwoULNjY2vXv3vnHjxsmTJysqKqKjoysrK+Vy+ciRIzkcTuPdKYo6cuQIQRCjR49Wr5TJZH/++WdiYqJcLre2th4zZoyXl1eTUB8+fHjs2LGnT5/SNG1iYuLv79+/f3+hUNisXw1o7WgAWp+ioqI+ffowlyiT+BBCgwcPrq2tVW+zcOFChNCaNWt8fHzU1/OzZ8+2bt2KEFq4cGFgYKB6/fXr12maPnfunJGREUIIx3E2m40Q4nK569ata3zqkJAQhNDevXstLS2ZfS0sLF4R6u7du0UiEbMll8tlFubNm8e8u2HDBoTQqlWrmuzVrVs3hJD64yQnJyOExo8fP3PmTHXM69evnzx5MkLoyJEjTXa/ePEiQigkJES95vr169bW1o1/YjiOL1iwoPFe8+bNw3EcIcThcHg8HrPx4sWwHuumAAAgAElEQVSL3+AXAtoSSOug1amvr+/atStCaPr06Xl5eTRN5+fnR0REIIQiIyPVmzFp3cDAYODAgQkJCeXl5ffu3ZPJZExaNzAw6NKlS3x8fGlpaVZWVllZWVZWlkAgIAjit99+k8vlJEkeOXLE0NAQw7BTp06pD8ukdQMDg5iYmJSUlKqqqpSUlJeFevLkSRzHhULh1q1b6+rqaJouLCzctm3b1q1bmQ2aldYNDAwsLCz279//9OnTwsLC3Nzcc+fOIYRGjhzZZPeoqCiE0MGDB5mXmZmZIpFIIBD89ttv5eXlJEkmJSX5+voihLZv385sc/PmTYSQt7d3Wloas6a4uDg2Nvbw4cPN+/WAVg/SOmh1Nm3ahBCaMmVK45UkSfr6+uI4npOTw6xh0rqZmRmTT9WYtM7lcvPz8xuv//jjjxFCs2fPbrxyx44dCCE/Pz/1Giat9+3b97VxUhTl7u6OEDpx4sTLtmlWWkcIxcfHN95MpVLZ2tqyWKzi4mL1SolEIhQKjYyM6uvrmTXh4eEIoW3btjXeNz8/n8/nOzk5MS/Xr1+PENqwYcNrPxdo6+CWKWh1Dh06hBCaMWNG45UEQYwfP56iqCtXrjReHxkZ+cKu4ZCQEKYzXe306dMIoVmzZjVeOWnSJBMTk5SUlJKSksbrp0+f/to4Hzx48OjRI3d399DQ0Ndu/CacnJyGDBnSeA2O45GRkSRJ7t+/X70yLi5OIpFMnDiR6fNRKBTHjx8XCoVME17Nzs6ub9++OTk5+fn5CCErKyuE0JEjR0pLS1skWtBqwS1T0Oqkp6cjhObPn890f6sVFRUhhAoLCxuvdHFxeeFBmqyXy+WFhYU8Hs/R0bHxejab7eLiUl5enpWVZW5u/trDNpaRkYEQ8vT0fO2Wb+iFJ42JiVm2bNnOnTtnz57NrNm5cydCKDo6mnmZk5Mjl8uFQuGwYcOa7Mv8JAsLC+3t7YcPH96tWzf1jdkBAwaEhIR0796d6W0H7Qn8RkGrI5PJUKM7pWqWlpZBQUE2NjaNV6pvVzbRZH1DQwNCyMDAoPFYGoZYLEYIyeXyNzlsY/X19erdW4SBgcHzK11dXXv37n3nzp20tDSEUH5+fkJCQufOnf39/ZkNmB9Xk3+BDE9Pz6CgIIFAgBDicDjXr1//+eef/fz8rl+//v333/fs2bNXr16ZmZktFT9oJaC1DlodIyOjmpqaLVu2ODk5tdQx9fX1WSxWWVmZQqFoMljw6dOnCCFjY+PmHtPExAQhxHRxvAyTbUmSbLK+qqrqzU8UHR1948aN2NjYVatW7dy5k6KomJgY9bvM2B6RSHT+/PlXH4fH433xxRdffPFFbW3tpUuX1q5de/HixcmTJycmJr55MKD1g9Y6aHV69+6NELp06VILHpMgCB8fH4qiEhISGq8vLCx89OgRj8fr3Llzc4/p7+9PEERKSgrTXn4hZtBhk46j8vLyV/8zaCIiIoLP5+/evVuhUMTGxrJYrClTpqjftbe3t7Kyys/Pz8nJecMDikSikSNHxsfHW1lZJScn19XVvXkwoPWDtA5anZkzZ2IYtnjx4ry8vCZvZWRkMN0pb4HpjF6yZIlSqVSvXLx4sUKhiIiIUA85f3PGxsbh4eFVVVVz586labrxW+pEyTyFdPjwYXXqpyhq3rx5FEW9+YlEItGoUaNKS0sXLVqUnZ09ZMgQ9Zh6hBCGYcxo92nTpjXpSqJpmum6QQg9fvy4yb8fhULR0NDAZrOf7+8CbRr8OkGr06dPn++++27RokU+Pj4TJ07s0qULi8XKzc29evXq9evXy8vL3yIFI4SmTZu2b9++69ev9+nTJzo6Wk9P78iRI6dPn7ayslq2bNnbhfrrr7/evHlz8+bNaWlp48ePNzExyc/Pj4+P79Wr16pVqxBCzs7OwcHB58+f79279+TJk0mSPHHiRHFxsZOT05s3rhFCMTExe/fuXblyJWp0s1Ttyy+/vHLlyrlz57y8vMaPH+/m5lZXV5eTk3Pq1CmBQMBk9h07dmzYsGHs2LGenp7W1tbFxcU7d+6sqKj45JNP9PT03u7jg1ZK1yMsAXixQ4cOdenSpfG16uzsPG/ePLlczmzw3//+VywWx8XFNdlx165dYrH4+aHiNE3X1tZOnTpVncUIghgxYgTzxJNaeHi4WCx+9OjRG8ZZUlLC/JNoHKf6QSFmg+DgYPW7/fr1y87O7t+/v1gsVo+4//vvv8VicUxMzMvOolKpOnXqJBaL7e3t1cPVG1MoFEuXLm3Siu/Zs+fGjRuZDY4ePdqjRw+CINQbiMXib775RqFQvOEnBW0FRv/7yyMArUppaWlBQQGPx7OysmqpMScymSw7O5skSUdHR0NDwxY5plwuz87OViqVFhYWjXOrWmFhYXFxsZmZma2tbYuc8YVoms7Nza2qqtLX17exseHz+c/HWVBQUFdXx8QJoxvbJUjrAADQrsD/agAAaFcgrQMAQLsCaR0AANoVSOsAANCuQFoHAIB2BdI6AAC0K5DWAQCgXYG0DgAA7QqkdQAAaFcgrQMAQLsCaR0AANoVSOsAANCuQFoHAIB2BdI6AAC0K5DWAQCgXYG0DgAA7QqkdQAAaFc6RFrPzMwsKyvTdRQAAKANHSKtr169+siRI7qO4i1JpVJdh9ChyeVyiqJ0HUXHpVKp6uvrdR1FG9Mh0jpCqO1O2Qo5Rbfg569bNE233T9eXWHpOoBWISEhITExsVevXgEBAbqOBQAA3klHaa2/wvYdO0MmfjD/UumwCdG7du/RdTgAAPBOIK2jTTv2SMLXkqOXSiasWbt9t67DAQCAdwKdMMjawjy14I6qcxCWn/q3wmTSJdVkF3yINcaCf3kAtDLbt29funSprqNoeRwOJykpSSgUtsjRIK2jX5Z9fy909JOTSx1c3A8fPXqHxlbfV0Vfpsc64pEueF8LDNN1hAAARn5+/qhRo2bOnKnrQFpY7969pVIppPUW4+jomJX+t1wu19PTQwh5IxTlihdI6b1Z9EdXVRiGwp2wSBfcWQTpHQDdE4vFTk5Ouo6ihREE0YJHg46GfzA5Xc1WgH3lgz8czzoUSNSrUN8TZPc/yTX3qXIYQQsAaN0grb+Gpxhb7k8UTmIv9ydSymn3OGXYOTIul1LAaGYAQKsEnTBvhMBQkDUWZE3UKok/86jYTGrGddUwGzzKFQ+0hs53AEArAmm9eURsFOWKR7nihVL6UC79n0RVRQOa6Ix94Ia7GUB6BwDoHnTCvCVrATbbC08dzTo9hEAIDTj1T+d7qVzXkQEAOjZI6++K6Xx/OpG9uhfxoJp2i1MGnyFjMykpqevIAAAdEqT1loFjKMAC2xxAFE1iT/XA43Ip+33KqMuqC4VQpggAoFWQ1luYHguNd8RPDGbdHcvyM8GWpKrs9pGzb6r+roD0DgDQBrhlqilWfGy2FzbbC39QTR/MocZcUOkRKMoVj3bDLfSQVCrdtWtXbW3t5MmTra2tdR0sAO1TZWXlsmXLrl27JpFIPD09586d26NHjzfffcOGDZcvX1a/nDdvXo8ePebOnfv06VP1yqCgoKlTp16+fHnDhg0RERFjxoxh1ickJKxbt27ChAnjxo1roU/zpiCta1xnQ2xJN2JRV3SjhN6VRXWKU3Y3QY8WBJfr2apE5it+7fXo3h0TExNdhwlAe5OZmdmnT5/y8nLm5f379+Pi4tasWfPpp5++4RESExOVSmVUVBTz0tbWFiEUHx8fGhrau3dvZqWDgwNCKDc39+LFi7m5ueq0vmbNmsuXL3t6ekJab7eYzvcAC2J1L2LbrfzPC5+qll1BCNXXFl+4+FdE+ARdBwhAu0LT9KRJk/7J6XZdkcgcPbxEkQ1ffPHFe++916VLlzc8jqur6+jRo5us7NGjx/Mrvb29q6qq7t275+3tXVFRcevWrcDAwHf+HG8D+ta1TY+F3u9qwlFKUWkWklU3FKR/+sBszi1Vcjl0vgPQYlJTU5OTkxFCKHg2+vY2mn0c/ecvhOEKhWLnzp1vfpyqqqrHjx8/fvw4KytLvbKkpCT7f2pra9XrIyMjd+3ahRDavXv3hAkTuFxui32e5oDW+j/qzu/nOHtznTy1cC6BQLBp3W8zPuunVDTMmPnp51/2359NT7mkqlehCHiyCYA3FniazKt78VvSxJx/lvzG/rPg6I9M7FFZ7uarOX8eePEAZDaO9g0iuhr//x/g0aNHb9y4gRASCASJiYnMyh9//HHt2rXM8rfffjt58mRmOTIysnv37j/++OOOHTv++OOPVatWvdPHe1uQ1v+B8/VlKX+xrRwkV47xu/VnmWr2NmbUlMlRUyarVCqmcttXPthXPnh6Fb0rixpwihRzMPXNVQDAy+wdyJKSL/6ae0toPPkXhBBChfeRcy+EEKorQ9XFCKGQTsbLhr+0YqKD8F+Nqg8++GDlypVNtlmzZs0Le8zNzc19fX1XrVqlUql8fX2b81FaEqT1fwj6hggQopUKRJHyv6/q+QTI028L/INwoYHmTtqkGifzZNOP3Qn1zVVPMTbBlnjfE+mzNRcFAG2VuR5C6MVfba0C+5iZmZWWlqKDX6LyPKRviq79gZRyhNCUcSOd9DX1hTgmJiY8PFxX7XQGpPV/wdgc0bAohJCqtlJVXlT/+A4hMibLC/l+gzA2RzsxqG+urulNnC+k/sggvk1TvmeJRbnio+xxNtwNAeAN8Hi8HTt2hISE0A1SdOb/m9sTJkwICwt7x4MrFAqZTMYsEwTRuA89LCwsPj6+WcMoWxwkiRcjREaG4z/ldxuA6wkaMtOURXmyO1fk924iLT40yiNQmB2+s48iL4I93hHf8pCy3KP85JrqWjE8uQrA6w0bNuzs2bM+Pj7MSyMjo6VLl8bGxr75Efh8fpOZGBBC+vr606dPt/yfKVOmIIS4XK6+vj5CiMPhBAUFiUQihJBAIHh+dy3A6A6QImbOnOnt7T1t2rR3OUj9gyTpzdPiiDmy1MscWzeOg0dLhfdqdXV1zOWCECqQ0kdy6Z2ZFFM2MsYV9zCEm6uaJZVK9fT0cBwaQLpBkqRSqVQnx8WLF7PZ7IULFzbrILW1tRKJxMrKSgMBtgwrK6s7d+6Ym5u3yNHgYn1TvM7+xh8uxgUihOOSK0dpUln31yGy/Jk2Y7BtVDaSR6DhZ1Weh8gVaVSRTJtRANDGiESi1pzTWxyk9WYT9g01ip6PMIyS1UmunVTVVEivn6Kkta/fs+V4irEl3YisCazNAcQzGe17VBlwgtzykKpVajMKAEBrBLdM3xJGsAxC30cIUZJqxZPHNE1xHTsri/L1fAK0f3N1ZQ/iXCEVl0PPT1L2MceiXPGR9jgH/mUDoGunT5+mafru3bvz58/X2knhT/9d4UJD8cQ5woAwhOHyu9cacu7XP0ypz0jS5s1VLoHC7PDYAUR2OHu8Ix6bSVnvhbLAAOhYdXV1fHx8SEhIbm7uw4cPtXZeSOsthm3laPzBIp57N0RRdZeOkFWlstTLioLH2ozBkIOiXPETg1l/j/5XWeDUchohtHHT5h4DBk+K+ejZM63eEgCgYzI0NGQGsJeWlmqzUCuk9ZbH6+xvOmMZy8ickktqjm9HlEp64zRZWaLNGJg5+a6Fsc4OI8RcNOEvld38Y5//uDap25y4CtOR4VO0GQwA7YBEIsnIyFAPV1c7f/68epmiKLn8X/NestnsI0eOfPPNN/r6+nfv3i0tLdVCqJDWNUjYN9R05nKE4WRZYV38bkomkd6Kp2QvqWGhGUxZ4MfjWe/VJyu6R6DOQWTIN3eSbtcotBkFAG1YRUXF0qVLL168yGKx9uzZ07iWwNGjR21sbJjl0tLScePGrVu3rvG+d+/e7dKli4ODQ35+fpcuXXbu3KlSqTQdMNwy1TwMMxj5MUKIktY2ZN2lJDV6PgHKp9k8794YS0s1AXAMRQT2PPbZfIlTD/zRZUN3f6cDcHMVgNcrKiqKiYnZvXu3qakpQsjV1XXKlCkXL14MDAyUSqV37txRV+g1MzOztLTs37+/et/6+vrVq1dLJBKE0O7duxFCYWFh27dv//jjjzUaM6R17cEFIqMp/0EIKYvypEnnEUKEsTldL+e6+iBM408VhYaGLntSsH3vclcnh1827RGaso/lU7GZ1IzrqmE2eJQrHmit+SAAaGsiIyM///xzJqczvL294+PjAwMDT548OXDgwMYbp6SkrFmzhiRJFouFEOLxeNu3b2+8gYeHh7r0o+ZAO00H2JYOJlN/0Ov6HiWtrTm9U1mcX5+RrHyWq+nzfjpjeuq1vw7Ebre2tjb4983Vr5JU9vvIr5NUj2pg7AxoM8p/X1K27ktE0+VbFpWt/wohVL7527INXyOEyjd9U75xAUKobMP88s0LEUJl678u3/Itoumydf8p/30Joumy3+ZVbPse0XTZb19UbP8BUVTpmi+Uhdnq4yckJGRlZYWGhjY+aVlZGbNw4cKFxrVfSkpKZDLZpUuXjhw58t13370sZpVKVVen2Z5YaK3rEq+TP6+TP0KoIed+7emdpp/9XJ9+m+PYmTDU3hx41oJ/5lxNr6Ljcqmh8Sp+ozlXAWjNxOGzEU0hDBNHfM4MKRZPnPO/hbnMNkaT5zHfho2mfIkwDGGYUeRXCMMRhhlFz/9nIWo+wgmE48bR8wmRkfr4Z86cGTRoUJOTJiQkzJkzByFUXV0tEAjU669cuRIUFBQcHFxVVfXHH3+8LGZTU9Nnz565u7u32E/hOZDWWwVh31Bh31CEUENehjTpvPH738rv3eB19sd5gtfu21I8xZin+F9zrnqKsShXfJIzLoSywKBVIkTi/y0YNV0wMG7Gwv8aUk1aVHV1deo7oow7d+6Ul5ePHTsWIUSS/5qL4/Lly9HR0QihhISEXr16vSxmLpcrlUrf9BO+FeiEaV0MR39iMvUHWqWU37tRG79HVVtJPkqmVS+eyUUTmCdXNwcQRZPZX/ngFwpp673KCRdVJ55QSkprUQDQKnTt2rWgoED9UqFQfPnll3/88QeHw0EI8Xi8xqUSk5OT/fz8EEIHDhwYN27c/fv3X3jMiooKS0tLjYYNab01wnkC4+gFhqOm0nJpw43TssTzyuL8hpz72i8LfDCQyItgh9phv6X/f1lgrcUAgG5FRUWVlpbeuXMHIVRcXDx79uxFixapx7p4eHhkZ//TEV9SUmJsbMzcKS0rKxMKhU0GsKtVVVW1VKXGl4FOmFaNZW4reH+RQF+//lFq7ckdBqM/QQjhAn22uZ3WYhBzUZQrHuWKP5HQ+7LpD6+qcAyFO2FTXHAXEYydAe0Zm80+fvx4fHx8QUGBUChcs2YN005nTJgw4eLFiy4uLgghoVD466+/MuvXr19fXV3t7+///AHlcrmZmZmm6zxDWm8beO7deO7dEELS66fqLh8xm/2LojCbbenQ+PaOptkJscZzrgacIG0EWKQLPskFN+VpLQoAtIogiJCQkBe+1alTp4MHDzIzEgsEAg+Pf+ZgcHNze9nR9u7d+44TP7yJNtMJU15efuXKlXv37j3/1sOHD1NSUlJSUtLS0rQfmJYJ+oZYfLMNFxrUZyRXbF2CKKr+QSLd8OKvexrCzLlaOIm93J9IKac94pRh58jYTEqmvVsAALQK06dPf/PpljIzM8VisZOTk0ZDQm2ltX7nzp2pU6cGBARkZGSYmpru2rWr8bvR0dEeHh48Hk8gEPzyyy+6ClLLDEdNRQjRSoX0Zrws9bLhuJnKJ5lcF2+Ev3RK9ZZFYCjIGguyJuQkcbKAis2k5t5SDbfFxzthw21xArpnQAdgZmY2fvz4N9zY0tLS1dVVo/Ew2kZa9/b2TkpKQghRFGVnZ1dRUWFsbNx4g1WrVpmYaG+sd+uBsTnGHy5CCJEVRTVndul59uB3G6iSVHHsNDgqtgk9FhrviI93xJ/J6LgcekUaNeM6NcYBi3bFDWtzo6bOzMzMHBkWsmH1KuaGEgDtiVAobPEt31Hb6IRRp4OioiKCIAwNDZtsMGrUqNDQ0JMnT2o9tNaCZWxpNnuVflA4WVFUtfeX+vTbyme5ZHmRNmOw4v9TNjJ+KCFgodEXVN4jP74hHlA67fSe6482b96izWAA6LBaXetp1qxZFPWvAdJr167FMAwhJJVKJ06cuG7dOoL4Vz/D9u3bnZ2ds7KyRo8ebW9v7+3trdWIWxmuq4/515sRQtJb8bXn9prOXEnVVrLMbHCBSGsxeIqxH/2J/3ZHBrMe0TF/IANLmWfogevpkVORCJ5sAkDDWkVal0gkFEWJRCKE0OTJk+l/j85mcrpcLh89evRHH30UFhbWZHdPT0+EkJeX14gRI27dutXB07qaoNdQQa+hCKGa6yfl925aLNjakJXGceiszTn5hg8dfCJurrzTUN6lNayZPzvuh7KRAGic9v62pk+f7u3tbWRk1LjqvEqliomJsbOzc3JyCg8Pb2ho6NmzZ69/QwgpFIoJEyaMGTMmKipKve/x48fr6uoUCgVTmb62tvbSpUudOnXS2idqKwxGfGTxzTZE03WXj1Zs+45WKhS5D7TzZNMfm9YtGNF9PHljz7qVfy0YyczJt+UhZbkH5uQDQFO011rv1q3blClTJk2apFQq1St3796dkpJSUFDAZrP79++/efPmWbNmPb/vnTt3ioqKtm7dunXrVoTQwYMHnZyctm7d2rVrVwMDg1GjRsnlcgzDoqKiAgICtPaJ2hgcN/n4O4SQqrq8+ugmrpOXfnCEqraSbemguXPy+fyFC75Wv2Tm5ItyxZ9K6cO59H8SVRUNaKIz9oEb7mYAQ2cAaBkYrd0Gk4uLy2+//TZ8+HDm5eDBg0NCQmbPno0Q2rlz5/r16xMTE1v8pEFBQQUFBeq5BFks1u7du/X02kZ9QolEoqEb6GT+Q/nhtbzACMLeA2EYbqCDoUQZNdifT4k9uQSfQBMdVJMdVWa81tWCl8lkPB5P048FgpchSVKpVKr/WpcvX75jxw4LCwvdRtXi7t27l5mZaWT0+qcLeTwem/2aO1Q67lvPyclRF6h0d3fPzdVIzXFbW1svLy910WSBQGBmZqaJE2mIvr6+Ro7r5S/2ikUISRPP1xz73fTTlQghwtAE19PSMCyEUA991MMG/bcnUzaS6H6G6m6KRbrgYxxaS9lIHMf19PQgretKk7T+9ddfjxkzRrchaYK+vr69vX1LHU3Hab2urk79CxMIBDU1NZo4C5/P9/DwCAoK0sTB2wFBj2BBj2CEUM2J7dLb5ywXx5IlT1iWDhihpcuDKRsZYEGs6U2cL6R2ZdKf3VAOscEjXbFhNjgLMir4H319faZKIngFHf/FmJqaVldXM8taKGwGXs0g7AOr/+7H2Jya0zvLfpuHKEr5NEsnZSNzI9hB1tiKNMpyL5SNBKB5dJzWfXx8mMdHEUKJiYk+Pj66jQcwTKb+YDZnNSWrq9z3S1XcWrpBTpY+1WYARlw01QO/FsZKHsVy0sc+vKrqfIhckqrKroX8DsBraK8T5tKlS+Xl5RKJ5Nq1a1KpNDg42NDQcMaMGSNGjAgICOByuatWrWoynSvQLVxoYP7lBoSQ4snjiu3f6wdH8H0CaJom9MVai8H+f2UjU8rp2EyqL5SNBOB1tJrWHz58+N5772VlZWVlZXXv3t3Q0LBv376bNm366aefVCrVypUrhw0bprV4wJvj2LlZLt6FaEp2J6H68AaTj5cQhqY4Xx/jam80kZ8J5mdC/NKLuPSMjs2kvr+j7GOOjXfExzni/FbxUB0ArYW2BzjqxMyZM729vbVQ5lgT6urqNDUS5q3QKhLDidrTsZLrJy2+2UbJJSwjc62VjVSTk4gpG3mzhH6+bCRFURiGMc8nvyOpVAojYXSoyUgY8CbgYgXNgxEshGGikGjL7/fiAlHNn1tKVk5HCJFlhdoMgykbeWIw6+5Ylp8JtiKNcthPzr6pSi2nF373X77IUGBgtHbDJm2GBEArAV9fwVvCWGyEkPFHS2ilgm6Ql2/9jmPjLJ70BVVTSRhp77EAKz422wub7YXfr6L3ZFGhf6SVrN9G/fcxUsi/Xthr/JhR7e/RFQBeDVrr4F1hbA7G1bOYv0UcMYcsfVq6enbdX4eoehklrdVmGF5ibJk/scOnjG/lioTGyMhGoW+xNamoTvn6fQFoT6C1DloMxuawLR0sl+yhaar+/q2q/auNor5m2zjjPIHWykb27t1LUJml3D8La5CIeeihXif7/cq+UDYSdCTNSOsqlaq2tlYgEDSeexuApnAcQ7ieTwDPsyeGE7Xn90uuHDGftwHjcHCBAWqJ25ivoK+vfzf51q5duzgcTlTUOgMDXrUCHc+ntjykpl9TDbfFo1zxQGsNBwGATjVjJMytW7d69+6dkZGhnmC7rYCRMLpFySW4nrDij/8qn2ZbLNyuklRrc+S7GlM2ckcmVdmcspEwEka3YCTMW2hGa10sFiOEeDx4CAQ0D1M7zPj9hZSsDtFU2W/z2Jb2xu9/S0lrcaGB1sKwEfxzczW9it6VRQ04RYo5WJQrHuOGm0PSAO1IM9og7u7u3bt337dvn+aiAe0bztdHOGGxYKs4Yg5ZUVy8fGrt2T00qaTqZdoMw1OMLfcnnk5kbw4gcupojzhl8BkyNpOSwM1V0C40r2994sSJixcvvn///qBBg5g56hj6+vpDhw7VQHigPcIwnK+P8/Utv9tDKxUNj/+u3LVCPHEOz70bYrGhbCQA76gZfev19fUv6+Hy8PDIyMhouahaGPStt3J0gxzhhCThWN1fcaaf/UyIjHA9gaZvrj6vsgEdyqViM6ncOjTOERvviAdYYNC3rlvQt/4WmtEy4nK52dnZL3zrtbN1APAKTG0Z/cDxgp6Dcb5+5e6ViiePzL/aTCsVOF97c3owZSOneuD5Enp/Nv3BVRULQ6NtiRh32tVQa6oFKX4AACAASURBVFEA8K6akdYxDHNyctJcKAAwd1CNor5W1VZhbE7JzzNZJlYmH39HKxowDldrYTQuG7ntAdnvFGUjoKFsJGgr3qYfMyUl5fbt23K53N7evmvXrs7Ozi0eFujgCJEYIWTx9RZVdZmqtrJk+SfCvqGi4VG0imSKFmiHnwnm4ada05dzpRiDspGgrWjetVlRUTFhwoS//vpLvQbDsOjo6C1btkA/DGh5GEaIzRBClt/toeSShux7Fdu+Nxw3U8+3H4bhSFv93QSGgqyxIGtCThJM2ci5t/55smmQFYZjCCGkVCrz8/NtbGxgBDDQueb9YURGRiYlJS1fvvz+/fvFxcWpqalz587dtWvXwoULNRQfAAghjM0hREZcly4WC7fzOvlLr58qWjxZUZBJKxXaDKNJ2cglqSr7/eTsm6qTqdm2Lh5+A0Ms7Jxu3rypzZAAeF4zRsI8ffrU1tb20KFDY8eObbx+yZIl69atKysra5Hy1poAI2HaH1VVKS40rD60viE33ezzXzEOT0OdM68eCcOUjVw7f7rM0J4e+h+UerT7vc1JVy9qIpKOCUbCvIVmtNYLCgoQQkOGDGmyftiwYRUVFRKJpCXjAuCVCLEZxuaIJ84x+fg7XE9YunpO+aZvEEI0qdVnipiykcGm9bTIAiGERGYZpdKdmRSUjQQ61Iy0bmJighBKS0trsj4tLY3P5wsEgpaMC4A3wzK1RhhmPneN4dgZVL2saPHkmuNbEUK0itRaDF/MnCo4tkAUGyPYHjl1xqcXC2n7/cqwc2RcLqWgtBYFAP9oxi1TV1dXX1/fKVOmbN68OTg4GMMwpVIZFxf35Zdfjh07Fp7XALqEEyxTa4SQ5eJYlaRGUZBZvnGBwciPBT0Ha+HkAQEB6amJ169f9/FZ4OnpiRCqVhBQNhLoSvPmMs3IyBg6dOiTJ084HI6xsXFpaalKperevXt8fLyxsbHmonxH0LfeAalqqxClqn98p/bUDuOYBRxHz7c7zjs+Zfp2ZSOBGvStv4XmDXDs1KlTenr6wYMHk5KS6urqHBwc/Pz8wsLCWCwYxAtaF2bku6BHMNexM65vWH1kY0Nmmsm0pYSBVtsfUDYSaF8zWuvp6enR0dF79uxxd3fXaEwtDlrrANG04mkW28qxfOMCjGCZTFv6hjVnWrYmDEWjGyX0rizqYA7V3RSLdMHHOOBCeOTj5aC1/haacbEqlcqUlJTGhRsBaDMwjGPrihEsk+k/ikJiEE0VLZ5cfXSzlqNgykZuDiCKJrNneeInn9DWe5UTLqpOPKFIuLkKWkgz0rqnp6eFhUVCQoLmogFA0zCCxbFzQzhhPn+LoOdgsqzw2cJw6Y3TWg6DR6AwO/xgIJEbwQ6yxlakUZZ7lZ9cU10rbsa9LgBeqBl94iwWa/369XPmzKmpqRk4cCAz3pFBEAS04kHbgvMEuJUjQsh83jqaJOVp12qO/S6eOIfr6qvNMJqUjfzwqorA0AQnLNIFdxbBzVXwNqDeemsHfetaoyzOxwUi6dXj9RnJRjHfsIwtkC7mMk0pp2MzqQM5lI0Ag7KR0Lf+FprRWmez2Zs3v7gv0tAQylGDNo9tYY8QEg2L4nr4EfqGFTuW0g31vPC5SLs5xc8E8zMhfulFXHpGQ9lI8BaacZlUVVUhhEJDQ62srDQWDwC6hmFcJy+EkNHkLxW56SSLXfbTDK6rj+GY6dqM4mVlI8c7YcNtcQJDCKG0tLSysrKAgAAoGwkaa0Zaz8rK+uSTT4qKijQXDQCtB8bmcN26klKp8cwVqtICVW1l2eo5wkHjhAFh2gyDKRs53hF/JqPjcugVadSM69QYB6x094KTfx5mmdiK5CV3k24aGBhoMyrQmjWjx9DZ2ZnFYj179kxz0QDQCuF8fa6TFyEyMpm5nOvQqSHz7+JlU+sfpmg5DCs+NtsLvxbGih9KcKj6A39sksxPrP7sXLlZ18OHD2s5GNCaNSOtm5qafvrpp//5z3+Y3hgAOhqWsSXbxoXr6ms05UuWqZXkytHS1XPI0qdaDsNTjC31Z3NwhEgFQkgul/3yANuZSdVC2UiAEGpWJ4xSqXz69Ont27cdHBx8fHwsLCzUb1lbW//6668aCA+A1ohj64oQEvYbybZ0xAWiqv2rKUm1UeRXzFzb2giAw1mwYMGK/3ZlGZhaG+gt+Xj88UJ67i1lH3MsyhUfaY9zoPJeB9aMtE7TdG5uLlM5QCaT5eTkqN+iKHhCDnQ8OM5180UIGY6Z1pB1F7HYZb99wbZyMhw3UwsnX7Tgq5gpE8vLy318fAiCGOf+/2Ujp11ThUDZyA6sGWmdw+EkJydrLhQA2iiMw+N17oEQMnp/ofLJY6peVrZmrqD3MOF7IzV6Xjs7Ozs7O/VLQw6KcsWjXHGmbOR/ElUVUDayQ4KvagC0GEJfzPPsifP4xu8vZFs5Kgoel676rP5BopbDsBFgs73w1NGs00MIhNCAU6TnIXJFGlUi13IgQDeandavXbsWFRXVo0ePiRMnMmuWL19+6dKllg4MgDaMZWbDdenCsXUzHDONMDCRJV0sW/+VsihPy2F4irHl/sTTiezNAUROHe0Rpww+Q8ZmUhK4udquNe+ptZ07d37wwQeenp5CoVA9gP3p06eXL18eOHCgBsIDoG1jpu9gWdjhfCGuJ6g5+QdZ9kw84TNcoL0aSkzZyAALYk1v4nwhtSuT/uyGcogNHumKDbPBWfCNvd1pxq+0vr7+888/nzVrVlpa2kcffaReHxwcfPPmTQ3EBkA7gREsnmdPwtBUPyhczzcAY7Ertv9QfWQjas7cZO8OykZ2EM1I6/fv36+url6wYAH277vrtra2tbW1Uqm0pWMDoL3BeXx+1/4YV89w7AyOvQdNqUp/mSVJOKblMJiykdfCWMmjWE762IdXVZ0PkUtSVdm1kN/bg2akdabWI0EQTdYXFxcTBMHhcFoyLgDaNcLAmO83ECNY4olzWSaWZFlh6eo59em3tByGvRD7ygd/NJ61awBR1YD6niC7/0muuU+V1Ws5ENCSmpHWO3furKend+DAAYRQ4wb7li1b/Pz82GyYuQuAZmNbOvA692CZWhuERGNcQX367fIt3yoLs7Uchp8JtqY3UTiJvdyfSCmnPeKUYefI2ExKRmo5ENACmnHLVCAQzJs3b/bs2bm5uSRJ1tXVHTly5Pfffz979uyxY9r+FglAO8NM30GrSJpUIIJdd/mIsiDTYNRUQl+stRheWzayqKjo2LFjFhYWI0aM0GYNetAszRsJs2TJEoTQzz//LJfLEUJjx441MjLaunVrWJhma9qVlpZ+++23zPLAgQMjIiIav5uZmbls2TKKor788ktPT0+NRgKARmEES8+nH0KIMDCS8wQYhlfFrcVYHMORHyG8af+n5rywbORgfuHBj3pT3sOJ4odBBw4f2bdLa/GAZmleWsdx/Pvvv583b15iYmJVVZWZmZm/vz+fz9dQcGq1tbV3795dt24dQsjMzKzxWyqVauTIkb///juPxxs7duzdu3ehlx+0A7ieUNBrCEJIFBQuf5CIECrf9A2vcw9NP7nahBUfm+2FzfbCH1TTX644IfUaQU9ai0jFqXlWcrkcJi1qnd5mthWRSBQUFNTiobwa03fv4OBgbGzceP3NmzednZ379u2LEOratetff/01dOhQLccGgOYQYjNh31CEkEHoB8qiXFVtZeXun4TvjdTz6qXNMDobYpE+JglnsiUUicpzlYgVfJ6IcqMinHER3FZrZVpd75jy30iSRAix2WwDA4PNmzf37dt306ZNjbcvKChQl8Wwt7d/8uSJDoIGQPPYNs58/yBCZKQ/aByG4w059yu2fad8mqW1AMaNGxdgb8CbZ83/uf+OLRu+8mVdK6Yd9yvDzpFxuZQCyv21Gq1rbsTHjx9HR0c3XuPu7r5jxw57e/sTJ04ghCoqKjw8PD7++GP1OEs2m82kfoSQUqmEHhjQ7vE8/BBCNKnU832PpijprbOK3HTR8GjCwPi1+74LFot15tihmpoaPp/PfHsOs/v/spHTr6mGQ9nI1kF7af3p06cpKSklJSWRkZGNu+SSk5OPHz9uYGAQHR3t5ub26gdWDQ0NMQxTKBTqI7i4uKxdu5ZZfvTo0ZgxYzT3EQBoPTAWm+83ECHEMragFfVIpao9E0srFaLh0RhLg90iTWbXg7KRrZCWOmGys7M7d+68YsWKTz75pK6uTr3+3LlzgwcP1tPTu3fvXo8ePRq/1djZs2c3b958/PjxqKiowMBAPT09uVzev39/hJCvry+O47/++uvGjRsrKyv79OmjnU8EQCuBC0TC90YSRmb8HsG4wABRVOWen7T/5CqUjWw9MForVSmYeTZqamqMjIxKSkrUo1kGDBgwfvz4mTNnIoTee++9SZMmTZs27fnd8/Pzjx8/XllZ6eXlNWrUKIIgSJLct29fZGQkQqiurm7Pnj0URU2ePPmFE/V+9NFHZmZmY8eOZV5yuVwvLy8NfdIWV1dXp6+vr+soOi6pVKqnp9e2xmgrnjxS5KTzew2p2rOK33Owlm+uMiga3Sihd2VRB3Oo7qZYpAs+xgEXNv9bBEmSSqUShtw0i5bSOqOqqqpxWidJksvlPnr0yMXFBSG0dOnSe/fu7d+/v8XPO3DgwOzsbBMTE+Ylm80+ceKEFsZltgiJRCIUCnUdRcclk8l4PF7bSuv/oGnycQolkxKWDg0JRzl9w1jWztqPol6FLpXg+3KJC0VYkCU90VEVbEG9edlISOtNcLnc1z7Sr8tbpqWlpRRFmZqaMi/Nzc3Pnj2riRN17tw5PDz8hd8DWj+apiGt6xCGYW2utf7//AYghGilQubpTyAVnXu3/mGKKHgiYWiitRCECI03QOPdUGUDOpRLrXlEzEikxzjgkS54gMXrO98hrb8FXaZ1FouFGs2DSpIkFJYBoMVhbI6g5xCEECWtJSuKKblE9ncCLavTD4rAOFythcGUjZzqgedL6P3Z9AdXVSwMTXDCIl1wZxHcXG1JumyDmJiYcDicwsJC5mVhYaGVlZUO4wGgfcMFIv0BY9iWDnqePRHBppUN1Uc2Sq6d0HIYTNnIx1A2UmN0mdZxHA8JCTl06BBCSKlUHjt2TNO1ZQAACCGWqbVoyCRcIOJ3G0DVVNIqsnLvKu2XBYaykRqivU6YPn36MAXCBg8ezOFwEhMTEUKLFi0KDg7OysrKyckRi8WjRo3SWjwAAI5DJ45DJ4SQnlcvsryIrCypO7dX0DeUY+uqtRheXTZSpVLt3r07LT1jWHBgSEiI1qJq07SX1tevX6/uRleXa/f19X3w4MGlS5cMDQ0HDRrE9LYDALRMr0tfhBBNKtm2rqryooYGufzeTf2BYwhDU+3F8KKykWZnf8i4fVnuM3rHtLlbf5JERIRrLZ62S6sDHHVl5syZ3t7ebXQkDIxb1622OG793VHSWsmNUzyXLmTZM7KqVH/AGIyrg7Eo96voAP+uNVO2IxtvlBQXWHr8wtGWHwDd/nSsixUA8CZwgUgUPJHj6Mlx7Ewr6ilpXe25fdLrp7Q8p7aXGOvu4UjcP4PkNcSDczcox+Az5M5Mqk6pzSjaHkjrAICXYplaG4R9SBiZcV26KMsKaRVZfWhd/YNErQWwafVPXZ6d0/vGNchI8mTnglme+MVC2h7KRr4S9GUDAF6P6+TJdfJECHGduygKMjl27rVndvF7DdH0zVUHB4frF+PVjyOFGULZyNeDtA4AaAa9ru/pIUSrSJaZjeLJIwzHZX9fFfYZTojNXr9zC4Gyka8GaR0A0GwYwRL2H4UQoqS1GIYpi/KUz3KUxU+EAWHavLlqI/hnTr70KnpXFjXgFCnmYFGueIwbbt6Byw1AWgcAvD1cIBINj0YIkWWF9Q9TycqShpz7GIst8A9GWhw+5CnGlvsTP3YnmLKRHnHKdykb2dbBLVMAQAtgmVobjp3BtnRgWzooch9QCnntmV31GcnajAHHUIAFtjmAKJrMnuWJn3xCW+9VTrioOvGEIjvSzVVI6wCAlsR18hJHzMF5Apa5bX36bVpRX3Nim/JptjZj4BEozA4/GEjkRrCDrLEVaZTlXuUn11TXitv/YzoIOmEAABrC7zaA320AolS4wKD+YTLG05MlXRT0HqrNJ1eblI388KqK6ABlIyGtAwA0CSf0B41DCFGyOppUNOSkEyIjZWE2v+cQnKe9qWzshdhXPthXPnhKOR2bSfU9QdoIsEgXfJILbsrTWhRaAmkdAKANOF/fIOxDhBBZVihLvkgW5amqy2lSyfcbqM2bq34mmJ8J8Usv4tIzOjaT+v6Oso85Nt4RH+eI89tLOoS+dQCAVrFMrcURcziOnXGRuP5BokpSLUk4Vv8oVZsxMGUjYwcQTyeyo1zxuFzKbp8y6rLqxBNKRaO0tLRR4VNCx028ffu2NqNqKZDWAQC6wXX2NoqeT4iMcIG+7NZZWkXWntunLMzRZgxM2cgTg1l3x7L8TLAVaZT9jqqeg4Yf5/Q+pR8YHDqqtLRUm/G0CEjrAAAd4/sNMoqej+EEhuOy5Iuq2qq6iwdV1eXajMGKj832wq+FsX6xz8KMbOj+U1HfaNLa58KtO9oMo0W0l84kAP6vvXsPauJO4AD+27wwCOH9sgiogFosaEGPm6MncBXU1gp3woFYnUGojgWrFa6MznVaRq/XWjrglA5q7VVPWxXFQeDE99RaweOitIJGRV7DKwYCgYSEPDb3Rzo57uZqK5hs2Hw/f+W3w+x+J5N8Z/lt9rcw1VGU88t/JITQaqVhWK5uqnOY9by6tZn3Qgyx4iOqE14M4Q91am6fJXyBvvP7bZ3zPqvSp8zirAvheFjvya+TgloHANvCETq5Jm8mhOhl3bqO+0TkoekV0JpRx4gYwuFa+uiurq5ny7/Of3e3Xq//y7Ev4hMCL3TT5a1G08XV9SGc1YEcgW1Pc6DWAcBG8byec0nbrtPpqJ4WdV2twD9Y23aX6+7jEBxu0ePGxsY2XIs1D1cFcFYFEIWWW9lBH3lIb/nOsMLfppeNRK0DgK1zmPOCw5wXCCHazvsjl04IAuep6msdQiL4voFWy+Dy38tGvtNgkF0jyUFU5lxOhLtt1btt/y8BADCOY2S85+Y9FF9Aq1Ujl0/SmlHl9SrDsNyaGfynU28t4IiTeOeWc90cSPJFQ9gp/Yff031qa6Z4EpytA8DUI0pIJ4QYx9T63g6VUuEY9bux1ibHiBhrLgsc5kaFuXHfXURuSI3lbXT4aV2oC7U+hJM+h+PM6LKRqHUAmKooB6FrSg4hRD/Qq7nXQHG5PI8ZhhG5MCzaaneumpaNjPHlfriEe7Gb/vtDY/5N3W/9qPUhnKRADp+JCRHUOgBMeTwPP48NOwkhY63Nqu9quC6e9Iic4+gsmBVmtQymZSNXBZDBMW5V54/P5PvDLM7rwZzf+Fr14irm1gGAPRxmh3lu3iMICKU1o4rqL+lR5WjDZf3jLmtmcHMg60M4F1fwbv+e97wrlVtnmHVcX9BgeKCw0rLAqHUAYCHHyHiv3L0cRye9vG+o8qBRr1PV1xpGBq2ZYeZ06q0FnNvJvJpELiEktkZvurgqtfDFVUzCAACbiRIzCCFGnVbbfk/X1eKcmKFtuTMtbAklsN6CvFZ+Jh/O1gGA/Si+wC1tu+uaHKLTjYqvqOrP6/o6Nff+RWjrPQ3vCc/k02j12W9u9ZsVGrt8VXt7+yQPhLN1ALAjXHdvj6z3CCHatrsjF49TXC7hCSgeXxAQarUM5ourj9XcE630nkZ63elPRxsf6rPPPhafyti4+bvLtZPZP87WAcAeCWY977X1Y4fQRYYh2eDJEsNQv6a5Xj/Qa80M3kKSG8apf433MueBPnwV8Z5D/3r9vbvNk9wtah0A7Jrji7E+eaVcV09t1yP5kQ+J0Th6+xtaNWzNDK+vTnD6Zh+pOyos356QkDDJvWESBgCAEEJEiRmixAyjQa+526BpuumakjP2oHHa/CiKL7D0oZOSkg7pdCfOVC9avfjtbVsnuTfUOgDAf1BcnntGHiHEoBhQ1f1D19fhGBVv6O91CIkglrypKDUlJTUl5ZnsCrUOAPB/cF08PDftJoRoOySK6i+c41N4PjMJTfOfm810tJ+BuXUAgCcRBM7zfnufcOFLur6OgS/36KSdY63NhiEZ07l+EmodAOAXcVy01HfXIb5PwFjLD7LSAqNBr5GIabWS6Vz/C7UOAPB0RAnpvrsOURzuaMOlgb/tMep1GonYaNAznetHmFsHAJgQinJ//R1CCK1UjFwu19z9pyhhrV7WLQiab9GLqz8LtQ4AMCkcJxevN/9KCNF1tQye+tTppdccgsONBj3fJ4CRPKh1AIBng+8f7JP/GSFE/f11xdnP3dcXUFw+R+TOFblZMwbm1gEAnjFhRIzvn78UBM7T3Gt4/PEWWq3Udt43jlnpaaeodQAAS3FeluZX+DVH6KS8Vvl43w5iNGrb7xHaYNGDYhIGAMDi3Nf9iRBCa0YVZz/nunm7pW7Vy3r4/nMscSzUOgCAlXCmOXptLSKE6Po6Bo58MP1XiY6LXzbqNDwPv2d4FNQ6AIC18X0DfXd+TghRN9UPnSxxW5vH8/TjCJ0400WT3znm1m1deHi4VqtlOoX92rRp0+XLl5lOYb+OHz++c+dOplNYkHBBtF/h19PmRap/uNH3QfYzuWcVZ+u2TiaT6fV6gcDiS4PC/6VUKkdGRphOYb+USqVSaXN351uCc/wa5/g1z2RXOFsHAGAV1DoAAKvYxSSMQqG4detWeXk500EmwmAwVFRUODg4MB3ETvX09Ny4ccNgsOwPjeGniMXitra2KfrltYTQ0NCIiIgn/41d1Lqvr29jY+MU/WQEBQVVVlZSjK4cZM9omm5ubu7s7GQ6iJ0aHBxUqVRT9MtrCUuXLv3ZWqeMRqN10gAAgBVgbh0AgFVQ6wAArIJaB3g6w8PD/f39TKewXyqVqqurCxexnwC1PpUolcro6Ghvb+/8/Hyms9ipV1999ZVXXklNTY2NjcVtStZXUlKycuXKnJycsLCw+vp6puPYKFwynUr0ev3Q0FBlZaVEItm7dy/TceyRVCr18fEhhGRnZy9ZsiQ7O5vpRHbq1KlTFRUVX331FdNBbJFd/MCRNXg8nqenJ9Mp7Jqp0wkhw8PDHh4ezIaxT3K5vK+v79y5czExMUxnsVGodYCnduLECZlMlpSUxHQQe3T+/Pny8vL29vYtW7YwncVGYW4d4OlUV1eXlJScOXOGw8HXhwHp6ekVFRX79+/Pzc1lOouNwufSJhiNxt27dycnJ0dFRTU1NY3fvmvXLn9//6CgoKKiIgYTsl5ZWVlaWtrixYtra2vHby8tLZ09e/aMGTPy8vIMBsOFCxd2795dVVXl4uLCVFRWqqqqyszMjI6OLi4uHr+9trY2PDzcy8srLS1tcHDQvEi1Wq3m8/lMJJ0CMAljK6RS6erVq3NyclQqlXnj4cOHT58+XVdXp1ar4+LiwsLCPvroo9bWVo1Gc/Xq1fPnz2N691np7u6OjY29deuWXC43b7x+/fp777135coVLy+vFStWlJWVFRYWikSixMREQkhWVtbmzZuZi8wqXV1dCxYs6Onp6erqMm/s7+9PTU09evRoXFxcVlZWXl4eRVEtLS1CobCrq+vQoUMMBrZpRrAlLi4u9fX15mFMTExZWZnp9fvvv79mzRqGctmLyMjIY8eOmYcbNmzIy8szvT527NjChQsZymUv3njjjR07dpiH+/bti4+PN72WSCRCoVCpVCoUiv7+foYCTg2YhLFpEonEvKxPRESERCJhNo+9wfvPLIlEsnDhQtPruXPn0jTd0dEhEonwT+qTodZtl9FoHBwcdHZ2Ng1FIhFubrQyuVw+/v3XaDTjp8jA0uRyuZOTk3koEokGBgYYzDNVoNZtF0VRbm5uw8PDpqFCofDy8mI2kr1xd3cf//4LhcLp06czG8muuLu7j7+Vd3h4GPdt/BKodZsWGhpq/mFMU1NTSEgIs3nsDd5/ZoWGht65c8f0+uHDhxRFzZw5k9lIUwJq3VY0NzeLxWKDwSCRSMRisemHXJmZmcXFxVKptK2t7cCBA5mZmUzHZK2WlhaxWDw6OtrW1iYWi00niRs3bjx8+PD9+/cHBgaKiorw/ltOT0+PWCyWyWRSqVQsFkulUkLI2rVrb968WVtbq1arCwsLU1NTx8/JwE/BmjC2IiUlpa2tzTysqanx8fGhabqgoODIkSNcLjc3N7egoIDBhOy2ffv2b7/91jwsKyuLiooihBQXF3/yyScajSY9Pb2oqIjHw2+CLWL//v0HDx40D7dt27Zu3TpCSHV1dX5+vlQqjYuLO3DgAC6W/hKodQAAVsEkDAAAq6DWAQBYBbUOAMAqqHUAAFZBrQMAsApqHQCAVVDrAACsgloHAGAV3DIHMHGPHj26cOHC8uXLg4KCqqurZTJZUFBQfHw807nAruFsHWDiTp48GRsbm5ycXFpaGhkZmZmZWVBQIBaLmc4Fdg21DjBBd+7cWbRoUXd3d29vb0ZGxowZMwgh06ZNe/DgAdPRwK6h1gEmKDg4eNmyZdeuXUtKSnJzcyOEaLXaxsZG8wN9ABiBWgeYIKFQyOVyr169GhcXZ9py6dIlHx+f+fPnMxsM7BxqHWDiRkdHGxoali5dahqWl5enp6cbjcaqqipmg4E9Q60DTFxdXV1gYKCfn595uHLlypqamoCAAGaDgT1DrQNMXG9vr+lpDyZZWVn19fU6nS4iIoLBVGDn8BgNAABWwdk6AACroNYBAFgFtQ4AwCqodQAAVkGtAwCwCmodAIBVUOsAAKyCWgcAYBXUOgAAq6DWAQBYBbUOAMAqqHUAAFZBrQMAsApqHQCAVVDrAACsgloHAGAV1DoAAKug1gEAWAW1p/IY7wAAAB9JREFUDgDAKqh1AABWQa0DALAKah0AgFVQ6wAArPJvJi1znvbRhnEAAAAASUVORK5CYII=",
"image/svg+xml": [
"\n",
"\n"
],
"text/html": [
""
]
},
"execution_count": 32,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"\"\"\"\n",
" fem(n, α, β, f; a=0, b=1)\n",
"\n",
"Approximate the solution to \n",
" -(α u′)′ + β u = f on (a,b)\n",
" u(a) = 0, u(b) = 0\n",
"using finite element method with hat functions\n",
"\"\"\"\n",
"function fem(n, α, β, f; a=0, b=1)\n",
" h = (b-a) / n\n",
" x = a .+ h .* (0:n)\n",
"\n",
" α_av = (α.(x[1:n])+α.(x[2:n+1]))/2\n",
" β_av = (β.(x[1:n])+β.(x[2:n+1]))/2\n",
" f_av = (f.(x[1:n])+f.(x[2:n+1]))/2\n",
"\n",
" K_d = [ (α_av[i]+α_av[i+1])/h for i ∈ 1:n-1 ] \n",
" K_dl = [ -α_av[i+1]/h for i ∈ 1:n-2 ] \n",
" K = SymTridiagonal(K_d, K_dl)\n",
"\n",
" M_d = [ h/3 *(β_av[i] + β_av[i+1]) for i ∈ 1:n-1 ]\n",
" M_dl = [ h/6 *β_av[i+1] for i ∈ 1:n-2 ]\n",
" M = SymTridiagonal(M_d, M_dl)\n",
" \n",
" F = [ (h/2) * ( f_av[i] + f_av[i+1] ) for i = 1:n-1 ]\n",
" \n",
" u = (K+M) \\ F\n",
" u = [0; u; 0] \n",
" return x, u\n",
"end\n",
"\n",
"a, b = 0, 2\n",
"α(x) = 1\n",
"β(x) = 1\n",
"f(x) = -8 + 16x^2 - x^4\n",
"u_e(x) = x^2*(4-x^2)\n",
"\n",
"x, u = fem(5, α, β, f; a, b)\n",
"x2, u2 = fem(10, α, β, f; a, b)\n",
"\n",
"P = plot(u_e, a, b; \n",
" title=L\"FEM solution to $-u''+u=f$\", \n",
" label=\"exact solution\", linestyle=:dot, color=:black)\n",
"plot!(x, u; \n",
" label=L\"n=5\", markersize=2, m=2)\n",
"plot!(x2, u2; \n",
" label=L\"n=10\", markersize=2, m=2)\n",
"\n",
"nn = [2^j for j∈2:12]\n",
"err_fem = zeros(length(nn))\n",
"for (i,n) ∈ enumerate(nn)\n",
" x, u = fem(n, α, β, f; a, b)\n",
" err_fem[i] = abs.( maximum(u - u_e.(x)) ) \n",
"end\n",
"\n",
"Q = plot( nn, err_fem;\n",
" title=\"error curves\",\n",
" xaxis=:log,yaxis=:log, m=2, xlabel=L\"n\", ylabel=\"error\",\n",
" label=\"FEM\" )\n",
"plot!( nn, nn.^(-2); linestyle=:dot, label=L\"O(h^2)\")\n",
"\n",
"plot( P, Q, layout=(2,1), size=(500,500) )"
]
},
{
"cell_type": "markdown",
"id": "a5464f54",
"metadata": {},
"source": [
"Here, we saw that this method converges with rate $O(h^2)$ in the max norm as the step size goes to zero.\n",
"\n",
"Today, we are going to see some of the ideas needed to prove that these methods converge. It will be easier to do this in the so-called *energy norm* instead of the max-norm: "
]
},
{
"cell_type": "markdown",
"id": "a66ec364",
"metadata": {},
"source": [
"## Error estimates\n",
"\n",
"This section is mainly based on Chapter 14.4 of @Suli: Theorem 14.6. \n",
"\n",
"Recall, we are approximating the solution to the weak formulation: find $u\\in V$ such that \n",
"\n",
"\\begin{align}\n",
" \\mathcal A(u,v) = (f,v) \\qquad \\text{for all } v\\in V\n",
" \\tag{weak}\n",
"\\end{align}\n",
"\n",
"with the approximation (projection): find $u_h \\in V_h$ such that\n",
"\n",
"\\begin{align}\n",
" \\mathcal A(u_h,v_h) = (f,v_h) \\qquad \\text{for all } v_h\\in V_h\n",
" \\tag{G$_h$}\n",
"\\end{align}\n",
"\n",
"(we now put subscript $h$ to distinguish $u$ and $u_h$ and $(\\text{G}_h)$ is for \"Galerkin\" with mesh size $h$) where\n",
"\n",
"\\begin{align}\n",
" \\mathcal A(u,v) &:= \\int_a^b \\left[ \n",
" \\alpha(x) u'(x) v'(x) + \\beta(x) u(x) v(x)\n",
" \\right] \\mathrm{d}x \\nonumber\\\\\n",
" %\n",
" (f,v) &:= \\int_a^b f v,\n",
"\\end{align}\n",
"\n",
"and $\\alpha(x) \\geq c > 0$ and $\\beta(x) \\geq 0$.\n",
"\n",
"Under these assumptions on $\\alpha, \\beta$, we can define a *norm*: \n",
"\n",
"::: {#thm-}\n",
"\n",
"The function \n",
"\n",
"\\begin{align}\n",
" \\| u \\|_{\\mathcal A} := \\sqrt{\\mathcal A(u,u)}\n",
"\\end{align}\n",
"\n",
"defines a norm on $V$ (recall $V$ is a space of functions with zero boundary conditions).\n",
"\n",
":::\n",
"\n",
"**\"Proof\".** \n",
"\n",
"*Positivity*: Since $\\alpha,\\beta\\geq0$, we have \n",
"\n",
"\\begin{align}\n",
" \\mathcal A(u,u) \n",
" %\n",
" = \\int_a^b \\left[ \\alpha (u')^2 + \\beta u^2\\right]\n",
" %\n",
" \\geq 0.\n",
"\\end{align}\n",
"\n",
"*Definiteness*: Let's do this for smooth $u$ (as mentioned before, $V$ is a bigger space but I don't want to introduce Sobolev spaces here). If $\\mathcal A(u,u) = 0$ then $\\int_a^b \\alpha(x) (u'(x))^2\\mathrm{d}x = 0$, then (since $\\alpha(x) \\geq c>0$), we have $u' = 0$ on $[a,b]$. Since $u\\in V$, we have $u(a) = 0$ and so \n",
"\n",
"\\begin{align}\n",
" u(x) = u(a) + \\int_a^{x} u'(t) \\mathrm{d}t\n",
" %\n",
" = 0.\n",
"\\end{align}\n",
"\n",
"That is, $\\|u\\|_{\\mathcal A} = 0$ implies $u = 0$. \n",
"\n",
"*Homogeneity:* Since $(\\lambda u)'(x) = \\lambda u'(x)$ we have \n",
"\n",
"\\begin{align}\n",
" \\| \\lambda u \\|_{\\mathcal A} \n",
" %\n",
" &= \\sqrt{\\mathcal A(u, u)} \\nonumber\\\\\n",
" %\n",
" &= \\sqrt{ \\int_a^b \\left[\\alpha(x) (\\lambda u'(x))^2 + \\beta(x) (\\lambda u(x))^2 \\right]} \\nonumber\\\\\n",
" %\n",
" &= |\\lambda| \\|u\\|_{\\mathcal A}.\n",
"\\end{align}\n",
"\n",
"*Triangle inequality:* This follows from the fact that $\\mathcal A$ is an inner product as we will see: since $\\mathcal A$ is symmetric ($\\mathcal A(u,v) = \\mathcal A(v,u)$), positive ($\\mathcal A(u,u) \\geq0$), and linear in each component (bilinear), we have\n",
"\n",
"\\begin{align}\n",
" 0 \\leq \\mathcal A(u+\\lambda v,u+\\lambda v) &= \\mathcal A(u,u) + 2\\lambda \\mathcal A(u,v) + \\lambda^2\\mathcal A(v,v) \\nonumber\\\\\n",
" %\n",
" &= \\|u\\|_{\\mathcal A}^2 + 2\\lambda\\mathcal A(u,v) + \\lambda^2 \\|v\\|_{\\mathcal A}^2\n",
"\\end{align}\n",
"\n",
"This is a quadratic in $\\lambda$ with at most one real root (since it is positive), and so the discriminant is less than or equal to zero: $D := 4\\mathcal A(u,v)^2 - 4 \\|u\\|_{\\mathcal A}^2 \\|v\\|_{\\mathcal A}^2 \\leq 0$. That is, we have the Cauchy-Schwarz inequality $\\mathcal A(u,v) \\leq \\|u\\|_{\\mathcal A} \\|v\\|_{\\mathcal A}$. We therefore have (taking $\\lambda=1$)\n",
"\n",
"\\begin{align}\n",
" \\mathcal A(u+v,u+v) &= \\|u\\|_{\\mathcal A}^2 + 2\\mathcal A(u,v) + \\|v\\|_{\\mathcal A}^2 \\nonumber\\\\\n",
" %\n",
" &\\leq \\|u\\|_{\\mathcal A}^2 + 2\\|u\\|_{\\mathcal A}\\|v\\|_{\\mathcal A} + \\|v\\|_{\\mathcal A}^2\\nonumber\\\\\n",
" %\n",
" &= (\\|u\\|_{\\mathcal A} + \\|v\\|_{\\mathcal A})^2,\n",
"\\end{align}\n",
"\n",
"and so we the triangle inequality $\\|u+v\\|_{\\mathcal A}\\leq \\|u\\|_{\\mathcal A} + \\|v\\|_{\\mathcal A}$. $\\square$\n",
"\n",
"It turns out that the solution $u_h$ to $(\\text{G}_h)$ is the orthogonal projection of $u$ (the solution to $(\\text{weak})$) to $V_h$:\n",
"\n",
"::: {#thm-EN}\n",
"# Energy norm\n",
"\n",
"Let $u$ solve $\\text{(weak)}$ and $u_h$ solve $(\\text{G}_h)$. Then, \n",
"\n",
"*(i)* Galerkin orthogonality: $\\mathcal A(u-u_h, v_h) = 0$ for all $v_h \\in V_h$. \n",
"\n",
"*(ii)* $\\|u-u_h\\|_{\\mathcal A} = \\min_{v_h \\in V_h} \\|u-v_h\\|_{\\mathcal A}$\n",
"\n",
":::\n",
"\n",
"**Proof.** \n",
"\n",
"*(i)* Using the linearity of $\\mathcal A$ and the fact that $V_h \\subset V$, we have\n",
"\n",
"\\begin{align}\n",
" \\mathcal A(u - u_h, v_h) &= \\mathcal A(u, v_h) - \\mathcal A(u_h, v_h)\\nonumber\\\\\n",
" %\n",
" &= (f,v_h) - (f, v_h) = 0\n",
"\\end{align}\n",
"\n",
"for all $v_h \\in V_h$. \n",
"\n",
"*(ii)* Using the bilinearity of $\\mathcal A$, we have\n",
"\n",
"\\begin{align}\n",
" \\|u-v_h\\|_{\\mathcal A}^2 &= \\mathcal A(u-v_h, u-v_h) \\nonumber\\\\\n",
" %\n",
" &= \\mathcal A(u-u_h + u_h - v_h, u-u_h+u_h-v_h) \\nonumber\\\\\n",
" %\n",
" &= \\mathcal A(u-u_h, u-u_h)\n",
" %\n",
" + 2 \\mathcal A(u-u_h, u_h-v_h) \\nonumber\\\\\n",
" % \n",
" &\\qquad + \\mathcal A(u_h - v_h, u_h-v_h) \\nonumber\\\\\n",
" %\n",
" &= \\|u-u_h\\|_{\\mathcal A}^2 + \\|u_h-v_h\\|_{\\mathcal A}^2.\n",
"\\end{align}\n",
"\n",
"Here, we have used the Galerkin orthogonality $\\mathcal A(u - u_h, u_h-v_h) = 0$. As a result, we have $\\|u-u_h\\|_{\\mathcal A} \\leq \\|u-v_h\\|_{\\mathcal A}$ for all $v_h \\in V_h$. $\\square$\n",
"\n",
"This theorem shows that $u_h$ is the best approximation (with respect to the energy norm) to $u$ in the space $V_h$. We next show how to use this to get an error estimate $\\| u - u_h \\|_{\\mathcal A}$ as a function of the mesh size $h$.\n",
"\n",
"Recall that the functions $\\phi_j$ are *cardinal* (that is, $\\phi_j(x_i) = \\delta_{ij}$). Therefore, the *interpolation* \n",
"\n",
"\\begin{align}\n",
" I_h u(x) := \\sum_{j=1}^m u(x_j) \\phi_j(x)\n",
"\\end{align}\n",
"\n",
"is a function that belongs to $V_h$ and satisfies $I_h u(x_j) = u(x_j)$. Therefore, by the theorem above, we have \n",
"\n",
"\\begin{align}\n",
" \\| u - u_h \\|_{\\mathcal A} \\leq \\| u - I_h u \\|_{\\mathcal A}. \n",
"\\end{align}\n",
"\n",
"We can obtain error estimates for finite element methods by finding an upper bound on the right hand side of this inequality. \n",
"\n",
"::: {#thm-}\n",
"\n",
"Suppose that $u$ is twice continuously differentiable on $[a,b]$. Then, piecewise linear finite elements converge with rate $O(h)$ in the energy norm:\n",
"\n",
"\\begin{align}\n",
" \\| u - u_h \\|_{\\mathcal A} \\leq C h \\| u''\\|_{L^2([a,b])}\n",
"\\end{align} \n",
"\n",
"where $\\| u''\\|_{L^2([a,b])} := \\sqrt{\\int_a^b |u''(t)|^2 \\mathrm{d}t}$ denotes the $L^2$ norm of $u''$.\n",
"\n",
":::\n",
"\n",
"\n",
"**Proof.** \n",
"\n",
"We apply @thm-EN and estimate the right hand side of $\\| u - u_h \\|_{\\mathcal A}\\leq \\| u - I_h u \\|_{\\mathcal A}$ where $I_hu$ is the interpolation defined above. (The proof that the right hand side of this is $O(h)$ is similar to a proof we had last semester).\n",
"\n",
"We will show that, on defining $v := u - I_hu$, we have \n",
"\n",
"\\begin{align}\n",
" \\|v\\|_{L^2([a,b])} &\\leq h \\|v'\\|_{L^2([a,b])} \n",
" %\n",
" \\leq h^2 \\|v''\\|_{L^2([a,b])}.\n",
" \\tag{*}\n",
"\\end{align}\n",
"\n",
"If this is true, then since $v'' = u''$, we have \n",
"\n",
"\\begin{align}\n",
" \\| u - u_h \\|_{\\mathcal A}^2 &\\leq \\| u - I_hu \\|_{\\mathcal A}^2 \\nonumber\\\\\n",
" %\n",
" &= \\int_a^b \\left( \\alpha(x) |v'(x)|^2 + \\beta(x) |v(x)|^2 \\right) \\mathrm{d}x \\nonumber\\\\\n",
" %\n",
" &\\leq \\Big( \\max_{[a,b]} |\\alpha| \\Big) \\|v'\\|^2_{L^2([a,b])} + \\Big( \\max_{[a,b]} |\\beta| \\Big)\\|v\\|^2_{L^2([a,b])} \\nonumber\\\\\n",
" %\n",
" &\\leq \\left[ h^2 \\Big( \\max_{[a,b]} |\\alpha| \\Big) + h^4 \\Big( \\max_{[a,b]} |\\beta| \\Big) \\right] \\|u''\\|^2_{L^2([a,b])}.\n",
"\\end{align}\n",
"\n",
"We conclude by taking square roots.\n",
"\n",
"Therefore, we only need to prove *Poincaré's inequality* $(*)$ (Poincaré's inequality is a bound of the form $\\|v\\|_{L^2} \\leq \\|v'\\|_{L^2}$ which is true e.g. when $v$ is zero on the boundary, or if $v$ has zero mean $\\int v = 0$). To prove this we need the Cauchy-Schwarz inequaltiy: $\\int fg \\leq \\|f\\|_{L^2} \\|g\\|_{L^2}$ (the proof of this is exactly the same as before:\n",
"\n",
"\\begin{align}\n",
" 0 &\\leq (f+\\lambda g,f+\\lambda g) \\nonumber\\\\\n",
" &= \\int |f|^2 + 2\\lambda \\int f g + \\lambda^2 \\int|g|^2. \n",
"\\end{align}\n",
"\n",
"This is a quadratic in $\\lambda$ and the discriminant being negative is equivalent to C-S.)\n",
"\n",
"Recall that $v := u - I_hu$. Since $I_hu$ is the piecewise linear interpolation of $u$ on $\\{x_j\\}$, we have $v(x_{k-1}) = v(x_k) = 0$. Moreover, since $v$ is twice continuously differentiable, there exists $\\xi_k\\in[x_{k-1},x_k]$ such that $v'(\\xi_k) = 0$. Therefore, we have \n",
"\n",
"\\begin{align}\n",
" v(x) &= v(x_{k-1}) + \\int_{x_{k-1}}^x v'(t) \\mathrm{d}t = \\int_{x_{k-1}}^x v'(t) \\mathrm{d}t \\nonumber\\\\\n",
" v'(x) &= v'(\\xi_{k}) + \\int^{x}_{\\xi_k} v''(t) \\mathrm{d}t\n",
" = \\int^{x}_{\\xi_k} v''(t) \\mathrm{d}t. \n",
"\\end{align}\n",
"\n",
"In the following, we will apply the the exact same argument to both $w = v$ and $w = v'$. On defining $y_k := x_{k-1}$ if $w=v$ and $y_k := \\xi_{k}$ if $w=v'$ and apply the Cauchy-Schwarz inequality, we have: for $x \\in [x_{k-1},x_k]$,\n",
"\n",
"\\begin{align}\n",
" |w(x)|^2\n",
" %\n",
" &= \\left|\\int_{y_{k}}^x 1\\cdot w'(t) \\mathrm{d}t\\right|^2 \\nonumber\\\\\n",
" %\n",
" &\\leq \\left(\\int_{y_k}^x 1^2\\mathrm{d}t\\right) \\left(\\int_{y_{k}}^x |w'(t)|^2 \\mathrm{d}t\\right) \\nonumber\\\\\n",
" %\n",
" &\\leq h \\int_{x_{k-1}}^{x_{k}} |w'(t)|^2 \\mathrm{d}t.\n",
"\\end{align}\n",
"\n",
"As a result, we have $\\int_{x_{k-1}}^{x_k} |w(x)|^2 \\mathrm{d}x \\leq h^2 \\int_{x_{k-1}}^{x_{k}} |w'(t)|^2 \\mathrm{d}t$ and, summing over $k$, we obtain $\\|w\\|_{L^2([a,b])}^2 \\leq h^2 \\| w'\\|_{L^2([a,b])}^2$. $\\square$\n",
"\n",
"## Potential presentation topic\n",
"\n",
"@Suli chapter 14.5 - *A posteriori* error analysis and mesh refinement. "
]
},
{
"cell_type": "markdown",
"id": "b1d5ab0c",
"metadata": {},
"source": [
":::{.callout-tip}\n",
"# TO-DO\n",
"\n",
"* Take a look at the suggestions for [Presentations/Posters/Papers](/Presentations.html) and sign up to present, \n",
"* The last few lectures have been more difficult: Read through the proofs/exercises from the notes and make sure you are up-to-date, \n",
"* (Sign-up to office hours if you would like)\n",
"\n",
"Previous reading:\n",
"\n",
"* Chapter 1 reading:\n",
" * @Burden sections 5.1 - 5.6, \n",
" * @FNC [zero stability](https://fncbook.com/zerostability/) \n",
"* Chapter 2 reading:\n",
" * Recap matrices: sections 7.1 - 7.2 of @Burden \n",
" * Jacobi & Gauss-Seidel: section 7.3 @Burden \n",
" * CG: section 7.6 @Burden\n",
"* Chapter 3 reading: \n",
" * Sections 11.1 - 11.2 of @Burden: shooting methods\n",
" * Lecture 16: @FNC [Galerkin method](https://fncbook.com/galerkin/)\n",
" * Today: Chapter 14.4 of @Suli: Theorem 14.6. \n",
"\n",
":::"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 1.11",
"language": "julia",
"name": "julia-1.11"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.11.6"
}
},
"nbformat": 4,
"nbformat_minor": 5
}