{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from IPython.display import display, Latex\n", "from IPython.core.display import HTML\n", "css_file = './custom.css'\n", "HTML(open(css_file, \"r\").read()) " ] }, { "cell_type": "code", "execution_count": 2, "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_CM6 : Schémas à un pas de type Runge-Kutta " ] }, { "cell_type": "markdown", "metadata": { "toc": true }, "source": [ "

\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Schémas de Runge-Kutta\n", "\n", "Les schémas de Runge-Kutta sont des méthodes à un pas qui approchent l'intégrale $\\int_{t_n}^{t_{n+1}} \\varphi(t,y(t))\\mathrm{d}t$ par une formule de quadrature en des points donnés toujours à l'intérieur de l'intervalle $[t_{n};t_{n+1}]$:\n", "\\begin{align*}\n", "t_{n,i} & \\stackrel{\\text{déf}}{=}t_n+c_ih \\qquad\\text{avec }c_i\\in[0;1]\n", "\\\\\n", "\\int_{t_n}^{t_{n+1}} \\varphi(t,y(t))\\mathrm{d}t\n", "&\\approx\n", "h\\sum_{i=1}^sb_i \\varphi(t_{n,i},y(t_{n,i})).\n", "\\end{align*}\n", "Le problème est que, si $c_i\\neq0$ et $c_i\\neq1$, alors $t_{n,i}$ n'est pas un point de la discrétisation et on ne connait pas $y(t_{n,i})$. \n", "L'évaluation des $\\varphi(t_{n,i},y(t_{n,i}))$ en les points intérieurs est alors approchée par une autre formule de quadrature:\n", "$$\n", "\\varphi(t_{n,i},y(t_{n,i}))\n", "\\approx\n", "\\varphi\\left(t_{n,i},y_n+h\\sum_{j=1}^{s}a_{ij} \\varphi(t_{n,j},y(t_{n,j}))\\right).\n", "$$\n", "Cela donne\n", "$$\n", "y_{n+1}\n", "\\approx\n", "y_n\n", "+\n", "h\n", "\\sum_{i=1}^s\\left[b_i \\varphi\\left(t_{n,i},y_n+h\\sum_{j=1}^{s}a_{ij} \\varphi(t_{n,j},y_{n,j})\\right)\\right]\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Une méthode de Runge-Kutta à $s\\ge1$ étages s'écrit:\n", "$$\\begin{cases}\n", "u_0=y(t_0)=y_0,\\\\\n", "u_{n+1}=u_{n}+h \\displaystyle\\sum_{i=1}^sb_iK_i& n=0,1,\\dots N-1\\\\\n", "K_i=\\displaystyle\\varphi\\left( t_n+hc_i,u_n+h\\sum_{j=1}^{s}a_{ij}K_j \\right) & i=1,\\dots s.\n", "\\end{cases}$$\n", "\n", "### Matrice de Butcher\n", "Les coefficients sont généralement organisés en deux vecteurs $\\mathbf{b}=(b_1,\\dots,b_s)^T$, $\\mathbf{c}=(c_1,\\dots,c_s)^T$ et une matrice $\\mathbb{A}=(a_{ij})_{1\\le i,j\\le s}$. Le tableau\n", "$$\n", "\\begin{array}{c|c}\n", "\\mathbf{c} & \\mathbb{A}\\\\\n", "\\hline\n", "&\\mathbf{b}^T\n", "\\end{array}\n", "$$\n", "est appelée *matrice de Butcher* de la méthode Runge-Kutta considérée. \n", "\n", "\n", "+ Si $a_{ij}=0$ pour $j\\ge i$ (i.e. la matrice $\\mathbb{A}$ est triangulaire inférieure stricte) alors la méthode est *explicite* car chaque $K_i$ peut être explicitement calculé en fonction des $i-1$ coefficients $K_1,\\dots,K_{i-1}$ déjà connus;\n", "+ dans les autres cas la méthode est *implicite* et il faut résoudre un système non linéaire de dimension $s$ pour calculer les $K_i$. L’augmentation des calculs pour les schémas implicites rend leur utilisation coûteuse; un compromis acceptable sont les méthodes *semi-implicites*: \n", " + si $a_{ij}=0$ pour $j> i$ (i.e. la matrice $\\mathbb{A}$ est triangulaire inférieure) alors la méthode est dite *semi-implicite* ou *diagonalement implicite* (diagonal implicit RK, DIRK), i.e. chaque $K_i$ est solution de l’équation non linéaire\n", " $$\n", " \\boxed{K_i}=\\varphi\\left(t_n+c_ih, u_n+ha_{ii}\\boxed{K_i}+h\\sum_{j=1}^{i-1}a_{ij}K_j\\right)\n", " $$\n", " Un schéma semi-implicite implique donc la résolution de $s$ équations non linéaires indépendantes. \n", " Si de plus tous les termes diagonaux sont égaux $a_{ii}=\\gamma$ pour $i=1,\\dots,s$, on parle de *singly diagonal implicit method* (SDIRK). *Remarque*: la méthode est implicite non pas parce que $u_{n+1}$ dépend de lui même, mais parce qu'au moins un $K_i$ dépend de lui même." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Ordre de convergence\n", "\n", "\n", "Soit $\\omega$ l'ordre de la méthode.\n", "
\n", "Proposition (consistance).\n", " \n", "La méthode de Runge-Kutta est consistante (i.e. $\\omega\\ge1$) ssi \n", "$\\begin{cases}\n", "\\displaystyle\\sum_{j=1}^s b_{i}=1&\n", "\\\\\n", "\\displaystyle c_i=\\sum_{j=1}^s a_{ij}&i=1,\\dots,s\n", "\\end{cases}$\n", "
\n", "\n", "**Remarque** En particulier, pour une méthode *explicite* on aura $c_1=a_{1,j}=0$ ainsi $K_1=\\varphi(t_n,u_n)$ et $c_2=a_{21}$ ainsi $K_2=\\varphi(t_n+c_2h,u_n+hc_2K_1)$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "Proposition (ordre).\n", " \n", "
\n", "\n", "(For ERK cf. page 135 de E. Hairer, S. R N0rsett, G. Wanner, *Solving Ordinary Differential Equations I. Nonstiff Problems*, Second Revised Edition, 2008)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Théorème (Barrières de Butcher)** \n", "Soit une méthode RK à $s$ étape d'ordre $\\omega$.\n", "- si la méthode est **explicite** alors $\\omega\\le s$\n", "- si la méthode est **implicite** alors $\\omega\\le 2s$\n", "- si $s\\ge5$ alors $\\omega\n", "\n", "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " Théorème de Kuntzmann et Butcher pour des méthodes d'ordre 2s\n", "\n", "Soit les conditions suivantes:\n", "+ $B_p$: $\\sum_i b_i c_i^{q-1}=\\frac{1}{q}$ pour $q=1,\\dots,p$\n", "+ $C_\\eta$: $\\sum_j a_{ij} c_j^{q-1}=\\frac{c_i^q}{q}$ pour $i=1,\\dots,s$ et $q=1,\\dots,\\eta$\n", "+ $D_\\zeta$: $\\sum_i b_i c_i^{q-1}a_{ij} =\\frac{b_j}{q}(1-c_j^q)$ pour $j=1,\\dots,s$ et $q=1,\\dots,\\zeta$\n", "\n", "Si les coefficients de la méthode de RK vérifient les conditions $B_p$, $C_\\eta$ et $D_\\zeta$ pour $p\\le\\eta + \\zeta + 1$ et $p \\le 2\\eta + 2$ alors la méthode est d’ordre $p$.\n", "\n", "Cf. E. Hairer, S.P. Norsett, G. Wanner, *Solving Ordinary Differential Equations I*, page 208\n", " \n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Remarque (le cas des systèmes)** Une méthode de Runge-Kutta peut être aisément étendue aux systèmes d’équations différentielles ordinaires. Néanmoins, l’ordre d’une méthode RK dans le cas scalaire ne coı̈ncide pas nécessairement avec celui de la méthode analogue dans le cas vectoriel. En particulier, pour $p \\ge 4$, une méthode d’ordre $p$ dans le cas du système autonome $\\mathbf{y}'(t)=\\boldsymbol{\\varphi}(t,\\mathbf{y})$, avec \n", "$\\boldsymbol{\\varphi}\\colon \\mathbb{R}^m\\to\\mathbb{R}^n$ est encore d’ordre $p$ quand on l’applique à l’équation scalaire autonome $\\mathbf{y}'(t)=\\boldsymbol{\\varphi}(\\mathbf{y})$, mais la réciproque n’est pas vraie." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### A-stabilité\n", "Rappelons la définition de schéma A-stable:\n", ">Soit $\\lambda\\in\\mathbb{C}$ avec $\\Re(\\lambda)=-\\beta$ où $\\beta$ est un nombre réel positif et considérons le problème de Cauchy\n", "$$\\begin{cases}\n", "y'(t)=\\lambda y(t), &\\text{pour }t>0,\\\\\n", "y(0)=y_0\n", "\\end{cases}$$\n", "où $y_0\\neq0$ est une valeur donnée. Sa solution est $y(t)=y_0e^{\\lambda t}$ donc $$\\lim_{t\\to+\\infty}y(t)=0.$$\n", ">Soit $h>0$ un pas de temps donné, $t_n=nh$ pour $n\\in\\mathbb{N}$ et notons $u_n\\approx y(t_n)$ une approximation de la solution $y$ au temps $t_n$.\n", ">Si, sous d'éventuelles conditions sur $h$, on a\n", "$$\n", "\\lim_{n\\to+\\infty} u_n =0,\n", "$$\n", "alors on dit que **le schéma est A-stable.**\n", "\n", "Dans le cas $\\lambda\\in\\mathbb{R}^-$, i.e. $-\\lambda=\\beta>0$, une méthode de Runge-Kutta à $s\\ge1$ étages pour $y'(t)=-\\beta y(t)$ s'écrit:\n", "$$\\begin{cases}\n", "u_0=y_0,\\\\\n", "u_{n+1}=u_{n}+h \\displaystyle\\sum_{i=1}^sb_iK_i& n=0,1,\\dots N-1\\\\\n", "K_i=\\displaystyle -\\beta \\left( u_n+h\\sum_{j=1}^{s}a_{ij}K_j \\right) & i=1,\\dots s.\n", "\\end{cases}$$\n", "\n", "On peut vérifier que, contrairement aux méthodes multi-pas, la taille des régions de stabilité absolue des méthodes RK augmente quand l’ordre augmente." ] }, { "attachments": { "A_stab_RK.png": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmoAAAJdCAIAAAANz8w9AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAgAElEQVR4nOzdeVxM+/8H8I9vIzTpqpRIG4m0EUmrqFBZK0nZ4yKkS103S2TrWq+Q9XZtyRZRkmWim63khiZEK2KqUZHOSJ3u/f1x7rdf39apzpkzM+f9fPRHpjPnvM+Y5tXncz6fz+n0zz//IAAAAAC0xX/oLgAAAACQPBCfAAAAQJtBfAIAAABtBvEJAAAAtBnEJwAAANBmEJ8AAABAm0F8AgAAAG0G8QkAAAC0GcQnAAAA0GYQnwAAAECbQXwCAAAAbQbxCQAAALQZxCcAAADQZhCfAAAAQJtBfAIAAABtBvEJAAAAtBnEJwAAANBmEJ8AAABAm0F8AgAAAG0G8QkAAAC0GYvuAkTh48eP6urqdFcBAABA1GpqalgsSpKu0z///EPFfsVKp06d6C4BAAAAPSiKOUa0PhFlL59EIP56YOwrwPDTRwh16sSIv5KbA6ePGPz+p7TtBNc+AQAAgDaD+AQAAADaDOITAAAAaDOmXPtkMsZe9iAw/PQBw98ADD99SkHrEwAAAGgzaWh9Yhh29erVnJwce3v7IUOGsNlsuisCAAAg5SR7SHdBQcHChQs5HE79B+/fv29lZVX/ERi5zuTTR4x/BeD04fTproI2lM7bkezO28mTJ3M4HDU1tf3798fHxzs4OCCErK2tuVwu3aWJESb/8hDgFWAyhv/vM/z0KSXBf5g8ePDA2toaIVRZWVnXYevo6MjhcBwcHG7fvl23JcP//gIMB+9/wFjQ+myar68vQigkJKT+xc6NGzcihBp05wIAAADkkuChQ8+fP+fz+YqKivUfTExMRAipqanRVBQAAABGkOD4RAipqKjU/yeXy92wYQNCaMmSJTRVBAAAgBEkuPO2Pi8vr969exsbGyOE/Pz81qxZQ3dFYgRuOAOvAJMx/H+f4adPKWmITxzHz549W1RURPxzwIABjbfp9L9QU+8qDMMaP9j4ERzHRb9ZRwqurKwUZjMh9ybOr1KTm2EYhhqh+vQl8VUS/emL5lVqMGyEaa8SiacvKa9Sg4966kjDkDwMw/h8/tevX2NiYojOWzU1tZycnLohRY1fRCk4awCEBCNvAaOI7ANf2n6vCgoKdHR0EELx8fHOzs7Eg/DxAZgM3v+AsWDiShtoa2sTV0AjIyPprgUAAIDUktT4xDAsKipq06ZNja9subm5IYRevHhBR10AAAAYQYJ7dUxMTDIyMkJCQoKDg+sexHG8c+fOCDpvgRjDMIzH48XHxyOE0tLS8vLymtvS0dFRSUnJwMBAV1dXW1u7fYeD9z9gLEo7byX496pu0b6SkhJiAiiO44cPH16+fLmamtr79+9ZrH9ntTL844Php4/E4BXAMOzRo0fnzp17/vz5i5cva3G8urqa+NGAgfoGRibNPfFK9Lm677t07aqspDx8+DBPT08HB4cGk55bQPvp0wtOn+GnjyA+m0Q0QBFCfn5+ysrKhw4dIqavNGiSMvwNBOjy9OnT8PDwhBs3ysvLvwkECKEp7p7DR1poaffTGzS4t3rfur/wWlbK55eUFL16wX2a9jj5Lif79StZWVm2vLzr1KlLly4dOnRoy0+H9z9gLIjPZvH5/EOHDtWlJkLIz8/Pw8MDblgGaESk5uWYmPKyMoSQz+LlQ81GmJlbamhpk7J/AYZlZjy7n5QYcfhAeVmpopKS69Spv/zyi66ubpPbw/sfMBbEZ+u4XO779+9HjRrV5L2y4eMDiACO4+Hh4Tt27vz44QNCyGfx8omu00zNzIVsYrbP40cP4i5fjDi8n9W58wgzs5CQEOK2ffXB+x8wFsRnR8HHB6AUn88PCgqKjIz8/v37FHfPuT/6jrCwav1p5MFxPPbyhbAd27Jfv+qjrn7yxIn6IQrvf8BYEJ8dBR8fgCI5OTmzZ89Oe/IEr6kJWLNhzoIlykKP6KHC40cPgn5a+uoFd9CgQXFxcUR3Lrz/AWNBfHYUwz8+GH76iJpXICcnZ86cOX+lp8vJsTfv3DveZbJcUxcOaJF487r/4vlfvnx2dnKKioqSl5dn8huA4e9/OH0E8dkRDH8DAXJhGLZixYqTp0798EOPkO17Jrl6UHp1s32I7txlPrO6Kyh8raiA9z9gJli0DwCxgON4WFiYoqJiRETExtBdT7MLXT28xDA7EUIsFsvVw4ubVzTE1AwhNG3aNBzH6S4KAKnCiGYZtD5Bx+Xk5FhZW5cUF09x99x14Kj4dNW2DMdxTcUuCKH+/fs/evRI+MUWAJAC0PoEgE44jk+ePNnAwOCHHkqJKc8OHj8jKdmJECIax4kpz8rKP2vr6OTk5NBdEQBSAuITgJY8ffpUSVk5NjY2eNvOxJRn+gZGdFfUHvoGRslPXqj31TQ0NLx27Rrd5QAgDSA+AWjWmjVrRo4cqdZbnZtXNH/RMvG8zCkkZRWVxJRnIyysp7q6QoIC0HEQn9Kv8b3XmaYdrwCfz9fX1w8NDZ01f1FiyjN6Z3OShcVinYm5bmFly6gEZfj7n+GnTymIT+kHw6ba+gpwOBxtHR3+p9Irt5I379wr0Y3OBhiYoAx//zP89CkF8QnA/wgJCXFxcRlmNjLpMVfEa++JRl2CTpkyFUYSAdBujJjRARNXgDBwHHd0dExKSvJZvHxD6C6paXT26S7z8WttgwdxHLcfOeRD4buC/HyYzQKkFaw61FEQn6BVfD7fwsIiNzf3dHSc/ThnusshU5PxiRAq5fOtTfV7KitlZWVJzd8KANQH8z4BoBaPx9PW0fn8pSIx5ZmUZWcLlFVUTl64mpubO2PGDLprAUDyQHxKPxh61/IrwOFwdPr169ZNLukxV0KndbbbCAurLbvC4q5dS01NpbsWqjD8/c/w06cUI3o1ofMWNOfatWtTXV0trGyPn4uRoLWE2qS5zts6o4Ybfix89/nzZ+jCBVIGOm8BoERoaKibm5uFle2ZmOvSmp3CiLx0DcMwNzc3ugsBQJJAfAKGWrBgwZo1a/wCgs7EXGd4q0tDS3vLrrAbN25IcRcuAKRjRK8mdN6CBhYsWBARERGwZsPKoGC6a6Fcq5236L/zWL5+KefxeKKpCgARgM5bAEiD4zijslNILBZr5/4jRUVFHA6H7loAkAyMaJZB6xMQcBy3tLRMS0s7EHHa1cOL7nJERJjWJ8F+5JAvn8s+FBZSXRIAogHLJnQUw+OT4aeP/vsK4Dg+YsSIp0+fXrmVLORqfL7zvJt8fPhIix6KSvoGRg0mujx+9ODE0YMIoYPHzzR+Vimfv/5nf+L7Jm+4jeN4SFBA6Sf+5h17SVykXvj4fPWCaz9yyO3btx0cHMg6Ou0Y/v6H00cQnx3B8DcQQAjhOD5o0KDc3FzhsxMh1Ke7TMsb2NjZ1x95dPlC1DKfWQihxnFVyuePGWnCLylWUe117c5DDS3txnvbE7pp17YQhFBqZm6TG7SP8PGJEBo13FCm0z+vXr0i6+gA0AiufQLQIe3LzjoBazakZubWfV25leyzeLmKai+E0L2kRL+Fc1rdQ/3svJPyvMlofPWCS2QnvYK37sjLy4Ol5AFoFcQnkHIdzE6EkHZ/XQ0t7bqvERZWm3fufZ77cYq7J0LoSvQ5AYa18PQG2dlkr6wAwzwnjWtHbaQbZT+WLd/9119/pbsQAMQdxCeQZh3PzhYsD/iF+CYz41lz2wiTnQihgGU/8kuKbezsya2wHVgsls/iZZFnzuA4TnctAIg1iE8gtSjNzvq6Kyg0+biQ2Zl48/qV6HP6BkYbf91NXZHCm+Y1+3tV1Y0bN+guBACxBvEp/Zi5ZjTV2Ynj+P5dvyKEGo+/JQiZnaV8/iz3iQihE+evyMt3J73OdtDQ0h4wUD8wMJDuQsjBzPd/HYafPqUYvVYZQzBw1DG52clJiG/wSEFuTvzVy69ecBFC3vMWNH5KXXYihCa5ejSXnTiO+873RggdiDitoaX9/m1BB0sly5yFi7dtWIPjuBQsZ8jA9399DD99Skn87wYADZDe7rwSfe5K9LnGj9vY2a9as6HJQ9RlJ0Io4vD+ia7TmtzsVMThe0mJNnb24raGg+P4CesCVnC53KFDh9JdCwBiCjpvgVShos9W38Boirsn8VX34IGI0+fjbjV3CGIcUE5RBTG/ZeHMaaV8foNtXr3grgtYoaLa6+AfTayxQC9ias3FixfpLgQA8QXxCaQHRdc7l678+eDxM8RXTlEFMT52mc+sx48eNPcUYjkFOTb72p2HCCF+SbHvfO/6Y1lxHCdmquw5+DuJCwyRaIq7Z0xMDN1VACC+ID6BlBDNOFs5NvtMzHUiQaeMtW3uamXdUkTEvcAQQveSEvft3Fa3Ae9DIdG7O8t9Yp/uMsSXuWF/4qfmhv37dJe5fCGKorMQhoOTS35BAY0FACDmID6BNODz+RoaGiKYo4IQYrFYdd2tE8ZYNjk/sv6Im9k+i4m43bUthBhtJBH0DYy+V1UVQIIC0AyIT+kn9SPXeTzeYAODGryWm1dEdXYSlFVUDkScRgjxS4pDggJa3rh+3HpOGkcsUdRbvW/9hQCJr9PRccRmp6PjUjNzx7tMpvIkWqGl3Q8h9OjRIxprIIXUv/9bxvDTpxTEp/ST7pHrPB5Pd8CAf/5BSY+5oryI6OrhRbQpIw7vb7VNqayiQkQjv6R4nudUhBCLxaq/ECDxpTdoMLG93qDBGlrajW/JIkrE0blciWkuN0e63/+tYvjpUwriE0gwIju7dZMTcXYSdh04SnyzzGdWq0vc2Y9zJgbu3ktK/OPIAcqLI8MUd8+kpCS6qwBATEF8AknF4XBozE5Ub1jQqxfcUxGHW91+37GTxDyWdQErJOIiqHJPldy8PLqrAEBMMeJGmHC/T+lz7dq1qa6uFla29W+3CZrUpvt91nf5QlTg8kUt308GAHFG6f0+4XMHSJ6QkJBt27aRmJ3v3xbw+SV/PX70/K8nbwvy8Jqagvy86urvTW6s1ltdSVkZIWRqZt5XU8thnIuCwg/iOXcTAEAdRjTLGN76lKbTx3F88eLFERERAWs2+AWuaXd2vn9bcCP+auylCzlvXldXV38T/NvA6qGobGHngBCyHevS3HO5fz0uL+UjhB4lcT6XlyKEOsvKslgs1V5qltajJrt7GhiZiFWatrv1+eoF137kkJKSEhVxOp22kqb3fzvA6SPKWp+MeGUZ/gaSGjiOGxkZZWVlbdkVNn/RsrY+vS4yX2ZmfBMIEEIjbceMHj9JQVFpgL4hW767uqZ2O6r68K6ghPfhY+G75FvxRKB2lpWVk2Pbj3Xy8J5jaWtHe99yu+Pz/dsCc8P++fn52traZBcFgChAfHYUxKcU4PF4BoaG5WVlbV0YIfP50+NHwq/HXfnyuRz9NzKHWdr209OnItjKPvFfPHvy4tmTqGPhn8tLu8mxB+gNXLVm42jHcXTlKMQnYCyIz46C+JR0165dm+7p2aZBtvm5OQd2/0qkZg9FZa+FS81txhgNMxdlhn14V5B8K/78icN5b7Jku3SdMNl1sd9KQxNR38ME4hMwFsRnR0F8Si7iYmfkmTMjRlodPxcjzEoCyXc4yxfO+fy5rKa62mvBMteZ8/UGN3E7a1H68K4g7sJpoj2q1rvP+i07JrpOE1mQQ3wCxqI0PmHeJxBfPB7PyMgoIiJi+crV5+NutZydAgw7djBMt5eC5+RxPXqq7D0Rnf7x2y/bfqM9OxFC6praiwPWJ78uOhmX1F1ReanPzMGaKscOhrW62II4EAgEdJcAgDiC+ARi6vDhwzo6OvxPpYkpz1YGBbewpQDDVvkuGKylsmH1ylHjJkQnpV9KSrdxcKJ9zE5jQ82tLiWlJzzJtnFw2rB6pZGOmviHqKysLN0lACCOID6ln8StGY1hmJWV1ZIlS5wmTk3NzNU3aLb5iOP4Kt8FgzSUz54+vmDF6pT8z9uPRIpDc7Nl6pra249EJjzJthw9dsPqlYPUleKvXqa7qGaJ4V8hbSJx739yMfz0KQXxKf0k67rv2bNne/fu/Srr9enouIPHzzTXYYvj+LGDYUY6amdPH//xp6CU/M+LA9bTu8Z6W9WFqPFw84Uzp9kOG5yfm0N3UVJIst7/pGP46VMK4hOICxzHJ0+ePG/+/CHDRiQ/eWE/zrm5LdOfPB6krrRh9UrL0WOTXn6UuOCsT11T+2j0zZNxSaWlZaPNjbdtXCPmfbkAAALEJxALly9fVlJWjo2N3R1+7HzcreZmpwgwbLKDjev40era/RKeZG8/EqnUU4IXxKkz1NwqkfvOfdbCA7u3mw7UgGYoAOIP4hPQjMfjGRsbu7m5DTE1S83MdfXwam7LYwfDBmkop6U+DNl79FJSevsWCRJbLBbrl22/JTzJrq39Z7S58bGDYXRXhNJSH3bp2hVmrQDQJIhPQBscx9esWdNfV/cjr+h0dNz5uFsaWtpNblnK59sOG7xh9UrHCa4p+Z9d3GaItlLRUdfUTuS+s3eesmH1yon2VrR35P7nP/ARAUDT4HdD+onn0DsOh6OkrBwaGrrUP/BpdmELVzpjLpwdaTzgE/9TeFTs9iORknuZU0gsFmv7kcjwqNi/HqcM1uxZzOPRVcnTtMfdu3en6+hkEc/3v8gw/PQpBfEp/cRt6F1OTo6VlZWjoyPRW7syKLi5qRE4js/1mPKT7wIjU7Mb6bk2Dk4iLpVGNg5OSS8/dpbtajlEL/P5U1pqKP3E79+vHy2HJpG4vf9FjOGnTymITyA6GIYtWLBAf/Dg12+yW+6tRQgV83imAzVuJcSF7D16NPqm1Dc6G1PqqRJz77m+sanzaIvbCfGiL+DPO7f7SX58AkARiE8gCjiOh4WFKSoqRkREbAzd1XJvLUIo+Q7H2lS/tvafhCfZUnyls1VKPVV+v3x7uIWNj7fb/t2/ivjo5WWlZmZmIj4oAJIC4hNQLiwsTElZ2d/f32WyW05RxfxFy1peyGZ3aMhM9wlGpmYJT7KlbHhtO7BYrKPRN3/8KWjn1o1b1q0W2XHfvy1ACNna2orsiABIFslejguIOQ6HM2fu3I8fPkxx9wzauLWFrloCjuPeU53vJSV6LVgWsGmnpC8XR6LFAesRQkd/C0UIrduyXQRH5H38ICsr27dvXxEcCwBJBB9PgBIcDmfWrFlFRUX6BkanL8W3sG5tHRzHbUz13+bnhR46xeQO2+bUJaiB8ZCpHpS/PpkZT2VYLBXh7q4KAANBfEo/Ed/uNCcnZ86cOX+lp2tq6Rw+dWGEhZUwzxJgmL3FkLf5eSfjkoaaC/UUBlocsD4/+/WKxfPkuys4OrlQeqwnKY8G6+tTegjRYPjtfhl++pSS+GufOI4/ePBgxYoVvXv33rRpE5fLpbsisSOyXx5iRoqhkdHrN9m7w4/9+SRTyOws5vGGDdL68vkLZGertoafGG5hs3CWB9UL+12JPmdjY0PpIUSD4eHB8NOnlGT/YRIVFeXt7d3gQTU1tVu3bhkZ/X9vIfz9RbWcnJzZs2enPXnyww89QrbvmeTqIfxly2Iez3KIXtdu7Jh7z6VjAVuq4Thub6RZ9Q17nJnX3OLA9fXpLvPxa22bDlHK5xv1U7t9+7aDg0N7ywSAfsSqERR9/ktw6xPHcSI7/fz87t+/X1lZGR8fb2xsXFRUtHLlSrqrYwoej2dpaak/eHDW6zd7D//xNLvQ1cMLspNSLBYr5t7zLl26TR0/iqJV/Z6lp7E6dzYxMaFi5wBIBwmOzwsXLiCEjI2Nw8LCrKys2Gy2s7Mzh8NBCHE4nIKCAprrk3bEGgj9dXWJ4Hye+7FNwYkgOztAqadK2KlLOW9erw9cQcX+n6c/6d69O4wbAqAFEhyfq1atQggtXLiw/oMqKipEd9OpU6foKYsBiOBUUVWNiIhY6h/YjuBEkJ0dNtTcyvfn4KhTf1CxINHVSxdsrK1J3y0A0kSCR95mZGRgGNZ4XlpmZiZCSElJiY6ixBGJl34xDFuxYkXU2bPfBIKANRsW+61q31p6kJ2kWOAflJ5yf/FcT25eEYmLGgowLPv1q5AN68naIb0YPvSB4adPKQlufaqoqGhrazdo9Dx48KCoqAghNHToUJrqEjs1NTUd30mDFmdOUcXKoOD2fWSX8vl2I4wgOzuOxWL9ejjym0AwfdJYEnebmfFMVlZWagYNkfL+l1wMP31KSXB8Nobj+MaNGxFCxsbG5ubmdJcjLjq4dg8RnMRytR0MToQQjuMT7C07dfoPZCcplHqqhB46xX3+LP7qZbL2eT8pkS0vLzUXPhm+dhXDT59S0hOfOI47OTkRQ4euXr3a4E3T6X+hpm6Dh2FY4wcbP4LjuOg360jB7d6MqI0IzhWBazoYnKjeukK/x9wWn+yMv3TWWLVz4y83O9PVi2ae/T38zcuGk4mJDeIvnW28NxzHf3QfR2xw9vfwuscFGHaPk/Cj+7jRhn1XL5p5j5NA1qBZF7cZpuaWPy3xIWuHdRc+Rf+Wa3IzWn7jmtyMltOHV6mtmzX4qKeOlMQnn8+vy8779+9ra2s32OCf/4WamgnEZrMbP9j4ERaLJfrNOlJwOzYjbpCiqtoLITTbZ3HHg5PgPMqcWFdIb3Dra/jRLvslNyHmfOgaf3c708O7NgvzFBzHfT1dUpLvIIR8fw6esWAp8fjT1AcjdXos9ZqUknyntKQ4Ieb8Uq9Jvp4uZAXer4cjK79WLPBy7/iuSvn87NevPD09kWjfci1sRstvXJOb0XL68Cq1dbMGH/XUkYb45PP5xsbGddlpZQXL1rQfEZw9evTw9/cfZT82NTN38869pAxLWeW7IDPjWXhUrNiuK5TwJLv+V3RSuu/PwQMGGyGEDu7Y9DT1QctPb5CdxBK1xOMrfaYjhLwWLEt4kp3+8Vt4VKyyaq+U5Du7ggNJqVypp0rQtr13E292fCmiZ+lp3eTkJk2aREphAEgxiY9PIjuLiorU1NRKSkogOzuCw+Go9+3r7+/v6DQxNTP34PEzrd4jRUgxF87GRJ/z/TnYxsGJlB1SQV1Tu/6X3mCjxQHrT8ffI35682p0C8+tn52hh07VZSdC6ObVi6UlxcqqvQI27VTX1GaxWDYOTtPnLkII3Yy9SFbx0+Yu0tDuN2VsRxdSuHTuzA8KCmzm3ZwcgLaS7Ph88OCBqqpqUVGRg4NDRkaG1Ax2IJcwFwBSU1N79+7t6Oio3FM1MeUZicGJEEq+w1mxeN7o8ZPqh4qkkGOzvRYsQy1GXf3sPBmX1OB2Me/zc0bajvnRP6j+9fgxzlMQQqUlxWWf+KTUyWKxdhyN4pcUHT8a3vrWzcBx/Er0uXnz5pFSkpig+gKYmGP46VNKggdl8fl8a2trhJCDg0NCQgIMMGtOyxcA2neDFOEV83g+3u7DLWy2hp8gd88ik/fmJUJohJVdkz9tkJ2Nu6ab/KMhImw7QkhZtReJQ6j0Bhs5TZ2+PWT9vB+Xtu/XIfv1K4TQ/PnzySpJHDB81iPDT59SEtz6XLHi3+XK9uzZU1hYWPC/+Hxy/qiXYsSMFP3Bg4kbpCSmPCM9O3Ecd7Q2le3Sde/JyxL69038pbNENHrOX9L4p61mZ2Mf3hWc/T08IeY8QuhH/yByq/Vbu0UgwNrdAE2IjemuoKCrq0tuVQBIJYn8REMIcbncs2f/nTZgbGzceIMZM2ZERUWJtiiJgeN4eHh4YGBgTU3Nll1hs30WU5RtbuNHfyopSXiSTeKaONSpPxGlorzsWdqjnNcvs19yEULKqr0aR2NFeVlddiKE2N0VWj2EsWrnuu/Do2JJvxKsrqndkQZoxOEDo+3syC0JAGklwfFJdwmSisPhTJkyBcOwKe6euw4cpS7Ytm1ck5b6MDwqVl1Tm6JDkCtoyezGDw4YbOS/bquFnWPjH4Wu8Se+UVbtVVpSvMjDKT71dauvJ7ExQijYf+GRCwmkz+HxWbE6Ieb8w+Qk2zFtWzbo1QtueVnp8uXLya0HAGnFiOUQYdVHAnGZ8+HDh/oGRifOXyFxcFBjtxPiF8x0X+j/i/gPF4q/dJYIztBD/95mgPvX46jfDyCElFV7NZmI9RuRJ+OS2N0V3O1MEUIjbcccjb7Z6hHfvORGhG0n+m+jk9JJT1A3O9PPn/jPcz8Q/xTyfp97QjcdOfBbxZcv5BYDAI3gfp+gQzp16oTj+Jo1a4xNTF6/yT4dHZeY8ozS7Czm8Xznew8bab2A7Gt7lHJxm0F8/bLtt6SXH4lmoov5wOYGxyqr9opOSh9qbqU32Mj352CEUErynforDTVHb7BR3UCqO9evkHcG/1q3fT+/pCj5Dkf4p+A4HnH4gMe0aaQXQzuGDz1l+OlTCuJT+iUkJCgpK4eGhi71D3yaXWg/zpnSw+E4PmWcrUQPF0IIKfVUibx+HyFUWlL8y+KZTU6mDAjZWddwXOAfRCywELrG/8O7glb3z2KxiMQ9f+IIiWUThppbKSqr7Ny6QfinpKellpeVLl26lPRiaMfwNdMZfvqUgviUZhiGWVlZOTk5DTE1S83MXRkULII8W+23+G1+3rFLNyViuFAL1DW1ie5cYZYHYrFYx/7bbTvT2boubnEcX71oprFq58aLFn0uK0MIDRhkQHLdCCGEFq9a+4KbIcAwIbePu3yxu4KCVN6nSHL/hiMFw0+fUhCfUissLEy5Z89XWa8PRJw+H3eL0t7aOjEXzl66cDZo216JWNW2VS5uM0bajkEIRf1+oPHC8Q0QNz9BCJWWFK9dOpd4kMVi6QwYiBDasvp/huQIMIy4vLokMJiCwtHkGXOrvgn27twqzMY4jkcc3i+VPbcAUAfiUwrxeDx9fX1/f3+nCVPSXua7eniJ5rjFPN5PvgtMzS3rlkqXAhv2/Nu5GuQ7p9X18OriNiHmfN00mNlLViKEsl9yz/4eTuzhzUtuyMpFqJn5MKSQY7Odpk6PjjotzMbpaakIIansuQWAOhCf0ubs2bM6/foVlwyAafQAACAASURBVPBPR8cdPH5GZD2oOI6PtR7Olu++9yRpN54UB+qa2kHb9iKEsl9yLwpxnXLvycvKqr0QQkFLZhMXQYkkQwiFrvE37dPNzc7U3c6UGHa7ae8x6ir3nL+ktPRTMY/X6pZS3HMLAHUgPqUHcaXTy8uLaHRSPUSogfWBK/glRVJwybOxaXMX1Q0LanWJWjk2e0/EeeJ7v9muxDfbj0QS929RVu2V/ZI70nZMeFRs+sdvlC6gbzTMvKa6esfmViYOET23PtK1UB8AIsCICZFMmPeZmprqOHbs14qKAxGnRdZbWyf5Dme2x+RVG36Vpm5bipR94ovsVuG/rvkp9vzJyq9fW5j3+fjRgyljbT9+/Ni7d2/RVCViTPj1bwGcPoJ5n6A5OI4vXbrU1tZ2iKkZN69I9NkpwLB5M6aamltOm7tIxIeWRCLLToTQuMnu379/b3kboudWWrMTMX7NdIafPqVgTLNk4/F4NjY2ubm5lC5d27IJYyy7dOn26+FIGCIvbgYaDqmprm5hg397bn18RFYSAFIDWp8SjBgl9PlLRWLKs/mLltGSXscOhmW9zNwaflyUjSogJDk2mxgJ3BwYcwtAu0FzQSLhOD5jxozo6GgbO/vj52LoGq2Tn5uzNTjIa8EySofACIkY5pr2IKnyawXv3duM9Md1P6qpru4sK0t831dLx2DocISQ8bCRyiqqisoq0jfWqb7JM+amJN/BcbzJv65gzC0A7caIq8pSdvGcz+dbWFjk5ubSMkqoDo7jpgM1amv/SeS+E3HDF8fx4o+Fd2/EZr/k5r3Jys56Wft37XfB/6ywY+ni0eRz32ZlfMjNqvunbFe5Tp1Qfz39HopK1g7jjYeN1B1kIE2B+uFdgdPwAbfuPzE0aSIj+3SX8fHx+f3330VfGAAiQOnQIWh9SpjU1NQx9vbdusklpjzTN6BzZZ/1gSs+lZREJ6WLJjvLPvFTkhPv3b7++MGfFRVfiLCU76FsbGXvbOGgptlfU8+wK7u7spq6jHD1lBS+RQh9yH2FVXxOT0p4m5XxcN0qhJAMq7McW77fgIHjpkyzdXDW7CfZ944m7hYXG3OxcXy+esFFDOi5lbK/ntuK4adPKUa8slLzBgoLCwsIDLSwsj34xxllFTqvNaY/eew6fjTVM1UEGJbx5NHVc6eSOTe+CTC8plpBqaehxRhTOye9oRbCJ6XwSgrffsh9lZv5V9aTB5mP7iKEOnfpOthoiNssH1tHFwm9vmus2tl4iOmNe2kNHt8TumnfrtBWh+ZKuuY6rhmC4adPaetTSnKlZVIQn3UXO30WL98Quove3wccxwdrqvTW0DrPeUxFJQIMS7oZd3jXlo+F76urBPI9lK0nTrdwctMaZNJVTnTdqrU4/iE362XavTsX/nj3JpPVWVahRw+XqZ5zlgWoqknSNA9j1c7d5ORyi782eNxAS9XG2urq1au0VAWACEB8dpSkxyeO44MGDaL9Ymcd1/GjUx4kJ738SG5rDMfxW7HRdampqWc4xmO+pbO7ghL9bb4qAfY26/mjhEs3Th9ECCmpqLpM9fTx/0Ui2qPEzb1TM3Pr3zaglM836qd2+/ZtBwcH2ioDgGIQnx0l0fHJ4/EGDx78+fPnK7eSR1hQsrx4m8RcOLvUZ2booVMubjPI2ue7vJyIfduvXTpX871KrFKzsVocz8l4TOSoDKuzhna/NaF7h1vZiXP/GBGf4RGRUz3+/7/sjyMHNq/9ubKyUpwrB6CDID47SnLjk8fj6Q4Y0K2b3PWkFNHccaxlpXz+SOMBRqZmR/97b8sOSvmTE7R03pfyMrym2n35OtvJM1X7apGyZ6pVCbAniXExh7d/yM2Sk1eYMX/xwp/WiOeQXSI+Fy5dEfLrnroHp08cW8IrzMrKav55AEg8iM+OktD4/PPPP51dXLp1k0t6zKV3oFCdMeYmRTzejfTcDuYE0U+7e+NqftFHTT3DGas2G1s5kD4OSDTevc68E32caIxO8vBevXWvuIWosWpnp6nTSz68jUv895bdAgzTVVPYu3fvihUr6K0NAErBmrdMdO3aNSdn52FmI1Mzc8UkO+sWGOpgPEQe3Wep2/OXxbN66QwMiUrcEZs2dNR4Cc1OhJDmQMO5a3efSP80dcnqmKgTtvp9NvgvFGBY688UIZ0BA1+/eln3z8yMZ7Kysl5e9F9HFwHiA5SxGH76lIL4FEfXrl2b6uo6fITFmZjrYtKUKebxtgYHOU2d3pEFhlL+5Ngba+1Yt8p0tMs+Tta649cHmlqSWCSNusqx3ZeuPZH+adLClTFRJ6wGqIRtWdvq7bVFRkNHt7r6e109cZcvdunaVUU8/iyjmiT2PJGI4adPKYns1Wwryeq8JbLTwsr2TMx1MRnWgeP4sIFaeG1tuxcYepeXs2zmlIKc15p6hst2HtccaEh6keKjSoCd+y34xumDcvIKPwVvm073jWiMVTsnPMl2Gj6gbvCtgZaqs9P4yMhIegsDgGrQecsgYpidCKEdW4KJW2G3oyQcx8O2rJ1sbVxWVrr6SMyO2DTpzk6EUFc59ty1u/dxsoaMGr/152WTLA3f5eXQW1I3OTZCKD83ByFUyueXl5XOnTuX3pIAkHQQn2JEPLMz8/nTI/t/8/05WG9wm9cIzPjrsbWeasS+HVOXrD6UnD901HgqKhRPqn21/HafXH0kpqysdKrtEHr7con5qa+zXiCE/rx7u0uXLhYWFnQVA4B0gPgUF+KZnTiOuzqN0dDut8A/qK1PXDXfc94U+559tPZxstyXrpXcwUEdMXTU+EPJ+fbTfSL27RhjqFFSxKOrkh6Kys//eoIQ4iTE9+3bly0e19QBkFwQn2KBw+GIYXYihBZ4uVd+rQiPimtTVe/ychxMtG5fu+QduC308iNJmcpJERkWa+7a3Tuupv2NkNPwAZFH99FShoWdw9uCPITQlehzEydOpKUGWjB86CnDT59SEJ/0y8nJmTxlihhmZ8yFs7cS4kIPnSLu2iGkyKP7ptoOwWv/3sfJGj9zCTMbnY1pDjQ8lJxv5jh5x7pVs1xsaenIramufv+2ACHk6ekp+qPTpaamhu4S6MTw06cUxCfNeDyeyZAh6n01xS07S/n8n1csGWk7RvjF+YgO2x3rVpk5Tj6UnM/wRmdjMiwWcTX0edojS92eIu7ItR3rkv0mKy31YZcuXYYNGybKQ9NLrH6tRI/hp08piE86YRg2cNCgbt3kLifcFbd3+dTxo2S7dN178rKQ25d94k+yNLh97dKyncf9dp+ERmdzho4af/ThO9mu3ZzM9P68FS/io3MS4hUVFcXtzQaAJIL4pA2O4yYmJixWZ/FZk6/Oto1rct683nviopCLNpQU8cYP0/38+UtIVKL1RAZ1DLaPgpLKoeT8gaYWP83z2Lp6mciOW1v79593bru4uIjsiABIMfgjlDampqa5ublXbiWLW3bm5+Yc2f+b14JlQ82FusHLk4d/Lpo+oRtbfte1dPG8TYoYkmGxgo7FxhzZfn7/lurv30P2HqP6iEPMLKq/V1V/r5oyZQrVxwKACSA+6bF06VIul3s6Ok4c7kFWH47j42zMNLT7BWzaKcz2f96K9587Td/MOuhYLHTYtokMi+W+dC1CKHr/Fl7hu4Pn4kXQpyrDYpmbm1N9FACYADpvaXDt2rWIP/4IWLPBfpwz3bU01KaZKpCdHee+dO3qIzEpyXdcbYeIYDhuV8YsdVuH4TM3GH76lIL4FDUej+fu7j5ipNXKoGC6a2moTTNVIDvJMnTU+JCoxIKc15MsDahOUL0BAyjdvxiSoPWuqcDw06cUxKdI4ThuYGgo313hTMx1umtpqJjHC1i+yGnqdGFmqkB2kmugqWVIVGJhQZ6HvRlFCdqrT1+EkL6+PhU7B4CBID5FatSoUeVlZdeTUsRt5gCO42Oth3ftJrdhz5FWN4bspAKRoDmvMmc5WVGRoMRbTlVVlfQ9A8BMEJ+iExYW9vDhw9PRccRNo8SK71wv4p4qrc5UeZXxFLKTIkSCvnievjlgCUWH6N27N0V7BoBpID5FJCcnJ2jNGp/Fy8VwuNDthPhbN+KDtu1t9Z4qJUW8ORPtIDupM9DUcu66PTFRJ3asX0XF/rt06ULFbgFgIPgEFAUcx02HDZOTY28I3UV3LQ0V83i+871NzS1nLFja8pYCDJtqYyLbTQ4WFaLU+JlLKr+UnT203dxmzKixsMRBR3Xq1InJw2cYfvqUgtanKLi5uX2tqBDPS55TxtkKszgfjuOejub/dPoPrI0gAlMXrdY3s145fzrt99mWAgxfM53hp08piE/Kpaam3rx1a8uuMPG85Pk2P0+YS56bA5YU5LwOPHgBslMEiDWJVDV0PBzMBRhGdzmSTdz+ZhUxhp8+pSA+qYXjuOPYsSNGWs32WUx3LQ3FXDgr5CXP65fPXbt0zn35uoGmlqKpDciwWBtO3xJUVsx0tqa7FgBAEyA+qUV02x7844y4/Q1IzPIU5pJnSRFvw8rFA00tiBXmgMgoKKkQU1nCtsArD4DYgfikUF23rbgtCo/juKO1adducsJc8pzuMKJL125Bx2JFUxuob6CppfvydScP7X2V8ZTuWgAA/wPikypEt62mlo4Ydtt6THD8VFISdfNhq5c8w3/dUFpStO54Agy1bZMqAVZL0tIHUxet7qXZb86k0SJYERcAIDz4TKTK1q1bv1ZUcB4+Fbdu25MRh1MeJAuzsO2rjKcnD+11X75Oc6ChSEqjUEnhWz+HQc39VFPPsO+AwXPX7qo/MOp+3LkDgfMQQueyvjV+SnT41uj9WxBChhaj62bBVgmwV2n34k/sy3x0FyE0fpbvGPd5HXz1ZFis1Ueu+DkM+mnutP2RMR3ZFTMxfOYGw0+fUtD6pASPxwsNDfVZvFzcRtvm5+asD/QXZmFbHMfnTXHopdlv6qLVoqmNRu/eZD6MvxA4yezd60xhtm8yO0sK38417bl90dTMR3d/6NkLIXTj9MGfJ5ststbuYEtUta/Wsp3HH/7JSfmT05H9MBPDw4Php08p8WoYSQ13d3c5tnzQxq10F/I/BBg2zsZMU6f/1vATrW68fa2/oLLi1yuPpazbdh8nS7WvVv1H3r3OfJl278SWlV8+FR8InLcjNq3lPdRlp/vydfWHUyVfjUQIWbp4TPkxUHOgYS2O52Q83uPn9eVTcU7G4w4OWrae6Jl0+ZTfHPeHOZ/ErT8DAGaC1if5UlNT/0pP37xzb6tXFkVsvK1Z5deKP67cafXz911ezqXI4+7L1zVIGqmkOdBw/Mwl7svXIYTevcmsErQ0z7K57KwSYMTjRHYihGRYrIGmlo4zFiKELu7f0vE6f9x8qEqAUbccLgCgTSA+yefq5qappePq4UV3If8j6KelOW9en4xLUurZ+jDg+VMd5LorMKHbto7t5JnEN2+znje3TXPZiRBKunwKIaSpZ9jgSiex28xHdyvK+B2sULWv1tx1e2IvnGnfUkTEyCN1dfUOlgEAIEB8kuzw4cMfP3zYub/1236J0u2E+KhTf/j+HDzU3KrVjSOP7ivhfWDaaNsPua+Ib7QGmTS5QQvZiRAaP3PJuaxvjTt+05P+vbErKas1OXouZCv8MGfS6HY8t/hjIUJoyJAhHS8DAIAgPsmF43hAQMAUd88RFq2nlMgQi8IPG2m9wD+o1Y1xHD+4a4uli4cUjLYV3uv0h4fXLkYIjZ/l21WuiS73uuwcP8tX+OUjanE85vB2hJChRXsCrzEZFmvVgfOlJUXXL59r3x7guikAZIHfJTKFh4djGCZWI4ZwHHcaNYJYFF6Yj87NAUsqP5d7/rRJBLXR4uj6JfUbgi9S//zyqZj4/oeevZzn+DV+Sl12IoQKc17V4riQ7fKYI9uJnf+4+VBH6/6vgaaWhhajN65cPHaSO2ShMBg+c4Php08p+PUjDY7jISGbprh7itVklQVe7kW8j9FJ6cKMYxJgWPzlC9I9YoiYkdmApp6hz8YwXeMRTeYikZ2aeobv3mRmProbc2S7MA3QutCdu24Pua/nj5sP+TkMOvfHwZk/NhH2oAGGhwfDT59S0HlLmvDw8PLyMrFqep6MOHwrIS700KlWF4UnLJo2XrZrtwnz/KkujEarj8Ts42Tt42Qt23nc0sWDeHCMx/zmspMQEpW4IzaN2D56/5bX6Q9bOEQtju9bNafuQun4mSSPlVXtqzV+lu/ujb/AzVgAoBHEJznEsOkp/AoJhHd5OZnP/pq7dleTF/+khnp/fdW+Wqp9tawnevrtPknMVzmxZWXMke3NPSUkKpGYtfnj5oPEegh7/LyaG0lbi+OhCyc9jL+AEJq7bg9F6+x7/rSpFq/ZvrYNf+iU8D7IyLD69u1LRT0AMBDEJznErenZphUSCGuXzesmr2Dh5E5lXWLHfena8bN8EULR+7fcj2t6PE7digdd5dhrI64hhL58Kt63ak7jLasEWOjCSUT/cEhUIuntzjpd5djuy9fFXjgjfAP0Y+E72S6ycLkUALJAfJJj7dq1YtX0FH6FBEJd05NRk1UIs1ZvJ9qUBwLnlRS+bXljzYGGRNxmPrp7I/J/BgRVlPFXjDUgVuzbx8mi+t6oE+b5t7UBCgAgEcQnCTgcDoZhywN+obuQf7VphQTCzuBABjY9CTIsFtGmRAgdXd96e3HW6u2aeoYIoRNbVtZfI3ffqjlfPhX/0LPXztg0EYy9akcDlJk6depEdwl0YvjpUwrikwRz5s7VNzDSNxBqeA7VbifEnz19QsgVEggCDHuUfGf8rCUMbHoS6rcpm+vCrSPDYgUcjCa+3+ozgVjk7+mfN4g+Wwsnt4wHiffjzjX4anktwPYZO+PHWrzm8pkIYTauKC9ji9kqkqLB8KGnDD99SknVx2VBQQGbzVYR7b2pnz59+vHDh4PHz4ryoM0hVkgwNbcUZoWEOtvX+ldXCaR7wG2rZq3e/ijhErFq/HD7iS2PnyJugXIgcN6XT8VH1/v67T55dvd64kc3Th9s8in7OFmkj8lSUFKxdPHYty3Yc75vq730z9IeDdDVJbcAAJhMeuIzKirK29t7xowZUVFRojxuUFCQck8VUzNzUR60STiOj7UeLvwKCXXPun39qvvyddI94LarnBwx7aSrnFyTGxBduFeO7kQIPUmMs57oqaKuWTezpTHriZ45GWkVZZ8QQu9eZ/YdMLjvgMEtF9ChE2iG50+bHsZfePIgaeQoByr2DwBojpQsSMHlco2NjRFCTcYndetuYBimoqq6dlPo/EXLqNh/m/w40+Pa1UvRSelCzvIkRB7dt2Pdqsa38QKS4udJZtjnT3cz3ze3gbFq54ySmtWLZn4uevfwYUszVgGQMsSlX4o+/6Xh2mdUVBSRnaK3devWbwKB58x5tBy9vpgLZ69dvRS0bW+bshMhdPzAbkOL0ZCdkstnY9iX8rJWb8PyKIljYtL0avgAgHaQ7Pjk8/mOjo7e3t50FXAgPHyKuyft9/Us5vF+8vVxmjp9xoKlbXriu7wcftHHacvXUVQYEAFd4xFd2d0j9jW77APhc3npoEGDRFMSAEwg2fG5YsUKDoeDEDpz5syZM2dEfPTU1NSvFRVzf/QV8XEbwHHc0dqULa8g/AoJdQ7t3CTfQ1nXeAQFdQERkWGxxs9acu3SOeKOnqABhs/cYPjpU0qy4xMhFBISUlJS4uVFw72pt23bptxThfZ7ky3wcv9UUhJ182FbF5TBcZyTEMvk+SpSY+yMH2u+V92KjW5uAyJZFRQURFiUuJCO4R3txvDTp5Rkf26eOnWKrkXIcBznJCb6rgig5eh1Yi6cJRaFV9fUbutzb8VGfxdgY2f8SEFdQKQUlFQ09QwP79ri7OrZ5AbEvbJHjybntqMAACTp8UnjAp4XL14UYNicBVQtaiqMukueQi4K38DhXVs09Qzr3/xSalQJMGJKSZOkcpzUjFWbdy31EGAY7VfiAWAIyY5P4TW4APDPP/80ns2CYZi8vHyDBxtvhuN4586dBw0apG9gpCzaJRoalOE0agRbXmHDniPteLoAwz4WvvcK2EJ6YaJBBGT2s5QvpfwXqUkVZZ8qykvLij8SP62pErTw3M7d2OiffxBCnWW7qPfTQwgNHmHzg3KvoXZOMjIsCQ1XfTObWrzm2G/bVqxr5b4FxBu4fW9+WjZrsmBaNoNXSSJeJZFd7mVKfDa+AND4ETabLcxmLBarsrKSmO5JbpFt4jvXS/j7YDeWdDOuukpg5zqb9MIoUlL4NvtZCvfR3XdZ3Pe5WZ3+859qQSXxo166Bqo6+l1VNKzn/7vWUg81je4qfZrcT7Wgsjj3BfF9WWHup4I3CKGrR3chhE6FBiKEOneVk1fooaE3eIjNOH0za009Q4m4NtxVjm1oMfrKuVNNxueztEeyXbpqa2sjhFgsVrvf/LRs1mTBtGwGr5Iwm9H+KtX/KaVRKgGfC2IoNjb2m0Aw2XU6XQXcToi/dSO+HbM86xzetcXQYrSYrzSU//JZyo1Lr9Lu57/KqGtQGjm6O473UuyjraIziK3YU7Zbm0+hl65Bg0fcNx9HCJV/fEuEa+GLNH5+FpGmMrJdFHoom9g4WrlMGzBkpDi/Yi5z/XYvm172id/krQJkZCR+nGD7ULdqikRg+OlTCuKzPTZt2kRjz60Aw4iFbds6y7P+Hj4Wvv8p8FdyCyNFRRmf+/DOjchDBa8za75hCKF+ZnaOSzdpD7VWVNduR1gKT7GPFkKol66B8TgPhFAtjleUfHjPTS18kZZ680LSpZP/YXX+QUnFxMZxnPdincFDqKukffTNbPCa6pio4z5+P9NdixhheHgw/PQpBfHZZjiO5+fnr9+6g64CJoyx7NxZ9tfDke3eA9Fz299oGIlVdVBFGT/+xL7kK1EV5Z9qa6p76Ro4+oZoD7XuqT2Qrr5TGRZLsY+WYh8t43Eezit3YuX8Dy/Tc1I5SRcOJ1062ZWtYGBuM37mEiNLe1rKa4zov7167lTj+OT+9bhP7960VAWAtIL4bLOkpKTv37/T1XN77GBY1svM8KhY4e/l2djZ38PFZMxtlQC7d/XM5UPbidTsZ2bnvHqftqk1pa3M9mErquhZjdOzGjfOL/RTwetXf8amXDj815142W5yls7TxKQ9auc6++j6pY3H35aX8lVVVemqCgCpBPHZZvv371fuqUJLz21+bs7W4CCvBctsHJzavRMcx7Necmes3ExiYe2Q//LZ2T3rXz15WFMlEOfUbEyGxeqla9BL18DOJ6g45wWRo0mXTiqq9pm+YoPNZC8ahxoZW9lXVwmSbsY1NwEUAEAWho4m6IjbHM48Ohbqw3F88lgbtrxCwKadHdlP9kvudwFmaudMVmFtUovj108eWGLbP8jVIufFM8elm35OyJuzP07PapxEZGcDRIgGxuf6HLmlqKV3eO2iWUOUDq9dXF7Co6UeBSUV9f6DEi6fb/D4g7u3zM3pv6ceANJEelqfWlpaM2bMsLS0pPQoOTk53wQCp0lTKT1Kk3ZsCf5UUhKdlN7BxSKuXYhEdCwdUCXA4v7Ym3AqXFBR3s/MznXzCU0TCxHXQBEZFkvTxGLO/rjqb9jDqH13j21LvhrV39DUd3tEb63+Ii7Gwtn9+vF9DR6s+FxOzFoBAJBFeuLTysrKyory5Wf/+OMPhNCAgfpUH6iB/NycI/t/8/05uN0zVerciI1u4S7QVKgSYCe2rkq+GvU3XmPk6G6/ZCMxxlX6yHZj2/kE2cwJfJF4OWHv6sCJw/sZDFny6++iDFEji9HR+7d8eFfQjnUcpRLDZ24w/PQpBZ23bRMTE2NjZy/ixQKJblsN7X4L/IM6uCsBhn35/NlmUnsW+WuHWhy/eGDrgpF9ky6dHDUvcO3dIvfNx6U1O+vIsFjG4zwC4rInrw0vLMgJnDg8eMZo3ttc0Rxda5AJQujujdgGj/fvL+p2sJhgeHgw/PQpBfHZBjiOv3v/fvqsuSI+LtFtGx4V1/HYznjySGRTVh5cO+9jrn7pwBb9URPW3i2y8wmSxKub7VY/RN/nvVnlMnR/4LwqAUb1cYnpKzevXGzwuKGhIdWHBoBRID7bgMvlCjDMzJzay6sNEN22Qdv2ktIXlxBzXr6HMtVTVnhvc5eM0t0fMHeA5Vj/y5num48zKjjrI0I0MD53/IrQB3HnFozse/3kAaoPOtx+YtaLDKqPAgDDQXy2wcWLFxFCGlraIjsi0W0r3/2HaXMXkbLDv1LuG1tROM2/Fsf3B84LmGD6vfq7z5FbTOiqFYYMi2U+bdHau0X6oyacCg1cZK1NaV+uzmCT798EZZ/49R9kw51YACAVxGcbxMTETHEX6XS6sJ1bP5WUHLt0k5SrrTiOFxfxTO3aP2e0ZTnP03zM1R/EnRvnty0gLltqBtaSRbYb233zcd/IlM7yPVa5DD27Z30tjlNxIOLyZ0pyYt0j3eTkVOi7OxAAUgnisw3evX/v4OQissMV83j7d2/3WrCs46NtCcSMT72h5KdaLY7v9HXfOGvsD701/C9nmk9bJBF3KaFFL10D38iU8StCrx7dtdhGh4pmaFc5tnwP5cazP5lJZLevEk8MP31KQXwKq6CgQIBh+gbkJJkwZkwZ3/FFEur7K+UeQkhZTZ2sHRJ4b3N9zPv8dSd+7PItvpEp0FvbKqIv1/9yZi1Cq5xNqLgaaj1x+suMdNJ3K4kYPvSU4adPKYhPYV29ehWJcMZnzIWzWS8zt+yPIHGSzM0rFw0tRpPbLrx+8sAql6H/ke1CNDpJ3LPUU+yjFRCXbea28FRoYPCM0eR25Ooam1VUfMH/u09dps5aAYA6EJ/CSkhI0DcwEs2MTxzHV/v7jrQd05G1bRv7UPiury5p8V+L4yGzxp4KDTRz9QmIy4ZGZzvIsFjOK3d6745+8zRlnpkaiUv9atnW2QAAIABJREFU6Q21+C7Aij8WEv+Ul5cna88AAALEp7C43ExLGzvRHGu13+LKrxUb9hwhcZ84jld8/qxrbEbK3spLeD85Gb9Ku+cWEuG8cidc6ewIPatxPyfksbp2WzpmYHrSdVL2qaDUEyGU9iCJlL0BABqD+BQKjuPln8uHmo0QwbHyc3Nios/7/hxM7qJrxR8Lq6sEpIwb4r3NXTHOsLLyq29kCnFnadBBbEWVgLhs7aFWu5Z5nt2zvuM7JEYPpSbfIf6poKDQ8X0CAOqD+BRKQUHBN4FANAsm+C+a102O3fH1+RogVnHr+Lih9KTrv7hZKvTSWBb1uJeuARmlAYQQkmGx5uyPGzUv8OrRXYfXLu74Do2t7POzXxPfOzlRNVsJAMaC+BRKfHw8Qqi3el+qD5R8h5OW+vDnLbtJv8jKe/cWIdTBXtb0pOu7lnmq65v6RqawFWEeIfnsfIKcV+1KunRy9RTzDg4m0hs6Mjc7i6zCJBfDZ24w/PQpBfEplDt37gwYqC+CcUP+i+f10xvk4kb+ku4Z6Y87eKOVK0d27lrmqT3UauZvMXCxkzrm0xb5HLn1NivjJyfjjiRoLw2df/5BAozyVXbFHMNnbjD89CkF8SmUV69eGRiZUH2UYwfDingfdxyNomLnuTmv+/TTa/fTb589evHAllHzAufsj4PspJqmiYXPkVsl7/M7kqDq/fWrqwTlpfzWNwUAtB3Ep1AKP3yger0hHMd3btnoNHU6WWsMNVBTXa2m2c7Jf+lJ149vDbCdG2DnQ/IVWdCcjicoMfj25tVohJCjoyPJ9QHAeBCfrePz+d8Egr4a1M5rPH40vPJrhd/aLVTs/MO7gnYv10dc7xw1LxCyU8TqEvQX15HtSNCucmyE0DdBJUJITk6O/PoAYDaIz9YVFhYihHr3IXmtu/pwHN8Tutlp6nRyJ6vUaXcHIJGd2kOtIDtpQSTo+zcv9q2a3Y6nG1qM/vMWORNJAQANQHy2Ljk5GVF8n7LjR8O/fC6nqOmJEErmXEcIqfZtWwOa9zZ3t58XMVaImrpA6zRNLNxCIlJvxrRjNgvVN3aVCAwfesrw06cUxGfr0tLSlHtS+DFEddMTIfS9qqqtTykv4a2eMkKxjzaMs6Wd8TiP0QvX3I87/+Ba226iojd05Lv8HIqqkhQMH3rK8NOnFMRn6/Ly8mzsKLzFNNVNT4RQ0o24Ns1aqcXxDd72rK5y8w8lQHaKAzufIP1RE8J/WdimG5zJ/6BYW1uL4F7ZAFAAPhlbl5X12mr0WIp2LoKmZzvs8fMseZ/vfzlT4tZGqP6GYeWfmvwRW7GnbLcmUqT849sWfoqV86u/CRBCst3kmnw1Mm5e6KGmIYJ7g0/dcIyXzV09ZcSBxCwhe2VV1DW/fxMghOBe2QCQDuKzdVXfq7T761K0cxE0PRFCPN5HWwsHITe+fvLAX3fi3UIiJPEmKlnJ8Zc2+LSwQT8zuwbd0XtdDRFCbiERjdfvfff8UcSisQgheSVV3zMpjfeWFBF699g2I0d3EcSnDIs1/1DCfs/hG2Y67Ir9S5heAUVVCse7AcBw0HnbCmLWCkV3ycZxfHvIehE0Pb+UfRJy0ifvbW7kjl/MPRZL9Frw8kqqRo7u9b/klVSJH+WlJUX+NLX6W+tr8TTIzsZNz3fPH909to304lvAVlTx2nGOl/fmwr4QUR4XANAYtD5bgWEYQkhevjsVO3+YnCQQYD4rVlOx8/pkWJ2F2awWxzfMGKOkoTvOL5TqkiilM8zWffPxBg9i5fw/j+9IvXA4Ly2pIP2+ntW4FvbQcnbW4vjNfUGpFw6TXnmrNE0sRi9cE/fHTjP7Sbomrdx+ruN3CAAANAdan624e/cuomyx+OUL5wwYbETRMkN1yj7xBV+/aOoZtrrlHj/PirKSmXsuSeVwIbaiyji/UKIZmnHjXAtbtpydxTkv9kwaSGRnXaNWlGzmBCr17bdl/oRWp/NK5f9jWzF85gbDT59SEJ+tqKioQAhRsVh8fm7Oly+f123fT/qeG/gmwBBCXdmtNKBznqc9v5/ovGqXJF7yFJIMi6UzzBYhVJL/qrlt6rKzn5ndytjXjfts753cVVlWIq+k6nPk1rgVNDTTZVismXsuVWEVe/w8RX90icPwmRsMP31KQXy2Ii0tbcBAfSr2HBIUICfHNhpmTsXO26oWxzfNdeoz2NR82iK6a6FQLY7n/5WMEBo2eV6TG9TPzuYmvMop9nQLiVgZ+1oEw4Wao9hHy3nVrmf3buc8T2t5SzUtqka9AcBwEJ+tyMvLo+JeKwIMS7pzy2vhUhHcBE2YFfuObVhW/Q1zDT5KdTF0wcr5xTkvjsy1riwrQQgNnTCz8TZ12YkQamGxCOeVO43HedDeLzp8qo+ius6muU4t///2MzQVWUkAMApcGmlFbe3fVOz2zMnfq79/95jb5mXY2qHVFft4b3OTr0aNXrhGarptubejubejm/yRvJLquBWhjad4Pr12Oi8tqe6fT2IixLwhTnTh7nU1vHx4+7Rla+kuBwDGgdZnK7iZXCpuVXZo7y6nqdOVqFwLUHj7V87u2v0HmzmBdBdCrV66Bj5Hbq2Mfd3knBwiO/uZ2fUzs0MIXd8dUJzzQrQFtpliHy1zj8Uxh7dXCZh+T2wARA/ikwaZz58W8T56zl9CdyEIIcR9mJj3In3Gr1G090aSqJ+Znf/lTOLLNzKFSMTinBdVlRUtnCZxvXPGjnPEeNpTfpOEmR5KL4clG7uwu2+dT+3NaCUaw4eeMvz0KQXx2RIcxxFCP/zQg9zdHj8SjhASk0FD+wPm99I1oHEUDBXYPXoq9tEivnrpGsz8LYZI0DOr3N89f9TkU+rGCsl2Y08PjUQIVZaVxG5bJsqy20G2G9t55c78l8+bWwu3Tz89EZckbhg+9JThp08piM+WFBYWfhMI9AYNJnGfOI7HxkT7/hwsgkFDrbp+8kBFWcmMHW27j4fEkWGxZv4WQ7QpIxaNxcr5jbcZOmFWXcNU08TC3GMxQoh7Ozrj5gVRltoOBvausmz5Q78saPKnQq42BQBoK4hPUXuYnIRVfh3jPEVkR2zubmW1OH5u70YjR3epGTHUAhkWa/a+WOL7k8sntrp93QILlzb4EGvKiy0ZFsvJf3vei2cv0+7RXQsADALx2RJciCkfbfXH4QOKyipUrzRUX9GH94YWoxs/fvPM4epvmP2SjSKrhF69dA2INmVxzotW25QyLNaC3+8Q3/++YIwwk39oZGDv2qOPVthPs+guBAAGoSQ+MQwrKCi4fv16VFTU9evXCwoK+PwmusvEX3x8PEJIQ0ubrB3iOH7vzzszfEQ6aOgbhjV5fyvmND3r1G9TNtmFWx+xNAFCqLKsJCZkoSjqay8ZFmvc8q1YxZf8l8/orgUApiAzPnEcv379uomJiby8vI6OjouLi7e3t4uLi46OjqqqqomJSVRUFBXtOQnyMDnpmwATZc8tQuhtXnbjB9OTrjOq6UmQYbGIYUEIoYQ9P7e6vfm0RcSYI+7t6DcPblJaWwf1N7eXZcv/scmf7kIAYApyRq9gGLZ79+5Lly65ubmFhoYOHjwYIfT161cul5uTk3Pp0qWMjIyMjAxvb29vb+8zZ854eXmRclyJczX6XA9FZVH23DbnyDrfXroG0tf01DAydwuJ6KGm0dwGmiYW3rujqyq/IISqv2Gy3dhuIRHEE5vc3n3TH7mP7yKEiKe09XAiQ1wBjdm8uKKML+TNtBmiU6dOTB59yvDTpxQ58SkQCOzt7YODgxs8bmRkhBAKDg7GMOzZs2cbN27kcDje3t6pqam7d+8Wh6GnLfv27Ru5O7wed8XF3ZvcfbZD/stnXz4Vu289TXch5CMmq7S8TYNblbV8Z1O2okoLGwhzOJExsHe9tMEnavf6xVtpuJOa2GJ4eDD89ClFTuetioqKlZVVCxuw2WwrK6vbt2+XlJSEhITs27dv2LBh4n9BNCkpycbOnqy9vX9b8OVz+bjJ7mTtsN2iD2yRU+ypbtDK3SKBZJFhsUYvXHM/7ryYD3QCQDqIeuStiopKcHBwSUmJgYGBsbGxmCdoRUWFMnnr6l25dB6JwWoJVQLs2X2O3fxfpGmZIUAwcfLCq6tSblyiuxAApJ9I45PL5RJ5qaKiEhUVtWTJEvFPUBJdPHNypO0Y2rus7109U1v93dDBld4yABWIVZZO/bqa7kIAkH4ijc/3798bGxt36tTpwYMHCKHg4OAlS5Z4eXkxYTgujuMfCgtHj59EdyHo8qHt/czsGt8FGkgHhyUhlV/KK8qY8lcpAHQRaXw6OztnZGQghNzd3YnIDA4OtrGx2bZtmyjLEN736mqyOm+zXnC/CTDbsTQv7V1Rxq8o/2ThKe5LuYJ20za1rq2pjj+xj+5CxAXD10xn+OlTStTXPt+8eYMQKioqKiwsJB4JDg6ePXu2iMsQUn5e/lCzEaTsKjbmIkJIXVOblL21W/yJfbU11dqm1vSWAagj243dz8wu+UoU3YWIC4YPPWX46VNKdPFZUFCwadMma2trhFB8fLy2tnbdj+p/L1aqvje9Wmw73LvDGWk7hqy9tVvylSgjR/fGN4sG0sTCc1lF+SfovwWAUhTGJ5/Pf/DgwYEDB7y8vDp16qSjo3Po0CE/P7/4+HhnZ2fqjiue3rx+RfuFT6Ln1ni8J71lAKr1N7eH/lsAqEZyfGIYdv369RUrVvTu3VtVVdXa2jo7O3vChAkZGRklJSU8Hi8sLIyB2fn+bcE3gWCYpS29ZaTejIGeWyaQYbH6mdk95sTSXQgA0oy0+IyKiiJWuw0KCnr58mVRURHxeFJSUo8ePfT19VVUmDvU80nqI4RQX61+9JZxN/pkL10D6LllgqETZn36WFglwOguBACpRU58btq0ydvb283NrbKy8vnz57dv387Pz58xYwZCKCMjw8XFRUNDg8nrxd9PSuyhqCzHpjO3anH8fd7rYZPn0VgDEJn+I0bXVAmyn6XQXQgAUouE+OTz+Rs2bEAIBQcHs/+bENra2lFRUZWVlSEhIQihoqIib29vDQ2NAwcOYJhk/EWM4/jff//dV4OEFU3TUh9a2Dl0fD8d8e5NZs03THso9NwyAltRpdsPSjciD9FdCP0YPnOD4adPKRLiU0VFxdjYWE1NrXHjks1mBwcHV1ZW7t+/X01NraioaPny5fLy8ps2ber4calWWFj4vaqqdx/1ju/qQ+F72md8Egu59dQeSG8ZQGSMx3m8SE2muwr6MXzmBsNPn1LkdN5GRkYihDQ0NIjlhBpgs9nLli17//79mTNn1NTUEEJEa5UhSvl8HMf79NWkt4zHnNh+Znawzi1z6Jo71NbW1lR/p7sQAKQTOfFpZGTE4/EyMjIUFBSa24bFYnl5efF4vPv37xsbG5NyXIlQUfGlprpatTcJrdiOKCspGmQ7gd4agCip6AyqqRIUFeTQXQgA0onMtoiKiooww2utrKyeP39O4nHFHOdmPKJ7vaEqAfa9sqK3HoP+agHEjUjfPEuluxAApJOoF+1joOd/PVFUpnnSTjmfhxBS1tSltwwgYkaO7mUlH+muAgDpBPFJubcFebQv14d9LkcIwV1WmEbDeGRZMY/uKmjG8KGnDD99SkF8Ug6vqaG7BFRZUW7k6E53FUDUFPto1+L0v/3oxfChpww/fUpBfDartLQUISQn19G1Dt68fkX7rJWa6uqe2nr01gBET0Vn0N+Mj08AKALx2Sw+n48QUpb8tQbxmpoucvJKffvTXQgQNWL0EACACtTGp8hW6eNyuVFRUVwul8Qjvn79mpT9dPrPfxR+6EHKrtrn779rK8tKNIzMaawB0KXbD0p0lwCAdKIqPvl8fu/evTt37mxiYrJp0yYul0s8vmnTpiaXVmi3FStWdOrUydjY2Nvb29jYuHPnzuTuv4Pevy0QVFb209OnsYavFRUIIdlucjTWAP6PvXuPiyn9HwD+tE4XcyrCpCGkcq/cIioKE1u5hUXl1lqXVcqutfvNLimLr5XVzWXdd6lIJFRorOi2Ca1GKqJEmm2EtU0undbvj7Pf+bVdpzqXmTmf9x/7ytkz53ye6TSfee5sMR/N/i6zAKglutJndnY2uelKbm5uQECAlZWVQCAICgpyc3OLiYmhKsNFRUWFhYUhhKysrMjFdRFC9vb28mwNEEJ/19YiGHbLVdDnDQBN6EqfkydPRggVFxdXVVUlJCQIhUKJRELm0bCwMHt7+/YvHC+VSj09PRFCCQkJd+7cIRfXFQqFCKEFCxa0vwhq46/Xr9gOAbAG+rw5PnOD48WnFTXps2EuxDAsMDBw2bJl2traLi4uycnJZB4ldzGjxPfff48QsrKyku+/jeP42bNnEUK5ubklJSVU3UjV1dbWwqwVztLR7YQQevv2LduBsIbjMzc4XnxaUZM+Hz16pKGh4efnl5iYKE+l69evRwht3bqV/CeO4y4uLlFRUR8+fKipqcHbvfllSkoKQuibb76pexDHcbICeuHChXZenxLKsMWphgaMr+Yuft+BCKHffoONPwGgGDUfrOQ+KmFhYa6urrq6ukOHDo2IiMjPz09OTkYIeXh41KueYlTs+5Gbm4sQsrW1rXecXHc3IyOj/bdov5JHRQih7j2MWYxBQ1MLOsA4rqgIFo4HgGLUpE8yY6WlpQUGBlpZWeXm5q5evdrKykpDQ6OyslIqlZqbm1M7nIeclNmohgmVRUUPChFFXxfa7G+CgA4wjnv3DrYtA4BqHyji7u4u/7mioqLRbk53d/e0tDRKbldcXExes7i4uN7/ioyMJO8lP8L0ewoAAEBpUJJ0GqKm9ikWiwcOHCj/J5/PJ7s5a2pqcnNzySopQig6Otre3l4gEPj5+VFyXwAAAIAVGlRVzvz8/Lp27bpx48amTpBKpdnZ2cePH4+OjkZUjAcjB2QXFxebmJjUPR4VFeXp6enu7h4VFSU/sw23Cw0NXbNmzbO/atsT5IE9oQHffJlbwea6o2PM+M5f7bSaMpfFGABbXj57HDLLouGfCXcQBMFu7wm7OF58Mk3Q1AZJ2ZjMnTt3pqamBgUFNXVCvSopVfdtSEnG3JKMBD3ZDgEhhN68fsl2CACwg8vJA3G++LSiLH1iGBYVFbV3714/P7/mZ2tgGGZpadn+O5INwvfu3at3nBxVpCQDiAYMHIwQKistYTGGD7Xvn+TCvAWO+ksK22UDQAsqZwTy+fzc3NyYmJhevXrJG07pM3v2bITQrl276h6UyWQikQgh5ODgQHcAiujIY3+lWZj3yWWvJE8QQsbGbE6dAkAtUfzByufznzx5MmHCBE9PT4FAEBER0f7F+ZqyaNEihJBIJJKvoEsQxM6dOxFCQqGQkgquetDW1nnz159sRwHYBC14AFCO+noJ2YobHh4ukUhWr16tq6vr5OTU1G5iYrE4KCjIw8OjDTcyMTEh58bY29uTt+jVq1dAQABC6MCBA+0viNrQ0tYu+i2Z7SgAO+6nJbEdAgDqia5mPR8fn+LiYjK9iUQi+W5iGhoaHh4eHh4eTk5OGhoaX3755aRJk9rc0vvLL7+QS/SRt5BIJEZGRpGRkZwdZNgoHdiqDHAYx9dM53jxaUVjk46JiUlUVNSBAweOHDmSkZEhlUrJXkmpVMrn8728vM6ePdvOlW8xDEtOThaLxWKxuKioyM3NTdnabAU9jbW0dSrKy3r2NmErBrL/9eWzxwY9+rAVA2CLODmW7RBYxvGFUzhefFrR3iOC47iPj4+Pjw99t7C0tKQja/bsScGcEwzDOnT46NnT0uE2du2/WttoaesghN5XV7EVAGBLrRLsWACAuoIxmU0i18Gvpm3oE2M0NDQ0dXh/PMxjOxDAtNcVZWyHAIDaoit9EgQRERHh5+cXFRVVd3n3ujuaKTlyrH/l8ybXpleQkaDn9csJVETUdppaWjCEhIP+kj7roKnFdhQAqCe6Gm8vX768evVq+T+trKyWLVs2b968UaNG2draikQicpMWLujStSvbISCern5ZQQ7bUQCmld/P/agDVlvznu1AAFBDdNU+Bw8ejBCKjIwMDw+Xb2FmaGhILhVE/pcjRoyyyUwRsRtDp27dXzx5CD1hXPMk9zc9A/a/vbGL40NPOV58WtGVPk1MTMgpJT4+Pnfu3KmoqCDzqEQiyc3NlUgkJSUlNN1a2Rj37vPqZSW7MegbdEPQE8Y94uTYHibmbEfBspoaNjdsYB3Hi08rGocOJSUlHTlyJCIiAiHE5/PJPFpVVZWWlpaQkMCdqZlj7cYjtpe97airh2npPBFnsRgDYJjspRQhZGYxku1AWMbxFZc4Xnxa0Zg+MQxLSkrKyMiou4g8juN2dnYuLi703VfZCHoYI4Qqylmu+XXpLoDRQ5xSWVrUQVOrW49ebAcCgHqiK33KZLKhQ4f26tWLz+dLpdJevXqJxWKa7qXkuvL5HXm8Z09L2Q1j5MSpMIOeUx7dTNHuiHfU1Wc7EADUE13pMz4+nuzjDAsLi46OlkgkVlZWTk5O8uXdOUUZ5q7YuXyC/tegB7hAnHx6kLVSbNsHgFqiK32S220mJCSkpaUFBgaSSxCIRCJ7e3uBQNDMrtpqqV//gYV5d9iNoe+Q4R20tB/euMpuGIAZ79/InpcU2rp8wnYgAKgtGkfeBgYG3rx5087ObuPGjeXl5cXFxfLBt+S+KErOxMREW0en/BkFfZZ2jhMe3S9ofhdxunXAMEPjPjkXjrEYA2BMye20jzDNkROnsh0I+zg+c4PjxacVjUOHNm7ciBDy8PAg04aJiYl88G1CAsstmQrS0tJ++uRx+6/zsesMhNCj+/ntv1R7jBZOf5SdArM/uaAoS9SpC1+H164tGdQDx9dM53jxaUVX+oyKivLz87O2tra1tXV2dq47yxPHcVUZeUsQ1EyZEvQ07sjDH+TfpeRqbea6xBchVJaXzW4YgG61BJEVs2+M82y2AwFAndHY9xkWFubq6rp69WqRSNS3b18/Pz9VWe2WchiGGXbvzvroIf0u/M6GgrtXzrAbBqDb85JChND4GW3ZhR4AoCC60ie5pC3Z2UkeCQsL09XVdXJyUqFV4yk02WVaUtxJtqNAw8ZNzorZB+236i3/2jkdXK/v4GFsBwKAOqMsfSYmJtb9J47jvr6+dTs7fX19EUIikcjV1VVXV5eq+6qKSZNdENtrDyGEpniuRNB+q9ZqCeLqga1jPp7FdiAAqDnK0ufx48frLWO7detW8geyszM0NPTDhw+5ubl1q6RKbtCgQTnZNyi5lLXNWJ2OvN+zMym5Wpv1HTxMv4shtN+qMfK7Efk9CSDODz3lePFpRWXj7dixY+s2zOJ4I6P+LC0tySophfelj7aWVvv3+yTxcFzQg/3FExBCE+YsyorZ9/4N59rPOeLulTP6XQyh5VaO42umc7z4tKIyfUokErJh1snJqd4u2QAh5DpzljJ0f7ou8dXU4RVcZz+RA8q9fyPLitnnssSH7UCUCMfXTOd48WlFy9AhkUjk6elpaGhILjCUnp5eb8UAchsWrpnu9glC6P49ltf+1e/C72U28NqRH9gNA9Ah58LxDlraHy9YxXYgAKg/KtNnVVXVhw8fiouLIyMjyc0+yQWG7O3tNTU1ySop2T964MABCu9Lnx49ehTm51F1NYuhw7W0dX5NPEvVBdvMfW3Q85LCP4ooKxpQElcPbh0yehyslgAAAyhLnwsWLCA7O01MTDw8PJKTk2tqauoteOvp6dm3b18NDY3c3Fyq7kurcePG5edRWVl0nOh08WwMhRdsG0vbSTq4XurPwWwHAqhUeifzzZ8vvDaEsB0IAJxAWfpsuJAQhmHyBW8rKirkVVIu+3Slz6P7BS8oGo7UHnN9N4qTY18+o2BJQqAk4ret7tStu6CPGduBAMAJNK55WxefzyerpFVVVeQEUG6yHe+opa19SQkqoFM8V/L0DTJPcLETWi39UZT3vKTQ54dD9Y5LSh+yEo/y4PjMDY4Xn1YMpU85HMdDQ0NVpRo6ZswYhFAldUOIMQwbaW1z8ug+qi7YZh0wzMFtAcxgURupPwfrdzG0tJ1U7/izR/dZiUd5cHzNdI4Xn1ZMp0/SjBkzWLlva3Xv3h0hVF1NZYJZ6fcVw+23fUz7PX1wr+HxeX4BH2Ga57bCJAeV9/LZY3Fy7KcBu9gOBAAOYSd9+vhw9yN7gtOUjjw8hsEK6ECrYaX3G9nsRYeHT/FcAT2gauBM0HKeXucxU2ChPgCYw076VBXGxsa4rt79gkaqbm2GYdg4h4nRh/ZSeM02W7BuG0/f4MreTWwHAtqu9E5m6e8ZX4QeZzsQALgF0mdzMAz7++/aP/98Re1lfb9e/7JSmpOVTu1l26ADhs3xXi9OjoU5oKorfttqw159G/Z6AgBoBemTBSOsR+vq6V2Kj2U7EIQQmuK5UgfXO71pKduBgLbIvRTz4umjdXtOsx2I8uL40FOOF59WkD5b0LdvX1ES9cvDLlm+KupgRLUS7HvaAcPWhkdXlj7MvcT+dBrQKu/fyE4HLLUY49Cr3yC2Y1FeHF8znePFpxWkzxZ00ten47Jr1n2r05F39eI5Oi5ez3ihC0Ko4mmT44MsbSeZDhmW+OM6mMSiWs5t9dHU7vhl2Ilmzql594axeJQTx9dM53jxaQXpswWmpqZ5Yuq3V+Ph+BBLqx++W0v5lRvqoMDfz/pDF2reVMMkFhXyR1GeODl26caQ5le4zRadRwjJlKCdAwA1A+mzBaNGjXpQmE/HlQO371KSAUQIIR0evth/uzg59n76JbZjAS2rJYgj3i4C0/6Osxcpcj7sHggA5SB9tkBfXx8hREcn5Qjr0UaCHnt3BFF+5Xo68nCE0FvZX82f5uS+vFf/IWeClstewketsosLXPZO9tfGoxfZDgQA7oL02YIJEyYghCrpWSTwOoJ+AAAgAElEQVRow/c//Hb917LSEjouLtelG1+7I6/RlRPq2Xzi2t9ETezGT2mNB7TT/fRL4uTY5YERBoaC5s+s/fc+uwAACkH6bAG5C1tVVQtVt7aZNuuTTp0NwrZ8R8fF20CHhy8PDC+981vWqZ/YjgU0TvZSGrl2Tr9hNoo021ZKyhgISclxfOYGx4tPK0ifLeDz+R15PGp3/ZTDMOxL/w1JcSfpnsEiMO5dlJutyJl2U+fZT5uXuPMrWEhBCdUSxOHPnbU64puOi9iORWVwfM10jhefVpA+W6ajrVPysIimi3su/gzT1Pxl7480XZ/UqVPn1y+eK3jyssAIw159j66eCp2gyuZSmP+LJ0Xbz2YpMpoaKdDhDQBoM0ifLRs4cEDR/UKaLs7D8U/mL9jzQxDdFdDXLxTNhR0w7PuT14i31Yc/d4bOM+WRcmhbVsy+r3bHKL4hdun9u9odebRGBQBnQfpsGU1TP+U27wiluwLq+PG0u5lXFT9fvws/MPLKn5Inx79woy8qoLjSO5lXD2x1nL14hKML27EAABCC9KkI+qZ+khiogGrr6KBWjsPsO3jYFyHHSu/8lnJoG01RAQX9UZR3aMVkgWn/lVtat8/d7ZSkXn360hQVABwH6bNl48ePRwg9eVxC3y3oroBO+Hg6av04zBGOLk7un109sBUyKItkL6WHP//YsFff4HO32vByTFOL8pBUC8eHnnK8+LSC9NkyY2NjhFD5MxrnAMgroDTNAVVw5YSGFvvvcJy9GDIoW2QvpeHzrbW1tXcl5So4XKiujISYgRZDEULV1dU0RKcaOL5mOseLTytIny0j5648fdLkkuuU2B62j745oIqvnNDQyi37IIOyQvZSGuExuvb92/+eyWxD7iT1H2KJEMrKyqI0NFXC8TXTOV58WkH6VIhxz550bFtWl3wOKE0VUC0tbUnpw7a9FjIo88jcSbytDr10t8XVhRpF7rFD7rfz+vVriuMDgPMgfSrE1NS0MJ/2ZQS8lnvzeHjglyvouLhpvwHPHt1v88shgzKJbLNtT+5ECL2sKEMI6ep3ojQ0AMA/IH0qxNnZmaaFh+rCMOxw9Jnfrv9KxzYs1rbjMxLatSE2ZFBmlN7JDJs7QktTM+JKQZtzJ0JIWlaq3ZHXpRufwtgAAHKQPhXCwODbf240Ucg3NFqz5BOC6vUK+pj1Qwi9rW7X3JiVW/Y5L/ZJ/fnHn1dPgxUV6FB6J/PQisn6nQ32pBTpd2lX5ivKzdbv1JmqwAAA9UD6VEj//v078nj3C+4xcK+zl69V/fXnpfhT1F52lJ0jQkjxpfuasth/x5dhkSU56ce/cINV/aiVcmjboRWT+wwa2rZxtvWU5N/pYdyb/Dk1NbXd0akqjs/c4HjxaQXpUyE4juvo6Ny5fZOBe/U1M58waYr/54uoXUWhew9jbR5+Pyez/Zca4ejy3ZGEp3dvRHiMhpXlKVFLED+vnnb9aLDj7MXb435rf+5ECBXcTJ8y8xPy52fPnrX/giqK42umc7z4tIL0qaiBAwZkZaQxc689RyK1tHWoHUOEYVjnLl0V3HelRYNHjQu9dFdXV2/PgjG5l9rVpQpePnscPK3f4zuZa8OjW7uuUFPIYbcjx4yj5GoAgIYgfSpq3rx5qSlXKO+SbBQPx78N2poUd/L+PSrHK5mZD7iXdZ2qqxkYCnYl5Y6c6Ho6YGnsBq/3b+hd8l5d5V6KCZll8ff7d8Hnb1G4nm3Zw/wOmGYfs/5UXRAAUA+kT0W5uroihMrLnjJzu2Wr/PiGRstmT6EwYdsLPy69f5fCIT8dMGzdntjVwUfzrsbvmmUBDbmt8v6N7OfV004HLB050fVQ1jPF91FRxMO7tzQ0NHg4Tv6zsJCuLYMA4CxIn4oyMTHh4Xh2VgZjdyTHEJ06+hNVFyRn0Ld25dsW2U2dt/vXwk4GXfYsGJP44zoYkauI++mXfpw5uOzerdXBR9ftiaWks7OugpvpQ4aOkP/zzdu31F4fAADpU1EYhvXs0YPutYfq6mtmPsV52rb1a6hah6hHbxOqRg/VY2AoCD53a5H/jptnjwRP61d6h/pbqA3ZS2nsBq/ItXO6dReEJ+fZTZ1H+S1qCeJu5lX5uCEAAB0gfbbC+PHjz8aeYPKOe45G6erpUzWGCMOw7kaC2ylJlFytng4Y5rLYJ0KU38mgy6EVk2M3eMG0lnpqCSLr1E8/OJvez7i0NCB0x7mb7ZzZ2RSygaHuuKHa2r/puJFK4PjMDY4Xn1aQPlth/vz5iJHFE+QwDDt4/NTtrIyE09GUXHDkGPvc9CuUXKpRBoaCXYl3lgaE3s+49IOzacqhbdCWSyq9kxk8rd+lsPUjJ7r+lFrs5L6cvnvdz8nU5uH9BlvKj7x/97akpIS+Oyozjs/c4HjxaQXpsxXGjh3bkcdjsvsTITR+onCE9Wj/zxe9eE5BZc7ZbV7Vq8rXL+itFzq5Lz+c9cxx9uJrR3bscDXLvRTD5ST6R1HengVjDq2YrK2lHXzh9ro9sTo8nNY73k5J0tXVg602AKAVpM9WwHHcwMCAye5P0vHTFzryeP9ZuaD9l7KyHqulw3sobsvGy63SAcNWbtm3MyHHuG+/+C3ewdP6cTCJvnz2OHaD154FY17/Ufbt4Qt7rxVRO7y2KRkJMS4z5zJwIwC4DNJn68z95BOGuz8RQjwcDw7fT0kTLg/HO3XufCctmZLAWiToY7b5RMqO8zeNTcw5lURL72TuWTAmZJZFYVrSIv8dBzOfWNpOYubW5IIJU+fW/7J19epVZgIAgCPUJ30SBOHn5xcREUHrXRYtWoQQYmD3lXrc5rpT1YRrP3HyxWN7KIlKQYI+ZkHRV3ecv9mzj1nc5pU7XM1SDm1Ty2UWagnifvql8PnWh1ZMfv1H2dKA0CPZEpfFPpTPS2lGw45PEmz5CQC11Cd9Llq0KCwsLCOD3o5JS0tLHo5npl2j9S6NOn76gq6e/n9WLmjnQgqzFy5D/6ujMImsie5MyLERTrv+884tE4xiN3ipzRQX2UtpyqFtQfYGJ/0X1L6p+vbwhSPZ5U7uy5lMnKSUM78IehhDx6ccx4eecrz4tFKH9EkQRFBQUHQ0NWNTm4dh2FArq8gjBxm4Vz08HD9x7tLtrIx2LqQweOgIbR5+OyWRqsBaRdDHbOWWfYd+K5v/ZdDjW9cPrZj83yl9Ug5tU9FZLu/fyHIvxexZMOYHZ9O046E2U9x2nL+591oRY0219bytlt3NvDrJZUa9431MzbOzqVnuWOVwfOgpx4tPKw1Vf3NLSkpmzJiRm5tL/tPd3T0qKqreORoaVBYzOjraw8OjSPJaviIak5bMnXk56XzSzQc9e5u0+SLTbS10+T2+O8JOBq2r+N7vF46EZl06S7x/281kwOjZy/rbfWzQow/bcbXg/RtZye20zBMRj7JTPsI0O3U1XLBuy5iPZzNf16yn8HbG5sXOojsldXfJtjLUdHab90pSSnfbDADKhqx805TmVDt9pqen29vbI4SMjIwsLCxEIhED6VMmk/ENDff/cnLSFMoW+FYcQRCDe3fT1NK5Ii5tcwNd4pkT/1m58Ojt53TPoFBQLUHcu3Ht4vG9uZkpNW9kHTt1GTN3pam1o1F/K62OShEh6Y+ivJKctFvxR/4oytPi6ep1MrCfNs91iS9Nqx+0wdEta29cjE3J+9e6jJA+AWfRmj5Vu4Pk8ePHCCF3d/fQ0NDk5GSRSMTATXEc7yEQ7I8IYSV9Yhh2KfWmwyiL4I3r/rN1V9su4jhlWgdMMz87dbjDx9SG1zYdMMzSdpKl7SQyj6bEHc86uffqga0Ioe7mQwZPmGFq7di1tzluwHSWqiWI5yWFJTlpBdcvPMpOQQhpdsS79+yzOviope1E5cmacheP7XHzWNLwuEFXft4tNelmBkBJqHb6tLS0LC4uNjExYfi+n332mb+/f7VMxkr7bV8zc791/ju3BtlNnDxO6NyGK/BwvEs3fsLRMCVJn3LyPIoQev1Cmnb+5O+plzKjIshUihCydJrTy2qMQQ8Tft+B+oY9KW8sffns8fvqqj8e5t1PS3p4M6X65XOEkGZHvJug54zlX435eHbv/hast9A2pbTwLkJo/qerGv4vy5Gj46IOMx4RAOpMST8IFGRpWX90PjOWLl0aEBCQmXaNlQooQmitf0B87Mlvvb3Opovr9nIpbvrcBYfCfnhbLVOS9tuG9LvwXRb7uCz2QQi9fiEtys0WZ6YU3ExLTk38u5aorXlPnmbpNAch1N/+n68RvSxt5FfADbrVbft9/0Yme/lc/k9pccHbqj8RQvfTkhBC8mSJENLU4enqdx403GaIjaP1pGldjajP03T4NfaIXieDQVbD2Q4EAE5QgQ8FStQbvf3hw4eGHaIymUxXV7fewYanEQRhaGhoZGTEVvst6eL1bEtTo/+sXLDnREIbOkEXrlzz894Q5Wm/bZ5+F/4IRxf5btJvq2UvpZKclKQ/K/+4dyO1tqbm3LbV5P+qeVut4DUxLR2Njz5CCOkbdOvSXTDcXmhmMdJk8NBugt6Gxso+dqmhWoJoquUWIaTfqXNtbS1BEBiGEQShqanZtoefldMaDVjx09C/u77aczVVfJcoLL6qvEuMzdXhSvpsqGFnMo7jDQ82PIJh2IcPH44cOfLpp5+y1X6LEOLh+JHouEVzZxwM2bbyqw2tfXmXbvxOBl2UsP1WETo8XNDHTLDYp97xWoKou5vpo7xbLyTP6p5gMcZRB9f730V4Sth52WZFuTdQEy23CCHT/oPev3v39OlTExMT8gGud4KCDz8rpzUaMCunwbukyGnK8y7RTQXSZ8OvEpGRkR4eHq26COXv7Ny5c719fC4mxM+a27pIKDR+otB94ZI9PwRNdJnZv8EqMy3y+3bzRr9lytx+21odMKxuxVEVK5Ftlpl0mqerBy23ANT9tKe1JqoOyyawAsfxPr17h/6wld0wNu8I7dPXdNnsKdWyVq+BN3n6J1o6vJQzv9ARGGDS22rZxWN73D/9vPnT2rleFQCgLhVIn8UNzJhRf1EVVoSHhz8ozGd+/du6MAw7e+n6+3dv1yye1drX8nB8kOXQ5Oj9dAQGmHTzynktHd6yL9Y3dQK5yEZCAtObBQGgxlQgfZo0gLPU3ViPo6OjgUGXpHNx7IbRXSDYczjydlbGvuDNrX2t99cbyx4WkBMegOqK27d9kIUVW93wAHCTCqRPpYVh2KxZbsFbA1lvE3Nydp09133PD0E5WemteuEYB6FeJ4NfY4/QFBhgQGnh3bKHBd7fBLAdiJLi+JrpHC8+rSB9tsu2bdu0tLTOnYlhOxC0PWyfef8Bfotmt3ZHs7mLl108tudttRpuH8YRZ/fv0OHhYxyEzZ9m2n9gamoqMyEpFZVel7T9OF58WkH6bBc+n29qasr6ACKEEIZhcRevvXv3prU7mi37Yr2WDu/CkRD6YgP0ef1CmpEQ47s+qMUzBwwZ+uzZsxZPAwAoSH3S54wZM4qLi0NDQxm+LzmA6EZm61pN6dCVzz+VcOXWb2nBG9cp/ioejo8dP/Hisb21bDdBgza4HL2/A6Y5y3Mp24EAwDnqkz5xHDcxMeHzmZ4ILxQK9fT1j+7fw/B9GzXCerTvV/+JOhiRKkpS/FXrgna8lf2VmRRLX2CADm+rZbHh30+f66ngoKH379/THRIA3KE+6ZNFm4OCzsaeePK4hO1AEEJorX/AwMEW3h7TFe8E7W1qPmTYyLh922kNDFDu5pXzHTBNv+8U6jsYP9n1Xn4+3SEBwB2QPing7e1tYNDlVJSyrD9wOf1Wp84GbuOGKt4JSs5gKbwN+0GqjFqCOLrlK4thI9u2ZwAAoJ0gfVJAPoOlDUv/0AHDsJQbYlnV61XzXRV8yRgHYVdDo/0bvGkNDFAoMyn2TdXrLREw6agFHJ+5wfHi0wrSJzVCQ0M1NTX3he1kO5B/dBcIdu05dDsrI/rgbgVfsm33EaiAqgp51bO3qXkrXlX7N30hKS2Oz9zgePFpBemTGjiOL1q0SHkqoAght7nus+e6b1u/RsG1FMY4CA0FPfdv8IYhuMov7qftra169htk8f7dW6m0ddOCAQBNgfRJGWWrgKLWr6VwOE5U8bQEhuAqOfmA21ZVPXFdPYSQTGm+3gGg6iB9UkZeAa1Umi/4GIZdvJ797t2bT2dOVGQYUW9T88FWw49u+QoWIVJm+zes0tTW+WYLrHQBAJsgfVIpNDTUwKDLhq/XsB3I/+PhuCgjp7T44bfeSxQ5f+ehk+/fvoFFiJRWxdPHGQkxX2zYAgvEA8AuSJ9UwnE8IGCj8swBJfU1M9+8IyQp7qQiw4gMjQSus+bGhn//+oWy1KFBXcGr5ujw8PmfrmI7EJXB8aGnHC8+rSB9Uszb2xvH8a98lrMdyL8sXrpy6ozZCg4j2hC8l6erH7Z2MQOBgVZJO3+i/HFR2M+xGIa19rUGXfkIoYKCAhriUmocH3rK8eLTCtInxTAMO3v2bGrKFWVYBbeuPUej+vQ1VWQYEYZhG4N33828mnPtIjOxAUW8rZZFrPMabDmsxc1VGkU29hYWFlIdFwAcBemTekKhsKex8VKP2azvA1oXhmEXrmQouCWLy6z5Q0eN3f3NZzCGSHmQI4Z+OgXfaQBQCpA+aZFy9epfr//85dA+tgP5F/mWLAdDtrV4cujPp99Uvd6/AfrYlELOtYswYggApQLpkxbm5ubTp0//7is/5ZnEQhphPXrF6i/2/BDUYidol278tZv+m5EQA024rHv9Qrr7m8/MB1ksWO7LdiwAgH9A+qRLdHS0gUGXVZ96sh1IfV9/F6TgWgoLlvuaD7LYvsINRuGyK2ztYuLdm+OJaWwHAgD4f5A+6YJh2O7dEakpV87ERLEdy79gGBZ38VpNzfv/rFzQ4snHE9N0OvJgFC6LLh7fezfz6qZdP0GzbdtwfOYGx4tPK0ifNHJ3d7e1td2wbo3yLIRL6srnH4qMvZ2VsS94c/Nn8nA87JfThbczY3dvYSY2UFdp4d2j33/pNHW2y6z5bMeiqjg+c4PjxacVpE96Xb58mah5r2zTQBFC4ycK3RcuUaQTdIyDcOrs+bHh38NmLAx7Wy373su5q6HR9v3H2Y4FAFAfpE964Tj+008/nY09oWxNuAihzTtCFewE3RC818R8QIDHJOgEZUwtQQR7f/L+TfVJ0Y02LJIAAKAbpE/akU24PksXKtsoXLITVJGZoBiGnUjO4unqfzV1BGxnxoy4n7bnZ6ftOBBlaCRgOxYAQCMgfTLh8uXLPBxf9amnUi2kgFozE5SH4+cy7hI1NduWTYcMSreLx/fGhn+//Iv/OEx2ZTsWAEDjIH0yAcfx+LNnszJSw3ZsZTuW+uQzQVNFSc2faWgk2L735/zstN3fLGUmNm4qvJ1x9Psv3TyWfL5uI9uxqAOODz3lePFpBemTIUKh8LPPPgveGpifJ2Y7lvrWb9o6cLCFt8f0FjtBHSa7bok4fPPKeRiIS5PC2xkBHpNMzAdsCN7LdixqguNDTzlefFpB+mROaGiomZnZ3KlOytYJihC68GuGrp6+Irtqu8ya7zprXmz495BBKUfmTmMT0zPXf6dpuFBtbS0dlwWAgyB9MgfDsNTU1DdvqpWwE5SH45dSs0uLHwZvXNfiyYEhB9w8lkAGpdbrF9IfVs4xNjE9l5FH31DbV69e0XRlALgG0iejBAIB2Qka6P8V27HU19fM3G+df9TBiBY7QRFkUKpVPH38xcdDP/xN0Jc7yZb5ly9f0nFxADgI0ifThELh+vXrD+0LV8KZoGv9AxTsBEWQQalTeDvDVzhQR0c78UYhffXON9UyhFB6unJtQwuA6oL0yYKAgAAHBwefpQuVbUtt9L9OUEX2BEX/y6DnD+46umUtzGZpG3l/Z/LvxV268em+3f0HD+i+BQAcAemTHSKRyMzMzGu+m7INI+Lh+IlzlxTcExQhFBhy4KvA7cnRB2A+aBuknT9B5k5a+zvrev/unVTJHjm6cXzmBseLTytIn+zAMCwzM/Pdu7eOoy2VbUF5+UzQ+/cUmmMzb8mKkKOn8rPTvp4xClb1U1AtQcTu3hKxzstp6mxmcmdFeRmmqVlTUyMSiei+l1Lh+MwNjhefVpA+WcPn84sePKiq+strvpuyDcT9+rugPn1Nl82eomBqd5jsGnUxvfJZ6Zcuw2Fl+RbVEsS2ZdPj9m7/fN2GnYdPMFPvfPa0VFNT06BL1xMnTjBwOwDUHqRPNgkEggvnz2emX/d0c1GqDIph2NlL12VVrwO/XKHgSwZZDU+596xz504BHpMuHocp/02qePr48/F9C29nhhw9xfy6QktX+qSmwbbbAFAA0ifLhEJh3JkzmenXlW09v+4Cwa49h5LiTiacjlbwJTwcP5eR5+axJHLH+u+9XKAht6GLx/f6Cgd+hFBS9n2G17MV37qB6+raO06SVVVxrfsTADpA+mTf1KlTgwIDd+/a8eO2ILZj+Re3ue6jbGz9P1+kyDwWEoZhgSEHNoceeCi++aXL8JxrF2mNUIW8rZZ97+USuWO909TZv959wvw+Ki8rpX1MTC2shr1//55r3Z8A0AHSp1Lw9/dfvdpnT2iwsmXQ6PiLis9jkXOZNf+q+ImRQLB9hVvY2sVvq5VrbBTz0s6fWDKi20Pxzc2hBxjr7Kznj/IyPT19Ho73GzAoPDyc+QDYwvGhpxwvPq0gfSqL7du3L1+2LHhroFJlUHIey+2sjFNHf2rtC89cy/n6+51Zl+J8Jg5IO8/R4SqvX0i/93KJWOc11HrMVfETl1nz2Yrk9m9pEyZ/jBBavGzl7ZwctsJgHseHnnK8+LSC9KlEdu3atXTpUmXLoCOsR0/+2HXb+jVlpSWtfe2C5b7pD6Rm/fr/9N3nX08fVfH0MQ0BKilyaspy296PxLf+u+/YscRUHo6zFQw5glpPTx8hNNbe4d3btzlcyqAA0AHSp3I5ePDg0qVLla0Vd8/RqG6Ght4e09owPJiH48cSru8+Hvf82WNf4cCwtYu5kETTzp9YYW8St3e7m8eSX8WlLFY6SS8rpQghu/ETEEKDhlgihH755Rd2QwJA1UH6VDoHDx708fYOC962yktZNmbBMCz+cuqTkkeK7MfSqDEOwrT7FV9/vzM3PdlXODB29xZ17RAtvJ2x1nV4xDovs37949NyA0MOsFjplCPXTODx/olk6crVMadOsRsSAKoO0qcy2r59++nTpy/En1ae+aB9zcxXrP4i6mCEgksRNYRh2ILlvtfvPZvntSJu7/YlI7rF7t6iTjVRMnEGeEx6+9er/aeSjiVc721qznZQ/yDXTOjK/2dNXUfh5OdSKUxfAaA9NLjQsayhoZLFvHDhgtusWWPtxkfGJbIyVrMegiAm2lhVVr64Ii5tZzzVMtn2b9dcOH2i5t1bW9e5M5ev6z3Agqo4GVZLEJlJsXH7tpc9LOhqaLRt95ExDkK2g/oXK0NNZ7d598U5qTn55JFqmczcSD8kJMTPz4/d2ACgFTnwmKbPf5XMK62loukTISQSiaZOm6arq5dyQyyvOrDoj/JyG0uzSS4zt/90vP1XIwjixOE9YVs3vq2W9TQb6LbyG+tJ03R47Dd1Kqji6ePr8cdjw7/vgGlaDLf2/nqjsiVOkpWhZmeDrhOdpkQcOiY/OG/aZKmkLD8/n8XAmKG6f/6UgOIj2tInNN4qNaFQWPzo0Zs31Y6jLfPz2thqSqHuAsG3QduS4k4qsqV2i8jm3Iyi5/tPJWlqfNjrv3zJiG5Ht6xV8lVz31bL0s6fWOs63Fc4MH5/sJvHkiu5j48lXFfO3El69bJy6EjrukfmLVxSXFwsU7LtCujA5eSBOF98WnHii4mqf/+SSqVjx459+PDh2cvXR4+1YzscNNFmqKS8/OLth9QOinnxXHpsX0j04b3VVX8hhD5euGqs82xzq9EdlKDhGiH0+oX0ofhWwtGwu5lXtXR4PYx7rfzqu8nT5yhDu3rzrAw1EUJZdx/26mMiP1gplVqaGkVFRbm7u7MWGQA0g8bb9lL19IkQIghixIgRYrH4q/UBX/ozvc54PdUy2VDznpYjRu2PvUTH9fNzc04c3iNKiP/rz5cIIVvXuSMcnfsPH2to3IeO2zWjliDKHhbcEMVnJsaWPSzQ0uHp6ev7fbt58vRPlGE8rYKsDDW1tHWKJH/Wy/QO1hb6urybN2+yFRgAdIP02V5qkD5Jn3322fHISOepM8MO/MxupScuJtp76YJte39xnU1j3aX0UVFc1JFL52IrJJL3b6tRnVTa1agnTbXS1y+k5SUPxJlXb1w6W3r/LkII09SymzDZedY8xynTVChrylkZavbpa5qZ+6De8cM/RWwNWP/nq1fKX4EGoG0gfbaX2qRPhFB0dPTSzz4z7tXnTNJVdgcTzfp4wm/p11PuPevSjfYwqmWy3JuZSXEn0369/OpFZc37dwihnmYDLW0nmluN6t3fQgfXa0PdtJYgKiVlb2V/ld6/W5Sb/bQo/27mVYQQpqmlpa3Tb+DgKTM/mfDx9J69TSgvEZOsDDXdF3rt3HOw3nGy/TY5OVkoVN5eWwDaA9Jne6lT+kQIFRUVjRg58q/Xr9ntCiUIYnBvvqBXn5OiG0xWXwiC+OPZ06sXz2Vdv3ov9/brP1+9e1NN/i/dzl2t7CYhhHqY9jfqbdbwtUW52a9fPEcIPS7ILXtYID+uqa3T2aBLD+PeU2Z+Yj5gcP8hQxn4TsCMstISZ+t+J+IvjZ/YSI50sLYYNKBffHw884ExRs3+/FsLio8gfbaH+j1AMpnM2tq6oKDgq/UBvuvWs9X4dvvmDbcpDl9t+sH9M29WAiCVlZa8fvXy1m+pD+6JH90vQAi9/vPVs6elDc/sY9qvY8eOCCHjPn2HDLfW1dMfZeeIEHQvSJUAACAASURBVFL1+mUzEk5H+3++qEjyutFm5x+3Be0JDYb2W6CuIH22l/qlT9L69et37txpYztuz+FIthpy/b/w/vngvtiU2/0HW7ISAGjeNysWJMWdfPZXbaP/Nz9PPGnMMGi/BeoK0md7qWv6RAhlZWVNnDRJW1snbP/RSVNcmA+AIIiRA/oQtbXtX4oIUI4giBE9OiKEmkqfCKEhfQzH2dupd/st4CxYNgE0ycbGpuKPP4x79lg4Z9qGdWuYXyAXw7DLaTdlVa+/9V7C8K1Bi8S3slo8Z+lKn8SkJCVZWhkAFQLpU+XhOJ6bmxsSEnL8yP7h/YyZX5you0AQuP1HqpYiAhTKSv21U2eD5s9xnu5G1NSkpKQwEhEA6kNtWzXrUuPG27qKiors7O0r/viDlfFEM4TjsrMymJnHAhQ0foCR89Tp0ceONNN4i6D9FqgvaLwFCjE3Ny97+tTf3z8seNukMcOePC5h8u7R8Rd19fSXzZkCzYBKIicr/dXLSq8VLQ+KJttvGQiJFeQHKGdxvPi0gvSpVjAM27p1a15e3p+vXthYmDHZG8rD8TNJv5YUFR4M2cbMHUHzLsXHdupsYDF0eItnku23WVktd5SqIi60PDWD48WnFaRPNSSvhh47/BOTvaEWQ4evWP3Fnh+CcrLSmbkjaApBEFEHI1ymzVTk5EFDLA26dD1w4ADdUQGgTjjRKciRvs+GioqKxo0bJ5FIZs6ZHxyxn5n1WsePHPxc+vxsuhg6QVmUKkryWzw7p/BJVz6/h16H5vs+EUI/bgs6tC/8RWUlM+EBwAzo+wRtZG5uXl5eHhISkhB/etTgvmdiohi46cXr2e/evfnPygXQCcqinZu+NulrqvhiGs7T3V6+eJGTk0NrVACoE0if6s/Pz+/ly5eDBg5Y672MgSFFPByPv3z9ZmZq8MZ1tN4INKWstOTR/YLvd4Qp/pJ+AwYhhHbv3k1bUACoG0if6k9DQwPH8fT09Avnz5eWPCKHFFXLZPTd0WLocL91/lEHI2AmKCvCtnzH4+G24x0VfwmGYUtXro45dYq2oFjD8aGnHC8+rVQ+fcpksvT0dD8/Pw0NDT8/v/T0dGgzrKempob8QSgUvnr1atOmTceP7Cfbcul7r9b6B4yysfX2mF5WWkLTLUCjXjyXJsWd/CZgc2sn/k6b9cm7t2+lUilNgbFF/vxzE8eLTyvVHlMTFRXl6enZ8HhCQoKLy/8vAMvZoUNNkclkTk5Ot3NyevfpuyP8J5p2PSMIYnDvbppaOkk3H6jiLtMqal/w5v27thU8qZS/54oMHUIIEQTR20Db399/69atNMcIAENg6FDjCIJYu3YtQsjd3T03N/fDhw9paWnu7u4IoaVLl7IdnVLDcTwjI+OuWPzXny9nTh6/ysuTjg5RDMNSb+XLqv5as3gWNAkwo1om2/ND0CfzF7Th+wqGYeMcJ8XExNARGADqR4XT5+XLlyUSiZGRUVRUlKWlJULIzs7ul19+QQhJJJL0dJh62AJyXG5UVNSVy4k0dYh2FwgOHIu59VsarKXAjLAt32Gampt3hLbt5ct91jwrL5fR2S8OgNpQ4fTp7++PEPr888/rHsQwzNfXF8EYQoW5u7u/qKwMCQn55dA+cyP9H7cFUZtEnZxdfb/6z54fghJOR1N4WdBQWWlJ1MGIDd9vb3NT+bARo95UV587d47awABQS6rdKVhSUmJsbFxviMTQoUNzc3N9fX1DQ//5Dg59n4qQyWR+fn5R0dE6Oh037wiZPmsuhYvOL18w90L86Z/Ppwy3oaWfFSCEvlmxIFWUdK9UWu8Xp2DfJ8nB2sLUpPfFixdpCBAApkHfZ5NMTEzqfVKUlJTk5uYihKZMmcJSUEpHwZHrOI4fPHhQWlExcED/NSs/HWrWg8KhuXuORvXpa7p4muOL5+o2sFNJ5GSlJ8Wd3B6yp51fembMnpty7RpVUSkDjs/c4HjxaaVW1TKCIJydnUUikZWV1a1bt+SfI1D7bK2ioqLFixffun1bV1cvcPuPlNRECYIYZt7zA9KIS70D6/lRiyCIOY7D/3r16s7Dsob/t1W1z/w88aQxw27fvj18eMtrzQOg5KD2qRB57kQIxcfH1/u41/g31NiXMplM1vBgwyMEQTB/WnsCbsNp5ubm6enpd8XiF5XP13ovG97PuP01UQzDrmT+/vaNzG3cUBiIS63gjetKix+evUxBrZFcfujUqVOI2UeumdNY+Ytr9DRWig/vUmtPq/dRTx8VqJZFRdVfqdXW1tbExKTuEZlMNnPmTDJ3pqWl2dn9q4MNap/tQW1N9I/ychtL8xE2tntOJDC8obe6un9PPMdxhM/ab9Zvany+ZqtqnwihVV6eBXd/z8/PpyhAAFhDa+1TBfJKw28QkZGRHh4e8n9KpVIrKyuJRIIay50I0icVKEyimWnX5k2fYj12HGTQ9iMIYpJlbw304feisqbezNamzzMxUWtXfVZVVQW/HaDqoPG2OfLcaWRkVFxc3DB3AkrIm3MH9O8nb85t26XG2jscijx9MzN11XxXaMVtp+CN62RVr89fSacw1Y2ysX337t2tW7eouiAAakkF0ueHBuRVz5KSEjJ3WllZ5ebm1mvRBSQKOwDqJVEHa4sbmW1ZnsLJ2RUyaPulipKiDkZ8G7Str5k5hZft1ccEIRQfH0/hNVnE8aGnHC8+rVS4VVMmk5mbm0skEqFQmJSU1My3b4433hIEQUcrnEgkWrhwoUQi6TdgUNsWzk1OSljqORtacdvmxXPpTDtLI4Hg16w7zZ/Z2sZbpF7dnzQ9/6qC48Xnet9nUzw8PKKjoxFCaWlpPXv2rPd/cRzn/2+vYI6nT1rJk+igIZbbdu1ubRKFDNo2ZJfn2zcy8SNJi2sMtSF9nomJ+tp35Z+vXsEvBag0SJ+NkEqlhoaGzZzg7u4uH7IL6ZNuIpFo8ZIlz8rKBg2xjDh0bNAQS8VfCxm0tQiCWDXf9XZWxtWsO4o027YhfT55XGJjYfbgwQNzcyqbhQFgGAwdasT9+/fZDgH8P6FQWPb0aXJy8p+vXkwaM6xVW7hAP2hrHQzZdjMz9cCxGGq7POsS9DRGCCUkJNB0fQDUACeqZVD7ZAxBELt37/72229lMtnMOfP9N20hx6G0KDPtmvtMF1xXD9Ykat6+4M17fghau37jWv8ABV/ShtongsVvgVqA2idQGRiG+fn5vXr1KiQkhNwHTcGa6Fh7hyxxEbkm0f17YvojVUlk7nRf6KV47myz8ROEGZmZdN8FANUF6VP9MT9ynUyi5D5o165ctrEwU2QftO4CwY27j7AOHeY4jkgVJTETqgqR586dew4ycDtH4eT3796pQXM6x2ducLz4tIL0qf7Yargmk2hFxR/+/v6hO7YqsploVz7/VuHjgYMtvD2m7wvezFioyo/h3IkQ6j9w8Lt3754+fcrM7ejD8Y4bjhefVpA+Ab0wDNu6devLly+XLl26O2THqMF9m199HsOwX7PuLP5s5cHQH5bPmULt3t0qivnciRDq2o2PELp79y5jdwRAtUD6BEwgNxN9WFQ0aOCAtd7LJo0Z1vxyRdt27T5wLOZ2VoazdT8ud4USBLF8zpQDIf9dtsqPydyJECKnkz58+JDJmwKgQiB9AuYIBAJyzb+//nw5c/L4VV6eldImd892cnbNEheRXaHRB3erQSdca714Ll013/VmZmrI3sOB239kPgAb23FJSdAJDUDjIH0Cppmbm5eXl5OjiixNjQ7/FNFUauwuENwqfDx1xuzgTV+vmu/64nmTuVb93L8nnmlneefmb8djL7jNdWclBkGPnq9fv2bl1gAoP0if6k85h975+fk9eVI6Z86coPXrJo0Zlp/XeAsthmH7j8fEXbp2N+em4+AeCaejuVANjT64e47jiE6dO2X8fn/8RCFbYViPGXs3L4+tu1NFOZ9/xnC8+LSC9Kn+ampq2A6hcTiOnzp1Ki8vj1yraMO6NU2lxhHWo++VSic7TwtYs3yecHRZaQmzkTLnxXPp8jlTtq1fM3XG7NTb+d0FAhaD6WzQRQ2+rCjt888MjhefVpA+1Z+SLyRrbm5e9vTppk2bjh/ZP7yfcTPV0KMxZ89cvFpW8sjZut9/13+hfoNyE05HOw7ucefmb7sPHd9/PIb1X1ynTp3/rq1V9QzK+tvILo4Xn1aQPgH7MAwLCAgoLi7md+vaYjW0oOzF2vUbY48d+HiEmdq05ZaVliyfM8X/80UjR48RP5Kw1dlZj9pM/QSADpA+gbIQCARisXjTpk1HD+wd3s+4qaX+MAxb6x9w77HUrF+/jX6fTbLsrdJLFFXLZPuCNztb98u9mbX70PHzV9Jb3IAMAKAMIH0CJUJWQwsK8jt0+MjGwuzwTxFNncnD8fNX0lNuiM3N+3l7THcYpHqjigiCSDgd/fEIs4Oh233WflNQ9kJJKp0AAEVA+gRKh+wNnT59+ubvvpk3bXIzfZx9zczjRamX026am/cLWLN8kmVvlUiiZOKcOMR4g+9Ss379buU/Xr9pK/RRAaBaIH2qP1UcuY5hWHx8/JHDh3+/dcPGwqyp8UQki6HD40WpV7PumJv32+C7dESPjvuCNyvn6NwXz6X7gjeP6NGRTJzXsu+ev5LelQ8btNFIFZ9/CnG8+LSC9Kn+VHfNaHd398LCQk2sw6Qxw5ppyCWRNdGcwifuC70Ohv7gbN1vtuOIVFGSMgzQJQgiJyt9+ZwpjoN7HAzd7r7QK6fwyfkr6fTtd00hVR86pLrPPyU4XnxacWIfadguW6URBOHu7h4bG7t05eqAbcGKNHISBHH+zKld2zeXFD8iamrGjJ+4cOWakWPHMzwqhyAI8a2sS/GxUQcjEEJGgh6fr/nKa7k3w+20bdsuGyH05HGJjYVZcXGxiYkJ1UEBwARat8vmRF6B9KkG9u3bt9rXd6zd+D2HIxVv7ayUSn/avev8mVN/SCRv31T3G2w5Z8HSkbbj+w+2pC/UstKS37Mzr19OSIo7iRDS1dObNnOO1wpvi6HD6btpMyB9As6C9NlekD7VQ05Ojv24cR078hJTfuvVx6RVr62WyS4lnDv8U4T4zu/v371FCI0ZP3HEGPshw6xN+w8y6MpvT8X0xXPp44f3nz0tladMnY68zp07T5v1ySfuC9nKmnKQPgFnQfpsL0ifaqO8vNx61KhnZWVnL18fPdauDVcgCKIgT5yZfv3c6ZjC/HsEQbx9U03+L2e3eQihYaPG6ht0QQj1G2SB6+rVfa2s6q8H+f/sfym+deNlpbQw786j+wXkkY48Xk/jXqNsbGfMmW9tM1Z5pm9C+gScBemzvTiePtWs+ARBWFpaFhQUfB8c+ukKn3ZerVomq3wuTb9+9dHDB5mp1xBCRfcL379/hxB6U13d8PyOPB75A66r28fEtGtXvp3jhAEDhwyxHKq0A2jbnD6vXEr8zHNOVVWVSk+qUbPnv7Wg+Ii29KnCfxVAQWq2ZjSGYWKxeOXKld995ffqxYsv/Te252o8HOfh+PyFXg3/F0EQ5WX/GnTatVu72nhVzp9/vvqoQweVzp1I7Z7/1uJ48Wml2n8YQBGq/vHXEIZhBw8e1NbWDv9xe1ZGWmRcIh1lxDCstT2saqbkYZGenl7L5yk39Xv+W4XjxacVzPsEqmr37t1HDh/OTL/u6eai/CsNqaKi+4VmpqZsRwGAkoL0CVSYu7t73JkzN29kerq5KMPyCGrm2q/JgwcPZjsKAJQUpE+g2qZOnZqUmHgr+zcbC7NKqZTtcNQHQRAvX1RaWtI4QRYAlQbpE6g8BweHogcP3rypdhxtCRmUKuSwqRkzZrAdCABKCtKn+uPCmtECgQAyKLWyszK0tbWNjY3ZDqS9uPD8N4PjxacVpE/1x5FZX2QG/fvvWsiglBAlJRgYGKjBuE2OPP9N4XjxaQXpE6gPgUBQWFCgoYHsRwx68riE7XBU29nYE66urmxHAYDygvQJ1Aqfz7+Xl9etaxcbC7Mbmelsh6OqyC8f8+fPZzsQAJQXpE+gbvh8fkFBgZmZ2czJ4yGDtk3yxQva2tpjx45lOxAAlBcnlkPk+KqP3EQQxMCBAx8+fNjmxeXVRhvWvHWwtvgI/V1QUEBTSAAwg9Y1b6H2qf64OfQOwzCog7ZNpVT6oDA/ICCA7UCowc3nX47jxacVpE/1x9k1oyGDtk38mZPa2trTp09nOxBqcPb5J3G8+LTiRKsmNN5yGbTitrbx1sHaoluXzhkZGfSFBAAzoPEWgLaDOmir3MhMf1CYHxQUxHYgACg7SJ9A/UEGVdzR/Xt69OwpFArZDgQAZQfpE3ACZFBFPHlccjb2xNfr1rEdCAAqgBOdgtD3CUjc7AdVvO9zlZfnlcuJLyor1WCtPgAQ9H2CdoKR6/J3AOqgzSCrnpuDgtQsd3L8+ed48WnFiWoZ1D5BXVyrgypY+5w3bfLt7N9evXqlZukTcBnUPgGgElkHHT58ONRB5a5cSkxNuXL27FnInQAoiBPVMqh9goYIgrC1tc3Ozo44dGzWXA+2w6FRi7VPgiCGmvXobsjPz89nLCoAGAC1TwCoh2FYRkbGnDlzfJYu/HEbp6c5Bvp/VS2r+vXXX9kOBABVAukTcBeGYadOnVq1alX4j9vnTZtMEATbEbHgRmb6oX3hX375pUAgYDsWAFQJJ1o1Od54y/HiIwXegejo6EWLF4+1G7/ncGRXPp+xwJjRTONtpVQ63nqIejfbcvz5h+Ij2hpvOfHOcvwBIgiC4+NBFHkHcnJy7MeN69iRl5jyW68+JozExZCm0idBEJ5uLreyfyspLuar3ZcGOY4//xwvPvR9gnbh8h8PSZF3YPjw4UUPHmhra9lYmJ2JiWIgKtZ5urlkpl+PP3tWjXMn4vzzz/Hi0wrSJwD/EAgEj0tKHBwcfJYu3LBujXp3hf64LSgz/XrcmTOwvC0AbQPpE4D/h2FYSkrKpk2bfj60z9PNpVIqZTsiWvy4LSh4a6CPt/fUqVPZjgUAVcWJTkGO932CNsjKynKaPPmv16+PxZ6fNMWF7XDapV7fJ5k7ly5devDgQRajAoAB0PcJANNsbGzKnz2ztbVdOGfaKi9PtWnI/XFbUMiOrV9//TXkTgDaCdKn+oM1o9v2DuA4np6eHhISknQ+bng/Y1Vf3o8giHnTJofs2BoUGLh9+3a2w2EOx59/jhefVpxo1YTGW9Ae5eXlEydOLCgomDlnfnDEfh6Osx1R6/TQ6yB+JJnlPKG05FFsbCz0dwLugMZbANgkEAjy8/NDQkKuXEowN9I/ExOlcm25lqZGkvKyu3fvQu4EgCqQPgFQiJ+fX3l5+dixY9d6L5s0ZphKtOUSBLFh3RqE0MCBA8ufPTM3N2c7IgDUBydaNaHxFlCoqKho3LhxEolknOOk4Ij9SrtE0ZVLib7Ll1TLqt69ewfPP+AmaLwFQImYm5uXl5dHRUX9fjvbxsJslZdnfp6Y7aD+5cnjklVengvnTOtuyC8uLmY7HADUE6RPANrC3d39RWVlSEhIctL5SWOGKUlzLpk4bSzMrl25vHfv3vz8fNhHBQCacKJVk+ONtxwvPqL5HSAIIiUlZfGSJc/Kyrp243stX7X4s8+Z37YlP08cHvzfs7EnDAy6BARs9Pb2li92yvEHAIrP8eIj2HGleTKZTCqVGhsbN7o+MscfIMCMnJycTZs2Xbp06d27d+McJy33WTNsxCi682i1THYxIT70h60PCvP19PU3BwXVTZwkeP4BZ0H6bBJBEJcvX961a5dIJCKPCIXCqKioejtIwMcHYIxMJouJidmyZcuz8vI31dWDhlh6en3m9PFUQc/Gv9u1TaVUeu1q8sljR1NTrmCamuZmZuHh4Y6OjvD1EYC6IH02KSIiYvXq1QghIyOjCRMmREdHkz9nZmaamJjIT4OPD8A8mUx27ty5nTt35hcUVMtkCKFxjpNsbO3tHScJevRsbTatlskelzzKzxPnZN84G3ui8rm0I4/Xp3dvNze3b7/9Fm92JQd4/gFnQfpsnEwm09XVRQhFRkZ6eHiQR2bOnCkSiYRCYXJysvxM+PgALCIIoqSkJCEh4eTJkwUFhTJZ1fv378n/NXPOfIRQ12784aNGN3xhTvaNyufSyufS1JQr5BEejnfu3Nl65MjVq1ePHTu2+awpB88/4CxIn40LCgoKCAiwsrK6c+eO/GB6erq9vT1CqKqqSv7hAh8fQHmQ/fRXr1598OBBSkoKQuj58+dPy8oanmncs2e3bt0QQtOnT+/evfuECRPqtqkoDp5/wFm0pk8V3ojczc3N3Nzc1taW7UCUHXx6KtU7gOM4juNeXl5sB8IVSvXbZx7Hi08rdXtnhw4dmpubW69KCg8Q4DJ4/gFnwapDLSspKYmIiCBzp5GR0fHjx9mOCAAAgDpTk6+ldfe0k48kqvt/1aOYALQBPP+As6D22TJfX9/w8HAjIyOEkKenp4eHR70tpTT+DTW2i6xMJmt4sOERgiCYP609AbNyGivvUqOnwbukyGnwLilyGrxLipzG+rtU76OePirwtbReVRIh5O3tbWdn1+jJiYmJrq6uCKGEhAQXFxfyIHz7BlwGzz/gLK5PXGn4DaJh82xdZA9o3amf8PEBuAyef8BZXE+fUVFR9Y7Y2tqamJg0tc5tVFSUp6cnqvOWcfzjg+PFR5x/B6D4UHy2o2AN1+d9NlXRNDc3l0gkgYGBGzdurHv8yJEjCCFfX18mglMFXP7jIcE7wGUc/+1zvPi0UuEvJvJaZk1NjbwCKpVKDQ0NEULFxcXyJVo4/v0LcBw8/4CzYORt4+bOnUsOtXV2dk5PTy8pKQkKCrKyskIIWVlZtW15MwAAAEARqv21ND09fc6cORKJpO5BX1/fL774AnZcAYAEzz/gLK4PHWqRWCyOi4tLTU394osvHBwcGm5DAR8fgMvg+QecBemzvTj+8cHx4iPOvwNQfCg+21GwBvo+Qbtw+Y+HBO8Al3H8t8/x4tMK0qf6Y2DxKmXG8eIDjv/24fmnD6RPAAAAoNUgfQIAAACtBukTAAAAaDVInwAAAECrcWJMM/ScAwAAZ8G8TwAAAEBZQOMtAAAA0GqQPgEAAIBWg/QJAAAAtBqkTwAAAKDVIH0CAAAArQbpEwAAAGi1Dps2bWI7BnoRBJGZmZmWloYQ0tXV1dLSYjsipkml0vj4+CdPnvTt2/ejjzj6hamkpCQxMdHU1JRrD4BMJhOJRNnZ2Xp6erq6utx8AGQyWWxsrJ6eXufOndmOhWlisfjKlSsIoe7du7MdCwvI4r9+/bpjx44Nt4Jurw/qq6qqytfXt155ExIS2I6LITU1NeHh4VZWVnWL7+vrm5uby3ZoTKuoqDAyMkIIFRcXsx0LcyoqKoRCYd3fvpGRUUVFBdtxMa2mpoZ8HyIjI9mOhTnFxcWBgYHkYy//7QcGBnLnAahXfIRQYGAgtbdQ5/QZGBgozxmRkZHyj5Lw8HC2Q2NCQkICWV4rK6vw8PC6D1NVVRXb0TGnqqpK/h2CU+lT/sD7+vrKv0caGRlx6rdfU1Mj/xzgVPqU//bd3d3rfvq5u7uzHRoT5J9+QqEwMjJS/vwLhUIK76K26bO4uLjhJyb5h2RkZMReXMwhc0bdj4yamhryIEe+QHz48CEhIaHuN1DupM/IyEjym5O8yMXFxeRvn/Lv4EpLXmSupU/5p1/duqb8oNp/f6qpqSH/6uv+xuUJlcLiq236JL9uWFlZ1T2Ym5vb8KlSS2Rfb8NnhXyGOPIFwt3dXV7/5lr6bLShRf7br6mpYSswxpBfIMjyNvwwVW/kk9+wosmR70/yTFn3Oa+pqSEPUth/p7bjCEJDQ4uLizMyMuoeFIvF5A88Ho+NoBgVGRkZGBhYr7e8U6dOCCGJRCKTyViKiznR0dEIofDw8Fu3brEdC6OkUin5w7x58+oed3BwQAhJJJKnT5+yEBazLly4gBByd3cvKiqaMGEC2+EwaurUqeHh4f7+/vWODxkyBCFUUFDARlDMcXFxqaioqKiowDBMflD+zJOfgZTAWj5FZZmYmNT9p0wmW7t2LUJIKBRSPwRLydjZ2dnZ2TU8vnv3boSQkZGR2r8DCKHAwMC1a9dyoaT1SCQS8gc+n1/3uPytyMjIqPfXoX6mTp26detWtS9mozw8PBoelMlk5BdKW1tbxiNiWr0nnyCI9evXkz/b2NhQdRe1rX3WlZiYOHToUF1dXYlEIhQKo6Ki2I6IHVKplPz7+fbbb9mOhQkbN27kYO5EdVpZGqo3EluNeXh4cDN3NmXnzp3kD1OnTmU3EiaJxWInJydNTc3o6GgjI6Pi4uK6VdJ24kT6dHV1lfd6Dh482MDAgN14WCGVSsmPTisrq5UrV7IdDqBdvVH7JLL5DnBNYmJiQEAAQig8PJxT3yq2bdsmEonIny0sLOrVStuJE+kzISEhNzeXHHYbFhbWq1cveecQR5C5UyKRGBkZiUQiCr9/AQCUXFRUlKurK0JIKBRy7auzt7d3cXFxeHg4+dGnq6ubnp5O1cVVPn1qNNCwbdbFxcXS0nLjxo1VVVUIIYlEcvLkSTaCpZ6Hh0e94jfs9khPT5fnztzcXGq/f7ErKiqq4QPAdlDss7S0RHV6QOvKy8tjPBzApqCgIE9PT4SQUChMSkri2ldnOzs7ExMTHx+f8vJysvmNwoX2VD59tgqO4+SElgMHDrAdC0PS09Pt7e0lEomVlZWa5U7QFD09vab+l7wXA6g9giCCgoLINtvAwMDk5GSu5c56vvnmG4SQSCSiat6Byr+b8rnAcmSGSExMfPXqlZOTU72EQQ67UpsPkdDQ0K1bt9Y9UnewDJk7kfp+8ZwxY0bDBwDInwGZTFb3eSAIgvyBC2MvOY4gCGdnZ7LbLzIystGxuOpKLBaLxWJLS0uyGUZO/thXIPomIAAACIpJREFUV1dTM6iQqgmkyobjE4c/1Fk5ITAwkAvT5JtHvhXcWTah4ZpTH5qYTs4F5KcBd5ZN+FBn0b60tDS2Y2EauWJGw+VByOEv9dbSaQ+1TZ+NLtrHnVWHqqqqyIGXvr6+bMeiFLiWPptZtI87SzbKcS19yldc4mDu/FBngaG6v3Fy4Au174nGh/99sqgfDw8Pcpqjr6+vjY3N9u3byfTp7u6u9lM//fz8wsLCEEJCobDR/s7Q0FBO9YOSQ4qKi4s5MmqfIIhevXqRo4d8fX27du1K9oFZWVndunVL/Zrxm0d+FHCkDVMmk5mbm5NDBRtdbsnW1tbHx4f5wJgk7/QVCoVeXl4XLlwgc4GRkdGTJ08oe/6pysNKqKKiot6eNVZWVgkJCWq/YvIHBb4ScaceRuJgqSsqKuotklC3MsopnKp9tjiwgwubrjTcrtHIyCg8PJzadkd1rn3KpaenP378uOEwIjVWUlLS/AnGxsacqoKQbwjXSk0QRH5+PrkIEaee/3qkUqlMJuPz+VxYhUomkzU/rx3Hce48CU0NI6IEJ9InAAAAQC1uzfsEAAAAKAHpEwAAAGg1SJ8AAABAq0H6BAAAAFoN0icAAADQapA+AQAAgFaD9AkAAAC0GqRPAAAAoNUgfQIAAACtBukTAAAAaDVInwAAAECrQfoEAAAAWg3SJwDqjCCIZv7ZFLFYPHToUA8PjxbPV/CCAKgf2HEFAPbV27/d0tJST0+v/RtsyWSymTNn/ve///3ll1/u3bsnEonCw8Nb3Co5PT3d3t6e/FkoFCYnJzdzspOTU1JSEqe2gQOABLVPAFhGEETnzp3Xrl3r6enp6el54cKFyZMn9+3bV1dXNygoqM2XLSkpMTc337Rp08CBA7/77ruKigqEkIODQ4svvHLliru7e0VFRVpamkgkIrcLbcqMGTNiYmLaHCQAqgtqnwAohaCgoICAAHltT14FLC4uNjExae3VyJenpaXZ2dkhhAiC0NTURAgp8vfu5OR09uxZsuLr4eExcODAjRs3NnWyVCq1srIqLy9vbYQAqDqofQKgFAoKChBCM2bMIP9pZ2dnZWWFEFq/fn1rLyUWi+vmToRQfn4+QsjX11eRlycnJ8sbjRcsWLB3795mTubz+eQdWxskAKoO0icASiE6OhohNHXqVPmRIUOGIITy8vJadR2pVDp58mRfX1957kQIxcXFIYTmzp3b2qgGDx4skUikUmkz50yYMIG8PgCcAukTAPaRtTcjo/9r7/5dkvviOIB/hjabguA4SBA2NDyXtqCpQIdaHIJCB8FBKCodHBpLB6NJqrFBkBQaGoRuQ9QgaBAFgf4BBjpckqY8kwd6hg/P4X799c3v9wlJ36/tns69XSF6cz5+zr1C12mllByo+/v7XU9RSr2+vnYOBgIBy7KSyaR9/OrqiogWFhb0NCnlV24sk8kQ0dPTU/9pvHQGGCuIT4Dh49Xb9va2HtE1W13OZUqpm5ubQCAQDAbD4bDT6bQXTpPJ5N3dXTabtbfsNhqNcrns8XgcDoeUMpFIuFyuycnJUqnU/66UUgcHB0T0/PzcfyYnPcBYQXwCDB+vDqempiqVSiKRcDqdp6enfr+/Wq3ag1Aptbq6mkqlkslkLpc7Pz+3LOvo6Ih/KqXktGtLXF47+ny+RqPhdrvf399XVlaI6OXlpf9d6ZZavr1e+pd2AUbWJwAMFW8pISK/38/tQkKIarXaOdPj8RBRs9nkw7OzMyIyTZMP4/E4EUUikbaz/H4/EZXLZcMweHLbiV21Wi0hBBHxL221Wr1m4j8JjCf80QMMmWmanFKftijVGalx5p2dnfFhNpvlxOVg02lXLBbtZ7VaLc7jSCSi85JD+u3trc9d8fX5Toioa5x/fn5Wq1UiMgzjv3xygJ8MxVuAIbu4uCCiUChERNPT05xthUKhbdre3p4QYnNzs1KpBAKBWCxWLBZzuRw/8ader1uWRUSLi4v2sx4fH4mIu2fX1tboz1ehhmHwnpNejo+PiSgejzscDsMwHh4euk7jcW4SBhgriE+AYVJKcd+N1+vlkfX1dSJKpVL2adwfZFlWNBotFAo7Ozu1Ws2+NeXj44OIPB5P2/Pz7u/viUgIwT20RHR5eUlE4XC4z12VSqVyuUxEsViMiJaXlzmGO6XTaSJaWlr6+kcGGBHDXv4CjDVOKSFE2wj9s37L32tms9le1+Faqy7tap0V3a9UbnlOPB7XFxdC6K8/dRGY67rUrdQMMPKw+gQYps4tK79+/eLMy+fzbZOvr6/th7wLxT4yOztrP2w0GpZlGYah16m6cku9O2YrlQpHeDAY5JGlpSW9l7RUKnG1mf4sPSORyP98tD3AjzTs/AYYX9zXQx3LSl5r2hd8uqXINM1ms9lsNk3TNAxDLyK5/6itwadzScrT4vE4PxS+611xp67f77cPcqILIYQQfEHdrNSrqwhgtCE+AYag1Wpls1neE0JEhmHYE1TXbw3D0BVUzlRN99wy7oBtSzK+vn2Q23eFEL2yk6/TealyuawTlH/E2dx/9wvACMMbVwCGQClVr9fbBu1vVrE/kE+PK6Vub29dLtf8/Hxbi5CUcnJysu31LHwR+4iUMp/P+3y+XuVWKSUXdTtf8yKlTKfTW1tbExMTUkq3272xsXFycvLvnxZgFCE+AUaE1+sNhUKBQOC7fxE//IiI8KJsGGdoHQIYEYeHh722l/xFyE4AhtUnwIhQSrlcrlqt9n2pxtk5PT2dyWSQnTDmEJ8AoyOXy83MzNgfp/B3RaPRubm53d3db7o+wA+C+AQYHUopIvrW1ScWnQAM8QkAADAwtA4BAAAMDPEJAAAwMMQnAADAwBCfAAAAA0N8AgAADAzxCQAAMDDEJwAAwMAQnwAAAANDfAIAAAwM8QkAADCw3zozWlU3YC9EAAAAAElFTkSuQmCC" } }, "cell_type": "markdown", "metadata": {}, "source": [ "Régions de stabilité absolue pour certaines méthodes RK explicites. \n", "Source: A. Quarteroni, F. Saleri, P. Gervasio *Calcul Scientifique: Cours, exercices corrigés et illustrations en Matlab et Octave.*\n", "\n", "![A_stab_RK.png](attachment:A_stab_RK.png)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercice : schémas RK à un étage\n", "Étudier les méthodes RK avec $s=1$:\n", "1. écrire le schéma **implicite** générique RK à un étage\n", "1. étudier son ordre de convergence\n", "1. écrire le schéma **explicite** générique RK à un étage\n", "1. étudier son ordre de convergence\n", "1. étudier la A-stabilité" ] }, { "cell_type": "markdown", "metadata": { "cell_style": "center" }, "source": [ "### Correction : écriture d'un schéma RK à un étage quelconque \n", "Ces méthodes ont pour matrice de Butcher\n", "$$\n", "\\begin{array}{c|cc}\n", "c_1 & a_{11}\\\\\n", "\\hline\n", " & b_1\n", "\\end{array}\n", "\\text{ avec }\n", "\\begin{cases}\n", "c_1=a_{11}\\\\\n", "b_1=1\n", "\\end{cases}\n", "\\quad\n", "\\implies\n", "\\begin{array}{c|cc}\n", "c_1 & c_1\\\\\n", "\\hline\n", " & 1\n", "\\end{array}\n", "\\quad\n", "\\implies\n", "\\begin{cases}\n", "u_0\t = y_0 \\\\\n", "K_1 = \\varphi\\left(t_n+hc_1,u_n+h c_{1}K_1\\right)\\\\\n", "u_{n+1} = u_n + h K_1 & n=0,1,\\dots N-1\n", "\\end{cases}\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "cell_style": "center" }, "source": [ "Par exemple, si on choisit $c_1=1$ on trouve le schéma implicite \n", "$$\n", "\\begin{array}{c|cc}\n", "1 & 1\\\\\n", "\\hline\n", " & 1\n", "\\end{array}\n", "\\quad\n", "\\implies\n", "\\begin{cases}\n", "u_0\t = y_0 \\\\\n", "K_1 = \\varphi\\left(t_{n+1},u_n+h K_1\\right)\\\\\n", "u_{n+1} = u_n + h K_1 & n=0,1,\\dots N-1\n", "\\end{cases}\n", "$$\n", "\n", "Notons que ce schéma n'est rien d'autre que le schéma d'Euler implicite car $u_{n+1} = u_n+h K_1$, donc si on remplace le $u_n+h K_1$ dans la définition de $K_1$ par $u_{n+1}$ on obtient\n", "$$\n", "K_1\n", "= \\varphi\\left(t_{n+1},u_n+h K_1\\right)\n", "= \\varphi\\left(t_{n+1},u_{n+1}\\right)\n", "\\quad\\implies\\quad\n", "u_{n+1} = u_n + h K_1=u_n + h \\varphi\\left(t_{n+1},u_{n+1}\\right)\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "cell_style": "center" }, "source": [ "### Correction : étude de l'ordre d'un schéma RK à un étage quelconque \n", "Si on cherche une méthode d'ordre $2$ alors il faut que $1\\times c_1=\\frac{1}{2}$: il existe une seule méthode d'ordre 2, elle est implicite et sa matrice de Butcher est\n", "$$\n", "\\begin{array}{c|cc}\n", "\\frac{1}{2} & \\frac{1}{2}\\\\\n", "\\hline\n", " & 1\n", "\\end{array}\n", "\\quad\n", "\\implies\n", "\\begin{cases}\n", "u_0\t = y_0 \\\\\n", "K_1 = \\varphi\\left(t_n+\\frac{h}{2},u_n+\\frac{h}{2}K_1\\right)\\\\\n", "u_{n+1} = u_n + h K_1 & n=0,1,\\dots N-1\n", "\\end{cases}\n", "$$\n", "\n", "Notons que $u_{n+1} = u_n+h K_1$, donc si on remplace le $u_n+\\frac{h}{2} K_1=\\frac{u_n+\\left(u_n+h K_1\\right)}{2}$ dans la définition de $K_1$ par $u_{n+1}$ on obtient le schéma RK1_M\n", "$$\n", "K_1\n", "= \\varphi\\left(t_{n}+\\frac{h}{2},u_n+\\frac{h}{2} K_1\\right)\n", "= \\varphi\\left(t_{n}+\\frac{h}{2},\\frac{u_n+u_{n+1}}{2}\\right)\n", "\\quad\\implies\\quad\n", "u_{n+1} = u_n + h K_1=u_n + h \\varphi\\left(\\frac{t_n+t_{n+1}}{2},\\frac{u_n+u_{n+1}}{2}\\right)\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "cell_style": "center" }, "source": [ "### Correction : écriture d'un schéma RK à un étage explicite \n", "Si le schéma est explicite alors\n", "$$\n", "\\begin{array}{c|cc}\n", "0 & 0\\\\\n", "\\hline\n", " & 1\n", "\\end{array}\n", "\\quad\n", "\\implies\n", "\\begin{cases}\n", "u_0\t = y_0 \\\\\n", "K_1 = \\varphi\\left(t_n,u_n\\right)\\\\\n", "u_{n+1} = u_n + h K_1 & n=0,1,\\dots N-1\n", "\\end{cases}\n", "$$\n", "\n", "Il existe un seul schéma explicite RK à un étage, il s'agit du schéma d'Euler explicite." ] }, { "cell_type": "markdown", "metadata": { "cell_style": "center" }, "source": [ "### Correction : étude de l'ordre d'un schéma RK à un étage explicite \n", "Si le schéma est explicite alors il est d'ordre 1." ] }, { "cell_type": "markdown", "metadata": { "cell_style": "center" }, "source": [ "### Correction : étude de la A-stabilité \n", "\n", "Une méthode de Runge-Kutta à 1 étage pour $y'(t)=-\\beta y(t)$ s'écrit:\n", "$$\\begin{cases}\n", "u_0=y_0,\\\\\n", "u_{n+1}=u_{n}+h K_1& n=0,1,\\dots N-1\\\\\n", "K_1=-\\beta \\left( u_n+hc_1K_1 \\right).\n", "\\end{cases}$$\n", "On trouve donc $(1+\\beta hc_1)K_1=-\\beta u_n$\n", "ainsi\n", "$$\\begin{cases}\n", "u_0\t = y_0 \\\\\n", "u_{n+1} = \\left(1-\\frac{(\\beta h)}{1+c_1(\\beta h)}\\right)u_n & n=0,1,\\dots N-1\n", "\\end{cases}$$\n", "Par induction on obtient\n", "$$\n", "u_{n}=\\left(1-\\frac{\\beta h}{1+c_1\\beta h}\\right)^nu_0.\n", "$$\n", "Par conséquent, $\\lim\\limits_{n\\to+\\infty}u_n=0$ si et seulement si \n", "$$\n", "\\left|1-\\frac{\\beta h}{1+c_1\\beta h}\\right|<1.\n", "$$\n", "Notons $x$ le produit $\\beta h>0$ (donc $x>0$) et $q$ la fonction $q(x)=1-\\frac{x}{1+c_1x}$ (avec $c_1\\in[0;1]$). \n", "\\begin{align*}\n", "|q(x)|<1\n", "\\quad\\iff\\quad\n", "-1<1-\\frac{x}{1+c_1x}<1\n", "\\quad\\iff\\quad\n", "\\begin{cases}\n", "\\frac{x}{1+c_1x}<2\\\\\n", "\\frac{x}{1+c_1x}>0\n", "\\end{cases}\n", "\\quad\\stackrel{\\substack{x>0\\\\c_1\\ge0}}{\\iff}\\quad\n", "\\frac{x}{1+c_1x}<2\n", "\\quad\\stackrel{c_1x\\ge0}{\\iff}\\quad\n", "x<2(1+c_1x)\n", "\\quad\\iff\\quad\n", "(1-2c_1)x<2\n", "\\quad\\iff\\quad\n", "\\begin{cases}\n", "x<\\frac{2}{1-2c_1}&\\text{si }1-2c_1>0\\\\\n", "\\forall x&\\text{sinon.}\n", "\\end{cases}\n", "\\end{align*}\n", "Conclusion:\n", "+ si $c_1\\ge \\frac{1}{2}$ le schéma est inconditionnellement A-stable,\n", "+ si $c_1<\\frac{1}{2}$ le schéma est A-stable ssi $\\beta h <\\frac{2}{1-2c_1}$.\n", "\n", "Pour nos exemples, on obtient:\n", "+ si $c_1=1$ (Euler Implicite) alors le schéma est inconditionnellement A-stable et d'ordre 1,\n", "+ si $c_1=\\frac{1}{2}$ alors le schéma est inconditionnellement A-stable et d'ordre 2,\n", "+ si $c_1=0$ (Euler Explicite) alors le schéma est A-stable ssi $h<\\frac{2}{\\beta}$ et il est d'ordre 1." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercice : schémas RK à deux étages\n", "\n", "Étudier les méthodes RK avec $s=2$:\n", "1. écrire le schéma **implicite** générique RK à deux étages\n", "1. écrire le schéma **semi-implicite** générique RK à deux étages\n", "1. écrire le schéma **explicite** générique RK à deux étages\n", "1. étudier l'ordre de convergences des schémas explicites\n", "1. étudier la A-stabilité des schémas explicites" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Correction : écriture d'un schéma RK à deux étages quelconque \n", "Les méthodes avec $s=2$ ont pour matrice de Butcher\n", "$$\n", "\\begin{array}{c|cc}\n", "c_1 & a_{11} & a_{12}\\\\\n", "c_2 & a_{21} & a_{22}\\\\\n", "\\hline\n", " & b_1 & b_2\n", "\\end{array}\n", "\\text{ avec }\n", "\\begin{cases}\n", "c_1=a_{11}+a_{12}\\\\\n", "c_2=a_{21}+a_{22}\\\\\n", "b_1+b_2=1\n", "\\end{cases}\n", "\\quad\n", "\\implies\n", "\\begin{array}{c|cc}\n", "c_1 & a_{11} & c_1-a_{11}\\\\\n", "c_2 & a_{21} & c_2-a_{21}\\\\\n", "\\hline\n", " & 1-b_2 & b_2\n", "\\end{array}\n", "\\quad\n", "\\implies\n", "\\begin{cases}\n", "u_0\t = y_0 \\\\\n", "K_1 = \\varphi\\left(t_n+hc_1,u_n+h(a_{11}K_1+ (c_1-a_{11})K_2)\\right)\\\\\n", "K_2 = \\varphi\\left(t_n+hc_2,u_n+h(a_{21}K_1+ (c_2-a_{21})K_2)\\right)\\\\\n", "u_{n+1} = u_n + h\\left((1-b_2)K_1+b_2K_2\\right) & n=0,1,\\dots N-1\n", "\\end{cases}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Voici un exemple, appelé méthode de Gauss:\n", "$$\n", "\\begin{array}{c|cc}\n", "\\frac{1}{2}-\\gamma & \\frac{1}{4} & \\frac{1}{4}-\\gamma\\\\\n", "\\frac{1}{2}+\\gamma & \\frac{1}{4}+\\gamma & \\frac{1}{4}\\\\\n", "\\hline\n", " & \\frac{1}{2} & \\frac{1}{2}\n", "\\end{array}\n", "\\qquad\\text{avec }\\gamma=\\frac{\\sqrt{3}}{6}\n", "$$\n", "\n", "Cette méthode est d'ordre 4 et on peut monter que c'est la seule méthode d'ordre 4 à deux étages." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Consistance\n", "True\n", "True\n", "True\n", "\n", "Ordre 2\n", "True\n", "\n", "Ordre 3\n", "True\n", "True\n", "\n", "Ordre 4\n", "True\n", "True\n", "True\n", "True\n" ] } ], "source": [ "from sympy import *\n", "gamma=sqrt(3)/6\n", "c=[1/2-gamma,1/2+gamma]\n", "b=[1/2,1/2]\n", "A=[[1/4,1/4-gamma],[1/4+gamma,1/4]]\n", "s=len(c)\n", "print(\"Consistance\")\n", "print(sum(b)==1)\n", "for i in range(s):\n", " print(sum(A[i]).simplify()==c[i])\n", "print(\"\\nOrdre 2\")\n", "print(sum([b[i]*c[i] for i in range(s)]).simplify()==1/2)\n", "print(\"\\nOrdre 3\")\n", "print(sum([b[i]*c[i]**2 for i in range(s)]).simplify()==1/3)\n", "print(sum([b[i]*A[i][j]*c[j] for i in range(s) for j in range(s)]).simplify()==1/6)\n", "print(\"\\nOrdre 4\")\n", "print(sum([b[i]*c[i]**3 for i in range(s)]).simplify()==1/4)\n", "print(sum([b[i]*c[i]*A[i][j]*c[j] for i in range(s) for j in range(s)]).simplify()==1/8)\n", "print(sum([b[i]*A[i][j]*c[j]**2 for i in range(s) for j in range(s)]).simplify()==1/12)\n", "print(sum([b[i]*A[i][j]*A[j][k]*c[k] for i in range(s) for j in range(s) for k in range(s)]).simplify()==1/24)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Correction : écriture d'un schéma RK à deux étages semi-implicite \n", "Les méthodes semi-implicites vérifient, de plus, $c_1=a_{11}$ ainsi la matrice de Butcher est\n", "$$\n", "\\begin{array}{c|cc}\n", "c_1 & c_1 & 0\\\\\n", "c_2 & a_{21} & c_2-a_{21}\\\\\n", "\\hline\n", " & 1-b_2 & b_2\n", "\\end{array}\n", "\\quad\n", "\\implies\n", "\\begin{cases}\n", "u_0\t = y_0 \\\\\n", "K_1 = \\varphi\\left(t_n+hc_1,u_n+hc_{1}K_1\\right)\\\\\n", "K_2 = \\varphi\\left(t_n+hc_2,u_n+h(a_{21}K_1+ (c_2-a_{21})K_2)\\right)\\\\\n", "u_{n+1} = u_n + h\\left((1-b_2)K_1+b_2K_2\\right) & n=0,1,\\dots N-1\n", "\\end{cases}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Par exemple, si on choisit $c_1=0$, $c_2=1$ et $a_{21}=a_{22}=b_1=b_2=\\frac{1}{2}$ on trouve le schéma semi-implicite \n", "$$\n", "\\begin{array}{c|cc}\n", "0 & 0 & 0\\\\\n", "1 & \\frac{1}{2} & \\frac{1}{2}\\\\\n", "\\hline\n", " & \\frac{1}{2} & \\frac{1}{2}\n", "\\end{array}\n", "\\quad\n", "\\implies\n", "\\begin{cases}\n", "u_0\t = y_0 \\\\\n", "K_1 = \\varphi\\left(t_{n},u_n\\right)\\\\\n", "K_2 = \\varphi\\left(t_{n+1},u_n+\\frac{h}{2}(K_1+K_2)\\right)\\\\\n", "u_{n+1} = u_n + \\frac{h}{2}(K_1+K_2) & n=0,1,\\dots N-1\n", "\\end{cases}\n", "$$\n", "\n", "Notons que ce schéma n'est rien d'autre que le schéma de Crank-Nicolson car $u_{n+1} = u_n+\\frac{h}{2}(K_1+K_2)$, donc si on remplace le $u_n+\\frac{h}{2}(K_1+K_2)$ dans la définition de $K_2$ par $u_{n+1}$ on obtient\n", "$$\n", "K_2 \n", "= \\varphi\\left(t_{n+1},u_n+\\frac{h}{2}(K_1+K_2)\\right)\n", "= \\varphi\\left(t_{n+1},u_{n+1}\\right)\n", "$$\n", "ainsi\n", "$$\n", "u_{n+1} = u_n+\\frac{h}{2}(K_1+K_2)=u_n + \\frac{h}{2}(\\varphi\\left(t_{n},u_n\\right)+\\varphi\\left(t_{n+1},u_{n+1}\\right))\n", "$$\n", "Il s'agit donc d'un schéma d'ordre 2." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Consistance\n", "True\n", "True\n", "True\n", "\n", "Ordre 2\n", "True\n", "\n", "Ordre 3\n", "False\n", "False\n" ] } ], "source": [ "from sympy import *\n", "c=[0,1]\n", "b=[1/2,1/2]\n", "A=[[0,0],[1/2,1/2]]\n", "s=len(c)\n", "print(\"Consistance\")\n", "print(sum(b)==1)\n", "for i in range(s):\n", " print(sum(A[i])==c[i])\n", "print(\"\\nOrdre 2\")\n", "print(sum([b[i]*c[i] for i in range(s)])==1/2)\n", "print(\"\\nOrdre 3\")\n", "print(sum([b[i]*c[i]**2 for i in range(s)])==1/3)\n", "print(sum([b[i]*A[i][j]*c[j] for i in range(s) for j in range(s)])==1/6)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Un autre exemple est la méthode dite DIRK (Diagonally Implicit Runge-Kutta) qui est d'ordre 3: \n", "$$\n", "\\begin{array}{c|cc}\n", "\\frac{1}{3} & \\frac{1}{3} & 0\\\\\n", "1 & 1 & 0\\\\\n", "\\hline\n", " & \\frac{3}{4} & \\frac{1}{4}\n", "\\end{array}\n", "\\quad\n", "\\implies\n", "\\begin{cases}\n", "u_0\t = y_0 \\\\\n", "K_1 = \\varphi\\left(t_{n}+\\frac{1}{3}h,u_n+\\frac{1}{3}hK_1\\right)\\\\\n", "K_2 = \\varphi\\left(t_{n+1},u_n+hK_1\\right)\\\\\n", "u_{n+1} = u_n + \\frac{h}{4}(3K_1+K_2) & n=0,1,\\dots N-1\n", "\\end{cases}\n", "$$\n" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Consistance\n", "True\n", "True\n", "True\n", "\n", "Ordre 2\n", "True\n", "\n", "Ordre 3\n", "True\n", "True\n", "\n", "Ordre 4\n", "False\n", "False\n", "False\n", "False\n" ] } ], "source": [ "from sympy import *\n", "c=[Rational(1,3),1]\n", "b=[Rational(3,4),Rational(1,4)]\n", "A=[[Rational(1,3),0],[1,0]]\n", "s=len(c)\n", "print(\"Consistance\")\n", "print(sum(b)==1)\n", "for i in range(s):\n", " print(sum(A[i])==c[i])\n", "print(\"\\nOrdre 2\")\n", "print(sum([b[i]*c[i] for i in range(s)])==Rational(1,2))\n", "print(\"\\nOrdre 3\")\n", "print(sum([b[i]*c[i]**2 for i in range(s)])==Rational(1,3))\n", "print(sum([b[i]*A[i][j]*c[j] for i in range(s) for j in range(s)])==Rational(1,6))\n", "print(\"\\nOrdre 4\")\n", "print(sum([b[i]*c[i]**3 for i in range(s)])==Rational(1,4))\n", "print(sum([b[i]*c[i]*A[i][j]*c[j] for i in range(s) for j in range(s)])==Rational(1,8))\n", "print(sum([b[i]*A[i][j]*c[j]**2 for i in range(s) for j in range(s)])==Rational(1,12))\n", "print(sum([b[i]*A[i][j]*A[j][k]*c[k] for i in range(s) for j in range(s) for k in range(s)])==Rational(1,24))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Correction : écriture d'un schéma RK à deux étages explicite\n", "Les méthodes avec $s=2$ explicites ont pour matrice de Butcher\n", "$$\n", "\\begin{array}{c|cc}\n", "0 & 0 & 0\\\\\n", "c_2 & c_{2} & 0\\\\\n", "\\hline\n", " & 1-b_2 & b_2 \n", "\\end{array}\n", "\\quad\n", "\\implies\n", "\\begin{cases}\n", "u_0\t = y_0 \\\\\n", "K_1 = \\varphi\\left(t_n,u_n\\right)\\\\\n", "K_2 = \\varphi\\left(t_n+hc_2,u_n+hc_2K_1\\right)\\\\\n", "u_{n+1} = u_n + h\\left((1-b_2)K_1+b_2K_2\\right) & n=0,1,\\dots N-1\n", "\\end{cases}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Correction : étude de l'ordre d'un schéma RK à deux étages explicite \n", "\n", "+ Un schéma RK à deux étage explicite est au mieux d'ordre 2.\n", "+ Pour avoir l'ordre 2 il faut que $b_2c_2=\\frac{1}{2}$.\n", "\n", "Conclusion: un schéma RK à deux étages explicite d'ordre 2 s'écrit \n", "$$\n", "\\begin{array}{c|cc}\n", "0 & 0 & 0\\\\\n", "\\frac{1}{2\\alpha} & \\frac{1}{2\\alpha} & 0\\\\\n", "\\hline\n", " & 1-\\alpha & \\alpha \n", "\\end{array}\n", "\\quad\n", "\\implies\n", "\\begin{cases}\n", "u_0\t = y_0 \\\\\n", "K_1 = \\varphi\\left(t_n,u_n\\right)\\\\\n", "K_2 = \\varphi\\left(t_n+\\frac{h}{2}\\alpha,u_n+\\frac{h}{2}\\alpha K_1\\right)\\\\\n", "u_{n+1} = u_n + h\\left((1-\\alpha)K_1+\\alpha K_2\\right) & n=0,1,\\dots N-1\n", "\\end{cases}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Voici deux cas particuliers:\n", "+ $\\alpha=\\frac{1}{2}$: schéma de Heun \n", "$$\n", "\\begin{array}{c|cc}\n", "0 & 0 & 0\\\\\n", "1 & 1 & 0\\\\\n", "\\hline\n", " & \\frac{1}{2} & \\frac{1}{2}\n", "\\end{array}\n", "\\quad\n", "\\implies\n", "\\begin{cases}\n", "u_0\t = y_0 \\\\\n", "K_1 = \\varphi(t_n,u_n)\\\\\n", "K_2 = \\varphi(t_{n+1},u_{n}+hK_1)\\\\\n", "u_{n+1} = u_n + h\\left(\\frac{1}{2}K_1+\\frac{1}{2}K_2\\right) & n=0,1,\\dots N-1\n", "\\end{cases}\n", "$$\n", "\n", "+ $\\alpha=1$: schéma d'Euler Modifié (ou du point milieu) \n", "$$\n", "\\begin{array}{c|cc}\n", "0 & 0 & 0\\\\\n", "\\frac{1}{2} & \\frac{1}{2} & 0\\\\\n", "\\hline\n", " & 0 & 1\n", "\\end{array}\n", "\\quad\n", "\\implies\n", "\\begin{cases}\n", "u_0\t = y_0 \\\\\n", "K_1 = \\varphi(t_n,u_n)\\\\\n", "K_2 = \\varphi(t_{n}+\\frac{1}{2}h,u_{n}+\\frac{1}{2}hK_1)\\\\\n", "u_{n+1} = u_n + h K_2 & n=0,1,\\dots N-1\n", "\\end{cases}\n", "$$\n", "\n", "+ $\\alpha=\\frac{2}{3}$: schéma de Ralston \n", "$$\n", "\\begin{array}{c|cc}\n", "0 & 0 & 0\\\\\n", "\\frac{3}{4} & \\frac{3}{4} & 0\\\\\n", "\\hline\n", " & \\frac{1}{3} & \\frac{2}{3}\n", "\\end{array}\n", "\\quad\n", "\\implies\n", "\\begin{cases}\n", "u_0\t = y_0 \\\\\n", "K_1 = \\varphi(t_n,u_n)\\\\\n", "K_2 = \\varphi(t_{n}+\\frac{3}{4}h,u_{n}+\\frac{3}{4}hK_1)\\\\\n", "u_{n+1} = u_n + \\frac{h}{3} (K_1+2K_2) & n=0,1,\\dots N-1\n", "\\end{cases}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Correction : étude de la A-stabilité des schémas explicites\n", "Le schéma donné appliqué à cette équation devient\n", "$$\\begin{cases}\n", "u_0\t = y_0 \\\\\n", "u_{n+1} = \\left(1-(\\beta h)+c_2b_2(\\beta h)^2\\right)u_n & n=0,1,\\dots N-1\n", "\\end{cases}$$\n", "Par induction on obtient\n", "$$\n", "u_{n}=\\left(1-(\\beta h)+c_2b_2(\\beta h)^2\\right)^nu_0.\n", "$$\n", "Par conséquent, $\\lim\\limits_{n\\to+\\infty}u_n=0$ si et seulement si \n", "$$\n", "\\left|1-(\\beta h)+c_2b_2(\\beta h)^2\\right|<1.\n", "$$\n", "Notons $x$ le produit $\\beta h>0$ (donc $x>0$) et $q$ le polynôme $q(x)=1-x+c_2b_2x^2$. \n", "+ Si $c_2b_2=0$, le schéma se réduit au schéma d'Euler explicite.\n", "+ Si $c_2b_2>0$, c'est une parabole convexe et le sommet est situé en $\\left(\\frac{1}{2c_2b_2}>0,1-\\frac{1}{4c_2b_2}\\right)$.\n", "$$ \n", "|q(x)|<1\n", "\\quad\\iff\\quad\n", "x<\\frac{1}{c_2b_2}\n", "$$\n", "En particulier, lorsque le schéma est d'ordre 2 on a $b_2c_2=\\frac{1}{2}$ et le schéma est A-stable ssi $h<\\frac{2}{\\beta}$.\n", "+ Si $c_2b_2<0$, c'est une parabole concave et le sommet est situé en $\\left(\\frac{1}{2c_2b_2}<0,1-\\frac{1}{4c_2b_2}\\right)$.\n", "$$ \n", "|q(x)|<1\n", "\\quad\\iff\\quad\n", "x<\\frac{1-\\sqrt{1-8b_2c_2}}{2b_2c_2}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exemples à 3 étages explicites\n", "\n", "Un schéma RK à $3$ étages explicite a pour matrice de Butcher\n", "$$\n", "\\begin{array}{c|ccc}\n", "0 & 0 & 0 & 0 \\\\\n", "c_2 & a_{21} & 0&0\\\\\n", "c_3 & a_{31} &a_{32}&0\\\\\n", "\\hline\n", " & b_1 & b_2 & b_3 \n", "\\end{array}\n", "\\text{ avec }\n", "\\begin{cases}\n", "c_2=a_{21}\\\\\n", "c_3=a_{31}+a_{32}\\\\\n", "b_1+b_2+b_3=1\n", "\\end{cases}\n", "\\quad\n", "\\implies\n", "\\begin{cases}\n", "u_0\t = y_0 \\\\\n", "K_1 = \\varphi\\left(t_n,u_n\\right)\\\\\n", "K_2 = \\varphi\\left(t_n+hc_2,u_n+hc_2K_1\\right)\\\\\n", "K_3 = \\varphi\\left(t_n+hc_3,u_n+h(a_{31}K_1+ a_{32}K_2)\\right)\\\\\n", "u_{n+1} = u_n + h\\left(b_1K_1+b_2K_2+b_3K_3\\right) & n=0,1,\\dots N-1\n", "\\end{cases}\n", "$$\n", "\n", "Voici trois exemples:\n", "\n", "+ schéma de Heun à trois étages\n", "$$\n", "\\begin{array}{c|ccc}\n", "0 & 0 & 0 & 0 \\\\\n", "\\frac{1}{3} & \\frac{1}{3} & 0&0\\\\\n", "\\frac{2}{3} & 0 &\\frac{2}{3}&0\\\\\n", "\\hline\n", " & \\frac{1}{4} & 0 & \\frac{3}{4} \n", "\\end{array}\n", "\\quad\n", "\\implies\n", "\\begin{cases}\n", "u_0\t = y_0 \\\\\n", "K_1 = \\varphi\\left(t_n,u_n\\right)\\\\\n", "K_2 = \\varphi\\left(t_n+\\frac{h}{3},u_n+\\frac{h}{3}K_1\\right)\\\\\n", "K_3 = \\varphi\\left(t_n+\\frac{2}{3}h,u_n+\\frac{2}{3}hK_2\\right)\\\\\n", "u_{n+1} = u_n + \\frac{h}{4}\\left(K_1+3K_3\\right) & n=0,1,\\dots N-1\n", "\\end{cases}\n", "$$\n", "\n", "Étant une méthode explicite à 3 étapes, elle est au mieux d'ordre 3. Vérifions qu'elle est effectivement d'ordre 3:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Consistance\n", "True\n", "True\n", "True\n", "True\n", "\n", "Ordre 2\n", "True\n", "\n", "Ordre 3\n", "True\n", "True\n" ] } ], "source": [ "#from sympy import *\n", "c=[0,1/3,2/3]\n", "b=[1/4,0,3/4]\n", "A=[[0,0,0],[1/3,0,0],[0,2/3,0]]\n", "s=len(c)\n", "print(\"Consistance\")\n", "print(sum(b)==1)\n", "for i in range(s):\n", " print(sum(A[i])==c[i])\n", "print(\"\\nOrdre 2\")\n", "print(sum([b[i]*c[i] for i in range(s)])==1/2)\n", "print(\"\\nOrdre 3\")\n", "print(sum([b[i]*c[i]**2 for i in range(s)])==1/3)\n", "print(sum([b[i]*A[i][j]*c[j] for i in range(s) for j in range(s)])==1/6)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "+ schéma de Kutta\n", "$$\n", "\\begin{array}{c|ccc}\n", "0 & 0 & 0 & 0 \\\\\n", "\\frac{1}{2} & \\frac{1}{2} & 0&0\\\\\n", "1 & -1 &2&0\\\\\n", "\\hline\n", " & \\frac{1}{6} & \\frac{2}{3} & \\frac{1}{6} \n", "\\end{array}\n", "\\quad\n", "\\implies\n", "\\begin{cases}\n", "u_0\t = 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+1},u_n+h(-K_1+2K_2)\\right)\\\\\n", "u_{n+1} = u_n + \\frac{h}{6}\\left(K_1+4K_2+K_3\\right) & n=0,1,\\dots N-1\n", "\\end{cases}\n", "$$\n", "\n", "Étant une méthode explicite à 3 étapes, elle a au mieux d'ordre 3. On vérifie ci-dessous qu'elle est d'ordre 3:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Consistance\n", "True\n", "True\n", "True\n", "True\n", "\n", "Ordre 2\n", "True\n", "\n", "Ordre 3\n", "True\n", "True\n" ] } ], "source": [ "# from fractions import Fraction\n", "# c=[0,Fraction(1,2),1]\n", "# b=[Fraction(1,6),Fraction(2,3),Fraction(1,6)]\n", "# A=[[0,0,0],[Fraction(1,2),0,0],[-1,2,0]]\n", "# s=len(c)\n", "# print(\"Consistance\")\n", "# print(sum(b)==1)\n", "# for i in range(s):\n", "# print(sum(A[i])==c[i])\n", "# print(\"\\nOrdre 2\")\n", "# print(sum([b[i]*c[i] for i in range(s)])==Fraction(1,2))\n", "# print(\"\\nOrdre 3\")\n", "# print(sum([b[i]*c[i]**2 for i in range(s)])==Fraction(1,3))\n", "# print(sum([b[i]*A[i][j]*c[j] for i in range(s) for j in range(s)])==Fraction(1,6))\n", "\n", "from sympy import *\n", "c=[0,1/2,1]\n", "b=[S(1)/6,S(2)/3,S(1)/6] # S() pour Symbol\n", "A=[[0,0,0],[1/2,0,0],[-1,2,0]]\n", "s=len(c)\n", "print(\"Consistance\")\n", "print(sum(b)==1)\n", "for i in range(s):\n", " print(sum(A[i])==c[i])\n", "print(\"\\nOrdre 2\")\n", "print(sum([b[i]*c[i] for i in range(s)])==1/2)\n", "print(\"\\nOrdre 3\")\n", "print(sum([b[i]*c[i]**2 for i in range(s)])==1/3)\n", "print(sum([b[i]*A[i][j]*c[j] for i in range(s) for j in range(s)])==1/6)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "+ schéma \n", "$$\n", "\\begin{array}{c|ccc}\n", "0 & 0 & 0 & 0 \\\\\n", "\\frac{1}{2} & \\frac{1}{2} & 0&0\\\\\n", "1 & -1 &2&0\\\\\n", "\\hline\n", " & -\\frac{1}{6} & \\frac{4}{3} & -\\frac{1}{6} \n", "\\end{array}\n", "\\quad\n", "\\implies\n", "\\begin{cases}\n", "u_0\t = 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+1},u_n+h(-K_1+2K_2)\\right)\\\\\n", "u_{n+1} = u_n + \\frac{h}{6}\\left(-K_1+8K_2-K_3\\right) & n=0,1,\\dots N-1\n", "\\end{cases}\n", "$$\n", "\n", "Étant une méthode explicite à 3 étapes, elle a au mieux d'ordre 3. On vérifie ci-dessous qu'elle n'est que d'ordre 2:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Consistance\n", "True\n", "True\n", "True\n", "True\n", "\n", "Ordre 2\n", "True\n", "\n", "Ordre 3\n", "False\n", "False\n" ] } ], "source": [ "from fractions import Fraction\n", "c=[0,Fraction(1,2),1]\n", "b=[-Fraction(1,6),Fraction(4,3),-Fraction(1,6)]\n", "A=[[0,0,0],[Fraction(1,2),0,0],[-1,2,0]]\n", "s=len(c)\n", "print(\"Consistance\")\n", "print(sum(b)==1)\n", "for i in range(s):\n", " print(sum(A[i])==c[i])\n", "print(\"\\nOrdre 2\")\n", "print(sum([b[i]*c[i] for i in range(s)])==Fraction(1,2))\n", "print(\"\\nOrdre 3\")\n", "print(sum([b[i]*c[i]**2 for i in range(s)])==Fraction(1,3))\n", "print(sum([b[i]*A[i][j]*c[j] for i in range(s) for j in range(s)])==Fraction(1,6))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "+ schéma Radau-IIa (d'ordre 5)\n", "$$\n", "\\begin{array}{c|ccc}\n", "\\frac{4-\\sqrt{6}}{10} & \\frac{88-7\\sqrt{6}}{360} & \\frac{296-169\\sqrt{6}}{1800} & \\frac{-2+3\\sqrt{6}}{225} \\\\\n", "\\frac{4+\\sqrt{6}}{10} & \\frac{296+169\\sqrt{6}}{1800} & \\frac{88+7\\sqrt{6}}{360}&\\frac{-2-3\\sqrt{6}}{225}\\\\\n", "1 & \\frac{16-\\sqrt{6}}{36} &\\frac{16+\\sqrt{6}}{36}&\\frac{1}{9}\\\\\n", "\\hline\n", " & \\frac{16-\\sqrt{6}}{36} &\\frac{16+\\sqrt{6}}{36}& \\frac{1}{9}\n", "\\end{array}\n", "\\quad\n", "\\implies\n", "\\begin{cases}\n", "u_0\t = y_0 \\\\\n", "K_1 = \\varphi\\left(t_n+\\frac{4-\\sqrt{6}}{10}h,u_n+h\\left(\\frac{88-7\\sqrt{6}}{360} K_1+\\frac{296-169\\sqrt{6}}{1800} K_2+\\frac{-2+3\\sqrt{6}}{225}K_3\\right)\\right)\\\\\n", "K_2 = \\varphi\\left(t_n+\\frac{4+\\sqrt{6}}{10}h,u_n+h\\left(\\frac{296+169\\sqrt{6}}{1800}K_1+\\frac{88+7\\sqrt{6}}{360}K_2+\\frac{-2-3\\sqrt{6}}{225}K_3\\right)\\right)\\\\\n", "K_3 = \\varphi\\left(t_n+h,u_n+h\\left(\\frac{16-\\sqrt{6}}{36}K_1+\\frac{16+\\sqrt{6}}{36}K_2+\\frac{1}{9}K_3\\right)\\right)\\\\\n", "u_{n+1} = u_n + h\\left(\\frac{16-\\sqrt{6}}{36}K_1+\\frac{16+\\sqrt{6}}{36}K_2+\\frac{1}{9}K_3\\right) & n=0,1,\\dots N-1\n", "\\end{cases}\n", "$$\n", "\n", "Étant une méthode implicite à 3 étapes, elle a au mieux d'ordre 6. On vérifie ci-dessous qu'elle est d'ordre 5:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Consistance\n", "True\n", "True\n", "True\n", "True\n", "\n", "Ordre 2\n", "True\n", "\n", "Ordre 3\n", "True\n", "True\n", "\n", "Ordre 4\n", "True\n", "True\n", "True\n", "True\n" ] } ], "source": [ "from sympy import *\n", "sq6=sqrt(6)\n", "c=[(4-sq6)/10, (4+sq6)/10, 1]\n", "b=[(16-sq6)/36,(16+sq6)/36,S(1)/9]\n", "A=[[(88-7*sq6)/360,(296-169*sq6)/1800,(-2+3*sq6)/225],\n", " [(296+169*sq6)/1800,(88+7*sq6)/360,(-2-3*sq6)/225],\n", " [(16-sq6)/36,(16+sq6)/36,S(1)/9]]\n", "s=len(c)\n", "print(\"Consistance\")\n", "print(sum(b)==1)\n", "for i in range(s):\n", " print(sum(A[i])==c[i])\n", "print(\"\\nOrdre 2\")\n", "print(sum([b[i]*c[i] for i in range(s)]).simplify()==S(1)/2)\n", "print(\"\\nOrdre 3\")\n", "print(sum([b[i]*c[i]**2 for i in range(s)]).simplify()==S(1)/3)\n", "print(sum([b[i]*A[i][j]*c[j] for i in range(s) for j in range(s)]).simplify()==S(1)/6)\n", "print(\"\\nOrdre 4\")\n", "print(sum([b[i]*c[i]**3 for i in range(s)]).simplify()==1/4)\n", "print(sum([b[i]*c[i]*A[i][j]*c[j] for i in range(s) for j in range(s)]).simplify()==1/8)\n", "print(sum([b[i]*A[i][j]*c[j]**2 for i in range(s) for j in range(s)]).simplify()==S(1)/12)\n", "print(sum([b[i]*A[i][j]*A[j][k]*c[k] for i in range(s) for j in range(s) for k in range(s)]).simplify()==S(1)/24)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exemples à 4 étages explicites\n", "\n", "Un schéma RK à $4$ étages explicite a pour matrice de Butcher\n", "$$\n", "\\begin{array}{c|cccc}\n", "0 & 0 & 0 & 0 & 0\\\\\n", "c_2 & a_{21} & 0&0&0\\\\\n", "c_3 & a_{31} &a_{32}&0&0\\\\\n", "c_4 & a_{41} &a_{42}&a_{43}&0\\\\\n", "\\hline\n", " & b_1 & b_2 & b_3 & b_4\n", "\\end{array}\n", "\\text{ avec }\n", "\\begin{cases}\n", "c_2=a_{21}\\\\\n", "c_3=a_{31}+a_{32}\\\\\n", "c_4=a_{41}+a_{42}+a_{43}\\\\\n", "b_1+b_2+b_3+b_4=1\n", "\\end{cases}\n", "\\quad\n", "\\implies\n", "\\begin{cases}\n", "u_0\t = y_0 \\\\\n", "K_1 = \\varphi\\left(t_n,u_n\\right)\\\\\n", "K_2 = \\varphi\\left(t_n+hc_2,u_n+h(a_{21}K_1)\\right)\\\\\n", "K_3 = \\varphi\\left(t_n+hc_3,u_n+h(a_{31}K_1+ a_{32}K_2)\\right)\\\\\n", "K_4 = \\varphi\\left(t_n+hc_4,u_n+h(a_{41}K_1+ a_{42}K_2+ a_{43}K_3)\\right)\\\\\n", "u_{n+1} = u_n + h\\left(b_1K_1+b_2K_2+b_3K_3+b_4K_4\\right) & n=0,1,\\dots N-1\n", "\\end{cases}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Voici deux exemples (le premier est le plus connus):\n", "+ schéma RK4-1 (\"La\" méthode de Runge-Kutta)\n", "$$\n", "\\begin{array}{c|cccc}\n", "0 & 0 & 0 & 0 & 0\\\\\n", "\\frac{1}{2} & \\frac{1}{2} & 0 & 0 & 0\\\\\n", "\\frac{1}{2} & 0 &\\frac{1}{2} &0&0\\\\\n", "1 & 0 & 0 & 1 & 0\\\\\n", "\\hline\n", " & \\frac{1}{6} & \\frac{1}{3} & \\frac{1}{3} & \\frac{1}{6}\n", "\\end{array}\n", "\\quad\n", "\\implies\n", "\\begin{cases}\n", "u_0\t = 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", "$$\n", "Il est bien d'ordre 4:" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Consistance\n", "True\n", "True\n", "True\n", "True\n", "True\n", "\n", "Ordre 2\n", "True\n", "\n", "Ordre 3\n", "True\n", "True\n", "\n", "Ordre 4\n", "True\n", "True\n", "True\n", "True\n" ] } ], "source": [ "from fractions import Fraction\n", "c=[0,Fraction(1,2),Fraction(1,2),1]\n", "b=[Fraction(1,6),Fraction(1,3),Fraction(1,3),Fraction(1,6)]\n", "A=[[0,0,0,0],[Fraction(1,2),0,0,0],[0,Fraction(1,2),0,0],[0,0,1,0]]\n", "s=len(c)\n", "print(\"Consistance\")\n", "print(sum(b)==1)\n", "for i in range(s):\n", " print(sum(A[i])==c[i])\n", "print(\"\\nOrdre 2\")\n", "print(sum([b[i]*c[i] for i in range(s)])==Fraction(1,2))\n", "print(\"\\nOrdre 3\")\n", "print(sum([b[i]*c[i]**2 for i in range(s)])==Fraction(1,3))\n", "print(sum([b[i]*A[i][j]*c[j] for i in range(s) for j in range(s)])==Fraction(1,6))\n", "print(\"\\nOrdre 4\")\n", "print(sum([b[i]*c[i]**3 for i in range(s)])==Fraction(1,4))\n", "print(sum([b[i]*c[i]*A[i][j]*c[j] for i in range(s) for j in range(s)])==Fraction(1,8))\n", "print(sum([b[i]*A[i][j]*c[j]**2 for i in range(s) for j in range(s)])==Fraction(1,12))\n", "print(sum([b[i]*A[i][j]*A[j][k]*c[k] for i in range(s) for j in range(s) for k in range(s)])==Fraction(1,24))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "+ schéma RK4-2 \n", "$$\n", "\\begin{array}{c|cccc}\n", "0 & 0 & 0 & 0 & 0\\\\\n", "\\frac{1}{4} & \\frac{1}{4} & 0 & 0 & 0\\\\\n", "\\frac{1}{2} & 0 &\\frac{1}{2} &0&0\\\\\n", "1 & 1 & -2 & 2 & 0\\\\\n", "\\hline\n", " & \\frac{1}{6} & 0 & \\frac{2}{3} & \\frac{1}{6}\n", "\\end{array}\n", "\\quad\n", "\\implies\n", "\\begin{cases}\n", "u_0\t = y_0 \\\\\n", "K_1 = \\varphi\\left(t_n,u_n\\right)\\\\\n", "K_2 = \\varphi\\left(t_n+\\frac{h}{4},u_n+\\frac{h}{4} 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_1-2K_2+2K_3)\\right)\\\\\n", "u_{n+1} = u_n + \\frac{h}{6}\\left(K_1+4K_3+K_4\\right) & n=0,1,\\dots N-1\n", "\\end{cases}\n", "$$\n", "Il est bien d'ordre 4:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Consistance\n", "True\n", "True\n", "True\n", "True\n", "True\n", "\n", "Ordre 2\n", "True\n", "\n", "Ordre 3\n", "True\n", "True\n", "\n", "Ordre 4\n", "True\n", "True\n", "True\n", "True\n" ] } ], "source": [ "from fractions import Fraction\n", "c=[0,Fraction(1,4),Fraction(1,2),1]\n", "b=[Fraction(1,6),0,Fraction(2,3),Fraction(1,6)]\n", "A=[[0,0,0,0],[Fraction(1,4),0,0,0],[0,Fraction(1,2),0,0],[1,-2,2,0]]\n", "s=len(c)\n", "print(\"Consistance\")\n", "print(sum(b)==1)\n", "for i in range(s):\n", " print(sum(A[i])==c[i])\n", "print(\"\\nOrdre 2\")\n", "print(sum([b[i]*c[i] for i in range(s)])==Fraction(1,2))\n", "print(\"\\nOrdre 3\")\n", "print(sum([b[i]*c[i]**2 for i in range(s)])==Fraction(1,3))\n", "print(sum([b[i]*A[i][j]*c[j] for i in range(s) for j in range(s)])==Fraction(1,6))\n", "print(\"\\nOrdre 4\")\n", "print(sum([b[i]*c[i]**3 for i in range(s)])==Fraction(1,4))\n", "print(sum([b[i]*c[i]*A[i][j]*c[j] for i in range(s) for j in range(s)])==Fraction(1,8))\n", "print(sum([b[i]*A[i][j]*c[j]**2 for i in range(s) for j in range(s)])==Fraction(1,12))\n", "print(sum([b[i]*A[i][j]*A[j][k]*c[k] for i in range(s) for j in range(s) for k in range(s)])==Fraction(1,24))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "+ schéma RK4-3 (règle 3/8)\n", "$$\n", "\\begin{array}{c|cccc}\n", "0 & 0 & 0 & 0 & 0\\\\\n", "\\frac{1}{3} & \\frac{1}{3} & 0 & 0 & 0\\\\\n", "\\frac{2}{3} &-\\frac{1}{3} & 1 &0&0\\\\\n", "1 & 1 & -1 & 1 & 0\\\\\n", "\\hline\n", " & \\frac{1}{8} & \\frac{3}{8} & \\frac{3}{8} & \\frac{1}{8}\n", "\\end{array}\n", "\\quad\n", "\\implies\n", "\\begin{cases}\n", "u_0\t = y_0 \\\\\n", "K_1 = \\varphi\\left(t_n,u_n\\right)\\\\\n", "K_2 = \\varphi\\left(t_n+\\frac{h}{3},u_n+\\frac{h}{3} K_1)\\right)\\\\\n", "K_3 = \\varphi\\left(t_n+\\frac{2}{3}h,u_n+\\frac{h}{3}(-K_1+3K_2)\\right)\\\\\n", "K_4 = \\varphi\\left(t_{n+1},u_n+h (K_1-K_2+K_3)\\right)\\\\\n", "u_{n+1} = u_n + \\frac{h}{8}\\left(K_1+3K_2+3K_3+K_4\\right) & n=0,1,\\dots N-1\n", "\\end{cases}\n", "$$\n", "Il est bien d'ordre 4:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Consistance\n", "True\n", "True\n", "True\n", "True\n", "True\n", "\n", "Ordre 2\n", "True\n", "\n", "Ordre 3\n", "True\n", "True\n", "\n", "Ordre 4\n", "True\n", "True\n", "True\n", "True\n" ] } ], "source": [ "from fractions import Fraction\n", "c=[0,Fraction(1,3),Fraction(2,3),1]\n", "b=[Fraction(1,8),Fraction(3,8),Fraction(3,8),Fraction(1,8)]\n", "A=[[0,0,0,0],[Fraction(1,3),0,0,0],[-Fraction(1,3),1,0,0],[1,-1,1,0]]\n", "s=len(c)\n", "print(\"Consistance\")\n", "print(sum(b)==1)\n", "for i in range(s):\n", " print(sum(A[i])==c[i])\n", "print(\"\\nOrdre 2\")\n", "print(sum([b[i]*c[i] for i in range(s)])==Fraction(1,2))\n", "print(\"\\nOrdre 3\")\n", "print(sum([b[i]*c[i]**2 for i in range(s)])==Fraction(1,3))\n", "print(sum([b[i]*A[i][j]*c[j] for i in range(s) for j in range(s)])==Fraction(1,6))\n", "print(\"\\nOrdre 4\")\n", "print(sum([b[i]*c[i]**3 for i in range(s)])==Fraction(1,4))\n", "print(sum([b[i]*c[i]*A[i][j]*c[j] for i in range(s) for j in range(s)])==Fraction(1,8))\n", "print(sum([b[i]*A[i][j]*c[j]**2 for i in range(s) for j in range(s)])==Fraction(1,12))\n", "print(sum([b[i]*A[i][j]*A[j][k]*c[k] for i in range(s) for j in range(s) for k in range(s)])==Fraction(1,24))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exemple à 5 étages explicite\n", "\n", "Un schéma RK à $5$ étages explicite a pour matrice de Butcher\n", "$$\n", "\\begin{array}{c|ccccc}\n", "0 & 0 & 0 & 0 & 0 &0 \\\\\n", "c_2 & a_{21} & 0&0&0&0\\\\\n", "c_3 & a_{31} &a_{32}&0&0&0\\\\\n", "c_4 & a_{41} &a_{42}&a_{43}&0&0\\\\\n", "c_5 & a_{51} &a_{52}&a_{53}&a_{54}&0\\\\\n", "\\hline\n", " & b_1 & b_2 & b_3 & b_4 & b_5\n", "\\end{array}\n", "\\quad\n", "\\implies\n", "\\begin{cases}\n", "u_0\t = y_0 \\\\\n", "K_1 = \\varphi\\left(t_n,u_n\\right)\\\\\n", "K_2 = \\varphi\\left(t_n+hc_2,u_n+h(a_{21}K_1)\\right)\\\\\n", "K_3 = \\varphi\\left(t_n+hc_3,u_n+h(a_{31}K_1+ a_{32}K_2\\right)\\\\\n", "K_4 = \\varphi\\left(t_n+hc_4,u_n+h(a_{41}K_1+ a_{42}K_2+ a_{43}K_3)\\right)\\\\\n", "K_5 = \\varphi\\left(t_n+hc_5,u_n+h(a_{51}K_1+ a_{52}K_2+ a_{53}K_3+ a_{54}K_4)\\right)\\\\\n", "u_{n+1} = u_n + h\\left(b_1K_1+b_2K_2+b_3K_3+b_4K_4+b_5K_5\\right) & n=0,1,\\dots N-1\n", "\\end{cases}\n", "$$\n", "avec \n", "$$\n", "\\begin{cases}\n", "c_2=a_{21}\\\\\n", "c_3=a_{31}+a_{32}\\\\\n", "c_4=a_{41}+a_{42}+a_{43}\\\\\n", "c_5=a_{51}+a_{52}+a_{53}+a_{54}\\\\\n", "b_1+b_2+b_3+b_4+b_5=1\n", "\\end{cases}\n", "$$\n", "\n", "Voici un exemple:\n", "+ schéma de Merson \n", "$$\n", "\\begin{array}{c|ccccc}\n", "0 & 0 & 0 & 0 & 0 &0 \\\\\n", "\\frac{1}{3} & \\frac{1}{3} & 0&0&0&0\\\\\n", "\\frac{1}{3} & \\frac{1}{6} &\\frac{1}{6}&0&0&0\\\\\n", "\\frac{1}{2} & \\frac{1}{8} &0&\\frac{3}{8}&0&0\\\\\n", "1 & \\frac{1}{2} &0&-\\frac{3}{2}&2&0\\\\\n", "\\hline\n", " & \\frac{1}{6} & 0 & 0 & \\frac{2}{3} & \\frac{1}{6}\n", "\\end{array}\n", "\\quad\n", "\\implies\n", "\\begin{cases}\n", "u_0\t = y_0 \\\\\n", "K_1 = \\varphi\\left(t_n,u_n\\right)\\\\\n", "K_2 = \\varphi\\left(t_n+\\frac{h}{3},u_n+\\frac{h}{3}K_1\\right)\\\\\n", "K_3 = \\varphi\\left(t_n+\\frac{h}{3},u_n+\\frac{h}{6}(K_1+K_2)\\right)\\\\\n", "K_4 = \\varphi\\left(t_n+\\frac{h}{2},u_n+\\frac{h}{8}(K_1+3K_3)\\right)\\\\\n", "K_5 = \\varphi\\left(t_n+h ,u_n+\\frac{h}{2}(K_1-3K_3+4K_4)\\right)\\\\\n", "u_{n+1} = u_n + \\frac{h}{6}\\left(K_1+4K_4+K_5\\right) & n=0,1,\\dots N-1\n", "\\end{cases}\n", "$$" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Consistance\n", "True\n", "True\n", "True\n", "True\n", "True\n", "True\n", "\n", "Ordre 2\n", "True\n", "\n", "Ordre 3\n", "True\n", "True\n", "\n", "Ordre 4\n", "True\n", "True\n", "True\n", "True\n", "\n", "Ordre 5\n", "False\n", "True\n" ] } ], "source": [ "from fractions import Fraction\n", "c=[0,Fraction(1,3),Fraction(1,3),Fraction(1,2),1]\n", "b=[Fraction(1,6),0,0,Fraction(2,3),Fraction(1,6)]\n", "A=[[0,0,0,0,0],[Fraction(1,3),0,0,0,0],[Fraction(1,6),Fraction(1,6),0,0,0],[Fraction(1,8),0,Fraction(3,8),0,0],[Fraction(1,2),0,-Fraction(3,2),2,0]]\n", "s=len(c)\n", "\n", "print(\"Consistance\")\n", "print(sum(b)==1)\n", "for i in range(s):\n", " print(sum(A[i])==c[i])\n", " \n", "print(\"\\nOrdre 2\")\n", "print(sum([b[i]*c[i] for i in range(s)])==Fraction(1,2))\n", "\n", "print(\"\\nOrdre 3\")\n", "print(sum([b[i]*c[i]**2 for i in range(s)])==Fraction(1,3))\n", "print(sum([b[i]*A[i][j]*c[j] for i in range(s) for j in range(s)])==Fraction(1,6))\n", "\n", "print(\"\\nOrdre 4\")\n", "print(sum([b[i]*c[i]**3 for i in range(s)])==Fraction(1,4))\n", "print(sum([b[i]*c[i]*A[i][j]*c[j] for i in range(s) for j in range(s)])==Fraction(1,8))\n", "print(sum([b[i]*A[i][j]*c[j]**2 for i in range(s) for j in range(s)])==Fraction(1,12))\n", "print(sum([b[i]*A[i][j]*A[j][k]*c[k] for i in range(s) for j in range(s) for k in range(s)])==Fraction(1,24))\n", "\n", "print(\"\\nOrdre 5\")\n", "print(sum([b[i]*c[i]**4 for i in range(s)])==Fraction(1,5))\n", "print(sum([b[i]*c[i]**2*(1-c[i]) for i in range(s)])+sum([b[i]*A[i][j]*c[j]*(2*c[i]-c[j]) for i in range(s) for j in range(s)])==Fraction(1,4))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exemple à 6 étages explicite\n", "\n", "Un schéma RK à $6$ étages explicite a pour matrice de Butcher\n", "$$\n", "\\begin{array}{c|cccccc}\n", "0 & 0 & 0 & 0 & 0 &0&0 \\\\\n", "c_2 & a_{21} & 0&0&0&0&0\\\\\n", "c_3 & a_{31} &a_{32}&0&0&0&0\\\\\n", "c_4 & a_{41} &a_{42}&a_{43}&0&0&0\\\\\n", "c_5 & a_{51} &a_{52}&a_{53}&a_{54}&0&0\\\\\n", "c_6 & a_{61} &a_{62}&a_{63}&a_{64}&a_{65}&0\\\\\n", "\\hline\n", " & b_1 & b_2 & b_3 & b_4 & b_5&b_6\n", "\\end{array}\n", "\\quad\n", "\\implies\n", "\\begin{cases}\n", "u_0\t = y_0 \\\\\n", "K_1 = \\varphi\\left(t_n,u_n\\right)\\\\\n", "K_2 = \\varphi\\left(t_n+hc_2,u_n+h(a_{21}K_1)\\right)\\\\\n", "K_3 = \\varphi\\left(t_n+hc_3,u_n+h(a_{31}K_1+ a_{32}K_2\\right)\\\\\n", "K_4 = \\varphi\\left(t_n+hc_4,u_n+h(a_{41}K_1+ a_{42}K_2+ a_{43}K_3)\\right)\\\\\n", "K_5 = \\varphi\\left(t_n+hc_5,u_n+h(a_{51}K_1+ a_{52}K_2+ a_{53}K_3+ a_{54}K_4)\\right)\\\\\n", "K_6 = \\varphi\\left(t_n+hc_6,u_n+h(a_{61}K_1+ a_{62}K_2+ a_{63}K_3+ a_{64}K_4+ a_{65}K_5)\\right)\\\\\n", "u_{n+1} = u_n + h\\left(b_1K_1+b_2K_2+b_3K_3+b_4K_4+b_5K_5+b_6K_6\\right) & n=0,1,\\dots N-1\n", "\\end{cases}\n", "$$\n", "avec \n", "$$\n", "\\begin{cases}\n", "c_2=a_{21}\\\\\n", "c_3=a_{31}+a_{32}\\\\\n", "c_4=a_{41}+a_{42}+a_{43}\\\\\n", "c_5=a_{51}+a_{52}+a_{53}+a_{54}\\\\\n", "c_6=a_{61}+a_{62}+a_{63}+a_{64}+a_{65}\\\\\n", "b_1+b_2+b_3+b_4+b_5+b_6=1\n", "\\end{cases}\n", "$$\n", "\n", "Voici un exemple:\n", "+ schéma de Butcher (d'ordre 5)\n", "$$\n", "\\begin{array}{c|cccccc}\n", "0 & 0 & 0 & 0 & 0 &0&0 \\\\\n", "\\frac{1}{4} & \\frac{1}{4} & 0&0&0&0&0\\\\\n", "\\frac{1}{4} & \\frac{1}{8} &\\frac{1}{8}&0&0&0&0\\\\\n", "\\frac{1}{2} & 0 &-\\frac{1}{2}&1&0&0&0\\\\\n", "\\frac{3}{4} & \\frac{3}{16} &0&0&\\frac{9}{16}&0&0\\\\\n", "1 & \\frac{-3}{7} &\\frac{2}{7}&\\frac{12}{7}&\\frac{-12}{7}&\\frac{8}{7}&0\\\\\n", "\\hline\n", " & \\frac{7}{90} & 0&\\frac{32}{90} & \\frac{12}{90} & \\frac{32}{90} & \\frac{7}{90}\n", "\\end{array}\n", "\\quad\n", "\\implies\n", "\\begin{cases}\n", "u_0\t = y_0 \\\\\n", "K_1 = \\varphi\\left(t_n,u_n\\right)\\\\\n", "K_2 = \\varphi\\left(t_n+\\frac{h}{4},u_n+\\frac{h}{4}K_1\\right)\\\\\n", "K_3 = \\varphi\\left(t_n+\\frac{h}{4},u_n+\\frac{h}{8}(K_1+K_2)\\right)\\\\\n", "K_4 = \\varphi\\left(t_n+\\frac{h}{2},u_n+\\frac{h}{2}(-K_2+2K_3)\\right)\\\\\n", "K_5 = \\varphi\\left(t_n+\\frac{3h}{4},u_n+\\frac{h}{16}(3K_1+9K_4)\\right)\\\\\n", "K_6 = \\varphi\\left(t_{n+1},u_n+\\frac{h}{7}(-3K_1+2K_2+12K_3-12K_4+8K_5)\\right)\\\\\n", "u_{n+1} = u_n + \\frac{h}{90}\\left(7K_1+32K_3+12K_4+32K_5+7K_6\\right) & n=0,1,\\dots N-1\n", "\\end{cases}\n", "$$" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Consistance\n", "True\n", "True\n", "True\n", "True\n", "True\n", "True\n", "True\n", "\n", "Ordre 2\n", "True\n", "\n", "Ordre 3\n", "True\n", "True\n", "\n", "Ordre 4\n", "True\n", "True\n", "True\n", "True\n", "\n", "Ordre 5\n", "True\n", "True\n" ] } ], "source": [ "from fractions import Fraction\n", "c=[0,Fraction(1,4),Fraction(1,4),Fraction(1,2),Fraction(3,4),1]\n", "b=[Fraction(7,90),0,Fraction(32,90),Fraction(12,90),Fraction(32,90),Fraction(7,90)]\n", "A=[[0,0,0,0,0,0],\n", " [Fraction(1,4),0,0,0,0,0],\n", " [Fraction(1,8),Fraction(1,8),0,0,0,0],\n", " [0,Fraction(-1,2),1,0,0,0],\n", " [Fraction(3,16),0,0,Fraction(9,16),0,0],\n", " [Fraction(-3,7),Fraction(2,7),Fraction(12,7),Fraction(-12,7),Fraction(8,7),0]]\n", "s=len(c)\n", "\n", "print(\"Consistance\")\n", "print(sum(b)==1)\n", "for i in range(s):\n", " print(sum(A[i])==c[i])\n", " \n", "print(\"\\nOrdre 2\")\n", "print(sum([b[i]*c[i] for i in range(s)])==Fraction(1,2))\n", "\n", "print(\"\\nOrdre 3\")\n", "print(sum([b[i]*c[i]**2 for i in range(s)])==Fraction(1,3))\n", "print(sum([b[i]*A[i][j]*c[j] for i in range(s) for j in range(s)])==Fraction(1,6))\n", "\n", "print(\"\\nOrdre 4\")\n", "print(sum([b[i]*c[i]**3 for i in range(s)])==Fraction(1,4))\n", "print(sum([b[i]*c[i]*A[i][j]*c[j] for i in range(s) for j in range(s)])==Fraction(1,8))\n", "print(sum([b[i]*A[i][j]*c[j]**2 for i in range(s) for j in range(s)])==Fraction(1,12))\n", "print(sum([b[i]*A[i][j]*A[j][k]*c[k] for i in range(s) for j in range(s) for k in range(s)])==Fraction(1,24))\n", "\n", "print(\"\\nOrdre 5\")\n", "print(sum([b[i]*c[i]**4 for i in range(s)])==Fraction(1,5))\n", "print(sum([b[i]*c[i]**2*(1-c[i]) for i in range(s)])+sum([b[i]*A[i][j]*c[j]*(2*c[i]-c[j]) for i in range(s) for j in range(s)])==Fraction(1,4))" ] } ], "metadata": { "colab": { "collapsed_sections": [], "default_view": {}, "name": "EdoExplicites.ipynb", "provenance": [], "version": "0.3.2", "views": {} }, "hide_input": false, "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.10.6" }, "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": true, "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": "", "title_sidebar": "Contents", "toc_cell": true, "toc_position": { "height": "calc(100% - 180px)", "left": "10px", "top": "150px", "width": "307px" }, "toc_section_display": true, "toc_window_display": false } }, "nbformat": 4, "nbformat_minor": 4 }