{
"cells": [
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Python version 3.8.10 (default, Nov 26 2021, 20:14:08) \n",
"[GCC 9.3.0]\n"
]
}
],
"source": [
"import sys #only needed to determine Python version number\n",
"print('Python version ' + sys.version)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# M62 TP 4 - Étude de la convergence"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"application/javascript": [
"IPython.notebook.set_autosave_interval(300000)"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Autosaving every 300 seconds\n"
]
}
],
"source": [
"%reset -f\n",
"%matplotlib inline\n",
"%autosave 300\n",
"\n",
"from matplotlib.pylab import *\n",
"from scipy.optimize import fsolve"
]
},
{
"cell_type": "markdown",
"metadata": {
"toc": true
},
"source": [
"Table of Contents
\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "WyG-bTRQE3f6"
},
"source": [
"# Étude de la convergence\n",
"\n",
"## Exercice\n",
"\n",
"Considérons le problème de Cauchy\n",
"\n",
">trouver la fonction $y \\colon I\\subset \\mathbb{R} \\to \\mathbb{R}$ définie sur l'intervalle $I=[0,1]$ telle que\n",
"$$\n",
"\\begin{cases}\n",
"y'(t) = y(t), &\\forall t \\in I=[0,1],\\\\\n",
"y(0) = 1\n",
"\\end{cases}\n",
"$$\n",
"dont la solution est $y(t)=e^{t}$. \n",
"\n",
"Estimer l'ordre de convergence des méthodes:\n",
"1. EE, AB$_2$, AB$_3$, AB$_4$, AB$_5$, N$_2$, N$_3$, N$_4$, EM, RK4$_1$, RK6$_5$, RK7$_6$ \n",
"2. EI, CN, AM$_2$, AM$_3$, AM$_4$, AM$_5$, MA$_2$, BDF$_2$, BDF$_3$, RK1_M\n",
"3. Heun, AM$_4$-AB$_1$, AM$_4$-AB$_2$, AM$_4$-AB$_3$\n",
" \n",
" \n",
"Remarque: puisque la fonction $\\varphi(t,y)=y$ est linéaire, toute méthode implicite peut être rendue explicite par un calcul élémentaire en explicitant directement pour chaque schéma l'expression de $u_{n+1}$. Cependant, nous pouvons utiliser le le module `SciPy` sans modifier l'implémentation des schémas (mais on payera l'ordre de convergence de `fsolve`).\n",
"\n",
"**Attention: les schémas multistep ont besoin d'initialiser plusieurs pas de la suite définie pas récurrence pour pouvoir démarrer. \n",
"Dans cette étude, au lieu d'utiliser un schéma d'ordre inférieur pour initialiser la suite, on utilisera la solution exacte (en effet, l'utilisation d'un schéma d'ordre inférieur dégrade l'ordre de précision).**"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" Rappels de la démarche pour le calcul de l'orde de convergence.
\n",
" \n",
"+ On écrit les schémas numériques :\n",
" + les nœuds d'intégration $[t_0,t_1,\\dots,t_{N}]$ sont contenus dans le vecteur `tt` (qui change en fonction de `h`) \n",
" + les valeurs $[u_0,u_1,\\dots,u_{N}]$ pour chaque méthode sont contenues dans le vecteur `uu`.\n",
"+ Pour chaque schéma, on calcule la solution approchée avec différentes valeurs de $h_k=1/N_k$. On sauvegarde les valeurs de $h_k$ dans le vecteur `H`. \n",
"+ Pour chaque valeur de $h_k$, on calcule le maximum de la valeur absolue de l'erreur et on sauvegarde toutes ces erreurs dans le vecteur `err_schema` de sort que `err_schema[k]` contient $e_k=\\max_{i=0,\\dots,N_k}|y(t_i)-u_{i}|$. \n",
"+ Pour afficher l'ordre de convergence on utilise une échelle logarithmique, i.e. on représente $\\ln(h)$ sur l'axe des abscisses et $\\ln(\\text{err})$ sur l'axe des ordonnées. \n",
" En effet, si $\\text{err}=Ch^p$ alors $\\ln(\\text{err})=\\ln(C)+p\\ln(h)$. \n",
" En échelle logarithmique, $p$ représente donc la pente de la ligne droite $\\ln(\\text{err})$.\n",
"\n",
" "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Schémas explicites"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "y69SGZjfIDo9"
},
"source": [
"### Schéma de Adam-Bashforth à 1 pas = schéma d'Euler progressif\n",
"$$\n",
"\\begin{cases}\n",
"u_0=y_0,\\\\\n",
"u_{n+1}=u_n+h\\varphi(t_n,u_n)& n=0,1,2,\\dots N-1\n",
"\\end{cases}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {
"colab": {
"autoexec": {
"startup": false,
"wait_interval": 0
}
},
"colab_type": "code",
"id": "_Bgo6mNyIQgu"
},
"outputs": [],
"source": [
"def EE(phi,tt,sol_exacte):\n",
" h = tt[1]-tt[0]\n",
" uu = [sol_exacte(tt[0])]\n",
" for i in range(N):\n",
" uu.append(uu[i]+h*phi(tt[i],uu[i]))\n",
" return uu"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "bJ2pbhejIQM2"
},
"source": [
"### Schéma de Adam-Bashforth à 2 pas\n",
"$$\n",
"\\begin{cases}\n",
"u_0=y_0,\\\\\n",
"u_{1}=y_1,\\\\\n",
"u_{n+1}=u_n+\\frac{h}{2}\\Bigl(3\\varphi(t_n,u_n)-\\varphi(t_{n-1},u_{n-1})\\Bigr)& n=1,2,3,4,5,\\dots N-1\n",
"\\end{cases}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {
"colab": {
"autoexec": {
"startup": false,
"wait_interval": 0
}
},
"colab_type": "code",
"id": "g38fKrIgSiBQ"
},
"outputs": [],
"source": [
"def AB2(phi,tt,sol_exacte):\n",
" h = tt[1]-tt[0]\n",
" uu = [sol_exacte(tt[0])]\n",
" uu.append(sol_exacte(tt[1]))\n",
" for i in range(1,N):\n",
" k1 = phi( tt[i], uu[i] )\n",
" k2 = phi( tt[i-1], uu[i-1] )\n",
" uu.append( uu[i] + (3*k1-k2)*h/2 )\n",
" return uu"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "nI8swyc6RxIR"
},
"source": [
"### Schéma de Adam-Bashforth à 3 pas\n",
"$$\n",
"\\begin{cases}\n",
"u_0=y_0,\\\\\n",
"u_{1}=y_1,\\\\\n",
"u_{2}=y_2,\\\\\n",
"u_{n+1}=u_n+\\frac{h}{12}\\Bigl(23\\varphi(t_n,u_n)-16\\varphi(t_{n-1},u_{n-1})+5\\varphi(t_{n-2},u_{n-2})\\Bigr)& n=2,3,4,5,\\dots N-1\n",
"\\end{cases}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "awcWzBp7SXvQ"
},
"source": [
"### Schéma de Adam-Bashforth à 4 pas\n",
"$$\n",
"\\begin{cases}\n",
"u_0=y_0,\\\\\n",
"u_{1}=y_1,\\\\\n",
"u_{2}=y_2,\\\\\n",
"u_{3}=y_3,\\\\\n",
"u_{n+1}=u_n+\\frac{h}{24}\\Bigl(55\\varphi(t_n,u_n)-59\\varphi(t_{n-1},u_{n-1})+37\\varphi(t_{n-2},u_{n-2})-9\\varphi(t_{n-3},u_{n-3})\\Bigr)& n=3,4,5,\\dots N-1\n",
"\\end{cases}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "heLmvMe_S0y6"
},
"source": [
"### Schéma de Adam-Bashforth à 5 pas\n",
"$$\n",
"\\begin{cases}\n",
"u_0=y_0,\\\\\n",
"u_{1}=y_1,\\\\\n",
"u_{2}=y_2,\\\\\n",
"u_{3}=y_3,\\\\\n",
"u_{4}=y_4,\\\\\n",
"u_{n+1}=u_n+\\frac{h}{720}\\Bigl(1901\\varphi(t_n,u_n)-2774\\varphi(t_{n-1},u_{n-1})+2616\\varphi(t_{n-2},u_{n-2})-1274\\varphi(t_{n-3},u_{n-3})+251\\varphi(t_{n-4},u_{n-4})\\Bigr)& n=4,5,\\dots N-1\n",
"\\end{cases}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "OldEmxFfTJfq"
},
"source": [
"### Schéma de Nylström à 2 pas\n",
"$$\n",
"\\begin{cases}\n",
"u_0=y_0,\\\\\n",
"u_{1}=y_1,\\\\\n",
"u_{n+1}=u_{n-1}+2h\\varphi(t_{n},u_{n})& n=1,2,3,4,5,\\dots N-1\n",
"\\end{cases}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "dADQEhyYTVQz"
},
"source": [
"### Schéma de Nylström à 3 pas\n",
"$$\n",
"\\begin{cases}\n",
"u_0=y_0,\\\\\n",
"u_{1}=y_1,\\\\\n",
"u_{2}=y_2,\\\\\n",
"u_{n+1}=u_{n-1}+\\frac{h}{3}\\Bigl(7\\varphi(t_{n},u_{n})-2\\varphi(t_{n-1},u_{n-1})+\\varphi(t_{n-2},u_{n-2})\\Bigr)& n=2,3,4,5,\\dots N-1\n",
"\\end{cases}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "LsgdqQnfTf66"
},
"source": [
"### Schéma de Nylström à 4 pas\n",
"$$\n",
"\\begin{cases}\n",
"u_0=y_0,\\\\\n",
"u_{1}=y_1,\\\\\n",
"u_{2}=y_2,\\\\\n",
"u_{3}=y_3,\\\\\n",
"u_{n+1}=u_{n-1}+\\frac{h}{3}\\Bigl(8\\varphi(t_{n},u_{n})-5\\varphi(t_{n-1},u_{n-1})+4\\varphi(t_{n-2},u_{n-2})-\\varphi(t_{n-3},u_{n-3})\\Bigr)& n=3,4,5,\\dots N-1\n",
"\\end{cases}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "p4f0txAsIwNG"
},
"source": [
"### Schéma d'Euler modifié\n",
"$$\n",
"\\begin{cases}\n",
"u_0=y_0,\\\\\n",
"\\tilde u = u_n+\\frac{h}{2}\\varphi(t_n,u_n),\\\\\n",
"u_{n+1}=u_n+h\\varphi\\left(t_n+\\frac{h}{2},\\tilde u\\right)& n=0,1,2,\\dots N-1\n",
"\\end{cases}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {
"colab": {
"autoexec": {
"startup": false,
"wait_interval": 0
}
},
"colab_type": "code",
"id": "O5rOYvtPI7TO"
},
"outputs": [],
"source": [
"def EM(phi,tt,sol_exacte):\n",
" h = tt[1]-tt[0]\n",
" uu = [sol_exacte(tt[0])]\n",
" for i in range(N):\n",
" k1 = phi( tt[i], uu[i] )\n",
" uu.append( uu[i]+h*phi(tt[i]+h/2,uu[i]+k1*h/2) )\n",
" return uu"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "_50Xo95mT9tC"
},
"source": [
"### Schéma de Runge-Kutta RK4-1\n",
"$$\\begin{cases}\n",
"u_0 = y_0 \\\\\n",
"K_1 = \\varphi\\left(t_n,u_n\\right)\\\\\n",
"K_2 = \\varphi\\left(t_n+\\frac{h}{2},u_n+\\frac{h}{2} K_1)\\right)\\\\\n",
"K_3 = \\varphi\\left(t_n+\\frac{h}{2},u_n+\\frac{h}{2}K_2\\right)\\\\\n",
"K_4 = \\varphi\\left(t_{n+1},u_n+h K_3\\right)\\\\\n",
"u_{n+1} = u_n + \\frac{h}{6}\\left(K_1+2K_2+2K_3+K_4\\right) & n=0,1,\\dots N-1\n",
"\\end{cases}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Schéma de Runge-Kutta RK6_5 (Butcher à 6 étages d'ordre 5)\n",
"$$\\begin{array}{c|cccccc} 0 & 0 & 0 & 0 & 0 &0&0 \\\\ \\frac{1}{4} & \\frac{1}{4} & 0&0&0&0&0\\\\ \\frac{1}{4} & \\frac{1}{8} &\\frac{1}{8}&0&0&0&0\\\\ \\frac{1}{2} & 0 &-\\frac{1}{2}&1&0&0&0\\\\ \\frac{3}{4} & \\frac{3}{16} &0&0&\\frac{9}{16}&0&0\\\\ 1 & \\frac{-3}{7} &\\frac{2}{7}&\\frac{12}{7}&\\frac{-12}{7}&\\frac{8}{7}&0\\\\ \\hline & \\frac{7}{90} & 0&\\frac{32}{90} & \\frac{12}{90} & \\frac{32}{90} & \\frac{7}{90} \\end{array}$$\n",
"\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Schéma de Runge-Kutta RK7_6 (Butcher à 7 étages d'ordre 6)\n",
"$$\n",
"\\begin{array}{c|ccccccc} \n",
"0 & 0 & 0 & 0 & 0 &0&0&0 \\\\ \n",
"\\frac{1}{2} & \\frac{1}{2} & 0&0&0&0&0&0\\\\ \n",
"\\frac{2}{3} & \\frac{2}{9} &\\frac{4}{9}&0&0&0&0&0\\\\ \n",
"\\frac{1}{3} & \\frac{7}{36} &\\frac{2}{9}&-\\frac{1}{12}&0&0&0&0\\\\ \n",
"\\frac{5}{6} & -\\frac{35}{144} &-\\frac{55}{36}&\\frac{35}{48}&\\frac{15}{8}&0&0&0\\\\ \n",
"\\frac{1}{6} & -\\frac{1}{360} &-\\frac{11}{36}&-\\frac{1}{8}&\\frac{1}{2}&\\frac{1}{10}&0&0\\\\ \n",
"1 & \\frac{-41}{260} &\\frac{22}{13}&\\frac{43}{156}&-\\frac{118}{39}&\\frac{32}{195}&\\frac{80}{39}&0\\\\ \n",
"\\hline \n",
" & \\frac{13}{200} & 0&\\frac{11}{40} & \\frac{11}{40} & \\frac{4}{25} & \\frac{4}{25} & \\frac{13}{200} \\end{array}$$\n",
"\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Schémas implicites"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "y69SGZjfIDo9"
},
"source": [
"### Schéma d'Euler régressif\n",
"$$\n",
"\\begin{cases}\n",
"u_0=y_0,\\\\\n",
"u_{n+1}=u_n+h\\varphi(t_{n+1},u_{n+1})& n=0,1,2,\\dots N-1\n",
"\\end{cases}\n",
"$$\n",
"avec $u_{n+1}$ zéro de la fonction $$x\\mapsto -x+u_n+h\\varphi(t_{n+1},x).$$"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {
"colab": {
"autoexec": {
"startup": false,
"wait_interval": 0
}
},
"colab_type": "code",
"id": "_Bgo6mNyIQgu"
},
"outputs": [],
"source": [
"def EI(phi,tt,sol_exacte):\n",
" h = tt[1]-tt[0]\n",
" uu = [sol_exacte(tt[0])]\n",
" for i in range(N):\n",
" temp = fsolve(lambda x: -x+uu[i]+h*phi(tt[i+1],x), uu[i])\n",
" uu.append(temp[0])\n",
" return uu"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "23PyYGzhQwuo"
},
"source": [
"### Schéma de Crank-Nicolson\n",
"$$\n",
"\\begin{cases}\n",
"u_0=y_0,\\\\\n",
"u_{n+1}=u_n+\\frac{h}{2}\\Bigl(\\varphi(t_n,u_n)+\\varphi(t_{n+1},u_{n+1})\\Bigr)& n=0,1,2,\\dots N-1\n",
"\\end{cases}\n",
"$$\n",
"avec $u_{n+1}$ un zéro de la fonction $x\\mapsto -x+u_n+\\frac{h}{2}(\\varphi(t_n,u_n)+\\varphi(t_{n+1},x))$."
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [],
"source": [
"def CN(phi,tt,sol_exacte):\n",
" h = tt[1]-tt[0]\n",
" uu = [sol_exacte(tt[0])]\n",
" for i in range(len(tt)-1):\n",
" temp = fsolve(lambda x: -x+uu[i]+0.5*h*( phi(tt[i+1],x)+phi(tt[i],uu[i]) ), uu[i])\n",
" uu.append(temp[0])\n",
" return uu"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "23PyYGzhQwuo"
},
"source": [
"### Schéma de AM-2\n",
"$$\n",
"\\begin{cases}\n",
"u_0=y_0,\\\\\n",
"u_1=y_1,\\\\\n",
"u_{n+1}=u_n+\\frac{h}{12}\\Bigl(5\\varphi(t_{n+1},u_{n+1})+8\\varphi(t_n,u_n)-\\varphi(t_{n-1},u_{n-1})\\Bigr)& n=1,2,\\dots N-1\n",
"\\end{cases}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [],
"source": [
"def AM2(phi,tt,sol_exacte):\n",
" h = tt[1]-tt[0]\n",
" uu = [sol_exacte(tt[0])]\n",
" uu.append(sol_exacte(tt[1]))\n",
" for i in range(1,N):\n",
" temp = fsolve(lambda x: -x+uu[i]+h*( 5*phi(tt[i+1],x)+8*phi(tt[i],uu[i])-phi(tt[i-1],uu[i-1]) )/12, uu[i])\n",
" uu.append(temp[0])\n",
" return uu"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "23PyYGzhQwuo"
},
"source": [
"### Schéma de AM-3\n",
"$$\n",
"\\begin{cases}\n",
"u_0=y_0,\\\\\n",
"u_1=y_1,\\\\\n",
"u_{2}=y_2,\\\\\n",
"u_{n+1}=u_n+\\frac{h}{24}\\Bigl(9\\varphi(t_{n+1},u_{n+1})+19\\varphi(t_n,u_n)-5\\varphi(t_{n-1},u_{n-1})+\\varphi(t_{n-2},u_{n-2})\\Bigr)& n=2,3,\\dots N-1\n",
"\\end{cases}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "23PyYGzhQwuo"
},
"source": [
"### Schéma de AM-4\n",
"$$\n",
"\\begin{cases}\n",
"u_0=y_0,\\\\\n",
"u_1=y_1,\\\\\n",
"u_{2}=y_2,\\\\\n",
"u_{3}=y_3,\\\\\n",
"u_{n+1}=u_n+\\frac{h}{720}\\Bigl(251\\varphi(t_{n+1},u_{n+1})+646\\varphi(t_n,u_n)-264\\varphi(t_{n-1},u_{n-1})+106\\varphi(t_{n-2},u_{n-2})-19\\varphi(t_{n-3},u_{n-3})\\Bigr)& n=3,4,\\dots N-1\n",
"\\end{cases}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "23PyYGzhQwuo"
},
"source": [
"### Schéma de AM-5\n",
"$$\n",
"\\begin{cases}\n",
"u_0=y_0,\\\\\n",
"u_1=y_1,\\\\\n",
"u_{2}=y_2,\\\\\n",
"u_{3}=y_3,\\\\\n",
"u_{4}=y_4,\\\\\n",
"u_{n+1}=u_n+\\frac{h}{1440}\\Bigl(475\\varphi(t_{n+1},u_{+1})+1427\\varphi(t_n,u_n)-798\\varphi(t_{n-1},u_{n-1})+482\\varphi(t_{n-2},u_{n-2})-173\\varphi(t_{n-3},u_{n-3})+27\\varphi(t_{n-4},u_{n-4})\\Bigr),& n=4,5,\\dots N-1\n",
"\\end{cases}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Schéma MS-2\n",
"$$\n",
"\\begin{cases}\n",
"u_0=y_0,\\\\\n",
"u_1=y_1,\\\\\n",
"u_{n+1}=u_{n-1}+\\frac{h}{3}\\Bigl(\\varphi(t_{n+1},u_{n+1})+4\\varphi(t_n,u_n)+\\varphi(t_{n-1},u_{n-1})\\Bigr)& n=1,2,\\dots N-1\n",
"\\end{cases}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "23PyYGzhQwuo"
},
"source": [
"### Schéma BDF2\n",
"$$\n",
"\\begin{cases}\n",
"u_0=y_0,\\\\\n",
"u_1=y_1,\\\\\n",
"u_{n+1}=\\frac{4}{3}u_n-\\frac{1}{3}u_{n-1}+\\frac{2}{3}h\\varphi(t_{n+1},u_{n+1})& n=1,2,\\dots N-1\n",
"\\end{cases}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "23PyYGzhQwuo"
},
"source": [
"### Schéma BDF3\n",
"$$\n",
"\\begin{cases}\n",
"u_0=y_0,\\\\\n",
"u_1=y_1,\\\\\n",
"u_{2}=y_2,\\\\\n",
"u_{n+1}=\\frac{18}{11}u_n-\\frac{9}{11}u_{n-1}+\\frac{2}{11}u_{n-2}+\\frac{6}{11}h\\varphi(t_{n+1},u_{n+1})& n=2,3,\\dots N-1\n",
"\\end{cases}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Schéma RK1_M:\n",
"$$\n",
"\\begin{cases}\n",
"u_0=y_0,\\\\\n",
"u_{n+1}=u_n+h\\varphi\\left(\\frac{t_n+t_{n+1}}{2},\\frac{u_n+u_{n+1}}{2}\\right)& n=0,1,2,\\dots N-1\n",
"\\end{cases}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Schémas predicteur-correcteur"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "23PyYGzhQwuo"
},
"source": [
"### Schéma de Heun\n",
"$$\n",
"\\begin{cases}\n",
"u_0=y_0,\\\\\n",
"\\tilde u = u_n+h\\varphi(t_n,u_n)\\\\\n",
"u_{n+1}=u_n+\\frac{h}{2}\\Bigl(\\varphi(t_n,u_n)+\\varphi(t_{n+1},\\tilde u)\\Bigr)& n=0,1,2,\\dots N-1\n",
"\\end{cases}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {
"colab": {
"autoexec": {
"startup": false,
"wait_interval": 0
}
},
"colab_type": "code",
"id": "1ewZyxhHRYxg"
},
"outputs": [],
"source": [
"def heun(phi,tt,sol_exacte):\n",
" h = tt[1]-tt[0]\n",
" uu = [sol_exacte(tt[0])]\n",
" for i in range(N):\n",
" k1 = phi( tt[i], uu[i] )\n",
" k2 = phi( tt[i+1], uu[i] + k1*h )\n",
" uu.append( uu[i] + (k1+k2)*h/2 )\n",
" return uu"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "23PyYGzhQwuo"
},
"source": [
"### Schéma AM-2 AB-1\n",
"$$\n",
"\\begin{cases}\n",
"u_0=y_0,\\\\\n",
"u_1=y_1,\\\\\n",
"\\tilde u=u_n+h\\varphi(t_n,u_n),\\\\\n",
"u_{n+1}=u_n+\\frac{h}{12}\\Bigl(5\\varphi(t_{n+1},\\tilde u)+8\\varphi(t_n,u_n)-\\varphi(t_{n-1},u_{n-1})\\Bigr)& n=1,2,\\dots N-1\n",
"\\end{cases}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "23PyYGzhQwuo"
},
"source": [
"### Schéma AM-3 AB-2\n",
"$$\n",
"\\begin{cases}\n",
"u_0=y_0,\\\\\n",
"u_1=y_1,\\\\\n",
"u_2=y_2,\\\\\n",
"\\tilde u=u_n+\\frac{h}{2}(3\\varphi(t_n,u_n)-\\varphi(t_{n-1},u_{n-1})),\\\\\n",
"u_{n+1}=u_n+\\frac{h}{24}\\Bigl(9\\varphi(t_{n+1},\\tilde u)+19\\varphi(t_n,u_n)-5\\varphi(t_{n-1},u_{n-1})+\\varphi(t_{n-2},u_{n-2})\\Bigr)& n=2,3,\\dots N-1\n",
"\\end{cases}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Schéma AM-4 AB-2/3/4/5\n",
"Écrire et implémenter les schémas AM-4 AB-2, AM-4 AB-3, AM-4 AB-4, AM-4 AB-5."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Convergence"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "cnwNf75iGe0F"
},
"source": [
"On initialise le problème de Cauchy"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {
"colab": {
"autoexec": {
"startup": false,
"wait_interval": 0
}
},
"colab_type": "code",
"id": "OLLu4aFJFENg"
},
"outputs": [],
"source": [
"t0, y0, tfinal = 0, 1, 3"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "ve4iOfOIGsYc"
},
"source": [
"On définit la solution exacte:"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {
"colab": {
"autoexec": {
"startup": false,
"wait_interval": 0
}
},
"colab_type": "code",
"id": "W3EcAN2eGz2j"
},
"outputs": [],
"source": [
"sol_exacte = lambda t : exp(t)"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "xpjn_ogYGo20"
},
"source": [
"On définit l'équation différentielle : `phi` est une fonction python qui contient la fonction mathématique $\\varphi(t, y)$ dépendant des variables $t$ et $y$."
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {
"colab": {
"autoexec": {
"startup": false,
"wait_interval": 0
}
},
"colab_type": "code",
"id": "df9F-MXWGm2a"
},
"outputs": [],
"source": [
"phi = lambda t,y : y"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Pour chaque schéma, on calcule la solution approchée avec différentes valeurs de $h_k=1/N_k$, à savoir $N_k=2$, $2^2$, $2^3$, ... $2^{10}$). On sauvegarde les valeurs de $h_k$ dans le vecteur `H`. \n",
"\n",
"Pour chaque valeur de $h_k$, on calcule le maximum de la valeur absolue de l'erreur et on sauvegarde toutes ces erreurs dans le vecteur `err_schema` de sort que `err_schema[k]` contient $e_k=\\max_{i=0,...,N_k}|y(t_i)-u_{i}|$ avec $N_k=2^{k+1}$."
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [],
"source": [
"H = []\n",
"\n",
"err_ep = []\n",
"err_AB2 = []\n",
"\n",
"err_er = []\n",
"err_CN = []\n",
"\n",
"err_em = []\n",
"err_heun = []\n",
"\n",
"for k in range(8):\n",
" N = 2**(k+3)\n",
" tt = linspace(t0,tfinal,N+1)\n",
" h = tt[1]-tt[0]\n",
" H.append(h)\n",
" yy = array([sol_exacte(t) for t in tt])\n",
" \n",
" # schemas explicites\n",
" uu_ep = EE(phi,tt,sol_exacte)\n",
" uu_AB2 = AB2(phi,tt,sol_exacte) # multipas explicite lineaire\n",
" uu_em = EM(phi,tt,sol_exacte) # Runge-Kutta explicite\n",
" # schemas implicites\n",
" uu_er = EI(phi,tt,sol_exacte) \n",
" uu_CN = CN(phi,tt,sol_exacte) # multipas implicite lineaire\n",
" # schemas predictor-corrector\n",
" uu_heun = heun(phi,tt,sol_exacte) # CN_EE\n",
" \n",
" # erreurs\n",
" err_ep.append(norm(uu_ep-yy,inf))\n",
" err_AB2.append(norm(uu_AB2-yy,inf))\n",
" err_em.append(norm(uu_em-yy,inf))\n",
" \n",
" err_er.append(norm(uu_er-yy,inf))\n",
" err_CN.append(norm(uu_CN-yy,inf))\n",
" \n",
" err_heun.append(norm(uu_heun-yy,inf))"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "SnKKU27oGyQb"
},
"source": [
"Pour estimer l'ordre de convergence on estime la pente de la droite qui relie l'erreur au pas $k$ à l'erreur au pas $k+1$ en echelle logarithmique en utilisant la fonction `polyfit` basée sur la régression linéaire. "
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {
"colab": {
"autoexec": {
"startup": false,
"wait_interval": 0
}
},
"colab_type": "code",
"id": "ySox-VsNGt8p"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"EE 0.92\n",
"AB2 1.92\n",
"EM 1.95\n",
"\n",
"\n",
"EI 1.12\n",
"CN 2.01\n",
"\n",
"\n",
"Heun 1.95\n"
]
}
],
"source": [
"print (f'EE {polyfit(log(H),log(err_ep), 1)[0]:1.2f}' )\n",
"print (f'AB2 {polyfit(log(H),log(err_AB2), 1)[0]:1.2f}' )\n",
"print (f'EM {polyfit(log(H),log(err_em), 1)[0]:1.2f}' )\n",
"print('\\n')\n",
"print (f'EI {polyfit(log(H),log(err_er), 1)[0]:1.2f}' )\n",
"print (f'CN {polyfit(log(H),log(err_CN), 1)[0]:1.2f}' )\n",
"print('\\n')\n",
"print (f'Heun {polyfit(log(H),log(err_heun), 1)[0]:1.2f}' )"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "4HPS2hE6G54k"
},
"source": [
"Pour afficher l'ordre de convergence on utilise une échelle logarithmique : on représente $\\ln(h)$ sur l'axe des abscisses et $\\ln(\\text{err})$ sur l'axe des ordonnées. Le but de cette représentation est clair: si $\\text{err}=Ch^p$ alors $\\ln(\\text{err})=\\ln(C)+p\\ln(h)$. En échelle logarithmique, $p$ représente donc la pente de la ligne droite $\\ln(\\text{err})$."
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {
"colab": {
"autoexec": {
"startup": false,
"wait_interval": 0
},
"base_uri": "https://localhost:8080/",
"height": 3818,
"output_extras": [
{
"item_id": 1
},
{
"item_id": 2
},
{
"item_id": 3
},
{
"item_id": 4
},
{
"item_id": 5
},
{
"item_id": 6
},
{
"item_id": 7
},
{
"item_id": 8
},
{
"item_id": 9
},
{
"item_id": 10
},
{
"item_id": 11
}
]
},
"colab_type": "code",
"executionInfo": {
"elapsed": 2188,
"status": "ok",
"timestamp": 1520423878951,
"user": {
"displayName": "Gloria Faccanoni",
"photoUrl": "//lh4.googleusercontent.com/-gY6sCpFtBJo/AAAAAAAAAAI/AAAAAAAABdo/a_W4-RMG5X0/s50-c-k-no/photo.jpg",
"userId": "116371262733782746288"
},
"user_tz": -60
},
"id": "oz1tVYNtG4-3",
"outputId": "9b89f5ec-83c5-4797-e057-8d6c20560051"
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAABJgAAAFSCAYAAAC+H2XFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdeVxVdfrA8c8DgoAi4ZIbCi65p1hYVpY6pVlWlpZazmhWWrY4lraY46RNavuvxZqyLLNwy9JRczQrsbTRNMMScQ8VV8QNBGT7/v74XhDxsi8XLs/79bovuGd97uFyvuc857uIMQallFJKKaWUUkoppYrLw9UBKKWUUkoppZRSSqnKTRNMSimllFJKKaWUUqpENMGklFJKKaWUUkoppUpEE0xKKaWUUkoppZRSqkQ0waSUUkoppZRSSimlSkQTTEoppZRSSimllFKqRDTBpColEYkRkZtcHYeriUgPEYnN8T5KRHoUYr1EEWlepsEppZQLVbZyQkSGiMi3ZbTtWSLykuP360VkhyvjUUpVTZXtvFzWRCRERIyIVHO8/6+IDHN1XEqVhCaYlEuJSDcR+VlETovICRFZJyJdXB1XZWWMaW+MiSjEcjWNMXvhwhsPpZSqaKpKOWGMCTfG9C6H/fxkjGld1HgcN0EtyzY6pVRlUFXOy+XNGHOLMeazgpbTRJ2VO0GnKgb9YyiXEZFawDJgFLAA8AauB865Mi6llFIVg5YTSilVseh52TkREUCMMZmujqUkRKSaMSa9vLZdlvsrqYocW0WmNZiUK7UCMMbMNcZkGGOSjTHfGmN+z1pAREaISLSIJIjINhG5Isf6oSLyu+PpyXwR8cmx3m0iEikipxxPWDrmmBcjIk871j0rIjNFpL6jWmqCiHwnIoE5lv9SRI449vOjiLTPMe9WR1wJInJQRMbl9WFF5AHHZzkpIitFJNgx/VoROS4iTRzvOznibpMj3vGO/ZwUkU9zftZc+8h+oiEiniLyvIjsccT3a459GBFpKSIjgSHAM45mc0sd8xuJyFciEicif4rI6Bz7uEpENonIGRE5KiJv5vtXVkqp4qsy5YSI3C8ia3O8NyLyqIjscqz7LxFpISL/c5x/F4iIt2PZHiIS6zjnH3fEPySP/eRuWt1ERL52nO/jRWR67nhE5EfH4lscZcWgQhzDZx2fN0FEdojIjc7/xEqpSqaqnZfXici7ju1sz3kuE5EIEZkiIuuAJKC5iAQ4Yjvs2PZLIuLpWN5TRF53nKf3An1z7S9CRB7K7ziKyOdAU2Cp43z8jGPZro5jdkpEtkiOLjMkV40nEZkkIl84fs+qBfSgiOwHfsjjWNQWew9ySOz9yOJcce4WW5ttiYg0yjHPiMhjIrIL2JXPtDYissqxjR0iMjDHNnxF5A0R2ef4O6wVEV8gq2w65TgW1ziWz+ue66IaTzmPeY6/9/+JyAlgkrNjoQpgjNGXvlzyAmoB8cBnwC1AYK759wAHgS6AAC2BYMe8GOAXoBFQG4gGHnHMuwI4BlwNeALDHMtXz7HueqA+0Nix7GagM1Ade2J9IUccDwD+jnlvAZE55h0Grnf8HghckcdnvRPYDbTF1hz8B/BzjvlTHPv1BX4HHs8xLwbYCjRxfNZ1wEuOeT2A2FzL3uT4/WngD6C14/h1Auo45hmgpeP3WVnbc7z3AH4F/ol9KtUc2Avc7Jj/P+Bvjt9rAl1d/V3Sl7705Z6vKlZO3A+szfHeAEscx6A9tnbA945zcgCwDRjmWLYHkA686YihO3AWaO2Yn32ez1luOD77FuD/gBqAD9Atn3ha5nif5zHEljsHgEaOZUOAFq7+PulLX/oq+asKnpfTgScBL2AQcBqo7ZgfAex3nKOrOZZZDHzoOKde6vi8DzuWfwTYzvlr+tWOc2u1HNt7qJDH8aYccTZ2/E1uxV7H93K8r5fH8pOALxy/hzhimO2I2TePY/ENMN9xvLyA7o7pfwGOO/5+1YF3gR9zrGeAVY7P6+tsmmO/B4DhjuN4hWOb7R3Lv+c4No2x341rHfvKir1ajv3lec+Vx/I5j3nW3/sJx7pOj4W+8n9pDSblMsaYM0A37D/6R0CcI+td37HIQ8CrxpiNxtptjNmXYxPvGGMOGWNOAEuBUMf0EcCHxpgNxj5Z+Qx7Yd41x7rvGmOOGmMOAj8BG4wxvxljzgGLsIVVVpyfGGMSHPMmAZ1EJMAxOw1oJyK1jDEnjTGb8/i4DwPTjDHRxla1nIp9ghPsmD8Je8PwC3AIeyLNabox5oDjs04B7s3ruObwEPAPY8wOx/HbYoyJL8R6XbAF0ovGmFRj+2r6CBic4zO3FJG6xphEY8z6QmxTKaWKrIqVE868Yow5Y4yJwj5o+NYYs9cYcxr4b84YHCYaY84ZY9ZgbwYGkr+rsDd6TxtjzhpjUowxawtYJ0t+xzADe/HfTkS8jDExxpg9hdyuUqoCq4Ln5WPAW8aYNGPMfGAHF9Y8mmWMiXJc39fGJt3GOM6px7AJ/Kxr6IGObWVd00/LZ78FHcec/gosN8YsN8ZkGmNWAZuwCafCmuSIOTn3DBFp6PhcjziOV5qjnAHbEuITY8xmx7EeD1wjIiE5NjHNGHMi17ZzTrsNiDHGfGqMSXf8Pb4C7hYRD2yy8O/GmIOO78bPjn05U9A9V0EOGWPedcRx0bFQBdMEk3Ipxz///caYIKAD9kL3LcfsJkB+F6RHcvyehK1NAxAMjHVUET0lIqcc22qUY/mjOX5PdvK+JmRXZX1ZbDOzM9gnAAB1HT8HYE/e+0RkTVbVTCeCgbdzxHMC+zSiseM4pGGfMHcA3jDGptFzOJDj9325PkteCjp+eQkGGuU6fs9jnxgBPIitHr1dRDaKyG3F2IdSShVKFSonnClUDA4njTFnc7wvTFnRBNhnitfHRJ7H0BizGxiDvak7JiLzcjaZUEpVblXsvHww13V57nNrzmv0YGztnsM5PsOH2JpMONbLfU2fl6JcxwcD9+Q6dt2AhoVcn5xxicgHjiZniSLyvCOWE8aYk07Wa0SOz2GMScTWnmrsbNt5TAsGrs4V/xCgAfZv5kPRjkWe91yF4CxWVQSaYFIVhjFmO+eTLGD/wVsUY1MHgCnGmEtyvPyMMXOLsa37gH7ATdgaRiGO6eKIeaMxph+24FiM7ewwr5gezhWTrzHmZwARaQy8AHwKvCEi1XOt3yTH702xtZwKUtjj5yyZ9WeuWP2NMbcCGGN2GWPuxX7mV4CFIlKjEPtRSqkScfNyoqQCc52LC1NWHACaSvFG4Mn3GBpj5hhjumEv9g22vFBKuZkqcF5uLCKS433uc2vO6+gD2FpXdXN8hlrGmKz+nw5z8TV9XvI7js6u3T/PdexqGGNedsw/C/jlWL5Bfts0xjxi7IjTNY0xUx3bry0ilzhZ7xD2PA+Aoxyqg23el1e8uacdANbkir+mMWYUtqlcCs6PhbPt5nfPlfUQJr9j4Wybqgg0waRcxtGZ21gRCXK8b4Jt+pXV5OpjYJyIXClWy0JWb/wIeERErnasV0NE+oqIfzHC9McWFPHYk9HUHPF7i8gQEQlw1EA6g20W4MwHwHhxdDAotgPAexy/C7ZgnomtHXQY+Feu9R8TkSARqY2tTTS/ELF/DPxLRC5zHIeOIlLHyXJHsX16ZPkFOCO2g1Zfx1OgDuIYflZE/ioi9YwdJeOUY528PrdSShVbFSsnSsNkxz6vxzY5+LKA5X/BljkvO46Bj4hcl8eyucuKPI+hiLQWkb84HpakYGsWaDmhlBuoguflS4HRIuLluHZvCyx3tqAx5jDwLfZhcS0R8RA7OEN3xyILHNsKEtsh+XP57De/45j7fPwFcLuI3Oy4bvcRO6BDkGN+JDDY8RnCgLvz2W9en+u/wPsiEujYzg2O2XOA4SIS6jjnT8U2XYwpwi6WAa1E5G+ObXuJSBcRaeu43/gEeFPsIESeInKNY19xQGauY5HnPZcxJg6b+PqrYzsPULxkqMqHJpiUKyVgO/LbICJnsQXTVmAsgDHmS2x/Q3Mcyy7Gtm3OlzFmE7Yd93TgJLajt/uLGeNsbLXPg9gOVXP3N/Q3IEZs9dtHsG2gncW0CPv0dp5j2a3YtswAo7HNzyY6quAOx56or8+xiTnYAmuv4/VSIWJ/E1uQfYstPGdiO9LLbSa2HfopEVlsjMkAbse2if8T++TgY+wTIIA+QJSIJAJvA4ONMSmFiEcppYqqypQTpeAI9rMcAsKxfWVsz2+FHOf7ltiOamOxndg6Mwn4zFFWDCzgGFYHXsaWH0ewN2jPF/eDKaUqlKp2Xt4AXIY9n00B7jb592k6FDtIzjbs51jI+aZqHwErsYMrbAa+zmsjBRzHacA/HOfjccaYA9gaW89jky4HsIP9ZN3rT8QmUk4Ckx3bLKq/Yfuu2o7tl2qMI87vHdv/CvvAogXn+5wqFGNMAtDbsd4hbLnxCrYsARiHHbhoI7bJ2yuAhzEmCXuM1jmORdcC7rnAfseexiYf2wM/FyVWVTC5uKsXpVRFIiIx2NENvnN1LEoppSoescNRf+HoD0UppVQpEJH7sdfg3Vwdi1KVhdZgUkoppZRSSimllFIlogkmpZRSSimllFJKKVUi2kROKaWUUkoppZRSSpWI1mBSSimllFJKKaWUUiVSzdUBlIW6deuakJAQV4fh9s6ePUuNGjVcHYaqZPR741q//vrrcWNMPVfH4WpaTpQP/X9XxaHfG9fScsLScqJ86P+7Kg793rhWfuWEWyaYQkJC2LRpk6vDcHsRERH06NHD1WGoSka/N64lIvtcHUNFoOVE+dD/d1Uc+r1xLS0nLC0nyof+v6vi0O+Na+VXTmgTOaWUUlWGiNwuIjNOnz7t6lCUUkoppZRyKxU+wSQizUVkpogsdHUsSimlKjdjzFJjzMiAgABXh6KUUkoppZRbcUmCSUQ+EZFjIrI11/Q+IrJDRHaLyHMAxpi9xpgHXRGnUkoppZRSSimllCqYq/pgmgVMB2ZnTRART+A9oBcQC2wUkSXGmG2lscO0tDRiY2NJSUkpjc1VOT4+PgQFBeHl5eXqUJRSqkxoOVF6ssoMpZRyJ1WxnNB7AKVUUbgkwWSM+VFEQnJNvgrYbYzZCyAi84B+QKESTCIyEhgJUL9+fSIiIi6YX7NmTerXr0/jxo0RkRLFX9UYYzh9+jRbtmwhMTExe3piYuJFx1mpguj3RlVUsbGx+Pv7ExISouVECRhjiI+PJzY21tWhKKVUqapq5UTO83mzZs1cHY5SqhKoSKPINQYO5HgfC1wtInWAKUBnERlvjJnmbGVjzAxgBkBYWJjJ3at8dHQ0QUFBVaIwKAv+/v4kJiYSFhaWPU1771fFod8b5Uoicjtwe8uWLS+al5KSUmVuGsqSiFCnTh3i4uJcHYpSSpWqqlZO6PlcKVVUFamTb2dnamOMiTfGPGKMaZFXcil7AwWMDlRVCoOyoMdOKeUOCurkW891pUOPo1LKXVW181tV+7xKqZKpSAmmWKBJjvdBwKGibEBHB1JKKaWUUkoppZQqfxUpwbQRuExEmomINzAYWOLimErdokWLEBG2b98OQExMDL6+voSGhtKpUyeuvfZaduzY4XTdL7/8kvbt2+Ph4cGmTZuKtN+QkBAuv/xyQkNDCQ0NZfTo0QDcf//9NGvWLHv6tddeW7IPqJQqO+HhEBICHh72Z3i4qyNSZcBZOSEiTJw4MXuZ48eP4+XlxeOPP+50G7/++iuXX345LVu2ZPTo0RhjCrXvnGVS1mv2bDseR0hICMePHy/hp1NKlSktJ9xezZo1L3g/a9asPMsCpZTKaePGzuzY8Sjnzh0us324JMEkInOB/wGtRSRWRB40xqQDjwMrgWhggTEmqojbzbeJXJGUUQE9d+5cunXrxrx587KntWjRgsjISLZs2cKwYcOYOnWq03U7dOjA119/zQ033FCsfa9evZrIyEgiIyN55513sqe/9tpr2dN//vnnYm1bKVXGwsNh5EjYtw+MsT9HjtSbB1cqx3KiefPmLFu2LPt91gOHvIwaNYoZM2awa9cudu3axYoVKwq9/6wyKes1dOjQ4n0QpVT50nIiTyJSQ0Q+E5GPRGRIeeyzPG7klFKqKM6ejeTIkZmsX9+MqKh7y+T85JIEkzHmXmNMQ2OMlzEmyBgz0zF9uTGmlaO/pSnF2G7pNJErowI6MTGRdevWMXPmzAtuHHI6c+YMgYGBTue1bduW1q1blygGpVQlk5gIK1fCo49CUtKF85KSYMIE18RV1ZVzOeHr60vbtm2za6/Onz+fgQMHOt3G4cOHOXPmDNdccw0iwtChQ1m8eHGJ4lJKVQLPPVelygkR+UREjonI1lzT+4jIDhHZLSLPOSb3BxYaY0YAd5RHfFk3chs2NC+3RFNcXBwDBgygS5cudOnShXXr1gEwadIkXn/99ezlOnToQExMDDExMbRt25YRI0bQvn17evfuTXJycpnHqZRyHWNSMeYccXHzy+T8VJFGkSux/EYHusCYMRAZmff89evh3LkLpyUlwYMPwkcfOV8nNBTeeivf3S5evJg+ffrQqlUrateuzebNm6lduzZ79uwhNDSUhIQEkpKS2LBhQ/7x57Jjxw4GDRrkdF5ERASXXHIJAD179sTT0xOAYcOG8eSTTwLw9NNP89JLLwHQvn17wvVJl1Kuk5QEP/8MERGwejX88gukp+e9/P795RZalVLBygmAwYMHM2/ePBo0aICnpyeNGjXi0KGLuyo8ePAgQUFB2e+DgoI4ePAgAOHh4bz22msXrdOyZUsWLlwIkF0mZXn33Xe5/vrr841bKeVCR47AG29AbKzz+e5bTswCpgOzsyaIiCfwHtAL27/rRhFZgu3b9Q/HYhmlsfNdu8aQmJhPOUHWjRwcPvwBhw9/iLd3fapXD8bDo7rT5WvWDOWyy/IvJ5KTky84R584cYI77rA5s7///e88+eSTdOvWjf3793PzzTcTHR1dwOfYxdy5c/noo48YOHAgX331FX/961/zXUcpVfkYk/vUZ8jMTOHw4Q9JSoqic+c1pbIft0owGWOWAkvDwsJGlGhDuW8aCppeSHPnzmXMmDGAvVGYO3cujz32WHZzBLBPpUeOHFmk5gytW7fOXj8/q1evpm7duhdNf+2117j77rsLvT+lVClKSbHJitWr7WvDBkhNBU9P6NIFxo2Dnj3hoYfgwIGL12/atPxjVuVeTgD06dOHiRMnUr9+/TwfKgBO+1vKGgVoyJAhDBmSf+uQnGWSUqoC278fXn0VPv4Y0tLAz+/iGkzgtuWEMeZHEQnJNfkqYLcxZi+AiMwD+mGTTUFAJPm04BCRkcBIgPr16xMREXHB/ICAABISEgBIS0slI6OwuSoDGFJTD5OefhY/v45Ol0pLS83efl58fX356aefst+Hh4ezefNmEhISWLVqFVu3nq/Qdfr0aQ4dOsS5c+fw8vLK3nZmZiaJiYkABAcH06JFCxISEujQoQM7duy4KIaUlJSLjkVpSUxMLLNtK/el35uiOgnkbiBWDfAE+nD69NBSO55ulWAqtAKeIBMSYps75BYcbGsVFEN8fDw//PADW7duRUTIyMhARHj00UcvWO6OO+5g+PDhAAwfPpzffvuNRo0asXz58jy3XdgaTEqpCiA11SaRsmoo/fyzTUp4eMAVV8Df/24TSt26gb//+fWmTbNNsHLePPj5wZQityau0gpd07UClhPe3t5ceeWVvPHGG0RFRbF06VIAMjIyuPLKKwFbhowaNYrYHDUZYmNjadSoEVC4GkxKqQpu1y54+WWYPRtEYOhQ2zxuwwYtJ6AxkPNpTCxwNfAOMF1E+gJL81rZGDMDmAEQFhZmevToccH86Oho/B1lc7t27+cbSESEZP8u4o2IJw0aDCc4eCLVqzco/Cdywj/H9YGPjw/e3t74+/tjjGHDhg34+vpesHyNGjXw8vLKXi81NTW7s3BfX9/s6X5+fiQmJl6w/ax9dO7cuUQx5yUiIoLcx1mpguj3pvBOnfqRbdseIz39JJmZpX8+yq1qJpgKMmVKqRfQCxcuZOjQoXz44YfZ07p3737BTQDA2rVradGiBQCffvppobZd2BpMSikXSEuDTZvO11Batw6Sk+1NQadOtm+lnj3h+ushv2RwVq2TCRPsU+umTe05qYDaKOpCpVbT1UXlxNixY+nevTt16tTJnubp6XlRGeDv78/69eu5+uqrmT17Nk888QRQuBpMSqkKautWmDoV5s8Hb2945BF4+unzNZSyEudVu5wQJ9OMMeYsMLxQGyjsg4hCbatsb+Ry6927N9OnT+fpp58GIDIyktDQUEJCQrIHiti8eTN//vlnmcahlHI9YzI5cOA19u6dgK9vczp2/C/R0fcTEHBtmZ6P3CrBVGoFQhncyM2dO5fnnnvugmkDBgxg6tSp2f1dGGPw9vbm448/drqNRYsW8cQTTxAXF0ffvn0JDQ1l5cqVhY4hZx9MHTt2zB56OmcfTAC//PIL3t7eRf2ISimw/SX99tv5hNLatbajboDLL4cRI2xC6YYbwNG3TqENGVLVbhQqrnIuJ7K0b98+39Hjsvz73//m/vvvJzk5mVtuuYVbbrml0HHk7oPpgQceYPTo0YVeXylVyjZtsueXxYuhZk0YOxaeegoaOLk50HIiFmiS430QcHFndfkorQcRNWqElvmNXG7vvPMOjz32GB07diQ9PZ0bbriBDz74gAEDBjB79mxCQ0Pp0qULrVq1Kpd4lFKukZZ2gu3bhxEfv4x69e6hdeuPqVatFl26/Fbm+xZnfTVUdmFhYSZrpJ0s0dHRtG3b1kURuYfcx1CrJqricLvvTUYGbNlyPqH0009w5oyd17atTSb17Andu0O9eq6NFRCRX40xYa6Ow9W0nCh70dHRHD161L3+31W5cLtyorh++skmllautDVcR4+2rxw1GMtCZSonHH0wLTPGdHC8rwbsBG4EDgIbgfuMMVFF3baWE+eV5efW/3dVHPq9yduZM78QFXUPqamHadHiTRo3fiy7L87Skl854VY1mJRSqsxlZtpmClkJpTVr4NQpO69VK7j33vMJJWdPl5VSSqm8GAOrVtnE0o8/2gcTL78Mo0ZBrVqujq5CEZG5QA+grojEAi8YY2aKyOPASmzvtZ8UJ7mklFKVjTGGgwffZc+ecXh7N6Jz53XUqtWl3OPQBJNSSoWH593UyRjYtu3ChFJ8vJ3XvDkMGGATSj16QOPGLvsISimlKrHMTFi61JY/Gzfa8uTtt+0Ion5+ro6uQjLG3JvH9OVA3qPjFKA0+2BSSqnykJ5+mh07HiIubiF16txGmzaf4eVVxK44SolbJZi0QFBKFVl4+IWdNe/bZy/ov/8ezp61I4IdO2bnBQfD7befTyi56dDPSimlyklGBnz5pe28+48/7IOLGTPsyHDVq7s6uiqpoD6YjDGl3tykInPH7lSUcicJCZFs23YPycl/0rz5qzRpMhYRD5fF41YJplIbHUgpVXVMmHDhSGAAKSnw6af2CXLv3uf7UWrWzDUxKqWUci9pafDFFzBtGuzaZfvs+/xzGDwYqrnV5blb8fHxIT4+njp16lSJJJMxhvj4eHx8fFwdilIqF2MMhw/PZNeux/HyqkNoaASXXNLN1WG5V4JJKaUKJTPTjsqzbJmtseSMCBw4YH8qpZRSpSElBT75BF55xTbL7twZFi6Eu+4CD9c9cVbn5dciIigoiNjYWOLi4so/MBfx8fEhKCjI1WEopXLIyDjLzp2jOHr0cwIDb6Jt23C8vS91dViAJpiUUlXFmTO249RvvrGvY8fsxXz16nDu3MXLN22qySWllFKlIzERPvwQXn8djhyBa6+Ff/8bbrlFy5oKJr8WEV5eXjTT2sxKKRc6e3YbUVH3kJQUTUjIJIKD/4GIp6vDyqaPSsrZokWLEBG2b98OQExMDL6+voSGhtKpUyeuvfZaduzY4XTdp59+mjZt2tCxY0fuuusuTmWNXKWUcm7PHttJaq9eULcu3H03LFoEN95o+146dgxmzry4A1U/P9vRqqoURKS5iMwUkYWujqU0OCsnRISJEydmL3P8+HG8vLx4/PHHnW5jwoQJNGnShJo1axZ5/zt37uTWW2+lZcuWtG3bloEDB3L06FEiIiIQEZYuXZq97G233UZERESR96FUlXHqFLz0EoSEwLhx0L69HTBi7Vq49VZNLimllCq0o0fD+fXXLqSlxdGx47eEhLxQoZJL4GYJJhG5XURmnD59utS2OSliUqltC2Du3Ll069aNefPmZU9r0aIFkZGRbNmyhWHDhjF16lSn6/bq1YutW7fy+++/06pVK6ZNm1aqsSlV6aWl2U65x42z/Vm0bAljxsDBg/bnmjUQFwdz5sB990GdOna0uBkzbAfeIvbnjBnnR5FTLiEin4jIMRHZmmt6HxHZISK7ReQ5AGPMXmPMg+UZX+fO9uuS+9W5c8m37aycaN68OcuWLct+/+WXX9K+ffs8t3H77bfzyy+/FHnfKSkp9O3bl1GjRrF7926io6MZNWpUdnOQoKAgpmjyVamCxcXZPv6Cg2HiRLjmGvjf/+C77+wgEZpYUkopVUgZGSns2PEw0dF/xd//SsLCIqld+yZXh+WUWyWYjDFLjTEjAwICSm2bk9dMLrVtJSYmsm7dOmbOnHnBjUNOZ86cITAw0Om83r17U83R8WPXrl2JjY0ttdiUqrSOH7cdpQ4eDPXq2c64333XNnF75x1bi2nbNnj1VbjhBuedpw4ZAjExtm+mmBhNLlUMs4A+OSeIfUTzHnAL0A64V0TalX9o9l7R2/vCad7ettVLSeRVTvj6+tK2bVs2bdoEwPz58xk4cGCe2+natSsNGzYs8v7nzJnDNddcw+233549rWfPnnTo0AGATp06ERAQwKpVq4q8baWqhIMH4cknbY2ladOgTx+IjISlS6FrV1dHpwqhLB5YK6VUcSUn7+G3367h8OEZNGnyLJ06/UD16o1cHVaeqmQfTGNWjCHySGShl+8xq0eBy4Q2COWtPm/lu8zixYvp06cPrVq1onbt2mzevJnatWuzZ88eQkNDSUhIIH3aAYwAACAASURBVCkpiQ0bNhS4v08++YRBgwYV9iMo5T6Mga1bbQfdy5bB+vU2MdSggW0C17cv3HQT+Pu7OlJVAsaYH0UkJNfkq4Ddxpi9ACIyD+gHbCvMNkVkJDASoH79+hc17QoICCAhIQGAZ5+tzh9/5P0MJjUV0tM9gfO1ENLTDZs2ZXD99c7XufzyTF55xUl/XznMmzePG2+8kYYNGxIQEMBPP/1EYGAgmZmZ9OvXj9mzZ2c3e6tduzYxMTHZMecl5/z58+fzzjvvXLRM8+bN+fzzz9m8eTPt27d3us2kpCTS09N56qmnmDx5Ml27diU9PZ2kpCSny6ekpJCYmKhN6FSRVcbvjc/hwzSdO5cGK1YgGRkc7dWL/ffdR1LTpnDypK1dqyoFHZVaKVVRxMV9zfbtwxHxpEOHpdSte5urQypQlUwwFSTmVAz7Tp8fWWrNvjUABAcEE3JJSLG3O3fuXMaMGQPA4MGDmTt3Lo899lh2EzmwF/8jR45kxYoVeW5nypQpVKtWjSFay0JVFcnJ9uI8K6m0f7+dfuWVtunBbbfBFVfoCDzurzFwIMf7WOBqEakDTAE6i8h4Y4zT9sPGmBnADICwsDDTo0ePC+ZHR0fj70hMenuDZz5N2n19oX5921evMba1S4MGgq9v3sWqtzf4+3vnOR/sg4gxY8bg7+/PkCFDWLJkCY899hgeHh7cddddTJ06lSZNmnDffffh7e2Nt7d3dsx5yTn/oYce4qGHHsonRm98fHycbtPPz49q1apx8803M3XqVCIjI6lWrRp+fn5Ol/fx8aFmzZrkPs5KFSQiIqLyfG+2b7c1lcLD7UnjwQfhmWdo0KwZDVwdm1JKqUopMzOVPXue4eDBt/H3v4r27Rfg4xPs6rAKpUommAqqaZSTTBbMC6bE+4yPj+eHH35g69atiAgZGRmICI8++ugFy91xxx0MHz4cgOHDh/Pbb7/RqFEjli9fDsBnn33GsmXL+P777xFtv6/c2cGDdrS3ZctsnxXJyVCjhu2w+5//tJ2jFqMJkKrUnJ30jDEmHnikUBvIZ/jpnN4qRDFx+DA0b25HHffxgV9/tRXpiqugcsLb25srr7ySN954g6ioqOzOtjMyMrjyyisBW4a8+OKLee4jPDyc11577aLpLVu2ZOHChbRv3541a9YUGOuECROyH3Yo5fbCw21/Svv32+bXU6bYzrqnToWFC23GefRoGDsWGjd2dbRKKaUqsZSU/URFDSQhYQONG4+mRYvX8PDI/wFlRaJXhuVk4cKFDB06lA8//DB7Wvfu3S/qR2nt2rW0aNECgE8//fSCeStWrOCVV15hzZo1+OUe9Uqpyi4zEzZutAmlb76B336z00NC7BPh226D7t3tnbyqqmKBJjneBwGHirKB0mz60LAhDB9uRx4fPrxkySUoXDkxduxYunfvTp06dbKneXp6ZteCLciQIUPyrf163333MW3aNL755hv69u0L2LKnca6b5t69ezNx4kQOHSrS4Veq8gkPh5EjISnJvt+3D4YNg4wMqFULxo+3g0jUq+faOJVSSlV68fHLiY7+G8ak0a7dl1x66d2uDqnI3CrBVNgn00XxQvcXSmU7c+fO5bnnnrtg2oABA5g6dWp2H0zGGLy9vfn444+dbuPxxx/n3Llz9OrVC7CduH7wwQelEp9SLnHmDKxaZZNKy5fDsWO2mdt118Err9j+lNq109F2VJaNwGUi0gw4CAwG7ivKBkq7nJg4EaKi7M+Syq+cyNK+fft8R4/L8swzzzBnzhySkpIICgrioYceYtKkSQWu5+vry7JlyxgzZgxjxozBy8uLjh078vbbbxMfH3/BshMmTKBfv36F+3BKVVYTJpxPLmXJyICAADsoxCWXuCQsVXbK4n5CKaXyk5mZTkzMP9m/fxo1anSiffsv8fO7zNVhFYsYU/LmXxVNWFiYyRppJ0t0dDRt27Z1UUTuIfcxrFR9JKjy56xJwZAhbAgP5+q4OJtU+vFHSEuzF+i33GJrKfXpA7Vruzp6tyUivxpjwlwdR0FEZC7QA6gLHAVeMMbMFJFbgbcAT+ATY8yU4mxfy4myFx0dzdGjR7WcUEVWIa4vjIHNmyEsj9OliK1564YqSzlR1pyVE6r0VYj/d1XpuNP35ty5w2zbdi+nT6+hYcMRtGz5Np6evq4OK1/5lRNuVYNJKVVBOGtScP/9MG4cVx85Yqe1a2eHcr7tNjvmu/blonIwxtybx/TlwPJyDkcpVVUcOWLLsFmz7IileWnatNxCUkop5Z5OnvyBbdvuJSMjkTZtZtOgwd9cHVKJ6R2dUqr0jR9/cZOC9HQ4dYpdTzzBZWPG2N6RlSpn2vRBKXWR1FRbq3bWLNtcOyMDunaFDz6wzbbHjLmwTPPzs7VylVJKqWIwJpN9+6YQEzMJP7/WhIb+QI0aBXeBUBnomN5KqdJx6hR8/jncdRccOOB8mXPnONi/vyaXSmhSxCRXh1BpGWOWGmNGBgQE5DW/nCNyT3ocVYWX1QRu9Gho1AgGDLBDQT79NERHw//+Bw8/DCNGwIwZEBxsm8UFB9v3+XSWr5RSSuUlNTWO33+/lZiYf3LppfdyxRW/uE1yCbQGk1KqJOLiYPFi+Oor+P57W0upcWPw94eEhIuX1yYFpWLymslM6jHJ1WG4HR8fH+Lj46lTpw6iHcsXmzGG+Ph4fHTER1URHT16vgncH39A9epw5522GXevXuDpefE6Q4ZoQkkppVSJnT69jqioQaSlHadVqw9p2HCE211zaoJJKVU0sbGwaJFNKv30k+3ktHlz25/SgAHQpQvMnXthH0ygTQpKyf7T+10dQqWWXxO5oKAgYmNjiYuLK//A3IyPjw9BQUHs27fP1aEo5bwJ3NVXw7//DYMGQWCgqyNUFYg2pVZKlTZjDLGxb7Jnz7P4+IRwxRX/w9+/s6vDKhOaYFJKFWzvXptQ+uor2LDBTmvXzo4SN2AAdOxomw5kyXrS62QUOSIiyj18dzD+u/G8vO7l7Pcy2R7vF7q/oLWZisAYsxRYGhYWNiL3PC8vL5o1a+aCqJRSpc4Y+O03m1SaMwfi46FhQxg3DoYNAx0xUuUhv3JCKaWKKi3tJNu3Dyc+/j/UrdufNm0+oVo15101uANNMJUjT09PLr/88uz3gwcP5rnnnqNHjx7s3buXffv2ZVeRu/POO/nuu+9ITEx0VbiqKjPG9kGRlVTassVOv/JKmyjq3x/atMl/G9qkoFRkZGYwK3IWs7bMAuC+y+9jzh9zMC9oHzdKKXWRY8fgiy+cN4G76SYdsVQppVS5SUj4laioezh37gAtW75F48aj3a5JXG5uVcqWVpXWV/fvp4u/Pz1zVJleffIkGxMSeKYEfcj4+voSGRnpdN4ll1zCunXr6NatG6dOneLw4cPF3o9SxZLV4enXX9uk0o4dtlbStdfCG2/YpFJIiKujrFK+3/s9T337FL8f/Z1rm1zL4kGLuTroaub8McfVoVVa2vRBKTeUmgrffHO+CVx6ujaBU0op5TLGGA4d+je7dz+Jt3d9QkN/IiCgq6vDKhduNYpcQaMDFVYXf38GbtvG6pMnAZtcGrhtG138/UsjTKcGDx7MvHnzAPj666/p379/me1LqWyZmbBuHYwda/tRCguDV16xHXW/9x4cPAhr18JTT2lyqRztOL6DO+bewU2f38SZc2eYf/d81g5fy9VBVwO2WZwqntIqJ5RSLpbVBO7vf7ejwPXvDxs32vJq2zZYvx4eeUSTS0oppcrUxo2d2bHjUc6dsxVE0tMTiI6+j127HiMw8CbCwn6rMsklcLMaTIU1ZtcuIgtoetbI25ubf/+dht7eHE5Npa2fH5NjYpgcE+N0+dCaNXnrssvy3WZycjKhoaHZ78ePH8+gQYMAuPHGGxkxYgQZGRnMmzePGTNm8K9//atoH0ypwkhPhzVrbE2lRYvg8GHw8rKj50ycCHfcAXXrujrKKik+KZ4X17zI+5vex7eaLy/f+DJ/7/p3fKpdOBqX9rmklKqyjh07Pwrc77+Dt/eFo8BpEzillFLl6OzZSJKStnH06KfUqXMHCQmbSEmJoVmzqTRt+iwiblWnp0BaCuchsFo1Gnp7s//cOZpWr05gKVyw5NdEztPTk27dujF//nySk5MJ0doiqjSdOwfffWeTSv/5j+3s1NcXbrnFdtLdty9ojQ6XSc1I5f2N7/Pimhc5fe40I64YwYs9X+TSGpe6OjSllHI9Z03grroK3n/fNoGrXdvVESqllKrCjEnFGIiLWwBAnTp30KDB/VUuuQRVNMFUUE0jON8sbmJwMP8+dIgXQkIu6JOpLAwePJi77rqLSZMmlel+VBVx9iysWGGTSsuWwZkzUKsW3H67bUrQpw/4+bk6yirNGMOSHUt4etXT7Dqxi94tevNG7zfocGkHV4fmtrQPJqUqkaxR4MLD7YORBg1sE7hhw+xIpkqVAS0nlFIlFR+/jG3bBtO58xpXh1Luql5KrRCykksL2rXjxWbNWNCu3QV9MpWV66+/nvHjx3PvvfeW6X6UGzt92l6I9+8P9erB3XfDypVwzz326W/W6Dr9+2tyycV+O/wbN86+kTvn30k1j2osv285K4as0ORSGdM+mJSqQMLDbf9+Hh72Z3g4HDtG0MKFEBoKV1wBH3wAf/mLrbl04IDtJ1CTS6oMaTmhlCqstLQTF7wX8cbDw5dGjR6hXbv5LorKtapkDaaCbExIYEG7dtk1lnoGBrKgXTs2JiSUqBZT7j6Y+vTpw8svv5z9XkQYN25c8QNX7i88HCZMgP37oWlTmDIFeveGJUvsyG/ffQdpadCwITzwgE0k3XCD9klRgRxKOMQ/fvgHsyJnUdu3Nu/d+h4jrxxJNQ/9GymlqpDwcBg5EpKS7Pt9+2zNpMxMWhoDXbrYwSYGD9YmcEoppSqcxMQ/2Lr1Tsc7Tzw8vGnQYDjBwROpXr2BS2NzJb2jceKZpk0vmtYzMLDETeQyMjKcTo+IiHA6PbGAjshVFePsYnzoUDsSHNinv6NH26RS1672ibCqMJLSknjj5zd4Zd0rpGakMvaasUy4YQKX+Fzi6tCUUqr8Pf/8+fIsS0YG1KrFL2+9xVXDh7smLqWUUqoAcXFfER09jGrVauHrexmBgb2qfGIpiyaYlKosnn324ovxzEzbr9Lq1dC5M4i4JjaVp0yTyZw/5jD++/HEnollQNsBvHLTK7So3cLVoSmlVPkyBjZsgAULbE1cZxISSGrWrHzjUkoppQrBmExiYiaxb9+/qFWrK+3bf0X16o1cHVaFogkmpSqyo0dt07cFC+DgQefLJCTYfipUhbN2/1qeWvkUGw9t5MqGVzKn/xyuD77e1WEppVT5MQY2bYL58+HLL21iydvbjmSanHzx8k5qkSullFKulp5+hujovxIfv5QGDR6gVav38fCo7uqwKpwqlWAyxiBaw6NYjDGuDqHqOH7cjvw2fz5ERNhaSm3aQECA7cQ7N70Yr3D2ntzLs989y8JtC2ns35jZd85mSMcheFTBoUorGh0dSKlyYAxs3mwfjixYADEx4OVl+wz817+gXz87umnOZt9gB5+YMsVlYSullFLOJCXtZOvWfiQl7aJly3dp3PgxzSvkocrc7fj4+BAfH6+JkmIwxhAfH4+Pj4+rQ3FfJ07AzJlw8812GOaHH7aj5Tz/PPz+O2zbZjs7zT3ym16MVyinU07zzKpnaPteW5bvWs7kHpPZ8fgO/tbpbyVOLr26f/9FI1muPnmSV/NqZqKc0tGBlCojxkBkpC23LrsMwsLgzTftA5JPP7U1cpcts30HBgTAkCEwYwYEB9vm3cHB9v2QIa7+JEoppVS2+Pjl/PrrVaSlHadTp+8ICnpck0v5qPA1mESkBvA+kApEGGPCi7OdoKAgYmNjiYuLK9X4qgofHx+CgoJcHYZ7OXUKFi+2T3dXrYL0dGjeHJ55BgYOhE6dLuxTKeuiO/cocnox7nLpmel89OtH/DPin8QnxTMsdBhT/jKFRv6l1ya7i78/A7dtY0G7dnQLCGDt6dPZ75VSyiWMga1bbY3bBQtg1y7w9IQbb4Tx4+HOO6FOnbzXHzJEyzCllFIVkjGG/ftf4c8/n6dmzU506LAYH59gV4dV4bkkwSQinwC3AceMMR1yTO8DvA14Ah8bY14G+gMLjTFLRWQ+UKwEk5eXF82000jlamfOwJIl9mJ85UpIS7NPbZ98EgYNsn0p5ZcR14vxCmfF7hWM/XYs2+K20T24O2/e/CZXNCz9PrF6BgbyeZs29P3jD2p5epIBLGjXrsSjWyqlVJFFRZ1v/rZ9ux21tGdPePppuOsuqFvX1REqpZRSxZaRcZbt2x8kLm4+9eoNok2bT/D09Ct4ReWyGkyzgOnA7KwJIuIJvAf0AmKBjSKyBAgC/nAsllG+YSpVChITYelSeyH+3//CuXMQFARPPGFrKl11lY7+VglFHYti7LdjWblnJS1rt2TRoEX0a92vzKrMfnfiBI/v2kVyZibJmZk806SJJpeUUuVn+/bzSaWoKFtu9egBf/879O8Pl17q6giVUkqpEktOjmHr1js5e/Z3mjd/mSZNntEmcUXgkgSTMeZHEQnJNfkqYLcxZi+AiMwD+mGTTUFAJPn0GSUiI4GRAPXr1yciIqLU41YXSkxM1OOcB4/kZOps2MClq1dTe/16PFNTOVenDnF9+3KsZ0/OtGtnn/gmJ8OaNa4Ot1xV9u/NydSTzIqZxbLDy/Cr5sejLR7lzkZ34nXEizVHSv9veQrbRngVUBeoga3W+eGBAzQ4cIDOpb5HpZRy2LnzfFLpjz9sUun662H6dBgwwPYZqJSb0cEglKq6Tp5cTVTUPRiTzuWXL6dOnT6uDqnSqUh9MDUGDuR4HwtcDbwDTBeRvsDSvFY2xswAZgCEhYWZHj16lF2kCoCIiAj0OOeQnGxrKM2fbzsyTUqC+vVhxAgYOJDq3boR5OFBVe/JqrJ+b1LSU3hnwztMWT+Fs6lneeyqx3ih+wvU8cunf5ESMMYw68gRxu3ZQ0JGBn+99FJWxMeztH17egYGsvrkyew+mLQmk1Kq1OzZcz6pFBlpp113Hbzzjk0qNSq9vuWUqoiMMUuBpWFhYSNcHYtSqnwYYzh48F12734KP79WdOjwH/z8LnN1WJVSRUowOat3ZowxZ4HhhdqAPnFQ5S0lxfaltGCB7VspMdH2PTF0qG3+dsMNtsNTVWkZY1i4bSHPfvcsf576k9ta3cZrvV6jTd02ZbbPnUlJPLxzJxGnTnFdrVp82Lo138TH80CDBtnJpJ6BgSxo146NCQmaYCoCLSeUcuLPP+HLL21Z9uuvdto118D//R/cfbdt1q2UUkq5oYyMFHbtGsWRI7OoU+cO2rb9nGrVark6rEqrIiWYYoEmOd4HAYeKsgF94qDKRWqqHfVt/nz4z39sx921a8Pgwbaj7h49oFpF+tdSxfXLwV94auVTrDuwjo71O7Lqb6u4qflNZba/c5mZvLJ/P1P27cPXw4MPWrViRMOGeIjQvkaNi5bvGRioyaUi0nJCKYf9+21Saf582LjRTrvqKnj9dZtUCtaRcpRSSrm3c+cOsXVrfxISNhAc/E9CQl5AJM9eeVQhVKS74I3AZSLSDDgIDAbuc21ISjmkpcH339unu4sWwalTcMkltrnAwIF2SGYvL1dHqUpoUsQkJvWYxIHTBxj//XjC/winfo36fHT7RwwPHY6nR9nVRvvp1Cke3rmT6KQkBtarx1stW9KwevUy259Syk2Fh8OECTaB1LQpTJlyfvTRAwdg4UJblq1fb6eFhcGrr9qkko62q5RSqoo4ffp/REX1Jz09gfbtv6Jevf6uDsktuCTBJCJzgR5AXRGJBV4wxswUkceBlYAn8IkxJqqI29WmD6p4nF2QDxoEq1fbC/Gvv4YTJ6BWLejXz87r1Qu8vV0duSpFk9dMJiMzg9f/9zrGGJ7v9jzPdXsO/+r+ZbbPk2lpPLt3Lx8dPkxw9ep8c/nl3FqnbPp1Ukq5ufBwGDnS9gEIsG+f7Qdw5UrYuxfWrbPTO3eGadPgnnugRQvXxauUUkq5wOHDM9m581GqVw/iiitWUbNmB1eH5DZcNYrcvXlMXw4sL8F2temDKjpnF+T33w+PPGL7VKpZE+64wyaVevcGHx+XhqtKX6bJ5LPIzwB46aeXuLfDvUy7cRrBl5RdExFjDPOPHWPM7t0cT0tjbFAQk5s1o4b22aWUKq4JE86XZVmSk+Hzz6FjR3jpJVvr9jLtuFQppVTVk5mZxu7dT3Lo0HsEBvaiXbt5eHnVdnVYbqUiNZFTyjWef/7iC/L0dMjMhK++gltuAV9f18SmytwTy59g+sbpF0ybu3Uureq0YlKPSWWyzz+Tk3l01y5WnDhBmL8//+3Ykc7+ZVdLSilVBezdax+QOCMCW7aUbzxKKaVUBZKaGkdU1D2cPr2GoKCxNG/+Mh4emg4pbW51RLWJnCo0Y+C332DuXNsszpnkZOivbXHdVVpGGm/87w0+2vwRl/hcwhu93+DBJQ9iXjBlt8/MTN6KjeWFmBg8RXi7ZUsea9wYT3E2iKZSSuXDGIiKsk24v/46/wRS06blF5dSSilVwSQkbGbr1rtISztG27ZfUL/+EFeH5LbcKsGkTeRUgaKjYd48+9q504725utrk0m56QW529p8eDMPLnmQyCORDGg7gHdveZeG/g15cMmDZbbPX86cYeSOHWw5e5Z+derw7mWX0USbWyqlisIYO+JbVlJp1y5bO+m66+DNN+3vuZvJ+fnZfgWVUkqpKujo0bns2PEgXl516Nx5Lf7+V7o6JLfmVgkmpZyKiTmfVNqyxV6A9+wJ48bZGkorVlzYBxPoBbmbSk5LZvKaybz+8+vUq1GPrwZ+Rf+252upvdD9hVLf55n0dP7x559MP3iQht7efN2+PXfVq1fq+1FKuan0dFi71iaUFi2C2Fj7cOQvf4GxY+3AEw0anF++Xr28R5FTSimlqghjMti7dzwHDrxGQEA32rdfiLd3fVeH5fbcKsGkTeRUtsOH4csvbRO4rKGYu3aFt9+2o+Y0bHh+2awLb70gd2trYtYwYukIdp3YxYOdH+S1Xq8R6Bt4wTKl3efS4rg4Ht+1i0OpqTzWuDFTmjWjVjW3Ou0qpcrCuXPw/fc2qfSf/8Dx43aAiT59YOpUuO02CAx0vu6QIVp+KaWUqtLS0k6ybdu9nDy5kkaNRtGy5Vt4eOjo3+XBre50tIlcFXfihO2Ue948iIiwnXR36mSHYh40CJo1y3tdvSB3W6dTTvPMqmeYsXkGzQOb8/3Q7/lLs7+U6T5jU1J4YvduFh8/TscaNfiqQweurlWrTPdZlYlIDeB9IBWIMMaEuzgkpYouMRH++1+bVPrmG0hIgFq1bDKpf3+bXKpRw9VRKqWUUhXa2bNRbN16Jykp+2jVagaNGmlqoDy5VYJJVUEJCbBkia2ptHKlbUrQsqWtjTR4MLRr5+oIlQst2bGEUd+M4kjiEcZdM47JPSfj5+VXZvvLMIb3Dx5kwp9/km4MrzRvzpNBQXh5eJTZPt2ViHwC3AYcM8Z0yDG9D/A24Al8bIx5GegPLDTGLBWR+YAmmFTlcOIELF1qk0orV9qaS3Xr2oci/fvbZnDVq7s6SqWUUqpSiItbzPbtf8PDowahoasJCLjO1SFVOZpgUpVPSgosX25rKi1bZjvoDgqCMWNsUumKK2w/S6rKOpp4lNErRrMgagEd63dk8aDFdGncpUz3uSUxkZE7dvBLQgI3BwbyfqtWNPf1LdN9urlZwHRgdtYEEfEE3gN6AbHARhFZAgQBfzgWyyjfMJUqosOHYfFim1RavRoyMqBJE3j4YZtU6tYNPD1dHaVSSilVaRiTyb59/yImZhL+/l1o3/5rfHyCXB1WleRWCSbtg8mNpaXZ/ijmzrWdnCYk2I5MH3jAJpWuvRa0lkiVZ4zh898/58mVT5KYmshLPV/imeuewcvTq8z2eTYjg8kxMbx54AB1vLyY07Ytgy+9FNEkZ4kYY34UkZBck68Cdhtj9gKIyDygHzbZFAREAnmeCERkJDASoH79+kRERJR63OpCiYmJepwBn0OHqLt2LfV+/JFa27YhxpDUpAlxgwZx/IYbSGjVyj4YMQZ++snV4bqcfm9UaROR5sAEIMAYc7er41FKlZ709AS2bx/K8eOLqV9/KK1afYinp47U7CpulWDSPpjcTGamvdCeNw8WLrSdnAYEwN13w7332pHgtMNk5RBzKoaHlz3Mt3u+5bom1/HR7R/Rtl7bMt3nivh4Ru3aRUxKCg81bMgrzZtT26vsklmKxsCBHO9jgauBd4DpItIXWJrXysaYGcAMgLCwMNOjR4+yi1QBEBERQZU8zsbAtm22ltLXX0NkpJ3euTO8+CL0749f27YEixDs2kgrpCr7vVFOFbHJtFOOBxMPisjCso5XKVV+kpJ2s3VrP5KSdtCy5Vs0bjxaH/K6mN6dq4rFGNi0ydZUmj8fDh0CPz+44w5bU6lPH+2PQl0gIzOD6b9MZ8IPExARpt8ynVFdRuEhZVej7WhqKmN272besWO08fNjTWgoN1xySZntT2VzdsVgjDFngeGF2oDWdFVlxRjYuNEmlBYtgp07ba2ka6+FN96Au+7Kf7AJpVReZlH4JtOewLRc6z9gjDlWPqEqpcrLiRMr2bZtMOBBp04rCQy80dUhKTTBpCqKrVttTaV582DPHvDygltusTWVbr9dR85RTkUdi+KhpQ+xPnY9t7S8hQ9u+4CmAU3LbH+ZxjDz8GGe2buXpIwMJoeE8GzTplTX5pnlJRZokuN9EHCoKBvQmq6q2MLD7QAS+/dD06YwZYrtjHvt2vNJpdhYW7O2Z0948kno1w8aNnR15KqK6Nz5fGW5nEJD4bffyj+elFWuMQAAIABJREFU0lKUJtPGmGnY2k7Fok2py582iVVFZzh37nN+//0zIAT4F1u2eAIRLo1KWZpgUq6zZ8/5pNLWrbYPpRtvhOeft096AwNdHaGqoFIzUpn20zSm/DSFWtVr8cVdX3Df5feVaZXY6LNnGblzJ2tPn6Z7QAAftm5Na7+yG5FOObURuExEmgEHgcHAfUXZgNZgUsUSHg4jR0JSkn2/bx8MGwaPPAKJieDjAzffbJNOt90GtWu7Nl5VJV1zjW2ZmZp6fpq3t61E54byajLtlIjUAaYAnUVkvCMRdRFtSl3+tEmsKoqMjCR27BjBsWNzqFfvHtq0+RRPT62IUJG4VYJJbxwqGGdPe3v0gAULbBO4jRvtctddB9On276V6td3aciq4tsQu4EHlzxIVFwU911+H2/d/Bb1atQrs/2lZGQwdf9+Xt6/n5qennzSujX3N2ig7bvLmIjMBXoAdUUkFnjBGDNTRB4HVmKbQXxijIkqyna1BpMqlvHjzyeXsmRk2L4CFy60zbe1pq1ysYkT4dNPL5zm6WmnuyGnTabzWtgYEw88UnbhKKVK28aNnalV6xpCQiZSvXpDUlL2s3XrnSQmRgIP0q7dR3o9XgG5VYJJbxwqEGdPe4cOtRfjAFdcAa++apsXNC27Jk3KfZxNPcs/fvgHb294m8a1GrPs3mX0bdW3TPe5+uRJHtm5k53JyQy59FLebNmSS729y3SfyjLG3JvH9OXA8nIOR1VFR47AkiWweDEcOOB8meRkGDCgfONSKg9RUbbGUkqKfe/tDcOHQ4MGro2rjJS4yXRe9IG1UhXD2bORJCVt4+jRTwkMvJnTp3/CmHQuv3wpf/xRQ5NLFZRbJZhUBeLsaW9mph0FbsMGaN3aNXGpSmnVnlWMXDaSmFMxPBr2KNNumkat6rXKbH/xaWmM27OHWUeO0NzHh287dqSXNnlxC3rjoPK1Y4dNKP3nP7B+ve24u3lz8PeHhISLl9cHJKoCOHYMnnrKPtsLCbEJptRUt669BKXQZDov+sBaqYrDmFSMgfj4/wBCvXqDqFnzCmCHq0NTedCeaVXpOXfOXpgPHJj3094zZzS5pArtRPIJ7l98P72/6I23pzc/3v8j7/V9r1STS6/u38/qkycBMMbw+ZEjNF+/ntn/z969x+dcvw8cf7232YbNNuzAMGfbHKep9O1bpIMK6aCUc6SISqmUhEoHyjfKIUQp9KMcQ0kaQpqEzQ7MYexgY+xo5/v9++NDDlEO++yz3buej8f3ofuz+76v6+ux+ey+3tf7eh87xqh69Yho106KS3ZEa71Kaz3Yw8PD6lREWWCzGYser70GQUEQGAijRhmfzt96CyIiIC4OZswwTjQ9X5UqxtZvISxis8Hs2cavVYsXG8Wk6GgYONAYa2kv3UtntkxvA5oppRKUUgO11kXA2S3T0cDiq90yLYQobzTHjy8+c3KcKKukg0lcn+Ji+OUXY6bSd99BRgZ4e4ObmzH49GKy2iuugNaab6O+ZdjaYZzMPcnrt77OmNvH4OrkWuKx2rm782hUFP9r1IgvU1JYf+oUTkoxq2lTBtauXeLxhBAWKygw7ltnO5WSk41Wjw4d4NlnoVu3v9+revUy/rx4ruDZ60KUsshIY8b8li1w220wc6ZRIwWj0LR3r/10L5X2lmnpdBXCeoWFJy94rJQzSjni5zeAgIAxbNsWY1Fm4t9IgUlcPa3h999p/Omn8PjjxpwKd3fj5LcnnjBOgvu//7twBhPIaq+4IklZSQxdPZQVsSu4odYNrOu9jtZ+rU2Ld5unJ494e9MnJgZnpXBzdGR58+Z0kq4luyQfHCqozExYu9YoKq1ZYzyuWtUYzt29O9x//7+fXNqrlxSUhOVOn4a334YPPzSmDsybZxxoeP4oklq1YONG63Is72SLnBDWysmJITKy65lHjjg4OP9VWHJxOduWKQWmssquCkzywcFkUVGwcKHRrXTwILUrVTKOY37iCeOX88qVzz1XVnvFVdJaM2fnHF7+6WXyi/OZeOdERrQfgZODef9M7Tt9mgExMWzNzKRp5crsy83l1Tp1pLhkx+SDQwWSnHxuSPfPP0NhodFh++ijRlGpUydwLfmuSCHM8sMPMHQoHDoE/fvDpElQs6bVWQkhRMk5eXIde/c+ioODM5UrN8HL666LCkuirLOrApN8cDBBfDx8841RWNqzx9jU36kTvPEGW729ubVLl8u/VlZ7xRWKOxnHU6ueIuxwGB3qd2B219k0rm5eobhYa6YmJPD6oUNUdnDgtbp1mX3sGGMCApiRlERHT086/ls3gxCi7ImJMba9LV9uDOkGaNQInn/eKCrdfLOxHU6IciQ5GUaMMJrDmzUzdnh26GB1VkIIUXK01iQmfkJc3AiqVm1By5YrcXUNsDotcQ3sqsAkSkhqKixZYnQqbdliXLv5Zpg6FXr0+GtiZFFYmHU5CrtQZCti8rbJjA0bi4ujC7O6zGJQ20GmHjsad/o0A2Jj+TUjgy41atDX15eh+/ezODiYjl5edPT05NGoqL8eC/sina52xmaD3383CkrLlxunwAGEhsI77xhFpeDgC/cPCVFOFBfDZ58ZM+jz842586+8Ai4uVmdm3+Q+IUTpstkK2b9/GMnJs6hR4wGCgr7GycnN6rTENZICkzBkZhq/nC9cCOvXG7/VNG9ubGvr2dM4plmIErTr2C4GrhzIzuSddA/szrT7plHb3byh2jat+TQxkVEHD+KsFF8GBtLH15dJR49eUEzq6OXF4uBgwrOypMBkh6TT1Q7k5184pPvYMXByMlo6hg83hnTXrWt1lkJcl1274Omnjfppp07GQYZNmlidVcUg9wkhSk9hYRp79z5CenoY9eqNokGDCSglB92XZ1Jgqsjy8oxhp4sWwfffG48DAozlsccfh5Ytrc5Q2KECWwGv//w6E7dMpEaVGizpsYSHgx42tWvpQG4uT8bEsCkjg/uqV2dWs2b4n1kCfuUSJxt29PKS4pIQZUlGxoVDurOyjCHd995rdCndd9+/D+kWohzIzoZx4+Djj6F6dfj6a2PUpTThCSHsTU5ONBERXcnPP0pg4Hz8/PpYnZIoAVJgqmiKioyV30WL4LvvjM4lb28YNMj4Debmm+W3GGGazfGbGbRjEEdzj9K/TX8+uvsjqlc2b6C2TWumJyby6sGDOCnFvGbN6OfnZ2oxSwhxFRYsuPxhEElJ54Z0b9hgDOn28TG6art3hzvukCHdwq6sWgXDhhk/Dk89Be+/bxSZhBDC3qSl/UBU1GM4OLjSpk0YHh7trU5JlBApMFUEWsP27cb2t8WLISUF3N3hoYeMotIddxjbC4QwwbiwcbzY/kVGrR/FjB0z8HP1Y13vddzV6C5T4x7KzeXJ2FjC0tPpXL06s5s2pY58GK3wZLZGGbJgAQwebJy7DsahEoMGGdve4uONvUEAjRvDCy8YRaWbbpIh3cLuJCQYc+iXLjWmE2zeDLfeanVWFZfcJ4Qwj9aahIQpHDjwElWrtjwzzPvvuwlE+SVVBXsWGWl0Ki1aZJxp6+IC999vFJXuuw8qV7Y6Q1EBjN84ns///JzEzERG3DyCu5zuMrW4ZNOamUlJvHLgAA5KMadZM56UriVxhszWKENGjz5XXDorL884ZKJdO6ObqXt3CAqSzlphl4qLYdo040ehqAjefRdeegmcna3OrGKT+4QQ5rDZCti//1mSk+dQs2Z3AgO/kmHedkgKTPbm8OFzRaWICHBwgDvvhDffhAcfBA8PqzMUFcTJ3JM8t/Y5ADxcPPh24LfcVOcmwkw8ffBwbi4DY2PZkJ7OXV5ezGnWjHrStSRE2ZKbaxwmER9/6a8rda57SQg79ccfxhDvP/6Azp2NQpOcpyKEsFcFBSfYu/cRMjI2Uq/e6zRo8LYM87ZTdlVgqhAtrZeaV3HnncaK78KFsG2b8bz27eGTT6BHD/D1tTZnUeH0+q4XCyMX/vV47/G93Pz5zYy9fSwd6FDi8bTWzEpOZuSBAyhgVtOmDKpVS7qWhCgrTpwwDpNYsQLWrTM6l5QytnBf7BKD94WwF1lZMGaM8Suajw/83/8Zv6rJ7UoIYa9ycqLODPNOJCjoa3x9e1mdkjCRXRWY7L6l9VLzKvr2BZvNeNyypdFf3bMnNGhgXZ6iwsrKz+LFH19kYeRCWvi0YH73+bSd1RY99tyHyJLuYIrPy2NQbCzrT52ik6cnnwcGEiBdS0JYLy6OOosXGx20W7YY9yp/f+jfHx54AJKTYejQC7fJValiLJwIYWe0hmXL4LnnjPn1Q4YY3+qenlZnJi5WIRashSglaWlriIrqiYNDFUJCNlKt2k1WpyRMZlcFJrv32mt/n1dhs0G1asYv7y1aWJOXEEDY4TAGrBjAkYwjvPqfVxnfYTwuTi6mxdNaMyc5mZcOHEADM5o04enataVrSQir2GwQHm50Ka1YAVFRNAZo1crovH3gAWjb9sJWDSeny58iJ4SdiI+H4cONU+JatTIO8b1JPmOVWXa/YC1EKTCGef+PAwdexs2tFS1arMTVta7VaYlSIAWmsq64GMLCjO1vR49e+jlZWVJcEpbJLczltZ9fY8r2KTSu3pjNAzZzS91b/vr62NvHlnjMo2e6ltadOsUdnp583qwZ9WVovbgCsjJdwvLyYMMGo6C0apXRleToCLfdBoMH85uPDzc//vjlX9+rlxSUhN0qKoIpU4wmPoBJk4zT4ipVsjYvIYQwk81WwL59Qzh2bC41az5EUNB8HB2rWp2WKCVSYCqLtIYdO4yi0jffwLFj4O4OVatCTs7fny/zKoRFtidsp9/yfsSmxTL8xuG81+k9qjpfeAMZ12FcicXTWjP32DFejIujWGumNWnCM7Vr4yBdS+IKycp0CTh5ElavNopKP/4I2dng5mZMKn7gAeOU0urVAcgzcai/EGXZ9u3GEO/du6FLF/j0UwgIsDorIYQwV0HBcfbufZiMjM0EBLxB/frjZZh3BSMFprIkNtYoKi1cCHFxxjm1998PTzxh/Ll06YUzmEDmVQhL5Bfl89bGt3h/y/vUqVaHn/v+zB0N7jA1ZkJeHk/t28cPJ09yu4cHcwMDaShdS0KUjkOHzm1927zZ6K6tVcvoPnrgAejYEWT2mRBkZMDrr8OMGVC7trEd7sEHZYi3EML+ZWdHEhnZlfz8ZIKCFuLr+w8dzMJuSYHJaomJRpfSwoWwc6fxG0jHjjBqFDz88IXTH89uI5B5FcJCu4/tpu/yvuxJ2cOTbZ5k8j2T8XD1MC2e1povjh1jRFwchVrzSePGDPX3l64lIcyktXF++tmiUkSEcb15c3j1VaOoFBoKDrIqKQQYPzJLlhhb4FJTjZlLb79tjMkUQgh7l5a2mqionjg6uhESsolq1W60OiVhESkwWeHUKWNJa+FCY76S1sYv6pMnw2OPGUtelyPzKoRFimxFfPDrB4zfOJ4aVWqw6vFVdGnaxdSYifn5DI6NZc3Jk/zXw4N5gYE0kq4lIcxRUAC//GIUlFauNBZAHBzg1lvho4+MolKjRlZnKUSZc+iQcSDiDz8Yc+xXrTJ+rRPlk8zqE+LKaa05evQjDh58BTe3EFq0WIGrax2r0xIWkgJTaTl9Gr7/3igqrVkDhYXQtCmMHQuPP278txBlVMyJGPot78fvib/zWPPHmHbfNGpUqWFaPK0181NSeCEujnybjSmNGzNMupaEKHnp6cY9acUKWLvWODSiShW45x6joHT//VCzptVZClEmFRYatde33jJm23/8MTz7rHE4oii/ZFafEFfGZss/M8x7Ht7ejxAY+IUM8xZSYDJVURGsX28UlZYtMwah1qpl9E0/8cTfj2sWooyxaRtTt0/ltZ9fo0qlKnzz8Dc81uIxU2Mm5efz9L59fJ+Wxn+qVWNeYCBNqlQxNaYQFcqRI+e2vm3caNyrfH2NDtoHHoBOnUA6BYX4R1u2GEO89+41ZixNnQp1ZNFeCFFBGMO8HyIj41cCAt6kfv2xMsxbAFJgKnlaw7ZtRlFp8WI4ftyYo9Szp1FUuu02Y5lLiDLu0KlDDFgxgI3xG+nStAuzu87Gz83PtHhaaxakpPBcXBy5NhuTGzXiuTp1cJQirChBFWLrw4IFf5/VFxx8rqi0a5fxvMBAeOklo6h0000yT0mISwgJOfcjc7F69YzdpF27lm5OQghhpezsCCIiulJYmEJw8Df4+Ji7+CzKFykwlZS9e8+dAHf4sHGaTrduRlGpc2dwcbE6QyGuiNaaOTvn8OK6F1Eo5nabS/82/VEmFnqOnelaWpmWxi1nupaaSteSMIHdb31YsODC00bj46FPH2PxQym45RaYONEoKsnWbCH+Vfv2EBVljCg7X+vW8Ouv4OZmTV5CCGGFEydWER39BI6O7rRps4lq1dpZnZIoY8p8gUkp1RAYDXhorR+xOp8LxMefOwFuzx6jM+nOO2H8eOjeXY4OEeVOUlYSg1YOYm3cWu5ocAdzu80lwDPAtHhaaxalpjJ8/35yiov5sFEjXpCuJSGuTXo6vPDCueLSWVpDjRrGp2QfH2tyE6KcGjMG5s698JqLizHQW4pLQoiKwhjmPYmDB0fh5taWli1X4OLib3VaogwytcCklJoLdAFStdYtzrveGZgCOAJztNbvX+49tNYHgYFKqW/NzPWKnThhnEO7cKGxdAXG8tYnn0CPHsYcCyHKGa01iyIXMWzNMPKK8vjk3k8Y2m4oDibupU4pKGDIvn0sO3GCm6tVY16zZgRWlcGAQlyV+Hhjj87KlcappEVFl37eyZNSXBLiKuXmwvTpxjDvs5ydYeBA8DNvx7gQQpQpNls+sbFPk5LyJd7ejxIYOA9HR9lpIC7N7A6mL4BPgflnLyilHIFpwF1AAhCulFqJUWx676LXP6m1TjU5x3+XnW3MrVi4ENatM36BDw425lr07AkNG1qdoRDX7HjOcYasHsJ30d/Rvk57vuz+JU1qNDEtnga+SUlh2P79ZBcXM7FhQ16sW1e6loS4ElrDzp1GQWnFCti927h+dp7SF19ASsrfX1evXqmmKUR5t2YNDBsGhw7BQw8Zj/PyjGb1MWOszk4IIUpHQUEqkZEPkpm5lfr1xxEQ8KapYzNE+WdqgUlrvUkpVf+iyzcCcWc6k1BKfQM8oLV+D6Pb6ZoopQYDgwF8fX0JCwu7qtf7rF9PwzlzcElNJd/Hh4P9+1Ps4YHP+vXU3LoVx7w88nx8SO3Rg5ROnchp2NCYZ3HkiPG/Cig7O/uq/55F2fLriV+ZvG8y2UXZDG4wmEfrPkpiRCKJJJoS7xTwYVERW6OjCQQ+AgIOHmTzwYOmxBPCLuTnG91JK1YYhaXERGMg93/+A5MmGfP+zs5TatnywhlMAFWqGAsiQoh/dfSosdN06VKjbrthA3TsCEOHwmefwYAB0r1k7yrEYRBCXIHs7D1nhnkfJzh4MT4+PaxOSZQDVsxg8geOnvc4Abjpck9WStUAJgAhSqnXzhSi/kZrPQuYBRAaGqo7dOhw5RktWAD/+99fv5C7pqQQ/MEHxtdq1DB+m3jiCVxvuYV6Dg7IOrAhLCyMq/p7FmVGel46z//wPPP3zifEL4T5D86nhU+Lf3/hFZp45Ajt3N3p6OX117Wxhw7x4dGjFDo58V6DBoysWxcnObVKiEs7edJomVi50hj2kpVlFIruuccY0H3ffeDt/ffX9epl/HnxKXJnrwshLqmwEKZOhbFjobjY+LEZOdLYEgdG19LevdK9VBHY/WEQQlyBEydWEBXVCycnD0JCNuPufoPVKYlywooC06V66vTlnqy1TgOeuaI3vtYVh9Gj/z4UFYxf3hMToVKlq3s/Icqwnw78xJMrnyQ5K5k3b3uT0beNxtnRuURjtHN359GoKBYHB9OialV67N3LxowMmlauzKjcXAYEmDc4XIhy6+DBc1vfNm82PuX6+cHjjxtdSp06GSeU/ptevaSgJMRV2LIFhgyBiAi4/35jrGaDBhc+p1Yt2LjRmvyEEKK0aK05cuQDDh16HXf3UFq0WI6LS22r0xLliBUFpgSg7nmP6wBJJfHG17zicLktbidOSHFJ2I3sgmxe+ekVZuyYQVDNIJYNWkZo7VBTYnX08mJxcDDdIyMp0prTNhsD/fyY2bQpv27aZEpMIcodmw127Di39S0y0rjeogW8+qrRqRQaamyHE0KUuBMnYNQo+PxzqFMHli0zfuxkvIgQoiIqLs5j377BpKR8hY9PT5o1m4ujY2Wr0xLljBUFpnCgiVKqAZAI9ASesCCPc+rVM07iudR1IezA5vjN9F/Rn0OnDvFS+5d4u+PbVK5k3g0jvbCQeceOkVlcDMDTtWoxs1kz0+IJUW7k5cHPPxsFpVWrIDnZmBr83/8aW7W7doVGjazOUgi7ZrPBvHlGHTcjA15+Gd58E9zcrM5MCCGsUVCQcmaY9zbq13+LgIA3ZJi3uCamFpiUUouADkBNpVQCMFZr/blSahjwI8bJcXO11ntLKN61bZGbMEGGogq7lFeUxxsb3mDytsk08GrAxv4b+W/Af02N+fOpUwyIiSExP58qDg68UKcOs5KTeczH54KZTEKUFKVUQ2A04KG1fsTqfP7mxAlYvdroVPrxR+Ne4+YG995rbH277z6oXt3qLIWoECIi4JlnYOtWuPVWmDHDaBoUQoiKKitrF5GR3SgsPEHz5t/i7f2w1SmJcszsU+Qev8z1NcAaE+Jd2xY5GYoq7NCOpB30XdaX6BPRPHPDM0y6exJuzuYtz+YWFzPq4EGmJiZSx9mZak5OLG3enI5eXtzp5fXXTCZZCxHnU0rNxThBNFVr3eK8652BKRgLEXO01u9f7j3OnEo6UCn1rdn5XrF9+87NU9q61WiZqFMH+vc3ikodOoCLi9VZClFhZGfDuHHw8cfg6Qlz50K/frIDVQhRsYSHh1CtWnvq1x+Di0stjh9fRnR0b5ycvAgJ+RV397ZWpyjKOSu2yJVNMhRV2InC4kLe2fQOEzZPwM/Njx96/cA9je8xNWZ4ZiZ9Y2KIOX2a4f7++FSqxH88PP7qWDo7kyk8K4sbTc1ElENfAJ8C889eUEo5AtOAuzDm9oUrpVZiFJsuPkn0Sa11aumk+g+Ki2H79nPzlGJijOtt2sAbbxiDXUJCZLiLEKVMa2O20vPPQ0ICDBoE779vHBIshBAVTU7OLk6fjiIlZR5Vq7YiK+t33N1vPDPMu5bV6Qk7YFcFpmveIieEnYhMjaTvsr78eexP+rbuy5TOU/B09TQtXqHNxrtHjvD24cP4OTuzrlUr7rrMVp+OXl509PIi7OBB0/IR5Y/WepNSqv5Fl28E4s50JqGU+gZ4QGv9Hka30zVRSg0GBgP4+voSFhZ2Va/3Wb+ehnPm4JKaSr6PD4f69qWoWjVqbt1Kjd9+w/nUKWyOjqS3aUPac89xon178v38jBdnZlbII6iys7Ov+u9ZiJL6vklKcmXq1CZs316Dhg2z+eSTfbRokUlExPXnKIQQ5ZXWBWgNWVm/A464ubW2OiVhR+yqwHTNW+SEKOeKbcV8tO0jxvwyBk9XT5Y9tozugd1NjRmTk0OfmBh2ZGXR29eXqY0b4yWnLoqS4Q8cPe9xAnDT5Z6slKoBTABClFKvnSlE/Y3WehYwCyA0NFR36NDhyjNasMAYwn1mVp9rSgpBkyYZX/PwMOYodeuGw733Ut3Dg+pAkyt/d7sVFhbGVf09C8H1f9/k58OkSca0AycnmDwZhg93w8lJtn4IIcSFiklO/pzTp2MJCal4C2Gi5NlVgUmIiijuZBz9lvdj69GtPBT0EDPvn4l3VW/T4tm05tPERF49eJAqDg4sCQ7mER8f0+KJCulS+8j05Z6stU4DnrmiN77WTtfRoy88COIsHx84ehScna/u/YQQptiwAYYOhdhYeOQRoy5cp47VWYnyRHZECHuVk3PhuVpKOaOUI35+AwgIGGNRVsLe2NVoQ6VUV6XUrIyMDKtTEcJU48LGYdM2pv0+jdYzWxN1PIoFDy3g2x7fmlpcOpqXx927d/N8XBx3eHoS2a6dFJeEGRKAuuc9rgMklcQba61Xaa0He3h4XN0Ljxy59PXjx6W4JEQZcOyYMUqzUycoLIS1a2HJEikuiat3zfcJIcqwtLQf2LmzPQBKVcLBoTK1ag3ippsO0rTpNFxc/CzOUNgLu+pgki1yoqIYv3E8W45uYf3B9XRu3Jk5XefgX83ftHhaa75OSWH4/v0Uac2spk0ZVKsWSgYWC3OEA02UUg2ARKAn8ISlGdWrB/Hxl74uhLBMcTHMnGk0Gebmwpgx8NprULmy1ZkJIUTZkJDwKXFxz+Pm1gqbrQBPzw4EBIyRopIwhV0VmISwd1prvtrzFQC/JfzGrC6zGNR2kKmFnhMFBTyzbx/fnTjBf6pV48ugIBrJb+6ihCilFgEdgJpKqQRgrNb6c6XUMOBHjJPj5mqt9/7D21xNvGvb+jBhAgwefOE2uSpVjOtCCEvs2AHPPAN//GF0Lk2fDk2bWp2VEEKUDTZbEXFxL5CUNI0aNboRFLQAJyc3q9MSdk4KTEKUE6/89AqTtk7663F2QTaDvx9MYlYi4zqMMyXm9ydOMCg2lpNFRbzfsCEj69bFUbqWRAnSWj9+metrgDUmxLu2TtdevYw/R482tsvVq2cUl85eF0KUmvR0eOMNo6Dk6wuLFsFjj4HcnoQQwlBUlMHevY9y6tQ66tZ9mYYN30MpR6vTEhWAXRWYZCifsFdr96/lqz1f4ezozIQ7JvDyTy+jx1525vF1yyoq4sUDB5iTnEzLqlVZ17o1rdxkxUOUf9d1n+jVSwpKQlhIa1i4EF56yRh/NmwYvP22cZCjEEIIQ27uQSIiupCbu59mzeZQq9ZAq1MSFYhdDfmWoXzC3uQU5DB09VDuW3gf3lW8CX8ViQF0AAAgAElEQVQqnJG3jDQ15ub0dFrv2MHnycm8Wrcu4TfcIMUlYTfkPiFE+RQTA3feCb17Gw2Ev/8OU6dKcUkIIc6Xnv4rf/xxIwUFKbRq9ZMUl0Sps6sOJiHsyfaE7fRZ1oe4k3GMbD+St+94G1cnVwDG3j62xOPl22y8eegQk44epb6rK5vatOFWT88SjyOEEEJcqdxcYzfqxInG2LPp041xaI6y00MIIS5w7Nh8YmOfwtW1Pi1bfk+VKk2sTklUQFJgEqKMKSwu5J1N7zBh8wT8q/mzod8GOtTvcMFzSnrm0u7sbPpERxORk8NTtWrxUaNGuDvJPw/C/shWaiHKjzVrjG1whw5Bnz4waZIxc0kIIcQ5Wts4dGgMR468i6fnHTRv/i2VKnlZnZaooOQTpBBlSOyJWPos60N4Ujh9W/dlauepeLia1/9frDWTjhzhzcOHqe7kxPctW3J/jRqmxRPCatc85FsIUWpSU1146CFYtgyCguCXX6BDB6uzEkKIsqe4+DTR0X05ceI7atV6iiZNpuHgUMnqtEQFZlcFJlmZFuWV1poZO2Ywct1IKleqzJIeS3gk+BFTYx7IzaVvdDRbMzN5uGZNZjZtSk1nZ1NjCiGEEJdTWAhTpsCYMTeiFLz7rjHQW25NQgjxd/n5SUREdCM7eyeNGk2mTp0XUHKcprCYXRWYZGValEdJWUk8ueJJfjzwI50bd2Zut7nUcq9lWjytNbOTk3kxLg4npfgqMJBevr5yQxIVgixECFE2hITArl2X/lr79qdYsKAmDRqUbk5CCFFeZGX9SUREV4qK0mnRYgU1a3a1OiUhADs7RU6I8ubbqG9pOaMlm+I3Mf2+6ax5Yo2pxaXk/Hy6RETw9L593FytGhHt2tHbz0+KS6LCkFPkhCgb2re/dGdS584wYUKkFJeEEOIyjh9fzp9/3opSDrRtu0WKS6JMkQKTEBbIyMug77K+9FjSg0Zejfjz6T8Z0m6IqYWeJamptAgPZ0N6OlMbN2Zd69bUdXU1LZ4QQghxOaNHg9YXXnN1hXnzQNY8hBDi77TWHDkykb17H6Jq1Ra0bfs7bm6trU5LiAvY1RY5IcqDsMNh9Fvej8TMRMbePpbR/x1NJUfzhvGdKixk+P79LEhNpZ27O/MDAwmsWtW0eEIIIcQ/2b0bhgwxZi4pZRSanJ3hySfBzw9iYqzOUAghyhabrYB9+57h2LF5eHs/SmDgFzg6VrY6LSH+RjqYhCgleUV5jFw3kju+vAMXRxe2PLmFcR3GmVpc+unkSVqGh/NNairj6tdnS0iIFJeEEEJYIjMTRoyAtm0hLg4+/hhcXIyvOTrCmDHW5ieEEGVRYWEau3ffzbFj8wgIeJPg4EVSXBJlll11MMnwVlFW7T62m97LehOZGsmQ0CFMumsSVZ3NK/ScLi7m1YMH+TQxkcAqVVjeogWh1aqZFk+I8kLuE0KUPq1hyRKjuJScDE8/DRMmQPXqEBsLn30GAwYY3UtCmEEp1R24H/ABpmmt11mckhBX5PTpWPbsuZ/8/ASCgr7G17eX1SkJ8Y/sqoNJhreKsqbYVszELRNpN7sdJ06fYM0Ta5h+/3RTi0u/Z2YSsmMHnyYm8ry/PztvuEGKS0KcIfcJIUrXvn1wzz3w2GNGAem332DGDKO4BEbX0q23SveSuDyl1FylVKpSKvKi652VUrFKqTil1Kh/eg+t9XKt9VNAf+AxE9MVosScOvUzO3feTHFxJm3abJDikigX7KqDSYiy5HD6Yfou68vmI5t5OOhhZnaZSc0qNU2LV2iz8U58PBPi46nl4sL61q3p5OVlWjwhhBDicnJz4b334IMPjOHdn3xizF1ydLzwebVqwcaN1uQoyo0vgE+B+WcvKKUcgWnAXUACEK6UWgk4Au9d9PontdapZ/77jTOvE6JMS0qaxb59Q6lSJZCWLb+ncuX6VqckxBX51wKTUqoJ8BqQq7V+1vyUhCjftNZ8uftLnlv7HEop5nefT+9WvU09IS4qJ4c+0dHszM6mr68vUxo3xrOSebOdhDif3CeEEOdbuxaGDYODB6FXL/jwQ9n+VtFdz31Ca71JKVX/oss3AnFa64Nn3v8b4AGt9XtAl0vEV8D7wFqt9c5/yHMwMBjA19eXsLCwq0lVXIPs7Gz5e75AMTAT+Ba4idOnx7B9+2HgsJVJlTnyfVN2XUkH01fAeOADAKVUC+AVrXVfMxMTojw6nnOcp79/mmUxy7gt4Dbmd59PgGeAafFsWjM1IYFRBw/i5ujIt82b87C3t2nxhLgMuU8IITh6FF54AZYuhcBA2LABOna0OitRRpT0fcIfOHre4wTgpn94/nDgTsBDKdVYaz3zUk/SWs8CZgGEhobqDh06XGN64kqFhYUhf8+GoqIsoqIe5+TJ1fj7P0+jRh/i4CAbji5Fvm/KriuZweSgtV6LUU5Fax0JtDA1KyHKodX7VtNyRktW71/NpLsmsaHvhhItLk08coRfTp366/GRvDxu2LGDEQcOcFf16kS2ayfFJWGVcnOfUEp1VUrNysjIsDoVIexGYaHRpRQUZHQvvfsu7N4txSVxgZK+T1yqLVxf7sla66la6xu01s9crrj01xvLfUJYIC8vnj///A8nT/5AkybTadLkYykuiXLpSgpMSUqpBpz5R/tMi6mciyjEGTkFOQz5fghdFnXBp6oP4U+FM/KWkTg6OP77i69CO3d3Ho2KYsPJk3x57BhBv//OrpwcRtaty8oWLfA7e9azEKWv3NwnZMi3ECVr82YICYGXX4Y77oCoKHjtNXB2tjozUcaU9H0iAah73uM6QNJ1vN9f5D4hSltGxm/88ceN5OUdoVWrtfj7D7E6JSGu2ZWURV8A5gB+SqkBQGcg8p9fIkTFsD1hO72X9ebAyQO8fMvLvN3xbVyczCn0dPTyYnbTptwbEUGB1jgpxdeBgfSSwRbCenKfEKKCSU2FV16BL7+EevVgxQro1s3qrEQZVtL3iXCgyZmiVSLQE3jiurMUopSlpHxDTEx/XFz8adMmjKpVg6xOSYjr8q8dTFrrwxg3geeAhsBGoI+5aQlRthUWFzL2l7H8Z+5/KCgu4Jd+vzDxrommFZcAfkhLY8j+/RRpowP81bp1pbgkygS5TwhRcdhs8NlnxoylhQuNbqWoKCkuiX92PfcJpdQiYBvQTCmVoJQaqLUuAoYBPwLRwGKt9V4zchfCDFprDh8eT3T041SrdiNt226X4pKwC1e0sfPMP+LfnvlfmaWU6gp0bdy4sdWpCDsWcyKGPsv6sCNpB/1a92NK5yl4uJrXRn26uJiXDxxgelIS9V1cqObkxHB/f2YkJdHJy4uOXl6mxRbiSpWX+4QQ4trt3AlDhsDvv0OHDjB9ujF3SYgrca33Ca3145e5vgZYUwKpXUA+TwizFRfnERv7JKmpi/D17UezZp/h4CCjLoR9uJIZTOWG7JkWZtJa8+nvnxLyWQiHTh3i2x7f8kX3L0wtLv2RlUXbHTuYnpTEwzVrkl1czNLmzXmrQQMWBwfzaFTUBYO/hRBCiJKWkQHPPQft2sHhw/D118YJcVJcEvZIPk8IMxUUpLB7d0dSUxfRoMF7BAbOk+KSsCsyml6IK5CUlcSAFQNYd2Ad9za+l8+7fU4t91qmxSuy2Xj/yBHGx8fjW6kS61u35o+sLJ719/+rY6mjlxeLg4MJz8qSLiYhhBAlTmtYtAhefNGYuTR0KLzzDnh6Wp2ZEEKUP9nZEUREdKGw8DjNm3+Ht/dDVqckRImTApMQ/2LJ3iU8/f3T5BXlMf2+6TwT+gzG4SfmOJCbS5/oaLZlZtLTx4fpTZrgVakSnS5RROooW+SEEEKYICYGnn3W6FQKDYXVq+GGG6zOSgjzyRY5YYa0tNVERfXE0bEaISGbcXeXf1CFfbKrLXJClKT0vHT6LOvDo98+SuPqjdn1zC6GtBtiWnFJa83nycm02bGDqJwcFgQFsSg4GK9KlUyJJ0RFpJTqqpSalZGRYXUqQpRJp0/D6NHQqpUxc2nGDPjtNykuiYpDtsiJkqS1JiFhChER3ahcuSk33PC7FJeEXZMOJiEu4ZdDv9BveT+SspIYd/s4Xv/v61RyNK/Qc7yggMH79rH8xAk6enryRWAg9VxdTYsnREWltV4FrAoNDX3K6lyEKGtWrYLhwyE+Hvr2hYkTwdfX6qyEEKJ8stkKiYt7jqSkmdSs+SBBQV/h6FjV6rSEMJUUmIQ4zxsb3iC3MJfJv02mSfUmbB24lRv9bzQ15uq0NAbGxHCqqIgPGzViRJ06OJi4BU8IIYQ4X3y8McR75UoIDoawMLj9dquzEkKI8quwMJ2oqB6cOrWeevVG0aDBBJSSzUPC/kmBSYgzdh3bxYTNEwAYGjqUiXdNpKqzeasMOcXFjDxwgJlJSbSsWpWfWrempZubafGEEEKI8xUUwOTJ8NZboBR88AGMGAGyM1tUZDKDSVyt8PAQqlVrT/36Y3BxqcXp03FERnYlN/cAzZrNo1at/lanKESpkQKTqPCKbcV8tO0j3tjwBgBrnljDvU3uNTXm75mZ9I6OJi43l5F16/JOgwa4OMiqhhBCiNIRFmacChcdDQ8+CB9/DPXqWZ2VENaTrdTiauXk7OL06ShSUubh5XUP6ekbUcqB1q3X4+l5m9XpCVGq5BOtqNDi0+NpNLURr65/lUJbIQD3LbwPNV4xLmxciccrstl46/Bhbtm5kzybjQ2tWzOpUSMpLgkhhCgVKSnQpw907Ah5efD997B0qRSXhBDiemhdgM2WR1raCoqLM6hevTOVKzexOi0hSp10MIkKa8GeBQxdMxSbtjHvgXn0a90Ph7cc0GO1KfH2nz5Nn+hotmdl0cvHh0+bNMFT9iEIIYQoBcXFMHOmcULc6dPwxhvw2mtQpYrVmQkhhL3RpKZ+Q35+AiEhG61ORohSJQUmUeGcyj3F0DVD+SbyG26pewtfPfgVDb0amhZPa83s5GRGxMXh7ODAN8HBPObjY1o8IYQQ4nzh4TBkCPzxB3TqBNOmQbNmVmclhBDlX3Hx6QseK+WMUo74+Q0gIGCMRVkJYZ1yUWBSSnUH7gd8gGla63UWpyTKqV8O/ULf5X05ln2Mdzq+w6u3voqTw7kfg7G3jy3ReCkFBTwVG8uqtDQ6eXryRWAgdVxdSzSGEEIIARASArt2Xfprfn6waBE89pgx0FsIcWky5Ftcqfz8ZCIjHzjzyBEHB+e/CksuLn6W5iaEVUwf/KKUmquUSlVKRV50vbNSKlYpFaeUGvVP76G1Xq61fgroDzxmYrrCTuUX5TNy3Ug6ze9ElUpV2PrkVkbfNvqC4hLAuA7jSizmqhMnaBkezrqTJ/m4cWPWtW4txSUhTKCU6q6Umq2UWqGUutvqfISwSvv24Oz89+stW0JMDPTsKcUlIf6N1nqV1nqwh4eH1amIMiw7ew87d95ETk4Urq4NqF37aW666SBNm06T4pKo0Eqjg+kL4FNg/tkLSilHYBpwF5AAhCulVgKOwHsXvf5JrXXqmf9+48zrhLhiESkR9F7Wmz0pe3jmhmf48O4Pqepc1bR42UVFvHjgALOTk2nj5sYvbdrQvKp58YQoz5RSc4EuQKrWusV51zsDUzDuC3O01u9f7j201suB5UopL+BDQLpcRYU0ZgzMnXvhNRcXWLcO5LOyEEKUjLS01URF9cTR0YOQkM24u4dYnZIQZYbpBSat9SalVP2LLt8IxGmtDwIopb4BHtBav4fxQeMCSikFvA+s1VrvvFQcpdRgYDCAr68vYWFhJfV/QVxGdnZ2mf57tmkbSxOXMuvgLNyc3Hi3xbu0d2tP+NZw02JGAe8CScDjQP/sbI6HhxNmWsTyp6x/34hS9wWyCCHEdcvJgY8/hoKCc9ecnWHgQGN7nBBCiOujtSYx8RPi4kbg5taGli1X4uLib3VaQpQpVs1g8geOnvc4AbjpH54/HLgT8FBKNdZaz7z4CVrrWcAsgNDQUN2hQ4eSy1ZcUlhYGGX17zkhM4H+y/vz86Gf6dq0K3O6zcGnqnmDtQttNt6Jj2dCfDx1XFwICwriNk9P0+KVZ2X5+0aUvtJahBDCXmkNK1bAc8/B0aPGNrjlyyEvDxwdja4mIYQQ18dmKyIu7nmSkqZTs2Z3goK+xtFRdigIcTGrCkyXmgBw2bPhtdZTgan/+qYylE8Ai/cu5pnvnyG/OJ9ZXWYxqO0glIlDJ/adPk3v6GjCs7Lo6+vL1CZN8HAqF/PzhSirSnwRAqTT1QrSsWiupCRXPvmkCb/9VoOGDbOZOnUfLVtmcvp0E1atqs3ddycRE7OfmBirM7068n0jhChLiooyiIrqycmTP1C37ss0bPg+Spk+yliIcsmqT8EJQN3zHtfB2FV0XbTWq4BVoaGhT13ve4nyJyMvg+Frh/PVnq+40f9Gvn7wa5rUaGJaPK01nyUl8eKBA1R2cGBJcDCP+JjXJSVEBWLKIoR0upY+6Vg0R34+TJoEEyaAkxN89BEMH+5GpUptAWjWDNLTYeZMf/z8yt/2Dfm+EVaSBWtxvtzcw0REdCE3N5amTWdTu/Ygq1MSokyzqsAUDjRRSjUAEoGewBMW5SLswOb4zfRZ1oejmUcZe/tYRv93NJUcK5kW71h+PgNjY1lz8iR3e3kxLzCQ2i4upsUTooIxZREC5IODKP/Wr4ehQ2H/fujRAyZPhjp1LnxOrVqwcaM1+QlR3smCtTgrI+M3IiMfQOsCWrX6AS+vTlanJESZZ3pvn1JqEbANaKaUSlBKDdRaFwHDgB+BaGCx1npvCcTqqpSalZGRcb1vJcqJguICXv/5dW7/4nacHJz4dcCvjOswztTi0vLjx2m5Ywcb0tP5pHFjfmjVSopLQpSsvxYhlFLOGIsQK0vijeX4aVFeJSUZ85XuusuYu/TDD7B48d+LS0IIIa5fSso37NrVAUdHd0JCtklxSYgrVBqnyD1+metrgDUlHEtWHCqQ6OPR9F7Wm53JOxkUMoj/df4fbs5upsXLKirihbg45h47Rls3N74OCiKoqgz3E+J6nFmE6ADUVEolAGO11p8rpc4uQjgCc0tiEeJMPOlgEuVKURFMm2YM6y4ogPHj4ZVXwNXV6syEEML+aK2Jj3+Hw4ffxMPjvzRvvhRn55pWpyVEuSGTiEW5o7Vmevh0Rv40kqqVqrLssWV0D+xuasytGRn0iY7mcF4er9erx9j69XF2kOF+Qlyv0lyEOPO+shAhyo1t22DIENi9Gzp3hk8/hUaNrM5KCCHsk82WT2zsIFJSvsbXtw/Nms3GwUF2KQhxNeyqwCQr0/bvWPYxBqwYwA9xP9C5cWfmdptLLfdapsUrtNkYf/gw7x05QoCrK5tCQviPbK0RotyS+4QoD9LSYNQomDPH2AL33Xfw4INg4oGoQghRoRUUnGDv3gfJyPiVBg3eoV691009hVoIe2VXLRgyW8O+LYteRovpLQg7HMa0+6ax5ok1phaXYnJyaL9zJxOOHKG/nx+7QkOluCREOSf3CVGW2Wzw+efGKXDz5sHIkRAdDQ89JMUlIYQwS05ODDt33kRW1g6Cg/+PgIDRUlwS4hrZVQeTsE9Z+VmM+HEEn//5OW1rteXrB78myDuoxN5/4pEjtHN3p6OXF2BswXshLo7pSUl4ODqytHlzHvT2LrF4QgghxMV27za2w23bBrfeCtOnQ8uWVmclRMUjna4Vy6lTPxMZ+TAODi60aRNGtWo3WZ2SEOWaXXUwySly9mfb0W20+awNc/+cy2u3vsa2gdtKtLgE0M7dnUejovjl1CmS8vO58Y8/mJqYSFs3NyLatZPikhB2RO4ToqzJzIQRI+CGG2D/fvjiC9i0SYpLQlhFOl0rjqSk2ezZ0xlX17q0bbtdiktClAC7KjDJDcF+FBYXMvaXsdw671aKbcVs7L+Rdzu9i7Ojc4nH6ujlxeLgYLpHRtLot9/YkZ3N8/7+/Na2LbVcZLCfEPZE7hOirNAaFi+GoCCYMgWeegpiY6FfP9kOJ4QQZtK6mLi4kezbNxgvrzsJCdlC5cr1rU5LCLsgW+REmbM/bT+9l/Xm98Tf6du6L1M7T8XD1bwPg5lFRcxPSSGzuBiAobVr83GTJqbFE0IIUbHt2wfDhsFPP0FICCxdCjfJwrkQQpiuuDiHqKhepKWtwN9/GI0a/Q8HB/lILERJkZ8mUWZorZm9czYjfhyBi6MLix9ZTI/mPUyNuSUjgz7R0RzOy6OKgwMv1KnDrORkHvH2/msmkxBCCFEScnPhvffggw/A1RU++cSYu+ToaHVmQghh//LzE4mI6Ep29m4aN55KnTrDrU5JCLtjVwUmGcpXfqXmpDJo5SBW7VvFnQ3v5IsHvsC/mr9p8QptNt6Kj+fd+Hh8KlXCw8mJpc2b09HLizu9vHg0KorFwcFSZBLCzsh9QlhlzRqja+nQIejVCz78EPz8rM5KCCEqhqysnUREdKW4OJOWLVdRo8Z9VqckhF2SGUzCcqv3rabljJasO7CO/93zP37s/aOpxaV9p0/znz//5J34ePr5+THE3/+v4hKcm8kUnpVlWg5CCGvIfUKUtqNH4aGH4P77wcUFNmyAr7+W4pIQQpSWEydW8Oef/0UpJ0JCtkpxSQgT2VUHkyhfcgpyGLluJDP/mEkr31b83PdnWvi0MC2e1prZycmMiIvD1cGBJcHBPOLjc8nndvTyku4lIYQQ16ywED7+GMaPB5sN3n0XXnoJnEv+rAohhBCXoLUmIWEyBw68jLt7O1q0WIGLi1T3hTCTFJiEJcITw+m9rDf70/Yzsv1I3rnjHVyczDuxLbWggEGxsaxKS+MuLy/mBQbiLyfECSGEMMGmTTB0KOzdC926GafE1a9vdVZCCFFx2GyF7N//LMnJs/H27kFg4Jc4Ola2Oi0h7J4UmESpKrIV8f6v7zN+43j83Pz4ue/PdGzQ0dSYq9PSeDImhoyiIj5u3Jjh/v44yBnQQlRIMoNJmCk1FV5+GebPh4AAWLHCKDAJIYQoPYWFp9i7twfp6T9Tr95oGjR4C6XsajKMEGWWXRWY5IND2Xbw1EH6LOvD1qNb6dmiJ9Pvm45XZfO2oZ0uLmbkgQPMSEqiVdWq/Ny6NS3c3EyLJ4Qo+7TWq4BVoaGhT1mdi7AfxcUwaxa8/jrk5Bh/jh4NVapYnZkQQlQsubkH2LPnfvLyDhIY+CV+fn2tTkmICsWuCkzywaFsGvvLWBp4NWD42uE4KkcWPLSAJ1o+YWrMHZmZ9I6OZl9uLiPr1uWdBg1wcZCVCyGEECXrjz9gyBAID4eOHWH6dAgMtDorIcS1kgXr8is9fTORkQ8Cmtat1+PpeZvVKQlR4dhVgUmUPWmn03hr01sA3B5wO/MfnE89j3qmxSvWmg+OHGHs4cP4OTuzvnVr7pBh3UIIIUpYejq88YZRUPLxgQUL4PHHQXZgC1G+yYJ1+XTs2NfExg7E1bU+LVuupkoVKRAKYQUpMAnTrDuwjv7L+wMw8c6JvNj+RRwdHE2Ldyg3lz7R0WzJzORRb29mNm2KV6VKpsUTQghRMYSEwK5dl/7a8OHw9tvg4VG6OQkhhACtbRw+PJb4+Hfw9OxI8+bfUqlSdavTEqLCkj1DosTlFeVx85ybuefre0jOTgbglfWv4PS2E+PCxpV4PK01848do/WOHUTk5PBVYCDfBAdLcUkI8TdKqa5KqVkZGRlWpyLKkfbtwdn5wmtKQY8eMHWqFJeEEMIKxcW5REU9QXz8O/j5DaRVqx+kuCSExaTAJErUnpQ9hM4KZXvidobfOJzTr58GQI/V6LGacR3GlWi8k4WFPBYVRb+YGNq4ubE7NJTefn4o2aMghLgErfUqrfVgD6kIiKvw0ktgs114zdXVKC4JIYQofQUFKezefQfHjy+mYcMPaNZsNg4Ozv/+QiGEqexqi5wM5bOOTdv4+LePee3n16heuTpre62lc+POpsZcf/Ik/WJiSC0s5L0GDXi5Xj0cpbAkhBCihGgNK1fCc89BURE4OBiFJmdnGDAA/PyszlAIISqe7OxIIiK6UFiYSvPm3+Ht/aDVKQkhzrCrDiZZmbZGQmYCd391Ny+te4l7G9/Lnmf2XFBcGnv72BKNl1dczEtxcdy1Zw/VnJzY3rYtowICpLgkhBCixBw6BN26QffuUK0aLF16bpucoyOMGWNtfkIIURGlpf3An3/egtYFtGmzSYpLQpQxdlVgEqVvyd4ltJrRim0J25jddTbLHluGd1XvC55TktviIrKzuXHnTiYnJPBs7dr8ccMNtHV3L7H3F0IIUbHl58OECRAcDL/8Ah9+CDt3woMPGl1LDg7SvSSEEKUhPDyE2Nih5OcbM10TE6cREXE/lSs3om3b36lWLdTiDIUQF7OrLXKi9GTmZ/J+zPv8uPFHbvS/ka8f/JomNZqYFs+mNVMSEhh18CBeTk6sbtmS+2rUMC2eEEKIimf9enj2Wdi3Dx55BP73P6hT59zXx4yBvXule0kIIUpDTs4uTp+OIiVlHpUrNyEnJ4IaNboSFLQQJyc3q9MTQlyCFJjEVdtyZAu9l/XmSPoRxtw2hjG3jaGSo3kntiXk5dE/Joaf09N5oEYNZjdrhvfFx/kIIYQQ1ygpyRjk/c030KgRrF0LnS8xRrBWLdi4sfTzE0KIikrrArSGnJwIwBFn59oUF2dJgUmIMkq2yIkrVlhcyJgNY7jti9tQKKa0mcJbHd8ytbi0JDWVVjt2sC0zk9lNm7KsRQspLgkhhCgRRUUwZQoEBsKyZTBuHERGXrq4JIQQwmrFJCfPJiqqp9WJCCEuQzqYxBXZl7aP3kt7E54UzoA2A/i488fs3LbTtHiZRUUM37+f+Skp3OjuztdBQTSpUsW0eEKIikFOGxVnbdsGQ4bA7t1wzz3w6acg3xZCCFE2ZGbuuOCxUs4o5YScTGYAACAASURBVIif3wACAmSfshBllRSYxD/SWjN752xG/DgCF0cXlvRYwiPBj5ga89f0dPrExHAkL4+xAQGMDgigkoM02wkhrp/WehWwKjQ09CmrcxHWSEuDUaNgzhzw94dvv4WHHgI5iFQIIcqG48eXEh3dGwClKqGU01+FJRcXOWFBiLJMCkziso7nHGfQqkGsjF3JnQ3v5IsHvsC/mr9p8QpsNsYfPsz7R45Q39WVX0NCaO/hYVo8IYQQFYfNBvPmwauvQno6jBwJb74JchCpEEKUDVprjh6dxMGDr1Kt2s0UFWXh6Xm7FJaEKEfsqsAkWx9Kztr9axmwYgCn8k4x+e7JPH/z8zgo87qIYnJy6B0dzR/Z2Tzp58fHjRvj7mRX355CCCEssnu3sR1u2za49VaYPh1atrQ6KyGEEGfZbIXs3z+U5OQ5eHs/RmDgPBwdK1udlhDiKtnVviOt9Sqt9WAP6Xq5ZrmFuQxfM5z7Ft6Hd1Vvdjy1gxHtR5hWXNJaMyMxkbZ//MGhvDy+a96czwMDpbgkhBDiumVlwYsvwg03wP79RgfTpk1SXBKiolFKBSmlZiqlvlVKDbE6H3GhwsJT7NnTmeTkOQQEvEFw8EIpLglRTsmnePGXP5P/pNfSXkSfiGbEzSN4t9O7uDq5mhYvpaCAgTExrD55kru9vJgXGEhtFxfT4gkhhKgYtIYlS2DECEhOhsGD4d13oXp1qzMTQlwtpdRcoAuQqrVucd71zsAUwBGYo7V+/3LvobWOBp5RSjkAs01OWVyVRP78cwi5uQcIDPwSP7++VickhLgOUmASFNuK+WjbR7yx4Q28q3qzrvc67mp0l6kxV504wcDYWDKLipjauDHP+vvjIBNWhRBCXKd9+2DYMPjpJwgJgaVL4aabrM5KCHEdvgA+BeafvaCUcgSmAXcBCUC4UmolRrHpvYte/6TWOlUp1Q0Ydea9RBmQkbEFeJaCAkdat16Pp+dtVqckhLhOUmCq4I5kHKHvsr5sjN/Iw0EP81mXz6hRpYZp8XKKi3kpLo7PkpNpXbUqv7RpQ/OqVU2LJ4QQomLIzYX33oMPPgBXV/jkE2PukqOj1ZkJIa6H1nqTUqr+RZdvBOK01gcBlFLfAA9ord/D6Ha61PusBFYqpVYDCy/1HKXUYGAwgK+vL2FhYSXxf0Fc0npgIjabN0VFH7Brlw0IszgnUV5kZ2fLz2cZJQWmCmxRxCKGrB5CsS5m3gPz6Ne6H6qEuogmHjlCO3d3Onp5/XVtRmIiYw4d4mRREa/UrctbDRrg4mBXY8CEEEJYYM0ao2vp0CHo1Qs+/BD85MAhIeyZP3D0vMcJwGV7FZVSHYCHABdgzeWep7WeBcwCCA0N1R06dCiBVMX5tNbEx7/F4cMT8PC4nYyMF+nQoZvVaYlyJuz/27vz8Kqqc4/j30UgDEmAMAZknjJACKRBpK0CtlRbRcGrCILIUKkoXOeBKgLWK1Zt9aoIpaIIFZVrQQFpoSoRwaFhiBASAsgYhsgQCAFCpnX/CFpEwOScs7PP8Ps8Tx49++yz1uvJ2q7sd68hNRVdn/5JCaYQdLTwKOOWjuPNjW/Sq0Uv5g6cS/sG7X1aR4+oKAZlZjI/IYHL69VjTHY2r+fm0rhGDT5OSqLPWYknERERT+zeDffcAwsXQlwcfPwx9O3rdlQiUgXO90TUXuhka20qGh7jutLSQrKzf8s337xJ06a3ERs7k5UrP3M7LBHxISWYQszKXSu5deGt7M3fyxN9nmDC5ROoXs33zaBvdDTzExK4cdMm6lSrRk5REVfWr8/fO3emfo0aPq9PREKXMSYeuBtoBHxkrZ3uckjisOJieP55mDKlfEHvqVPLd4sLD3c7MhGpIjlAy7NetwD2+aJgY0x/oH+HDh18UZycUVR0kIyMgeTnr6Zt26do1eoRn82cEBH/oflJIaKotIgJH06gz+w+hIeFs3rUaib2nuhIcgnODH8tLOR4aSk5RUUMbNSIj7p1U3JJRL7HGPOaMeYbY0zGOcevNsZkG2O2GWMeuVgZ1tosa+0dwCAgxcl4xX0rV0K3bvDww9CvH2RmwiOPKLkkEmLSgI7GmLbGmHBgMLDIFwVbaxdba8fUq1fPF8UJcOLEZtatu4yCgrUkJMyndesJSi6JBCklmELA5kOb6TWrF0+vfprR3Uez/nfr6dnCuS11jhQXMygzk5HZ2QD89yWX8OmxY6zIy3OsThEJWLOBq88+cNbuQL8GEoAhxpgEY0yiMWbJOT9NznzmOmAV8FHVhi9VJTcXbrsNeveGkydh0SJ47z1o08btyETEScaYt4DPgVhjTI4xZrS1tgQYBywDsoD51tpNbsYp55eX9xHr1/eitLSAbt1SadLkJrdDEhEHaYpcELPWMmPNDO5ffj91atRh4c0LGRA3wNE6PzxyhNs2b+ZAURER1arxXpcu/LJBAwY0avTdmkx9tf6SiJyh3YGCmy92eSkthSVLmvPqq20pLAxj2LA9DB26i1q1ytCvMDhpdyA5m7V2yAWOL+UiC3Z7SlPkfGf//lls2XIHtWvHkpi4hNq127gdkog4TAmmIJVbkMvoRaP5YOsHXNX+Kl6//nWaRTVzrL7C0lIe3bGDP+fkEFu7NjdecgkDGjX6Lpn07ZpMacePK8EkIj9GuwMFCW93eVmzBsaOLf/nlVfCtGkQF9caaO2zGMX/aHcgcZO1djGwOCUl5Xa3YwlU1paxffsE9ux5hujoX9G583yqV9eUQ5FQ4PcJJi3eWnlLtixh1PujyD+dz4tXv8hdl95FNePcbMiMggJuycpi44kT3Nm8Oc+2b0+dsLAfnNc3OlrJJRGpCMd2B9KT6cCQlwePPQbTp0PTpjBvHgweDFqyQ0TEv5WWniQr61YOHVpA8+Zj6dDhRao5tOariPgfR9dg0uKtVetE0QnGLhlL/7f60zyqOWvHrGV8z/GOJZfKgP/NySFl7Vpyi4pYkpjItE6dzptcEhGpBMd2B9Lirf7NWpg7F+LiYMYMGD8eNm+GIUOUXBIR8XenT+8nPb03hw4tpH375+nYcZqSSyIhxukrfjbwMjDn2wNnLd7aj/KbiDRjzCIgDJh6zudHWWu/ObN46yNnypLzWLNvDUMXDGXr4a080OsBnrzySWpWr+lYfftOn+ZhYM22bVzbsCGzYmNpoi18RMQ3vtsdCNhL+e5At/iiYI1g8h/du0N6+vnf69kT/vnP8nNERKqS+gnPFBRsYOPGaykuPkKXLu/TqFF/t0MSERc4mmDS4q3OK7WlvLX7LWbvmk2D8AY81/U5ksOT+XzV547V+SnwHFBoLfcaQ//Dh8n87DMyHatRgokWb5WzndkdqA/QyBiTA0yy1s4yxny7O1AY8JqvdgfS2hr+o1cvyMyEoqLvH+/dGz7+GKppn1sRcYH6ico7fHgpmZk3ExZWj+7dPyUqSk8HREKVG2MWtXirD0xOncyIbiO4deGtrNq9ikGdBzHjmhlE13ZujaOCkhLu3raN1w4cIDkykrsLChge5N+z+J4Wb5WzaXeg0PXYYzBr1veP1aoFb7+t5JKISKDIyXmZbdvuJjIyicTExdSseYnbIYmIi9xIMDm2eGuosNYy5ZMpPP/F8wDMHTiXoYlDMQ4uUPFlfj5DMzPZXljIhFatmNymDZ+tXOlYfSIiTtCTaf+wfTv8939/f/RSeDiMGgUxMe7FJSIiFVNWVsLXX9/H3r0v0bDhdSQkzCMsLMLtsETEZW48I3Rs8VZjTH9jzMxjx475oji/lHcqjyF/L3/g37VpV7664yuGdR3mWHKppKyMJ3bu5Gfr1lFiLZ9068ZT7doRrsfLIiJSSadPw5NPQufO8MknMGlS+aglgLAwmDjR3fhERELhfsJbJSXHyci4nr17X6JFi/vp0mWBkksiAriTYPpu8VZjTDjli7cu8kXBwb470Ij3RtDgmQa8s+kdAFbtXkXb/23L5NTJjtS3/dQprkhPZ9LOnQxu0oSvevTg8vr1HalLRKQq6MbBPf/6FyQmlieR+veHrCyYPBlGjiyfEjdypEYviYj7gv1+wluFhbtZv/7nHDmyjE6dZtChw3OU7+EkIuLwFLmqXrw1WNfWKCotYuLHE5nz1Rw6NezE3wb+jUtfvRQ76YIzC71ireWNAwcYv20bYcC8+HiGNG3qSF0iIlVJU+Sq3r59cN998M470KFD+e5wV131n/cnToRNmzR6SUTE3+XnryEjoz+lpSfp2vUfNGjQz+2QRMTPOL2LXJUu3hqMNw5ZB7MYumAo6w+sZ0zyGP581Z+JCHduCOrh4mLu2LKFdw8epHe9esyJj6fVt/MXREREKqikBN59twVz5pSvtTRlCjz00H+mxH2rWbPy6XIiIuK/Dh5cQFbWMMLDm5KU9BEREQluhyQifsiNRb6lAqy1TF8znfuX309keCTv3fwe18dd/937k3pP8nmdHx45wm2bN3OwuJg/tmvH/S1bEubgwuEiIhKcPvsMxo6FDRs6cPXV8PLL0L6921GJiPy4YJ0R4SlrLXv2PMv27Q9Tt+5ldOnyPuHhTdwOS0T8VFCt1Bwsa2vkFuRy7VvXctfSu+jTpg8bx278XnIJYHKfyT6rr7C0lPu2baPfhg3UrV6dL5KTeahVKyWXRCToBEs/4a8OHYLf/hZ+9jM4cgSmTMlg6VIll0QkcGgNpv8oKysmO/t2tm9/mMaNbyYp6WMll0TkooIqwRQMHcKSLUtInJ7Ixzs+5qVfv8TSW5YSE+ncqqcZBQVcum4dz+fkcFfz5qz9yU9IjopyrD4RETcFQz/hj8rK4NVXITYW3ngDHnywfBHvK644hJ5ViIgEnuLiPDZsuJoDB2bRuvVjJCTMIyysttthiYif0xQ5P3Gy+CT3L7ufGWtnkNQ0iTdveJPOTTo7Vl+Ztby0dy8Pf/019apX54PERH7TsKFj9YmISHBKTy+fDvfFF3D55fDKK9Cli9tRiYiIp06d+poNG66hsHA7cXFvEBMz3O2QRCRABFWCKVDnTK/dt5ahC4aSfTibB3o9wJNXPknN6jUdq2/f6dOM2LyZf+Xl0b9hQ16NjaVJeLhj9YmISPDJz4fHH4eXXoKGDctHLt16KxqxJCISwI4dW01GxgCsLSMp6UPq17/C7ZBEJIBoipyLSstKeXrV01w26zIKigr4aPhHPPurZx1NLi08eJCuaWmsOnaMGZ068X6XLkouiUjI0BpM3rMW3nkH4uLgxRdhzBjYvBmGD1dySUQkkOXmziM9/UqqV48mOfkLJZdEpNKCKsEUSHYf282Vc65kwkcTGBg3kA1jN3Bl2ysdq6+gpITRmzdzw6ZNtKlVi/UpKfyueXOM7gZEJIQE2oMIf7NlC/zqVzB4MDRrVj4tbvp0aNDA7chERHwjFB9EWGvZuXMKWVlDqVu3F8nJX1CnTke3wxKRAKQEkwvmbZxH1+ldWbd/HbOvn807N75Dg9rO/XX+ZX4+3das4fUDB/h9q1Z8lpxMbJ06jtUnIiLB5dQpmDgREhMhLQ1efhn+/W+49FK3IxMR8a1gfxCRltad7Ow7OX16PwClpYVkZQ1j587JxMSMIClpOTVq6KmBiHhGazBVoaOFR7lr6V3M2ziPn7b8KXMHzqVddDvH6ispK+Op3bt5YudOWtSsySfdunF5/fqO1SciIsHngw9g/HjYsQOGDYNnn4UY5zY3FRERB504kc7Jk5nk5r5O48aDOXFiEwUFabRt+xStWj2i2Q0i4pWgSjBZaxcDi1NSUm53O5Zzrdy1klsX3sre/L080ecJJlw+gerVnPv6vz51iluzsvg8P59hTZvycseO1KseVL9uERFx0O7dcPfd8N57EB8PK1ZAnz5uRyUiIt6ytghrITd3NgDR0b8iJmaEkksi4jVNkXNYUWkRv//o9/SZ3Yca1WqwetRqJvae6FhyyVrL7P376bZmDVknT/JWfDxz4+OVXBIRITTX1qisoiL44x/Lk0rLl8PTT0N6upJLIiLBKi/vQzIzB7sdhogEAWUdHJR9KJuhC4aydv9aRncfzQtXv0BkeKRj9R0uLuaOLVt49+BB+tSvz5y4OFrWquVYfSIigcafR7r6g9RUuPNOyMqCAQPghRegdWu3oxIREScYUwNjqhMTM5LWrSe6HY6IBAElmBxgreUva//Cfcvuo3aN2iwYtICB8QMdrfPDI0e4bfNmDhYX88d27bi/ZUvCNMxVREQqIDcXHngA/vY3aNMGFi+Ga691OyoREfEVa8vYvv33Z14ZqlWrSUzMKFq3nkjNmlpYT0R8I6gSTP6wyPfBEwcZvWg0i7cspl+7fsweMJvmUc0dq6+wtJTf79jB8zk5xNepw5LERLpHRTlWn4iIBI/SUpgxAx59FE6ehMcegwkTQBuNiogEj9LSk2RlDefQob9TvXpDGje+iTZtJimxJCI+F1QJJrenPizdupRR74/iaOFRXrjqBcb3HE8147tlrp7ZvZseUVH0jY4GIKOggOsyMthRWMi4Sy7hj+3aUScszGf1iYhI8EpLg7FjYe1a+MUvYNo0iI11OyoREXf5wwNrXzp9+gAZGddx/Pga2rf/My1a3KPFvEXEMVrk2wdOFZ9i3NJxXDPvGppENCHt9jTuvuxunyaXAHpERTEoM5OPjhzhf3NySF6zhp2FhUxt25aXOnZUcklERH5UXl75Oks9e8K+ffD22/Cvfym5JCIC5Q+srbVj6tWr53YoXiso2Mi6dT05cWITXbq8R8uW9yq5JCKOCqoRTG5IP5DOLX+/haxDWdx72b089YunqFXdmYW1+0ZH80rHjvx640aKrSXcGN5NSOCGxo0dqU9ERIKHtTB3bvlaS4cPw913w5QpULeu25GJiIivHT78TzIzBxEWFkX37p8SFZXsdkgiEgKUYPJQmS3jT5/9iUc/fpTGEY1ZPmw5/dr3c7TOhQcPMnbLFqy1ADzUsqWSSyIi8qM2bSqfDvfpp3DZZbB8OXTr5nZUIiLihL17X2Hr1vFERnalS5fF1KrVwu2QRCREKMHkgT3H9jD8veGk7kzlhvgbmHntTBrWaehYfQUlJdyzbRuzDhygU61alAHjLrmE6fv2cWV09HdrMomIyMUF29oa5+reHdLTz/9egwbw17/CqFFQTRPkRUSCjrWlfP31A+TkvEDDhv2Jj59H9eqRboclIiEkqP7ENMb0N8bMPHbsmGN1vJPxDl1ndCVtbxqzrpvFuze962hy6cv8fLqvXctrBw5wS5MmHCkt5e+dO/NE27bMT0hgUGYmK/LyHKtfRCSYBNPaGufTqxeEh//weFwcZGfDb3+r5JKISDAqKSkgI2MgOTkv0KLFPXTpslDJJRGpckH1Z6aTNw75p/MZvnA4g/8+mNiGsaTfkc6o7qMcWyivpKyMP+zcyc/WraO4rIxPunUjKTKS+QkJ341Y6hsdzfyEBNKOH3ckBhERCSwTJ8K53VLNmrBiBTRq5E5MIiLirMLCHNLTL+fw4Q/o2HEaHTo8jzHa/EdEqp6myFXA6t2rGbZwGLuP7WZS70k8evmj1Air4Vh9O06dYlhWFp/l5zO0SROmdepEverVubx+/R+c21dT5EREBCgshFdfhZKS/xwLD4fRoyEmxr24RETEOcePr2fjxmspLT1OYuIHNGx4tdshiUgIU4LpIopLi3nikyd4atVTtKnfhlUjV9GrZS/H6rPWMjc3l3Fbt1INmBcfz5CmTR2rT0REgsPy5TBuHGzdCtddB8uWwenTEBZWPqpJRESCz6FDi8jMHEKNGo3o3n01kZGJbockIiEuqKbI+dLWw1v5+es/58lPn2R40nDSf5fuaHLpSHExgzMzuW3zZrpHRvJVjx5KLomIyEXt3Qs33wxXXVX+evlyeP/9/yzkPXKkRi+JiAQbay179jxPRsYAIiI6k5z8pZJLIuIXNILpHJNWTKJlvZbc8897CA8L5/9u+j9uTLjR0To/zstjeFYWucXFTG3blgdbtSLMobWdREQk8JWUwEsvweOPl//7H/4ADz5Yvt4SlI9a2rRJo5dERCrL33cbLSsrYdu28ezbN4NGjf6L+Pg5hIXVcTssERFACabvOXTyEE+sfAKAX7T9BbMHzKZF3RaO1Xe6rIzHduzgT3v20Kl2bd5PTOQnUVGO1SciIoFv9Wq4807YsAF+85vyRFO7dt8/p1kz+OQTd+ITEQlk1trFwOKUlJTb3Y7lXCUlx9i0aRB5ectp2fJh2rV7CmM0IUVE/IcSTGccP32cbjO6AfBcv+e4t9e9VHPwf9iZJ04wNCuL9IIC7mjenOfatyciTLs9iIjI+R06BA8/DK+9Bi1bwoIFMGDAD3eNExGR4HPq1E42bryWU6eyiY19lWbNRrsdkojIDyjlDUxOnUzdp+uy9/heAB741wOEPRHG5NTJPq/LWsvLOTn8ZO1a9p4+zaIuXZjeqZOSSyIicl5lZfDXv0JsLMyZAw89BJmZMHCgkksiIqEgP/9L1q3rSVHRXrp2Xabkkoj4raAaweTpnOnJfSYzuc/k8jKmGOwk6/vggAOnTzMqO5t/HDnCbxo04LW4OJqGhztSl4hIKDHGRAArgUnW2iVux+Mr69fD2LHw5ZdwxRXwyivQubPbUYmISFX55pt32bz5VsLDm5OY+AEREXFuhyQickFBNYLJWrvYWjumXr16bofyA4sPHSJxzRpWHD3Kyx07siQxUcklEQl5xpjXjDHfGGMyzjl+tTEm2xizzRjzSAWKehiY70yUVe/YMbj7bkhJgR07ykcupaYquSQiEiqsteza9TSZmTcRGZlMcvIXSi6JiN8LqhFMvjCp9ySflneitJT7t23jL/v30y0ykjfj40mIiPBpHSIiAWw28DIw59sDxpgwYBrQD8gB0owxi4AwYOo5nx8FdAUygVpVEK+jrIW334b77oPc3PLRS08+CdHRbkcmIiJVpaysiC1bxnLgwGs0aTKE2NjXCAsL+C5OREKAEkzn+HaqnC+sPX6cWzIz2XrqFA+2bMkf2ralZrWgGjQmIuIVa+1KY0ybcw5fCmyz1m4HMMa8DVxvrZ0KXHtuGcaYvkAEkACcMsYstdaWnee8McAYgKZNm5KamurD/xLv7d5dhxde6Mj69dHExuYzefJWYmOP89VXbkfmuYKCAr/7nsX/qd1IKCsuzmPTpv/i6NEVtG79OG3aTMZowT0RCRBKMDmg1Fqe2b2bx3fuJCY8nI+Skuirx88iIhV1CbDnrNc5QM8LnWytfRTAGDMCOHS+5NKZ82YCMwFSUlJsnz59fBSud06ehP/5H3j2WahTp3ydpTFj6hIW9hO3Q/Naamoq/vI9S+BQu5FQderU12zYcA2FhTuIi5tLTMwwt0MSEakUJZh8bFdhIbdmZfHpsWMMatyYGZ06EV2jhtthiYgEkvM9qv3R3RestbN/tGAPN4NwyuLFMH487NoFw4fDM89A06ZuRyUiIlXt6NFVZGQMACAp6UPq17/c5YhERCpP87V86M3cXLqmpZFeUMAbcXG8nZCg5JKISOXlAC3Pet0C2OeLgv1lM4hdu+D66+G66yAionwB7zfeUHJJRCQU5ea+yVdf/YIaNRqSnPyFkksiErCUYPKBo8XF3JKZybCsLBIjIvgqJYXhMTGaLy0i4pk0oKMxpq0xJhwYDCzyRcHGmP7GmJnHjh3zRXGVVlQETz8N8fHw4YflI5bS06F3b1fCERERF1lr2bFjMllZw6hX76ckJ39OnTr+McJWRMQTSjB5aeXRoyStWcP8b77hD23akNqtG21r13Y7LBGRgGCMeQv4HIg1xuQYY0Zba0uAccAyIAuYb63d5Iv63BzBtGIFJCXBhAnw619DVhY8+CBooKuISOgpLS0kK2sYu3ZNISZmBF27LqNGjQZuhyUi4hWtweShorIyJu3cyR9376Z97dqsTk6mZ926boclIhJQrLVDLnB8KbC0isNxxIED8MAD8Oab0K4dfPAB/OY3bkclIiJuKSo6SEbGQPLzV9O27VO0avWIZj6ISFBQgskDm0+cYFhWFmsLChgdE8MLHToQWV1fpYiIv6vKRb5LS2H6dHj0USgshMcfh0ceAQ1yFREJXSdObGbjxmsoKtpHQsJ8mjS5ye2QRER8RlPkKsFay4y9e0leu5YdhYUs6NyZV+PilFwSEQkQVTVF7t//hksvLd8hrmdPyMiAKVOUXBIRCWV5eR+zfn0vSksLSEpaoeSSiASdgEgwGWMijDFrjTHXuhXDN0VFXJ+RwditW/l5vXps7NGDgY0buxWOiIh4wOlFvo8cgTvugMsuK58aN38+LFsGHTs6Up2IiFSQ2/cT+/e/zoYNVxEe3pzk5C+pV+8yN8IQEXGUowkmY8xrxphvjDEZ5xy/2hiTbYzZZox5pAJFPQzMdybKH7f08GES09JYfuQIz7dvzz+7dqV5zZpuhSMiIh5yagSTtTB7NsTGwquvwj33lC/ifdNNoGU1REQ8F+j3E9aWsX37BLKzR1G/fl+Skz+jdu02VR2GiEiVcHpu12zgZWDOtweMMWHANKAfkAOkGWMWAWHA1HM+PwroCmQCtRyO9QdOlpby0NdfM23fPhIjIvgwKYnEyMiqDkNERPzYxo1w552wahX89Kfwyivlu8WJiIhPzCZA7ydKS0+xefNwDh58l2bNfkfHji9RrZq2DhWR4OVogslau9IY0+acw5cC26y12wGMMW8D11trpwI/GLJqjOkLRAAJwCljzFJrbdl5zhsDjAFo2rQpqampXsW+FfgfYBdwI3D7iRMcXrMG70oNLgUFBV5/zxJ61G4kWBQUlK+r9PzzUL8+zJoFI0ZAtYCYfC4iEhgC937iCPAokA2MZf/+m9i/f7UX5QUn/V0onlC78V9urE59CbDnrNc5QM8LnWytfRTAGDMCOHS+zuDMeTOBmQApKSm2T58+PROFUwAACOdJREFUFQ7omd276REVRd/oaMqs5bk9e/j99u3UCQtjeefO9GvQoMJlhZLU1FQq8z2LgNqNuMvTXeS6d4f09PO/d/vtMHUqNGzofXwiIlIhfnc/kZbWnbp1e9GmzUSKiw+zceNtFBcfIj5+AY0bD6hwOaFGfxeKJ9Ru/JcbCabzrUZhf+xD1trZvg+lXI+oKAZlZvJyhw7M2L+f1KNHCTeGOXFxSi6JiAQRa+1iYHFKSsrtlflcr16QmQlFRd8/fsMNMHOmDwMUEZGKcOx+wtMHESdOpHPyZCYHDszCWkuNGtF0776SqKifVKocEZFA5sZA/hyg5VmvWwD7fFGwp7sD9Y2OZnZsLEOyslh99CiRYWH8IzGRAdolTkREgIkTf7hYd+3aMG2aO/GIiIQ4x+4nvNkMwtoirC0CiikpOca+fbM4fXq/L8ISEQkIbiSY0oCOxpi2xphwYDCwyBcFe9MhXNOoEdc1bEgxcG+LFlypkUsiInJGs2YwahSEhZW/Dg+HkSMhJsbduEREQpRj9xO+Yu1p9u//C5mZg90ORUSkyjiaYDLGvAV8DsQaY3KMMaOttSXAOGAZkAXMt9Zu8lF9Ho1gAliRl8fq/Hwmtm7N9H37WJGX54uQRETEj3jTT0ycCDXObP4TFlb+WkREnFXV9xO+YEw41arVpnnzO0hIeMftcEREqozTu8gNucDxpcBSB+rzaG2NFXl5DMrMZH5CAn2jo+lbv/73XouISHDwtJ+A8lFMI0fCX/6i0UsiIlWlqu8nPF2Dqfyz4RgTRkzMSFq3nkjNmuooRCS0aDNlIO348e8lk/pGRzM/IYG048ddjkxERPzJxInw859r9JKISLDydMmNiIhuNGv2W3r23E6nTtOUXBKRkOTGLnKO8fSJw0OtWv3gWN/oaI1eEhGR72nWDD75xO0oRETE3/Tosd7tEEREXBdUI5i8WeRbRERERERCmzdr9YmIhLqgSjCJiIhcjG4cRETkYvTAWkTEc0GVYNKNg4iIXIxuHEREREREnBFUCSbdOIiIiIiIiIiIVL2gSjCJiIiIiIh4SjMiREQ8pwSTiIiIiIgImhEhIuKNoEow6YmDiIiIiIiIiEjVM9Zat2PwOWPMQWCXF0XUA3yRpfK0nMp8riLnenvOhd5rBBz6kXL9ha9+p1VRhyflVPYzvmg3nrQZULtxoo7KlNHaWtvYy/oCnvqJSp/j6fu63n1fh7+0mYqcF+z9RKC0mcqWo36CgO8n/O3vwou9HyjXOwTONa9+wn8ESpupbDkX7iestfo55weY6WY5lflcRc719pwLvQescft3VdW/06qow5NyKvsZX7QbT9rMmffUbnxcR1XEqR9nvvNA6Sc8fV/Xu+/r8Jc24027qMB7AdFuAqXNVFWs+vGf352//V14sfcD5Xr35e/U6TrUT/jPT6C0GV+WE1RT5HxoscvlVOZzFTnX23N89X24qSr+G9xsN5X9jC/aTbC3GQicdhMs33cgCbV+wtv3A0GwX+++bjMVOS/Y+4lAaTO+LEcqTn8XVr4Ofxco17z6Cf8RKG3GZ+UE5RQ5qRrGmDXW2hS345DAonYjEjp0vYsn1G5EQoeud/GE2o3/0ggm8cZMtwOQgKR2IxI6dL2LJ9RuREKHrnfxhNqNn9IIJhERERERERER8YpGMImIiIiIiIiIiFeUYBIREREREREREa8owSQiIiIiIiIiIl5RgklERERERERERLyiBJM4whgzwBjzV2PM+8aYX7kdjwQGY0w7Y8wsY8y7bsciIs5SPyGeUD8hEjrUT0hlqY9wnxJM8gPGmNeMMd8YYzLOOX61MSbbGLPNGPPIxcqw1r5nrb0dGAHc7GC44id81G62W2tHOxupiHhL/YR4Qv2ESOhQPyGVpT4iOBhrrdsxiJ8xxlwBFABzrLVdzhwLA7YA/YAcIA0YAoQBU88pYpS19pszn/sT8Ka1dl0VhS8u8XG7eddae2NVxS4ilaN+QjyhfkIkdKifkMpSHxEcqrsdgPgfa+1KY0ybcw5fCmyz1m4HMMa8DVxvrZ0KXHtuGcYYAzwN/EOdQWjwRbsRkcCgfkI8oX5CJHSon5DKUh8RHDRFTirqEmDPWa9zzhy7kPHAL4EbjTF3OBmY+LVKtRtjTENjzAyguzFmgtPBiYhPqZ8QT6ifEAkd6iekstRHBBiNYJKKMuc5dsH5ldbaF4EXnQtHAkRl281hQH9AiAQm9RPiCfUTIqFD/YRUlvqIAKMRTFJROUDLs163APa5FIsEDrUbkdCh6108oXYjEjp0vUtlqc0EGCWYpKLSgI7GmLbGmHBgMLDI5ZjE/6ndiIQOXe/iCbUbkdCh610qS20mwCjBJD9gjHkL+ByINcbkGGNGW2tLgHHAMiALmG+t3eRmnOJf1G5EQoeud/GE2o1I6ND1LpWlNhMcjLUXnMIoIiIiIiIiIiLyozSCSUREREREREREvKIEk4iIiIiIiIiIeEUJJhERERERERER8YoSTCIiIiIiIiIi4hUlmERERERERERExCtKMImIiIiIiIiIiFeUYBIREREREREREa8owSQiIiIiIiIiIl5RgknEIcaYXxpj5rodh4iI+Cf1EyIicjHqJyTQKMEk4pwkYL3bQYiIiN9SPyEiIhejfkICihJMIs5JAmKMMZ8aYw4YY37pdkAiIuJX1E+IiMjFqJ+QgKIEk4hzkoBD1trLgTuBoS7HIyIi/kX9hIiIXIz6CQkoSjCJOMAYUwNoADx35lB14Kh7EYmIiD9RPyEiIhejfkICkRJMIs5IAL6y1paded0VyHAxHhER8S/qJ0RE5GLUT0jAUYJJxBlJwFdnve4KbHApFhER8T/qJ0RE5GLUT0jAUYJJxBlJfL8D6IKeOIiIyH+onxARkYtRPyEBx1hr3Y5BREREREREREQCmEYwiYiIiIiIiIiIV5RgEhERERERERERryjBJCIiIiIiIiIiXlGCSUREREREREREvKIEk4iIiIiIiIiIeEUJJhERERERERER8YoSTCIiIiIiIiIi4pX/B6pkEr9Z4uvfAAAAAElFTkSuQmCC\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"figure(1,figsize=(20,5))\n",
"\n",
"subplot(1,3,1)\n",
"loglog(H,err_ep, 'r-o',label='AB-1=EE')\n",
"loglog(H,err_AB2, 'g-+',label='AB-2')\n",
"loglog(H,err_em, 'c-x',label='EM')\n",
"xlabel('$h$')\n",
"ylabel('$e$')\n",
"title(\"Schemas explicites\")\n",
"legend() \n",
"grid(True)\n",
"\n",
"subplot(1,3,2)\n",
"loglog(H,err_er, 'r-o',label='AM-0=EI')\n",
"loglog(H,err_CN, 'b-v',label='AM-1=CN')\n",
"xlabel('$h$')\n",
"ylabel('$e$')\n",
"title(\"Schemas implicites\")\n",
"legend() \n",
"grid(True)\n",
"\n",
"subplot(1,3,3)\n",
"loglog(H,err_heun, 'y->',label='Heun')\n",
"xlabel('$h$')\n",
"ylabel('$e$')\n",
"title(\"Schemas predicteur-correcteur\")\n",
"legend() \n",
"grid(True)\n",
"\n",
"show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"colab": {
"collapsed_sections": [],
"default_view": {},
"name": "EdoExplicites.ipynb",
"provenance": [],
"version": "0.3.2",
"views": {}
},
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.10"
},
"latex_envs": {
"LaTeX_envs_menu_present": true,
"autoclose": true,
"autocomplete": true,
"bibliofile": "biblio.bib",
"cite_by": "apalike",
"current_citInitial": 1,
"eqLabelWithNumbers": true,
"eqNumInitial": 1,
"hotkeys": {
"equation": "Ctrl-E",
"itemize": "Ctrl-I"
},
"labels_anchors": false,
"latex_user_defs": false,
"report_style_numbering": true,
"user_envs_cfg": true
},
"toc": {
"base_numbering": 1,
"nav_menu": {},
"number_sections": true,
"sideBar": true,
"skip_h1_title": true,
"title_cell": "Table of Contents",
"title_sidebar": "Contents",
"toc_cell": true,
"toc_position": {},
"toc_section_display": true,
"toc_window_display": false
}
},
"nbformat": 4,
"nbformat_minor": 4
}