{ "cells": [ { "cell_type": "code", "execution_count": 101, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "execution_count": 101, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from IPython.core.display import HTML, display\n", "css_file = './custom.css'\n", "HTML(open(css_file, \"r\").read())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# 2021 CR : le sujet comporte trois exercices indépendents." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercice 1 : étude d'un schéma\n", "\n", "On étudie le schéma de résolution du problème de Cauchy $y'(t)=\\varphi(t,y(t))$ et $y(t_0)=y_0$ défini par la formule de récurrence\n", "\n", "$$\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{1}{2}h,u_n+\\frac{1}{2}hK_1\\right)\\\\\n", "K_3 = \\varphi\\left(t_{n}+h,u_n+hK_1\\right)\\\\\n", "u_{n+1} = u_n + h\\left(aK_1+bK_2+cK_3\\right) & n=0,1,\\dots N-1\n", "\\end{cases}\n", "$$\n", "où $a,b,c \\in [0;1]$.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ">**Q1 [3 points]** \n", "Pour quelles valeurs du triplet $(a,b,c)$ retrouve-t-on \n", "- la méthode d'Euler explicite (EE) ?\n", "- la méthode d'Euler modifiée (EM) ?\n", "- la méthode de Heun (H) ?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- Si $(a,b,c)=(1,0,0)$ on retrouve la méthode d'Euler explicite (EE):\n", "$$\n", "u_{n+1} = u_n + h K_1 = u_n + h \\varphi\\left(t_n,u_n\\right)\n", "$$\n", "- Si $(a,b,c)=(0,1,0)$ on retrouve la méthode d'Euler modifiée (EM):\n", "$$\n", "u_{n+1} = u_n + h K_2 = u_n + h \\varphi\\left(t_{n}+\\frac{1}{2}h,u_n+\\frac{1}{2}h\\varphi\\left(t_n,u_n\\right)\\right)\n", "$$\n", "- Si $(a,b,c)=(\\frac{1}{2},0,\\frac{1}{2})$ on retrouve la méthode de Heun (H):\n", "$$\n", "u_{n+1} = u_n + \\frac{h}{2} (K_1+K_3) = u_n + \\frac{h}{2} \\left(\\varphi\\left(t_n,u_n\\right)+\\varphi\\left(t_{n}+h,u_n+h\\varphi\\left(t_n,u_n\\right)\\right)\\right)\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ">**Q2 [2 points]** \n", "Le schéma peut être vu comme une méthode de type Runge-Kutta. Écrire la matrice de Butcher associée. \n", "Le schéma est explicite, semi-implicite ou implicite ?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{array}{c|ccc} \n", "0 & 0 & 0 & 0 \\\\ \n", "\\dfrac{1}{2} & \\dfrac{1}{2} & 0 & 0 \\\\ \n", "1 & 1 & 0 & 0 \\\\ \n", "\\hline \n", " & a & b & c \n", "\\end{array}\n", "$$\n", "Le schéma est explicite (la matrice $\\mathbb{A}$ est triangulaire inférieure à diagonale nulle). " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ">**Q3 [1 points]** \n", "Sans faire de calculs, indiquer quel est l'ordre maximale du schéma." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Soit $\\omega$ l'ordre de la méthode et $s=3$ le nombre d'étages. Étant un schéma explicite, $\\omega\\le s=3$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ">**Q4 [3 points]** \n", "Étudier théoriquement l'ordre du schéma en fonction de $a,b,c$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "+ Si $\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", "alors $\\omega\\ge1$\n", "\n", "+ Si $\\omega\\ge1$ et $\\displaystyle\\sum_{j=1}^s b_j c_j=\\frac{1}{2}$\n", "alors $\\omega\\ge2$\n", "\n", "+ Si $\\omega\\ge2$ et $\\begin{cases}\\displaystyle\\sum_{j=1}^s b_j c_j^2=\\frac{1}{3}\\\\\\displaystyle\\sum_{i=1}^s\\sum_{j=1}^s b_i a_{ij} c_j=\\frac{1}{6}\\end{cases}$ alors $\\omega\\ge3$ \n", "\n", "Il est donc\n", "- au moins d'ordre 1 ssi $a=1-b-c$ \n", "- il est d'ordre 2 ssi, de plus, $\\frac{b}{2}+c=\\frac{1}{2}$ donc ssi $\\begin{cases}a=c\\\\b=1-2c\\end{cases}$ \n", "- il ne peut pas être d'ordre 3 car la deuxième condition n'est jamais satisfaite.\n", "\n", "On peut vérifier nos calculs comme ci-dessous:" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Étude de la consistance (ordre 1)\n" ] }, { "data": { "text/latex": [ "$\\displaystyle \\sum_{j=1}^s b_j=a + b + c\\text{ doit être egale à 1}$" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/latex": [ "$\\displaystyle i={}0\\quad \\sum_{j=1}^s a_{ij}-c_i=0\\text{ doit être egale à 0}$" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/latex": [ "$\\displaystyle i={}1\\quad \\sum_{j=1}^s a_{ij}-c_i=0\\text{ doit être egale à 0}$" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/latex": [ "$\\displaystyle i={}2\\quad \\sum_{j=1}^s a_{ij}-c_i=0\\text{ doit être egale à 0}$" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "\n", "Étude de l'Ordre 2\n", "\n" ] }, { "data": { "text/latex": [ "$\\displaystyle \\sum_{j=1}^s b_j c_j=\\frac{b}{2} + c\\text{ doit être egale à }\\frac{1}{2}$" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "\n", "Étude de l'Ordre 3\n", "\n" ] }, { "data": { "text/latex": [ "$\\displaystyle \\sum_{j=1}^s b_jc_j^2=\\frac{b}{4} + c\\text{ doit être egale à }\\frac{1}{3}$" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/latex": [ "$\\displaystyle \\sum_{i=1}^s\\sum_{j=1}^s b_ia_{ij}c_j=0\\text{ doit être egale à }\\frac{1}{6}$" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import sympy as sym\n", "sym.init_printing()\n", "\n", "from IPython.display import display, Math\n", "\n", "sym.var('a,b,c')\n", "\n", "cc=[0,sym.S(1)/2,1]\n", "bb=[a,b,c]\n", "A=[[cc[0],0,0],[cc[1],0,0],[cc[2],0,0]]\n", "s=len(cc)\n", "print('Étude de la consistance (ordre 1)')\n", "calc=[sum(bb)]\n", "display(Math(r'\\sum_{j=1}^s b_j='+sym.latex(calc[0]) +'\\\\text{ doit être egale à 1}' ))\n", "for i in range(s):\n", " calc.append(sum(A[i])-cc[i]) \n", " display(Math(r'i={}'+str(i)+'\\quad \\sum_{j=1}^s a_{ij}-c_i='+sym.latex(calc[-1])+'\\\\text{ doit être egale à 0}'))\n", "\n", "print(\"\\nÉtude de l'Ordre 2\\n\")\n", "calc=sum([bb[i]*cc[i] for i in range(s)]) \n", "display(Math(r'\\sum_{j=1}^s b_j c_j='+sym.latex(calc )+'\\\\text{ doit être egale à }\\\\frac{1}{2}')) \n", "\n", "print(\"\\nÉtude de l'Ordre 3\\n\") \n", "calc3_1=sum([bb[i]*cc[i]**2 for i in range(s)])\n", "display(Math(r'\\sum_{j=1}^s b_jc_j^2='+sym.latex(calc3_1)+'\\\\text{ doit être egale à }\\\\frac{1}{3}' ))\n", "calc3_2=sum([bb[i]*A[i][j]*cc[j] for i in range(s) for j in range(s)])\n", "display(Math(r'\\sum_{i=1}^s\\sum_{j=1}^s b_ia_{ij}c_j='+sym.latex(calc3_2)+'\\\\text{ doit être egale à }\\\\frac{1}{6}' ))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ">**Q5 [2 points]** \n", "Étudier théoriquement la A-stabilité du schéma pour les valeurs du triplet $(a,b,c)$ qui donnent l'ordre maximal. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Soit $\\beta>0$ un nombre réel positif et considérons le problème de Cauchy\n", "$$\\begin{cases}\n", "y'(t)=-\\beta y(t), &\\text{pour }t>0,\\\\\n", "y(0)=1.\n", "\\end{cases}$$\n", "Sa solution est $y(t)=e^{-\\beta t}$ donc $$\\lim_{t\\to+\\infty}|y(t)|=0.$$\n", "\n", "Le schéma appliqué à ce problème de Cauchy s'écrit\n", "$$\n", "\\begin{cases}\n", "u_0\t = y_0 \\\\\n", "K_1 = -\\beta u_n\\\\\n", "K_2 = -\\beta \\left(u_n+\\frac{1}{2}hK_1\\right)=-\\beta \\left(u_n-\\frac{1}{2}\\beta h u_n\\right)\\\\\n", "K_3 = -\\beta \\left(u_n+hK_1\\right)=-\\beta \\left(u_n-\\beta h u_n\\right)\\\\\n", "u_{n+1} = u_n + h\\left(aK_1+bK_2+cK_3\\right) & n=0,1,\\dots N-1\n", "\\end{cases}\n", "$$\n", "d'où\n", "$$\n", "u_{n+1} = \\left( 1-(a+b+c)\\beta h + \\left(\\frac{b}{2}+c\\right)(\\beta h)^2 \\right)u_n\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "L'ordre maximal (i.e. l'ordre 2) est obtenu pour $a+b+c=1$ et $b=1-2c$, on trouve alors \n", "$$\n", "u_{n+1} \n", "= \n", "\\left(\\dfrac{1}{2}(\\beta h)^2 -(\\beta h)+1\\right)u_n\n", "$$\n", "Vérifions ces calculs:" ] }, { "cell_type": "code", "execution_count": 41, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$\\displaystyle K_1=- \\beta u_{n}$" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/latex": [ "$\\displaystyle K_2=\\frac{\\beta u_{n} \\left(\\beta dt - 2\\right)}{2}$" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/latex": [ "$\\displaystyle K_3=\\beta u_{n} \\left(\\beta dt - 1\\right)$" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/latex": [ "$\\displaystyle u_{n+1}=- \\frac{u_{n} \\left(2 a \\beta dt - b \\beta^{2} dt^{2} + 2 b \\beta dt - 2 \\beta^{2} c dt^{2} + 2 \\beta c dt - 2\\right)}{2}$" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "Avec les paramètres qui donnent l'ordre 2 on obtient\n" ] }, { "data": { "text/latex": [ "$\\displaystyle u_{n+1}=\\frac{u_{n} \\left(\\beta^{2} dt^{2} - 2 \\beta dt + 2\\right)}{2}$" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import sympy as sym\n", "sym.init_printing()\n", "from IPython.display import display, Math\n", "\n", "# pour ne pas effacer l'affectation de \"h\", ici je vais l'appeler \"dt\"\n", "sym.var('u_n,u_np1,beta,dt,x, a,b,c')\n", "\n", "K1 = -beta*u_n\n", "display(Math('K_1='+sym.latex(K1)))\n", "K2 = -beta*(u_n+sym.S(1)/2*dt*K1).factor()\n", "display(Math('K_2='+sym.latex(K2)))\n", "K3 = -beta*(u_n+dt*K1).factor()\n", "display(Math('K_3='+sym.latex(K3)))\n", "u_np1 = (u_n+dt*(a*K1+b*K2+c*K3)).factor()\n", "display(Math('u_{n+1}='+sym.latex(u_np1)))\n", "print(\"Avec les paramètres qui donnent l'ordre 2 on obtient\")\n", "display(Math('u_{n+1}='+sym.latex(u_np1.subs(a,1-b-c).subs(b,1-2*c).simplify())))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On note $x=\\beta h$ et on étudie la fonction $q\\colon \\mathbb{R}^+\\to\\mathbb{R}$ définie par $q(x)=\\frac{1}{2}x^2-x+1$ car le schéma est A-stable ssi $|q(x)|<1$. \n", "Il s'agit d'une parabole convexe de sommet en $x_s=1$ et $q(x_s)=\\frac{1}{2}>-1$ donc $q(x)<1$ pour $00$ pour tout $x$.\n", "On conclut que le schéma est A-stable pour tout $h<\\dfrac{2}{\\beta}$.\n", "\n", "On peut vérifier les calculs avec `sympy`:" ] }, { "cell_type": "code", "execution_count": 89, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$\\displaystyle q(x)=\\frac{x^{2}}{2} - x + 1$" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/latex": [ "$\\displaystyle q(0)=1$" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/latex": [ "$\\displaystyle \\lim_{x \\to \\infty}\\left(\\frac{x^{2}}{2} - x + 1\\right)=\\infty$" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "On cherche les points stationnaires dans R^+\n" ] }, { "data": { "text/latex": [ "$\\displaystyle q'(x)=x - 1$" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "q est une parabole dont le sommet est\n" ] }, { "data": { "text/latex": [ "$\\displaystyle q'(x)=0 \\iff x\\in\\left[ 1\\right]$" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/latex": [ "$\\displaystyle x=1\\text{ est un minimum et}$" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/latex": [ "$\\displaystyle q(1)=\\frac{1}{2}$" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaUAAAEYCAYAAAD8hukFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAd8klEQVR4nO3de3SV5YHv8d+bbBIIuRAgIQkbCBAIIZAESLwiR0EuDTRW9ChWRxQtwuiyl2ln2WnHU20d0aOrtWqlGasD1YKXsUARqGhhigjEyE1AJUACCQnkToAQyOU5fyA5UkISIHu/7977+1krS3f2k/f5ZS3dv7y357WMMQIAwAmC7A4AAMA5lBIAwDEoJQCAY1BKAADHoJQAAI5BKQEAHINSAgA4BqUEAHAMV2cGWZYVK+l6SQmSTknaJSnfGNPiwWwAgABjtbeig2VZN0l6TFJvSdsklUvqLmm4pKGS3pX0vDGmzvNRAQD+rqNS+r+SXjTGHGrjPZekGZKCjTH/7bmIAIBA0W4peXxyy3pNZ4ut3Bgzqo33b5S0XFLh1996zxjzpPcSAgC8qVMXOliW9UfLsqK+8TrRsqyPumD+/5I0rYMxG4wxGV9/UUgA4Mc6e/Xdx5K2WJaVbVnW9yR9IOk3Vzq5MebvkqqvdDsAAP/Q6cN3lmWNl7ROUqWkMcaYI1cwb+ukRUVFmjFjhnbt2nXBoPXr1+u2225T/4GDdHzyLzRzdLSe/6cb2txgbm6ucnNzJUmnTp3S7t27ryCes/1ixW69nV+sXb+YqqAgy+44ANBZHX5gdfbw3T9Jek3SvTp7yG2VZVnpVxStE8aOHauDBw9q57atGhAZrPfW5V907Ny5c5Wfn6/8/Hz16NHD09Fstae0TinxkRQSAL/T2cN3t0kab4xZYoz5qaR5OltOHhUZGanw8HBJ0oRRA2V6DVB5RYWnp3U0Y4y+KKtTSnyE3VEAoMt1qpSMMd8xxpR/43WepKs9luprR44c0bnDixFnaqSQHqoz3T09raOV1JzS8dNNGhkf1fFgAPAx7a7oYFnWzyX9zhhzwcUIxpgzlmVNlBRmjFl5OZPfddddWr9+vSorK+V2u/XEE0+osbFRkjRv3jy9++67euWVV+RyueTqO0jKekifHz6mpNjA3UvYXXr2PuWRCZE2JwGArtfRMkOfS/qLZVkNkrZKqtDZFR2GScqQ9KGk/7jcyZcsWdLu+4888ogeeeQRSVJTc4tG/+ID7Sg+plvHuC93Sp93sPqk3L16KLlf4BYzAP/VUSndboy53rKsf9XZJYbiJdVJekPSXGPMKU8HPMcVHKRR/SP1+eFj3prSkT7ZV6Ww0GD1CAm2OwoAdLmOSmmcZVmDJN0t6aZ/eK+Hzi7O6jWj+/fSn/IOqqm5Ra7gwFvg3BijHSW1mjoyzu4oAOARHX2yL5S0RtIISfnf+Prs6396VfqAKDU0tmjv0RPentoRiqrqVVvfqIyBveyOAgAe0W4pGWN+a4xJkfSaMWbIN74GG2OGeCljqzT32Q/jzw/XentqR9heXCNJyhhAKQHwT529JHy+p4N0RmKfMEV0d2lHSWCeV9p+qFZhIcEazkUOAPyUT52YsSxLae4o7SwJzD2lbcW1SnNHKZiVHAD4KZ8qJensIbwvy46robHZ7ihe1dDYrC/K6pQxINruKADgMT5XSunuKDW1GH155LjdUbxqd2mdGpsN55MA+DWfK6XRX1/sEGiH8LYXn/19x3DlHQA/5nOllBDVXTcmx6ikxqu3SNlue3GtEqK6q19kYK/9B8C/+VwpWZalkOAgfbD7Sh7n5Hu2F9dwfxIAv+dzpSRJmYnRKqqqV8Xx03ZH8YrKE6dVXH2K80kA/J5PltK4QWevQPvsYI3NSbxj+6Gz55O48g6Av/PJUhrVP0ohriB9dvCCJ2r4pe3FtQoOsjS6P89QAuDffLKUQl3BSusfpfxA2VMqrtWIuAhWBgfg93yylCRpXGK0dh0+5vc30ba0GO0oruV8EoCA4LOllDmotxqbjXb6+Tp4BeXHldi3p64d2sfuKADgcT5bSucudsj38/NKWwqr9fnhY0p3s6cEwP/5bCn17hmiITE9tdXPzyttKaxWfFR3uaN72B0FADzOZ0tJkjIHReuzgzUyxtgdxSOMMdpyoFpXD+4ty2JlcAD+z8dLqbdq6hu1v+Kk3VE8orDypCpPnNbVQzifBCAw+HQpjW29idY/zyttKTz7e101uLfNSQDAO3y6lIbG9FR0WDflF/nneaUtB6rUNzxUQ/r2tDsKAHiFT5eSZVkaNyhah6rr7Y7S5Ywx2lLI+SQAgcWnS0mSrhnSR1sKq3XkWIPdUbpUSc0plR1r0NVDOHQHIHD4RSlJ0sZ9lTYn6VqcTwIQiHy+lEbGRyo6rJs27vezUjpQpV5h3TQ8NsLuKADgNT5fSkFBlq4b2lef7Kvyq/uV8oqqlZXYW0FBnE8CEDh8vpQk6bqkPjpS16ADlf5xv9KRYw06WFWvqzl0ByDA+EUpXT+0ryTpEz85r7SlsErS/z9fBgCBwi9KaVCfMPXv1UMb91XZHaVLbCmsVkSoSynxkXZHAQCv8otSsixL1w3to0/2V6q5xffPKx2qOqnvjElQMOeTAAQYvyglSbo+qa/qGpq0u9S3n69UXF2vj/dVaXDfcLujAIDX+U0pXTf03P1Kvn0I7+8FFZKkCcP72pwEALzP1lKaM2eOYmNjNWrUqDbfN8bo0UcfVVJSktLS0rR169aLbis2sruG9wvXJz5+v9KGvZVKiOquoTHsKQEIPLaW0n333ac1a9Zc9P3Vq1eroKBABQUFys3N1fz589vd3nVD++rTomr56u1KTc0t2ri/UjcMi2G9OwABybLjhtM7f7+pddKGhgZ9vutzZWVmXTBu79696tWrl2JjYyVJeXl5ysjIUEhISJvbrak/o71HT6jmT/+iaT973UPpPefE6SbtLq3TsNhw9e7Z9u8I+LO3HrrW7gjwrA7/2nb0OaXTZ04rNDS09XVoaKhOnz7d5tiysjId+HK3ZIwaG5u8FbFL1dY3SpIie3SzOQkA2MOWPSVJrZMWFRVpxowZ2rVr1wWDpk+frp/+9KcaP368JGnSpEl69tlnNW7cuItu+IdvbdebP/2uyg/s8UBsz7r1dxtljLTs4evtjgIAnuDbe0put1vFxcWtr0tKSpSQkNDuzwzp21PHTjWq/LhvPcriWH2jdhTXasLwGLujAIBtHF1KOTk5Wrx4sYwx2rx5s6KiohQfH9/uz0xK6SdJWvdluTcidpmP91UoM7G3bkqmlAAELpedk991111av369Kisr5Xa79cQTT6ix8ex5lXnz5ik7O1urVq1SUlKSwsLC9PrrHV+8kBIfoW7BQfrwi3LdmTXQ079Cl/nr7qPaX35Cae5edkcBANvYWkpLlixp933LsvTyyy9f0jYty1Jkd5c2FFSoobFZ3bsFX0lErzjT1KJ1X5brW6PjWFoIQEBz9OG7yxXRvZsaGlt85kbaLYVVOn66SVNGxtkdBQBs5ZelFB7qUs+QYH34hW+cV/pg91H16Bas8cNYWghAYPPLUrIs6YZhMfrbF+WOfxqtMUZr9xzVhOF9feJQIwB4kl+WkiRNSonVkboG7S6tsztKuz4/fExH6ho4dAcA8uNSumlErKLDummjw59G+8HuowoOsjRxRKzdUQDAdn5bSn3DQzUiLlJv5Rc7+hDe2j1HlZUYrWjWugMA/y0lSZqRHq8DFScdewivqPKEvjp6nEN3APA1vy6l7FHxcgVZ+suOUrujtOn9nUeUEhehKan97I4CAI7g16UU3TNENwzrq7/sKFVLi/MO4f1lZ6nCQl1yR4fZHQUAHMGvS0mScjISVHqsQVsP1dgd5TxfHTmuL48cV056+wvMAkAg8ftSmjwyTqGuIK1w2CG8FTsOKzjIUvbo9heYBYBA4velFB7q0s0p/bTq8zI1NbfYHUfS2Rtml28v1fVJfRUTEdrxDwBAgPD7UpKkb6cnqPLEGX2yv8ruKJKkT4uqFR7q0p2ZbrujAICjBEQp3Zgco+uG9tG6r5yxFt5bn5aopOaUbuKGWQA4T0CUUvduwRoaE643txxSbf0ZW7Mcb2jUqs/L9O30eIWF2PrkEABwnIAoJUm666qBOtPUove2HrY1x8qdZTrV2Kw7MgfYmgMAnChgSmlkQqQyBvTSkrxDti07ZIzRpv1VunZIH2UM4AmzAPCPAqaUJOm7Vw1UQfkJ5R+0556l/IM1WrGjVNlp8bIsnjALAP8ooEppRnq8rhnc27Z7ll7dcEC9wrrp9rFcdQcAbQmoUgoLcWm0O0p/2nJIJTX1Xp37YNVJfbDnqO65epB6hPAwPwBoS0CVkiTdf/1gWZL+8HGhV+d9fWORXEGW7r12kFfnBQBfEnCllNCrh3LSE/TWp8U6Vt/olTmP1Tfq7fxi5aT3V2xkd6/MCQC+KOBKSZK+N2GI6s80640tRV6Z7095B1V/plkPjB/slfkAwFcFZCmlxEfq3msH6b82HlRdg2f3lurPNGnZ9lLdffVAjUyI9OhcAODrArKUJOmOzAGqOHFav/+f/R6d5/WNRfrqyHHN5Io7AOhQwJbSqP5RuiUjQX/4uFBH6xo8Mkdt/Rkt/J/9ujklVuMGRXtkDgDwJwFbSpL0L5OT1dxi9JsP93pk+6+s36+R8ZH68dRkj2wfAPxNQJfSwD5hmn9jkrYX12pHcW2Xbnt36TH954YDGtYvXCPiOJcEAJ0R0KUkSQ/eMFi19Y368Ts7dLqpuUu22dTUon9ftkvRYSH6yZQRXbJNAAgEAV9Kkd276emZo1VQfkK//aigS7b5nx8fUFOL0f/59khFhXXrkm0CQCAI+FKSpBuTY3VHpluvrN+vzw5WX9G2th6q0fMf7JU7uoe+nZ7QRQkBIDBQSl/72fSRmjyyn/75za06cuzyrsarPHFa//ruTsVFddfTM9NYCRwALhGl9LWoHt30g5uH60RDk+5/PU/HLvEJtccbGnX/659KMnrl7rGK6sFhOwC4VJTSN6TER+qlu8eqd88Q3fOHLSo/3rk9pmP1jXrsv3fqi7I6/Vt2ika7eYAfAFwOW0tpzZo1Sk5OVlJSkhYsWHDB++vXr1dUVJQyMjKUkZGhJ5980uOZbkqO1Zzxg7Wv/KRm/m6jPi2qanf8juJa3fq7jdq4r1Kv3DNOE0f083hGAPBXLrsmbm5u1sMPP6y1a9fK7XYrKytLOTk5Gjly5HnjbrjhBq1cudKr2Sal9NNbD12jhev3638v3Kyc9Hjde22ixg6MVlCQpZYWo52Hj+m9z0r0xpaDmpTST89MSFNWYm+v5gQAf2NbKeXl5SkpKUlDhgyRJM2aNUvLly+/oJTskubupWduT9PgmJ768Iujun3hJg2PDdfp5hbFhIcq/2CNxif11b3XJuqHNw/n0m8A6AK2ldLhw4c1YMCA1tdut1tbtmy5YNymTZuUnp6uhIQEPffcc0pNTW1ze7m5ucrNzZUkVVRUdEnGiO7d9JOpI/S9G4Zo3VflKq09pb1HT2hEXITuuWaQ/tfwGEX3DOmSuQAANpaSMeaC7/3jJdRjx47VwYMHFR4erlWrVuk73/mOCgravsF17ty5mjt3riQpMzOzS7P2CgvRrWNY5RsAPM22Cx3cbreKi4tbX5eUlCgh4fybTSMjIxUeHi5Jys7OVmNjoyorK72aEwDgPbaVUlZWlgoKClRYWKgzZ85o6dKlysnJOW/MkSNHWveo8vLy1NLSoj59+tgRFwDgBbYdvnO5XHrppZc0depUNTc3a86cOUpNTdXChQslSfPmzdO7776rV155RS6XSz169NDSpUtZJQEA/JjV1rkdL/DopJmZmcrPz/fkFACAS9fhXgUrOgAAHINSAgA4BqUEAHAMSgkA4BiUEgDAMSglAIBjUEoAAMeglAAAjkEpAQAcg1ICADgGpQQAcAxKCQDgGJQSAMAxKCUAgGNQSgAAx6CUAACOQSkBAByDUgIAOAalBABwDEoJAOAYlBIAwDEoJQCAY1BKAADHoJQAAI5BKQEAHINSAgA4BqUEAHAMSgkA4BiUEgDAMSglAIBjUEoAAMeglAAAjkEpAQAcg1ICADiGraW0Zs0aJScnKykpSQsWLLjgfWOMHn30USUlJSktLU1bt261ISUAwFtsK6Xm5mY9/PDDWr16tfbs2aMlS5Zoz549541ZvXq1CgoKVFBQoNzcXM2fP9+mtAAAb7CtlPLy8pSUlKQhQ4YoJCREs2bN0vLly88bs3z5ct17772yLEvXXHONamtrVVZWZlNiAICnWcYYr086bdo0U1BQoLq6Og0aNEiSVFVVpZMnT2rgwIGt4/bt26e4uDiFh4dLkvbu3Su3262wsLALtllRUaHKykpJ0unTp5WRkeGF36TrVVRUKCYmxu4Yl4389iK/vcjfvs8+++yvxphp7Q4yxtjxZd5++23zwAMPmHMWL15sHnnkEfNN2dnZZsOGDa2vJ06caPLz801HwsLCOhzjVOPGjbM7whUhv73Iby/yd6jDfrDt8J3b7VZxcXHr65KSEiUkJFzyGACA/7CtlLKyslRQUKDCwkKdOXNGS5cuVU5OznljcnJytHjxYhljtHnzZkVFRSk+Pt6mxAAAT3PZNrHLpZdeeklTp05Vc3Oz5syZo9TUVC1cuFCSNG/ePGVnZ2vVqlVKSkpSWFiYXn/99U5tu2/fvp6M7lFz5861O8IVIb+9yG8v8l85Wy50kOTRSTMzM5Wfn+/JKQAAl87qaAArOgAAHINSAgA4hl+V0rlli3bt2tXmskVOVlxcrJtuukkpKSlKTU3VCy+8YHeky9Lc3KwxY8ZoxowZdke5ZLW1tbr99ts1YsQIpaSkaNOmTXZHuiS//vWvlZqaqlGjRumuu+5SQ0OD3ZHaNWfOHMXGxmrUqFGt36uurtbkyZM1bNgwTZ48WTU1NTYmbF9b+X/yk59oxIgRSktL06233qra2lobE15cW9nPee6552RZVut9n97mN6X0zWWLUlNT21y2yMlcLpeef/55ffHFF9q8ebNefvlln8p/zgsvvKCUlBS7Y1yW73//+5o2bZq+/PJL7dixw6d+j8OHD+u3v/2t8vPztWvXLjU3N2vp0qV2x2rXfffdpzVr1pz3vQULFmjSpEkqKCjQpEmTHP3HZVv5J0+erF27dmnnzp0aPny4nn76aZvSta+t7NLZP47Xrl173iIG3uY3pfTNZYssy2pz2SIni4+P19ixYyVJERERSklJ0eHDh21OdWlKSkr0/vvv68EHH7Q7yiWrq6vT3//+dz3wwAOSpJCQEPXq1cvmVJemqalJp06dUlNTk+rr6x1/T9+ECRPUu3fv8763fPlyzZ49W5I0e/ZsLVu2zI5ondJW/ilTpsjlOntR8zXXXKOSkhI7onWoreyS9MMf/lDPPvusLKvD6xE8xm9K6fDhwxowYEDra7fb7XMf6ucUFRVp27Ztuvrqq+2Ockl+8IMf6Nlnn1VQkO/9Z3XgwAHFxMTo/vvv15gxY/Tggw/q5MmTdsfqtP79++vHP/6xBg4cqPj4eEVFRWnKlCl2x7pkR48ebb0XMT4+XuXl5TYnunyvvfaavvWtb9kdo9NWrFih/v37Kz093dYcvvfpcRFtXdpuZ9tfrhMnTui2227Tb37zG0VGRtodp9NWrlyp2NhYjRs3zu4ol6WpqUlbt27V/PnztW3bNvXs2dPRh47+UU1NjZYvX67CwkKVlpbq5MmTeuONN+yOFbCeeuopuVwu3X333XZH6ZT6+no99dRTevLJJ+2O4j+l5A9LEjU2Nuq2227T3XffrZkzZ9od55Js3LhRK1asUGJiombNmqW//e1vuueee+yO1Wlut1tut7t17/T222/3qed3ffjhhxo8eLBiYmLUrVs3zZw5U5988ondsS5Zv379Wp8EUFZWptjYWJsTXbpFixZp5cqVevPNN33mD+P9+/ersLBQ6enpSkxMVElJicaOHasjR454PYvflNI3ly0yxrS5bJGTGWP0wAMPKCUlRT/60Y/sjnPJnn76aZWUlKioqEhLly7VxIkTfeov9bi4OA0YMEBfffWVJOmjjz7SyJEjbU7VeQMHDtTmzZtVX18vY4w++ugjn7pQ45ycnBwtWrRI0tkP91tuucXmRJdmzZo1euaZZ7RixYo2n2bgVKNHj1Z5ebmKiopUVFQkt9utrVu3Ki4uzvthOrNqqwe+POL99983w4YNMyEhIeZXv/qVp6bxiA0bNhhJZvTo0SY9Pd2kp6eb999/3+5Yl2XdunVm+vTpdse4ZNu2bTPjxo0zo0ePNrfccouprq62O9Ilefzxx01ycrJJTU0199xzj2loaLA7UrtmzZpl4uLijMvlMv379zevvvqqqaysNBMnTjRJSUlm4sSJpqqqyu6YF9VW/qFDhxq32936//BDDz1kd8w2tZX9mwYNGmQqKio8MXWH/cAyQwAAb2GZIQCA76CUAACOQSkBAByDUgIAOAalBABwDEoJAOAYlBIAwDEoJQCAx3z66adKS0tTQ0ODLMvqaVnWbsuyLnyQ09coJQCAx2RlZSknJ0c///nPJelZSW8YY3ZdbLzLa8kAAAHp8ccfV1ZWliRlSnq0vbHsKQEAPKq6ulonTpyQpAhJ3dsbSykBADxq7ty5+uUvfylJb0p6pr2xlBIAwGMWL14sl8ul7373u5K0QFKWZVkTLzaeVcIBAN7CKuEAAN9BKQEAHINSAgA4BqUEAHAMSgkA4BiUEgDAMSglAIBjUEoAAMeglAAAjmHLKuHV1dW68847VVRUpMTERL399tuKjo6+YFxiYqIiIiIUHBwsl8vFKg0A4Ods2VNasGCBJk2apIKCAk2aNEkLFiy46Nh169Zp+/btFBIABABbSmn58uWaPXu2JGn27NlatmyZHTEAAA5jSykdPXpU8fHxkqT4+HiVl5e3Oc6yLE2ZMkXjxo1Tbm5uu9vMzc1VZmamMjMzVVFR0eWZAQCe57FVwi3L+lBSXBtv/SwqKmpZbW1t6zeio6NVU1NzwcDS0lIlJCSovLxckydP1osvvqgJEyZ0ODerhAOAI3W4SrjHLnQwxtx8sfeSk5NVVlam+Ph4lZWVKTY2ts1xCQkJkqTY2FjdeuutysvL61QpAQB8ky2H73JycrRo0SJJ0qJFi3TLLbdcMObkyZM6fvx4679/8MEHGjVqlFdzAgC8y5ZSeuyxx7R27VoNGzZMa9eu1WOPPSbp7OG67OxsSWfPO40fP17p6em66qqrNH36dE2bNs2OuAAAL+HJswAAb+HJswAA30EpAQAcg1ICADgGpQQAcAxKCQDgGJQSAMAxKCUAgGNQSgAAx6CUAACOQSkBAByDUgIAOAalBABwDEoJAOAYlBIAwDEoJQCAY1BKAADHoJQAAI5BKQEAHINSAgA4BqUEAHAMSgkA4BiUEgDAMSglAIBjUEoAAMeglAAAjkEpAQAcg1ICADgGpQQAcAxKCQDgGJQSAMAxKCUAgGNQSgAAx6CUAACOQSkBABzDllJ65513lJqaqqCgIOXn51903Jo1a5ScnKykpCQtWLDAiwkBAHawpZRGjRql9957TxMmTLjomObmZj388MNavXq19uzZoyVLlmjPnj1eTAkA8DbLGOP1Se/8/SYjSdt3bNfQIUMVERFxwZi6ujoVFRUpLS1NknTo0CFJ0sCBAzvc/tr/mKPJ//ZaV0YG4AVvPXSt3RHgWVZHA1zeSHE5Tp8+rdDuoa2vQ0NDVXe87qLjy8rKVFpWevZFt+6ejgcA8ACP7SlZlvWhpLg23vqZMWaZJN1444167rnnlJmZecGgd955R3/961/16quvSpL++Mc/Ki8vTy+++KJH8gIAPM6+PSVjzM1X8vNut1vFxcWtr0tKSpSQkHDFuQAAzuXYS8KzsrJUUFCgwsJCnTlzRkuXLlVOTo7dsQAAHmRLKf35z3+W2+3Wpk2bNH36dE2dOlWSVFpaquzsbEmSy+XSSy+9pKlTpyolJUV33HGHUlNT7YgLAPASW66+k2TLpAAAW3V4Tsmxh+8AAIGHUgIAOAalBABwDEoJAOAYlBIAwDEoJQCAY1BKAADHoJQAAI5BKQEAHINSAgA4BqUEAHAMSgkA4BiUEgDAMSglAIBjUEoAAMeglAAAjkEpAQAcg1ICADgGpQQAcAxKCQDgGJQSAMAxKCUAgGNQSgAAx6CUAACOQSkBAByDUgIAOAalBABwDEoJAOAYlBIAwDEoJQCAY1BKAADHoJQAAI5BKQEAHINSAgA4BqUEAHAMSgkA4BiUEgDAMVw2zWvZNC8AwMHYUwIAOAalBABwDEoJAOAYlBIAwDEoJQCAY1BKAADH+H8TjMT7eP5fyAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "Conclusion: -1-1\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAACoAAAAOCAYAAABZ/o57AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAB0UlEQVRIDb2W7U0CQRCGD0IBqB1oBxAqEDuA2IF0oH/5Z6ADtAKFDqAEpAOxA0MH+Dzn7uVyhFNAb5JhZudj592vOWrb7TYZDodvSZI04Q38GWQLKWlbw/rPg1yTc4NeGTVCJUH1KT6LldHv0UfwAP0pZ++hP8dxVbIOCHdqlgdZVjzELcti/sNXZ1KPc37g5KuwwAPTjg8XqDvqHTyE3gl2gZVRg51ZHVHNhU3JvUQ+wo47cBfuwy7ee6x/QZy2BHmHGMBpHuOx9kiMtT/AbsQF3MRmfFL351Ai2eLtkHcV5ATpo+zi82rot4MIOiVsPsprOLN9e9JFmGv3GRE3hlPASG3HATUxkEAE5mN0V8+QWYdgvHOl8G9CblFMMbyGeVIfujveQvZieyom/XZs0QxMCYjS+chzNz32dPcKwV7NzqlAndNdPZUEKbUB7T3Ok1dq+RdA9x1lvthPejyVOUCzj04+qZ4fVKEDJO5eVg6bx+uCbzNjTsHfPRVok/nkfeSHpAjMFiYo20+e7AY+TP0ZMfYzvq75pyRSCLLn2cxNEISTLmDv4oQYW48++5vf/eh/wb5zbNi8Y5K90dhX+AOWlvizPzfoLsq2JMVeas31Fx2vogF49TMKAAAAAElFTkSuQmCC\n", "text/latex": [ "$\\displaystyle \\text{True}$" ], "text/plain": [ "True" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "Notons que q(x)>0 pour tout 0-1\")\n", "display(sym.solve(q>-1))\n", "print(\"Notons que q(x)>0 pour tout 00))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ">**Pour la question 6 on va fixer $(a,b,c)$. Pour cela, écrivez votre nom et prénom dans la cellule ci-dessous et vous obtiendrez un choix pour les paramètres.**" ] }, { "cell_type": "code", "execution_count": 61, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$\\displaystyle a=\\frac{1}{3},\\quad b=\\frac{1}{3},\\quad c=\\frac{1}{3}$" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "from IPython.display import display, Math\n", "import sympy \n", "sympy.init_printing()\n", "\n", "nom=\"Faccanoni\"\n", "prenom=\"Gloria\"\n", "\n", "DEN=list(range(3,11))\n", "idx=sum([ord(c) for c in nom+prenom])%len(DEN)\n", "#for idx in range(len(DEN)):\n", "my_c=sympy.S(1)/DEN[idx]\n", "my_b=1-2*my_c\n", "display(Math(r'a='+sympy.latex(my_c)+',\\\\quad b='+sympy.latex(my_c)+',\\\\quad c='+sympy.latex(my_b)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ">**Q6 [2 points]** \n", "Implémenter le schéma avec les valeurs $a,b,c$ trouvées avec vos nom et prénom et approcher la solution du problème de Cauchy suivant avec $N+1$ points (quelles valers de $N$ peut-on choisir? justifier la réponse)\n", "$$\n", " \\begin{cases}\n", " y'(t)=-50y(t), &t\\in[0;1]\\\\\n", " y(0)=1\n", " \\end{cases}\n", "$$ \n", " Vérifier ensuite empiriquement l'ordre de convergence." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Pour la A-stabilité il faut que $h$ soit $<\\frac{2}{50}$ donc $N>\\frac{1-0}{h}=\\frac{50}{2}$" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmEAAAFRCAYAAAA8Z3p7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3de5xddX3v/9cntyaBEHJjRhPCBAooUEggQEUtUY8/gUeV+mh/PdB56BFsU7B4jqd6KkqRWo2lHrwiCvFSqp1T4KhHEbnY2sdYK6DBGm6maA4wIQKGhNxDEpJ8zx9rT7Iz7JnsmVlrr71nv56Px37s2Wutvdcn+0uYd77ftb7fSCkhSZKkxhpXdgGSJEntyBAmSZJUAkOYJElSCQxhkiRJJTCESZIklcAQJkmSVAJDmCRJUgkMYZIkSSUwhEkaVEQ8FRGnl12HJI1FhjCpwSJiRkSkiLhvwPabIuJTw/ic34iIL0dEX0RsjYifRcT5NY67KCJWRcT2iPi/EfHaeusEXg78R701DfI5MyPi/1TO3xcRfzTSYyPiHyLimYjYEhG/iIg/HrC/KyLujIiNEfFsRHwuIiYMOOb4iNgZEf9Q4/yDflcRcUVEPBARuyLi5gHv66185rbK47F6666n5pHWPVTNzVx3Pe+VxgJDmNR4C4FngZMi4mUDtq8cxudMAJ4CzgWmA1cDt0VEV/8BEfFG4G+BS4BpwO8Aj9f5+b8FPJFS2jGMmmq5AdgNdADdwBci4uQRHvs3QFdK6QjgLcBHI+KMqv2fB9YBLyP7Ps8F3lXjHCsGnriO7+pp4KPAVwap/YqU0uGVx4kD9g1Vdz01j7TuQ9XcrHUP+V5prDCESY23EHgA+CeyX2xExHiy0POzej8kpbQ9pfRXKaUnU0r7Ukp3AE8A1aHkw8Bfp5Turxzzq5TSr+o8xanA/42Iz0TEcxHxdOUXZ90i4jDg94GrU0rbUkr/BtwOvG0kx6aUHk0p7ep/WXkcV/UxC4DbUko7U0rPAncD+0NcRFwEbAK+X6PcIb+rlNI3U0rfAjYM5zuoo+4hax5N3aOpucy663ivNCYYwqTGW0TW4/Ut4Pcq214BjAdWAUTEHRGxaZDHHbU+NCI6gBOARyuvxwOLgTkRsToi1laGjKbUWeeplfffSdYzdRPw/qrz1VPjCcDelNIvqj73QQb8sh7OsRHx+YjYQTZM+kylvn6fAS6KiKkRMRc4nywcEBFHAH8NvHfgiXP4rgD+JiLWR8SPImJJjXMMVvegNbdr3UO9VxpLDGFS4/UPO34XeG1ETKtseySl9CJASul3U0pHDvL43YEfGBETgR7g71NK/ddwdQATgT8AXls5xyLgL+us87eAZSmle1JK+4CfV++ss8bDgc0DPncz2fDTQHUdm1J6V2Xba4FvAruqdv+ALLRtAdaS9Th+q7LvI8CXU0pP1Tj3aL+r9wPHAnOB5cB3IqK6h26ouoequV3rHuq90phhCJMaKCJ+A3glsDKltBH4CVkPQn/v2Eg+cxzwNbJrqa6o2vVC5fn6lNIzKaX1wCeBC+r4zABOAb5TtfkUBgSxOmwDjhiw7Qhg62iOTSntrQxXzgMur9Q8DriHLCgcBswGZgB/GxELgf8EDHbjw4i/q0o9P04pbU0p7Uop/T3wo1rvHVj3UDVX/kxtV3cd75XGjJfcySKpUKeQ/QLqvwC5f0iyA/g//QdFxF1kPQS1/DCldH7luAC+XHn/Bf09aQAppY0RsZbsOp7hWlB5Xl21bRFVPR111vgLYEJEHJ9S+mVl32lUhkwHGM6x/SZw4BqlmcDRwOcq1zHtioi/I7sw/WmgC1iTfWUcDoyPiJNSSqeP8ruqJQFRR91D1fwXwJI2rHvI9472Dyg1lZSSDx8+GvQA/pgsoPS/Pobs4uPngdeM4PNuBO4HDh9k/1+T3V12FFlPxQ+Bj1Ttvxm4ucb7fg+4d8C2p4AzRlDjLcA/kvWYvJpsiPHk4R5b+TNcROWXMvAmYDtwYdX7HweuJAsLR5IF2x5gKtBZ9bgO+DowZxjf1QRgMtkdg1+r/Nx/njdVve6u1HViPXUPVnNl36jqHqzmyr6mrLue9/rwMVYepRfgw0c7PYDPkQ3BVG9bCewDpg3zs44h60nYSTaU1//orjpmItlUApvIpsX4LDC5av/3gT+p8dlXA1+oej0LeBH4jRH8mWeS9aBtB9YAfzRg/13ABw91LDCH7DqkTWTXIT08sHaya4t6gY3AeuB/A0fVqOmvgH8YsO1Q39VfceAOwf7HX1XqWkE2bLqJLBS/sd666615JHUPVnNVXU1Z96He68PHWHlESnn1YktqJRExiezuw1NT1TCmJKkxDGGSJEkl8O5ISZKkEhjCJEmSSmAIkyRJKoEhTJIkqQQtN1nr7NmzU1dXV+Hn2b59O4cddljh51H9bJPmY5s0J9ul+dgmzakR7fLTn/50fUppTq19LRfCurq6eOCBBwo/T29vL0uWLCn8PKqfbdJ8bJPmZLs0H9ukOTWiXSKib7B9DkdKkiSVwBAmSZJUAkOYJElSCVrumjBJklSuF198kbVr17Jz586ySxmV6dOns2rVqlw+a/LkycybN4+JEyfW/R5DmCRJGpa1a9cybdo0urq6iIiyyxmxrVu3Mm3atFF/TkqJDRs2sHbtWhYsWFD3+xyOlCRJw7Jz505mzZrV0gEsTxHBrFmzht0zaAiTJEnDZgA72Ei+j8JCWER8JSLWRcQjg+yPiPhsRKyOiIci4vSiahmWnh7o6uLc178eurqy15IkqamMHz+ehQsX7n9ce+21uX32ypUrufPOO3P7vMEUeU3YzcDngK8Osv984PjK42zgC5Xn8vT0wNKlsGMHAdDXl70G6O4uszJJklRlypQprFy5spDPXrlyJQ888AAXXHBBIZ/fr7CesJTSvwLPD3HIhcBXU+Z+4MiIeFlR9dTlqqtgx46Dt+3YkW2XJElNbfPmzZx44ok89thjAFx88cV88YtfBODyyy9n8eLFnHzyyVxzzTX737NixQrOOeccTjvtNM466yw2b97Mhz70IW699VYWLlzIrbfeyvbt27n00ks588wzWbRoEd/+9rdzqbfMuyPnAk9VvV5b2fbMwAMjYimwFKCjo4Pe3t5CCjp3zRpqjeimNWv4QUHnVP22bdtWWNtrZGyT5mS7NJ+x1ibTp09n69atpdbwwgsvcOqpp+5//ed//uf8/u//Ph//+Md529vexuWXX85zzz3HRRddxNatW7nyyiuZOXMme/fu5c1vfjPnnXcexx13HH/4h3/I3/3d33HGGWewZcsW9u3bxwc/+EH+/d//nU984hMAXHPNNbzqVa/iM5/5DJs2beJ1r3sdZ5999kvWndy5c+ew2rnMEFYz79Q6MKW0HFgOsHjx4lTYOk/z52dDkAPE/Pmu+dUEXHut+dgmzcl2aT5jrU1WrVq1f2qH97wH8h4VXLgQPv3poY+ZMmUKDz300Eu2X3jhhXz3u9/lfe97Hw8++OD+Ont6eli+fDl79uzhmWeeoa+vj4jg5S9/+f626T928uTJTJo0af/r3t5e7r77bm644QYAdu/ezcaNG+ns7Dzo3JMnT2bRokV1/znLDGFrgaOrXs8Dni6plsyyZaSlS4nqIcmpU2HZsvJqkiRJddu3bx+rVq1iypQpPP/888ybN48nnniC6667jhUrVjBjxgze8Y53sHPnTlJKdd3VmFLiG9/4BieeeGKutZYZwm4HroiIW8guyN+cUnrJUGRDdXcTwFNvv4q5+9Yw7pj5WQDzonxJkmo6VI9Vo33qU5/ila98JR/72Me49NJLue+++9iyZQuHHXYY06dP59e//jV33XUXS5Ys4YQTTuDpp59mxYoVnHnmmWzdupUpU6Ywbdq0g4Zb3/SmN3H99ddz/fXXExH87Gc/G1aP12AKC2ER8Y/AEmB2RKwFrgEmAqSUbgTuBC4AVgM7gEuKqmVYurt5y3XdTJmynnvvnV12NZIkqYYXXniBhQsX7n993nnncemll/KlL32Jn/zkJ0ybNo3f+Z3f4aMf/Sgf/vCHWbRoESeffDLHHnssr371qwGYNGkSt956K+9+97t54YUXmDJlCv/8z//M6173Oq699loWLlzIBz7wAa6++mre8573cOqpp5JSoqurizvuuGPUf4bCQlhK6eJD7E/AnxV1/tHo6IAnn5xUdhmSJGkQe/furbm9ei3IT37yk/t/vvnmm19y7NatWznzzDO5//77X7JvxYoVB72+6aabRljp4Jwxv4bOTti40RAmSZKKYwirobMTnn9+EqnmvZqSJEmjZwiroaMD9uwZx8aNZVciSZLGKkNYDf3Tfvz61+XWIUmSxi5DWA39IezZZ8utQ5IkjV2GsBo6OrJnQ5gkSSqKIawGhyMlSWpu48ePZ+HChZxyyim8+c1vZtOmTQA8+eSTnHLKKfuP++IXv8jpp5/Oxia80NsQVsOMGTBhwj57wiRJalJTpkxh5cqVPPLII8ycOXP/uo7Vvva1r3H99dfzve99jxkzZpRQ5dAMYTVEwIwZuw1hkiTloacHurpg3Ljsuacn149/1atexa9+9auDtt12221ce+21fO9732P27OZcAafMtSOb2syZu3n22clllyFJUmvr6YGlS2HHjux1X1/2GnJZm3nv3r18//vf553vfOf+bX19fVxxxRX87Gc/o7P/GqMmZE/YIGbO3D36a8IKTv6SJDW9q646EMD67diRbR+F/rUjZ82axfPPP88b3/jG/fvmzJnD/Pnzue2220Z1jqIZwgYxY8aLoxuO7E/+fX2Q0oHkbxCTJLWTNWuGt71O/deE9fX1sXv37oOuCZs6dSp33XUXN954Iz1N/HvXEDaImTN3s24dDLI+6KEVlPwlSWop8+cPb/swTZ8+nc9+9rNcd911vPjii/u3z5kzh7vvvpsPfvCD3HPPPbmcK2+GsEHMnLmbfftgw4YRfkBByV+SpJaybBlMnXrwtqlTs+05WbRoEaeddhq33HLLQdsXLFjA7bffzqWXXsqPf/zj3M6XFy/MH8SMGbuBbMLWo44awQfMn58NQdbaLklSu+i/+P6qq7KOiPnzswA2yovyt23bdtDr73znO/t/fuSRR/b/fNppp73kzslmYU/YIGbOPBDCRqQByV+SpJbQ3Q1PPgn79mXPOdwVORYYwgbRH8JGfIdkdzcsXw7HHJNNPHbMMdlr/8OTJEk4HDmo6uHIEevuNnRJkqSa7AkbxNSpe5kyxUW8JUmqJaVUdglNZSTfhyFsEBHZQt4u4i1J0sEmT57Mhg0bDGIVKSU2bNjA5MnDW2nH4cghdHTYEyZJ0kDz5s1j7dq1PPfcc2WXMio7d+4cdnAazOTJk5k3b96w3mMIG0JnJ/zyl2VXIUlSc5k4cSILFiwou4xR6+3tZdGiRaWd3+HIITgcKUmSimIIG0JHB6xfD1WrIEiSJOXCEDaEzs7sed26cuuQJEljjyFsCP0hzCFJSZKUN0PYEPpDmHdISpKkvBnChtDRkT03TQjr6YGuLhg3Lnvu6Sm7IkmSNEJOUTGE/hDWFMORPT2wdCns2JG97uvLXoNLI0mS1ILsCRvC1KlwxBFN0hN21VUHAli/HTuy7ZIkqeUYwg6haWbNX7NmeNslSVJTM4QdQtNM2Dp//vC2S5KkpmYIO4TOzibpCVu2LBsfrTZ1arZdkiS1HEPYITTNcGR3NyxfDsccAxHZ8/LlXpQvSVKL8u7IQ+jshM2bYedOyGmh9ZHr7jZ0SZI0RtgTdgjOmi9JkopgCDuEppuwVZIkjQmGsEOwJ0ySJBXBEHYIrh8pSZKKYAg7hKOOyp4NYZIkKU+FhrCIOC8iHouI1RFxZY390yPiOxHxYEQ8GhGXFFnPSEyaBDNnOhwpSZLyVVgIi4jxwA3A+cBJwMURcdKAw/4M+HlK6TRgCfCJiJhUVE0j1TQTtkqSpDGjyJ6ws4DVKaXHU0q7gVuACwcck4BpERHA4cDzwJ4CaxqRppmwNU89PdDVBePGZc89PWVXJElSWykyhM0Fnqp6vbayrdrngFcCTwMPA/8tpbSvwJpGpGnWj8xLTw8sXQp9fZBS9rx0qUFMkqQGKnLG/KixLQ14/SZgJfB64DjgnyLihymlLQd9UMRSYClAR0cHvb29+Vc7wLZt2/af58UXj+NXv3o5vb0/LPy8jfDb730vk3fsOHjjjh3sfO97uX/uwJzcPKrbRM3BNmlOtkvzsU2aU9ntUmQIWwscXfV6HlmPV7VLgGtTSglYHRFPAK8AflJ9UEppObAcYPHixWnJkiVF1bxfb28v/ef58Y/h61+HxYuXcPjhhZ+6eOvW1dw8ed06GvHdjlR1m6g52CbNyXZpPrZJcyq7XYocjlwBHB8RCyoX218E3D7gmDXAGwAiogM4EXi8wJpGZMxN2Dp//vC2S5Kk3BUWwlJKe4ArgHuAVcBtKaVHI+KyiLiscthHgHMi4mHg+8D7U0rri6pppMbchK3LlsHUqQdvmzo12y5JkhqiyOFIUkp3AncO2HZj1c9PA/9fkTXkYcytH9ndnT1fdRWsWZP1gC1bdmC7JEkqXKEhbKwYc8ORkAUuQ5ckSaVx2aI6zJmTTac1ZnrCJElS6QxhdRg/HmbPNoRJkqT8GMLqNOYmbJUkSaUyhNXJ9SMlSVKeDGF1GpPrR0qSpNIYwurU3xOWBi68JEmSNAKGsDp1dsLu3bB5c9mVSJKkscAQVqcxN2Fr3np6oKsrm8ujqyt7LUmSBmUIq9OYW7ooTz09sHQp9PVl47V9fdlrg5gkSYMyhNVpTM6an5erroIdOw7etmNHtl2SJNVkCKuTPWFDWLNmeNslSZIhrF4zZsCECYawmubPH952SZJkCKvXuHHZxfkOR9awbBlMnXrwtqlTs+2SJKkmQ9gwOGv+ILq7YflyOOYYiMiely/PtkuSpJomlF1AK+nogGeeKbuKJtXdbeiSJGkY7AkbBhfxliRJeTGEDUN/CNu3r+xKJElSqzOEDUNHB+zdCxs2lF2JJElqdYawYXDCVkmSlBdD2DA4YaskScqLIWwYXMS7QVwMXJLUBpyiYhgcjmyA/sXA+9ei7F8MHGDu3PLqkiQpZ/aEDcMRR8DkyfaEFcrFwCVJbcIQNgwR2ZCkIaxALgYuSWoThrBhcsLWgrkYuCSpTRjChsn1IwvmYuCSpDZhCBsmhyML5mLgkqQ24d2Rw9TZCevXw549MMFvrxguBi5JagP2hA1TZyekBM89V3YlkiSplRnChskJWyVJUh4MYcPkhK2SJCkPhrBhcv1ISZKUB0PYMDkcKUmS8mAIG6bDDoPDD3c4sqW4ILgkqQk5ycIIOGFrCxlqQXCnwZAklciesBFwwtYW4oLgkqQmZQgbAdePbCEuCC5JalKGsBFwOLKFuCC4JKlJGcJGoKMDNm6EXbvKrkSH5ILgkqQmVWgIi4jzIuKxiFgdEVcOcsySiFgZEY9GxA+KrCcv/XOFrVtXbh2qgwuCS5KaVGF3R0bEeOAG4I3AWmBFRNyeUvp51TFHAp8HzksprYmIo4qqJ0/VE7YefXS5tagOLgguSWpCRfaEnQWsTik9nlLaDdwCXDjgmD8CvplSWgOQUmqJviUnbJUkSaNVZAibCzxV9XptZVu1E4AZEdEbET+NiLcXWE9uXD9SkiSNVpGTtUaNbanG+c8A3gBMAe6LiPtTSr846IMilgJLATo6Oujt7c2/2gG2bds26Hl27w7gXO677wl+8zf7Cq9FmaHaROWwTZqT7dJ8bJPmVHa7FBnC1gLVV0zNA56uccz6lNJ2YHtE/CtwGnBQCEspLQeWAyxevDgtWbKkqJr36+3tZajzzJgBU6YsYMmSBYXXosyh2qQhenqyiV7XrMmmuVi2rK2vN2uKNtFL2C7NxzZpTmW3S5HDkSuA4yNiQURMAi4Cbh9wzLeB10bEhIiYCpwNrCqwptx0dDgc2Xb6l0Dq64OUDiyB5FqUkqQRKCyEpZT2AFcA95AFq9tSSo9GxGURcVnlmFXA3cBDwE+AL6WUHimqpjw5YWsbcgkkSVKOCl3AO6V0J3DngG03Dnj9P4H/WWQdRejshAceKLsKNZRLIEmScuSM+SPkcGQbcgkkSVKODGEj1NkJW7fC9u1lV6KGcQkkSVKODGEj5FxhbcglkCRJOSr0mrCxrH/W/F//Go49ttxa1EAugSRJyok9YSNUvX6kJEnScBnCRsgQJkmSRsMQNkJz5mSXBRnCJEnSSBjCRmjCBJg92wvzNUo9PdDVBePGZc/Ovi9JbcML80fBWfM1Kv3LIPXPwt+/DBJ48b8ktQF7wkaho8MQplFwGSRJamuGsFHo7HQ4UqPgMkiS1NYMYaPQPxyZUtmVqCW5DJIktTVD2Ch0dMDOnbBlS9mVqCW5DJIktTVD2Ci4dJFGxWWQJKmteXfkKFRP2HrCCeXWohblMkiS1LbsCRuF/vUjvUNSkiQNlyFsFByOlCRJI2UIG4VZs2D8eHvC1CScfV+SWsohrwmLiBXAQ8DD/c8ppeeKLqwVjBsHRx1lCFMTcPZ9SWo59fSEXQj8b2AScBnwZET0FVpVC3HCVjUFZ9+XpJZzyJ6wlNLTwNPA3QAR8UrgDwquq2W4fqSagrPvS1LLOWRPWEQcNH13SmkVcHJhFbUY149UU3D2fUlqOfUMR94aEWsj4ocR8fmI+CTwiqILaxWdnbBuHezbV3YlamvOvi9JLeeQISyl9KqU0jzgEuCfgEeB3y26sFbR2QkvvggbN5Zdidqas+9LUsupe8b8lNJqYHWBtbSk6glbZ80qtxa1OWffl6SW4jxho+SErZIkaSQMYaNUvX6kNKY4+askFcoFvEfJEKYxyclfJalw9oSN0vTpMGmSw5EaY5z8VZIKZwgbpQgnbNUY5OSvklQ4Q1gODGEac5z8VZIKZwjLQUeHw5EaY5z8VZIKZwjLgT1hGnOc/FWSCufdkTno7ITnnoO9e2H8+LKrkXLi5K+SVCh7wnLQ0ZGtHbl+fdmVSE3KOcck6SUMYTlwrjBpCP1zjvX1QUoH5hwziElqc4awHBjCpCE455gk1WQIy0H/It7eISnV4JxjklSTISwH9oRJQ3DOMUmqyRCWg8MPh8MOM4RJNTnnmCTVZAjLiRO2SoNwzjFJqqnQEBYR50XEYxGxOiKuHOK4MyNib0T8QZH1FMkJW6UhdHfDk09mc7k8+aQBTJIoMIRFxHjgBuB84CTg4og4aZDj/ha4p6haGsEQJjWQ845JGgOK7Ak7C1idUno8pbQbuAW4sMZx7wa+AawrsJbCORwpNYjzjkkaI4oMYXOBp6per61s2y8i5gJvBW4ssI6G6OyEDRtg9+6yK5HGOOcdkzRGFLl2ZNTYlga8/jTw/pTS3ohah1c+KGIpsBSgo6OD3t7evGoc1LZt24Z1ns2bXwacyLe/fR9z5uwqrK52Ntw2UfHKaJNz16yp/T+XNWv4gf99AP5daUa2SXMqu12KDGFrgaOrXs8Dnh5wzGLglkoAmw1cEBF7Ukrfqj4opbQcWA6wePHitGTJkqJq3q+3t5fhnGfzZvjkJ+HYY1/FGWcUV1c7G26bqHiltMn8+dkQ5AAxf77/fVT4d6X52CbNqex2KXI4cgVwfEQsiIhJwEXA7dUHpJQWpJS6UkpdwNeBdw0MYK3CCVulBnHeMUljRGEhLKW0B7iC7K7HVcBtKaVHI+KyiLisqPOWxRAmNYjzjkkaIwqdJyyldGdK6YSU0nEppWWVbTemlF5yIX5K6R0ppa8XWU+RXD9SaqA85x1zugtJJSnymrC2MnkyTJ9uT5jUUvqnu+i/27J/uguwZ01S4Vy2KEdO2Cq1GKe7kFQiQ1iOOjoMYVJLWbNmeNslKUeGsBx1dnpNmNRS5s8f3nZJypEhLEcOR0otxukuJJXIEJajjg7YsgVeeKHsSiTVxekuJJXIEJaj/rnCHJKUWkie012AU15IqpshLEdO2Cq1uf4pL/r6IKUDU14YxCTVYAjLUf+ErYYwqU055YWkYTCE5cjhSKnNOeWFpGEwhOXoqKOyZ3vCpDbllBeShsEQlqOJE2H2bEOY1Lac8kLSMBjCctbR4XCk1LbynvLCOy2lMc0FvHPmhK1Sm+vuzmeeMRcXl8Y8e8JyZgiTlAvvtJTGPENYzvqHI1MquxJJLc07LaUxzxCWs87O7B+r27aVXYmkluadltKYZwjLmbPmS8qFd1pKY54hLGf9s+Z7h6SkUSlicXHvtpSaindH5syeMEm5yetOS/BuS6kJ2ROWM0OYpKbk3ZZS0zGE5WzWrKyn3+FISU3Fuy2lpmMIy9n48dkakvaESWoq3m0pNR1DWAGcsFVS0/FuS6npGMIK4PqRkpqO61pKTce7IwvQ2Qk//3nZVUjSAK5rKTUVe8IK0Nnp0kWSxjDvtJRyYQgrQEcH7N4NmzaVXYkkFcA7LaVcGMIK4Fxhksa0Iu609BoztSFDWAEMYZLGtLzvtOy/xqyvL7uOo/8aM4OYxjhDWAFcP1LSmJb3nZZeY6Y25d2RBbAnTNKYl+e6ll5jpjZlT1gBZsyAiRMNYZJUF2fzV5syhBUgwglbJaluRVxj5kX+agGGsIK4dJEk1SnPa8y8yF8txBBWEEOYJA1Ddzc8+STs25c9e5G/2oAhrCAOR0pSCbzIXy3EEFaQzk5Ytw727i27EklqI04kqxZiCCtIZ2cWwDZsKLsSSWojTiSrFmIIK4gTtkpSCZxIVi2k0BAWEedFxGMRsToirqyxvzsiHqo87o2I04qsp5GcsFWSSpLXRf7gNWYqVGEhLCLGAzcA5wMnARdHxEkDDnsCODeldCrwEWB5UfU0miFMksaAvK8x8/oyVSmyJ+wsYHVK6fGU0m7gFuDC6gNSSvemlDZWXt4PzCuwnoZyOFKSxoA8rzHz+jINUGQImws8VfV6bWXbYN4J3FVgPQ01bRpMmWJPmCS1tDyvMfP6Mg1Q5ALeUWNbqnlgxOvIQthrBtm/FFgK0NHRQW9vb04lDm7btm2jPs+RR57Ngw9uobd3VT5Ftbk82kT5sk2ak+2Ss7lz4eabD942zO932/e6UAsAAAxvSURBVLZtpDVrav9iXLOGH9hepSj770qRIWwtcHTV63nA0wMPiohTgS8B56eUak7okFJaTuV6scWLF6clS5bkXuxAvb29jPY8XV0AU1iypCOHipRHmyhftklzsl2aT29vLzF/fjYEOUDMnz/y9urpyXrS1qzJrlNbtmx0NyK0mbL/rhQ5HLkCOD4iFkTEJOAi4PbqAyJiPvBN4G0ppV8UWEspOjocjpQkVTiHmQYoLISllPYAVwD3AKuA21JKj0bEZRFxWeWwDwGzgM9HxMqIeKCoesrg+pGSpP2cw0wDFDpPWErpzpTSCSml41JKyyrbbkwp3Vj5+Y9TSjNSSgsrj8VF1tNonZ3ZjPkvvlh2JZKkptDMc5g5fUbDOWN+gTo6sh7i554ruxJJ0piT5xxmDm2WwhBWICdslSQVJs9rzBzaLIUhrECGMElSYfK8xqyI5Zkc3jykIqeoaHvOmi9JKlR3dz5TUgwyfcaolmdauvRA71r/8CY4hUYVe8IK1B/C7AmTJDW1vKfPcHizLoawAk2dCkccYQiTJDW5vKfP8M7NujgcWbCODocjJUktIK+hTch3eHMMD23aE1YwJ2yVJLUd79ysiyGsYCecACtWjO4GE0mSWkqL3Ll57utfX+rwpiGsYH/5l9nzu96VzX8nSVJbyGt1gDwnpYWDJqaNkiemNYQVrKsLPvIR+O534etfL7saSZJazBi+c9MQ1gD/9b/C6afDu98NGzeWXY0kSS2k2e/cHAVDWANMmABf/GK2huSVV5ZdjSRJLSbPhc/zHt4cBUNYg5x+Ovz3/56F9x/+sOxqJElqU3kPb46CIayBPvzhrBd16VLYtavsaiRJakNVw5spj+HNUTCENdBhh8EXvgD/8R9w7bVlVyNJUpuqDG/+4F/+ZfTDm6NgCGuw88+Hiy+Gj30MVq0quxpJklQWQ1gJPvWprFfsT/80u8ZQkiS1H0NYCTo64Lrrsgv0v/zlsquRJEllMISV5JJLYMkS+B//A555puxqJElSoxnCShIBN90EO3fCe95TdjWSJKnRDGElOuGEbG3J226DO+4ouxpJktRIhrCS/cVfwEknZQt8b9tWdjWSJKlRDGElmzQpW9Loqafg6qvLrkaSJDWKIawJnHMOXH45fPazsGJF2dVIkqRGMIQ1ib/5m2zqij/5E3jxxbKrkSRJRTOENYnp0+Fzn4MHH4RPf7rsaiRJUtEMYU3krW+FCy+Ea66Bxx8vuxpJklQkQ1gTich6w8aPz64RS6nsiiRJUlEMYU1m3rzs+rDvfQ/+1/8quxpJklQUQ1gTuvxyOPvsbCb9DRvKrkaSJBXBENaExo+H5cth0yZ43/vKrkaSJBXBENakTj01W9z75pvhX/6l7GokSVLeDGFN7Oqr4bjj4E//FF54oexqJElSngxhTWzKFLjpJli9Gj760bKrkSRJeTKENbk3vAH+y3+Bj38cHn647GokSVJeDGEt4Lrr4MgjsyWN9u4tuxpJkpQHQ1gLmD0bPvUp+PGP4cYby65GkiTlYULZBag+3d3w1a/CBz4AEyfCa14Dr3gFjDNGS5LUkgxhLSIi6wV7/euzuyUBZsyAc86BV786e5x5ZnYxvyRJan6GsBZy7LHwxBPwy1/Cj34E//Zv2fN3v5vtnzgRzjjjQCh79avhqKPKrVmSJNVWaAiLiPOAzwDjgS+llK4dsD8q+y8AdgDvSCn9e5E1tboIOOGE7HHJJdm29evh3nuzQPajH8H118MnPpHt+83fPBDIXvMaOPFEhzAlSWoGhYWwiBgP3AC8EVgLrIiI21NKP6867Hzg+MrjbOALlWcNw+zZ8Ja3ZA+AnTvhpz89EMruuAP+/u+zfTNnHjyEefLJ2RDm5MlZwJMkSY1RZE/YWcDqlNLjABFxC3AhUB3CLgS+mlJKwP0RcWREvCyl9EyBdY15kycfCFkAKcEvfnHwEOYdd9R+3+TJWSjrD2b1/tz/PG5c9og48Kh+PZx9/Y9+jzwym+efz34eKjAOti/vkGlohYcfnsWWLWVXoYFsl+ZjmzSn5547jCVLyjt/kSFsLvBU1eu1vLSXq9Yxc4GDQlhELAWWAnR0dNDb25t3rS+xbdu2hpynkY49Nnu8/e2wadNEHnnkCJ59dgq7d49j167ssXv3gceuXePZtWsc27ePY+PGA/t37Rpfdcw49uxp1PjmKQ06j+r3W2UXoJpsl+ZjmzSj88+fw3HH9ZZ2/iJDWK1+gjSCY0gpLQeWAyxevDgtaUBs7e3tpRHnKdPv/V4+n7N3bzYE+sILsG9f1vPW/9z/qH49nH3VVqxYwZlnnvmS7dUG2zfUe0Yi789rVQ888ACLFy8uuwwNYLs0H9ukOf3yl8+U+ru+yBC2Fji66vU84OkRHKMmN348HHZY9ijS889v59RTiz2HhmfLlm2cfnrZVWgg26X52CbNacuWXaWev8hxpBXA8RGxICImARcBtw845nbg7ZH5bWCz14NJkqR2UFhPWEppT0RcAdxDNkXFV1JKj0bEZZX9NwJ3kk1PsZpsiopLiqpHkiSpmRQ6T1hK6U6yoFW97caqnxPwZ0XWIEmS1IyctlOSJKkEhjBJkqQSGMIkSZJKYAiTJEkqgSFMkiSpBIYwSZKkEhjCJEmSShCpxRbBi4jngL4GnGo2sL4B51H9bJPmY5s0J9ul+dgmzakR7XJMSmlOrR0tF8IaJSIeSCm52moTsU2aj23SnGyX5mObNKey28XhSEmSpBIYwiRJkkpgCBvc8rIL0EvYJs3HNmlOtkvzsU2aU6nt4jVhkiRJJbAnTJIkqQRtHcIi4ryIeCwiVkfElTX2R0R8trL/oYg4vYw6200d7dJdaY+HIuLeiDitjDrbyaHapOq4MyNib0T8QSPra1f1tEtELImIlRHxaET8oNE1tps6/v81PSK+ExEPVtrkkjLqbCcR8ZWIWBcRjwyyv7Tf9W0bwiJiPHADcD5wEnBxRJw04LDzgeMrj6XAFxpaZBuqs12eAM5NKZ0KfASvtShUnW3Sf9zfAvc0tsL2VE+7RMSRwOeBt6SUTgb+/4YX2kbq/LvyZ8DPU0qnAUuAT0TEpIYW2n5uBs4bYn9pv+vbNoQBZwGrU0qPp5R2A7cAFw445kLgqylzP3BkRLys0YW2mUO2S0rp3pTSxsrL+4F5Da6x3dTzdwXg3cA3gHWNLK6N1dMufwR8M6W0BiClZNsUq542ScC0iAjgcOB5YE9jy2wvKaV/JfueB1Pa7/p2DmFzgaeqXq+tbBvuMcrXcL/zdwJ3FVqRDtkmETEXeCtwYwPranf1/F05AZgREb0R8dOIeHvDqmtP9bTJ54BXAk8DDwP/LaW0rzHlaRCl/a6f0IiTNKmosW3graL1HKN81f2dR8TryELYawqtSPW0yaeB96eU9mb/wFcD1NMuE4AzgDcAU4D7IuL+lNIvii6uTdXTJm8CVgKvB44D/ikifphS2lJ0cRpUab/r2zmErQWOrno9j+xfJsM9Rvmq6zuPiFOBLwHnp5Q2NKi2dlVPmywGbqkEsNnABRGxJ6X0rcaU2Jbq/X/Y+pTSdmB7RPwrcBpgCCtGPW1yCXBtyuaHWh0RTwCvAH7SmBJVQ2m/69t5OHIFcHxELKhcFHkRcPuAY24H3l65c+K3gc0ppWcaXWibOWS7RMR84JvA2/wXfUMcsk1SSgtSSl0ppS7g68C7DGCFq+f/Yd8GXhsREyJiKnA2sKrBdbaTetpkDVnPJBHRAZwIPN7QKjVQab/r27YnLKW0JyKuILuTazzwlZTSoxFxWWX/jcCdwAXAamAH2b9gVKA62+VDwCzg85Welz0ujFucOttEDVZPu6SUVkXE3cBDwD7gSymlmrfpa/Tq/LvyEeDmiHiYbBjs/Sml9aUV3QYi4h/J7kSdHRFrgWuAiVD+73pnzJckSSpBOw9HSpIklcYQJkmSVAJDmCRJUgkMYZIkSSUwhEmSJJXAECZJklQCQ5gkSVIJDGGS2l5EzIuI/1x2HZLaiyFMkrJlZE4vuwhJ7cUZ8yW1tYh4Ddkai5uArcBbU0pPlFuVpHZgCJPU9irrK77PdRUlNZLDkZIEJwKPlV2EpPZiCJPU1iJiFrA5pfRi2bVIai+GMEntbgHwdNlFSGo/hjBJ7e4/gNkR8UhEnFN2MZLahxfmS5IklcCeMEmSpBIYwiRJkkpgCJMkSSqBIUySJKkEhjBJkqQSGMIkSZJKYAiTJEkqgSFMkiSpBP8PovyDur5u7TQAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAf0AAAFRCAYAAAB+EnQdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deXxcdb3/8dcnaZO0TVfapnvSNGloKS3QUtQrPyjqlcvVi1xF0eJyVYosZRFQoagsooigoiyKgF4RLShcQWUTbfVyVWjLZgtkaZomaUK6L2mayTKf3x+ZYhqTdtJm5syceT8fjzweOUvmfD6dnrwz3znzPebuiIiISPhlBV2AiIiIJIdCX0REJEMo9EVERDKEQl9ERCRDKPRFREQyhEJfREQkQyj0RUREMoRCXySDmdlPzOxrQdchIsmh0BdJU2b2KTP7u5m1mNmbZna3mY0Kuq6BYma5ZnafmW00sz1m9pKZ/dtB9r/YzFabWcTMftLL9p+ZWaOZ7TazCjP7bEIbEElBCn2RNGRmVwDfBK4CRgJvAwqB35tZTi/7D+rn4/dr/wQZBNQBp9DV45eBh82sqI/9G4CvAff3sf0bQJG7jwD+A/iamc0fyIJFUp1CXyTNmNkI4Hpgqbs/5e7t7l4DfJiu4D83tl+NmX3RzF4F9prZIDM73sxejL1yfgjI6/a4ve0/ycweMbMtZrbBzC45RG3DzWyTmZ3SY/1UM3MzOyrePt19r7tf5+417h51998CG4Beg9rdH3X3XwPb+ti+zt0j+xdjXzPirUckDBT6IunnHXSF9aPdV7p7M/Ak8J5uqz8K/Dswiq7z/dfAA8AY4JfAB3s8dvf9o8BvgFeAycC7gMvM7L0Hqe0KYK27/6lHbXXAXuBYADP7rZnt7OPrt709sJkVADOBdQc5/kGZ2V1m1gK8ATQCTxzuY4mkI4W+SPoZC2x1945etjXGtu/3PXevc/d9dL0FMBj4bmx04FfAqh4/333/E4Fx7n6Du7e5ezXwI+Cc3ooys2zgAuDe2PI4MyvutksHMATA3d/n7qP6+HpfL489GHgQ+G93f+MQ/z59cvcLgeHAyXT90RQ5+E+IhItCXyT9bAXG9vG++8TY9v3qun0/CdjkB95ac2OPn+++fyEwqfurcOAaoKCPuuYA44GnY8ufB5YAmNkQusJ2c59d9cHMsuganWgDLu7vz/fk7p3u/hwwha4/UkQyhkJfJP38la5XqP/ZfaWZDQP+DfhDt9XdA74RmGxm1m3dtB6P3X3/OmBDj1fhw939jD7qmgzscPfdseXT+UfInwLsAF6K1fqkmTX38fVkt54MuI+uPzQ+6O7tfRz7cAxC7+lLhlHoi6QZd99F14V83zez081scOyK9l8C9XS9Ku7NX+kaYr8kdpHefwILD3KoF4DdsYv7hphZtpnNMbMT+9h/OzDCzKab2UeBHGB27GOE19H1tkI01sO/uXt+H1/dP5Z3NzALeH/sLYc+xXrKA7KBbDPL2z8aYmbjzewcM8uP9fFeuq5f+OPBHlMkbBT6ImnI3W+ha6j9VmA38Dxdr8zf1e0K9Z4/00bX6MCn6HrV/RF6XAzYY/9O4P3AcXRdNb+VrvfrR/bxI6uA5cDLwGfo+ljcO4DKWH3f7EeLmFkhcH7s+G92GwlYHNv+pJld0+1HrgX2AV+i6xMM+2LroGsE4wK6/ijaQde/22Xu/lh/ahJJd3bg23siIiISVnqlLyIikiEU+iIiIhlCoS8iIpIhFPoiIiIZQqEvIiKSIVLhTloJNXbsWC8qKkroMfbu3cuwYcMSeoxkCEsfEJ5ewtIHhKeXsPQB4elFfRxozZo1W919XG/bQh/6RUVFrF69OqHHWLlyJaeeempCj5EMYekDwtNLWPqA8PQSlj4gPL2ojwOZWc/ptd+i4X0REZEModAXERHJEAp9ERGRDKHQFxERyRAKfRERkQyh0BcREckQCn0REZEMEfrP6YuIiKSiB5uaWFZdTW0kwrTcXM4FTk3wMfVKX0REJMkebGpiSXk5GyMRHNgYiXBrbH0iKfRFRESSbFl1NS3R6AHrIrH1iaTQFxERSbLaSKRf6weK3tMXERFJkn2dndywcSPex/ZpubkJPb5CX0REJAme3b6dz1VUsL61lf83YgSrmpvZ122IPxe4qbg4oTVoeF9ERCSBtrS18YnXX+c9r75Klhl/mDePP51wAj8qK6MwNxcDCnNzuRJYXFCQ0Fr0Sl9ERCQB3J2fNjVxRVUVuzs7ubawkGXTppGXnQ10BXz3kF+5cmXCa1Loi4iIDLDKlhY+V1HBH3fu5B0jRnBPWRnHDBsWdFkKfRERkYHSFo1ya10dN9TUkJuVxd2lpSyZNIkss6BLAxT6IiIiA+Ivu3axpLycdS0tfGjcOG4vKWFSgq/G7y+FvoiIyBHY1dHB1dXV/KChgSm5uTw+Zw7vHzs26LJ6pdAXERE5DO7Oo1u3srSykqa2Ni6ZPJkbp09n+KDUjdbUrUxERCRF1bW2cnFlJY9v28Zx+fk8PmcOC0aMCLqsQ1Loi4iIxKnTnTs3bWLZhg10uvOt4mIumzKFQVnpMe2NQl9ERCQOrzQ3c155Oav27OH0MWO4q7SU6UOGBF1Wvyj0RUREDqKls5Pramr4dl0dRw0ezM9nzeKc8eOxFPkYXn8o9EVERPrw9PbtXFBRwYbWVj4zYQK3zJjBmMGDgy7rsCn0RUREetjc1sblVVX8fPNmyoYMYeVxx3HKqFFBl3XEAr/ywMyuM7NNZvZy7OuMPvY73czKzazKzL6U7DpFRCT83J37Gxs5+oUX+OWWLXy1sJBXTjwxFIEPqfNK/zvufmtfG80sG7gTeA9QD6wys8fd/bVkFSgiIuFW0dLC+RUVrNy5k5NHjuSHM2cyKwXmyx9IqRL6h7IQqHL3agAzWw6cCSj0RUTkiLRFo3yztpabNm4kLyuLe2bO5DMTJ6bMfPkDKVVC/2Iz+wSwGrjC3Xf02D4ZqOu2XA+clKziREQknJ7buZMlFRW83tLCR8aN47slJUxIsfnyB5K5e+IPYvYsMKGXTcuAvwFbAQduBCa6+6d7/PzZwHvd/bOx5Y8DC919aR/HWwIsASgoKJi/fPnygWqlV83NzeTn5yf0GMkQlj4gPL2EpQ8ITy9h6QPC08vh9NEM3AP8BigALgPeNvCl9ctAPR+LFi1a4+4Let3o7inzBRQBa3tZ/3bg6W7LVwNXx/OY8+fP90RbsWJFwo+RDGHpwz08vYSlD/fw9BKWPtzD00t/+ohGo/5wU5NP+L//86wVK/zzlZW+p709ccX1w0A9H8Bq7yMTAx/eN7OJ7t4YWzwLWNvLbquAUjObDmwCzgE+lqQSRUQkBDa2tnJRRQW/276dE/Lz+e2xxzJ/+PCgy0qqwEMfuMXMjqNreL8GOB/AzCYB97r7Ge7eYWYXA08D2cD97r4uqIJFRCR9dLrz/fp6rt2wAQe+PWMGSydPTpv58gdS4KHv7h/vY30DcEa35SeAJ5JVl4iIpL+X9uzhvPJy1jQ3c8aYMdw1cyaFeXlBlxWYwENfRERkoO3t7OSrGzbwnfp6xg0ezEOzZ3P2uHFpOV/+QFLoi4hIqDy5bRsXVFSwMRJhycSJ3FxczOg0ni9/ICn0RUQkFJra2risqorlmzcza+hQ/ve443hnSKbPHSgKfRERSWtR4N6GBq6qrqals5Pri4r44rRp5GbghXqHotAXEZG09cbevVwOvFpRwSkjR/LDsjLKhg4NuqyUpdAXEZG0E4lG+cbGjXyjtpZc4L6yMv5rwoSMv1DvUBT6IiKSVv68cyfnV1TwRksLHxs/ng9t3sxZEycGXVZa0BseIiKSFna0t3NeeTmnvPwyrdEoTx57LA/Ons3ooAtLI3qlLyIiKc3deWjzZi6rqmJreztXTZ3KV4uKGJadHXRpaUehLyIiKatm3z4urKzkye3bWTB8OE/NnctxGTZf/kBS6IuISMrpiEa5fdMmvrJhAwZ8t6SEiydPJlsX6h0Rhb6IiKSUNbH58l9qbub9Rx3FHaWlTMvg+fIHkkJfRERSQnNHB1+uqeF79fUU5OTwq2OO4T/HjtXH8AaQQl9ERAL3u23buLCigtpIhAsmTeIbxcWMHKSIGmj6FxURkcA0RiJcWlXFL7ds4ZihQ/m/44/nHSNHBl1WaCn0RUQk6aLu/KixkS+uX09rNMrXpk/nqqlTydF8+Qml0BcRkaR6be9elpSX83+7d7No1Ch+MHMmMzVfflIo9EVEJClaOzv5em0tN9fWMjw7mx+XlfFJzZefVAp9ERFJuJU7dnB+RQUV+/ZxbkEB354xg3E5OUGXlXEU+iIikjDb29u5av167n/zTYrz8nhm7lzeM2ZM0GVlLIW+iIgMOHfnF7H58re3t/OladP4cmEhQzVffqAU+iIiMqCq9+3jgooKntmxg5OGD+fZefOYm58fdFmCQl9ERAZIezTKd+rrua6mhkFmfL+khAs0X35KUeiLiMgRW7V7N+eVl/PK3r18YOxYvl9SwhTNl59yFPoiInLY9nR0cO2GDXx/0yYm5uTw6DHHcNa4cUGXJX1Q6IuIyGF5fOtWLqqsZFMkwoWTJvH14mJGaL78lKZnR0RE+qUhEuGSykoe2bqVOcOG8cvZs3mb5stPC4GHvpldB5wHbImtusbdn+hlvxpgD9AJdLj7gmTVKCIiXfPl/7ChgS9VV9PmzjemT+eKqVMZrPny00bgoR/zHXe/NY79Frn71oRXIyIiB1jb3MySigr+uns37x49mh/MnMmMIUOCLkv6KVVCX0REUtC+zk6+tnEjt9TVMTI7m58efTTnFhRovvw0lSpjMheb2atmdr+Zje5jHweeMbM1ZrYkmcWJiGSiP+7YwdzVq/l6bS2Lx4/njYUL+bhukJPWzN0TfxCzZ4EJvWxaBvwN2EpXqN8ITHT3T/fyGJPcvcHMxgO/B5a6+5/7ON4SYAlAQUHB/OXLlw9MI31obm4mPwSzTYWlDwhPL2HpA8LTS1j6gL572QXcDTwNTAYuB+Ynt7R+CctzMlB9LFq0aE2f1725e8p8AUXA2jj2uw64Mp7HnD9/vifaihUrEn6MZAhLH+7h6SUsfbiHp5ew9OH+z71Eo1H/aWOjj33uOR+0cqVfs369t3R0BFNcP4TlORmoPoDV3kcmBv6evplNdPfG2OJZwNpe9hkGZLn7ntj3/wrckMQyRURCbf2+fXyuooJnd+zg7SNGcM/MmcwJwatnOVDgoQ/cYmbH0TW8XwOcD13D+cC97n4GUAD8T+x9pEHAz939qWDKFREJj/ZolNvq6rh+40ZyzLirtJTzJ00iS+/bh1Lgoe/uH+9jfQNwRuz7amBeMusSEQm714BL1qzh73v38p9jx/K90lIm5+YGXZYkUOChLyIiybW7o4Nrqqu5C5jc0cFjc+bwH2PHBl2WJIFCX0Qkg/x6yxYurqykoa2Ns4CfnHgiwzVffsZIlc/pi4hIAtW3tnLW2rWctW4dYwcP5m8nnMBSUOBnGD3bIiIh1unO3Zs2cc2GDXS4883iYi6fMoXBWVmsDLo4STqFvohISL3a3MyS8nKe37OHfx09mrtnzqRY8+VnNIW+iEjI7Ovs5IaNG7m1ro7Rgwbx4KxZfHT8eE2fKwp9EZEw+f327XyuooLq1lY+PWECt8yYwVGDBwddlqQIhb6ISAhsaWvjivXreaCpiZlDhrBi3jxOHd3X/cskUyn0RUTSmLvz06YmrqiqYndnJ18uLOSaadPIy84OujRJQQp9EZE0VdnSwucqKvjjzp38y4gR3FNWxuxhw4IuS1KYQl9EJM20RaN8q66OG2tqyMvK4gczZ3LexImaL18OSaEvIpJG/rJrF0vKy1nX0sLZ48Zxe0kJEzVfvsRJoS8ikgZ2dXRwdXU1P2hoYEpuLr+ZM4f3ab586SeFvohICnN3Ht26laWVlTS1tXHZlCncUFREvqbPlcOg/zUiIimqrrWViysreXzbNo7Pz+c3xx7L/OHDgy5L0phCX0QkxXS6c8emTVy7YQNRd26dMYNLJ09mUJbukSZHRqEvIpJCXt6zhyUVFazas4d/GzOGu0pLKdJ8+TJAFPoiIimgpbOT62pq+HZdHUcNHszy2bP58Lhxmi9fBpRCX0QkAA82NbGsupraSIRxgwcTdWdrRwfnTZzIN4uLGa358iUBFPoiIkn2YFMTS8rLaYlGAdjc3o4B106bxo3FxcEWJ6Gmq0JERJLsmurqtwJ/PwceaGoKpiDJGAp9EZEkenHPHmojkV639bVeZKAo9EVEkmBneztLKys5cc2aPn/xTtN0upJgek9fRCSB3J2fb97MFVVVbGlv56LJk5k7bBiXVlUdMMQ/NCuLm/R+viSYQl9EJEFe37uXCysrWblzJwuHD+eJuXM5ITaj3pDs7Leu3p+Wm8tNxcUsLigIuGIJO4W+iMgA29vZydc2buS2ujrys7N7vfXt4oIChbwknUJfRGQAPbZ1K5dUVlIbifCpCRP4ZnEx43Nygi5LBEiRC/nMbKmZlZvZOjO7pY99To/tU2VmX0p2jSIiB7Nh3z7+4+9/5wNr1zJi0CD+fNxx/PjooxX4klICf6VvZouAM4G57h4xs/G97JMN3Am8B6gHVpnZ4+7+WnKrFRE5UCQa5ba6Or62cSNZwK0zZnDJ5MkM1s1xJAUFHvrABcDN7h4BcPfNveyzEKhy92oAM1tO1x8KCn0RCcwfduzgoooKyvft40PjxvGdGTOYkpcXdFkifUqFP0VnAieb2fNm9iczO7GXfSYDdd2W62PrRESSrjES4WOvvca7X3mFDneePPZYfnnMMQp8SXnm7ok/iNmzwIReNi0DbgL+CFwKnAg8BBR7t8LM7Gzgve7+2djyx4GF7r60j+MtAZYAFBQUzF++fPkAdvPPmpubyc/PT+gxkiEsfUB4eglLHxCOXjqBhyIRHszNpR1YDHwUSNd37cPwnID66GnRokVr3H1Bb9uSMrzv7u/ua5uZXQA8Ggv5F8wsCowFtnTbrR6Y2m15CtBwkOPdA9wDsGDBAj/11FMPv/g4rFy5kkQfIxnC0geEp5ew9AHp38vfdu3igspKXgbeO3o0d5SWUjJ0aNBlHZF0f072Ux/xS4Xh/V8DpwGY2Uy6/mje2mOfVUCpmU03sxzgHODxpFYpIhlpW3s7S8rLeftLL7GlrY3rgCfnzk37wJfMlAqhfz9QbGZrgeXAJ93dzWySmT0B4O4dwMXA08DrwMPuvi6wikUk9KLu3N/YSNnzz3N/YyNXTp3K6wsXcgpg3SbZEUkngV+97+5twLm9rG8Azui2/ATwRBJLE5EM9WpzMxdUVPCX3bt558iR3FVayrEheM9YJPDQFxFJFXs6Oriupobb6+sZPXgwPy4r4xMTJhwwfa5IOlPoi0jGc3d+tWULl1VV0djWxpKJE/l6cTFjBg8OujSRAaXQF5GMVtnSwsWVlTyzYwfH5+fz6Jw5nDRiRNBliSSEQl9EMtK+zk5urq3l5tpa8rKy+H5JCRdMnky2hvIlxBT6IpJxnty2jYsrK6lubWXx+PHcOmMGE3Jzgy5LJOEU+iKSMepaW7msqopHt27l6KFD+eO8eSwaPTroskSSRqEvIqHXHo3y3fp6rq+pIQp8Y/p0Pj91Kjm6E55kGIW+iITan3fu5MKKCta1tHDmUUfx3ZISioYMCboskUAo9EUklDa3tXHV+vX8tKmJorw8Hp8zh/ePHRt0WSKBUuiLSKh0unNPQwPXbNjA3s5Olk2bxjWFhQzNzg66NJHAKfRFJDRW797NBZWVrN6zh9NGjeLO0lKOHjYs6LJEUoZCX0TS3s72dpZt2MDdDQ0U5OTw81mzOGf8eN0YR6QHhb6IpC1352dNTVy5fj1b29tZOnkyN0yfzshB+tUm0hudGSKSll7bu5cLKyr4065dnDR8OE/Nncvxw4cHXZZISlPoi0ha2dvZyY01NdxWX8/w7GzumTmTz0ycqDvhicRBoS8iacHdeWzrVi6pqqIuEuHTEyZwc3Ex43Jygi5NJG30O/TNbBjQ6u6dCahHROSfVO/bxyWVlfxu+3aOHTaMX8yezb+MHBl0WSJp55Chb2ZZwDnAYuBEIALkmtkW4AngHnevTGiVIpKRItEo36qt5abaWgaZcduMGSydPJnBmj5X5LDE80p/BfAscDWw1t2jAGY2BlgE3Gxm/+PuP0tcmSKSaZ7dvp2LKiup2LePs8eN49szZjAlLy/oskTSWjyh/253b++50t23A48Aj5jZ4AGvTEQyUkMkwuerqnhoyxZKhgzhqblzee+YMUGXJRIKhwz9/YFvZqV0vdrf5+4X9baPiMjh6ohGubOhgS9v2EBbNMr1RUV8YepU8jR9rsiA6c8bYw8AvwROBjCzOWb204RUJSIZ5a+7drFgzRouq6riX0aOZN3ChXylqEiBLzLA+hP6We7+JNAJ4O5rgTkJqUpEMsK29nbOKy/nHS+9xLaODh455hieOPZYZujWtyIJ0Z+P7DWY2XTAAaxrUmudmSLSb1F3fvzmm3xx/Xp2dXZy1dSpfKWwkHxNnyuSUP05wy4D7gUmmNl/AacDaxNSlYiE1ivNzVxQUcFfd+/m5JEjuau0lDn5+UGXJZIR4g59d68xs9OBDwDzgD8B9yeqMBEJl73A5VVVfK++nqMGD+YnRx/NJwoKdCc8kSSKZ3Iec3cHcPcO4Fexr173ERHpzt15eMsWLgK219dz/qRJ3DR9OmMG65O+IskW1+Q8ZvYI8Ji71+5faWY5wDuBT9I1gc9PDrcIM1sKXAx0AL9z9y/0sk8NsIeuCwk73H3B4R5PRJKjoqWFiyoreXbHDkqBJ044gYUjRgRdlkjGiif0Twc+DfwidiHfTrou4MsCngG+4+4vH24BZrYIOBOY6+4RMxt/kN0XufvWwz2WiCTHvs5Ovl5byy21tQzJyuKO0lKOrqxU4IsELJ7JeVqBu4C7YjPvjaVrgp6dA1TDBcDN7h6JHW/zAD2uiATgd9u2sbSykg2trZxbUMC3iouZkJvLykrdokMkaP26a4W7t7t74wAGPsBM4GQze97M/mRmJ/Z1eOAZM1tjZksG8PgiMgBqW1s5a+1a3vf3v5OXlcWKefN4YNYsJuTmBl2aiMRYvNffmVku8EGgiG4jBO5+Qxw/+ywwoZdNy4CbgD8Cl9J1F7+HgOKeFwaa2SR3b4gN//8eWOruf+7jeEuAJQAFBQXzly9ffsj+jkRzczP5IfjIUVj6gPD0kg59tNN1Ze/+6Tk/AXwI6HmZXjr0Eo+w9AHh6UV9HGjRokVr+rrurT+h/xSwC1hDbFY+AHe/7UiKiz3uze6+Mra8Hnibu285yM9cBzS7+62HevwFCxb46tWrj6TEQ1q5ciWnnnpqQo+RDGHpA8LTS6r38aedO7mwooLXWlr4wNixfLekhMI+7oSX6r3EKyx9QHh6UR8HMrM+Q78/k/NMcffTj7iaf/Zr4DRgpZnNBHKAAy7WM7NhdE0DvCf2/b8ChxxhEJGB8WBTE8uqq6mNRJiWm8sXpk3jb7t380BTE0V5efxmzhzeN3Zs0GWKyCH0J/T/YmbHuvvfB7iG+4H7zWwt0AZ80t3dzCYB97r7GUAB8D+xSTwGAT9396cGuA4R6cWDTU0sKS+nJRoFYGMkwkWVlWQB1xYWcvW0aQzVjXFE0kJ/Qv+dwH+ZWTUQAQxwd597JAW4extwbi/rG4AzYt9X0zULoIgk2bLq6rcCv7sJOTncOH16ABWJyOHqT+ifTizoE1SLiKSg2kik1/WNbW1JrkREjlQ80/Duofeg3/8HgGbbEAmhTnfua2zs8y/9afoonkjaiWdynuHJKEREUseq3bu5sLKS1Xv2UDZkCBsjEVq7DfEPzcripuLiACsUkcPRr8l5RCTctra1saS8nJNefJH6SIQHZ83i9YULubesjMLcXAwozM3lnrIyFhcUBF2uiPRTf97TF5GQ6nTn3sZGrqmuZldHB5dPmcJXi4oYMajrV8TiggKFvEgIKPRFMtwLu3dzUWwo/5SRI7mjtJQ5IZjdTET+mUJfJENtbWvjmg0buLexkQk5OTw4axYfHT+e2HwYIhJCCn2RDHOooXwRCS+d5SIZ5IXdu7mwooI1zc0ayhfJQAp9kQywta2Nqzds4L7YUP7PZ83iHA3li2Qchb5IiHW686OGBq7ZsIHdHR18fsoUvqKhfJGMpTNfJKSe372bi2JD+aeOGsUdpaUcM2xY0GWJSIAU+iIhs38o/97GRibm5PCLWbP4iIbyRQSFvkhodB/K39PZyRWxq/KHayhfRGL020AkBDSULyLxUOiLpDEN5YtIfyj0RdJQz6H8K6dO5SuFhRrKF5GD0m8IkTTzfGyCnRebm1kUG8qfraF8EYmDQl8kTWxpa+Pq6mrue/NNJmkoX0QOg0JfJMV1unNPQwPLNJQvIkdIvzVEUthrwBVr1mgoX0QGhEJfJAVtaWvjS9XV3A9Mamtj+ezZfHjcOA3li8gRUeiLpJBOd34YG8pv7uzkI8CPFi7UUL6IDAj9JhFJEX/btYuLKit5sbmZ00aN4vulpWxetUqBLyIDRr9NRAL21lB+7Kr87kP5m4MuTkRCRaEvEpCeQ/lXTZ3Kl3VVvogkkH67iATgr7Gh/JdiQ/l3lJYyS1fli0iCBR76ZvYQUBZbHAXsdPfjetnvdOB2IBu4191vTl6VIgNjc2wo/8e9DOWLiCRa4KHv7h/Z/72Z3Qbs6rmPmWUDdwLvAeqBVWb2uLu/lrRCRY6AhvJFJBWkzG8c63qp82HgtF42LwSq3L06tu9y4Ey65i4RSWkayheRVJEyoQ+cDDS5e2Uv2yYDdd2W64GTklKVyGHqPpQ/OSeHh2bP5mwN5YtIgMzdE38Qs2eBCb1sWubuj8X2uZuuV/O39fLzZwPvdffPxpY/Dix096V9HG8JsASgoKBg/vLlywemkT40NzeTn5+f0GMkQ1j6gGB76QQeB+4H9gFnA58AhhzGY+k5ST1h6QPC04v6ONCiRYvWuPuCXje6e+BfdI04NAFT+tj+duDpbstXA1fH89jz58/3RJVjrWIAABQaSURBVFuxYkXCj5EMYenDPbhe/rJzpx+/apWzYoW/66WX/LXm5iN6PD0nqScsfbiHpxf1cSBgtfeRiakyvP9u4A13r+9j+yqg1MymA5uAc4CPJas4kUPpOZT/8OzZfEhD+SKSYlIl9M8BftF9hZlNouujeWe4e4eZXQw8TddH9u5393UB1ClygE53ftDQwLWxq/K/OHUq1xYWkq+r8kUkBaXEbyZ3/1Qv6xqAM7otPwE8kcSyRA7wYFMTy6qrqY1EmJaby6cmTOCxbdt4ubmZd48ezfdLSjhaV+WLSApLidAXSXUPNjWxpLyclmgUgI2RCNdv3Mjo7GwN5YtI2lDoi8RhWXX1W4HfXf6gQZw9fnwAFYmI9F9W0AWIpIPaSKTX9fV9rBcRSUUKfZGDaIxE+Pjrr9PXbBbTcnOTWo+IyJFQ6Iv0oi0a5dbaWma+8AIPb97MB446iiFZB54uQ7OyuKm4OKAKRUT6T6Ev0sOz27czb/Vqrqqu5pSRI1l34on8z7HH8qOyMgpzczGgMDeXe8rKWFxQEHS5IiJx04V8IjG1ra18vqqKR7ZupTgvj9/MmcP7xo59a/viggKFvIikNYW+ZLzWzk5uq6/npo0bAbixqIgrp04lLzs74MpERAaWQl8y2u+2bePSykrWt7bywbFjua2khMK8vKDLEhFJCIW+ZKT1+/ZxWVUVv922jaOHDuWZuXN5z5gxQZclIpJQCn3JKC2dnXyjtpZbamvJycriW8XFXDJlCjlZuqZVRMJPoS8Zwd15ZMsWrli/ntpIhHMLCvhmcTGT9Dl7EckgCn0Jvdf37mVpZSV/2LmTucOG8bNZszh51KigyxIRSTqFvoTW7o4Obqip4fZNm8jPzub7JSV8btIkBmkoX0QylEJfQsfd+T3wsRdeoLGtjc9MmMDXi4sZn5MTdGkiIoFS6EuovNLczMWVlTwHLMjN5ddz5rBwxIigyxIRSQkKfQmFHe3tfKWmhrs2bWL0oEFcCXzzhBPI0j3uRUTeojc3Ja1F3bm3oYGZL7zAXZs2ccGkSVScdBL/Dgp8EZEe9Epf0tYLu3dzcWUlq/bs4V9GjOCOuXM5bvjwoMsSEUlZCn1JO1va2ri6upr73nyTCTk5PHD00SwuKMD0yl5E5KAU+pI2OqJRftjYyLUbNtDc2ckVU6bwlaIiRgzSf2MRkXjot6Wkhed27uTiykpe2buXd40axfdKS5k9bFjQZYmIpBWFvqS0xkiEL1RX87OmJqbm5vLL2bP54LhxGsoXETkMCn1JSe3RKN/btInrampoi0ZZNm0aVxcWMkz3uBcROWwKfUk5z27fziVVVbze0sIZY8Zwe0kJJUOHBl2WiEjaU+hLyqhtbeXzVVU8snUrxXl5/GbOHN43dmzQZYmIhIZCXwLX2tnJbfX13LRxIwA3FhVx5dSp5GkoX0RkQAUe+mb2EFAWWxwF7HT343rZrwbYA3QCHe6+IGlFSsL8bts2Lq2sZH1rKx8cO5bbSkoozMsLuiwRkVAKPPTd/SP7vzez24BdB9l9kbtvTXxVkmjr9+3jsqoqfrttG0cPHcozc+fynjFjgi5LRCTUAg/9/azrM1gfBk4LuhYZOA82NbGsupraSIRpubl8taiI6tZWbqmtJScri28VF3PJlCnk6B73IiIJlzKhD5wMNLl7ZR/bHXjGzBz4obvfk7zS5HA82NTEkvJyWqJRADZGInymvBwHFo8fzy0zZjApNzfYIkVEMoi5e+IPYvYsMKGXTcvc/bHYPncDVe5+Wx+PMcndG8xsPPB7YKm7/7mPfZcASwAKCgrmL1++fCDa6FNzczP5+fkJPUYyDHQf5wBNvawfDTw6YEfpnZ6T1BOWXsLSB4SnF/VxoEWLFq3p67q3pIT+oZjZIGATMN/d6+PY/zqg2d1vPdS+CxYs8NWrVx95kQexcuVKTj311IQeIxkGuo+slSvp7X+XAdEE/3vpOUk9YeklLH1AeHpRHwcysz5DP1XeSH038EZfgW9mw8xs+P7vgX8F1iaxPumHqDv//eab9DVR7jQN6YuIBCJVQv8c4BfdV5jZJDN7IrZYADxnZq8ALwC/c/enklyjxGHNnj2886WX+NQbb1CUm0tejznyh2ZlcVNxcUDViYhktpS4kM/dP9XLugbgjNj31cC8JJcl/bC1rY1rNmzg3sZGxg0ezP1lZXxywgR+sXnzAVfv31RczOKCgqDLFRHJSCkR+pK+OqJRftDQwJdratjT0cGlU6ZwXVERI2P3uF9cUKCQFxFJEQp9OWx/3rmTpZWVvLp3L6fF7nF/jO5xLyKSshT60m/1ra1cVV3N8s2bmaZ73IuIpA2FvsQtEo3ynbo6vrZxIx3ufLmwkC9Nm8ZQ3RhHRCQtKPQlLk9s28ZlVVVU7tvHmUcdxbdLSigeMiToskREpB8U+nJQVS0tXL5+Pb/dto2ZQ4bw1Ny5vFc3xhERSUsKfenV3s5Ovr5xI7fW1ZGTlcUtxcVcqhvjiIikNYW+HMDdeXjLFq5cv576SIRzCwr4ZnGxbowjIhICCn15SzVw/SuvsHLnTo7Lz2f57Nn8y8iRQZclIiIDRKEv7Ghv56s1NdwJjGpu5u7SUs6bNIlsfQRPRCRUFPoZLOrO/Y2NXL1hA9vb23kfcP9JJ3HU4MFBlyYiIgmgq7Iy1PO7d3PSiy9yXkUFZUOGsHr+fC4HBb6ISIjplX6GaWpr40vV1fzkzTeZmJPDz2bN4mPjx2NmrAy6OBERSSiFfoZoj0a5Y9MmrqupYV80yhemTuXawkKGD9J/ARGRTKHf+Bngjzt2sLSyktdaWjh9zBi+W1JC2dChQZclIiJJptAPsdrWVq5Yv55fbdnC9Lw8Hpszh/cfdZRujCMikqEU+iHU2tnJt+rq+EZtLQA3FhVx5dSp5OnGOCIiGU2hHyLuzuPbtnF5VRUbWls5e9w4bp0xg2l5eUGXJiIiKUChn6YebGpiWXU1tZEI03JzuXjyZJ7dsYOnd+zgmKFD+cO8eZw2enTQZYqISApR6KehB5uaWFJeTks0CsDGSISrqqvJM+O7JSVcOGkSg3VjHBER6UGhn4aWVVe/Ffjdjc3J4dIpUwKoSERE0oFeDqah2kik1/Wb+lgvIiICCv20srWtjfPLy/E+tk/T7W9FROQgFPppoCMa5Y76ekpfeIH7Ghs5ffRohvR4z35oVhY3FRcHVKGIiKQDhX6KW7ljByesWcPSqirm5+fz6okn8uS8efyorIzC3FwMKMzN5Z6yMhYXFARdroiIpDBdyJei6lpbuXL9eh7esoXC3FweOeYYzho79q3Z9BYXFCjkRUSkXxT6Kaa1s5Nb6+r4em0tDlxfVMRVU6cyRLPpiYjIEQo89M3sOOAHQB7QAVzo7i/0st/pwO1ANnCvu9+c1EITrOdseh+KzaZXqNn0RERkgAQe+sAtwPXu/qSZnRFbPrX7DmaWDdwJvAeoB1aZ2ePu/lqyi02EN/bu5dKqKp7RbHoiIpJAqRD6DoyIfT8SaOhln4VAlbtXA5jZcuBMIK1Df3dHBzfU1HD7pk0My8ri9pISLtBseiIikiCpEPqXAU+b2a10fZrgHb3sMxmo67ZcD5yUhNoSIurOT998ky9VV7O5vZ3PTJzITdOnMz4nJ+jSREQkxMy9r6leBvAgZs8CE3rZtAx4F/And3/EzD4MLHH3d/f4+bOB97r7Z2PLHwcWuvvSPo63BFgCUFBQMH/58uUD10wvmpubyc/Pj2vfN4DvAa8Ds4GlwNGJK61f+tNHqgtLL2HpA8LTS1j6gPD0oj4OtGjRojXuvqDXje4e6Bewi3/88WHA7l72eTvwdLflq4Gr43n8+fPne6KtWLHikPs0RSL+mddfd1uxwguee87/u7HRO6PRhNfWH/H0kS7C0ktY+nAPTy9h6cM9PL2ojwMBq72PTEyF4f0G4BRgJXAaUNnLPquAUjObDmwCzgE+lqwCj0R7NMqdmzZxXU0Ne6NRrpg6lS8XFjJiUCr804uISCZJheQ5D7jdzAYBrcSG5c1sEl0fzTvD3TvM7GLgabo+sne/u68LrOI4/WHHDi6prOS1lhbeO3o03y0p4ehhw4IuS0REMlTgoe/uzwHze1nfAJzRbfkJ4IkklnbYavbt44r163l061aK8/J4bM4c3n/UUW/NpiciIhKEwEM/TPZ1dnJLXR0319aSBXxt+nSumDKFPM2mJyIiKUChPwAceGTLFq6oqmJjJMJHxo3jWzNmMFWz6YmISApR6B+hdXv3ciXw4rp1HDtsGCtnzeKUUaOCLktEROSfKPTj9GBTE8uqq6mNRJiWm8uyadNY19LCHZs2MRS4o7SU8ydOZJBm0xMRkRSl0I/Dg01NLCkvpyUaBWBjJMKSyq5PFn5u0iROb2jgzMmTgyxRRETkkPSyNA7LqqvfCvzuJubkcPfMmYwMoCYREZH+UujHoTYS6XX9m21tSa5ERETk8Cn04zAtN7df60VERFKRQj8ONxUXM7THBXpDs7K4qbg4oIpERET6T6Efh8UFBdxTVkZhbi4GFObmck9ZGYsLCoIuTUREJG66ej9OiwsKFPIiIpLW9EpfREQkQyj0RUREMoRCX0REJEMo9EVERDKEQl9ERCRDKPRFREQyhEJfREQkQyj0RUREMoS5e9A1JJSZbQE2JvgwY4GtCT5GMoSlDwhPL2HpA8LTS1j6gPD0oj4OVOju43rbEPrQTwYzW+3uC4Ku40iFpQ8ITy9h6QPC00tY+oDw9KI+4qfhfRERkQyh0BcREckQCv2BcU/QBQyQsPQB4eklLH1AeHoJSx8Qnl7UR5z0nr6IiEiG0Ct9ERGRDKHQP0xmttTMys1snZnd0sc+o8zsV2b2hpm9bmZvT3adhxJPH7H9ss3sJTP7bTLr649D9WJmU81sRey5WGdmlwZR56HE+X/r9Ng+VWb2pWTXGA8zu87MNpnZy7GvM/rY7/JYr2vN7BdmlpfsWg+mH32k9Pkebx+xfVP6fI+nl3Q43/vxf2vAzvdBR/LDmcrMFgFnAnPdPWJm4/vY9XbgKXf/kJnlAEOTVmQc+tEHwKXA68CIpBTXT3H20gFc4e4vmtlwYI2Z/d7dX0tqsQcRTx9mlg3cCbwHqAdWmdnjqdRHN99x91v72mhmk4FLgNnuvs/MHgbOAX6SpPriddA+YlL6fI+Jpw9I8fM95lC9pPz5HnOoc2RAz3e90j88FwA3u3sEwN0399zBzEYA/w+4L7ZPm7vvTGqVh3bIPgDMbArw78C9Saytvw7Zi7s3uvuLse/30PVLbXJSqzy0eJ6ThUCVu1e7exuwnK4/FNLVIGCImQ2iKygbAq6n39LkfI9Lmpzvh5Qm53s8BvR8V+gfnpnAyWb2vJn9ycxO7GWfYmAL8OPYMNm9ZjYsuWUeUjx9AHwX+AIQTV5p/RZvLwCYWRFwPPB8Emrrj3j6mAzUdVuuJ3V/mV1sZq+a2f1mNrrnRnffBNwK1AKNwC53fybZRcbhoH2QHuc7HLoPSI/zHeLrBUjp8x0O3ceAnu8K/T6Y2bOx9xh7fp1J1yuT0cDbgKuAh83MejzEIOAE4G53Px7YCyT9vdcj7cPM3gdsdvc1ya69pwF4TvY/Tj7wCHCZu+9OWgP/OP6R9tFbX4F8DOcQvdwNzACOoyvQb+vl50fT9aplOjAJGGZm5yaxhf11HFEfpMf5Hs/zkS7nezzPyf7HSeXzPZ4+BvZ8d3d99fMLeAo4tdvyemBcj30mADXdlk8Gfhd07YfRxzfo+suyBngTaAF+FnTth9NLbP1g4Gng80HXfATPyduBp7stXw1cHXTth+irCFjby/qzgfu6LX8CuCvoeg+jj5Q/3+PsIy3O93h6iW1L6fM9zudkQM93vdI/PL8GTgMws5lADj1ukuDubwJ1ZlYWW/UuINUuIImnj6vdfYq7F9F1gdUf3T3pr8TicMheYq+Y7wNed/dvJ73C+ByyD2AVUGpm063rgrFzgMeTWmUczGxit8WzgLW97FYLvM3Mhsaen3fR9d5ryoinj3Q43+PsIy3O93h6SYfzPc5zZGDP96D/uknHL7p+Ef8s9gS9CJwWWz8JeKLbfscBq4FX6fplPjro2g+nj277nwr8Nui6D7cX4J10DYu9Crwc+zoj6NoP8//WGUAFXSMBy4Kuu49eHgD+Hvv3fhyY2Ecv1wNvxHp+AMgNuvbD7CPVz/e4+ui2fyqf74fsJU3O93j/bw3Y+a4Z+URERDKEhvdFREQyhEJfREQkQyj0RUREMoRCX0REJEMo9EVERDKEQl9ERCRDKPRFREQyhEJfRPrFzJrj2GdI7IZB2bHld5vZAz32yTGzP1vX3fVEJAkU+iKSCJ8GHnX3ztjyPOCl7jt4121C/wB8JMm1iWQshb6I9JuZFZnZ62b2IzNbZ2bPmNmQbrssBh7rtjwPmGBm/2tmb5rZu2Prfx3bV0SSQKEvIoerFLjT3Y8BdgIfhK5he6DY3Wu67TsP2OruJwMX8o+gXwucmLSKRTKcQl9EDtcGd3859v0aum4NCjCWrj8CADCzwcAY4NbYqkH7t8eG/9vMbHgyChbJdAp9ETlckW7fd9IV5gD7gLxu22YDr7h7NLY8lwNvIZoLtCaqSBH5B4W+iAwod98BZJvZ/uCfB7zSbZe5dN1KFDM7Ctji7u3JrVIkMyn0RSQRnqHrfubQFfqvdts2h3+80l8EPJHEukQymrl70DWISMiY2fHA593944fY71HgancvT05lIplNr/RFZMC5+0vAiv2T8/QmdpX/rxX4IsmjV/oiIiIZQq/0RUREMoRCX0REJEMo9EVERDKEQl9ERCRDKPRFREQyhEJfREQkQyj0RUREMsT/B7vvumoOg47JAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "%reset -f\n", "\n", "from matplotlib.pylab import *\n", "from scipy.optimize import fsolve\n", "\n", "# def RK(tt,phi,y0): # comme suite geometrique\n", "# x=beta*h\n", "# q=x**2/2-x+1\n", "# uu = [y0]\n", "# for i in range(len(tt)-1):\n", "# uu.append( q*uu[i] )\n", "# return uu\n", "\n", "\n", "a,b,c = 1/3, 1/3, 1/3\n", "\n", "\n", "def RK(tt,phi,y0):\n", " uu = [y0]\n", " h = tt[1]-tt[0]\n", " for i in range(len(tt)-1):\n", " K1=phi(tt[i],uu[i])\n", " K2=phi(tt[i]+h/2,uu[i]+h/2*K1)\n", " K3=phi(tt[i+1],uu[i]+h*K1)\n", " uu.append( uu[i]+h*(a*K1+b*K2+c*K3) )\n", " return uu\n", "\n", "\n", "t0, y0, tfinal = 0, 1, 1\n", "beta=50\n", "\n", "phi = lambda t,y : -beta*y\n", "sol_exacte = lambda t : exp(-beta*t)\n", "\n", "figure(figsize=(10,5))\n", "N=int(50/2)+1\n", "tt = linspace(t0,tfinal,N+1)\n", "h = tt[1]-tt[0]\n", "\n", "if h>=2/beta:\n", " print('Attention: h > limite de stabilité')\n", " \n", "yy = [sol_exacte(t) for t in tt] \n", "# uu = RK(tt)\n", "uu = RK(tt,phi,y0)\n", "plot(tt,yy,'b-',label=(\"Exacte\"))\n", "plot(tt,uu,'ro',label=(\"RK\"))\n", "title(r' $N$={}, $h$={}'.format(N,h))\n", "xlabel('$t$')\n", "ylabel('$u$')\n", "legend() \n", "grid(True)\n", "\n", "\n", "H = []\n", "err = []\n", "N = 50\n", "for k in range(7):\n", " N+=100\n", " tt = linspace(t0, tfinal, N + 1)\n", " h = tt[1] - tt[0]\n", " yy = [sol_exacte(t) for t in tt]\n", " uu = RK(tt,phi,y0)\n", " H.append(h)\n", " err.append(max([abs(uu[i] - yy[i]) for i in range(len(yy))]))\n", " \n", "omega=polyfit(log(H),log(err), 1)[0]\n", "figure(figsize=(8,5))\n", "plot(log(H), log(err), 'c-o')\n", "xlabel('$\\ln(h)$')\n", "ylabel('$\\ln(e)$')\n", "title(r'Ordre $\\omega$={:1.2f}'.format(omega))\n", "grid(True);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercice 2 : étude d'un schéma multipas\n", "\n", "On notera $\\varphi_k\\stackrel{\\text{déf}}{=}\\varphi(t_k,u_k)$.\n", "\n", "Soit la méthode multipas\n", "\n", "$$\n", "u_{n+1}=u_{n-1} +h\\left(b_{-1} \\varphi_{n+1}\n", " +b_0\\varphi_{n}+b_{1} \\varphi_{n-1}\\right).\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ">**Q1 [2 points]** \n", "Pour quelles valeurs des paramètres $b_0, b_1, b_2$ est-elle zéro-stable?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Le premier polynôme caractéristique est\n", "$$\n", "\\varrho(r)=r^{p+1}-\\sum_{j=0}^p a_jr^{p-j}=r^2-a_0r-a_1=r^2-1=(r-1)(r+1)\n", "$$\n", "dont les racines sont \n", "$$\n", "r_0=1,\\quad r_1=-1.\n", "$$\n", "La méthode est donc zéro-stable pour tout choix des paramètres car\n", "$$\n", "\\begin{cases}\n", "|r_j|\\le1 \\quad\\text{pour tout }j=0,\\dots,p\n", "\\\\\n", "\\varrho'(r_j)\\neq0 \\text{ si } |r_j|=1.\n", "\\end{cases}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ">**Q2 [3 points]** \n", "Pour quelles valeurs des paramètres la méthode est convergente d'ordre 4?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* Une méthode consistante d'ordre $\\omega\\ge1$ et zéro-stable est convergente d'ordre $\\omega$. \n", "* On a $p=1$: c'est une méthode à $q=p+1=2$ pas avec $a_0=0$ et $a_1=1$. \n", "* La méthode est implicite si $b_{-1}\\neq0$.\n", "\n", "La première barrière de Dahlquist affirme qu'un schéma implicite à $q=2$ étapes consistante et zéro-stable peut être au mieux d'ordre $\\omega=q+2=4$. \n", "\n", "Pour que la méthode soit consistante d'ordre 4 il faut que\n", "$$\n", "\\begin{cases}\n", "\\displaystyle\\sum_{j=0}^p a_j=a_0+a_1=1,\n", "\\\\\n", "\\displaystyle\\sum_{j=0}^p (-j)a_j+\\sum_{j=-1}^p b_j=1,\n", "\\\\\n", "\\displaystyle\\sum_{j=0}^p (-j)^{2}a_j+2\\sum_{j=-1}^p (-j)^{1}b_j=1,\n", "\\\\\n", "\\displaystyle\\sum_{j=0}^p (-j)^{3}a_j+3\\sum_{j=-1}^p (-j)^{2}b_j=1,\n", "\\\\\n", "\\displaystyle\\sum_{j=0}^p (-j)^{4}a_j+4\\sum_{j=-1}^p (-j)^{3}b_j=1\n", "\\end{cases}\n", "$$\n", "Dans notre cas cela reviens à imposer\n", "$$\n", "\\begin{cases}\n", "0+1=1\n", "\\\\\n", "{}-1+\\big(b_{-1}+b_0+b_1\\big)=1\n", "\\\\\n", "1-2\\big(-b_{-1}+b_1\\big)=1\n", "\\\\\n", "{}-1+3\\big(b_{-1}+b_1\\big)=1\n", "\\\\\n", "1-4\\big(-b_{-1}+b_1\\big)=1\n", "\\end{cases}\n", "$$\n", "On cherche donc $b_{-1}$, $b_0$ et $b_1$ tels que\n", "$$\n", "\\begin{pmatrix}\n", "1&1&1\\\\\n", "2&0&-2\\\\\n", "3&0&3\\\\\n", "4&0&-4\n", "\\end{pmatrix}\n", "\\begin{pmatrix}\n", "b_{-1}\\\\b_0\\\\b_1\n", "\\end{pmatrix}\n", "{}=\n", "\\begin{pmatrix}\n", "2\\\\0\\\\2\\\\0\n", "\\end{pmatrix}\n", "$$\n", "Il s'agit d'un système surdeterminé qui admet l'unique solution\n", "$$\n", "b_{-1}=b_1=\\frac{1}{3},\\qquad b_0=\\frac{4}{3}\n", "$$" ] }, { "cell_type": "code", "execution_count": 66, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOIAAAAyCAYAAABBLgf4AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAJEklEQVR4Ae2c6ZXUOBDHh3kTAEcEO2TAEQGQwbIbwUIG8PjGNx5kMGwEHBmwE8EsZABEAEwGs/+fV7Lls9VW222rq95z67COqr9UUulwX7u6ujoaopcvX57q/Rlp5H80lNbeGQKGQIWA9OWOQud6Xsn/pnrT9h23o6oYZX6m0Fc93/Q8rt6YzxAwBDYhIP35ojTozQv5v+phUuuka30zojK9Vg4U8an8bztzryhSMjA6IcvTFbG9U1YdBh9U6F35L3da+IEUNgZD5bkueJgZUUSwZ2Kr0Ukt5AJK+FBelPCt/KtXQicWHZAR6qDIdYK/JfRPPff09I7KBwXMFsKmYqj8l3oeqMrveljmtZZ4fabpc/hU5ixmD8nBoHKQJNnpBI9dW747SBAShd4FhpQhNpjUHsqPdVajPkVkRvynlnKlASc0IJgpttI2zIhtPxCiXzVqKaI6LvYslEvH/VMy5WJe/98y9rtWBPza8HZTgJYiBglYU6yapICYpMXRy6oFMeazR2BIEVctvJTwVAKwPvKj0KrlMebzRiBbRVSzZXHsknf3M+k8AlkqombBJxLQTFLfyuYuHoHsFNGZpNfNJF183zMGAwQ6D/SD92v0sja8L0XkAD8kzm5OXfw3ucVZaZjA/IbAvhDIThGlYJx/ts5AFf+LeLl2Z3Zfvc3q7UUgO9O0V9KjI85H/RnpQLKsX91y0t3MWspphZsEw+xmxGYbaAZk0wZzFeJ6ESbrhdzBz1KK1Jn8OJmRxt/o+KA4jnU+ybXLDhHtPDWGSYoo5phh6NhcJv6pcOvGQISMkyYRT7Pcl10yFuJtVnN8yViM7WxTY5hkmoo5Dsy5Sc7o2lqXjRV6jfkMi6rVDIsKi1hfkiIGlbAj+SkIH7LXsKha37CosAh9rTV6lyKeuhxRl741+vl1x0HPiGBmWFR9zbCosPA+YeJ1qrVp2KWIXrFiZ7jCNA0q8fUeomtYVK1uWFRYhD4+Tvc6Vsa3/ipDCsV5W/TGi9J/VnrWiBeu1Pty3yn+owsnOSqH/8zh/G+WTZcUZg2LCr2psahq2s637/6k+n8Xx2xw1u5C12ZEJfL3MxnNNpLSM8WyDsCc5W81OBL4Sw/b497EVXAcufIph13ZRZNhUTXP1FhUNW3nW0J/Eg9MUBwZnclf6sixF0WRKOEfeh7IH/vpkJ9iyVPYv87Fj+YXpLjX7nmG6+M3ua6sG3Lvbkq7gPexWHAP9okeLIloyhGLaOE7EgqPrXFcCobiA+sOZfwsPxPZ0QkCEaEH2/U3x6y8UcTM+aUjD2UWNxD0ji8hfsgtDtDlcqjOQXLsrFsoeBQ3SqRyh/+oNShIaa8FwVRvDBaA7hUWjLYi8ZsNFlsJ3kgsHEbjuBQMxQemKfswWI9nJ4GMY77Ip1PV1oIq1Hc0v2ZkBuQfrArSe9Z7KCIXsGNnXp99o6syd6lcG+sLEmzEQrwx2DFwldZCkH/n3iVjkSLsnDjOgOFNsGBGZJS9LRfT9LtczEw6TAxh43qF8+n9LQ4UjveM/E2Fo05GtWa8olZLg1isVqpxjG/EQn0DS4klB9bYv3r8AM6XM/zrHH9zQj+hPz1XOLZPKvnySfKwYYPMhb4de5b1Arv1vZ5z+QEylkplUj6UjnUm0y4gEu4iZt9iJOh6uY84x3tq1UNYpJY9W/6psQjKpzNiMd1T3Bse+Vm6YK4RZh2F+fZCTzYkuZAZq6ic9EpFREolQBlREgCKIczScEeTfO8dgEP5UcI+JS3zqRwW5Fd6ttrYKAuI9FCPkv5KrGcsFlFcwqOeXLCg/Rn0C6tIcoUXz2mLVwEo3F8uB7ggfmvvXBgOMSYemAWZ7Zmsylm+poiuAEzUOy6Di+p1OKp4pLTsipKPkQxl9nTpPQ0XsDeCq7LITzpMl8koqIfODm9jaBMWY8os8wQ8rh4LycKH2bRtcebshVQcinkpt+ygCtNxL3yaFNfVOXl/2sCj//PucPA5OunI5K+qsQvo/R3JihkUMP2asJVGgnvAMXVDcEnbDLfyE6EyRn3RoXzU6QcFlIswA0anTNSjZ/QmivIOYqG6kwkexxSifEvFAiULZz/C5UDj+GZT76MeJoeoPjOEkcqYBcMBHpCxJUeXIvqZis67CwLosnIBgR9gfT27qKNWhsqGdxb4XhGP5EfJ2K1lIwAzsovYKOh715V+8XGSZ5FYiC8GByyQsFMyQ4ZXK+krfuCkLcv2lH82SsCwxqMrh7hW3z+upZwgoMpZgN+Sy2yEbcwMihk3JbEj90T1hTOcb9DOhb/S0mF/TMmU6mB0Z6EOD4VfYTCZkhaJhQRGEZuDXjOONsNUBaPyP4b2gOPWGIrfIfrZfNk1IzbTJIcFXAlicmFxBTDKXrqnyCEeaNDC3/OD4jJoTEYqn5HQsBAIwgIl84NjgbniajeoHF6tpY+LnxPHMf2pkCn2ZxZFjGVmV+lcI98Iy1Ocnx3ZVGqR3k+qhK0KZ4owLNKBHoPhtrVObppuy9AU6QUkaw1MwtqN9ynqWnqZhkV6C02BYZYzoodagLEdjhKyCYB5Ue7IyX9QZFikN/eUGHZ9j8imBd8k8lnTXnap0iFrlyBZME25cDC0a9rOmGGMYZHeqGMwVJ5e3ToI0xTYBQI7dGzgcOkAQA6WDIv0pt81hlkqokDi8BeztEneNMVcPQgyLNKbeQ4Ms1REQc/dVD66POiZz3VBwyJdFyfHMNfNGkxQPsPCDemeC9TOr8IEGfoNi/RGnRzDXBWxddgrpWSz5roe/4lWevOsowTDIr2dJsdwaNeU+6CtWw3pMs1TgnhnHRjyf6pw76XvebjaTy2GRTruu8BQZdAH+VfC1olESxFhWRn43xdMO77AMDIEDIEdICB9YgOR9SYfPddm2b7NGtZQD5UYU87IEDAEdoOA361/1yyuTxG9tnZ+qdAsxMKGgCEwjICb1NAn/jyMW1416lREl5BL0PwPqb8sXctoAUPAEIhDwCnhuUsd7luUBXQqIm+VmVmRp/jfReKMDAFDYDsEpEeYo6wLWebdVbj1UTAldm7W8MKTMp7KX3w6JL9t3nhgzDUENiAgfWFzhpnwlfyDn9n9Bz9EOpmdf7OzAAAAAElFTkSuQmCC\n", "text/latex": [ "$\\displaystyle \\left\\{ b_{0} : \\frac{4}{3}, \\ b_{1} : \\frac{1}{3}, \\ b_{m1} : \\frac{1}{3}\\right\\}$" ], "text/plain": [ "{b₀: 4/3, b₁: 1/3, bₘ₁: 1/3}" ] }, "execution_count": 66, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import sympy as sym\n", "sym.init_printing()\n", "\n", "sym.var('b_m1,b_0,b_1')\n", "\n", "eq2 = sym.Eq(b_m1 + b_0 + b_1, 2)\n", "eq3 = sym.Eq(2 * b_m1 - 2 * b_1, 0)\n", "eq4 = sym.Eq(3 * b_m1 + 3 * b_1, 2)\n", "eq5 = sym.Eq(4 * b_m1 - 4 * b_1, 0)\n", "sym.solve([eq2, eq3, eq4, eq5])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ">**Q3 [3 points]** \n", "Vérifier empiriquement l'ordre de convergence sur le problème de Cauchy\n", "$$\n", "\\begin{cases}\n", "y'(t) = y(t)\\ln(1+y(t)), &\\forall t \\in I=[0,1],\\\\\n", "y(0) = 2,\n", "\\end{cases}\n", "$$\n", "après en avoir calculé la solution exacte." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- On définit la solution exacte (calculée en utilisant le module `sympy`) pour estimer les erreurs:" ] }, { "cell_type": "code", "execution_count": 81, "metadata": {}, "outputs": [], "source": [ "%reset -f\n", "%matplotlib inline\n", "\n", "from matplotlib.pylab import *\n", "\n", "# variables globales\n", "t0 = 0\n", "tfinal = 1\n", "y0 = 2\n", "\n", "phi = lambda t,y:y*log(1/y)" ] }, { "cell_type": "code", "execution_count": 82, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAANoAAAAyCAYAAAA6CsU0AAAACXBIWXMAAA7EAAAOxAGVKw4bAAANHklEQVR4Ae2d7XUdNRCGb3xcgDEdQAcBKiDpAEIFQAfh+F/++ZAOIBUk0AFQQYI7IB3EuIPwPopGaLVarbT37v1wNOestSuNRprRjGak1V4/eP/+/aZDl8AxSeDZs2d/qD8/Kb05pn6lfVH/LpT3m65vdX+XlsfPZ/FDv+8SOLQEpLAo7m/HbmTIyRvXT7r9W/cY3SScT5ZkCkTsobJf6PpM1ys9/5hB61ldAoskIH3CyN4q/XURgUIlr7vQ/0L3Re9TIDMqEq0bXc7YVPj5CMFnPFgSOoow8Sbu8vcpwj2/S6BFAtKlH4T/o9IvWuqVcEULL4NjuNX1pS4cxSfK35mhiZ4D0fxZN58p/dZnDZLm0FGEHnkKfw4o9YcugYUSkE4RIaGoWSVdSHYjune6cAhEXi+X0qmppzbwag+VMmGMoNnQROGxLtz7zmeFUe96xsciATY/rqVTb0+cYYztZ/GBJx3AEkPDo3VvNhBjf1gqAe8BLpU+X0rjWOqJB5ZSTBZ45wEUN0NU0Vw6ld/pwsCIc691degS2IUEUMqdb37somMLaeDV/pDt4NmCh570aELCc/2ti/cZXMw4Zqndo0kYHbaTgHTqqSgQZt2biVs8YRu8/zNb0e1mkzU0IcM8W6EYWLBKPXPPdmZfn0kQHbaWAJsUv99DfWLi+MbbkRNS1tBUgjVeCDF16X195sTW/2wrAekWSxCWJvfGm0UysYgv7EBOGdoTVTJkV99bJ4Jhh6hDl8C2ErgSAbbfj/qY1RImxRMRH/aDx3YwMjRvUISOqUFhfBuVDwzQUel/ugTaJXDfoyOWXrzAxjnl12heZvHajCzen7nZR5WJPx0Bj9uTLoFqCUh3CBtzk3k1jRNAfOX7yIQyNjQJwdxeMCQMS7gI5w2VBI+Vlxrih5Id/RX90H6JZC1eicYxlbXw04J7TDyqL075lO4zOvrUy+ByH7LQ2GBHXDiosaGRKeAozFdCfqqLjRGMigq8WGRL9hddq4FvA8OuAdwzfTp5aOQbfk+Vd3SJ9dmqkzUCUht8CUAYZxsT7ll59gzaWoBjcnq86FDxWr2CrgSA98TIefFXBcLF0Bi4dJe0qv4xIC3hm37X8i48BpwDtkQKf+p5p+cKRbMa1Pa/Qn6j1M321RVPDFH84aRwVg+OytDUIeL2v5RmT3D7cl6ij0JXlZH/tVLc9UmB52uSb5jxOFvzLjrQ4KzqQQzN84Gh/ar7sCsHj/cNxB9ek+jv8dmRMccMUApLie2ZkfnsIQXqUf8UYY5veNoV7znZ7VNmtvb+Z5+NHqgt29N4eGyG9kSzQCn8I9TIfjng61Efr3hqMMc3/NwX3s3QVl+fHYESGI+fH42hyUBYm1nHpmTErM4J6SmgvnvfN4VwbPmVfNPto+ZdfJgBzYnYdv3mxnqOztGXSya2jLk8X6u3agTPQozK5938pkLwVL7shdJ4ncCMPdruFQ4KRiwPPQaTj+vYRXqtNP20gvrQCW3pfu+gfrXwnuWbTovO3nj3fSaEtZCOcWOHbjAmUZ8wFMaDgw3wS6TBqfXJz/mFA4AL3H1IPoq/F6sZmsR3JaFzKBlPxW5XrPx4HfJj+FIPo/WZ6jPQ7JKB/0gpijkFKAnKWQWiRXvV+J4oh6rjCSLXVgvvWb4hqnZW4z3utNphR5LJa7DJpHwMh4nNTWgeD8MKPwegPDY2GOfaw8H2PuvQa0V1ey/AhLKOR/MD8tqzgWGkQs3N4sx0KZ4n4RLq3MQZmXvq14YwG/Vz57teC3if4xs2d857IjuMDENJwzlesRCNMNEh+ytd6RqZyQA8dhEH3k95OYDfjwnQyYtzCWdnP+woWg+8BBkMMwq813UiWbxImkfsfpfgxY/UKa3PwEVRDj2QrbzP8Q1fq/GuccKbMTnZxEh7DhhDXdx/p8vGk+dtAH43olsa6xF94e9MT0fEGzPUF9Pz2pqXGFprpVniJkSlKAiKH8JG5TGw5NXMfkJzgwK+rQdc3sQfBrFpACfoLM4+Qd4tAijJjTEDCLUJJ/mEyvApW/21itrbuZ7C0L7gfOWGWMukH4pifLnPI5yLnegPdTYSdjBO3ceDbdUwSOhUgWistUaj/VreS3xDZxXeIezBwkVkNwWGw3shJk02ssjDSPlkP0ykep4DNz4T4zdX91TLb9c2NAbCBsmElFufUWYDZ3hxOlijaJDYGAH/JkbSPR4tbS9B+f9RdHa+RvufulPCtC853kt8Q24V3q2fkoFNhLQzCM1V5oxc+azhAJ5HO5GupP6PecKDRx/1Xd4K0+nk2VYk5isPFM0PHIPFzlUKGM1XaaZ/DgYkGsy8HKZNjQxUjm7l8inbN9TyXuKbPu+D96/VDi/NHyZCIiR8rnyLJDASPBg7ke5bK6UlT5iQc4/vfGZrvRyt6jz6W4vcgltJ8+3aHo3dKMIMBgwB2zsWG7i4ny/1YDNnnM89dH4RHXfiWmn6/szwMeI1vZS1U5PW8l7im3Z2wrtkhhGxa8jrhI2ekbX7TRjd49WYpHjGmAAUk+cwVtzroozzkgGUx6TC+EyNS8DVjdGP81a9V784dE4fB5NfoVEmEb65rOGnQMYVMaHcFg8VqyEGh3dgCH3r39oXPQzukdKpQ8O8B+OXZZu9kurQRxbqZsx6PB5QvyZ5V9livuFQ9ffCu9ohEsGgXIipFCXC46In8DdrbKpD2I+R884uGLGeVwHfXtPXIHRE9TDOrb4IEQ3kw3vG52f6MwlCtJmOCoNwDyK6/tHFII9A+YQYYebTPTTwSOm2flyXwVrqkZj5qX9wWMD7NnzD7+q8iyeihY3SsI7TPYrI6wzy4IF13hyYV8nqzVzllnL1C527Uop8sgCOrpEeKw9vxv8CgMZSMB7fFQ0N6mrICVi36exDPoRudeWAmetlVIBn5KVmGKiozN2qjN0rF/unZaVn1aMf1GvZ/SqR3LasifelfNPJPfL+Rs0hY7xXDpgg2cUtgupbtJKNaoqV2wsx/rk+lfSYutBYCmZoVWs0Zqn0NAANT+Vbp5hF7ItswjnCitRYDTdO2RZ3oUWcOXOPQJZ6whnSi4qX8L6Ebzq3F941dngvjAMP8Z1S29T4VPfM+t8r34xIj0XAq7m1YhFr+0I2eOb0YlKPVRfH8K+ueO3a0isztJvzilpYfM5AyC95p8myUptiigHFZT/VNbsYBU/0wLeQpER+L2XqSzPvqtPEN4yozl55p49qdjIMaxAuBon+rAbqK1FFjU4U9djTeKJ0SbTkdtHVl6FHUwYWiKukg8xYGBihwrWujcrpFDMEMxi4pZP0Kl4GdEw1Z40M6sKtwlvWk/3WauH7xHlnve9+yVc8Y7xVIFz0rvaLEDxVzkFsRKdFj6EBrSWGRjuuD+e6ceAbJ2TjPyK6mUCpbYA4ZD2T1p6k/0C4/+0SGEuA9R6AIrZ4f7exIT3EW7k1P0Q84HXINyA0za7PGvW46YsQa1xt4IiYGJwNnVGgTDLCexXyPGBwdnLA8kix8Np4PK7X77sE0Dd0B0/GWq8KVIfI6rVHRv9uk4qpB0OnU5ykSpUeQwOjaQUmEcBNJM7Q9EC4yDZn6h6D66NGBFP5EUq/7RIoSuCVSk0Zi4i+0F4j8Ij3Sr0VtCwCA+dSF8Zcgho9xtlgtK2A4dNnFx2aodFxFx4aNSFAHEuOO7+Zyrd6Pe0SqJQAERSTe5WxCc8ZjcdHN4NTUB7ejryBDut5ElQnq9+ZCjUGm6nmwtgwGZz5Bml0YFB6xvg2Kk877wQT53saoHfoEqiSgNcfZntea7QA+OlyBp1k15aQ1ICQD72eglo9hsZcCDpoQ/1wRwWVGSYD82ggOhcX1QjrMFWMf2s/5INLmRI8X4cugVYJsGT5QTpUMoiUJrqW09XUIYBT0staPcajpe2lfUqfeQUy+GkHPNqdMulk6JQ3Htyx7Q7FvyURGvYC4rRAPJOoWocugXkJSG+Y8dG/q3nsgDFQetHAM3GlERk66d5jhZrDm1o95iV9tX77/mBLg/eN575t3DGn7J8q5U3/S11YPOcVyQuxpu4hwCmPuZP0QuvQJTArAdOna+kURjcH4Nd8EYIOsw6cAmt3To8x4rnTJXEb0B39/krx9H5cu993CawlARkYh895P4uSNoHqEH4+UornGYDy9vpVhNozzxp+Jcw6dGY3Pe0SOKAEvlfbrNXC8iXXF5W3fhGCEbZ4o7RZDB8atUDkx2deI8/cDa1WhB1vNQlIMVkDccwvXqLk2mPjjZDQ4IVuJr8IEV3WgOwhFA3YiMWpr0PdsHMYl6f3wsMg2Q3NnnTphpZKrD8fRAJSUM6sskVf8iB4GHA4cI5RslcwF26y/zBnwEIZAXWqvKH6QMhI+EpbWehrtKxYeuahJCBlZb3GxkjWMyzpl2jh0ap/mkD4bACONjRybQuX3Xk8a/FfhnVDy0mv5x1UAlJetur5Bqx6W/0QHVb/LtTuX7qKRkbf/gPY+FwBQmFyJAAAAABJRU5ErkJggg==\n", "text/latex": [ "$\\displaystyle \\frac{d}{d t} y{\\left(t \\right)} = y{\\left(t \\right)} \\log{\\left(\\frac{1}{y{\\left(t \\right)}} \\right)}$" ], "text/plain": [ "d ⎛ 1 ⎞\n", "──(y(t)) = y(t)⋅log⎜────⎟\n", "dt ⎝y(t)⎠" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAHQAAAAbCAYAAACtOKuoAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAFEElEQVRoBe2a63EUORCADeUAriACIAMOMoAMgAzAGUDxz/8oLgMgAgoysDPgIAPIAJcz8H2fLIlZjeaxszM2y01XafXqbrXUL2nsGxcXFwcr7PcJHB8f32UHr6iPbu73Vlbp4wk8ov5u+8bqofFI9rTCK18i+hHllHKyKnSCIuMhPoT0jHIeWbyhfk15w3wai1PLVqz3nXLPVQ6XXerP4s6hGdreUcxX/zR3R/+E/gPqV83xKW14hJzYQ/sJHD3ygPovqmxAq0J7Tq05xcE9of+Bcod2PsAGjor+0ehPbsJfPobRMaCR/SsidC9mUSiM7kYhegUYi9fL5BomlZtlP1Ge0q4pU6lUgl46CuCjZxmivczYPmfsPfW28A2Cx9C+oD7dWaEwMim7mTHWqeKfUDbCFbS/O7xFwB/I/blLUOY8WMsggKsCv1I0kEBDrcFsrVDoNrx5p2cLzAxDt6k7N8p8BvBS3Nea9gkMa0H2sUKzV433hKLySjB0f2POIp4GM9q7S2bN/mSFRkFfU1cvAc5TvH0ZrjLQ1zuPnM+D+9EI77wuUdmPxh2A9n0alo29X86GXw3kDDwNWzxvxlt7Z+BU/OwScrUqLwJdoNBu6KyCIJ30YxN/hcWVDumdPlOqgDJMOzlK0Q/hl9o91kBjfsv8mDRVo+8cm+yhcHyGQH1W9Rgc8855uXqkk35fvPQ5e3iEvBrpBjCmMg2d2yhHA9EzM0DfpfyMM6YxyUNZ3PAytAE3n622Ioz0zyh9RlEhu/oh9usN9A4r61VPqTXSn1GS987H9thKHvIygkmrYfthYmcICoWxDI3nfm34Sj8fcpz7QK0QCfS+1iUBHJVoGJWfwt5nzNvbF+ryZiu9fPJatK8FkE15fUIkJXkOfjw4TwLF9iwpYk5eSb5UJw8NlxsWSo/n5iHrRTnhR8IH1K38Cb1KOo18DFEqrAu8ZGgAgwAf1xqF22BmGGwaYWPqVxMc+cq/+YQwHKrg6oWP8V6Ap/TyvUXxbDXovmjVy2+bycO4+JdIpALKS0zNG7XoEq+5rjReDPpAer14EJBxFs8oF4Kv6/tc8NbdlFdlTH5GRF7yK6MSQ8vCIey9uKTN6I1lLHdz5ZiWl8MR7RKkGbJIc6iGcZ2QooxPLC83wm2K78dWSgmzv/mPHhoUQ60SPOAcbhkzdDg2enPQiJ8sn2YnDBlFJ+GME8HwkPnKPWnGPWywOmz0zDfmnaDgOO6GveElD07oZzRUXA2kOYAmGwFtPaDJVxTp5TMI0C6WQ1k8pZtBOfYBoalQvap8itTyp/sSryv/beRPlJGeOKVR6KHlevJuATwWyaEs5PqG2BYoN2UobbTornugqVA35yEHYDN6mqV2mCrnYUBs/2RFwUMv9Ftl7WD+Zq5UcpvbsiOtr1VRZsdneegvK36be1OhXtF9b7oR32PhL+DUOXTSTvCRhu/LGsjnHXzCB3jqrvzUZSw1nouMIZsfBUwH7jl9q7VfM+JFZJibaee/oMRN+pbUk1rAuAeQ324thJ4BaMOliToZTQ/2OrXNCdwUmYP1M5R/nwtA21Cph5XPlUuEy1+teqol68XSrzDzCQSFwtOLi2E0gX+vMxzVcl/AYc7njfmx63KUeG3UEV+6/DzaQFg7O51AyqF6zC0O2ce1YdAcWMud5WI+dcyl3mzHgk+QqZ49do3/LV5nDh17ItHjRv1bSTSYz9SjnitjZVjxfp3Af+gIw+kMuX0qAAAAAElFTkSuQmCC\n", "text/latex": [ "$\\displaystyle y{\\left(t \\right)} = e^{C_{1} e^{- t}}$" ], "text/plain": [ " -t\n", " C₁⋅ℯ \n", "y(t) = ℯ " ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAHUAAAAVCAYAAAB48KHmAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAF7UlEQVRoBe2Z65EUNxCAhy0C4JEBZMAjAiCDg4sAyMCu+3f/KMgALgIMGQARGMgAHIHPm8H5+4QkRlppZnZZfFCmq7SSWq1Wq1+SZi+cnZ0NNRwfH18Cd0s87bf1+K/+Mg2gu2uUz8uov1AtmSMN1DcoH1v8V/WCEP0G7l3Ev6/Hf/WXaSDqUcVvCzqCNuhCNKTO8pT2G4pBmOHCOFIZPGDkFeUy7XWm6jQis0cM36OM6V3oBUWPekz9e4fFXtGsoxJPKK77lv79vS6wkBnrqsfb9b7pK1fShZnw1D74j9QZ6GvUNfWLjOw0oHkTafNe60i9DcFnCMcGarKLC/8VB+/TzwXcKf3n1C74KdJ894o1TUc3WWirlLdPwVjfqDmiTsYL7Olr0OfUOrlFOTXmB9p3A1H8of+MpjTymgN5FBmhNuocg8GFKBrrMeUO7WeUwgnov46M3Mh5nMlGwHnBUxbWoWsQr84yoCcNr+7MjjXIwzlbw9ZGZQXPW1PHTYQq0ka1uoKaQs4taip5/qvuA/bcSptG4yfG6ujT6Q0UAyBD5CGvmj7T9BpbGZUF9BxD/SHtIjobCxgtfzTwPy2KPReKrzfCuGdpz4k13tTR1jKevB7U68z1L84RpPG4IQ9wBUvpNQ23agVqpaEW7SSO9TyXvfgU6WtyUmOQ+SpOx0zn/HXar8AXRwR9o8q13IOG9Lhxrnv3Eui8FnhhLHglIubki0zCxTqch4y3sp685NmK/IrN125tVAVffx0uWkmhiwyFkPJpCVownevAR5lUrCl/Z4CPyvNIuEc7RxNtjXSD4uVksE2lEfMLANw/9L2lvqb09ANJkHGRfiSGl2u5t+JS5VgEnU8Hm4K/GbwyJlilDgvI3FDvCWVqEZZE6RfK0S/8PTceUT6M0LNN6FWiCva2+C2gQTVKNmhkpkJ976lg4YhSp0kjJigeumYkOhHQAT12lkKSKThUY5K8tMsUaA91m+wzrOj42JW5xWdJL9SDNzBeK6W5IHRuMABtFeZ7VlzGh8EFP8xfLyDrksT1Vc6fNRFjKZsc1mM79NXRIllZ1+DReXpp2eXV9aS+mC/NZcohbZ9HBxpV5BPKe4oppsdksQfKGF7Zw+j7ftQbFzkEdPuGJMuUwlOkqmwdfawHx3Z6XrQ2Am8d/Aq15+UULHUSs4t7DEfESo4wV+npzDwR14CQdqBLCmqQZJRfU1IEZOQ5NpIzjQ1Vi5NodG6z1Ql7MC2bwax7GWzMR8efWmOAjw5/nTpHKG2dqKVXeU0GU+TnBdZvBsFG9UVJQ+TcTHsMeqoe5njvDBhgLJ2R/8MAMum0awQyMoo7Afh0EdF4gv2NG3EYmf/RMVrGCTNZy4jX4euLkTptOY2Rmpwt8Gj8+BXQ7wHuL0Bt1ITfqJmU8r8blkkhBH2VodKeOL7BYEcEvPRWb5/pE+COnIY7THwHPz/VjbOITuhXsXQBUnYj8yF12scp/dQG3QX5quQNYL7G1nF8mpnix3AXXCtQ0qfEMW2rXUTzYqPKiYW9PXoou2mv/WmjMtXYtQeC/jaApw6kt5oWJwE6I8Hz5ZaE9FWi54wOmZzCfpJbRdtPBh1sU0APxS0dnDLoEC3lSy+8pKSID4jRj/pyPbNdDWMnG48ZKOlYHOMn2/W/NHrtAYL3HteTzOYG4Wua0SG+C/+59ZeMI5vK13ghTVNfom8a1GHUz6Rhofdt6SuiZyiG54H5OsDUh47ABLoNm20VqfOi/NwUKMjIGKjzuUvbqLYY7RrXC85UtKpko2vrCGPOGMx68toaVtWMja8T1fhOXZTh7U4Bw9XbNsUb248GpnhlNSpboKHq87CgY653jd5ttqDtdeBhlMqjuLd06M0kOl2GOlL1UBVuCs7emql3bMDL82jv5+2O4nSnIafnt5eTI+pDap1cuEpRef6RsSStGs2erXPvUEiaoOPMRjqyKJNfAQvdFmeq7CHUSyQy1bykvzfjwu9/A1GPBsdUqt7QB/RmsNbnzEwbeWt0bbVho38BxUpgN3qgreoAAAAASUVORK5CYII=\n", "text/latex": [ "$\\displaystyle \\left\\{ C_{1} : \\log{\\left(2 \\right)}\\right\\}$" ], "text/plain": [ "{C₁: log(2)}" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAJEAAAAbCAYAAAB8822dAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAGjklEQVRoBe2a3XEUORCA1y4HwJkMTAbAZWBncJDB4Qxw+c1vLi4DIAIOMoAMOMjAzuBczsD3fbIkZjWanZllf7zcdJUsTau71Wq1Wi2t9+7u7mYTbNYCFxcXR4x4Rn262ZHXM9r+esROUnsscEz/VQ/NznRPTrThpSL6vGbIM8oT2n9sePjqcOpBMToOgjiHTLs3HWfZFhtrsAhXlCddA9L3iL6PlGPae110q8Aj/xVyrqm/JHlx/PP4nZzrT/C3DZq3fIfj+CAhp3p5C2BMDW106YKP0IRFotZB8mLUGKC5pbhAaz3yoi4n1C8KPd7E8QOa9lsa3yhNx3dO0p1NTlRYb5lPDHkN39Ak2XzoH8eB7xXlne0K3FRwq0a9QaAOUoJ6ZcenUzpxTynfJab+QtGJHq3EiRB0RNGQC2Eo3UIhW+rUWAxtiDc6hGgCrssBFmnpIrj7PUbyEbKIodkHnzmVkUwdHlMuwfk9oxbngqujUUP8V8opfSfUJXhc1pxfXHD0kqH4Vv/jn3YilHBSOlCvE0Gjs5nE/UV7ZwB9XRzD+QvaaSeas4x2IvjHRK05G8H7GYS7v3k0No8ZdTJX+SQjte83v6VvcQnAeQQH50u4VNNXzkunMm8Kc0901DroyX4DMbqJUG8Xj6mD0n0CoEuTdxfuErxH2e8akeJGcLe7oBsDxnzKYEaOHL1o6wQurhtZqDnG4X1X66+0Ny1sgYjjus7Pii4/Hf9oaSdCuLvznLqaUNpP8Raishn4NgoZXuXfFTCPuUFnnd/F9Agpd+u65/KcAVy0EsSlhFfnDrZFPxfevKXrhND+NXmg7wFe106Zz2jXaIMT/sxxpvBaUnavAbuGRpe3yxcmnIgfeK3BPUa6FmQT6pujqEcJ2thjRbDfTf2SYoSq5UF0BXAuNXmhMzqQr+pBRvyeUTdt4NjX+4FjuT8vEbhoNzq4E7ktxUc++TsnUfJs+dsjxAiUAd3dBBsDxjMf8Tg1wgSg7SJa0joYkTwB3lHysReIiz/0K0/eFkS5bnTzK29kzl3nLI8/+a+WikQIdSJNj+SzBUaiRbmS/O6YZICWgAeE8B3FSKTR3BQ6/yVlLcA4yg9OStuFDDcoao+VpIdj6zTNo8bE+hs0VAHU9e/EH3HNyk1eu1mbrKuDdYaKHAPFaXixplMGz3uVUom8sLHvPXV+kKKtl/ogNpcP8a3jOGHl2XY3OJGv9M3dyPjWSE4gy+V7K4AO6ntO+TcqoB0M5eq+E4CuOrjr4g0y6B1x2llnmVsrJwXONfK5odVn/yKAR5v5lnRyEAlDggzCCONNJDsRbaNFDqG0hecUFZ4D+HUakznpvUksOpN9y3ASvYAcxxpE2xBm6O91UGiUm4xviJ+BM3yfU0YbV/4tQYiS6J4dn7bO4y3STdEC+lwrna4WjVr0BULnDBHyAAEaLCVmLnp57okrz1e9sKQDlUGesCAZ027I78R7AR2Dsr2EIwk0Hiwa2RtNU18da6NX+JGqt8jRX4cwH3JD6EhGVZ3HG3LnZqDPuXtEeuPMDghfJ0DnqWUUCinNAR96azKgUeey4NagJe4Q3KIB5VmUDzmECuiM2wQNLmj81/fN8Ar8me9y48Tuh1uhszbvs3trAvB5dI9ZC/OsvP5GovBB7cIrKB9l4IxS4gYbFB7p0w6n2Ql9jtjJuMKO4OzoPJevrVD+zojCBsEPhihc0h40mMwfzCOawjSyCXSKVInco0hnqYE8M3iy49F2pzflSiK/cnoBXiNGkNtL/INgUE4EeTrKf3BOrVEWaDqR0aO8ttfyIQeQriufmcuHcACTbOlLRzQSleOBagMy1pITMZLj+yNmC9SbMvpoaAn6HyCaTqRBXdgAGNCdb6ktoA7xeyBs/8nOgQyjjZl/bTH8LaZ0rLa09WLyDSMNE3UWb5lggAWaTmQG73uQxkuZvSLyseRHhA/UPmzVQDk+kIUfWam78o0uB63JXAsO3XzZ9ah1zukfwPyubZy16PArCO3899ho2GPq2q+3M/AaPf9rxBhjwBsSb+rq+8UYWRPt9i2wrwospu8E+YmbtseQkaS82kuewN277I41Wsk/wS9ggeBEzMPk1yMqQXi1xplquUygoc+nAPOdrgQ7yZqrI718+SlhjmD62DkLpJzIyHDIwvrg5hFjTlPLhcoJ+ixgbuSNbCh4XV82gg0dY6LboAU6c6KhOsTIMuhfXqOTfqIedLUfqsNEt10L/AcKp4OBi6apdQAAAABJRU5ErkJggg==\n", "text/latex": [ "$\\displaystyle y{\\left(t \\right)} = e^{e^{- t} \\log{\\left(2 \\right)}}$" ], "text/plain": [ " -t \n", " ℯ ⋅log(2)\n", "y(t) = ℯ " ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import sympy as sym\n", "sym.init_printing()\n", "\n", "t = sym.Symbol('t')\n", "y = sym.Function('y')\n", "edo= sym.Eq( sym.diff(y(t),t) , y(t)*sym.log(1/y(t)) )\n", "display(edo)\n", "solgen = sym.dsolve(edo)\n", "display(solgen)\n", "\n", "consts = sym.solve( sym.Eq( y0, solgen.rhs.subs(t,t0)) , dict=True)[0]\n", "display(consts)\n", "solpar=solgen.subs(consts).simplify()\n", "display(solpar)\n", "\n", "sol_exacte = sym.lambdify(t,solpar.rhs,'numpy')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- On calcule la solution approchée pour différentes valeurs de $N$:" ] }, { "cell_type": "code", "execution_count": 85, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Multipas ordre=4.04\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfkAAAFACAYAAAC/cwVZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3de7yVY/7/8dcnKpPMz1DKqcMMOhBlh8LQzimGmIwptoZhbCQTg6HZvuNrqGESohz65lT2yKFJETqwK4dCWwedndJByVm7RIfP749rmbbau9autda91r3fz8djPfa+D+71+Vjdfbqudd3XZe6OiIiIxE+NqAMQERGR9FCRFxERiSkVeRERkZhSkRcREYkpFXkREZGYUpEXERGJqZ2jDiDV6tWr502aNEnLtVevXs2uu+6almtnWlxyiUseEJ9clEf2iUsucckDUp9LaWnp5+5ef/P9sSvyTZo0Ydq0aWm59sSJE+nQoUNarp1pccklLnlAfHJRHtknLrnEJQ9IfS5m9nFF+9VdLyIiElMq8iIiIjGlIi8iIhJTKvIiIiIxpSIvIiISUyryIiIiMaUiLyIiElMq8iIiIplQXAxNmkCNGrTr1i1sp1nsJsMRERHJOsXFUFgIa9YAsMunn4ZtgIKCtL1t5C15MzvHzOaY2UYza1tuf4GZzSj32mhmraOMVUREZLsUFf23wP/XmjVhfxpFXuSB2UAXYHL5ne5e7O6t3b010B1Y5O4zoghQRERkhyxeXLX9KRJ5d727zwMws62ddi7wREYCEhERSaWFC6F2bVi7dstjjRql9a2zoSWfjK6oyIuISC757jv4+9+hVSswg5o1f3q8Th3o0yetIZi7p/UNAMxsAtCwgkNF7j4qcc5E4Fp3n7bZf3sUMMTdW23l+oVAIUCDBg3yhg8fnqrQf6KsrIy6deum5dqZFpdc4pIHxCcX5ZF94pJLLuWxx5tvcuA99/CzTz5hxUkn8eFll7H7O+/wyyFDqL1yJd/Vq8eiwkJWnnhiSt4vPz+/1N3bbnHA3bPiBUwE2law/y7gb8leJy8vz9OlpKQkbdfOtLjkEpc83OOTi/LIPnHJJSfyWLLE/eyz3cG9eXP3V16p8LRU5wJM8wpqYlZ315tZDeAcID1NcxERkVRYtw7694fmzeGFF6BvX5g5E/LzIw0r8iJvZr81s6VAe2CMmY0td/g4YKm7fxhNdCIiItvw+uuQlwfXXhuK+pw50Ls31KoVdWTRF3l3H+nu+7l7bXdv4O6nlDs20d3bRRmfiIhIhT7/HC6+GI49Fr75Bp59Fp57Dpo2jTqy/4q8yIuIiOSUjRthyBBo1gyGDoXrr4e5c+HMM6OObAuRPycvIiKSM2bOhMsvhylT4Ljj4L774OCDo46qUmrJi4iIbMuqVfCXv4Tv3t9/Hx57DCZOzOoCD2rJi4iIVM4dnnkGrroKli+HSy8NE9jssUfUkSVFLXkREZGKvP8+nHoq/P730KBB6KK///6cKfCgIi8iIvJTa9fCzTfDIYeEwn7PPfDWW3DUUVFHVmXqrhcREfnRuHFwxRWhFX/uuWGCm733jjqq7aaWvIiIyLJl0LUrnHIK1KgB48fDv/+d0wUeVORFRKQ6W78e7r4bWrSA0aPhlltg1ixI0cIxUVN3vYiIVE9TpoRn3mfODAPsBg6EX/4y6qhSSi15ERGpXr74AgoL4eijw+8jRsCYMbEr8KAiLyIi1cXGjfDII2GluIcfhmuugXnzoEsXMIs6urRQd72IiMTf7Nmha/611+CYY8Lz7q1aRR1V2qklLyIi8VVWBtddB61bh1b7ww/D5MnVosCDWvIiIhJH7jByJPTqBUuXwp/+BLfdBnvuGXVkGaWWvIiIxMuHH8Lpp8PZZ4ei/sYb8H//V+0KPKjIi4hIXHz/Pdx6a1gZbvJkuPNOmDYN2rePOrLIqLteRERy38svQ48esHAhnHMO3HUX7Ltv1FFFTi15ERHJXcuXw3nnhRnqNmyAl16Cp55SgU9QkRcRkdyzYQPce2945n3ECLjppvCY3CmnRB1ZVlF3vYiI5Ja33grPvL/zDpx8cpiO9sADo44qK6klLyIiOWHnVatCcW/XLnTTP/lk6J5Xga+UWvIiIpLd3GHYMI7s1Qu+/TY8+37zzfDzn0cdWdZTkRcRkew1d24YNT9pEt+1bEmtkpIwe50kRUVeRESyz+rVYW33/v1ht91g8GCm/+pXdFCBrxJ9Jy8iItll1Cho2RJuvx26d4cFC+CSS6CGSlZV6f+YiIhkh0WLoHNnOOus8H37q6+GBWXq1486spylIi8iItH64YeweEzLlvDKK9CvX3g87thjo44s5+k7eRERic7EiWFg3bx50KUL3H037L9/1FHFRuQteTM7x8zmmNlGM2tbbn9NM3vMzN41s3lm1jvKOEVEJIU+/TR8356fD2vXwpgxYeY6FfiUirzIA7OBLsDkzfafA9R291ZAHnCpmTXJbGgiIpJSGzbAffdBs2ZhMpsbb4Q5c+C006KOLJYi765393kAZrbFIWBXM9sZ+BnwA/BtZqMTEZGUKS0NM9a9/TZ07Lip2EvaZENLvjLPAKuB5cBi4A53/zLakEREpMq+/hp69oQjjoAlS+Df/4YJE1TgM8DcPf1vYjYBaFjBoSJ3H5U4ZyJwrbtPS2wfA/QALgR+AbwKnOruH1Zw/UKgEKBBgwZ5w4cPT0MWUFZWRt26ddNy7UyLSy5xyQPik4vyyD6R5eLOXi+/zAH33UfNb75h2Zln8tFFF7FhO2PRZ1K5/Pz8Undvu8UBd8+KFzARaFtuexDQvdz2w8Dvt3WdvLw8T5eSkpK0XTvT4pJLXPJwj08uyiP7RJLLvHnuHTu6g/sRR7hPm7bDl9RnUjlgmldQE7O5u34x0NGCXYF2wPyIYxIRka1ZswaKiuDQQ8Oz7vffD1OmQF5e1JFVS5EXeTP7rZktBdoDY8xsbOLQIKAuYfT928Aj7j4rojBFRGRbxoyBgw+Gvn2hWzeYPx8uuwx22inqyKqtbBhdPxIYWcH+MsJjdCIiks0WLw7Lvz77LLRoESa4Of74qKMSsqAlLyIiOWrdujAFbYsWMHZsmJp2xgwV+CwSeUteRERy0Kuvhmfe58yBM8+EAQOgceOoo5LNqCUvIiLJ++wzuPBCOO44KCuD0aNDN70KfFZSkRcRkW3buBEefDBMYFNcDDfcEFrxZ5wRdWSyFequFxGRrZs+PXTNv/kmdOgAgwaFZWEl66klLyIiFfv2W7jqKmjbFj76CIYNC+u9q8DnDLXkRUTkp9zhqafg6qthxYrQir/1VvjFL6KOTKpIRV5ERDZZuDAsJjN+fJilbtSosLCM5CR114uICHz3Hfz979CqVfjufeDA8FMFPqepJS8iUt299FJovX/wAZx3HvTvDw0rWjhUco1a8iIi1dXSpfC738Gpp8LOO8PLL4fH41TgY0NFXkSkulm/Hu68M0xHO2ZMGFQ3cyZ07Bh1ZJJi6q4XEalOXn8devSAWbPgN7+Be++Fpk2jjkrSRC15EZHq4PPP4U9/gmOPha++gpEj4bnnVOBjTi15EZE427iRhmPGhO/ev/kGrrsujKKvWzfqyCQDVORFROJq1iy4/HKav/EG/PrXcN99cMghUUclGaTuehGRuFm1Cq65Bg4/HBYuZP7118OkSSrw1ZCKvIhIXLjDM8+EUfN33gkXXwwLFrCiUycwizo6iYCKvIhIHHzwAZx2GpxzDtSvD1OmhKVh99gj6sgkQiryIiK5bO1a+Mc/4OCDw+Nxd98Nb78N7dpFHZlkAQ28ExHJVePHwxVXwHvvQdeuoYt+n32ijkqyiFryIiK55pNPoFs3OPnk8D382LEwfLgKvGxBRV5EJFesXw8DBkDz5vDss3DzzfDuu6HYi1RA3fUiIrlg6lS4/HKYMQNOOSUsBXvAAVFHJVlOLXkRkWz25Zdw6aVw9NHw2Wfw9NPw4osq8JIUFXkRkWzkDo8+Cs2awUMPwdVXw7x5YXpaPfMuSVJ3vYhItpk9O3TNv/YatG8PDzwAhx4adVSSg9SSFxHJFmVl8Ne/Qps2MHcuDBkSCr0KvGwnteRFRKLmHkbL9+oFS5aE6Whvuw3q1Ys6MslxkbfkzewcM5tjZhvNrG25/bXM7BEze9fMZppZhwjDFBFJj48+gjPOgC5dYPfdQ8t9yBAVeEmJyIs8MBvoAkzebP8lAO7eCjgJ6G9m2RCviMj2KS6GJk2gRg1o3DjMM9+yZVghrn9/KC2FY46JOkqJkci76919HoBtOVq0JfBy4pyVZvY10BZ4K6MBioikQnExFBbCmjVhe/Hi8DrySBgxAvbbL9r4JJayuWU8EzjTzHY2s6ZAHrB/xDGJiGyfoqJNBb68Tz9VgZe0MXdP/5uYTQAaVnCoyN1HJc6ZCFzr7tMS2zsD/YB84GOgJvDgj+dvdv1CoBCgQYMGecOHD09HGpSVlVG3bt20XDvT4pJLXPKA+OSiPLZkGzZw3IknUtHT7W7GpFdeScn7VEafSfZJdS75+fml7t52iwPunhUvYCLQdivH3wBabus6eXl5ni4lJSVpu3amxSWXuOThHp9clMdmpk51b9PGPYyh3/LVuHFq3mcr9Jlkn1TnAkzzCmpi1nbXm1kdM9s18ftJwHp3nxtxWCIiyfniizAdbfv2oUv+yiuhTp2fnlOnDvTpE018Ui1EXuTN7LdmthRoD4wxs7GJQ3sB75jZPOB6oHtUMYqIJG3jRnjkkbBS3I/T0c6fD/fcA4MHh1H1ZuHn4MFQUBB1xBJj2TC6fiQwsoL9i4BmGQ9IRGR7zZoFPXrA66+HR+Huu++ns9UVFKioS0ZF3pIXEcl5q1bBX/4Chx8OCxaElvzkyZqOViIXeUteRCRnuYelX6++GpYvD8/B9+0Le+wRdWQigFryIiLbZ+FCOOUU6NoVGjSAKVPCanEq8JJFVORFRKriu+/g73+HVq3gzTfh3nvh7bfhqKOijkxkC+quFxFJ1gsvQM+eYVGZggK44w5oWNE8XyLZQS15EZFtWbwYfvtb+M1voHZteOUVePxxFXjJeiryIiKV+eEHuP12aNECxo0La7zPnAn5+VFHJpIUddeLiFRg9xkzwjPv8+bBWWfB3XeHCWxEcoiKvIhIeStWwHXX0frxx8Pa7889B6efHnVUIttF3fUiIgAbNsCgQWE62qeeYlH37jBnjgq85DQVeRGRt96CI48MI+ePOAJmzWLRRRdtuaCMSI5RkReR6uvLL+Gyy6BduzBj3fDhYYBdMy2bIfGgIi8i1Y87PPpoKOZDhkCvXmGluK5dwwpxIjGhIi8i1cu778Jxx8Ef/wgHHgilpXDXXfDzn0cdmUjKqciLSPWwahVcey20aRMei3voIXjtNTjssKgjE0kbPUInIvHmDs88E1aKW7YMLrkE/vlP2HPPqCMTSTu15EUkvt57Dzp1gt//HurXDyvFDR6sAi/Vhoq8iMTPd9/BTTfBIYfA1KkwYEBYKa5du6gjE8koddeLSLy8+GJ43v3DD+G888JKcXvvHXVUIpFQS15E4mHJEjj7bDjtNKhZE15+GYqLVeClWlORF5Hctm4d9OsXVop78UXo2xdmzYKOHaOOTCRy6q4Xkdw1eTJcfjnMnQudO4fv3ps0iToqkayhlryI5J5PP4ULLoDjj4fVq2HUqPBSgRf5CRV5EckdGzbAffeF6WifeAKKija14kVkC+quF5Hc8PbboWu+tBROOAEGDgzLwopIpdSSF5Hs9tVX0KMHHHVUmLHuiSdg/HgVeJEkqMiLSHZyh6FDQ9f8gw/Cn/8cVorr1k0rxYkkSd31IpJ9Zs8OrfdXXw2z1I0bB61bRx2VSM5RS15EskdZGVx3XVgpbs6csNb766+rwItspyoXeTPb1cx2SlUAZtbPzOab2SwzG2lmu5c71tvM3jezBWZ2SqreU0SyjDuMGBEmtLnjjvB43IIFcPHFUENtEZHttc27x8xqmNl5ZjbGzFYC84HlZjYnUaAP3MEYxgOHuPuhwEKgd+J9WwLdgIOBTsB9qfzHhYhkifffD1PR/u53YXW4118PLfh69aKOTCTnJfNP5BLgV4Ti29Dd93f3vYBfA1OB28zs/O0NwN3Hufv6xOZUYL/E72cCw939e3f/CHgfOHJ730dEsszatXDzzWGluNdfh7vvhmnT4Oijo45MJDaSGXh3oruv23ynu38JjABGmFnNFMVzEfBk4vd9CUX/R0sT+0Qk1730Ulgp7oMPwmj5/v1hn32ijkokdszdkzsxdMv3Br5z9yuq9CZmE4CGFRwqcvdRiXOKgLZAF3d3MxsETHH3xxPHHwJecPcRFVy/ECgEaNCgQd7w4cOrEl7SysrKqFu3blqunWlxySUueUB8ctlaHrVXruSAQYOoP3kya/bfn/d69eKrvLwMR5icuHweEJ9c4pIHpD6X/Pz8Undvu8UBd0/qRWhVnwrMSmwfAgxN9r/fxrUvAKYAdcrt6w30Lrc9Fmi/rWvl5eV5upSUlKTt2pkWl1zikod7fHKpMI8ffnDv1899113dd9nF/dZb3deuzXhsVRGXz8M9PrnEJQ/31OcCTPMKamJVhq3WcPcXgQ2JfxzMThT6HWJmnYDrgc7uvqbcodFANzOrbWZNgQOBt3b0/UQkw159NTwSd911kJ8f5povKoLataOOTCT2qlLkP0kUWwcwMwN+loIYBgK7AePNbIaZPQDg7nOAp4C5wEvAFe6+IQXvJyKZsHIlXHghHHccrFoVVol77jlo2jTqyESqjarMeHcVMARoaGZ/JDzWNntHA3D3A7ZyrA/QZ0ffQ0QyaMMGeOAB6N07LAPbu3doue+6a9SRiVQ7SRd5d1+U6Fo/CzgMmAQ8nK7ARCQHlZZyeM+eYY75/HwYNChMcCMikdhmkTczS3ypj4fn2Z9JvCo8R0Sqoa+/Dq31++9nl913h+JiOPdcLSQjErGkJsMxsyvNrFH5nWZWy8w6mtljhNHxIlLduMOwYWGluAcegJ49eXPoUDjvPBV4kSyQTJHvRBhR/4SZfWJmc83sI+A94FzgLnd/NI0xikg2mjMndMn/4Q9hMN20aXDPPWyIyXPMInGwze56d18L3EeYO74mUI8wIc7X6Q5ORLJQWRnccgvceSfsthsMHqyFZESyVJXWk/cwve3yNMUiItnMHUaOhKuugiVL4KKL4LbboH79qCMTkUokXeTNrDZwNtCk/H/n7v9IfVgiErni4jCYbvFi2HvvsCrcrFnQqhU88QQcc0zUEYrINlSlJT8K+AYoBb5PTzgikhWKi6GwENYkJqH85JPwKiiARx+FnavUCSgiEanKnbqfu3dKWyQikj2KijYV+PJee00FXiSHVGWkzBtm1iptkYhIdli6FD7+uOJjixdnNhYR2SFVKfLHAu+Y2QIzm2Vm75rZrHQFJiIZtm4d3HEHNG9e+TmNGlV+TESyTlX63ToBRmKBGhGJkVdfhcsvD8++n346nHDCll32depAHy0lIZJLkpnWdhUVF/YfC/7PUx2UiGTIypVhCdihQ6Fx47BSXOfO4Vj9+ptG1zdqFAp8QUG08YpIlSQzGc5umQhERDJowwZ48MFQxFevhr/9Lfxep86mcwoKVNRFcpyGyYpUN2+/HbrmS0uhY8ewUtzWvocXkZyleShFqouvvgrF/aijwjPvTzwBEyaowIvEmIq8SNxt3BgmsGnWLMwz36tXWO+9WzetFCcSc+quF4mzd9+FHj3CJDbt28O4cdC6ddRRiUiGqCUvEkerVsE110CbNjBvHjz0UCj0KvAi1Ypa8iJx4g5PPw1XXw3Ll8Mll0DfvrDnnlFHJiIRUEteJC4WLoRTToGuXaFBA5gyJTwmpwIvUm2pyIvkujVr4H/+JywB++abcO+94TG5o46KOjIRiZi660Vy2fPPw5VXwqJFcP750K8fNGwYdVQikiXUkhfJRR9/DGedBWecEWapKymBYcNU4EXkJ1TkRXLJ99+HgXQtWsD48XD77TB9OnToEHVkIpKF1F0vkitefhmuuAIWLIAuXeCuu7T0q4hslVryItnuk0/g3HPhxBNh/Xp44QUYMUIFXkS2SUVeJFutXw8DBoS55UeOhP/9X5g9G049NerIRCRHqLteJBu98UaYjnbmTOjUKTwWd8ABUUclIjkm8pa8mfUzs/lmNsvMRprZ7on9e5pZiZmVmdnAqOMUyYjPP4eLL4ZjjoEvvgjd8i+8oAIvItsl8iIPjAcOcfdDgYVA78T+tcD/ANdGFZhIxmzcyN7PPRdWihs6FP761zDnfJcuWilORLZb5EXe3ce5+/rE5lRgv8T+1e7+GqHYi8TXO+/A0UfT7M47w6x1M2aER+Pq1o06MhHJcZEX+c1cBLwYdRAiGfH112G2uiOOgEWLmPe3v4VJbQ4+OOrIRCQmzN3T/yZmE4CKpuIqcvdRiXOKgLZAFy8XlJldCLR1955buX4hUAjQoEGDvOHDh6cw+k3KysqoG5PWVVxyyck83GkwYQK/uv9+an7zDcvOPJNFF13E15B7uVQgJz+TCsQlD4hPLnHJA1KfS35+fqm7t93igLtH/gIuAKYAdSo4diEwMNlr5eXlebqUlJSk7dqZFpdcci6POXPcjz/eHdyPPNK9tPS/h3Iul0ooj+wTl1zikod76nMBpnkFNTHy7noz6wRcD3R29zVRxyOSFmVlcP31cNhhMGtWWAJ2yhQ4/PCoIxORGMuG5+QHArWB8RZGEU9198sAzGwR8HOglpmdBZzs7nOjClSkytzDRDa9esHSpXDRRXDbbVC/ftSRiUg1EHmRd/dKHwB29yYZDEUktT74IAyse/FFOPRQePJJOProqKMSkWok8u56kdhZuxZuvjmMkn/ttbCQTGmpCryIZFzkLXmRWHnxxdB6/+AD6NYN+veHffaJOioRqabUkhdJhSVL4Oyz4bTTYOedYcIEeOIJFXgRiZSKvMiOWLcO+vWDFi1CK75v37CozAknRB2ZiIi660W226RJYaW4uXPhzDPh7ruhSZOooxIR+S+15EWqasUK6N4dOnSANWvguefg2WdV4EUk66jIi1SmuDgU7ho1ws9hw2DgQGjeHJ56Cm68EebMgdNPjzpSEZEKqbtepCLFxVBYGFrqAB9/DBdcECa3OfFEGDQIDjoo2hhFRLZBRV6kIkVFmwr8j9yhXj0YN05rvItITlB3vUhFFi+ueP8XX6jAi0jOUJEX2dzMmVCrVsXHGjXKbCwiIjtARV7kR99+C1ddFVaGq1Vry0Jfpw706RNNbCIi20FFXsQ9zE7XrBnccw9cemkYaPfww9C4ceieb9wYBg+GgoKooxURSZoG3kn1Nm8eXHEFlJRA27YwejQccUQ4VlCgoi4iOU0teameVq+GG26Aww6D6dPh/vth6tRNBV5EJAbUkpfqxR1GjYJevcII+gsvhNtvh732ijoyEZGUU5GX6uPDD8MysC+8AK1awauvwrHHRh2ViEjaqLte4m/tWvjHP6BlS5g8Ge68E0pLVeBFJPbUkpd4e+kl6NkTPvgAunaF/v1h332jjkpEJCPUkpd4WrIEzj4bTj0Vdt4Zxo+H4cNV4EWkWlGRl3j54Qf417+gRQt48cUwec3MmWFRGRGRakbd9RIfkyZBjx4wdy507gwDBmiNdxGp1tSSl9y3YgV07w4dOoSV40aPDo/JqcCLSDWnlrzkrvXr2fc//4HHHgsj6G+8EXr3DnPMi4iIirzkqKlT4fLLOXDGDDj5ZLj3XjjooKijEhHJKuqul9zyxRdwySXQvj189hlzbropPCanAi8isgUVeckNGzfCkCGhmD/yCFxzDcybx2cdOoRV4kREZAvqrpfsN316GDU/dSr8+tcwaFCYllZERLZKLXnJXt98A3/+c1gC9sMPwwC7SZNU4EVEkhR5kTezfmY238xmmdlIM9s9sf8kMys1s3cTPztGHatkiDsUF0OzZjBwIFx2GcyfD3/4g7rmRUSqIPIiD4wHDnH3Q4GFQO/E/s+BM9y9FXABMCyi+CST5s6Fjh3h/POhUSN4663QPf+LX0QdmYhIzom8yLv7OHdfn9icCuyX2D/d3T9J7J8D7GJmtaOIUTKgrAyuvx4OOyxMQ/vAAzBlSuiqFxGR7ZJtA+8uAp6sYP/ZwHR3/z7D8Ui6ucPIkdCrFyxdCn/8I9x+O9SvH3VkIiI5z9w9/W9iNgFoWMGhIncflTinCGgLdPFyQZnZwcBo4GR3/6CS6xcChQANGjTIGz58eIozCMrKyqhbt25arp1p2ZDLz5Yt44B77mHPt96i7Je/ZOFVV/FtFQfVZUMeqRKXXJRH9olLLnHJA1KfS35+fqm7b9n16e6RvwjfuU8B6my2fz/C9/THJHutvLw8T5eSkpK0XTvTIs1lzRr3m25yr13bfbfd3O+6y33duu26lD6T7KM8sk9ccolLHu6pzwWY5hXUxMi7682sE3A9cLy7rym3f3dgDNDb3V+PKj5JsRdegCuvDI/EdesG/fvDPvtEHZWISCxFPvAOGAjsBow3sxlm9kBif0/gAOB/EvtnmNlekUUpO2bxYujSBX7zG6hZEyZMgCeeUIEXEUmjyFvy7n5AJftvBW7NcDiSCsXFUFQUCvv++0O7dvD882GQXd++YUraWrWijlJEJPYiL/ISM8XFUFgY1nWHUOgXL4a8PBgxAho3jjY+EZFqJBu66yVOioo2FfjyPv9cBV5EJMNU5CV11q+Hjz+u+NjixZmNRUREVOQlRd54Y+uz0zVqlLlYREQEUJGXHfXZZ3DxxXDMMfDFF2HVuDp1fnpOnTrQp0808YmIVGMq8rJ9Nm6EBx8MK8UNHQrXXQfz5sGAATB4cPj+3Sz8HDwYCgqijlhEpNrR6HqputJS6NEjrBB3/PFhlbiDD950vKBARV1EJAuoJS/J++oruOIKOOKIMMDu8cehpOSnBV5ERLKGirxsmzs89ljomn/gAejZE+bPD611s6ijExGRSqi7Xrbu3XdD1/xrr4WZ68aOhTZtoo5KRESSoJa8VGzVqjD9bJs2YUDdkCHw+usq8CIiOUQtefkpd3j6abj6avjkE7jkEq5zXlkAAAxXSURBVPjnP2HPPaOOTEREqkgtedlkwQI4+WTo2hUaNICpU8PjbyrwIiI5SUVewlzzRUXQqhW8/TYMHBh+HnVU1JGJiMgOUHd9dTd6dJil7uOPoXt36NcvtOJFRCTnqchXU7ssXw5nnBHWeT/4YJg0CY47LuqwREQkhVTkq5vvv4d//Ysjbr0VatYMLfdevcLvIiISKyry1cnYsWEim/ff54vjj2evxx+H/faLOioREUkTFfnqYOnS8EjcM8/AgQfC2LHMrVWLvVTgRURiTaPr42zdutAd37x5+O79llvCDHYnnxx1ZCIikgFqycfVpElhMZk5c8IAuwEDoGnTqKMSEZEMUks+blasCI/CdegAq1eHR+RGj1aBFxGphlTk42LDhjCJTbNm8NRTcOONm1rxIiJSLam7Pg6mTg0rxU2fDiedFIr9QQdFHZWIiERMLflc9sUXYQGZ9u1h5crQgh87VgVeREQAFfncUlwMTZpAjRph0ZjGjeGRR+Daa8NysOecA2ZRRykiIllC3fW5orgYCgvDYjIAX34Zin3fvnD99dHGJiIiWUkt+Vxxww2bCvyPNm6E+++PJh4REcl6aslnO/fQil+6tOLjixdnNh4REckZkbfkzayfmc03s1lmNtLMdk/sP9LMZiReM83st1HHmnFz5kB+fnjuvVatis9p1CizMYmISM6IvMgD44FD3P1QYCHQO7F/NtDW3VsDnYAHzax69DyUlcF110Hr1mEa2sGD4aGHoE6dn55Xpw706RNNjCIikvUiL5ruPq7c5lTgd4n95b+A3gXwTMYVCfewiMzVV8OyZfCnP8E//wn16oXjZlBUFLroGzUKBb6gINqYRUQka0Ve5DdzEfDkjxtmdhTwMNAY6O7u66MKLO0WLgzLwI4fD23ahGLfrt1PzykoUFEXEZGkmXv6G8hmNgFoWMGhIncflTinCGgLdPHNgjKzFsBjwHHuvraC6xcChQANGjTIGz58eIozCMrKyqhbt25Kr1lj7VoaFxez/5NPsrFWLT66+GKWde4MO+2U0vfZXDpyiUJc8oD45KI8sk9ccolLHpD6XPLz80vdve0WB9w98hdwATAFqLOVc0oI39Fv9Vp5eXmeLiUlJam94OjR7k2auIP7+ee7L1+e2utvRcpziUhc8nCPTy7KI/vEJZe45OGe+lyAaV5BTYx84J2ZdQKuBzp7ue/hzazpjwPtzKwx0AxYFEmQqfbRR9C5c3jVqQMTJ8KwYdCwos4OERGR7ZMN38kPBGoD4y1MyTrV3S8DjgVuMLN1wEagh7t/Hl2YKfD999CvXxgwt9NO4fdevaBmzagjExGRGIq8yLv7AZXsHwYMy3A46TNuXBhY99578LvfwV13wX77RR2ViIjEWOTd9bG3dGlYOOaUU8L22LHw9NMq8CIiknYq8umybl3ojm/eHJ5/Hm65JUxsc/LJUUcmIiLVROTd9bE0aRJccUWYlvaMM2DAAGjaNOqoRESkmlFLPpVWrAjzzHfoAKtXw+jR4aUCLyIiEVCRT4X16+Hee6FZM3jqKbjxxk2teBERkYiou35HTZ0Kl18OM2bASSfBwIFw0EFRRyUiIqKW/Hb7/POwgEz79vDZZ6EFP3asCryIiGQNFfmq2rgxLP3arBk8+ihcey3MmxcekwuT+YiIiGQFFfmtKS6GJk2gRg1o0oRGQ4eGlvull8Ihh4Qu+n79YLfdoo5URERkC/pOvjLFxVBYCGsS0+l//DFNH3kkFPShQ+H889VyFxGRrKYiX5miok0FPsEA/t//C4/JiYiIZDl111dm8eKK9y9bltk4REREtpOKfGUaNarafhERkSyjIl+ZPn3CWu/lbKhdO+wXERHJASrylSkoCI/KNW4cBtg1bsyCa68N+0VERHKAivzWFBTAokXh2fhFi1h54olRRyQiIpI0FXkREZGYUpEXERGJKRV5ERGRmFKRFxERiSkVeRERkZhSkRcREYkpFXkREZGYUpEXERGJKXP3qGNIKTP7DPg4TZevB3yepmtnWlxyiUseEJ9clEf2iUsucckDUp9LY3evv/nO2BX5dDKzae7eNuo4UiEuucQlD4hPLsoj+8Qll7jkAZnLRd31IiIiMaUiLyIiElMq8lUzOOoAUiguucQlD4hPLsoj+8Qll7jkARnKRd/Ji4iIxJRa8iIiIjGlIr8VZva/ZrbMzGYkXqdVcl4nM1tgZu+b2Q2ZjrMqzOxaM3Mzq1fJ8avNbI6ZzTazJ8xsl0zHmIwk8tjdzJ4xs/lmNs/M2mc6xmRtK5fEOTuZ2XQzez6TsVXF1vIws/3NrCTxWcwxs15RxJiMJP5sZf39bma3mNmsxN9b48xsn0rOy+r7vQp5ZP39nmwuiXNTdr+ryG/bXe7eOvF6YfODZrYTMAg4FWgJnGtmLTMdZDLMbH/gJGBxJcf3Bf4MtHX3Q4CdgG6ZizA528ojYQDwkrs3Bw4D5mUitqpKMheAXmRpDpBUHuuBa9y9BdAOuCIb75Mk7pFcud/7ufuh7t4aeB74++Yn5Mj9vs08EnLhfk82F0jh/a4iv+OOBN539w/d/QdgOHBmxDFV5i7gr8DWBmLsDPzMzHYG6gCfZCKwKtpqHmb2c+A44CEAd//B3b/OXHhVss3PxMz2A34DDMlUUNthq3m4+3J3fyfx+yrCX2D7Zi68pG3r88iJ+93dvy23uSuV55PV93syeeTK/Z7sZ5Lq+11Fftt6JrpYHjazX1RwfF9gSbntpWThX15m1hlY5u4zKzvH3ZcBdxBaMcuBb9x9XIZCTEoyeQC/BD4DHkl0eQ0xs10zE2HykswF4G5C4dmY/qiqrgp5/Hh+E6AN8GYaw6qyJPPIifsdwMz6mNkSoIAKWo25cL/DtvMgR+53SCoXSPH9Xu2LvJlNSHwftfnrTOB+4FdAa8JN0L+iS1SwL5JHFraRSxFb7x4i8Y+YM4GmwD7ArmZ2fvoj3yKOHcqD0Do5HLjf3dsAq4FIvjtNwWdyOrDS3UszEnDlcezoZ/LjdeoCI4CrNmvZZEQK8siV+x13L3L3/YFioGcF/30u3O/bzIPcud+T+UxSf7+7u15JvIAmwOwK9rcHxpbb7g30jjrezWJsBawEFiVe6wn/em+42XnnAA+V2/4DcF/U8W9HHg2BReW2fw2MiTr+7czln4TW4iJgBbAGeDzq+KuaR+LcmsBY4C9Rx70Dn0fW3+8V5Na4kr+7svp+r0IeWX+/VyGXlN/v1b4lvzVmtne5zd8Csys47W3gQDNrama1CANXRmcivmS5+7vuvpe7N3H3JoQ/RIe7+4rNTl0MtDOzOmZmwAlk0QCWZPNIbC8xs2aJXScAczMb7dZVIZfe7r5f4pxuwCvunvHWVmWSzSPx5+khYJ673xlBqFtVhXsk6+93ADM7sNxmZ2B+Badl9f0OyeWRC/c7JJ1Lyu93Ffmt+5eZvWtms4B84GoAM9vHzF4AcPf1hG6XsYQb5Cl3nxNVwFW1WS5vAs8A7wDvEv585MQMU+XzSLgSKE58dq2BvtFEVnUV5JKTNsvjGKA70NG28UhqtsnR+/22RDfxLOBkwmjtXLzft5lHQi7c78nmklKa8U5ERCSm1JIXERGJKRV5ERGRmFKRFxERiSkVeRERkZhSkRcREYkpFXkREZGYUpEXERGJKRV5EdkqMytL4pyfmdkkC0uxYmYnmtmwzc6pZWaTLax4JiIZoCIvIqlwEfAfd9+Q2D4MmF7+BA9Ls74MdM1wbCLVloq8iGyTmTUxs3lm9n9mNsfMxpnZz8qdUgCMKrd9GNDQzF41sxVmdmJi/7OJc0UkA1TkRSRZBwKD3P1g4GvgbAjd8MAv3X1RuXMPAz53918DPdhU2GcDR2QsYpFqTkVeRJL1kbvPSPxeSlh+GaAeoegDYGY1gT2AOxK7dv7xeKI7/wcz2y0TAYtUdyryIpKs78v9voFQvAG+A3Ypd6wlMNPdNya2D+WnyzTXBtamK0gR2URFXkR2iLt/BexkZj8W+sOAmeVOORSYBWBmewKfufu6zEYpUj2pyItIKowDjk38fhiJop5wCJta8vlA2tbOFpGf0nryIrLDzKwN8Bd3776N8/4D9Hb3BZmJTKR6U0teRHaYu08HSn6cDKciiVH4z6rAi2SOWvIiIiIxpZa8iIhITKnIi4iIxJSKvIiISEypyIuIiMSUiryIiEhMqciLiIjElIq8iIhITP1/wiE3yNr/x6EAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "from scipy.optimize import fsolve\n", "\n", "def multipas(phi, tt, y0):\n", " h = tt[1] - tt[0]\n", " uu = [y0]\n", " uu.append(sol_exacte(tt[1]))\n", " for i in range(1,len(tt) - 1):\n", " temp = fsolve( lambda x : -x+uu[i-1]+h/3*(phi(tt[i+1],x)+4*phi(tt[i],uu[i])+phi(tt[i-1],uu[i-1])) ,uu[i])\n", " uu.append( temp )\n", " return uu\n", "\n", "H = []\n", "err_mp = []\n", "N = 10\n", "for k in range(7):\n", " N+=20\n", " tt = linspace(t0, tfinal, N + 1)\n", " h = tt[1] - tt[0]\n", " yy = [sol_exacte(t) for t in tt]\n", " uu_mp = multipas(phi, tt, y0)\n", " H.append(h)\n", " err_mp.append(max([abs(uu_mp[i] - yy[i]) for i in range(len(yy))]))\n", "\n", "print ('Multipas ordre=%1.2f' %(polyfit(log(H),log(err_mp), 1)[0]))\n", "\n", "figure(figsize=(8,5))\n", "plot(log(H), log(err_mp), 'r-o')\n", "xlabel('$\\ln(h)$')\n", "ylabel('$\\ln(e)$')\n", "grid(True)\n", "show()\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.5" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": true, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": {}, "toc_section_display": true, "toc_window_display": false } }, "nbformat": 4, "nbformat_minor": 4 }