{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from IPython.display import display, Latex\n", "from IPython.core.display import HTML\n", "css_file = './custom.css'\n", "HTML(open(css_file, \"r\").read()) " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Python version 3.10.6 (main, Nov 14 2022, 16:10:14) [GCC 11.3.0]\n" ] } ], "source": [ "import sys #only needed to determine Python version number\n", "print('Python version ' + sys.version)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# M62_TD2 : Exercices schémas multipas et Runge-Kutta" ] }, { "cell_type": "markdown", "metadata": { "toc": true }, "source": [ "

\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Rappels schémas multipas\n", "\n", "
\n", "\n", "Considérons une méthode numérique à $q=p+1$ étapes (=pas) linéaire.\n", "Elle s'écrit sous la forme générale\n", "$$\n", "u_{n+1} \n", "= \n", "\\sum_{j=0}^p a_ju_{n-j}\n", "+\n", "h\\sum_{j=0}^p b_j\\varphi(t_{n-j},u_{n-j})\n", "+\n", "hb_{-1}\\varphi(t_{n+1},u_{n+1}),\n", "\\qquad\n", "n=p,p+1,\\dots,N-1\n", "$$\n", "où les $\\{a_k\\}$ et $\\{b_k\\}$ sont des coefficients donnés et $p\\ge0$ un entier.\n", " \n", "En notant $\\varphi_{n-j}\\stackrel{\\text{déf}}{=}\\varphi(t_{n-j},u_{n-j})$ et en écrivant explicitement les séries, on a\n", "$$\n", "u_{n+1} \n", "= \n", "a_0u_{n}+a_1u_{n-1}+\\cdots+a_pu_{n-p}\n", "+\n", "h b_0\\varphi_{n}+h b_1\\varphi_{n-1}+\\cdots+h b_{p}\\varphi_{n-p}\n", "+\n", "hb_{-1}\\varphi_{n+1}.\n", "$$\n", "\n", "
\n", "\n", "1. La méthode est **consistante** ssi\n", "$$\n", "\\begin{cases}\n", "\\displaystyle\\sum_{j=0}^p a_j=1,\n", "\\\\\n", "\\displaystyle-\\sum_{j=0}^p ja_j+\\sum_{j=-1}^p b_j=1\n", "\\end{cases}\n", "$$\n", "\n", "1. La méthode est **zéro-stable** ssi\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", "$$\n", "où $\\varrho$ est le premier polynôme caractéristique défini par\n", "$$\n", "\\varrho(r)=r^{p+1}-\\sum_{j=0}^p a_jr^{p-j}\n", "$$\n", "\n", "1. **Consistance + zéro-stabilité = convergence**\n", "\n", "1. Une méthode convergente est d'**ordre** $\\omega\\ge1$ si, de plus,\n", "$$\n", "\\sum_{j=0}^p (-j)^{i}a_j+i\\sum_{j=-1}^p (-j)^{i-1}b_j=1 \\quad i=2,\\dots,\\omega.\n", "$$\n", "\n", "1. 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}\\xrightarrow[t\\to+\\infty]{}0.$\n", " Si, sous d'éventuelles conditions sur $h$, on a $|u_n|\\xrightarrow[n\\to+\\infty]{}0$ alors on dit que **le schéma est A-stable.**\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercice : consistance\n", "\n", "On notera $\\varphi_k\\stackrel{\\text{déf}}{=}\\varphi(t_k,u_k)$.\n", "La méthode multipas suivante est-elle consistante?\n", "$$\n", "u_{n+1}=u_{n-1}+ \\dfrac{h}{4}(3\\varphi_{n}-\\varphi_{n-1})\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Correction\n", "\n", "C'est une méthode à q=2 pas explicite : \n", "- $p=1$\n", "- $b_{-1}=0$\n", "- $a_0=0$ et $a_1=1$\n", "- $b_0=\\frac{3}{4}$ et $b_1=-\\frac{1}{4}$\n", "\n", "La méthode n'est pas consistante car\n", "$$\n", "\\begin{cases}\n", "\\displaystyle\\sum_{j=0}^p a_j=1,\n", "\\\\\n", "\\displaystyle-\\sum_{j=0}^p ja_j+\\sum_{j=-1}^p b_j=-1+\\frac{3}{4}-\\frac{1}{4}\\neq1\n", "\\end{cases}\n", "$$\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercice : convergence\n", "\n", "On notera $\\varphi_k\\stackrel{\\text{déf}}{=}\\varphi(t_k,u_k)$.\n", "La méthode multipas suivante est-elle convergente?\n", "$$\n", "u_{n+1}=-4u_{n}+5u_{n-1}+ h(4\\varphi_{n+1}+2\\varphi_{n})\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Correction\n", "\n", "C'est une méthode à 2 pas implicite : \n", "- $p=1$\n", "- $b_{-1}=4$\n", "- $a_0=-4$ et $a_1=5$\n", "- $b_0=2$ et $b_1=0$\n", "\n", "La méthode est consistante car\n", "$$\n", "\\begin{cases}\n", "\\displaystyle\\sum_{j=0}^p a_j=-4+5=1,\n", "\\\\\n", "\\displaystyle-\\sum_{j=0}^p ja_j+\\sum_{j=-1}^p b_j=-\\big(0\\times a_0+1\\times a_1\\big)+\\big(b_{-1}+b_0+b_1\\big)=-5+4+2+0=1.\n", "\\end{cases}\n", "$$\n", "\n", "Le premier polynôme caractéristique est\n", "$$\n", "\\begin{aligned}\n", "\\varrho(r) &=r^{p+1}-\\sum_{j=0}^p a_jr^{n-j}=r^2-\\big(a_0 r^1+a_1r^0\\big)\n", "\\\\\n", "&=r^2-(-4r+5)=r^2+4r-5=(r-1)(r+5)\n", "\\end{aligned}\n", "$$\n", "dont les racines sont $r_0=1$ et $r_1=-5$.\n", "La méthode n'est pas zéro-stable car $|r_1|>1$, donc la méthode ne converge pas." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercice : convergence\n", "\n", "On notera $\\varphi_k\\stackrel{\\text{déf}}{=}\\varphi(t_k,u_k)$.\n", "La méthode multipas suivante est-elle convergente?\n", "$$\n", "u_{n+1}=-u_{n}+u_{n-1}+u_{n-2}+ 4h\\varphi_{n}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Correction\n", "\n", "C'est une méthode à 3 pas explicite : \n", "- $p=2$\n", "- $b_{-1}=0$\n", "- $a_0=-1$, $a_1=1$ et $a_2=1$\n", "- $b_0=4$, $b_1=0$ et $b_2=0$\n", "\n", "La méthode est consistante car\n", "$$\n", "\\begin{cases}\n", "\\displaystyle\\sum_{j=0}^p a_j=-1+1+1=1,\n", "\\\\\n", "\\displaystyle-\\sum_{j=0}^p ja_j+\\sum_{j=-1}^p b_j=-(0\\times(-1)+1\\times1+2\\times1)+(0+4+0+0)=-3+4=1.\n", "\\end{cases}\n", "$$\n", "\n", "\n", "\n", "Le premier polynôme caractéristique est\n", "$$\n", "\\varrho(r)=r^{p+1}-\\sum_{j=0}^p a_jr^{n-j}=r^3-\\big(-r^2+r+1\\big)=r^3+r^2-r-1=(r-1)(r^2-1)=(r-1)^2(r+1)\n", "$$\n", "dont les racines sont $r_0=1$ de multiplicité $2$ et $r_1=-1$.\n", "La méthode n'est pas zéro-stable car $\\varrho'(r_0)\\neq0$, donc la méthode ne converge pas." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$\\displaystyle \\left(r - 1\\right) \\left(r + 1\\right)^{2}$" ], "text/plain": [ "(r - 1)*(r + 1)**2" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import sympy as sym\n", "sym.var('r')\n", "rho=r**3+r**2-r-1\n", "sym.factor(rho)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercice : schéma à deux pas d'ordre 4\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 d'ordre $\\omega=q+2=4$. \n", "\n", "Construire un schéma implicite à $q=2$ étapes consistante et zéro-stable d'ordre $\\omega=q+2=4$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Correction\n", "\n", "Un schéma implicite à $q=2$ étapes s'écrit\n", "$$\n", "u_{n+1}=a_0u_n+a_1u_{n-1}+h\\big(b_{-1}\\varphi_{n+1}+b_0\\varphi_n+b_1\\varphi_{n-1}\\big)\n", "$$\n", "\n", "La méthode est d'ordre 4 ssi\n", "$$\n", "\\begin{cases}\n", "\\displaystyle\\sum_{j=0}^p a_j=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", "ssi\n", "$$\n", "\\begin{cases}\n", "a_0+a_1=1,\n", "\\\\\n", "{}-\\big(0a_0+1a_1\\big)+\\big(b_{-1}+b_0+b_1\\big)=1\n", "\\\\\n", "\\big(0^2a_0+1^2a_1\\big)-2\\big(-1b_{-1}+0b_0+1b_1\\big)=1\n", "\\\\\n", "{}-\\big(0^3a_0+1^3a_1\\big)+3\\big((-1)^2b_{-1}+0^2b_0+1^2b_1\\big)=1\n", "\\\\\n", "\\big(0^4a_0+1^4a_1\\big)-4\\big((-1)^3b_{-1}+0^3b_0+1^3b_1\\big)=1\n", "\\end{cases}\n", "$$\n", "ssi\n", "$$\n", "\\begin{pmatrix}\n", "1&1&0&0&0\\\\\n", "0&-1&1&1&1\\\\\n", "0&1&2&0&-2\\\\\n", "0&-1&3&0&3\\\\\n", "0&1&4&0&-4\n", "\\end{pmatrix}\n", "\\begin{pmatrix}\n", "a_0\\\\a_1\\\\b_{-1}\\\\b_0\\\\b_1\n", "\\end{pmatrix}\n", "{}=\n", "\\begin{pmatrix}\n", "1\\\\1\\\\1\\\\1\\\\1\n", "\\end{pmatrix}\n", "$$" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWIAAAAyCAYAAABmpdE5AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAMR0lEQVR4Ae2d7ZXVNhCGL5wtgJAKsnTARwUhHeSjgpAOyOEf/zihA5IKSOggoQISOghUENgONu/jlXxlW/aVr+XvmXO0kmVppHklj8cjXe2t6+vrQxc9f/78UvdfUUbpb7rK2j1DwBAwBAyBGwSkL+8r9VbhhdIvu3C53XVTlZ/q/r8KHxS+6ypr9wwBQ8AQMASOCEh/vtcVevOZ0v8qYNRG6VabRaxKv6gGivgnpX+N1rZMQ2AlCGgOY50wl39aSZezd9Nh8IcYP1D6KnsDO2HYF0eVvyNosIxRxGCPYVuhi8qVu1DBx0qihH9V2pRwDCTLWxsCKCAslF2RUwK/SehPCg8VWq2yXQHTU9ghOKrulcLXavKjAm7ehou3zTXxM/1U5d1aD8hvtA0ENI8xKnZJkh0l8J17ll/vEoQMQg/FkfrqBkbtY6X5OqtQmyLGIv6rUtIuDIEVIuAmPQ+BfYqvcPw21mX/IkS/VqihiDVx8WdANnFvcLC/60bgB81pc6+tewy30nvvG75XF6ihiIMC+JSMDIHVIiAFjEui2Hq5WiGs47tAoEsR7wIAE3KbCEgJX0oy/KPeCtmmoCbVJhAwRbyJYTQhIgjYtssIKJa1TARMES9zXKxXAxCQFfxE1c0lMQBDqzotAqaIp8XbWhsZAeeSuGMuiZGBNvZZEYj+oCNrC8bMEJgWAXzDj6SI+QFHSOzdvHT5HxQXe+XDApY2BOZCwBTxXMhbu6MgIAXL/vfGHnjlfyZfsZ2ZMgryxnQIAuaaGIKe1V0TAuyP93vk19TvnH390jG7m5PpDnllx9Es4h3Ooj2JLAuYRTvcFRA/L8Vl8U5x57GERemN/HEyI43/RdcfymNb35+K7ccuieM8Jo6bVcQCjdPjoP8U+CXLL27ykbcLkrz4RUc7bUv8sTDhz2Eyn3Td+MWQ8mcl9WmS81KWjIX6Nqk7ZslYDJmMY+K4SUUswP4R4BzG/AbgFaMw/lH8jcKmN/g7WSc5bUttXQlXMAXvvxV2S4bFcegNiyMWqanN+Yg1CdhDyvalQgkDhJsYXG9+bymyKkx92haW959gbXQwLI6TwLA4YhGm7oYXpGOK2PvTsHbWSHyGxc6dfad8fIRYx0aZEBCe3u/Y2KmQqYnVsDEsjkNlWByx8Clh4nVqQwfFFLF/sNZq4dD/2IFF3iXh5fP4WDwMAQ65Zl+un2TDuK27tmFxHD/D4ohFmMJIbOigi7CESz9TzIOVbOGoLFY0i2MoOxbHcAN8q/xsK9Pixf/OYx9o6+KL7jXeNKpTp8ZnQb1A32u1uwj5+/Y7U3km1Xth4A9ff6Tr17ouXUND2hGfk+M+hH/muobFEdBRsTg20z8185x6oR6za+WJQrljpWIR64b3oTb+lUebuKoD4CzW/Kw0AeULH3YppChGFe0mxwdlx+p8F3kl22WdZemT74T6tiT5fbcmid244AdkbPi3Woz9jwpMNPIGkeOfMu6D2slR2bA4ojg2FseW+qfmnlNqHwMFBfxK6fIZue1FUSbK83uFr5X2n/H+djRWOZQa25dQwGEd0o3PVZVBOROeEkeZRjJVFsX6heIHkdt9s/xm7L71GuXVn77ys4jIm5AXVzKpfE75k9tNKMhLCGLO0MeDi0l/yzWkvCWM+01nxvubisXW5kAM0SQsYhVT8zSnVouj+s5XPcqYnVwYMocLBCJDAd/FV7ouHiilUwhlCiClie0qMRCVT1OVYTfDf4oLd4ViFs7YUJ5kfatcSr9ivmHXpYO3lnGd5KI+8gO4n6Bg3osS5e/FM0Nhxg63RH1skK944eneFONeiqL2rsuLEwmVvXWiSJ/bKVhMOgeWjEUfYOtlJdfqcZQMHNPKOhxfj68uAiG7lFhQrJLEgq74ksWUhxCTu77Yh9LiP5kWpHL4e1HEHMQSWtO+SO9YfK4UqBdTdD4vS1uuc8nyq1+86FBapaXoeKw5ir1w/cuGXSrQ6ON+08zNX+GbU7mGrE+lT2Ix9RxYMhanwOy6vzEcCwMRi/hKQt9TjGvio2I+M1EanaQyKDZCXeGinA66XypopVHMlK0rQdrm7VbPV9bZRLu0V6dCYGWW/aoX6HPdR/4+fFdWFpy9wvVd97/i4kU75bj79ueKO7GYq1MztXsSC80NvpRwNfI1/reCf4Fzch774J8qD/3AfML1eVInqdxqSPLg0kXmQt/e9j3XDfwWvyu8VRogU6muRItPNCqLDzsn4IUSjhFWuFeQsfvn5CFgbFGPQa98Rqtvbf3q026K/H34TVY2t/yOHy9iPrt4iNrwHWPcB+OWAY9yLkSwGNy/KRmMiUXAm2eVL6aHyntJUBqXJZ/rXOPyxNB7prAZklzIzJdxafSWihgpVQBlzEMCQJ2ksjxoFetTeTDHwuUNB536STFKuO1hLRjwR3zxQ18rnFzkUhkG75Pi0gWgNG2gIH5UKMjlfVZ8kqevE8aql0P+kGVrmr4qJMlfY+IXJhsvO3iq7Nnyu3ZYBwhfesyb38W7vmbgipdR9nEvOZ+ZyIDHuVgk9Zj+KZwzB5L4h4UmwILxx+i7r8CifjhfmJds8fLE+SXlC85nnhtPiWOsj2ofKxhrH2OltPIrithVxEVx31VwWa0Rnw18SvhdEACGRXyXPMXwglBaMQL0kyCLF/Up5xV8jFeYh/XLS6DYoaH0bwrl24eCAU8mOP04h1LlP4d3WSfoa5L8Ko9FgVLk8w8qrpXnrw9Ke0yHyM+LzePMWNMOL3NPtBEj8B5j3GNtJeVlwOMUFkn9aCsU9C9pDrTxSckP2jp3bnRiIf5+R9Uj9ee175PyUcxXiksFpWsU1ztfZmgcyDY6ji19Lf4hgfoRvnwOF5HC3oeKQvXpSLFSmaGM6lTJU6MeeNwUIcjUq1/XeRXX4pF8spcDO1QIrTxVtrSco4U6Ml07FVld8VheB6fTt3rKn9Q+PEeSvxBIvCcd9xBFtc1c83MAxc81L+bWOa17Z+Ohurx0knBXubOI/p1TUfWWigVKNrR+uS4VpOs3i/lvFDAOk3TFKYzEZzIcI31BxoYcMUXsLRUmb04C8LITAoM0APv2crbVhxcW/Zs+FTZWdmz5Jx93jSdzlwUer4gPSvPCZZcOC0Fd4z02HpNOH8m6SCzUL14OWNyhUnrEGAUAoSP8i5OxLMczKDNJciCORR8dD9INnXd7EinUiDqBI/5Lxd5dgPXAJ8xs5IDJua+4VRa1xZsdJz0LD0Va17hvZqMp5Fcbc4w7Lpgnajv82vEPNPhHaWw8xH+OObBILDQAKOL6C7Gex5jhquA5Kf/H4JpwVL9j9KmeGbOI62WyXQvAEsxsTIcx4mFFUYxOaoe34C7ln2HcsbJwFRAKUh94oN1VazTqfJhpDiwVC5SsfzkWAyJ8WNspyeHVcPesDMdSnq7EpIq4qyNz3NOATqKE55Atpc2tyi+5eMC/CDFQnreO/QJyeLtIbxEPw6IxzGdlnItjamOTuSZSO2TlDIHcCOghwteIW4gtQ5XV6txtLZ2fYZFnhHLjuGuLOM+QGJelIqCHhe1QKGEWgfhEL1fkld4VGRZ5hnssHG9dX1fPSFFDrLJ+VuBYw9lWKfPAZlwMgRsENJdxTbC3+tSuic1DZljkGeK+OKp8q24110SeMTEuC0dADwEr9Cze8aMTHojdkmGRZ+hz4miKOM+YGJcFIaAHhM3/uCXq5F0TuCt2QYZFnmEeG0dTxHnGybgsCwHOD+HQ7V1bvm5IDIs8c3NUHG2xLs8gGZdlIYALgmM4iUN66C4q+1fDAhtMGxZ5BnVUHE0R5xkk47IsBBo/nJFSZrHujoI/onNZPR6vN4ZFHmxHxbFr1wTnQDR+1ZJHJuNiCIyLgOYufuBw/l7quvPQn3F7NB93wyIP9kNxVH3mIP+VvLEjraGI6bIqsKeNTztOYDMyBAwBQ8AQGIiA9CkLyPiaOfS+YmG3LdbhQ3uswnzKGRkChoAhYAgMR8Dv1inPYPYs2xSx19atJ1V5BhYbAoaAIWAIdCPgjFr0Kf+ujV95ViiqiF1BDsThP2/4w1IqFe3CEDAEDAFD4DQCTgm/dSXDdYuyclQRc1eVsYoJ/BKp9cQqyhoZAoaAIWAINBGQ7sQdgV8YN+8DXTcOhadWdLGOG55U8VLpQhErbYt3HhiLDQFDwBDoQED6ksU5LOEXSnceufs/FmzC2EyCG64AAAAASUVORK5CYII=\n", "text/latex": [ "$\\displaystyle \\left\\{ a_{0} : 0, \\ a_{1} : 1, \\ b_{0} : \\frac{4}{3}, \\ b_{1} : \\frac{1}{3}, \\ b_{m1} : \\frac{1}{3}\\right\\}$" ], "text/plain": [ "{a₀: 0, a₁: 1, b₀: 4/3, b₁: 1/3, bₘ₁: 1/3}" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import sympy as sym\n", "sym.init_printing()\n", "\n", "sym.var('a_0,a_1,b_m1,b_0,b_1')\n", "\n", "eq1 = sym.Eq(a_0 + a_1, 1) # eq1 = a_0+a_1-1\n", "eq2 = sym.Eq(-a_1 + b_m1 + b_0 + b_1, 1)\n", "eq3 = sym.Eq(a_1 + 2 * b_m1 - 2 * b_1, 1)\n", "eq4 = sym.Eq(-a_1 + 3 * b_m1 + 3 * b_1, 1)\n", "eq5 = sym.Eq(a_1 + 4 * b_m1 - 4 * b_1, 1)\n", "sym.solve([eq1, eq2, eq3, eq4, eq5])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# monter en ordre pour un 2 pas\n", "# import sympy as sym\n", "# sym.init_printing()\n", "\n", "# sym.var('a_0,a_1,a_2,b_m1,b_0,b_1,b_2')\n", "\n", "# eq1 = sym.Eq(a_0 + a_1+ a_2, 1) # eq1 = a_0+a_1-1\n", "# eq2 = sym.Eq(-a_1+2*a_2 + b_m1 + b_0 + b_1+b_2, 1)\n", "# eq3 = sym.Eq(a_1 +4*a_2 - 2 * b_m1 - 2 * b_1-4*b_2, 1)\n", "# eq4 = sym.Eq(-a_1-8*a_2 + 3 * b_m1 + 3 * b_1+12*b_2, 1)\n", "# eq5 = sym.Eq(a_1 + 16*a_2- 4 * b_m1 - 4 * b_1+8*b_2, 1)\n", "# sym.solve([eq1, eq2, eq3, eq4, eq5])\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On reconnait la méthode MS$_2$.\n", "\n", "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-1)(r+1)\n", "$$\n", "dont les racines sont $r_0=1$ et $r_1=-1$.\n", "La méthode est zéro-stable 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", "$$\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercice : convergence\n", "\n", "On notera $\\varphi_k\\stackrel{\\text{déf}}{=}\\varphi(t_k,u_k)$.\n", "Soit la méthode multipas\n", "$$\n", "u_{n+1}=(1-\\gamma)u_{n}+\\gamma u_{n-1}+ \\frac{h}{4}((\\gamma+3)\\varphi_{n+1}+(3\\gamma+1)\\varphi_{n-1}).\n", "$$\n", "1. Pour quelles valeurs de $\\gamma$ est-elle consistante?\n", "1. Pour quelles valeurs de $\\gamma$ est-elle zéro-stable?\n", "1. Pour quelles valeurs de $\\gamma$ est-elle convergente?\n", "1. Pour quelles valeurs de $\\gamma$ est-elle convergente d'ordre $2$? et d'ordre $3$?\n", "1. Pour $\\gamma=\\frac{1}{2}$ vérifier empiriquement la convergence sur le problème de Cauchy\n", "$$\n", "\\begin{cases}\n", "y'(t) = -4y(t)+t^2, &\\forall t \\in I=[0,1],\\\\\n", "y(0) = 1\n", "\\end{cases}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Correction\n", "\n", "C'est une méthode à 2 pas implicite : \n", "- $p=1$\n", "- $b_{-1}=\\frac{\\gamma+3}{4}$\n", "- $a_0=1-\\gamma$ et $a_1=\\gamma$\n", "- $b_0=0$ et $b_1=\\frac{3\\gamma+1}{4}$\n", "\n", "\n", "1. La méthode est consistante pour tout $\\gamma$ car\n", "$$\n", "\\begin{cases}\n", "\\displaystyle\\sum_{j=0}^p a_j=1-\\gamma+\\gamma=1 \\text{ pour tout }\\gamma,\n", "\\\\\n", "\\displaystyle-\\sum_{j=0}^p ja_j+\\sum_{j=-1}^p b_j=-\\big(0\\times(1-\\gamma)+1\\times \\gamma\\big)+\\big(\\frac{\\gamma+3}{4}+0+\\frac{3\\gamma+1}{4}\\big)=1 \\text{ pour tout }\\gamma\n", "\\end{cases}\n", "$$\n", "\n", "1. Le premier polynôme caractéristique est\n", "$$\n", "\\varrho(r)=r^{p+1}-\\sum_{j=0}^p a_jr^{p-j}=r^2-(1-\\gamma)r-\\gamma=(r-1)(r+\\gamma)\n", "$$\n", "dont les racines sont $r_0=1$ et $r_1=-\\gamma$.\n", "La méthode est zéro-stable ssi\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", "$$\n", "donc ssi $-1\\le -\\gamma\\le1$ et $-\\gamma\\neq1$, i.e. ssi $-1< \\gamma\\le 1$.\n", "\n", "1. La méthode est convergente ssi elle est consistante et zéro-stable donc ssi $-1< \\gamma\\le 1$.\n", "\n", "1. La méthode est d'ordre au moins 2 pour tout $\\gamma$ car\n", "$$\n", "\\sum_{j=0}^p (-j)^{2}a_j+2\\sum_{j=-1}^p (-j)^{1}b_j\n", "{}=\n", "\\big(0^2a_0+1^2a_1\\big)-2\\big(-1b_{-1}+0b_0+1b_1\\big)\n", "{}=\n", "\\gamma-2\\big(-\\frac{\\gamma+3}{4}+\\frac{3\\gamma+1}{4}\\big)\n", "{}=1\n", "$$\n", " Pour que la méthode soit d'ordre 3 il faut que\n", "$$\n", "1\n", "{}=\n", "\\sum_{j=0}^p (-j)^{3}a_j+3\\sum_{j=-1}^p (-j)^{2}b_j\n", "{}=\n", "{}-\\big(0^3a_0+1^3a_1\\big)+3\\big((-1)^2b_{-1}+0^2b_0+1^2b_1\\big)\n", "{}=\n", "{}-\\gamma+3\\big(\\frac{\\gamma+3}{4}+\\frac{3\\gamma+1}{4}\\big)\n", "{}=\n", "2\\gamma+3\n", "$$\n", " donc ssi $\\gamma=-1$. Cependant cette valeur n'appartient pas au domain de zéro-stabilité. \n", " On conclut que la méthode est convergente d'ordre $2$ pour $-1< \\gamma\\le 1$ et non convergente sinon.\n", "\n", "1. On définit la solution exacte (calculée en utilisant le module `sympy`) pour estimer les erreurs:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%reset -f\n", "%matplotlib inline\n", "\n", "from matplotlib.pylab import *\n", "# rcdefaults()\n", "rcParams.update({'font.size': 16})\n", "\n", "# variables globales\n", "t0 = 0\n", "tfinal = 1\n", "y0 = 1\n", "\n", "phi = lambda t,y : -4*y+t**2\n", "\n", "gamma=1/2" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAKgAAAArCAYAAADsb8PCAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAHO0lEQVR4Ae2c4XHVOBDHnUwKAK6DXAdwVHChA+AqADqA4VPyLZPrIKSCO+jgoIKQdBA6gEsHuf/PSB5blv1kP8t27u3OKLJXK2m9+mu1kp23d3d3V6TSycnJY8leKB0q/a37N6l1Tc4sMMYCB0MqCZDXkn+iHFT/M6SuyZoFxlhgf2glgfPI1fk8tK7JmwWGWmAwQNXBM6VvAurt0M5M3iww1AJjAIoHNe851NImP8oCvTGovCSboTOlb0rflQAmG6VTJSOzQHYLdALUxZofpQGbIgBaKPcbI/Og2YfGOsAC0SVeQHygMsD5zoMTYRFAvRbP4s/SHPYntwWiAFWnLOsPBMQPgQIWfwYGsdu8FugC6Et121jGnVclJvXLfF7NrHWzgCzQikEdEFniQyAC2kLlDeDCM0q3gLPve1eDCQ+9Et/Cpp+2aPxtAbRWWm6Mavecf/ImCZA+51p5KEOxUb8FzmS36hWxrs8lfqX0a3+13SxtLfEyGDMZL+lndyEegOR46asS9Ey8rOBU+1X/P7uM/02Vi9dehPtaOhPLeyLePxQP+/4vKXWMYnItgDoLvVD+VBXeKvlzUDzoI3jKmfXZyPWROmAMLjrdF8J7+ol+X3Qeree2Y7k35Gum0VoOqKgHwlszOd6lVnNGuFUenjqkNrGYnHTGATxXvoolXnrgGN4or8KQscZRG1uP5aoAqgdic/ZF+ZOYUVw58VorxFAZ/N+Vr2az0acvz6dywOBfhqxCb+l0I53YX7CKjibVn2Qsu5b40YptWRFv0hc+ELsRm/6I9EM96q+JOvXVAPIc6MuburWAc8pQaZKxXBtAX2qw+pbpzi+pXD3qM3PXQlF9pSPg5C0dKwGhCXE0vMVI/ePNmShTTZZJxnI1AJWBiFc2nQzgkT4pdRH1y/PaLoGZ+S19HRDx9ueAwgGDeDu2Ksyp7h/Spc85JOuidiYby4PkXgcKSkk82Wslgv+r+sO7sgvl9TgHb9N6CSAZBpmAnfbwMgwqcdul8j+V14n6tDOJoesNp14n6EuszLOQV6R6W29KqsYGXqjv3pMZlS82ltkAKhu914OxjDGbLpTqoMHLwa/Tb7ppxZ+qD+g+u3aOlAPALiLAB9BJpLboL1neNdq7gdikr8ofJik3k5D0YdITZvStXouNZRaA6mGJZy6djQFUuHzFvCWzNJRzTZQZda7rjMg19ZNjOemZ02ul6Bt5hNlZHCl1HuktPZZZACoT8y8hHkx4y9PA7HitkPdIvNtArn5Lnb74E1m8AEBfA6Xom6yn7MlzfVEa8nwvauPQ6ktlhGCtVSsQXHQsD6Rk+v8dB5qHt2prD57yEmjKGSQMWi3v4uFd4bXiTfGipDrI4xnDD1hC+U0gD+Wz3A/QN7l/tYlNo+fDyY3UBNUe9uSTyr6lvXD9ki8ylgC0BFVN9ykv2QSFHzjzoMQ83sP6/lieAWKMqIORKlDrGuMyaHWiPu0kkepPHoO6jlP1TdIzkxAA5Y0dG8464UA49oKP9/TL/yJjeVDXLMM1RghnaCz+pGvkkI9RI56T0fwxRghyPGjYX6y9kqd2csWgqfp26pa7QM9ebj7DfsT/Vzw2pfUTFsQWGcv9UMGJ7xtg0UPjWUixpRqwPe3ovwKe2sBLMsNDcFKVJTDGp2xOStV3Tp1S+8K+pJAWGcvcHpTlgfNOXnt9V/IfRFRLdc0Kf+k6XG58Me1wsE1QXygPzz+9HODP5RV9Hyl5qr4pbc0iI5sS7vgV7Ej3jEX9rHmZseRrprnS8fHxmdJVV38qu1F63FXex1e9Q+r3yVjZdGMtW88ylr0eVLOIgJlDdmbWoB8LU128JjOx3HkqZ9nAA75S6iLq4AHHeEFmOPWNJrbAkmPZG4NKMXbgAAxwNeJGAKd0o+SXhdAsbGRYtj0B9A+S7zzLVBnHUYM/nHA6UK86zvKdWj6JBRYby14PyqNp0InroDBuhA84f1AYITya/wKf2JMYMmwjUq1g90j8w044lYifxnjd1PZ3XW65sdwUl7lYoxXbiX+eK+ZTu8STbzfpRjlyyKfImsx0MWiqLbcdy41f1Mvr8dXNV+UND6V7Psz4pJzZZWQWyGKBxhIvsLFks9HgzItjIZZkNkqnSoXKWdYBKjEpsn2fvqnYyCywnQUqD+rAR+wX/lgYoHyo8lvfla4Jmj8qz/ma1Hdn+Q5bYJ9nF9DwiIAz9cfCGq/yaMPILJDDAiVA1TDLOsdG4TEN3jO28+7i59DR2txhC3iA8s1mA4jOqxJnts4/xWvxd9iG9ugZLbDvgNg6iFefgLZQeQO4YuE9G3zXBmwjs8CkFvAelEYbX6vovoozBUB++QKvCVV8bihT5stgGZkFJrMAHpTdOV6yApkDHcdL/jeE6r/kcZ8/JZvMcNbQPBYoj5kESJZ43pVfKv2ixDt0gMvmCR4fsJbfWSoHyLxaZNdf6D7cWME2MgtMYoH/AMXi2HRj1MYPAAAAAElFTkSuQmCC\n", "text/latex": [ "$\\displaystyle \\frac{d}{d t} y{\\left(t \\right)} = t^{2} - 4 y{\\left(t \\right)}$" ], "text/plain": [ "d 2 \n", "──(y(t)) = t - 4⋅y(t)\n", "dt " ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQUAAAAuCAYAAAAhigzQAAAACXBIWXMAAA7EAAAOxAGVKw4bAAALzklEQVR4Ae2d7XEUORCGF5cDML4IDjIAEwEmA+AiADKA8j/+uSADjggMZAAXgTkygAzOOAPufWRJaGakmdnd+Vq7u2qsb6n1jtTqljTrW79+/VoZGQJzIfD69esDtX3i27/j3WeKv5yLp5ve7v5NB8D6PzsCbyQAXgQu5H8n/7967oY4c6dFYG/a5qw1Q6CBwHMJguMk9o38dxR3L4kz74QImFCYEGxrKosAWsLXbIpFzoLALdtTmAV3a7SAgDQENIXHcs18KGC0SbTXvD6q7H35W/drbE9hE4StzEYIaDCyqch+wSP5f9Qr8QP3seLv19OWEu7qw1L4hA/P63t5L/Qc6QkbufKWyYRCGRtLGR4B9g4YmAzSCmkAE4+W0LmSVQpOHyj2YXpW2lsUpmgET8gl/0s5vfZpTCiAmNFUCDxSQz/8YI1tKoxAeCWX9JUP4za0iVhoPk+2D/OxM3zLJhSGx9RqLCPAKvspTfYCgGNIhEJYydh8fJXmW5C/0YcF8TYIKyYUBoHRKikhoInOJGKSs5+ARnBPcWx4nct9K5c9hrDXIO8VKS3eXQhxc7k9+jAXa6O0a6cPo8BqldYR0MRiA/Gj3Fv1tF0J73IfxDt7CuzZ3JafvYYi7RVTLMEQGBYBbPFvw1Y5eW3XoQ+doJlQ6ITIMgyEAGbEl4Hqmqua69CHTuxMKHRCtBsZpBJimzeIeD3Y8rOR2g/7CZ9nY2LLhq9DH/pCYEKhL1ILzqcBmz2D9gOZjb25j/ZYYVfiI2oKnrcFo9pg7Tr0odGpXIQJhRwqOxSnycVgLU160pZgx1dscfHMpuOs2ssGr/g69KFXtwcRCnrJvV5w33y9OLdMK+GJWs7ZflyBgQWc9bDTzEOYz5PDHQCipqZDNegEl/hwpoTcJQirdXDY9T784TtLP1ppvzW1R6JeLqorL7y0WqW1MFj52IXzaaPtEeCz47d6mGiRFOZdcBkITeFUbuXCUMw4nYeLSO/Ex3OalLuL738n+yCsuRMCMRYgjoUZH5/l/u1ian+2uqegSlEDH8jtfftMeREilyWGavxNHhRfTDAGL+riZcKAA1HpaEUv5Pbuc1LHYF61D/bf5HJtGJ5/yq3cAVCYn9XqPJcejCmr6FogsLGm4AfiidzsF20+vfFFnOJZ2dgR/6AnnXSzAyp+EFj8NNipnicpf/Kj4XAdF4nr1HK5g5DqRbXvUyfCCCGAYDrEX2JAafCJ8F0UxiV+LX45CGwsFNQFBjGTpEQMSgZv44s4xVGO8ou4yqqJw0qLmgW/DxVu2LuK+6QH7YE8FRte4a3It0fdfQl+76pcECSEVz7M9WHMBQSN41Nh0jHdGv1SfCupjHuPcrOqZmthS+xEYIn4biMUnqpDbZOaQd74Ig6UGGB6UHexey87kRs/wz9qgsn+Zwc/CA76XVyhx2fV4cfkjhNc/MA7Jg/7BwFPfqTkXA8UtJ+r0Hp/ESg8RuMgsDh89zbppwYe9mzXxGCFadvgovzTTdofsoz6wmrLqtrnF4QvlO/DkO1vW5f4RxgEjYFTBt4NhDaGNkH6mdwgLEgzMgSKCDhNQQMGacXgYXXB3o+qok97L9f9WIOvCS2goUIrD4IA7YH6WL3qX8QpKhLlqSe2FVMm8ohfeHSnJ/K3CbDAEYKszWQK+Squ6gYPVuv/fAI4D6IlqW7wa2CoeDSJNk3Os2KOIVBFIJgPbBgySFll+PmmdJCxmofVJ5Q+kqcxOVSeif7F13Mst81O/q68CJFOUj201StvUhk786kgS5KiN0yaRl9ijsSj+lhto9qeJBW9KgPf1M/GpSsrF83kRM+sJxhq38gQaCCw7wdosD2ZxKjIKeW0Ala+er56ma7JQ3lW6k4Sj2HyduZdM0MQdn20hEbV4gscEJqcCjROYRRH/7jvT3qKB4JiZ78DEO9G1xiBffWNzcAwYBngp7X+MoDrcYeKY9UsEWW6JhqqOJNqTqIfK/W/a3/E8ah8B3pcv+Wy2tNPqNSPoIFQDjMF4mYZdx4a5pdLnemP+IHX0J+Ui4BRTjB3amOqd+t/QaY6KvcvUubwKx382SwuvQey1SlqbvWENLwO/218Km2x+Kb9xY+mEAY5AwJQo+mgNAY+cb0HsMqQP6yQ8haJwebaLuYYPwFtBX47Sf1Cq0B4OAGqMC6TImgbuTqccFSexd/gE4+5Sb/y/eM4c6M+qFzrhM6Btm6c2mAcNTS1devJ5R+Kf9WzM/juJ0BgfzPI04nKoOYCTNAkQva2yUSZlcpEQSJ/XGFDBXKZjNTTSSpfkrJtZTtXMRWGx+eqn0HfpS1wc7NL+8nxE0yzXJrFGQKLQyAVCqzu9YmR20+gE+Qr7QdQJgoRTaTKCkthT2gK9fZCWsVVHVkpW8m0WYCjPE5d4LG4Eqp98p3qWZfoH+ZCg8BFzyZCplGXRRgCQyKwl1RWmaAasKz4PLkNMSb9g6Rs6o2TXXWgDZRu0qHuReGRVjCVX/zRZzQkzvcRDhUCAz1OIMhNNahKvpYAZZ3mFPKoHrQmNJ9Z+x74MdcQqCOQagocj3EfgYHMeXr4t13RDEgKn8nP7b4cUU+fL+KYLGNpADm+snHqL9eXbysRwYAADJP/Qn6+KNv42FBlubmJEABTjmAhwrP3+4oV+2sINBGIQkEDlckQz/X9QMYur2gQVKE44nG5nFRZ8RQmPyZEkZTHmR5ycwKnWG6sBPFB30eZqKq7aJaM1Z9drFc4sanN0e0o72FbTMTXgeo4SeohzKKxlDHMnArYwRthFroKfwoTHxa6I/lZ/LijFOexEwqKcGquXLeDK5dKUaef6SkRZWAiMFLKl4uHKcobLR+BS7HIMzahecaBOXZjG9TPBKuMdYURCmh+2+wNbY0vPKg/TOzIn/zsk3H0zdGr408uAgEtPi7a8jMPucXM//d0AmRPERAVnDnf1Z/3clB9i51VGkeX7BfQUG/y+SkXjz57F15QRvqhB0BZPZxf4XAXYUGcbseK+sQN1VHf1dJxE38skLkfxmXRTLWHtcEeCF/44xSNeRwoaAgpf4zXKDjIqDIs0AimuB0QzAcS+D6fQc1eAtIkVKpgkTA3qCxKnmLO3wnv5K0w9jtpd3zCBzMpqGG7w/jCOBWOmA1br5Yjd4s5weI3qnDcog9oWBUMheulnnqVx4r4rvj6D+8w1zkNY3H7EcyHokZQrzUNqwIaxg58qafTdiafypO/sU+R1mv+G4XAXxoPFdV3gb0/F0+McTaiUceZgBArL4vcrCR+mNRslkdSXNAaUv7Ixz5g4D/m954D3K1+jq1eo4UNgXUQ0OBkkeD0h6v2qOe4cbN7nbrGziu+0IjDRENDRHvAZt9oQR2TX/GERoAwYB+kU7tRHk7G0BLc7dOwpzAmj1a3IdBAQAMQdRxNcye0RvGJsAoTDA2BibeojVHxiBbgtHHP21e5rUQZZeBdRFPYhEIrZJY4IgKYkWGSjdjMMFWLV7QE1G40BNRwJhL2edAeFJyXxAtXBfgNVATYmR5OFbr4QwNCW4vmvwmFed/jjWxdA5Dd8tTWXTQOnl+O7Nj7wMRhYz2YOVz4c7b4kjohnjBrEGLu2DTHm/LwDhommwmFHFoWNxoCGoissJzt74TZ4IHAXIjqNXHin0mH1oBAwJSYjcQLZgNmQJ2C+dDgT/kRzJw4Nk4O9+u1WNgQGBkBhAJfnKK2puRsWx/P6lWZhGnGKf3ig0mPEGPVrZDi4BPhcFhJmD7g7lCIl/pRY5YT5cOk4Pc7g7azkp/3gnt1JJktaZGGwAgIaNBhj/NUSPE/FcFFqThQKxlmCogfNkN53Bl+hg2ERqM/mXxjRiGwwK4uuI58o5E/5UH4IpTrQhdB4fZ4TFPwqJkzOwJMLp4lEoIK25z/CRInnvyo4BxLzm0K1Sf4SjwxycGTDV3Hs1y0ATQ0BEh9T+dYcW6z0e4pCCGj+RDwg5PBGuxe1PHzMEDn46zasvhhhT3Rc5GkICjiKpzET+4VH+CXallgWvkgSnncfYQCc5xcuG+f/geDxz/qvGWFFwAAAABJRU5ErkJggg==\n", "text/latex": [ "$\\displaystyle y{\\left(t \\right)} = C_{1} e^{- 4 t} + \\frac{t^{2}}{4} - \\frac{t}{8} + \\frac{1}{32}$" ], "text/plain": [ " 2 \n", " -4⋅t t t 1 \n", "y(t) = C₁⋅ℯ + ── - ─ + ──\n", " 4 8 32" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAF4AAAAyCAYAAADV5GxPAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAGv0lEQVR4Ae2c63EVORBGr10bgGEzMBkAG8GaDGAzWJMBlH/Z/1yQgdkIWMgAyADIAMjASwbec3RHqnlo5mq4D9sz7irdaUmt16eeVksje+/q6moxRGdnZ4fkXygD/2RIdu554PMQDD4RzuFfD+GxP5RJ4RfkfyN8Jzwbkr3LC4r5tcLpBOy+EVTaLO31aTyFXlFC4J/Dv8mWvoWJFRjPq64f8BScV6R/zA2HdLX4HeER/M+cTDsNOetV863bcipug7LAI3iE1AfCG/jYyUbB2xipABHkNCb4p4xFYJ/Bv3dcPAXuH8Il4TFB8O+RXgQ8srGOH7CfKdcx0X2m5mVVOHXQ+ATomDEcA4RgR4qafhITBJjgRDj+tzF9zNM6kNdSHME7cQ3qA16Njx1qFLjlEW2wgCTNrQDa1rDipIlng35rxIjQkYMqLXWuLXNb44xNZbpX7z9pUfuD51bP2wAfbfuDdl19Gq+c9m3SBOhqok7Ezh2IjsZPGulqcACuzRX0Pwian8+EndJcgRdsw6IyNV94Jq9mFzMwZGp20f61twHgupCuZ+/g4/q29X7NCniAfWjIoBpNjeZnJzQr4EH0i2GXmt03i3Oz8ZqUjwDfdpXdnUo727vMDfiwI19ivPxlEvTjte26lO0JUeh3f6D7hFx+yBz7MyvgAdazJ7fw9c3SIaA9Ia2h7cQ9v5Gi3XfxdUP0gefah4YbB55OqT2eiXgwVNeQ0GHyHaja1dE+0rdOtCvADZBzjSK31WPwjQJPZz1G9rDpnKBfnICHf0pQ09Qgd4uzpo0AD6Bqua+m2vwn8bA5qSNL2nuCb4EyKzWuXnYCvOtDg/YbsWVEYKSkrcvo4K+H/noGHvp3QK+VdHI8co2HR7Ws6bGMM2KoYjYop/GaAskPISuJyjUbbkoapqWn4CXp//bkTTVZRYyYpjHmNF4b/R1AV5oDZHw7tOvKh683qeY8o6bXPYq81IpU2vJ75tr1rGhmU9mudwv6q8ORqAF8bTCdT1WpRJOJX6iKQKB+zcyQKWrWnolR3tfWCY+bnozUzUmivyqk7ucFvP0OlIAnUfD+Irg4ltpgNx9SibYvJWu/tHNAOCa4lS8iZLWbfv98VFTgBgjRVxVU8D2uCGdFe6enp2qQA1cT/64GBruakP0PKcHbWy29/LoV6+dpB6Lt06/vfKUpqbNUhvqHLxCVVlSTo86iccciyKuorokX9cXVhW8sWaazYucqqRr1TQqmhrjPr1V6rshG02hnFEgbbbxZWXAt9+mQdjdq2w/48Co0ZXtjYQGmTLJdvZJ87UEugD4gM9ksxq4r7ZURTfnruo3XDunqfSoEEtG0A4223rQOUZ+vV1jdO5kzSKjGL0Zpc5mAd/wICL7mIx4QmdxLyGs6PNPwklDDXbIQaUeEADpPF8W1iXpcU64IxQvy2o2uUQH9dB2LN/LSG1+38bF6vRuBFLQSX96jAK9MWMZNVwQ4TCBpGz0Moz5NoxMevxrB3miKl8MaJ5o54CPY+vKRHxyZYCDg27ITor24Jo1uj7KHFIp91TEwrtI0xlrJRaV5jIyK9JL0pLXES0iN75TJAR99eDs1KQI0xyR4EfgFvLbXI+t0ywDeyXDDkzaS8JpM/fDO2T3pWUI2YhgxTXL7idsxQ6cOq8Gc0HTgiWsLt0lFdyfpgCCnybFD9E3t/0koWv8sUyPflgblNL4hsK0IA1EL4qu8rWba9frKC54hEP1wzahi6aF58DyofUNYc+R3BRWlo8WpdAFzbcAX9G3jIoAlcCV3J5XzKkiaoFZnDlrx0dFZAd9GB2DV7GBW4JPXAd/32S9sLsnvLJbtulfFZwk8wMVzouK7k1UZF92NmMe5Aq/GBq0FUL2akruTLqruWQb/qAyZIro2r6aodzsQAkiPtLXlvXcnkXFT6ceePhM0uqezAh7gRt+dpIwu6H2eyacfjXKmwKyAZ/yj7k4CtmboAc+k6fC6ktr6tWhuwGtStNNtN9EjASkdGyDjAuxRdnsxdTIuFV6H5ra4tkFcAKxA6penu5OkqdEupl5w1b7XycPDtRfYIeDv11ubAg9gpXcnPWUV/M5RN2ljfPheDDvA07m4hVYLJkeMT3OSTEpugMj88ulnq76IYdu0LfpsvB3zlYoFW/XdRQsROKrk3rbl+4CPttCTwzv6BQQqpRU/P+h3zFMW+ErQBeQFvIvPHY1AoALd+6RSckWX0eVvFnizKKzWG9zRtVd2Re4ogwBYaV7cL2ims/+5w2LZ/95hRiQqOoQPwMNvdPcW25jKE3z0/dX0c/hBl/N/8IZuItEOiEUAAAAASUVORK5CYII=\n", "text/latex": [ "$\\displaystyle \\left\\{ C_{1} : \\frac{31}{32}\\right\\}$" ], "text/plain": [ "⎧ 31⎫\n", "⎨C₁: ──⎬\n", "⎩ 32⎭" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQcAAAAuCAYAAAAlf9ztAAAACXBIWXMAAA7EAAAOxAGVKw4bAAALgElEQVR4Ae2d7XEUORCGF8oBGF8GkAHYEZzJ4OAiADKA8j/+UZABRwQuyACIAI4MIAN8zoB7H5UkJI3mY73rnZ61umqsb03rldRqtTTrW79+/Vo1Gkbg5cuXh8px5nPd9e4TxV8Ol2ypDQFbCGjM/qvnQckV8Yp7JPdHSDsInuYOIvBaoD0LOeR/Kz9g3gtxzW0IWEdA4/a5eLxf8ql4Fr/7cqNgIM/tMmMLVxF4KuBOk5TX8t9VXAfoJE/zNgTMIODHbzb5E+YY29+SsPM2zaFEpB5Ga/haT2qxDQHbCEgwOM1AXP6Tcqp4tsiM7b/0XCrMoncu1wmKJhyExhgJrAxU5QfQHwHEsfItfRwBYYkW9l7PA/nN23LEY5hYNI7JR5jt5yciSpq5fWi+b/TAZySF0SReyEVzeCX3Q0yUpwmHFA35BRAAYk94KH9HDVMcgxhJ2zHqKG52GuN/dgYTBjyv7xR1oedYDxPMPHm+mVSpHYox8VFxGPXcJPP5Zm2feICvbNJXAGZMd4RaEw5dpJCiDFIGbEYC2q0OirS8uvXynzXGQEB4oiE8ghX5q8YyA2zWWHiqSFZjhEGYeGFycaoVhMNW26d3MYlR/cfomfKi2TJej/D3FVAa44UtBbxm1IRDBocLPNRfgM3AUhigWS1IX/kwbi/w5JuBqvzPwMc+v5I9OeMjjhGNAybYtbZZ9fNeN/4mvuhQ+e6pXBAohFc+/EUuQixqDQqTjqGd96yacACFnJCkYTVwKQILwcDxJcIBMCFUyhfOZ+tPh39b7C2fG40BtIQ7aUsUh/oOMU5MkHhikruJDkMKM47RerAvBMHGcfwXPRBazyvn058mHASCgGJCMdmd5JTLmS/GMaTrG7nYIEjDjaS0uOeMkTN4JvA/A1c355Uef1Zn1PnSeD0ZCJVljDFBf/pCTFwWpDCRffT6jupAKAStA8Np2BIhzOCbdE4q4ruacBAiAoSV4JNcpP+p3ACigi49WyVcpKE/Y/wbYnWvWBHuaJEsLCd6WKG/6rkSqS7qYaJi0HSrva8fYbGxhqq6EFodweXfVV3kmnAQYgkhFKIalsQvxbt0/peCs+PTT6wwkVlYuJocTyumNkZlUPc/6mEFT8cfAoP4WagJhxx2OiOzN+TJ5kNL5988wH0MalJ/0INK/l7uHe/vy17GozFAhyrHqQ30hx5U/3AK4iJ3+acJB4+2OoH9XpDgu+yDrbxr6fxvBYQdVSKsnVFabrrK83a2FQjodYW0y6/6sG+ZodtmOJmfETpopQ6Kklp+BMZSaOn8LwVn+MQwzRZim+MjnBiYwaEJh99dke3X1fHsIdEklkJL538pOMMn24ewjUj5PvaBuMCkiQN+7sqwjeiQH4ed+F1EbEU4qAGTJtHUfLtoeOUdR4pzF5rEp9tiyC3VxkoxM1FL5z9MDtphnTg9yAyFGissJowbjIoIj5KG2scxqNP8QiHGoB5sEbONwVub/tiLGoABhRuFo4Y85QEA7hCY2lvRIeIJAUdncL+BcOfYh3irtFT+xbfDW7gyNphcTAaENMY4s30g3uDXXf2WCzF+Oh9eTW2f8jGPECDf9UAIh1nnyUbCQcwjLU/kTj6H9SBw1dRsx7uuaX8aAjccgSsLB01upPxnuQ9qGPp0DDedrxuVRvyfcmvqV626FtcQaAjsGIFNbA7sk8L5bI1t1C5UrYtKIuXCxyCV5BZlBQEJ8FM9XK1dPO1TW3bRGZsIh8cCe2hrgPW883UjjfLlKI/20cg2AvTRvvTTPrXl2kfNlYSDJjW2hrFPldEchoyUlH987S1sL2gINASuhMABpfwKjurIV2Bc7ogagU97Jze1zKIVdM5ylQeBMPR1o5IjUZ564rtiSvM0BBoCsyPghIO4ONPE5tNQNAJ+1iqdsKzuxKd0rEDH3qDyTPjerxvTCuTnyAZhMkqql3dNyptU9k3lUoGWJDVvQ6AhMIbAgSYQ98TD1U1W8tKAWNMS2LuV+dJ3UYbz6iGiPAbLURKP1U9KRwu2DA2BhsCVEUBzwGgYJjJawquiNlbsMo5bbJdFvjRImSF7A3mxOZgwdKn9G//bL9Vxi0b1kdJp62c967Q5fts/UO9k3od4VFqfduZuLCq9JqAnaWcqO5nHvnaW8fvUlrJtVsLxnoPAZkJzJTR+bqo4tAruJPCDqkGArOT/T3EM3JrdgcFPOvcbOumKd6Q03uc+b/VRzTGIgPqJLSW/Kzjrbb1tQLNPbdkGHmN1HCQZ2J+zElwmcUxgbjNGweDT2BL0rYCUWalMFAzycxU0rZcslKeeUVLZvlVtqOykVW2ogpbWELjJCKTCgf1/eTxZszeAF/n67AWZvUETm5WH/KWAQV0t36eoLqmOmkrbzdhiGgINga0hkN5zyCaqJiQaAE/29Zl/MxP9xPtLJ0561YF2EH/qusjItetSYBRZWrAh0BCYC4FUc+DjKe4zcK35p57wH6Tj9iBh8lz+8DVdEu281PNW9bgrt3L79qoInqYRlOi1cEPACAJROGgSYxOI9wK8kGDfnmkU8K044nH5/Dpb/RUmP1uLXlIetyWRWxM8veXmTBCvGGf5Vt+kQBNfaGlnCUaEMfiawVi80O8BP6dVKlz7zJl84Utf7tRcEFb5bKwpbi/IKi5OOIg5tIVTue4LS7l0HCv/kwH0KUNHh84eyNpJouMpvyRCU7I8OJlkWV8ojHDAGDx2rDzUD5dK5NmI4EMVMMEjj/Jjj+J3G+IvNsuPYEDzjAuM/IwVbu4OnoApzxhtpS1jL1knXW2ygEuV5ds+lk46T3K4W5JivHdQKY1blNgT6MzJ5PNTLr2FObn8HBnFKz/EYZbEH4KcI+eSEO6pNlGmj4ZVNzdet9FX8PhUdTHWAgWtJuUxLDohz0plWEyY2H1b2Zh3yKN6ttWWodesmzY7Ln0MB+HgwBd4z/VwbIjkDmpdX1ni2YaQfx0if1w91ik4R17hwHbC3IpTYIF9KK60RZqVIFpXhqOwJVwStqjvSjssEhAkaEFrLUZFHRaDZnEJ24peDWEITTpXD/twhEqf4TFWQT4FyN+xY8RM9jx/i99MHbbHorv+Th9wsoSKHiYdq/C6wvtamieemNzZfw5TXNAiUh7Jhy0rtKHkpxQaZfqiwpZxiTckF4XojphVxyHM+JVhrpijtuNGo+2O2Jj0GvGFyh0mG1of2gT7+SsJ/kkv3SCT+EJDQChgKxndtigPH+qxHR28pr4BSyaKWsLFaQ4mUDHGhDoJ9RXNaBFajvhEY2CysYdFY4BvXFMkHtmmIRi4J4NK/VXPIPky9MeUre5gXVYTLeISbA5WMZuTL7Y/oyvanAym7xavaA2o4mgMqOZMJvbuQZtQcH4SPxyDv9GDBnauh1OIMR7RitDgRreu87fwahyobeZwacKh0pfqKFbfdB9cyWUnyvPLMR+2EbY+GCfD9oeLbSb36eKLLQ8CzR251hBVHvrB7HauxvOmcVZwacKh6El1DCsuVvFFbCc8+2wfMpXbDzC0CAQDavysJH4wMrKlKClsKzo8Kj9C+kiu9ZOYsk2Tw5ZxaTaHbjciHE7UaeWZOgMbgxjxrGTZZOxWs5sY8cHkR5ixAmekOPhkdT7KEuYJuHsY4if+JMAQG8rHVuOe3KABreSnb3CXJLiHmkmaWVyacCi6TgOP/TpPRornNyq4RBMHa5ZhpoD4wWjKg+CqTRqER6c9M7CL8AK/Uogde14ij8qDIEZAlwIYgbEYO5Bv15hjFpcmHMa67nc6k4zHIiGw2Ldn/yhIYdRyjjNrQmPX7Sgn+kp8MdnBFOMvk4Q4tAO0MwRJaffhiv++GSXN4tLuOTAiB8gPUAZs2BOjpn+xNkjFD6vtmZ4LPYEQGHFFDpFzueIFDFPNC1yzD6+Ux91n6OERi371P6z15F9EtFVc/gdJK7WZdiYTRQAAAABJRU5ErkJggg==\n", "text/latex": [ "$\\displaystyle y{\\left(t \\right)} = \\frac{t^{2}}{4} - \\frac{t}{8} + \\frac{1}{32} + \\frac{31 e^{- 4 t}}{32}$" ], "text/plain": [ " 2 -4⋅t\n", " t t 1 31⋅ℯ \n", "y(t) = ── - ─ + ── + ────────\n", " 4 8 32 32 " ] }, "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) , phi(t,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": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Multipas 1.9692586606652251\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAqcAAAFKCAYAAAAg1/bFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABCNElEQVR4nO3de5yWc/7H8dens5oiRQo1Tp1ZlLYN20SRQ7QO65CfWLaVSGUTIsecV9IBbTmkWZGwG0mRwsrqoFSUSiepdFJN02Fqvr8/vvcw3Waa+56577nu+5738/G4H1f3dbo/H3fX3cf3+72+lznnEBERERFJBOWCDkBEREREJI+KUxERERFJGCpORURERCRhqDgVERERkYSh4lREREREEoaKUxERERFJGBWCDiCZ1K5d26Wnp8fl3Dt27KBatWpxOXdpSpU8IHVySZU8IHVySZU8IHVyUR4Fmz179kbn3GExO6FIBFScRiE9PZ1Zs2bF5dzTpk0jIyMjLucuTamSB6ROLqmSB6ROLqmSB6ROLsqjYGa2MmYnE4mQuvVFREREJGGoOBURERGRhKHiVEREREQShopTEREREUkYKk5FREREJGGoOBURERGRhKGppERERCRwc+bMObdChQr3OeeOQI1nqSrXzH52zs3bu3fvIy1atFhS0E4qTkVERJJdZib070/bVaugfn0YOBC6dAk6qojNmTPn3MqVKw9NT0/fc9BBB20pV66cCzomiT3nHDk5ORW2bdt2xrp1696bPXt2zxYtWkwK30//ZyIiIpLMMjOhWzdYuRJzDlau9O8zM4OOLGIVKlS4Lz09fU+1atV2qjBNXWZGpUqV9tauXXtrenr63ooVK95R0H4qTkVERJJZ//6Qnb3/uuxsvz5JOOeOOOigg3YFHYeUnmrVqmU7544paJuKUxERkWS2alV06xNTObWYli1mBmAFbdOYUxERkWSUmwuvvAJm4Aqo6+rXL/2YRGJALaciIiLJZt48+OMf4S9/gWOPhSpV9t9etaq/KUokCak4FRERSRZbt0KvXnDqqbB4Mbz4ol+OHAkNGuDMoEEDGDEiqe7WT0XPPvtsLTNrYWYtvv7668rh2999993qedvfeeed6sU59+LFiyvlrevTp0+9//znP785z6WXXpp+5JFHnli8LIKh4lRERCTROefvvm/cGJ591t+Nv3gxXH89lCvnC9EVK5g+dSqsWKHCNIFUq1Ytd9SoUbXC17/yyiu1qlWrlhurzxk0aFDdDz/88DfF6YMPPrj2jTfeWBqrzykNKk5FREQS2cKF0K4dXHMNHH00fPklPPccHHpo0JFJBM4999wtb775Zq3c3F/r0KysLJs0adIhHTt23BLvz2/WrNnu008/fWe8PyeWVJyKiIgkou3boW9fOPlk+PpreOEFmDEDWrYMOrLk8Pzzh1Kv3omUK9eCevVO5PnnA6nmr7322k1r166tNHny5LS8dWPGjKm5b98+u/TSS/crTlu1atWoVatWjcLPceSRR5546aWXphf2GWbWAmDIkCF184YK9OnTpx78tlt/8eLFlcysxWOPPXbYjTfeeNShhx76u4MOOuiUdu3aHZ9/mADAiBEjarZu3bphzZo1f1e1atVTmjRp0nTIkCG/aQV+6KGHDj/22GObValS5dQaNWqc3Lx58yajR48+JPL/SvvT3foiIiKJxDkYNw769IE1a+CGG+Cxx6B27aAjSx7PP38ovXs3YNcu3wi3dm0levduAMBNN20uzVCOO+64PS1btsx65ZVXanXs2DELIDMzs9Y555zzc/Xq1WPSrf/hhx8uat++feNLL710U/fu3TcApKen7znQMc8880zdpk2bZg8fPnzF+vXrKzz88MNHnnvuuQ0XL168sHLlyg7g+++/r9y5c+ctDRs2XFeuXDk3bdq06r17926wc+fOcnfccccGgOeee+7QBx544OhevXr92LZt26zs7Oxy8+bNO2jTpk3FrjFVnIqIiCSKxYvhllvgww/hlFPgzTehdeugowrOX/5yNAsWVI36uHnzqrFnz/5zaO7aVY7bbkvnxRcPi+pczZtn8+KLq6OOIZ+rr75644ABA47Ozs5etWHDhgozZsyoMW7cuAKfK18cZ5999g6AevXq7cn7c1GqVau2b8qUKUvLly8PQJMmTXade+65jYcPH16rd+/eGwEee+yxdXn779u3jwsuuGD7unXrKo4aNeqwvOJ0xowZaQ0bNsx+6qmn1ubte8UVV2wtST7q1hcREQnajh1w991w4okwcyYMGeKXZbkwLYnwwrSo9XHWtWvXLTk5OTZ27NhDRo0adWitWrVyLrroom1BxJKnU6dOW/IKU4BzzjlnR506dXK++OKLannr5s+fX7lTp07HHH744SdVqlSpRaVKlVq8/vrrtVesWPHL3GWnnXbajkWLFlXt2rXr0e+880717du3l7i2VMupiIhIUJyDd97x00OtWgXXXgtPPAF16gQdWWIobotlvXonsnZtpd+sr1t3D19+ubikYUWrZs2aue3bt/95zJgxtX744YdKl1xyyab8hWEQ6tSpkxO+rnbt2jlrQ//dtm7dWq5jx44Nq1Spknvffff90LBhw92VK1d2Q4cOPWzcuHG/jDHp0aPHpl27dtno0aMPGzNmzOEVKlRwbdu23TpkyJDVjRo1OuDQgsKo5VRERCQIy5bBBRfAJZdAjRrwySf+iU8qTEtuwIA1VKmy/3jOKlVyGTBgTUAR0bVr103Tp08/eMmSJQfdcMMNmwrap3Llyrk5OTm/ad3dunVrzBsT169fXzF83caNGyvWrVt3D8DUqVPTfvzxx0rDhw9f2aNHj80dOnTY8cc//jF77969+8VXrlw5+vbtu3H+/Pnfrlu3bu7QoUOXz5s3r9rll19+bHFjU3EqIiJSmnbuhPvug2bN4NNP4emnYc4cOPPMoCNLHTfdtJlBg1ZSt+4ezHyL6aBBK0v7Zqj8OnfuvO3888/ffPXVV29o2bLlroL2Ofroo/esWLGi8q5du34pAN9///20HTt2FFmvVaxY0e3cuTPium7ChAk19+3b98v7yZMnV1u/fn3F1q1b7wDI+8yKFSv+8mzcDRs2lJ8yZcohhZ3zsMMO2/fXv/51S6dOnTYvWbLkoEhjCadufRERkdLy7rvQsycsXw5XXQVPPQX16gUdVWq66abNQRaj4SpUqMCECROWH2ifq666avNrr71W+4orrki//vrrNy5btqzysGHD6qSlpe070HEAxx133K4PP/zw4LfffntbrVq19tavXz8nPT39N133eXbs2FG+Q4cOx3fr1m3DTz/9VOGhhx46skGDBrtvvvnmTQBnnXVWVlpa2r6ePXvWv+eee37Mysoq98QTT9StWbPm3qysrF/GJFx11VUN0tLS9rVp02bHEUcckfPtt99WGT9+fK0zzjij2GNq1XIqIiISbytWwMUXQ6dOULkyfPQR/OtfKkxlP506ddr+xBNPrJw7d261K6+88oQxY8bUfuWVV5bXqFGjyOJ08ODBK6tWrZp75ZVXHt+2bdsmzz777AFnJejVq9faY489dlf37t3T+/XrV79Zs2bZkyZN+i5vGql69ertzczMXLZv3z677rrrjnvggQeOvPbaazdedtll+xX8bdq0yZo3b16122+/vf7FF1/c8B//+EfdSy65ZPPrr79+wEL8QNRyKiIiEi+7d/vW0YEDwQwef9zf/FTpt/fqSGrp2bPnpp49exY4tjTPhRdeuN05Nzv/ur59+27s27fvxvzr1qxZM7+oc59zzjk7Fi5c+G34Z4wfP35FQZ9dqVIlN3LkyB9Gjhz5Q2HxXXTRRdsvuuiib8LXP/300z/m/fnWW2/ddOuttx4wz2glZcupmR1pZi+a2Toz221my83s0QiOe9nMXAGvZ0ohbBERKUsmT/ZTQ91zj7/xadEiuOMOFaYiRUi6llMzSwf+CywHegLrgXTg+AhPsQG4KGzd2oJ2FBERidrq1dC7N4wfDyecAJMmwbnnBh2VSNJIuuIUeB5YA7RzzuUN9J0exfF7nHNfxD4sEREp0/bsgWeegQcfhH374KGHoG9fP8ZUJEE0atRoT/hQgkSTVMWpmR0HnAtcm68wFRERCdbHH0OPHvDtt3DRRTB4MKSnBx2VSFJKtjGnp4eWO81sSmi86RYzG21mtSI8x+FmttHM9prZd2bWz8yCfUyDiIgkpx9/hKuvhrPOgl27YMIE+Pe/VZiKlECyFad5c268CHwHnAf0Ay4APjCzovKZC9wO/Bk/7nQ68CjwQjyCFRGRFJWTA4MGQePG8NZbflL9hQvhwguDjixZ5ebm5gby3HsJhnMOwBW0zUIbA2Fm7YEpEew63TmXYWZ3AwOBCc65X25qMrMrgLHA+c6596OMYRDQC2jonFtSwPZuQDeAOnXqtBg7dmw0p49YVlYWaWlpcTl3aUqVPCB1ckmVPCB1ckmVPCB1cokmj4O//poTBg8m7fvv2fT737P01lvZeeSRcY4wMrH+Ptq1azfbOdcyZicsxNdff/35cccdd3C1atV2xvuzJDFkZWVV/f7779eddNJJZ4VvC3rM6edAkwj2yw4t8+bRCi9oJ4eWpwBRFafAa/jitCXwm+LUOTcCGAHQsmVLl5GREeXpIzNt2jTide7SlCp5QOrkkip5QOrkkip5QOrkElEe69f7qaBGj4b69eGtt6jVuTO1LHEa/JL1+9i7d+8DK1asGJqens5BBx20q1y5csG1nEncOOfIycmpsG3btrR169bl5uTkPFHQfoEWp865bGBRFIcszDu0kO25xQgj71dFF4KIiPzWvn3w3HN+vtLsbLjrLujfH6pVCzqylHHqqad+MGfOnFuWLVt2n3PuCJJv2KFEJtfMtjjnPtm7d++jLVq0+E2jIATfchqtL4B1QEdgaL71HUPLmcU459X4wrQ4x4qISCr74gu4+Wb46ito3x6GDoVGjYKOKiWdeuqpHwAfBB2HBC+pilPn3F4zuxN42cyeB97CT74/EJgGTM3b18w+Aho4544PvW8AvIofm7oUqAz8CbgOeME5t6z0MhERkYS2cSPceSeMGgX16sHrr8Pll/tHkIpIXCVVcQrgnHvFzHLxd+lfD2wGxgB3uf3v7irP/vltD+3bD6iDby39Fv+UqeGlELqIiCS63FwYOdJ33W/bBn//OwwYANWrBx2ZSJmRdMUpgHPuVXwr6IH2yQh7vxnoHL+oREQkqc2a5bvwZ86Etm1h2DBo1izoqETKnKQsTkVERGJm82ZOGDTIT6B/+OEwZoyfWF9d+CKB0N1wIiJSNuXmwksvQaNG1Hv3XejZExYvhi5dVJiKBEjFqYiIlD1z58KZZ8Jf/gINGzLrhRfgmWfg4IODjkykzFNxKiIiZcfWrXDbbdCiBSxZ4ltOP/2UHccfH3RkIhKiMaciIpL6nIPMTH/3/U8/wU03wcCBULNm0JGJSBgVpyIiktoWLIAePeCTT6BVK3jvPd9yKiIJSd36IiKSmrZv9y2lJ5/sC9QRI2DGDBWmIglOLaciIpJanIM33oA+feDHH+HGG+HRR6F27aAjE5EIqOVURERSx6JFcM45cOWVUKeObyn95z9VmIokERWnIiKS/Hbs8I8cPekk/4SnoUP9snXroCMTkSipW19ERJKXc/DOO9CrF6xaBV27whNP+Cc9iUhSUnEqIiLJaelS/1Sn99+HE0+ETz+FM84IOioRKSF164uISHLZuRPuuw+aN4fPPoNBg2DOHBWmIilCLaciIpI83n3Xt5YuXw5XXw1PPgn16gUdlYjEkFpORUQk8S1fDhddBJ06QZUqMHWqf+KTClORlKPiVEREEtfu3fDww9C0qS9In3gC5s6Fdu2CjkxE4kTd+iIikpg++ABuucXf+HTZZX5s6VFHBR2ViMSZWk5FRCSxrF7ti9GOHcHMF6njxqkwFSkjVJyKiEhi2LMHHn8cGjeGiRNh4ECYP98/8UlEygx164uISPCmToUePfzjRzt39l346elBRyUiAVDLqYiIBOfHH+Gqq+Dss33L6bvvwttvqzAVKcNUnIqISOnLyYGnn4ZGjXwxet99sGABXHBB0JGJSMDUrS8iIqXrk098F/6CBXD++fDss3DccUFHJSIJQi2nIiJSOtavh2uvhbZtYft2eOcd342vwlRE8lFxKiIi8bV3Lwwd6rvwx46Fu++Gb76Biy/2U0WJiOSjbn0REYmfGTPg5pv9U506dPBFasOGQUclIgksqVpOzew6M3MHeB0RwTk6m9lXZrbLzFaa2T1mVr404hcRKTM2bIAbboA2bfyf33jDT6avwlREipBsLafvAX8IW2fABOB759y6Ax1sZucC44FRQB/gFOARoDrQL+bRioiUNfv2wciRcNddflxp374wYACkpQUdmYgkiaQqTp1zG4AN+deZ2ZlALeC+CE7xGPCZc65b6P3HZpYG3GNmg4oqbkVE5ABmzfJd+DNnQkaG78Jv1izoqEQkySRVt34hugJ7gLEH2snMjgZOBsaEbXoVqAicF4/gRERS3ubN0L07tGoFq1dDZqZ/4pMKUxEphqRqOQ1nZgcBlwPvOuc2FbF73q/kgvwrnXPLzSwbaBqHEEVEUlduLkdMnAiXXw5btsBtt8EDD0CNGkFHJiJJLKmLU6AzUAN4JYJ9Dw0ttxSwbUu+7SIiUpS5c+Hmm2k8YwacfjoMGwa/+13QUYlICjDnXHAfbtYemBLBrtOdcxkFHD8JOBWo55zbW8RndcF36Td2zi0O27YGmOScu6GA47oB3QDq1KnTYuzYA44eKLasrCzSUuCGgVTJA1Inl1TJA1Inl2TOo3xWFse8+CJH/vvf5NSowTfXXcfPnTpBueQeJZbM30l+sc6jXbt2s51zLWN2QpEIBN1y+jnQJIL9ssNXmFldoD0wpKjCNGRzaFlQC+kh+bbvxzk3AhgB0LJlS5eRkRHBR0Vv2rRpxOvcpSlV8oDUySVV8oDUySUp83AOxozxd99v2ADdu1PpoYf4ed685MulAEn5nRQgVfKQsi3Q4tQ5lw0sKubh1wDliaxLH2BhaNkMmJG30szSgarAN8WMQ0QktS1YAD16wCef+JueJk6EU08NOioRSVHJ3A9zLfC1c25uJDs751YB84AuYZuuAXKA92ManYhIstu+HW6/HU4+2ReoI0b4Jz6pMBWROAq6W79YzOxUoDlw+wH2+Qho4Jw7Pt/qu4F3zewF4DX8JPz3AIM1x6mISIhz8PrrvjBduxZuvBEefRRq1Qo6MhEpA5K15bQrsBfIPMA+5Qkrvp1zE4HLgNbAB0Bv/BOi7oxPmCIiSWbRIujQAa66Co44wreUjhihwlRESk1Stpw6524Dbitin4xC1r8FvBWHsEREkteOHfDww/CPf0C1ajB8OHTrBuXLBx2ZiJQxSVmciohIjDgHb78NvXr5pztddx08/jgcfnjQkYlIGZWs3foiIlJSS5bAeefBpZdCzZrw6afw0ksqTEUkUCpORUTKmp07YcAAaN4cPv8cnnkGZs+GM84IOjIREXXri4iUKRMmQM+esGIFdOkCTz4JdesGHZWIyC/UcioiUhYsXw4XXeRfVavCxx/7Jz6pMBWRBKPiVEQkle3aBQ89BE2bwtSpvqV07lzQIy5FJEGpW19EJFV98AHccgssXQp//rOfJuqoo4KOSkTkgNRyKiKS7DIzIT0dypXzy2ef9Xfgd+zo102e7J/4pMJURJKAWk5FRJJZZqafLD87279fuRJuuw0qVoRHHoE+faBy5WBjFBGJgopTEZFk1r//r4VpfocfDnfdVfrxiIiUkLr1RUSS2apVBa//8cfSjUNEJEbUcioikoz27oVhwwrfXr9+6cUiIhJDajkVEUk2M2ZAy5bQq5d/ylOVKvtvr1oVBg4MJDQRkZJScSoikiw2bYK//hXatIGNG+HNN2HePBg5Eho0ADO/HDHCP/1JRCQJqVtfRCTR5ebCiy/CnXfC1q3Qty8MGABpaX57ly4qRkUkZag4FRFJZHPnQvfu8MUXcOaZMHy478oXEUlR6tYXEUlEW7f6+UpbtIDvv4fRo2H6dBWmIpLy1HIqIpJInIOxY/3k+evX+1bThx+GmjWDjkxEpFSoOBURSRSLFkGPHjB1qr8bf8IEvxQRKUPUrS8iErTsbLj7bjjpJJgzB557zo8xVWEqImWQWk5FRIL0n/9Az56wciV07QpPPOEfPSoiUkap5VREJAjLl0OnTnDxxVC9OnzyCbz8sgpTESnz1HIqIlKadu+m/pgxkJkJ5cvDU0/5ltOKFYOOTEQkIag4FREpLVOmwC23cOx338Fll8GgQXDUUUFHJSKSUNStLyISbz/+CFdeCeecA7m5zHv8cRg3ToWpiEgBVJyKiMTL3r2+dbRxY3jnHXjgAZg/ny2tWgUdmYhIwkqq4tTMrjMzd4DXEUUc/3Ihxz1TSimISFnx3//6pzv16eMfO7pwIQwYAFWqBB2ZiEhCS7Yxp+8BfwhbZ8AE4Hvn3LoIzrEBuChs3doYxCYiAhs2QL9+8NJLcPTR8NZb0LkzmAUdmYhIUkiq4tQ5twFfXP7CzM4EagH3RXiaPc65L2Idm4iUcbm5MHIk3HknbN/uC9R774Vq1YKOTEQkqSRVcVqIrsAeYGzQgYhIGTVnDnTvDl9+CW3bwvDh0LRp0FGJiCSlpBpzGs7MDgIuB951zm2K8LDDzWyjme01s+/MrJ+ZlY9jmCKSqn7+GW69FU47zT/hacwY+PhjFaYiIiWQ7C2nnYEawCsR7j8XmA0sBKoAfwIeBU4Abox9eCKSkpyDf/0Lbr/djzG9+WZ46CE45JCgIxMRSXrmnAvuw83aA1Mi2HW6cy6jgOMnAacC9Zxze4sZwyCgF9DQObekgO3dgG4AderUaTF2bHxGD2RlZZGWlhaXc5emVMkDUieXVMkDEiOXqitWcMLgwdScO5dtjRvzXe/eZDVsGNU5EiGPWEmVXJRHwdq1azfbOdcyZicUiUDQxWlVoH4Eu2Y751aFHVsXWA0Mcc71LkEMrYD/AVc751470L4tW7Z0s2bNKu5HHdC0adPIyMiIy7lLU6rkAamTS6rkAQHnsmOHbx39xz+genV47DG48UYoF/3oKH0niUd5FMzMVJxKqQu0W985lw0sKubh1wDlibxLvzB587sEV6WLSOJyzk+gf9ttsHo1XH89PP44HHZY0JGJiKSkZL4h6lrga+fc3BKe52p8YTqzxBGJSGr5/nu48EK45BI/nvSzz+DFF1WYiojEUVLeEGVmpwLNgdsPsM9HQAPn3PGh9w2AV/FTTi0FKuNviLoOeME5tyzOYYtIsti1C558Eh55BCpUgKef9nflV0jKn0wRkaSSrL+0XYG9QOYB9inP/vltBzYD/YA6+NbSb4GewPD4hCkiSWfyZOjRA5YuhSuu8GNMjzwy6KhERMqMpCxOnXO3AbcVsU9G2PvN+KmnRER+64cfoE8fGDcOTjjBF6kdOgQdlYhImZPMY05FREouJ8e3jjZuDBMm+Dvy589XYSoiEpCoWk7NrDXQEWgN1AMOAjYCi4HpwDvOuS2xDlJEJC4+/dRPoL9gAVxwAQwZAsccE3RUIiJlWkQtp2bW1czmA5/jJ6yvCizBzw+6Bfg9MBJYY2Yvm5l+3UUkcf30E1x3Hfzxj7Btm58qasIEFaYiIgmgyJZTM5sHHA6Mxk/fNNcVMHO/mR0MXAh0ARaa2fXOuddjHK+ISPHt2wf//CfcdZefVP+uu6B/f6hWLejIREQkJJJu/ZeA551zuw60k3NuK/7u+Uwz+x1wRAziExGJjdmzoXt3mDkT2rWDYcOgSZOgoxIRkTBFdus7554pqjAt4Jh5zrkPih+WiEiMbNnip4Y67TT/hKfMTPjoIxWmIiIJqlh365t3pJmlxTogEZGYcA5Gj/Z34T//vJ9Ef9EiuPpqMCv6eBERCURUxamZVTGzYUA2sAr42cxOjkdgIiLFtnAhZGRA165w7LEwaxYMHgwHHxx0ZCIiUoRoW04fwd/w1B9/81O5vHOY2dtm1iO24YmIRCErC+64A04+2U8P9c9/wn//C6ecEnRkIiISoWiL08uBu5xzTwOTw7ZNAS6LSVQiItFwDsaP9+NIn3zSt5guXgw33gjl9KwREZFkEu2v9qHAd4VsWw40Llk4IiJRWroUzj8fLrsMatXyLaUjR0Lt2kFHJiIixRBtcboAyChkWzagAV0iUjp27YL774fmzX1B+swzfmxpmzZBRyYiIiUQ1eNLgRHAEDObBbwbtu13wI8xiUpE5EDef9/ffb9sGVx1FTz1FNSrF3RUIiISA1EVp865UWbWBngL+ARwQDMzOwG4Gxgb+xBFRLzKP/3ku+/Hj4dGjeDDD+Hss4MOS0REYijallOcczeY2XT8HfsGvBLa9DHwQAxjExHxcnLgmWdoNWCAn6N04EC4/XaoXDnoyEREJMaiLk4BnHOjgdFmlg4cCaxxzq2IYVwiIt4nn8DNN8PChWxp04bamZmQnh50VCIiEifFKk7zhArSFTGJREQkv/XroW9fePVVaNAA/v1vFtSoQYYKUxGRlFbk3fpm9m8zi3gG69BTpPqY2U0lC01EyqR9+2D4cD+mdOxY6N8fvvkGLroo6MhERKQURNJyugr4wszmApnAZ8DXzrm9eTuYWT2gFdAJuARYA/wl5tGKSGqbORO6d4fZs/2NTsOG+SJVRETKjCJbTp1ztwJNgS+B+4GZwC4z22xma81sF7Aafwd/M6AXcJJz7st4BS0iKWbLFl+U/v738OOP8NprMGWKClMRkTIoojGnzrllwK1mdjvwB+D3QD2gCrAJWAR84pxbGa9ARSRFZGb6rvpVq+Doo6FjR3j7bdi0CW67DR54AGrUCDpKEREJSLTznO4BpodeIiLRycyEbt0gO9u/X7UKRoyA44+HyZPh5JMDDU9ERIIX7eNLRUSKr3//XwvT/PbsUWEqIiJAMaaSMrOuwFVAfXy3fn7OOXdcLAITkRTjnG8pLcjq1aUbi4iIJKyoilMzuxf/FKgFwFxgdxxiEpFUs3Qp3HqrL1ALUr9+6cYjIiIJK9qW0xuAwc653vEIRkRSzK5d8Pjj8OijUKkSXHMNvPXW/l37Vav6x5GKiIgQ/ZjTWsCEeAQSKTOrZWaDzex7M9tpZsvNbKiZHRbh8Z3N7Csz22VmK83sHjMrH++4RcqcyZPhxBPh/vuhc2dYtMg/7WnECP/EJzO/HDECunQJOloREUkQ0bacTgd+B0yNQyxFMjMD/gM0BAYA3+LnYH0IaGFmbZwrrN8QzOxcYDwwCugDnAI8AlQH+sU3epEyYs0a6NMH3ngDTjjBF6kdOvy6vUsXFaMiIlKoaIvTXsBbZrYJmAhsDt/BOZcbg7gKcwLQBvibc25EaN00M8sFnsMXrYsPcPxjwGfOuW6h9x+bWRpwj5kNcs6ti1fgIilv714YMgQGDICcHHjwQejbF6qE3zcpIiJSuGi79b8DmgMvAeuBnLDXnphG91uVQsttYet/Di0LzcfMjgZOBsaEbXoVqAicV/LwRMqozz+HFi18i+mZZ8LChXDvvSpMRUQkatG2nD4IFNptXgoWAp8A95rZUvyTqZriu/jfd859e4Bjm4WWC/KvdM4tN7Ps0HlEJBqbNkG/fjBqFBx1FIwfD3/6kx9PKiIiUgzRPiHq/jjFEennOzM7H9/aOTPfpveAy4s4/NDQcksB27bk2y4iRcnNhZde8oXp1q2++37AAEhLCzoyERFJcnaA+4fi/+Fm7YEpEew63TmXETrmX0AGfr7Vb4EmoT/PBjoVNubVzLrgu/QbO+cWh21bA0xyzt1QwHHdgG4AderUaTF27NiIcotWVlYWaSnwD3uq5AGpk0us86i2bBkNBw3i4IUL+fnEE1nSuzc7jjkmZuc/EH0niSdVclEeBWvXrt1s51zLmJ1QJBLOuQO+gFxgX4SvvUWdL+zcVYHGEbzqh/a/AD+s4Oyw83QIrb/4AJ91XmifPxSwbQfwZFHxtmjRwsXLxx9/HLdzl6ZUycO51MklZnls2+Zc797OlS/vXO3azr30knO5ubE5d4T0nSSeVMlFeRQMmOWi+HddL71i8YqkWz9u40ydc9n4caOROjG0nBm2/svQsgnw70KOXRhaNgNm5K00s3R8kfxNFHGIlB3OwZtvQq9esHYtdOsGjzwCh2okjIiIxF6RxakLeJxpmLypnloBH+Zb//vQck1hBzrnVpnZPKALMDLfpmvwMw28H8M4RVLDkiVwyy1+rtKTT/Y3PLVuHXRUIiKSwqKdSipobwE/AqPNrLuZtTOz7sBoYDXwdt6OZvZR6I7+/O4G2prZC2aWYWa9gXvwj2TVHKcieXbt8k92OvFEmDEDBg+GmTNVmIqISNxFO5VUoJxz28ysNXA/cAdQF1iLf6Tq/c65rHy7lycsP+fcRDO7DLgPuA4/V+sjgB7sLZJn0iTfWrpsGVx5JTz9NNStG3RUIiJSRiRVcQrgnFsN/Oau+gL2yyhk/Vv4FlgRye+HH6B3bz++tGFDmDIF2rcPOioRESljkq1bX0RiLSfHt442aQLvvgsPPwxff63CVEREApF0LaciEkP//S907w7z58P558OQIXDssUFHJSIiZZhaTkXKoo0b4YYb4Iwz4Oef4e23faupClMREQmYilORsiQ3F0aOhEaNYPRouOMO+OYb6NwZzIKOTkRERN36ImXG3Lm+C/+LL+DMM2H4cGjePOioRERE9qOWU5EUV37HDv90pxYt/PRQr7wC06erMBURkYSkllORVOUcvPEGrXr0gM2b4W9/g4ED9dhRERFJaCpORVLRd9/5ifSnTGHPCSdQeeJEaNUq6KhERESKpOJUJJXs3AmPPgqPPw5VqsCQIcxu0oQMFaYiIpIkNOZUJFW8/74fR/rQQ3DZZbB4sW89LV8+6MhEREQipuJUJNmtXg2XXuon0a9YET76CDIz4Ygjgo5MREQkaipORZJVTg489ZR/7OjEif5mp3nz4Kyzgo5MRESk2DTmVCQZffaZn7N0wQK48EJ49lk45pigoxIRESkxtZyKJJMNG+D66/0k+tu2wTvvwH/+o8JURERShopTkWSQmwsjRvjHjo4ZA3fe6R87evHFeuyoiIikFHXriyS6r77yXfj/+x+0besfO9q0adBRiYiIxIVaTkUS1dat0LMntGwJy5fD6NHw8ccqTEVEJKWp5VQk0TgHY8dCnz6wfr1vNX34YahZM+jIRERE4k7FqUgiWbwYevTwc5W2aOFvdjrttKCjEhERKTXq1hcJQmYmpKdDuXJ++dJLcO+9cNJJMGsWDB3qx5iqMBURkTJGLacipS0zE7p1g+xs/37lSrjhBt+df8018OSTerqTiIiUWSpORUpb//6/FqZ5nIPDD4dXXw0mJhERkQShbn2R0rZqVcHrN2wo3ThEREQSkIpTkdL06adQoZAOi/r1SzcWERGRBKTiVKQ0bNgA110Hf/wjVK8OlSvvv71qVRg4MJDQREREEomKU5F4yv/Y0cxM6NfPd+uPGgUNGvhHjzZo4Pfp0iXoaEVERAKXdMWpmdUys8Fm9r2Z7TSz5WY21MwOi+DYl83MFfB6phRCl7Jm7lw4/XT429/gxBP9+8ceg2rVfCG6YoUvXlesUGEqIiISklR365uZAf8BGgIDgG+BpsBDQAsza+Occ0WcZgNwUdi6tbGOVcqwbdtgwAAYMgRq1YJXXoH/+z/fSioiIiIHlFTFKXAC0Ab4m3NuRGjdNDPLBZ7DF62LizjHHufcF3GMUcoq52DcOOjdG9au9S2mjzyix46KiIhEIdm69SuFltvC1v8cWiZbPpIqli6Fjh3hiiugTh2YMQOee06FqYiISJSSrZhbCHwC3GtmLc0szcxa4bv433fOfRvBOQ43s41mttfMvjOzfmZWPq5RS+ratQvuvx+aN/cF6eDB8OWX8PvfBx2ZiIhIUkqqbn3nnDOz84FXgZn5Nr0HXB7BKeYCs/FFbhXgT8Cj+OECN8Y0WEl9kydDjx6+1fTKK+Hpp6Fu3aCjEhERSWpW9P1Dcfxws/bAlAh2ne6cywgd8y8gA3gAf0NUk9CfZwOdnHO5UcYwCOgFNHTOLSlgezegG0CdOnVajB07NprTRywrK4u0tLS4nLs0pUoeUHgulTZs4Pjhwzl82jSyjzqKJbfdxpaWLQOIMDJl4TtJNqmSB6ROLsqjYO3atZvtnEvcHzhJTc65wF5AVaBxBK/6of0vABxwdth5OoTWX1yMGFqFjr2qqH1btGjh4uXjjz+O27lLU6rk4VwBueTkODdokHNpac5Vruzcgw86t3NnEKFFJaW/kySVKnk4lzq5KI+CAbNcgHWCXmXzFWi3vnMuG1gUxSEnhpYzw9Z/GVo2Af4dZRh58/sE14QsiW/GDOjeHebN8zc+DR0Kxx0XdFQiIiIpJ9luiFoXWrYKW59398maYpzzanxhGl7wisDmzdCtG7RpAxs3wptvwsSJKkxFRETiJKluiALeAgYCo83sIXyra2PgPmA18Hbejmb2EdDAOXd86H0D/I1UY4GlQGX8DVHXAS8455aVXhqS8JzjiEmT4PLLYcsW6NPH35VfvXrQkYmIiKS0pCpOnXPbzKw1cD9wB1AX/3SnCcD9zrmsfLuXZ//8tgObgX5AHXxr6bdAT2B43IOX5LFgAXTvTuPPPoM//AGefx5OOinoqERERMqEpCpOAZxzq4EbItgvI+z9ZqBzfKKSlJCVBQ8+CIMGQY0aLPr732n8+ONQLtlGv4iIiCQv/asr4hy8/TY0bQpPPgnXXguLF7PuggtUmIqIiJQy/csrZdvy5dCpE1xyCRxyCHz2GYwaBbVrBx2ZiIhImaTiVMqm3bvhkUegWTOYNg3+8Q+YPRtOPz3oyERERMq0pBtzKlJiU6fCzTfD4sVw6aXwzDNw1FFBRyUiIiKo5VTKknXr4Jpr4OyzISfHz1f65psqTEVERBKIilNJffv2wbBh0LgxjBsH997rp4s677ygIxMREZEw6taX1DZrln/s6KxZvsV02DBo1CjoqERERKQQajmV1PTzz3DLLdCqFfzwA7z2GkyZosJUREQkwak4ldTiHGRm+i78557zBeqiRXDllWAWdHQiIiJSBHXrS+pYtAh69PB347dq5W94OvXUoKMSERGRKKjlVJJfdjb07w8nnQRz5vgW088/V2EqIiKShFScSvLIzIT0dP9I0fR0//699/xE+o884rvuFy2Cm26C8uWDjlZERESKQd36khwyM6FbN99KCrByJXTt6qeJatLEP+WpbdtAQxQREZGSU3EqyaF//18L0zz79sEhh8DcuVCpUhBRiYiISIypW1+Sw6pVBa/fulWFqYiISApRcSqJb+NGqFq14G3165duLCIiIhJXKk4lceXmwsiRfuL87GyoEDYKpWpVGDgwmNhEREQkLlScSmL6+ms44wz461+heXOYPx9efhkaNPCT6TdoACNGQJcuQUcqIiIiMaQboiSxbN8O998PgwdDzZq+IL32Wl+QNmumYlRERCTFqTiVxOAcvPUW3HYbrFnjp4169FE49NCgIxMREZFSpG59Cd7338MFF8Bll0Ht2jBjBrzwggpTERGRMkjFqQRn9254+GHfXf/ppzBoEMyaBa1bBx2ZiIiIBETd+hKMjz6CHj1g8WK4/HJfmB55ZNBRiYiISMDUciqla906f1NT+/awdy+8/z688YYKUxEREQFUnEpp2bcPhg2Dxo3hzTdhwAA/PVTHjkFHJiIiIglE3foSf7NmwU03wezZvsV02DBo2DDoqERERCQBJV3LqZnVNrMXzWyDme00s/+Z2blRHN/ZzL4ys11mttLM7jGz8vGMucz6+We45RZo1cpPD/XaazB5sgpTERERKVRSFadmVhmYCnQE7gAuAVYD75pZRgTHnwuMB2YC5wGDgXuAR+ITcRnlHGRm+i78557zBeqiRXDllX4yfREREZFCJFu3/uXAiUA759w0ADObBMwDngBaFXH8Y8Bnzrluofcfm1kacI+ZDXLOrYtP2GVH1VWrfNf91Klw2mkwcSKcemrQYYmIiEiSSKqWU6A1sBOYnrfCOeeAycBpZlboLd9mdjRwMjAmbNOrQEV8S6oU186dcM89tLzhBj+2dPhwP5m+ClMRERGJQrK1nO4DckIFaX67Q8vmwJpCjm0WWi7Iv9I5t9zMsoGmMYuyrJk40XfdL1/OTx06cMSrr0KdOkFHJSIiIkko2YrTxUANM2vinPs23/o/hJYHet5l3rYtBWzbUsSxUpAffoBevWD8eD++dOpUFplxhApTERERKSb7bSNkKX64WXtgSgS7TnfOZZjZIcB3wArgBmAt0A14ECgPXOmce72Qz+qC79Jv7JxbHLZtDTDJOXdDAcd1C30GderUaTF27NjIkotSVlYWaWlpcTl3rNm+fRw5fjzHvPQS5Oay8tprWf3nP+MqVkyqPIqSKrmkSh6QOrmkSh6QOrkoj4K1a9dutnOuZcxOKBKBoFtOPweaRLBfNoBz7mczuxR4Bfg6tG0ZcD/wEL5YLczm0LKgFtJD8m3fj3NuBDACoGXLli4jIyOCcKM3bdo04nXumPr8c+jeHb7+Gi64AIYM4dhjjuHY0OakySMCqZJLquQBqZNLquQBqZOL8hBJHIEWp865bGBRlMd8ambHAcfjW0u/A/rib5Sac4BDF4aWzYAZeSvNLB2oCnwTTRxlzqZN0K8fjBoFRx0Fb78NF1+sqaFEREQkppLtbn3A36HvnFvinFuELyz/CrzqnMs6wDGr8FNOdQnbdA2QA7wfr3iTWm4uvPgiNGoEr7wCffvCt99C584qTEVERCTmgu7Wj5qZPQrMBjbiW0/74ovLu8L2+who4Jw7Pt/qu/ET9r8AvAacgp+Ef7DmOC3A/Pm+C/+//4UzzvAT6jdvHnRUIiIiksKSseW0DvAMfm7T+0PL051z4WNGyxNWfDvnJgKX4edL/QDojX861J1xjTjRZWZCejqUK+eXL77oW0hPOcU/2WnUKJg+XYWpiIiIxF3StZw65/4S4X4Zhax/C3grljEltcxM6NYNsrP9+5Ur4cYb/SNIb7wRHnsMatUKNkYREREpM5KuOJUY69//18I0j3N+Ev1//jOYmERERKTMSsZufYmlVasKXv/TT6Ubh4iIiAgqTsu2adOgfPmCt9WvX6qhiIiIiICK07Lpp5/g2muhXTs45BCoXHn/7VWrwsCBgYQmIiIiZZuK07IkNxdeeMHPWTp2rB9vunKlvxu/QQM/b2mDBjBiBHQJnw5WREREJP50Q1RZ8dVXfs7S//0PMjL8nKWNG/ttXbqoGBUREZGEoJbTVLdtG/TqBS1bwvLl8OqrMHXqr4WpiIiISAJRy2mqcg7GjYPevWHtWrjpJj+OtGbNoCMTERERKZRaTlPR0qVw3nlwxRV+vtIvvoDhw1WYioiISMJTcZpKdu+GBx/0jxn9/HMYPBi+/BJatQo6MhEREZGIqFs/VXz4Idx8MyxZ4ltMn34a6tULOioRERGRqKjlNNmtWwdXXw0dOvhxph984KeJUmEqIiIiSUjFabLatw+GDvVzlo4fD/fdB/PnwznnBB2ZiIiISLGpWz8ZzZrl776fPdu3mA4bBiecEHRUIiIiIiWmltNk8vPPcMst/ganH3/03fcffKDCVERERFKGitNk4Bz8619+4vznnoNbb4Vvv/U3PpkFHZ2IiIhIzKhbP9EtXgw9esBHH8Fpp8HEiXDqqUFHJSIiIhIXajlNVDt3woABcNJJfozpsGEwY4YKUxEREUlpajlNRJMm+dbS77+HLl3gqafgiCOCjkpEREQk7tRymkjWrIHLL/ePHq1Y0XfljxmjwlRERETKDBWnQcjMhPR0KFfOL199laPGjfM3PL37Ljz8MMybB2edFXSkIiIiIqVK3fqlLTMTunWD7Gz/fuVK6NqV453zLaZDh8KxxwYbo4iIiEhAVJyWtv79fy1M8zjHnoMPptJ772lqKBERESnT1K1f2latKnB1xW3bVJiKiIhImafitLTVr1/g6t2HH17KgYiIiIgkHhWnpW3gQKhadf91Vavy/Y03BhOPiIiISAJJuuLUzGqb2YtmtsHMdprZ/8zs3AiPfdnMXAGvZ+Ic9q+6dIERI6BBA9+N36ABjBjBT+3bl1oIIiIiIokqqW6IMrPKwFSgNnAHsA64AXjXzDo456ZFcJoNwEVh69bGMs4ideniX/lNm1aqIYiIiIgkoqQqToHLgROBdnmFqJlNAuYBTwCtIjjHHufcF3GLUERERESKLdm69VsDO4HpeSuccw6YDJxmZkcGFZiIiIiIlFyyFaf7gJxQQZrf7tCyeQTnONzMNprZXjP7zsz6mVn52IYpIiIiIsWRbN36i4EaZtbEOfdtvvV/CC0PLeL4ucBsYCFQBfgT8ChwAqDb5UVEREQCZr9thCzFDzdrD0yJYNfpzrkMMzsE+A5Ygb8Rai3QDXgQKA9c6Zx7PcoYBgG9gIbOuSUFbO8W+gzq1KnTYuzYsdGcPmJZWVmkpaXF5dylKVXygNTJJVXygNTJJVXygNTJRXkUrF27drOdcy1jdkKRCARdnFYFCp6Vfn/ZzrlVoWPOBF4BjgltWwa8DDwEtHXOfRJlDK2A/wFXO+deO9C+LVu2dLNmzYrm9BGbNm0aGRkZcTl3aUqVPCB1ckmVPCB1ckmVPCB1clEeBTMzFadS6gLt1nfOZQOLojzmUzM7Djge31r6HdAXf6PUnGKEkffM0OCqdBEREREBkm/MKfDLHfpLAMwsDfgr8KpzLqsYp7saX5jOjF2EIiIiIlIcgXbrF4eZPYq/qWkjvvW0L5ALnO6c25xvv4+ABs6540PvGwCvAmOBpUBl/A1R1wEvOOe6R/DZG4CVscwnn9r4nJJdquQBqZNLquQBqZNLquQBqZOL8ihYA+fcYTE8n0iRkrHltA7wDHA48BPwNnBf/sI0pDz757cd2Az0C53DAd8CPYHhkXxwPC9QM5uVCuN6UiUPSJ1cUiUPSJ1cUiUPSJ1clIdI4ki64tQ595cI98sIe78Z6ByHkEREREQkRpJtEn4RERERSWEqThPHiKADiJFUyQNSJ5dUyQNSJ5dUyQNSJxflIZIgku6GKBERERFJXWo5FREREZGEoeI0AGa2wsxcAa/OER7f2cy+MrNdZrbSzO4xs/JxDjuSuK4K5fFDhPuXN7PeZrbAzHaY2Voze9vMTop3rEXEFVUeoWMOMrP7zWyJme02s/Vm9q6ZVYpnrBHEFXUu+Y491syyQ8cfH4/4oogl4jzMrIaZDTCzz81sk5n9HPpz51IItajYivN3KyGudzOrbmZvmNnS0PX6s5n9z8yuifD4hLjeS5pH6BwJcb3HIpd850qY610k6e7WTyEfAPeHrVtc1EFmdi4wHhgF9AFOAR4BquOnyQqEmR0CDALWRXHYQ/iYHwWm4ufnuwf42Mx+55yLuqAqqeLkYWYVgffxj9R9FPgGOAzogJ/SLBDF/E7yGw5sBQ6KVUzFUYw86gM3Ay/h/47lAlcBb5vZLc65YfGIsyjF/LuVSNd7JWAv/u/4Cvxc0VcAr5rZYc65QUUcnyjXe4nySLDrvaTfSX4Jcb2LgMacBsLMVgCfOeeK83+3XwHbnHNt860bgP+Rr++cK24hUiJmNgJoAKwF2jvnjorgmB+Bac65q/Ota4yff/Ym59wL8Yr3ADEVJ487gbuBZs651XEOMWLFySXfsVfjC6lHQ8sTnHNL4xJo0bFElYeZVcM/SC47bP1H+Dzqxy3YA8dVnL9bCXm952dmM4A059yJReyXcNd7flHkkZDXe36R5pJv/4S53kVA3fpJxcyOBk4GxoRtehWoCJxX2jEBmNnpwDVAjygPrQRsC1v3c2hZ6n83S5DHzcC4RPqHqgS5YGY1gaeBv/Pr9xGI4uThnNsRXpiGzALqxSq2aBQnj0S93guwCciJYL+Eut4LEGkeCXe9FyDSXBLqehfJkwg/CGVVp9D4nt1m9kWE4+GahZYL8q90zi0HsoGmMY6xSKEurhHAk8X4P+3hwDVmdnFonOCxoXU/AK/HONQDKm4eZlYfOBr43sz+aWbbQmMDPzKzk+MUblExleQ7AXgCWOScezW2kUUnBnmE+yOwKAbniUoJ8ki46x3AvApmVsvMugHn4p/aV5SEud6heHkk4vUeiqu43wkkyPUukp/GnAZjAjATWI5/lOot+PFw/+ecC28lye/Q0HJLAdu25Ntemvrhxzk9Gu2BzrkBZrYbeItf/0fpOyCjgMfRxltx88hrieuH/06vDJ3nAWCamZ3knFsVsygjU+zvxMzOAK7Fj20MWrHzCBf6B7s1vvWytBU3j0S83sG3/g4J/TkHuM05N7qogxLseofi5ZGI1zsU8ztJsOtd5FfOOb1K8ALaAy6C17QDnKM8/odudRGf1SV0rkYFbFsDjCrNXIDjgZ1Ax3zneBn4IcLP6x46/gEgA7gM3/W6DKiXDHkAbULnWgdUzbf+aGA38HiyfCf4bteFwGP51l0XOv/xyZJHAZ+dAewCRpckhwC+j4S63vMddxjQEuiIb/ncB/wtgs9LiOu9JHmQYNd7CXOJ2/Wul14lfanltOQ+B5pEsF9BY+AAcM7tM7NxwONmVtc5t7aQXfNaFwpqMTkk3/biijaXZ/F33X4RuhMZ/A+ehd7vds7tLOgEZnYofuD9k865+/Ktn4q/67Qv0Dv6FIBSzAM/tgvgvy7fOEfn3GozW0TJWyRKM5de+L9bz+Y7tmpoWd3MqjvntkcV/a9KM49fmNlpwH9C57ohypgLUpp5JNr1DoBzbgOwIfR2kplVBZ4ysxedcwWOc0yw6x0oXh4k3vWe9/nFyaUX8bveRUpExWkJhX6gYjGOzfJOeYB9FoaWzYAZvxxolo7/UfmmJAEUI5em+LuPtxSwbQswGP8DWJCG+O6wmWExbDazZUT2A12gUs7je3xrUEHfm+GnMSq2Us6lKXAEvlUu3BxgHv4GnaiVch4AmNmJ+Cnb5gKXHuAf6YiVch6Jdr0XZhbQFT9EqbDpoBLpei9MJHkk2vVemEhyidv1LlJSKk4TgJlVAC4HVrkDTA3jnFtlZvPw3X0j8226Bj/O6P24BvpbVwJVwtbdCbTA53OgeQvz8myFb9kCfmlhOR7/41haip2Hcy7HzN4D/mhm1ZxzO+CXGycaAf+OT8iFKsl38hi+yzm/jvjxddcQwTy8MVSSPDCzE4Ap+GLiwkhaWeOkJH+3Eu16L0xbIAv46QD7JNL1Xpgi80jA670wkXwniXS9i+wv6HEFZe2Fnwx8LH4Qejv8P16f4v9P/MqwfT8CloatOx//f+cv4Mdt9caPp3sy6NxC8b1MAePpCsllAn6c1oPA2cCf8f/HvwdomUR5NMX/QzAN6IQvOhYA64E6yfSdFLDPdSTIGLRI8wAOx3cVbwYuwN8Ilf9VORnyCK1LmOsd+Bv+wQZd8MXPJaHfMgf0iyCXhLjeY5BHwlzvJc2lgPMlzPWuV9l+qeW09C3H/+P5JH68Tza+q6ujc+6DsH3LE9a67ZybaGaXAffhf0jW458YMzC+YZfYb3LBP8nkdnzBfjt+DsQ5wBnOuVmlG17ECvpOvjGzs4DH8VPi5AAfA52dc+tLP8SIFfSdJKPwPPK60gHeLWD/Y/DFa6JJ9Ot9PnAx8BT+t2sjfgL9C51z74Xtm8jXe4nySLDrvaTfiUhC0hOiRERERCRhaBJ+EREREUkYKk5FREREJGGoOBURERGRhKHiVEREREQShopTEREREUkYKk5FREREJGGoOBURERGRhKHiVEREREQShopTEYkrM7vfzIr9tA8zG2JmE8LWNTUzZ2YdDnBcbzP72sz0OycikkT0oy0iCcvMjsM/P/yBsE2nhpYHeuzl8/hHBXeNQ2giIhInKk5FJJH1AuYV8Oz1FsAy59yWwg50zu0ERgN/j194IiISaypORaRU5XXzm9kJZvaemWWZ2UozG5C/C97MKgPXAP8q4DQtgJlm9n9mNsfMdprZN2bWLmy/sUBTM2sTv4xERCSWVJyKSFDeBqYCnYF38F33+bvgWwOHAJ/mP8jMDDgZaAN0AR4G/oz/PRsd9hlzgW1Ax5hGLiIicVMh6ABEpMz6h3PupdCfPzSzs4CrgLx1rQEHfB12XEOgOjDFOXdp3kozOxoYZmYHhbr0cc7lmtnXoXOJiEgSUMupiATlvbD3C4D6+d7XA7Y55/aE7dcitLw7bH3t0P47w9ZvCJ1LRESSgIpTEQnK5rD3u4Eq+d5XCa0Ldyqwwjm3OGz9Kfy2lRVgJ3BQcYMUEZHSpeJURBLVJqBmAetbAHMKWH9KIesPBTbGMC4REYkjFacikqgWARXN7Ki8FaGboU4Bvsq/o5nVBBqErw85BghvZRURkQSl4lREEtUnoWWrfOuOAw7mty2kp4SW+603s0PwN1B9goiIJAUVpyKSkJxzK4AvgU75VufdDFVQcbob+CZs/QXAHvy0VSIikgTMuWI/8lpEJK7M7DpgMFDXOZddjOPfBzY65/4v1rGJiEh8qDgVkYRlZuWB+cCLzrmnojz2ZOALoLlzbmkcwhMRkThQt76IJCzn3D7gL0DUrabAEcD1KkxFRJKLWk5FREREJGGo5VREREREEoaKUxERERFJGCpORURERCRhqDgVERERkYSh4lREREREEoaKUxERERFJGCpORURERCRh/D/RSnmI3hgTUQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "from scipy.optimize import fsolve\n", "\n", "def multipas(phi, tt, sol_exacte):\n", " h = tt[1] - tt[0]\n", " uu = [ sol_exacte(tt[i]) for i in range(2) ]\n", " for i in range(1,len(tt) - 1):\n", " temp = fsolve(lambda x: -x + (1-gamma)*uu[i]+gamma*uu[i-1] + h/4 * ( (gamma+3)*phi(tt[i+1],x)+(3*gamma+1)*phi(tt[i-1],uu[i-1]) ), uu[i])\n", " uu.append(temp[0])\n", " return uu\n", "\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(tt)\n", " uu_mp = multipas(phi,tt,sol_exacte)\n", " H.append(h)\n", " err_mp.append(max([abs(uu_mp[i] - yy[i]) for i in range(len(yy))]))\n", "\n", "print (f'Multipas {(polyfit(log(H),log(err_mp), 1)[0])}')\n", "\n", "figure(figsize=(8,5))\n", "plot(log(H), log(err_mp), 'r-o', label='Multipas')\n", "xlabel('$\\ln(h)$')\n", "ylabel('$\\ln(e)$')\n", "legend(bbox_to_anchor=(1.04, 1), loc='upper left')\n", "grid(True)\n", "show()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercice : convergence\n", "\n", "On notera $\\varphi_k\\stackrel{\\text{déf}}{=}\\varphi(t_k,u_k)$.\n", "Soit la méthode multipas\n", "$$\n", "u_{n+1}=\n", "\\left(\\frac{1}{3}\\gamma^2+\\frac{1}{2}\\gamma\\right)u_{n}+\n", "\\left(\\frac{2}{3}\\gamma^2+\\frac{1}{2}\\gamma-1\\right)u_{n-1}+\n", "h\\left[\\left(-\\frac{1}{6}\\gamma^2-\\frac{1}{4}\\gamma+2\\right)\\varphi_{n}+\n", " \\left(-\\frac{1}{6}\\gamma^2-\\frac{1}{4}\\gamma\\right)\\varphi_{n-1}\\right].\n", "$$\n", "1. Pour quelles valeurs de $\\gamma$ est-elle consistante?\n", "1. Est-elle convergente pour ces valeurs?\n", "1. Quel ordre de convergente a-t-elle?\n", "1. Pour les valeurs de $\\gamma$ qui donnent la convergence, vérifier empiriquement l'ordre de convergence sur le problème de Cauchy\n", "$$\n", "\\begin{cases}\n", "y'(t) = -10\\big(y(t)-1\\big)^2, &\\forall t \\in I=[0,1],\\\\\n", "y(0) = 2.\n", "\\end{cases}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Correction\n", "\n", "C'est une méthode à $q=2$ pas explicite : \n", "- $p=1$\n", "- $b_{-1}=0$\n", "- $a_0=\\frac{1}{3}\\gamma^2+\\frac{1}{2}\\gamma$ et $a_1=\\frac{2}{3}\\gamma^2+\\frac{1}{2}\\gamma-1$\n", "- $b_0=-\\frac{1}{6}\\gamma^2-\\frac{1}{4}\\gamma+2$ et $b_1=-\\frac{1}{6}\\gamma^2-\\frac{1}{4}\\gamma$\n", "\n", "\n", "1. La méthode est consistante ssi\n", "$$\n", "\\begin{cases}\n", "\\displaystyle\\sum_{j=0}^p a_j=\\gamma^2+\\gamma-1=1,\n", "\\\\\n", "\\displaystyle-\\sum_{j=0}^p ja_j+\\sum_{j=-1}^p b_j=-\\big(0a_0+1a_1\\big)+\\big(b_{-1}+b_0+b_1\\big)=\\gamma^2+\\gamma-1=1\n", "\\end{cases}\n", "$$\n", "donc ssi $\\gamma=-2$ ou $\\gamma=1$.\n", "\n", "1. 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=\n", "\\begin{cases}\n", "r^2-\\frac{1}{3}r-\\frac{2}{3}&\\text{si }\\gamma=-2,\n", "\\\\\n", "r^2-\\frac{5}{6}r-\\frac{1}{6}&\\text{si }\\gamma=1,\n", "\\end{cases}\n", "$$\n", " dont les racines sont \n", "$$\n", "\\begin{cases}\n", "r_0=1,\\ r_1=-\\frac{2}{3}&\\text{si }\\gamma=-2,\n", "\\\\\n", "r_0=1,\\ r_1=-\\frac{1}{6}&\\text{si }\\gamma=1,\n", "\\end{cases}\n", "$$\n", " La méthode est zéro-stable ssi\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", "$$\n", " Cette condition est vérifiée pour les deux cas.\n", "\n", " La méthode est convergente ssi elle est consistante et zéro-stable donc ssi $\\gamma=-2$ ou $\\gamma=1$.\n", "\n", "1. La méthode est d'ordre au moins 2 pour les deux choix de $\\gamma$ car\n", "$$\n", "\\sum_{j=0}^p (-j)^{2}a_j+2\\sum_{j=-1}^p (-j)^{1}b_j\n", "{}=\n", "\\begin{cases}\n", "\\frac{2}{3}+2\\frac{1}{6}=1&\\text{si }\\gamma=-2,\n", "\\\\\n", "\\frac{1}{6}+2\\frac{5}{12}=1&\\text{si }\\gamma=1,\n", "\\end{cases}\n", "$$\n", " La méthode n'est pas d'ordre 3 pour aucun des deux choix de $\\gamma$ car\n", "$$\n", "\\sum_{j=0}^p (-j)^{3}a_j+3\\sum_{j=-1}^p (-j)^{2}b_j\n", "{}=\n", "\\begin{cases}\n", "{}-\\frac{2}{3}+3\\frac{-1}{6}\\neq1&\\text{si }\\gamma=-2,\n", "\\\\\n", "{}-\\frac{1}{6}+3\\frac{-5}{12}\\neq1&\\text{si }\\gamma=1,\n", "\\end{cases}\n", "$$\n", "\n", "1. On définit la solution exacte (calculée en utilisant le module `sympy`) pour estimer les erreurs:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%reset -f\n", "%matplotlib inline\n", "\n", "from matplotlib.pylab import *\n", "# rcdefaults()\n", "rcParams.update({'font.size': 16})\n", "\n", "# variables globales\n", "t0 = 0\n", "tfinal = 1\n", "y0 = 2\n", "\n", "phi = lambda t,y : -10*(y-1)**2" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "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) , phi(t,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": null, "metadata": {}, "outputs": [], "source": [ "from scipy.optimize import fsolve\n", "\n", "def multipas(gamma, phi, tt, sol_exacte):\n", " h = tt[1] - tt[0]\n", " uu = [sol_exacte(tt[0])]\n", " uu.append(sol_exacte(tt[1]))\n", " for i in range(1,len(tt) - 1):\n", " uu.append( (gamma**2/3+gamma/2)*uu[i]+(2*gamma**2/3+gamma/2-1)*uu[i-1]\n", " +h*(-gamma**2/6-gamma/4+2)*phi(tt[i],uu[i])+h*(-gamma**2/6-gamma/4)*phi(tt[i-1],uu[i-1]) )\n", " return uu\n", "\n", "\n", "H = []\n", "err_mp1 = []\n", "err_mp2 = []\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_mp1 = multipas( 1, phi, tt, sol_exacte)\n", " uu_mp2 = multipas(-2, phi, tt, sol_exacte)\n", " H.append(h)\n", " err_mp1.append(max([abs(uu_mp1[i] - yy[i]) for i in range(len(yy))]))\n", " err_mp2.append(max([abs(uu_mp2[i] - yy[i]) for i in range(len(yy))]))\n", "\n", "print ('Multipas $\\gamma$= 1 ordre = %1.2f' %(polyfit(log(H),log(err_mp1), 1)[0]))\n", "print ('Multipas $\\gamma$=-2 ordre = %1.2f' %(polyfit(log(H),log(err_mp2), 1)[0]))\n", "\n", "figure(figsize=(8,5))\n", "plot(log(H), log(err_mp1), 'r-o', label='Multipas b=1')\n", "plot(log(H), log(err_mp2), 'b-d', label='Multipas b=-2')\n", "xlabel('$\\ln(h)$')\n", "ylabel('$\\ln(e)$')\n", "legend(bbox_to_anchor=(1.04, 1), loc='upper left')\n", "grid(True)\n", "show()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercice : convergence\n", "\n", "On notera $\\varphi_k\\stackrel{\\text{déf}}{=}\\varphi(t_k,u_k)$.\n", "Soit la méthode multipas\n", "$$\n", "u_{n+1}\n", "=\\left(b^2-\\frac{9}{2}b+\\frac{11}{2}\\right)u_{n}\n", "+\\left(b^2-\\frac{7}{2}b+\\frac{3}{2}\\right)u_{n-1}\n", "+h\\left(\\frac{b}{2}(b-2)\\varphi_{n}\n", " +\\frac{b}{2}\\left(b-\\frac{10}{3}\\right)\\varphi_{n-1}\\right).\n", "$$\n", "- Pour quelles valeurs de $b$ est-elle consistante?\n", "- Est-elle convergente pour ces valeurs?\n", "- Quel ordre de convergente a-t-elle?\n", "- Pour les valeurs de $b$ qui donnent la convergence, vérifier empiriquement l'ordre de convergence sur le problème de Cauchy\n", "$$\n", "\\begin{cases}\n", "y'(t) = y(t)(1-y(t)), &\\forall t \\in I=[0,1],\\\\\n", "y(0) = 2.\n", "\\end{cases}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Correction\n", "\n", "C'est une méthode à $q=2$ pas explicite : \n", "- $p=1$\n", "- $b_{-1}=0$\n", "- $a_0=b^2-\\frac{9}{2}b+\\frac{11}{2}$ et $a_1=b^2-\\frac{7}{2}b+\\frac{3}{2}$\n", "- $b_0=\\frac{b}{2}(b-2)$ et $b_1=\\frac{b}{2}\\left(b-\\frac{10}{3}\\right)$\n", "\n", "\n", "- La méthode est consistante ssi\n", "$$\n", "\\begin{cases}\n", "\\displaystyle\\sum_{j=0}^p a_j=2b^2-8b+7=1,\n", "\\\\\n", "\\displaystyle-\\sum_{j=0}^p ja_j+\\sum_{j=-1}^p b_j=-\\big(0a_0+1a_1\\big)+\\big(b_{-1}+b_0+b_1\\big)=\\frac{5}{6}b-\\frac{19}{6}=1\n", "\\end{cases}\n", "$$\n", "donc ssi $b=3$. Vérifions ces calculs ci-dessous:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import sympy as sym\n", "sym.init_printing()\n", "\n", "sym.var('b')\n", "p=1\n", "bm1=0\n", "a0=b**2-sym.Rational(9,2)*b+sym.Rational(11,2)\n", "a1=b**2-sym.Rational(7,2)*b+sym.Rational(3,2)\n", "b0=sym.Rational(1,2)*b*(b-2)\n", "b1=sym.Rational(1,2)*b*(b-sym.Rational(10,3))\n", "\n", "cond1=a0+a1\n", "cond2=-(0*a0+1*a1)+(bm1+b0+b1)\n", "display(cond1)\n", "display(cond2.factor())\n", "display(sym.solve([sym.Eq(cond1,1),sym.Eq(cond2,1)],b))\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Dans la suite on considerera $b=3$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from IPython.display import display, Math\n", "display(Math('a_0='+sym.latex(a0.subs(b,3))))\n", "display(Math('a_1='+sym.latex(a1.subs(b,3))))\n", "display(Math('b_0='+sym.latex(b0.subs(b,3))))\n", "display(Math('b_1='+sym.latex(b1.subs(b,3))))" ] }, { "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\n", "\\stackrel{b=3}{=}\n", "r^2-r\n", "$$\n", " dont les racines sont \n", "$$\n", "r_0=1,\\ r_1=0\n", "$$\n", " La méthode est zéro-stable ssi\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", "$$\n", " ce qui est bien vérifié.\n", "\n", " La méthode est convergente ssi elle est consistante et zéro-stable donc ssi $b=3$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- La méthode est d'ordre au moins 2 car\n", "$$\n", "\\sum_{j=0}^p (-j)^{2}a_j+2\\sum_{j=-1}^p (-j)^{1}b_j\n", "{}=\n", "\\big(0^2a_0+1^2a_1\\big)+2\\big(-(-1)b_{-1}+0b_0+1b_1\\big)\n", "\\stackrel{b=3}{=}\n", "1\n", "$$\n", " La méthode n'est pas d'ordre 3 car\n", "$$\n", "\\sum_{j=0}^p (-j)^{3}a_j+3\\sum_{j=-1}^p (-j)^{2}b_j\n", "{}=\n", "{}-a_1+3\\big(b_{-1}+b_1\\big)\n", "\\stackrel{b=3}{\\neq}1\n", "$$" ] }, { "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": null, "metadata": {}, "outputs": [], "source": [ "%reset -f\n", "%matplotlib inline\n", "\n", "from matplotlib.pylab import *\n", "# rcdefaults()\n", "rcParams.update({'font.size': 16})\n", "\n", "# variables globales\n", "t0 = 0\n", "tfinal = 1\n", "y0 = 2\n", "\n", "phi = lambda t,y : y*(1-y)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "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) , phi(t,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": null, "metadata": {}, "outputs": [], "source": [ "from scipy.optimize import fsolve\n", "\n", "def multipas(phi, tt, sol_exacte):\n", " h = tt[1] - tt[0]\n", " uu = [sol_exacte(tt[0])]\n", " uu.append(sol_exacte(tt[1]))\n", " for i in range(1,len(tt) - 1):\n", " uu.append( uu[i]+h*3/2*phi(tt[i],uu[i])-h/2*phi(tt[i-1],uu[i-1]) )\n", " return uu\n", "\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, sol_exacte)\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 b= 3 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', label='Multipas b=3')\n", "xlabel('$\\ln(h)$')\n", "ylabel('$\\ln(e)$')\n", "legend(bbox_to_anchor=(1.04, 1), loc='upper left')\n", "grid(True)\n", "show()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercice : convergence predictor-corrector\n", "\n", "On notera $\\varphi_k\\stackrel{\\text{déf}}{=}\\varphi(t_k,u_k)$.\n", "Soit les méthodes multipas\n", "\\begin{align*}\n", "u_{n+1}&=u_{n}+\\frac{h}{12}\\left(23\\varphi_{n}-16\\varphi_{n-1}+5\\varphi_{n-2}\\right)&\\text{[P]}\n", "\\\\\n", "u_{n+1}&=u_{n}+\\frac{h}{24}\\left(9\\varphi_{n+1}+19\\varphi_{n}-5\\varphi_{n-1}+\\varphi_{n-2}\\right)&\\text{[C]}\n", "\\end{align*}\n", "- Étudier consistance, ordre et zéro-stabilité de la méthode P\n", "- Étudier consistance, ordre et zéro-stabilité de la méthode C\n", "- Soit la méthode PC suivante:\n", "$$\n", "\\begin{cases}\n", "\\tilde u=u_{n}+\\frac{h}{12}\\left(23\\varphi_{n}-16\\varphi_{n-1}+5\\varphi_{n-2}\\right)\n", "\\\\\n", "u_{n+1}=u_{n}+\\frac{h}{24}\\left(9\\varphi(t_{n+1},\\tilde u)+19\\varphi_{n}-5\\varphi_{n-1}+\\varphi_{n-2}\\right)\n", "\\end{cases}\n", "$$\n", " Calculer empiriquement l'ordre de convergence de la méthode PC sur le problème de Cauchy\n", "$$\n", "\\begin{cases}\n", "y'(t) = \\frac{1}{1+t^2}-2y^2(t), &\\forall t \\in I=[0,1],\\\\\n", "y(0) = 0,\n", "\\end{cases}\n", "$$\n", " dont la solution exacte est $y(t)=\\frac{t}{1+t^2}$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Correction schéma P\n", "\n", "P est une méthode à $q=3$ pas explicite : \n", "- $p=2$\n", "- $b_{-1}=0$\n", "- $a_0=1$ et $a_1=a_2=0$\n", "- $b_0=\\frac{23}{12}$, $b_1=\\frac{-16}{12}$ et $b_2=\\frac{5}{12}$\n", "- La méthode est consistante d'ordre $\\omega=3$ car\n", "$$\n", "\\begin{cases}\n", "\\displaystyle\\sum_{j=0}^p a_j=1,\n", "\\\\\n", "\\displaystyle\\sum_{j=0}^p -ja_j+\\sum_{j=-1}^p b_j=-a_1-2a_2+\\big(b_{-1}+b_0+b_1+b_2\\big)=1\n", "\\\\\n", "\\displaystyle\\sum_{j=0}^p (-j)^{2}a_j+2\\sum_{j=-1}^p (-j)^{1}b_j=a_1+4a_2-2\\big(-b_{-1}+b_1+2b_2\\big)=1\n", "\\\\\n", "\\displaystyle\\sum_{j=0}^p (-j)^{3}a_j+3\\sum_{j=-1}^p (-j)^{2}b_j=-a_1-8a_2+3\\big(b_{-1}+b_1+4b_2\\big)=1\n", "\\\\\n", "\\displaystyle\\sum_{j=0}^p (-j)^{4}a_j+4\\sum_{j=-1}^p (-j)^{3}b_j=a_1+16a_2-4\\big(-b_{-1}+b_1+8b_2\\big)\\neq1\n", "\\end{cases}\n", "$$\n", "Vérifions ces calculs ci-dessous:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$\\displaystyle \\omega\\ge1$" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/latex": [ "$\\displaystyle \\omega\\ge2$" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/latex": [ "$\\displaystyle \\omega\\ge3$" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/latex": [ "$\\displaystyle \\omega\\le4$" ], "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", "p=1\n", "bm1=0\n", "a0=1\n", "a1=0\n", "a2=0\n", "b0=sym.Rational(23,12)\n", "b1=sym.Rational(-16,12)\n", "b2=sym.Rational(5,12)\n", "\n", "\n", "aa=[a0,a1,a2]\n", "bb=[b0,b1,b2]\n", "\n", "cond=[]\n", "cond.append(sum(aa)) \n", "cond.append(sum([-j*aa[j]+bb[j] for j in range(len(aa))])+bm1)\n", "if cond==[1,1]:\n", " display(Math(\"\\omega\\ge1\"))\n", " for n in range(2,9):\n", " cond.append(sum( [ (-j)**n*aa[j]+n*(-j)**(n-1)*bb[j] for j in range(len(aa)) ])+n*bm1)\n", " if cond[n]==1:\n", " display(Math(\"\\omega\\ge\"+str(n)))\n", " else:\n", " display(Math(\"\\omega\\le\"+str(n)))\n", " break" ] }, { "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^3-a_0r^2-a_1r-a_2\n", "{}=\n", "r^3-r^2=r^2(r-1)\n", "$$\n", " dont les racines sont \n", "$$\n", "r_0=1,\\ r_1=r_2=0\n", "$$\n", " La méthode est zéro-stable ssi\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", "$$\n", " ce qui est bien vérifié.\n", "\n", " La méthode est convergente car consistante et zéro-stable." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Correction schéma C\n", "\n", "C est une méthode à $q=3$ pas implicite : \n", "- $p=2$\n", "- $b_{-1}=\\frac{9}{24}$\n", "- $a_0=1$ et $a_1=a_2=0$\n", "- $b_0=\\frac{19}{24}$, $b_1=\\frac{-5}{24}$ et $b_2=\\frac{1}{24}$\n", "- La méthode est consistante d'ordre $\\omega=4$ car\n", "$$\n", "\\begin{cases}\n", "\\displaystyle\\sum_{j=0}^p a_j=1,\n", "\\\\\n", "\\displaystyle\\sum_{j=0}^p -ja_j+\\sum_{j=-1}^p b_j=-a_1-2a_2+\\big(b_{-1}+b_0+b_1+b_2\\big)=1\n", "\\\\\n", "\\displaystyle\\sum_{j=0}^p (-j)^{2}a_j+2\\sum_{j=-1}^p (-j)^{1}b_j=a_1+4a_2-2\\big(-b_{-1}+b_1+2b_2\\big)=1\n", "\\\\\n", "\\displaystyle\\sum_{j=0}^p (-j)^{3}a_j+3\\sum_{j=-1}^p (-j)^{2}b_j=-a_1-8a_2+3\\big(b_{-1}+b_1+4b_2\\big)=1\n", "\\\\\n", "\\displaystyle\\sum_{j=0}^p (-j)^{4}a_j+4\\sum_{j=-1}^p (-j)^{3}b_j=a_1+16a_2-4\\big(-b_{-1}+b_1+8b_2\\big)=1\n", "\\\\\n", "\\displaystyle\\sum_{j=0}^p (-j)^{5}a_j+5\\sum_{j=-1}^p (-j)^{4}b_j=-a_1+32a_2+5\\big(b_{-1}+b_1+16b_2\\big)\\neq1\n", "\\end{cases}\n", "$$\n", "Vérifions ces calculs ci-dessous:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import sympy as sym\n", "sym.init_printing()\n", "from IPython.display import display, Math\n", " \n", "p=1\n", "bm1=sym.Rational(9,24)\n", "a0=1\n", "a1=0\n", "a2=0\n", "b0=sym.Rational(19,24)\n", "b1=sym.Rational(-5,24)\n", "b2=sym.Rational(1,24)\n", "\n", "\n", "aa=[a0,a1,a2]\n", "bb=[b0,b1,b2]\n", "\n", "cond=[]\n", "cond.append(sum(aa)) \n", "cond.append(sum([-j*aa[j]+bb[j] for j in range(len(aa))])+bm1)\n", "if cond==[1,1]:\n", " display(Math(\"\\omega\\ge1\"))\n", " for n in range(2,9):\n", " cond.append(sum( [ (-j)**n*aa[j]+n*(-j)**(n-1)*bb[j] for j in range(len(aa)) ])+n*bm1)\n", " if cond[n]==1:\n", " display(Math(\"\\omega\\ge\"+str(n)))\n", " else:\n", " display(Math(\"\\omega\\le\"+str(n)))\n", " break" ] }, { "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^3-a_0r^2-a_1r-a_2\n", "{}=\n", "r^3-r^2=r^2(r-1)\n", "$$\n", " dont les racines sont \n", "$$\n", "r_0=1,\\ r_1=r_2=0\n", "$$\n", " La méthode est zéro-stable ssi\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", "$$\n", " ce qui est bien vérifié.\n", "\n", " La méthode est convergente car consistante et zéro-stable." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Correction schéma PC\n", "\n", "Il est d'ordre $4$ car le prédictor est d'ordre $3$ et le corrector d'ordre $4$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On définit la solution exacte pour estimer les erreurs:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%reset -f\n", "%matplotlib inline\n", "\n", "from matplotlib.pylab import *\n", "# rcdefaults()\n", "rcParams.update({'font.size': 16})\n", "\n", "# variables globales\n", "t0 = 0\n", "tfinal = 1\n", "y0 = 0\n", "\n", "phi = lambda t,y : 1/(1+t**2)-2*y**2\n", "\n", "sol_exacte = lambda t : t/(1+t**2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On calcule la solution approchée pour différentes valeurs de $N$ pour les trois schémas:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from scipy.optimize import fsolve\n", "\n", "\n", "def multipas_P(phi, tt, sol_exacte):\n", " h = tt[1] - tt[0]\n", " uu = [sol_exacte(tt[0])]\n", " uu.append(sol_exacte(tt[1]))\n", " uu.append(sol_exacte(tt[2]))\n", " for i in range(2,len(tt) - 1):\n", " uu.append( uu[i]+h/12* ( 23*phi(tt[i],uu[i])-16*phi(tt[i-1],uu[i-1])+5*phi(tt[i-2],uu[i-2]) ) )\n", " return uu\n", "\n", "def multipas_C(phi, tt, sol_exacte):\n", " h = tt[1] - tt[0]\n", " uu = [sol_exacte(tt[0])]\n", " uu.append(sol_exacte(tt[1]))\n", " uu.append(sol_exacte(tt[2]))\n", " for i in range(2,len(tt) - 1):\n", " temp = fsolve ( lambda x : -x+uu[i]+h/24* ( 9*phi(tt[i+1],x) + 19*phi(tt[i],uu[i])-5*phi(tt[i-1],uu[i-1])+phi(tt[i-2],uu[i-2]) ), uu[i])\n", " uu.append(temp)\n", " return uu\n", "\n", "def multipas_PC(phi, tt, sol_exacte):\n", " h = tt[1] - tt[0]\n", " uu = [sol_exacte(tt[0])]\n", " uu.append(sol_exacte(tt[1]))\n", " uu.append(sol_exacte(tt[2]))\n", " for i in range(2,len(tt) - 1):\n", " u_tilde = uu[i]+h/12* ( 23*phi(tt[i],uu[i])-16*phi(tt[i-1],uu[i-1])+5*phi(tt[i-2],uu[i-2]) )\n", " uu.append( uu[i]+h/24* ( 9*phi(tt[i+1],u_tilde) + 19*phi(tt[i],uu[i])-5*phi(tt[i-1],uu[i-1])+phi(tt[i-2],uu[i-2]) ) )\n", " return uu\n", "\n", "\n", "H = []\n", "err_p = []\n", "err_c = []\n", "err_pc = []\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_p = multipas_P(phi, tt, sol_exacte)\n", " uu_c = multipas_C(phi, tt, sol_exacte)\n", " uu_pc = multipas_PC(phi, tt, sol_exacte)\n", " H.append(h)\n", " err_p.append(max([abs(uu_p[i] - yy[i]) for i in range(len(yy))]))\n", " err_c.append(max([abs(uu_c[i] - yy[i]) for i in range(len(yy))]))\n", " err_pc.append(max([abs(uu_pc[i] - yy[i]) for i in range(len(yy))]))\n", " \n", "print (f'Multipas P ordre = {polyfit(log(H),log(err_p) , 1)[0]:1.2f}' )\n", "print (f'Multipas C ordre = {polyfit(log(H),log(err_c) , 1)[0]:1.2f}' )\n", "print (f'Multipas PC ordre = {polyfit(log(H),log(err_pc), 1)[0]:1.2f}')\n", "\n", "figure(figsize=(8,5))\n", "plot(log(H), log(err_p), 'r-o', label='Multipas P')\n", "plot(log(H), log(err_c), 'b-o', label='Multipas C')\n", "plot(log(H), log(err_pc), 'c-o', label='Multipas PC')\n", "xlabel('$\\ln(h)$')\n", "ylabel('$\\ln(e)$')\n", "legend(bbox_to_anchor=(1.04, 1), loc='upper left')\n", "grid(True);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Rappels schémas Runge-Kutta\n", "\n", "\n", "
\n", "\n", "Une méthode de Runge-Kutta à $s\\ge1$ étages s'écrit:\n", "$$\\begin{array}{c|c}\n", "\\mathbf{c} & \\mathbb{A}\\\\\n", "\\hline\n", "&\\mathbf{b}^T\n", "\\end{array}\n", "\\quad\\implies\\quad\n", "\\begin{cases}\n", "u_0=y(t_0)=y_0,\\\\\n", "u_{n+1}=u_{n}+h \\displaystyle\\sum_{i=1}^sb_iK_i& n=0,1,\\dots N-1\\\\\n", "K_i=\\displaystyle\\varphi\\left( t_n+hc_i,u_n+h\\sum_{j=1}^{s}a_{ij}K_j \\right) & i=1,\\dots s.\n", "\\end{cases}$$\n", "\n", "
\n", "\n", "1. Soit $\\omega$ l'ordre de la méthode.\n", " - Si $\\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}$ alors la méthode est consistante (i.e. $\\omega\\ge1$)\n", "\n", " - Si, de plus, $\\displaystyle\\sum_{j=1}^s b_j c_j=\\frac{1}{2}$ alors $\\omega\\ge2$\n", " - Si, de plus, $\\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", " - Si, de plus, $\\begin{cases}\\displaystyle\\sum_{j=1}^s b_j c_j^3=\\frac{1}{4}&\\\\\\displaystyle\\sum_{i=1}^s\\sum_{j=1}^s b_i c_i a_{ij} c_j=\\frac{1}{8}\\\\\\displaystyle\\sum_{i=1}^s\\sum_{j=1}^s b_i a_{ij} c_j^2=\\frac{1}{12}\\\\\\displaystyle\\sum_{i=1}^s\\sum_{j=1}^s\\sum_{k=1}^s b_i a_{ij}a_{jk} c_k=\\frac{1}{24}\\end{cases}$ alors $\\omega\\ge4$ \n", "\n", "1. La méthode est **zéro-stable**\n", "\n", "1. **Consistance + zéro-stabilité = convergence**\n", "\n", "1. 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}\\xrightarrow[t\\to+\\infty]{}0.$\n", " Si, sous d'éventuelles conditions sur $h$, on a $|u_n|\\xrightarrow[n\\to+\\infty]{}0$ alors on dit que **le schéma est A-stable.**\n", "\n", " Une méthode de Runge-Kutta à $s\\ge1$ étages pour $y'(t)=-\\beta y(t)$ s'écrit:\n", " $$\\begin{cases}\n", " u_0=y_0,\\\\\n", " u_{n+1}=u_{n}+h \\displaystyle\\sum_{i=1}^sb_iK_i& n=0,1,\\dots N-1\\\\\n", " K_i=\\displaystyle -\\beta \\left( u_n+h\\sum_{j=1}^{s}a_{ij}K_j \\right) & i=1,\\dots s.\n", " \\end{cases}$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercice\n", "Dans la littérature on trouve beaucoup de schémas de Runge-Kutta. Ci-dessous une petite liste de ces méthodes. Pour chaque méthode de Runge-Kutta indiquée, il est indiqué le nombre d'étages $s$, son ordre de convergence $\\omega$ et si le schéma est explicite, semi-implicite ou implicite.\n", "\n", "Exercice: choisir une méthode et vérifier si l'ordre de convergence est correct." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "| \t| Matrice de Butcher \t| s \t| $\\omega$ \t| Exp/Imp/Semi-Implicite \t|\n", "|--------------\t|-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------\t|---\t|----------\t|----\t|\n", "| Gauss-1 \t| $\\begin{array}{c|cc} \\frac{1}{2} & \\frac{1}{2}\\\\ \\hline & 1 \\end{array}$ \t| 1 \t| 2 \t| I \t|\n", "| EE \t| $\\begin{array}{c|cc} 1 & 1\\\\ \\hline & 1 \\end{array}$ \t| 1 \t| 1 \t| E \t|\n", "| EI \t| $\\begin{array}{c|cc} 0 & 0\\\\ \\hline & 1 \\end{array}$ \t| 1 \t| 1 \t| I \t|" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "| \t| Matrice de Butcher \t| s \t| $\\omega$ \t| Exp/Imp/Semi-Implicite \t|\n", "|--------------\t|-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------\t|---\t|----------\t|----\t|\n", "| SDIRK order 3 | $\\begin{array}{c|cc} \\gamma & \\gamma & 0\\\\ 1-\\gamma & 1-2\\gamma & \\gamma\\\\ \\hline & \\frac{1}{2} & \\frac{1}{2} \\end{array}$ avec $\\gamma=\\frac{3\\pm\\sqrt{3}}{6}$ \t| 2 \t| 3 \t| SI \t| \n", "| Gauss-2 \t| $\\begin{array}{c|cc} \\frac{1}{2}-\\frac{\\sqrt{3}}{6} & \\frac{1}{4} & \\frac{1}{4}-\\frac{\\sqrt{3}}{6}\\\\ \\frac{1}{2}+\\frac{\\sqrt{3}}{6} & \\frac{1}{4}+\\frac{\\sqrt{3}}{6} & \\frac{1}{4}\\\\ \\hline & \\frac{1}{2} & \\frac{1}{2} \\end{array}$ | 2 \t| 4 \t| I \t|\n", "| DIRK-2 \t| $\\begin{array}{c|cc} \\frac{1}{3} & \\frac{1}{3} & 0\\\\ 1 & 1 & 0\\\\ \\hline & \\frac{3}{4} & \\frac{1}{4} \\end{array}$ \t| 2 \t| 3 \t| SI \t|\n", "| Hammer & Hollingsworth \t| $\\begin{array}{c|cc} 0 & 0 & 0\\\\ \\frac{2}{3} & \\frac{1}{3} & \\frac{1}{3}\\\\ \\hline & \\frac{1}{4} & \\frac{3}{4} \\end{array}$ \t| 2 \t| ? \t| SI \t|\n", "| Randau-IIa-2 \t| $\\begin{array}{c|cc} \\frac{1}{3} & \\frac{5}{12} & -\\frac{1}{12}\\\\ 1 & \\frac{3}{4} & \\frac{1}{4}\\\\ \\hline & \\frac{3}{4} & \\frac{1}{4} \\end{array}$ \t| 2 \t| 3 \t| I \t|\n", "| CN \t| $\\begin{array}{c|cc} 0 & 0 & 0\\\\ 1 & \\frac{1}{2} & \\frac{1}{2}\\\\ \\hline & \\frac{1}{2} & \\frac{1}{2} \\end{array}$ \t| 2 \t| 2 \t| I \t|\n", "| Heun \t| $\\begin{array}{c|cc} 0 & 0 & 0\\\\ 1 & 1 & 0\\\\ \\hline & \\frac{1}{2} & \\frac{1}{2} \\end{array}$ \t| 2 \t| 2 \t| E \t|\n", "| EM \t| $\\begin{array}{c|cc} 0 & 0 & 0\\\\ \\frac{1}{2} & \\frac{1}{2} & 0\\\\ \\hline & 0 & 1 \\end{array}$ \t| 2 \t| 2 \t| E \t|\n", "| \t| $\\begin{array}{c|cc} 0 & 0 & 0\\\\ \\frac{2}{3} & \\frac{2}{3} & 0\\\\ \\hline & \\frac{1}{4} & \\frac{3}{4}\\end{array}$ \t| 2 \t| 2 \t| E \t|\n", "| Ralston \t| $\\begin{array}{c|cc} 0 & 0 & 0\\\\ \\frac{3}{4} & \\frac{3}{4} & 0\\\\ \\hline & \\frac{1}{3} & \\frac{2}{3} \\end{array}$ \t| 2 \t| 2 \t| E \t|\n", "| \t| $\\begin{array}{c|cc} 0 & 0 & 0\\\\ \\frac{3}{4} & \\frac{3}{4} & 0\\\\ \\hline & 0 & 1 \\end{array}$ \t| 2 \t| 1 \t| E \t|\n", "| \t| $\\begin{array}{c|cc} 0 & 0 & 0\\\\ \\frac{2}{3} & \\frac{2}{3} & 0\\\\ \\hline & \\frac{1}{2} & \\frac{1}{2} \\end{array}$ \t| 2 \t| 1 \t| E \t|" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "| \t| Matrice de Butcher \t| s \t| $\\omega$ \t| Exp/Imp/Semi-Implicite \t|\n", "|--------------\t|-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------\t|---\t|----------\t|----\t|\n", "| Gauss-3 \t| $\\begin{array}{c|ccc} \\frac{1}{2}-\\frac{\\sqrt{15}}{10} & \\frac{5}{36} & \\frac{2}{9}-\\frac{\\sqrt{15}}{15} & \\frac{5}{36}-\\frac{\\sqrt{15}}{30} \\\\ \\frac{1}{2} & \\frac{5}{36}+\\frac{\\sqrt{15}}{24} & \\frac{2}{9}&\\frac{5}{36}-\\frac{\\sqrt{15}}{24}\\\\ \\frac{1}{2}+\\frac{\\sqrt{15}}{10} & \\frac{5}{36}+\\frac{\\sqrt{15}}{30} &\\frac{2}{9}+\\frac{\\sqrt{15}}{15}&\\frac{5}{36}\\\\ \\hline & \\frac{5}{18} &\\frac{4}{9}& \\frac{5}{18} \\end{array}$ \t| 3 \t| 6 \t| I \t|\n", "| Radau-IIa-3 \t| $\\begin{array}{c|ccc} \\frac{4-\\sqrt{6}}{10} & \\frac{88-7\\sqrt{6}}{360} & \\frac{296-169\\sqrt{6}}{1800} & \\frac{-2+3\\sqrt{6}}{225} \\\\ \\frac{4+\\sqrt{6}}{10} & \\frac{296+169\\sqrt{6}}{1800} & \\frac{88+7\\sqrt{6}}{360}&\\frac{-2-3\\sqrt{6}}{225}\\\\ 1 & \\frac{16-\\sqrt{6}}{36} &\\frac{16+\\sqrt{6}}{36}&\\frac{1}{9}\\\\ \\hline & \\frac{16-\\sqrt{6}}{36} &\\frac{16+\\sqrt{6}}{36}& \\frac{1}{9} \\end{array}$ \t| 3 \t| 5 \t| I \t|\n", "| Lobatto IIIa \t| $\\begin{array}{c|ccc} 0 & 0 & 0 & 0 \\\\ \\frac{1}{2} & \\frac{5}{24} & \\frac{1}{3}&-\\frac{1}{24}\\\\ 1 & \\frac{1}{6} &\\frac{2}{3}&\\frac{1}{6}\\\\ \\hline & \\frac{1}{6} & \\frac{2}{3} & \\frac{1}{6} \\end{array}$ \t| 3 \t| 4 \t| I \t|\n", "| Lobatto \t| $\\begin{array}{c|ccc} 0 & 0 & 0 & 0 \\\\ \\frac{1}{2} & \\frac{1}{4} & \\frac{1}{4}&0\\\\ 1 & 0 &1&0\\\\ \\hline & \\frac{1}{6} & \\frac{2}{3} & \\frac{1}{6} \\end{array}$ \t| 3 \t| 4 \t| SI \t|\n", "| Heun3 \t| $\\begin{array}{c|ccc} 0 & 0 & 0 & 0 \\\\ \\frac{1}{3} & \\frac{1}{3} & 0&0\\\\ \\frac{2}{3} & 0 &\\frac{2}{3}&0\\\\ \\hline & \\frac{1}{4} & 0 & \\frac{3}{4} \\end{array}$ \t| 3 \t| 3 \t| E \t|\n", "| Kutta \t| $\\begin{array}{c|ccc} 0 & 0 & 0 & 0 \\\\ \\frac{1}{2} & \\frac{1}{2} & 0&0\\\\ 1 & -1 &2&0\\\\ \\hline & \\frac{1}{6} & \\frac{2}{3} & \\frac{1}{6} \\end{array}$ \t| 3 \t| 3 \t| E \t|\n", "| \t| $\\begin{array}{c|ccc} 0 & 0 & 0 & 0 \\\\ \\frac{1}{2} & \\frac{1}{2} & 0&0\\\\ 1 & -1 &2&0\\\\ \\hline & -\\frac{1}{6} & \\frac{4}{3} & -\\frac{1}{6} \\end{array}$ \t| 3 \t| 3 \t| E \t|\n", "| \t| $\\begin{array}{c|ccc} 0 & 0 & 0 & 0 \\\\ \\frac{1}{2} & \\frac{1}{2} & 0&0\\\\ 1 & -1 &2&0\\\\ \\hline & \\frac{1}{6} & \\frac{2}{3} & \\frac{1}{6} \\end{array}$ \t| 3 \t| 3 \t| E \t|\n", "| \t| $\\begin{array}{c|ccc} 0 & 0 & 0 & 0 \\\\ 1 & 1 & 0&0\\\\ \\frac{1}{2} & \\frac{1}{4} &\\frac{1}{4}&0\\\\ \\hline & \\frac{1}{6} & \\frac{1}{6} & \\frac{2}{3} \\end{array}$ \t| 3 \t| 3 \t| E \t|\n", "| \t| $\\begin{array}{c|ccc} 0 & 0 & 0 & 0 \\\\ \\frac{2}{3} & \\frac{2}{3} & 0&0\\\\ \\frac{2}{3} & 0 &\\frac{2}{3}&0\\\\ \\hline & \\frac{1}{4} & \\frac{3}{8} & \\frac{3}{8} \\end{array}$ \t| 3 \t| 3 \t| E \t|\n", "| \t| $\\begin{array}{c|ccc} 0 & 0 & 0 & 0 \\\\ \\frac{2}{3} & \\frac{2}{3} & 0&0\\\\ \\frac{2}{3} & \\frac{1}{3} &\\frac{1}{3}&0\\\\ \\hline & \\frac{1}{4} & 0 & \\frac{3}{4} \\end{array}$ \t| 3 \t| 3 \t| E \t|\n", "| \t| $\\begin{array}{c|ccc} 0 & 0 & 0 & 0 \\\\ \\frac{2}{3} & \\frac{2}{3} & 0&0\\\\ 0 & -1 &1&0\\\\ \\hline & 0 & \\frac{3}{4} & \\frac{1}{4} \\end{array}$ \t| 3 \t| 3 \t| E \t|\n", "| \t| $\\begin{array}{c|ccc} 0 & 0 & 0 & 0 \\\\ \\frac{1}{3} & \\frac{1}{3} & 0&0\\\\ 1 & -1 &2&0\\\\ \\hline & 0 & \\frac{3}{4} & \\frac{1}{4} \\end{array}$ \t| 3 \t| 3 \t| E \t|\n", "| \t| $\\begin{array}{c|ccc} 0 & 0 & 0 & 0 \\\\ \\frac{1}{3} & \\frac{1}{3} & 0&0\\\\ \\frac{2}{3} & \\frac{1}{3} &\\frac{1}{3}&0\\\\ \\hline & \\frac{1}{3} & \\frac{1}{3} & \\frac{1}{3} \\end{array}$ \t| 3 \t| 1 \t| E \t|" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "| \t| Matrice de Butcher \t| s \t| $\\omega$ \t| Exp/Imp/Semi-Implicite \t|\n", "|--------------\t|-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------\t|---\t|----------\t|----\t|\n", "|Lobatto \t| $\\begin{array}{c|cccc} 0 & 0 & 0 & 0 & 0\\\\ \\frac{5-\\sqrt{5}}{10} & \\frac{5+\\sqrt{5}}{60} & \\frac{1}{6} & \\frac{15-7\\sqrt{5}}{60} & 0\\\\ \\frac{5+\\sqrt{5}}{10} & \\frac{5-\\sqrt{5}}{60} &\\frac{15+7\\sqrt{5}}{60} &\\frac{1}{6}&0\\\\ 1 & \\frac{1}{6} & \\frac{5-\\sqrt{5}}{12} & \\frac{5+\\sqrt{5}}{12}& 0 \\\\ \\hline & \\frac{1}{12} & \\frac{5}{12} & \\frac{5}{12} & \\frac{1}{12} \\end{array}$ \t| 4 \t| 6 \t| I \t|| \t| $\\begin{array}{c|cccc} 0 & 0 & 0 & 0 & 0\\\\ \\frac{1}{4} & \\frac{1}{8} & \\frac{1}{8} & 0 & 0\\\\ \\frac{7}{10} & -\\frac{1}{100} &\\frac{14}{25} &\\frac{3}{20}&0\\\\ 1 & \\frac{2}{7} & 0 & \\frac{5}{7} & 0\\\\ \\hline & \\frac{1}{14} & \\frac{32}{81} & \\frac{250}{567} & \\frac{5}{54} \\end{array}$ \t| 4 \t| 5 \t| SI \t|\n", "| RK4-1 \t| $\\begin{array}{c|cccc} 0 & 0 & 0 & 0 & 0\\\\ \\frac{1}{2} & \\frac{1}{2} & 0 & 0 & 0\\\\ \\frac{1}{2} & 0 &\\frac{1}{2} &0&0\\\\ 1 & 0 & 0 & 1 & 0\\\\ \\hline & \\frac{1}{6} & \\frac{1}{3} & \\frac{1}{3} & \\frac{1}{6} \\end{array}$ \t| 4 \t| 4 \t| E \t|\n", "| RK4-2 \t| $\\begin{array}{c|cccc} 0 & 0 & 0 & 0 & 0\\\\ \\frac{1}{4} & \\frac{1}{4} & 0 & 0 & 0\\\\ \\frac{1}{2} & 0 &\\frac{1}{2} &0&0\\\\ 1 & 1 & -2 & 2 & 0\\\\ \\hline & \\frac{1}{6} & 0 & \\frac{2}{3} & \\frac{1}{6} \\end{array}$ \t| 4 \t| 4 \t| E \t|\n", "| Règle 3/8 \t| $\\begin{array}{c|cccc} 0 & 0 & 0 & 0 & 0\\\\ \\frac{1}{3} & \\frac{1}{3} & 0 & 0 & 0\\\\ \\frac{2}{3} &-\\frac{1}{3} & 1 &0&0\\\\ 1 & 1 & -1 & 1 & 0\\\\ \\hline & \\frac{1}{8} & \\frac{3}{8} & \\frac{3}{8} & \\frac{1}{8} \\end{array}$ \t| 4 \t| 4 \t| E \t|\n", "| \t| $\\begin{array}{c|cccc} 0 & 0 & 0 & 0 & 0\\\\ \\frac{1}{2} & \\frac{1}{2} & 0 & 0 & 0\\\\ 0 &-1 & 1 &0&0\\\\ 1 & -1 & \\frac{3}{2} & \\frac{1}{2} & 0\\\\ \\hline & \\frac{1}{12} & \\frac{2}{3} & \\frac{1}{12} & \\frac{1}{6} \\end{array}$ \t| 4 \t| 4 \t| E \t|\n", "| \t| $\\begin{array}{c|cccc} 0 & 0 & 0 & 0 & 0\\\\ 1 & 1 & 0 & 0 & 0\\\\ \\frac{1}{2} &\\frac{3}{8} & \\frac{1}{8} &0&0\\\\ 1 & -2 & -1 & 4 & 0\\\\ \\hline & \\frac{1}{6} & \\frac{1}{12} & \\frac{2}{3} & \\frac{1}{12} \\end{array}$ \t| 4 \t| 4 \t| E \t|\n", "| \t| $\\begin{array}{c|cccc} 0 & 0 & 0 & 0 & 0\\\\ \\frac{2}{3} & \\frac{2}{3} & 0 & 0 & 0\\\\ \\frac{1}{3} &\\frac{1}{12} & \\frac{1}{4} &0&0\\\\ 1 &-\\frac{5}{4} & \\frac{1}{4} & 2 & 0\\\\ \\hline & \\frac{1}{8} & \\frac{3}{8} & \\frac{3}{8} & \\frac{1}{8} \\end{array}$ \t| 4 \t| 4 \t| E \t|\n", "| \t| $\\begin{array}{c|cccc} 0 & 0 & 0 & 0 & 0\\\\ \\frac{1}{2} & \\frac{1}{2} & 0 & 0 & 0\\\\ 1 &0 & 1 &0&0\\\\ 1 &0 & 0 & 1 & 0\\\\ \\hline & \\frac{1}{6} & \\frac{2}{3} & 0 & \\frac{1}{6} \\end{array}$ \t| 4 \t| 3 \t| E \t|\n", "| \t| $\\begin{array}{c|cccc} 0 & 0 & 0 & 0 & 0\\\\ \\frac{1}{2} & \\frac{1}{2} & 0 & 0 & 0\\\\ \\frac{1}{3} &\\frac{1}{3} & 0 &0&0\\\\ \\frac{2}{3} &\\frac{1}{3} & 0 & \\frac{1}{3} & 0\\\\ \\hline & 0 & -2 & \\frac{3}{2} & \\frac{3}{2} \\end{array}$ \t| 4 \t| 3 \t| E \t|\n", "| \t| $\\begin{array}{c|cccc} 0 & 0 & 0 & 0 & 0\\\\ \\frac{1}{2} & \\frac{1}{2} & 0 & 0 & 0\\\\ \\frac{1}{3} &\\frac{1}{3} & 0 &0&0\\\\ \\frac{2}{3} &\\frac{1}{3} & 0 & \\frac{1}{3} & 0\\\\ \\hline & 0 & -1 & 1 & 1 \\end{array}$ \t| 4 \t| 2 \t| E \t|" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "| \t| Matrice de Butcher \t| s \t| $\\omega$ \t| Exp/Imp/Semi-Implicite \t|\n", "|--------------\t|-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------\t|---\t|----------\t|----\t|\n", "| Merson \t| $\\begin{array}{c|ccccc} 0 & 0 & 0 & 0 & 0 &0 \\\\ \\frac{1}{3} & \\frac{1}{3} & 0&0&0&0\\\\ \\frac{1}{3} & \\frac{1}{6} &\\frac{1}{6}&0&0&0\\\\ \\frac{1}{2} & \\frac{1}{8} &0&\\frac{3}{8}&0&0\\\\ 1 & \\frac{1}{2} &0&-\\frac{3}{2}&2&0\\\\ \\hline & \\frac{1}{6} & 0 & 0 & \\frac{2}{3} & \\frac{1}{6} \\end{array}$ \t| 5 \t| 5 \t| E \t|" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "| \t| Matrice de Butcher \t| s \t| $\\omega$ \t| Exp/Imp/Semi-Implicite \t|\n", "|--------------\t|-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------\t|---\t|----------\t|----\t|\n", "| Fehlberg | $\\begin{array}{c|cccccc} 0 & 0 & 0 & 0 & 0 &0&0 \\\\ \\frac{1}{4} & \\frac{1}{4} & 0&0&0&0&0\\\\ \\frac{3}{8} & \\frac{3}{32} &\\frac{9}{32}&0&0&0&0\\\\ \\frac{12}{13} & \\frac{1932}{2197} &-\\frac{7200}{2197}&\\frac{7269}{2197}&0&0&0\\\\ 1 & \\frac{439}{216} &-8&\\frac{3680}{513}&-\\frac{845}{4104}&0&0\\\\ \\frac{1}{2} & -\\frac{8}{27} &2&-\\frac{3544}{2565}&\\frac{1859}{4104}&-\\frac{11}{40}&0\\\\ \\hline & \\frac{25}{216} & 0&\\frac{1408}{2565} & \\frac{2197}{4104} & -\\frac{1}{5} & 0 \\end{array}$ \t| 6 \t| 4 \t| E \t|\n", "| Butcher \t| $\\begin{array}{c|cccccc} 0 & 0 & 0 & 0 & 0 &0&0 \\\\ \\frac{1}{4} & \\frac{1}{4} & 0&0&0&0&0\\\\ \\frac{1}{4} & \\frac{1}{8} &\\frac{1}{8}&0&0&0&0\\\\ \\frac{1}{2} & 0 &-\\frac{1}{2}&1&0&0&0\\\\ \\frac{3}{4} & \\frac{3}{16} &0&0&\\frac{9}{16}&0&0\\\\ 1 & \\frac{-3}{7} &\\frac{2}{7}&\\frac{12}{7}&\\frac{-12}{7}&\\frac{8}{7}&0\\\\ \\hline & \\frac{7}{90} & 0&\\frac{32}{90} & \\frac{12}{90} & \\frac{32}{90} & \\frac{7}{90} \\end{array}$ \t| 6 \t| 5 \t| E \t|\n", "| Nystrøm \t| $\\begin{array}{c|cccccc} 0 & 0 & 0 & 0 & 0 &0&0 \\\\ \\frac{1}{3} & \\frac{1}{3} & 0&0&0&0&0\\\\ \\frac{2}{5} & \\frac{4}{25} &\\frac{6}{25}&0&0&0&0\\\\ 1 & \\frac{1}{4} &-3&\\frac{15}{4}&0&0&0\\\\ \\frac{2}{3} & \\frac{2}{27} &\\frac{10}{9}&-\\frac{50}{81}&\\frac{8}{81}&0&0\\\\ \\frac{4}{5} & \\frac{2}{25} &\\frac{12}{25}&\\frac{2}{15}&\\frac{8}{75}&0&0\\\\ \\hline & \\frac{23}{192} & 0&\\frac{125}{192} & 0 & -\\frac{27}{64} & \\frac{125}{192} \\end{array}$ \t| 6 \t| 5 \t| E \t|" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "| \t| Matrice de Butcher \t| s \t| $\\omega$ \t| Exp/Imp/Semi-Implicite \t|\n", "|--------------\t|-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------\t|---\t|----------\t|----\t|\n", "| Butcher \t| $ \\begin{array}{c|ccccccc} \n", " 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\\\ \n", " \\frac{1}{2} & \\frac{1}{2} & 0&0&0&0&0&0\\\\ \n", " \\frac{2}{3} & \\frac{2}{9} &\\frac{4}{9}&0&0&0&0&0\\\\ \n", " \\frac{1}{3} & \\frac{7}{36} &\\frac{8}{36}&-\\frac{3}{36}&0&0&0&0\\\\ \n", " \\frac{5}{6} & -\\frac{35}{144} &-\\frac{220}{144}&\\frac{105}{144}&\\frac{270}{144}&0&0&0\\\\ \n", " \\frac{1}{6} & -\\frac{1}{360} &-\\frac{110}{360}&-\\frac{45}{360}&\\frac{180}{360}&\\frac{36}{360}&0&0\\\\ \n", " 1 & -\\frac{41}{160} &\\frac{22}{13}&\\frac{43}{156}&-\\frac{118}{39}&\\frac{32}{195}&\\frac{80}{39}&0\\\\ \n", " \\hline \n", "\t\t\t\t\t\t\t\t& \\frac{13}{200} & 0 &\\frac{11}{40} & \\frac{11}{40} & \\frac{4}{25} & \\frac{4}{25} & \\frac{13}{200}\n", "\t\t\\end{array}$ \t| 7 \t| 6 \t| E \t|" ] } ], "metadata": { "hide_input": false, "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.6" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": true, "title_cell": "", "title_sidebar": "Contents", "toc_cell": true, "toc_position": {}, "toc_section_display": true, "toc_window_display": false }, "vscode": { "interpreter": { "hash": "916dbcbb3f70747c44a77c7bcd40155683ae19c65e1c03b4aa3499c5328201f1" } } }, "nbformat": 4, "nbformat_minor": 4 }