{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Arvo:400,700,400italic' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=PT+Mono' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Philosopher:400,700,400italic,700italic' rel='stylesheet' type='text/css'>\n",
       "\n",
       "<style>\n",
       "\n",
       "/*\n",
       "Largeur de la page\n",
       "*/\n",
       "div.container { \n",
       "    width:100% !important; /*!important; was necessary to override many styles*/\n",
       "}\n",
       "\n",
       "#notebook_panel { /* main background */\n",
       "    background: #888;\n",
       "    /*color: #f6f6f6;*/\n",
       "}\n",
       "\n",
       "div.cell { \n",
       "    width: 100%; /* set cell width   */\n",
       "    /* margin-top:1em; margin-bottom:1em; */\n",
       "}\n",
       "\n",
       "div #notebook { /* centre the content */\n",
       "    background: #fff; /* white background for content */\n",
       "    width: 95%;\n",
       "    margin: auto;\n",
       "    padding-left: 1em; /* padding-left: 20%; padding-right:20%; */\n",
       "}\n",
       "\n",
       "/*\n",
       "#notebook li { /* More space between bullet points */\n",
       "margin-top:0.8em;\n",
       "}\n",
       "*/\n",
       "\n",
       "/* draw border around running cells */\n",
       "div.cell.border-box-sizing.code_cell.running { \n",
       "    border: 3px solid #111;\n",
       "}\n",
       "\n",
       "/* Put a solid color box around each cell and its output, visually linking them together */\n",
       "div.cell.code_cell {\n",
       "    background-color: #f4f3e0; \n",
       "    border-radius: 10px; /* rounded borders */\n",
       "    padding: 1em;\n",
       "    margin-top: 1em;\n",
       "    /* overflow:auto; */\n",
       "}\n",
       "\n",
       "\n",
       "/*\n",
       "div.text_cell { width: 105ex }\n",
       "*/\n",
       "\n",
       "div.text_cell_render{\n",
       "    font-family: 'Arvo' sans-serif; /* 'Times New Roman','Palatino','Arial' */\n",
       "    line-height: 130%;\n",
       "    font-size: 115%;\n",
       "    width:100%;\n",
       "    margin-left:auto;\n",
       "    margin-right:auto; /* padding-left:3em; padding-right:3em; */\n",
       "}\n",
       "\n",
       "\n",
       "/* Formatting for header cells */\n",
       ".text_cell_render h1 {\n",
       "    font-family: 'Philosopher', sans-serif;\n",
       "    font-weight: 400;\n",
       "    font-size: 30pt;\n",
       "    line-height: 100%;\n",
       "    color: rgb(12,85,97);\n",
       "    margin-bottom: 0.1em;\n",
       "    margin-top: 0.1em;\n",
       "    display: block;\n",
       "    text-align:center;\n",
       "}   \n",
       ".text_cell_render h2 {\n",
       "    font-family: 'Alegreya Sans', sans-serif;\n",
       "    font-weight: 700;\n",
       "    font-size: 24pt;\n",
       "    line-height: 100%;\n",
       "    color: rgb(171,165,131);\n",
       "    margin-bottom: 0.1em; /* -0.4em; */\n",
       "    margin-top: 0.1em;\n",
       "    display: block;\n",
       "}   \n",
       "\n",
       ".text_cell_render h3 {\n",
       "    font-family: 'Philosopher', serif;\n",
       "    margin-top:12px;\n",
       "    margin-bottom: 3px;\n",
       "    font-style: italic;\n",
       "    color: rgb(95,92,72);\n",
       "}\n",
       "\n",
       ".text_cell_render h4 {\n",
       "    font-family: 'Philosopher', serif;\n",
       "}\n",
       "\n",
       ".text_cell_render h5 {\n",
       "    font-family: 'Alegreya Sans', sans-serif;\n",
       "    font-weight: 300;\n",
       "    font-size: 16pt;\n",
       "    color: grey;\n",
       "    font-style: italic;\n",
       "    margin-bottom: .1em;\n",
       "    margin-top: 0.1em;\n",
       "    display: block;\n",
       "}\n",
       "\n",
       ".text_cell_render h6 {\n",
       "    font-family: 'PT Mono', sans-serif;\n",
       "    font-weight: 300;\n",
       "    font-size: 10pt;\n",
       "    color: grey;\n",
       "    margin-bottom: 1px;\n",
       "    margin-top: 1px;\n",
       "}\n",
       "\n",
       ".CodeMirror{\n",
       "    font-family: \"PT Mono\";\n",
       "    font-size: 100%;\n",
       "}\n",
       "\n",
       "</style>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%reset -f\n",
    "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": "markdown",
   "metadata": {},
   "source": [
    "# M62 : Correction Examen 2019 Session 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Choix d'un schéma"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ">**Q1 [5 points]**  \n",
    "Considérons le problème de Cauchy\n",
    "$$\\begin{cases}\n",
    "y'(t)=y(t)-2e^{-t},\\\\\n",
    "y(0)=a\n",
    "\\end{cases}$$\n",
    "1. Calculer l'unique solution exacte à ce problème de Cauchy en fonction de $a$ pour $t>0$.\n",
    "2. Tracer la solution exacte sur l'intervalle $[0;4]$ pour $a=1$, $a=1.1$ et $a=0.9$.\n",
    "3. Soit $a=1$. Parmi les schémas EE (Euler Explicite), EI (Euler Implicite) et EM (Euler modifié), lequel faut-il choisir? Justifier la réponse en 5 lignes au plus."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Correction Q1**  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Il s'agit d'une EDO linéaire d'ordre 1, on peut donc calculer la solution à la main:\n",
    "$$\n",
    "a(t)y'(t)+b(t)y(t)=g(t)\n",
    "$$\n",
    "avec $a(t)=1$, $b(t)=-1$ et $g(t)=-2e^{-t}$ donc\n",
    "+ $A(t)=\\int \\frac{b(t)}{a(t)}\\mathrm{d}t=\\int -1\\mathrm{d}t= -t$,\n",
    "+ $B(t)=\\int \\frac{g(t)}{a(t)} e^{A(t)}\\mathrm{d}t\n",
    "       =\\int -2e^{-t} e^{-t}\\mathrm{d}t\n",
    "       = e^{-2t}$,\n",
    "\n",
    "d'où $y(t)=ce^{-A(t)}+B(t)e^{-A(t)}=c e^{t}+ e^{-t}$.\n",
    "\n",
    "En imposant la condition $a=y(0)$ on trouve l'unique solution du problème de Cauchy donné:\n",
    "$$\n",
    "\\boxed{y(t)=(a-1)e^{t}+e^{-t}.}\n",
    "$$\n",
    "\n",
    "Sinon, on peut utiliser le module de calcul symbolique mais **attention**, pour ne pas avoir des mauvaises intéractions entre le module `matplotlib` et le module `sympy`, on utilisera l'import de sympy avec un alias:\n",
    "```\n",
    "import sympy as sym\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "EDO:\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAALgAAAArCAYAAADLwUIqAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAHyUlEQVR4Ae2c63EUORSFxy4HYNgI1mTAIwJMBkAGu2SwlP/xz8VmAI7ACxkAEdiQwTqDBWfgPZ8sqVoadbd6xppRe3SretR6654+un2llr13c3OzyJV37949VtkzXUe6/lH8TW7dVq4hsA0EDqZ0KkL/UPknCpkVX6bUbWUbAttAYH9qpyL3sa3zdWrdVr4hsGkEJhNcA3yh60pEv970YFt/DYGpCKxCcCx4s95TkW7lt4LAoA8uK81i8r2uK13/6YLYLDRPdTVpCFSPQC/Bra/9SRqwqITgC4VuYdksePWPtg0QBJIuioh8qDzI/daRm8ISiP5Dac3/NnC0n9oRSBJcg8YtORSRP0YKNP87AqRF60agj+CvNezADbFWHZ/cuSl1a9ZG1xAQAksEt0TGRYmJDOkXyg+IT1qThkBNCIij33VhjJcJ3hmoWVh24ux/8yUTkr90DXTy221DYOsIiJcY58cKDX9TFpwFJFbazABGrMIvFbA9eElc8sI1cBst86s+/BiGesgtN9RGbXm5OuWWq02/guNhnWgMMX30bRO+Ut6ZwPtL4W+6znVhwd/btA+6Lyq2H2Zh/CZJ9Xuk8rxV/k5lzi1tl3Vf9VkJM4whh/8wxteKs1FyvjflNKEqbEQ0OAb5TOHb3A5VlsmIYvHOT24TVZSbs+4aOyRzz+yp7n8SV7q3qKVBVl/f1cepws/0VR3BNTB8qG8KnzDAWGw+Siy5Scoj/blC3KzZidVtlrpr7JD7g0Le9EZ0jxXF8PCsNrI5oX446fpAoeFAn4tiBrilH0AZcoHwsQAT6xAL9ag/13Pqc9Z9CXeRDOv9p54HHw0f6MoW1WPNR5tj8kZljRurEG7wFvcGrkYL/ksD7AVDeZD4WOGjlOZK/6X03xV6JVPlakxj7LpmqbvF/VC4eusJxkqH3Licj3Sfs56i2kqi9nlb4Nq+0sVYjvZXaqlQJQ0KIMZAYJYa/6pnGNR/3ZNXbfI90B0XZOgYNYQrLRi9C9vJicKrYi6KnUG8nuiUjfdg8ac4M/sPhV1Li/+25KupDKTG7TCzUiH7nNS/UBjvnFCfdoL+FN+YaEyMc6d0l87svKUEV2Oh/GChaTGChJxSReAJLk2XDyZjwg9vd1wWsD+nrWIEVwcn6oABY5XPdHnC2QGwrReDwsqbQQaicpD2q20L98QvZIKCt5F/FTAhskRt0V92edsoB87isXf7m4Xu3QGXuBdGkJv1kttZMd0oHbzBHVfCEN+WPYnLmgqZP7atYP1VhOB2sO5VARnjBSFpwYy2OmD54rI2ywR99bplqA+oWaKxBoBkVRooNCfdB9S4qyzesp+FiX/L6p5nwzEQLG2XA5A+Ph6ipPWkCME1JHwxN3j84dNomCgTp1Hkoa6hVxT1hvxv2sAHZ6JsS2ahu54PGH3TNQUrb3HHwFX7WGiwiN90pCOcVmVRiPAx8YviS+6pyV3j50CN5v/fiJGO1NYeRRQakiqEkADYdU+YwaRNUkZtUcfNft32ytgk6a14Fxlz0d2OM/mtYV0c1DY+8EOFKVfSGCnleau+bn9D9SG4IeVQoTXymL3xyhqffKF+nYXvNo97AZFTAjDU8xND91iB66gw9WknS1Qfi2LazqpwW2jMB6dU9bpP0De7qPDk+bIl6C237o3LqNDtkDn3NbvdVQserFoxsx6KxURmVhuSSmFmOv9AyJEUAAwYCmMJ/G/VAUjKx+1jwR2QcRtLcbVzpz54p4Pqde+M9U5uhSWLSvahg0Wl0nhW7i3Os8ElWRKeqa4xF3Sp3lBCaYKjDIQzYgF4qohTlpnuyE0ZyPqMm4R44qoOVpoDVikweO3GpE80Vzxpp3TXs2BCs6hkt8v52Q7kY6U5lyT1xZPnSXrOl0vXZlZYmuDMZE4lMnD2O3noENCdSjzXfVeIA1JKaIuzDlj9RQewuCzuRimrHPc1FN813dkBgeTm+UTAeIOj5/ZRF64lnGBLFyFe5JnV+KkepbNX6wYe+yOQAJjVePIzfrdsjfca987qXup5DFpwAY5PxUcaiLOpf7bJzGY2rzKjsZrUn6vssu5FnlmWBRfR2UrEqnqfV/f4Tcljq+uOVG27DwG4NFmiOkzC4LhmVsXKCu2y7iUexf5YowIcnxbx23O3UbO1Bql+2vhdBmwxxQuVsfYpv4rVH2t30/m7rPudYz1KcPXI9ly8l81A+tLJW0s0qdhZ4VOu+9I12J4t588FDxauPHOXdS/xaEZdFAGOG3KpMLCOirMg4pwBfm+ThkCVCASLTJEVl4OFDr4v23q4JSw0T3UtlO+24Mw+NHlKY1svdWyVKk0aAltFwFtwS17IGv+zTUgd/5UGX6Y+qU7Jz/xbBaZ1fj8QMD64iIpFhtyc3+7uXHCf+mebwWfz+wFF0+I+IuAWmbglfE1yn9CdrljvePeEvL50V6+FDYEqEHAE58x2QGRr1fHJg0PofelVaNMG0RCIENi3hMVFCYisOKRfKD8gvpKw3kG6bYPkJg2BqhBwFpxBdX1v4t7PFoE5xog1R3w6EfIUuDySmjQEqkEAC85HFay0J6klLduDl3ak3f8ilTq26k+L2fItaAhUgYDZJhShcVE4VHWhi8Po57ogPotP0jjja0iskInAZ3F2XRaKxwtTkps0BKpA4H/d3k++Dst9QAAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$\\displaystyle \\frac{d}{d t} u{\\left(t \\right)} = u{\\left(t \\right)} - 2 e^{- t}$"
      ],
      "text/plain": [
       "d                    -t\n",
       "──(u(t)) = u(t) - 2⋅ℯ  \n",
       "dt                     "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Solution générale:\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAJ8AAAAYCAYAAADkri+AAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAFkklEQVRoBe2a63HUSBCAbZcD4CCDIwMDEQAZABlgMjjK//zPxWXAEQGPDIAIKMgAMgCcge/7Bs0yK43eWq3WbFfNtjTqbnX39PT0jPbw6urqYA/L8cD5+flntHkM/rYcrabTJLXvcB980zl2rCQG5gYyfoIPx8paIn/ZvqMlKvkH6/QA279cY/vX7Jss8xHVf9Nal4qudNd4ACqm6RM6n9Ee0S5pH2iv6b8WgVhn3yTBh/B/cNY38DtwI0Bj9J+A/20k/AMf4hPrvYsuftxF95TtOx5rBAKdrbfAnYIJug80g++U9t/Y92+CH72svU5pD2lmogjv1ZkWMhX4eXwwET5BjlnvusKafaMyXzFIH8F3ct4qnjubH3K9tiRzb/99cDq4OTGz9qGPWfyMdkEz0Fb6ce1EMyDN3i+4n2zyIEuZb8F/gRcN6GgQveig5DNow7jn7Bub+VTgZYMSOtQs8SNDI5/81jpbB5xjtntLU18nRaXeou8dzeCTZuoMtcoKhS7W0BUdeG8jwBN8Dp5sYpRfWOilH/pAxb6jPtwZ2ictRqqgteAqe0QZBZ/8DvoS4CNK3KXdQaemQTdAL6FZy+QTGHAbGZ8KOWbeofL151J8WpgTUMW+wcGH812C2hzkLGzahMj/JKi2xR9sMQM7M59yXZkoJdXM4m9KfVPcuhLc5v3Wmu502/SY4p1zyqjYF5ZdDHWmaLTR+Zn7tZTNvbO9PDBmtcrSA60B51KqTJcnNxfyfwKXNyXyK2ftfdzPBuikjp1369A6YXRkL+A9+sOM9r1g1NfP6Q9BBjbbLqkEqdW10L8XytkXgg8pZzzUEWazV7RVMNBnUD4CPwan4BJVGQToDCh3tMp6AG6qDb5KQ2sF5PiuTrSJsC/wlfVOHofLOOAVW8qE3iPPYGlalits8Ki38v1sFnjBZtoz2tQ7ZkQOhzl1PS6cEGsNA8VlJQX7cs52Jpdpu/ClNPKbeVoBPWOQtNL2JHCSCE3lwS+KzC966QdLB3d2lV0/fdr3vnie+tGAtH8xMLeux1juhiA6RSdelLyhk8p9ktykhSXDmwzI1zagLmEO3jZBOw7wQVv9GnSE7gYt2A02e2mnUGdHzKjyubwLt2ieGVbKlvC0ww+8yo3vTjmiPbnJ2rYSbETXVLn02swXHakhOjBdcp219vVyEjLliTOey1poC+BaxgkfmH3VtxWwyyxpkIbJyr3YAY3ZMycjTEJoyvVujrZzH/JywXVQ6OIxzZD3bUTXOqOOkwfWRuVjkeBUDImZMSEPS27doGmEjlgFLderjJEIkd/BbwX462Z6E2/bTJdXHU+R3+Xb9D3o2rJ5Tp9Y1uSeLa1vNl3T4DNTlYPMei8EEE534/EGfFl4ywwgTw7W6kR41jJGwmDm67rcZWd6ImvopccsYVMFrs0W2CBdrvxoe6/2ucxWQL/QhgRzRdZEHbPqmgafLw71gobglBPQXVpchj2DioEniYF6z4sMrIIKHrObWSXnZAv0csBnxG2uC73M9mZ9P215eBztDS/l3izuZPKDf2p/eN7hp/IVBzn6xH7bkmBWXdPgc8v/CseogGdRBqPB4TdMC+XXtBS89/wuB8p6CZ8Z5QBcl1Ec2E1lNF/dCdDPz2Z+U9VWd6AxyCwJDMrBxyHw+n3YkkO/erQkeL91u3+p8vt3bl3H/rFAZ67Orn6b0X6FoWFDAvawdacBGywrDNyt21LoMnTDMes4HI18m7N56Aw2m8i/h2k9YNaOmXtayRNLG5X51IWZFg9QXaY7ATxmPZdla6mdhcIOJ5/lgzWy5cV3+uvKDB7vIXogrfliX18cinWY+gSSxyZDM2Zf/TZGT5A54QbXgxtTbEcEj8582llkAI8NWmc8NG5eLPA7Z8od8eVezZ4e+B8TAEgwJPzZ9wAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$\\displaystyle u{\\left(t \\right)} = C_{1} e^{t} + e^{- t}$"
      ],
      "text/plain": [
       "           t    -t\n",
       "u(t) = C₁⋅ℯ  + ℯ  "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Prise en compte de la condition initiale u(0)=a:\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAGoAAAAVCAYAAACuVXuDAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAD2ElEQVRoBe2Z61EUQRCAgSIA0AwwA5EIPDIAQ4AMpPjHPwozECJAyADMAMlAjACKDM7vG3e2Zmdv9+bkVqS4ruqbR/d0z/ZjXrc8Ho+Xcjg6Olqj74P91K9z+qI9fwtg5w2kvgdvqd/lGlbyDpg+0/e96r/J6Yv2MBaonKODTqhfgSZLDctpRkHcgXIBrlN/rLk6KpWwPcjbYMqvolPQKNmnPOgQ8aq6sYMZo303qaf2atgB2hUdj5S7kbAaK1W5RXnXJyTyw2PmHYLH4G46hvoO+JX+EXgCvlrADmbGGfgAup0YvNPgFgaTpobcUTWhq1IpNipU+JG2QhtA3yVolsnzqvc47GDmhMygbnCbVTPDyswj/uxfRobp23JSIk9nmr6tjTHhWVQLLTBTRmF0lzEjorHUdegy1b910BbdM1qg2FE4yWXM1HUPuyzQYya5Tz0Z0PcTIdeU+6XCqvkaWM7jHnTO7p1fKF8c5I5y4+s6jUQjFRkfgyinb2ksMhZynJNB4nJbBIwZwRhPV2Hppc+T1IjS02jXNxbJ/wdMBtabVE+9RzF5jfEJ7HJEPIWUZFOqI9SRvwbugT9axJ4O+DWq14XNHraaBJ+O1UkH1NP90XrRibYW9nwVbay9os2XVmnoIJcIS/eerlNa8HD28QyZDPCpKEQupfuaUS5oyJkgyikc5Leo+zTjV39vkDHGuXnZn2WO2uzJK0c6V+QZUOv0nVEeUh7rKDu9C7m0GYU3YDAw7RQ8HBR9AOONBCM4fABtS59G6gihPRS4KjSCDb3O20B0+esE+PzuosztFDI/gg5yzvrkekW5VDRi3IO8nE2C8PHwOXgabClzGtO86ejUIWLuEJ23BL3hQPv+R2CeBrQHN++pYc7BUclkNa7L1CRwSRF6swLB8pmhzwnp3uQ8vHzH7PbkVxJszzl/X4gemacZHiB3VOxvlQzy471h+2jo+14D6BuBwUmUtYIG0180kOV+MwanHkIqvUZg7Qj6DCyD76ZSv01f7siKNHjxttIQ9vsp2h5Sen48T2mtOh/o05CbXHjhpYwOUegFtLk/viLTyNKw0dBUe8FgchN26dAw56AZ5Zzt6zrVQhoG0OspVIgHKm3lN4XH60CZ8jOTo5SFAp0T97Mp4udDRue7UknV/MLbWjZmUl/GMkyTOT1Zd/HSN8wnLKSWWiDPqNaNuFRQHx8R5Z5hFpr6G7Tdy+4pX+RzDnMfGtZQELeVoCv/41CD+q7mJa73chhGL37mbgHsrpN+gd6f6kt7w1FqhaizPBR4MjmnvXAYhhgaKru76mj/lt1/A7cQaJAyB7uWAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$\\displaystyle \\left\\{ C_{1} : a - 1\\right\\}$"
      ],
      "text/plain": [
       "{C₁: a - 1}"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Solution du problème de Cauchy:\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAMkAAAAYCAYAAABQpNmPAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAF+0lEQVR4Ae2a63EVORCFrykHAEsGkIHZjWAhAyADIIOl+Od/lDcDIAIeGQAR8MgAMoAlA+/5hDSl0eg1b0HdrpI1I6mlc1pqqTW+J5eXl4ejHC0wxQLn5+efpHdP+dcp+i3r+NxOjk7S8lS1i02L6KrQ/af8pF2U05CF3K5M6+aodbTA4bZs8Pk3tUOP22InibzvhlLx2K1tt6fxfwWMNfZZgwd9auxHSneVfii9U3qp8l/eYVLcFnESdf6PDPVV+RvlWVEbvPRM+b/ZhjtVjuGyE8TqYde0tfrmPvJUeXHOqwE30jDkNjvcUofsKNdrjaV27DwH5Q8bsUkHQ5hGcekUd3wQZjacL0rcEXqisjVtfabBTP+9QX+Plx63WSeJnZj3ym/FbGPr2XHu6LkXiumd8r+Vc2TvLhZrksvuAD0AFusLFX1X+lOJSb2WsqXKF7W1+iMaeK38mvJmRfiwy0UFwEdqa9ZnjNtpRQe5JgB4lmmAMYlhmcxQ0EOf+LYFKXFpAaPBoIlkY7nHi54JdVkMOVna1t1Oq/E5wbiPjrqTqL1ZG8qf54DPqbOY7ozsY8Btbrh1v0ASgNxVBqeF1UMfI7cgJS4tYJyEYQVb3xSQDxbME+W9KKESJPPeytz7kAfcJjuJDE/8XjIOu0XuYof+fR/hHs+VXPaAtuSYS9qak+mm7Ma9ki9bg01wSeAb9zXgZsItkcSjIYwXfdJ77wjU+2uVPwiMwSkxuLipDY5BCEWfhFpcLNH/oDz8ooU+/fTG0/vWEuXigxB2uBCSsdi+KeH8dyOcVNykVNtanJg7Tgh4IqyLxyo3zqCc0Gr3MLmEE+BjJcbNOIk6eqJKjMDpwIWwW7Qqw3lYDCYG1rMTLox4XU/Ujsl4Z/u6rZwFmJIvqsCpiqJ+GKuqrdfZZ+mFuL3q7jHKxdWqD8bF0W/p2Zyeyt/qHX7Plczice0bzatsLS5wxdb83MTcM5QTp+M0j5WakC1xnloDuPiSBR1esimLXcrYbcK2vgFTen4b9NmhiyKca+5cSS4alzochE3EOIgFy3P0vmXrW8uKthY/5gLn52uPP+c4DuVNyNY4T8WaiXYG4X7wNLAEBgrLaPKHUm4HRS93H6EPFhqLcG/JcSHEuiobdaerBVvD7yA9+L2nD6tXk3W7eE3jyjY1tnaRAXz5aoZcV3qr90FobWoLf6SXigCw+UH1sc2vFAEsjjNHg5PELHTlTDoT2S0GlbGzUDbKQNJDx+1KekxKbnEmlTauYOPo8R/B76C22Df6f6SNedTY2ji+MId3x8lQ1VfMCbALoT2fjqeMtTjOHMFTr5LYPQwfIAIhd9J4zU2ohTPEBBLodYtLz+xOLBhf0CcMKIp0UztSTre0IzldMAy4gNmWh6GG+SKn+o6f66jhvNbWLvRumIqBthlO30nY+UNn4F5hFoIWBBf4V8p/GIg/QyV0YtK7j0gHZ+O4D/tnd6O8KOojuiMVFesagCHFhR5CjB0/yw1nDNvUjbxdqxpbw4HwaiDwVCqFzwO9lQo2xXnFI9GbZBnkTHV89XELm+/izkFQo/wvHiLSTYh02MFS/5ElDHH9R7rZrCjKxfJlk+gcSGU4PLb5aNENfnJjy7fK3KLG5jmpsTX3LxMFuI6YPyVO8RbmycHaFKd/kvB574UMAgC+j+M0GPZCZVziXir5wjtffWJCX8+kx+lzUJ6KO5mQNU8Ihq+RHBfCUOyCDViQtOUkcXZxl0gVbSfC42zvFjW/pWLOuGR390oPUdHW6CnhFBfS45MxwnsLc/QTjf4Kz6Y45/7AEUNO+hIjouzOTCj/qNpdhGMyl93BFwC0ZmvgChMn8tSLe4HxstV+uDWlZ3acqbsMpw36rcgcLq1wSOFozdbgJHT3w/cU9t3LZ50koNeO4P751LvT5JhJh1OEcIywpRmZwqUZ8Akgrdo6AbfJ4rknCaSI2cfG5bSfegIx5loyhctaWJbqt1VbL8Vv9X5mnyQgtLtV1Y/91JYL8Bvl1SfP6lbwBhjDxVNr8rF1WzdptAio/wEZK4Zd8EGGjAAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$\\displaystyle u{\\left(t \\right)} = \\left(a - 1\\right) e^{t} + e^{- t}$"
      ],
      "text/plain": [
       "                t    -t\n",
       "u(t) = (a - 1)⋅ℯ  + ℯ  "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import sympy as sym\n",
    "sym.init_printing()\n",
    "\n",
    "sym.var('t,a')\n",
    "u=sym.Function('u')\n",
    "f=u(t)-2*sym.exp(-t)\n",
    "edo=sym.Eq(sym.diff(u(t),t),f)\n",
    "print(\"EDO:\")\n",
    "display(edo)\n",
    "\n",
    "solgen=sym.dsolve(edo,u(t)).simplify()\n",
    "print(\"Solution générale:\")\n",
    "display(solgen)\n",
    "\n",
    "t0=0\n",
    "u0=a \n",
    "consts = sym.solve( [solgen.rhs.subs(t,t0)-u0 ], dict=True)[0]\n",
    "print(\"Prise en compte de la condition initiale u({})={}:\".format(t0,u0))\n",
    "display(consts)\n",
    "\n",
    "solpar=solgen.subs(consts).collect(t)\n",
    "print(\"Solution du problème de Cauchy:\")\n",
    "display(solpar)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dd3hc5Zn38e+ZrqnqXbIky3KRLRt3Q+gkIbSQQOglgcQQ0tksyS6bTXs3Pdm0DcRhSYAk9E3okNCbwb1XWZbVuzQzml6e948ZCzewLEuakXR/rmsuSTNnztw6Hus355zn3I+mlEIIIYRIN7pUFyCEEEIciwSUEEKItCQBJYQQIi1JQAkhhEhLElBCCCHSkgSUEEKItCQBJYQQIi1JQAkhhEhLhuEspGlaPnAaUAwEgG3AOqVUfAxrE0IIMYVpH9RJQtO0s4FvAtnARqALsAA1wHTgMeDnSinP2JcqhBBiKjleQP0U+I1SqukYjxmAiwC9UurxEb24pmUC9wBzAQXcpJRaPZJ1CSGEmFw+MKCGFtK0SqXU/uPdd8Ivrmn3AW8ope7RNM0EWJVSAyezTiGEEJPDcANqg1Jq4RH3rVdKLRrxC2uaE9gMVCnpWCuEEOIIHzhIQtO0WUAt4NI07ZOHPOQkcS7qZFQB3cAfNU2bD6wHvqKU8n3AcyTIhBBi4tFG8qTjDTOfSeI8UyZw8SG3hcDnRvKChzAk13OXUuoUwEdiQMZhNE1bqWnaOk3T1tXW1p7kSwohhJgohnuIb8VoD17QNK0QeEcpVZH8+XTgm0qpC9/vOYsXL1br1q0bzTKEEEKMvdHfg9I07T80Tct6v3DSNO0cTdMuGskLK6U6gGZN02Ym7zoX2DGSdQkhhJh8jneh7lbgaU3TgsAGEueMLMAMYAHwIvCDk3j9LwF/SY7gawA+cxLrEkIIMYkcL6AuV0qdpmnaHSQu0i0CPMCfgZVKqcDJvLhSahOw+GTWIYQQqRaJRGhpaSEYDKa6lJSyWCyUlpZiNBpHZX3HC6hFmqZNA64Fzj7isQwSbY+EEGJKa2lpweFwUFFRgaaN6HTLhKeUore3l5aWFiorK0dlnccLqLuB50kMCT90dIJGYsh31ahUIYQQE1gwGJzS4QSgaRo5OTl0d3eP2jo/cJCEUurXSqnZwL1KqapDbpVKKQknIYRImsrhdNBob4NhTbehlPr8qL6qEEIIcRwyH5QQQkxhr7/+OgsXLsRgMPDYY4+lupzDSEAJIcQUVl5ezp/+9CeuueaaUV/3znYPL+7oHPHzhzVhoRBCiPR36aWX0tzcTDAY5Ctf+QorV6487nMqKioA0OlGd3+l3xdm5QPriMYU580pGNE6JKCEEGKUXfn70Z3W7uFbVgxruXvvvZfs7GwCgQBLlizhsssu47bbbmP37t1HLXv77bdzww03jGqdB0Vjcb744AY63SEevmX5iNcjASWEEJPEr3/9a/72t78B0NzczN69e3n44YfHvY4fPbeLt+p7+cnldZxSnjXi9UhACSHEKBvuHs9oevXVV3nxxRdZvXo1VquVs846i2AwyJVXXjmue1D+cJR73mzi06dWcMXispNalwSUEEJMAm63m6ysLKxWK7t27eKdd94BGNc9KH84Sr8/wrLKbO68cPZJr09G8QkhxCRw/vnnE41Gqaur41vf+hbLlw/v3M/atWspLS3l0Ucf5ZZbbmGk8+5FYnEO9PrRaxq/u3YhRv3Jx4vsQQkhxCRgNpt57rnnTvh5S5YsoaWl5aReOx5XNPf5Mep1ZNtM5NjNJ7W+g2QPSgghxIgppWjpDzAYipJrN2EyjF6sSEAJIYQYsU5PiIFAmEKnhUyraVTXLQElhBBiRPp9Ybq8QbKtJvIco3NY71ASUEIIIU7YYChKy0AAu9lAcVbGmHRzl4ASQghxQkKRGAd6fZj0OsqzrejGaKoRCSghhBDDFo3FaXcH0dCoyLViGIXh5O9HAkoIIaaw4U63ceedd1JWVobL6WQwFKUix4rZoB/T2iSghBBiChvudBsXXnQRjz33Kgooz7ZiNY/9ZbQSUEIIMUlceumlLFq0iNraWlatWjWs51RUVFBXV/eB020opSibOZ+MzBx0GjgzjKNV8geSThJCCDHKJtt0Gx2eIP3+xLVO40kCSgghJomxmG6jZzBEtzdEjt08Jtc6fRAJKCGEGGWTZbqN3sEQbQMBcu1milyWMbnW6YNIQAkhxCQw2tNt9PnCtA4EcFqMFKYgnEAGSQghxKQwmtNtDPjDtPb7uer8M4YuxL3jjjsoLS3F7/dTWlrKd77znbH9hQBNKTXmLzJaFi9erNatW5fqMoQQ4jA7d+5k9uyTn6AvHXgCEQ70+rGa9FTk2tDrTmzP6X22xYh2v2QPSgghBACDwQgH+vxYjDoqcq0nHE6jTQJKCCEEg6EoXd4QZoOOylwb+g+4Lmq8pL4CIYQQKTUYjNLY4yMSU1Tm2sa0v96JkFF8QggxhQ0GIzT2+jHpdVTm2TCmSTiBBJQQQkxZQ+GUPKyXTuEEElBCCDEleYMRur2htA0nSINzUJqm6TVN26hp2tOprkUIIaYCtz9MY6+fWFxR4jBw3TVXU11dzbJly2hsbDzmcx5++GHq6uqora3ljjvuGJc6Ux5QwFeAnakuQgghpoI+X5imPj8ZRj2VeTbu+9MfycrKor6+nq997Wt84xvfOOo5vb29/Ou//isvvfQS27dvp7Ozk5deemnMa01pQGmaVgpcCNyTyjqEEGIyON50Gz2DIVr6/djMhsRoPZ2OJ554ghtvvBGAyy+/nJdeeokjGzg0NDRQU1NDXl4eAOeddx6PP/74mP8+qT4H9UvgDsCR4jqEEGLUpON0G9G4IhqLo9M0TAbdULPY1tZWysrKADAYDLhcLnp7e8nNzR1ab3V1Nbt27aKxsZHS0lL+/ve/Ew6HR/V3PJaUBZSmaRcBXUqp9ZqmnfUBy60EVkJi5kchhBDHdqzpNh566CHa3UG8wSg2k57irAx0hzR+PVa7uyMbw2ZlZXHXXXdx5ZVXotPpOPXUU2loaBjbX4bU7kGdBlyiadoFgAVwapr2Z6XUdYcupJRaBayCRC++8S9TCCFOTLpMt+EPBLj4E5dTv3cPep122Ei9g3tQpaWlNDc3U1paSjQaxe12k52dfdT6L774Yi6++GIAVq1ahV6vH/PfKWUBpZT6N+DfAJJ7UF8/MpyEEEIMz7Gm22h3B/nBb/+XIpeFXLv5mFNmXHLJJdx3332sWLGCxx57jHPOOeeYy3V1dZGfn09/fz+/+93veOSRR8b8d0r1OSghhBCj4Pzzz+fuu++mrq6OGTU1zF+4hFAkRlm2lSyr6X2fd/PNN3P99ddTXV1NdnY2Dz300NBjCxYsYNOmTQB85StfYfPmzQD853/+JzU1NWP7CyHTbQghxElLp+k2AuEYHZ4g/lCU8hwrDotxXF9/NKfbkD0oIYSYJDyBCE19fvQ6jao8Gxmmif0nfmJXL4QQAqUUvb4w7QMBLEY9FTk2jIZ06MNwciSghBBiFCiljjm4YDxet20gSK8vhNNipCw7dRMNjvYpo4kfsUIIkWIWi4Xe3t5R/wN9PNF4nOb+AL2+EHkOM9NyUhtOvb29WCyWUVun7EEJIcRJKi0tpaWlhe7u7nF7zUgsTu9gopuDw2JgwGtgoG3cXv6YLBYLpaWlo7Y+CSghhDhJRqORysrKcXu957e18y+PbCbDZODu6xayuOLoC2snAwkoIYSYIOJxxX+/uIffvFzP/LJMfn/dIgpdo3dILd1IQAkhxAQw4A/z/57ZyRObWrlicSnf+/hcLMaxbzeUShJQQgiR5ra0DPD5P2+gyxvkh5+Yx2WLSlMyYnC8SUAJIUSaUkrx53eb+P5TO8hzmHn01lNZUJaZ6rLGjQSUEEKkIU8wwn/+fRtrGvs4tTqH/75iAVm29++pNxlJQAkhRJrZ2NTPlx/aSNtAkH+/YDafObUCXYqub0olCSghhEgT8bji96838PN/7KbAaeGRW1awaFpWqstKGQkoIYRIA52eID96bhd/29jKBfMK+eEn63BljG8n8nQjASWEECn21OY2/uPv28ixmfjBJ+Zy9dLyKTFK73gkoIQQIkXc/gjfemIbT25uY35ZJr+4Yj7T8+ypLittSEAJIUQKvLm3h68/upmewRC3f7iG286ajkEv/bsPJQElhBDjyBuM8LMXdvPUlnayrEZW3XAqdaVT59qmEyEBJYQQ4+SV3V3c+X9bafcE+dp5Naw8o2rStys6GRJQQggxxvp9Yb7/9A7+b2MrM/LtPP75U1lYPnWHjw+XBJQQQowRpRTPbu3g209uY8Af4cvnVPOFc6oxG2SvaTgkoIQQYgw09fr5zye34QtFKXJl8MDNy5hd5Ex1WROKBJQQQoyiUDTGqtca+O0r9Rh0Grd/uIYbT62QEXojIAElhBCj5O36Hv7jiW00dPu4cF4R37pozqSeUHCsSUAJIcRJaur184Nnd9LY6yMeV/zpM0s4a2Z+qsua8CSghBBihAZDUX73Sj33vLkfvabxhbOn89kPVWIxyZ/W0SBbUQghTlA8rvjbxlZ+/PwuurwhPnlKCXecP0sO540yCSghhDgBb+/r4b63GnlhRyfzyzK5+/pFck3TGJGAEkKIYdjR5uHHz+/itT3dFLss/PaaU7hgbtGUnEhwvEhACSHEB2ju8/Pf/9zD3za14rQY+fcLZnHDigppUTQOJKCEEOIYur0hHlzTxG9frkfT4JYzpvP5M6fjsk7tSQTHkwSUEEIcos8X5vev7+P+tw+QYdJx+aJSvnRuNUWujFSXNuVIQAkhBIlg+vM7B/j9a/vwR2J8fH4xXz53BlUygWDKSEAJIaa0Lk+QP7zRwF/ebWJuiZOzZ+Xz1fNmUJ3vSHVpU54ElBBiSmrp97Pq9QYeWttMNBbnkvnF3HZ2NTUFEkzpImUBpWlaGXA/UAjEgVVKqV+lqh4hxNSwvc3NPW/sp3XAz8amAS5bWMqtZ06nIteW6tLEEVK5BxUF/kUptUHTNAewXtO0fyqldqSwJiHEJKSU4o29Pax6vYE363uwmvRcs6ycX155CsWZMvghXaUsoJRS7UB78nuvpmk7gRJAAkoIMSqCkRhPbmrl4XUtrD/QT77DzB3nz+TapdNkuPgEkBbnoDRNqwBOAd49xmMrgZUA5eXl41qXEGJiancH+PM7B3hwTTN9vjAX1RXxk8vr+PiCYpnNdgJJeUBpmmYHHge+qpTyHPm4UmoVsApg8eLFapzLE0JMEEop1h3o509vNfL89g7iSnHe7AI+c2oFK6bnoGnSkmiiSWlAaZpmJBFOf1FK/V8qaxFCTExuf4THN7Tw4JomegfDaBrc/KFKrl8+jbJsa6rLEychlaP4NOB/gZ1KqV+kqg4hxMSjlGJD0wB/efcAz2xpJxSNM78sk8+dXsVFdUVYzSk/OCRGQSr/FU8Drge2apq2KXnfvyulnk1hTUKINNbhDvJ/G1t4bH0LrgwjezsHuXxRKdcsK6e22JXq8sQoS+UovjcBOSgshPhAwUiMf+7o5LH1Lbyxt5u4gqUV2XxqcQkXzCvGJntLk5b8ywoh0k4srnh3fy9PbW5ja6ubba0eil0WvnB2NZctLJWLaqcICSghRFpQSrGt1cMTm1p5aksbnZ4QVpOe65dP498+NpsVVTkyOeAUIwElhEgZpRTb2zw8t62dvV2D/GN7J0a9xlkz87lkfjHnzS4gwyTXLU1VElBCiHGllGJLi5tnt7Xz3NYOmvr86HUaH19QzI8+OY+PzS2SLg8CkIASQoyDYCTG6oZeXtrZyfoD/exs92LQaZxWncsXzp7Oh+cUkm0zpbpMkWYkoIQQY6J3MMQru7t5cUcnr+/txh+OkWHUc0ZNLreeMZ2zZubLnpL4QBJQQohREYsrtra6eW13N9vb3PxzZydKQYHTzCdOKeG8OQWsqMrBYpRzSmJ4JKCEECPW6Q7w+t4eXt/bwxt7uxnwR9A0mFfi4mvn1XD2zHzmljilD54YEQkoIcSw9fvCvLu/l7f39bJ6Xy9ZNhNr9veRazdz7qwCzqjJ5fQZeXI+SYwKCSghxPtyByKsP9DH2/WJUNrZ4UEpsJr0LKnI5rzZ+Xz74jnMLnTKNUpi1ElACSGGtLsDrG3sZ11jH2v297G708spZZlsa/OwqDyL28+rYcX0HOpKMzEZdKkuV0xyElBCTFHhaJwd7R42NfWzpcXN2sY+mvsDANhMehZOy+Jjc4tYMT2butJMGdwgxt2ECqgOT5AnN7dRV+JiWo5VTrwKMUzxuOJAn58dbW7WHxhgU3M/29o8hKNxIDHS7qyafGYWOlhSkc3sIgcGvewhidTSlJo4k9Raimaowht/CYDDYmBeiSt5c1Jbksm0bKscBxdTXiyuaOzxsa3NzbZWN1tb3Wxv9eANRXGY9UTjMLfEySnlWSwoy+SU8kyKXBmpLltMbiP6wzyhAqp6Tp36n0f/Qb8vzNbkf7yd7R7Ks63s6/ZhNemZWehgdpGTOUVOZhc5mVlgx26RiwHF5OQNRtjb6WV7u5ed7R52tHnY3eEl32HmQJ8fk0HH7EIHc0tc1JW6mFvioibfjtEgh+vEuJr8AWUumqGKbvwlmVYji8qzWDgti/mlLiwmHQ1dfna0exL/Sds9eINRAJZVZtM6EKCmwMGMAjs1+Q5qChxMz7dhNU2oI5xiCvOFouzrHqSxx8/2djd7Orzs6RykdSDA0sps1uzvw2kxMDv5wWxOkZPaEic1BQ6McqhOpN7kD6ia2vnqu/c+xcbmftYf6Gdft48lFVlsaBqgpsDBgrJMFpS5mF+aidVkYHenl+Y+Pxua+qnvGmRf9yCRWOL3nZFvxx+OMT3fTlWujel5NqryEt8XuixyfkuMu3hc0e4J0tjj40Cvjz2difdsfdcg7e4gAEsqstjUPMD0PDszCxMftmqLnVTn2ynJzJD3rUhXkz+gzEUzVPXnfsPiimxWTM9hTpGTaEyxqbmfjc0DbG4ewJPcc7KZ9FxQV4TLYmRe8tBGWWYGTf0B9nZ6aXcH2dQ8QEPPIA3dPvzhGAA5NhP+cIxpOVam5VipyLExLcdGRa6VkswMijMz5BOpGLFQNEZrf4Dm/gBNfX76fWG2tbpp7PVxoNdPKDloYUa+nbaBANPz7UzPs1M99DXxfpT3oJhgJn9Azaidr676r7+wuqGX+q5BagrstA0EWTQti6WV2SytyMJhMbKzw8PmZje9gyFe2NE5NFLJZtJTW+xifpmLmYWJwyDV+XaMeo0OT5CGbh+t/QF2d3o50OujsddPU5+fcDROXYmLLa1udBoUOi2UZlspy7Iys8BOptVEcWYGxZkWijMzZDjuFOYPR2kbCNI2EBi6dXqD7O/x09znp8MT5ND/ctOyrRgNOipybFTmWqnItVGZY6Mi10aR7MmLyWPyB9TixYvVunXrAOjyBtlwoJ+36nuHLigEsJv01CVHJi2uyKauxEWXN5QcyZQYWBGOxtnW5gHAqNeYnmdnTpGTOcWJ0KopdJBrNwOJwy4dniBNfYmwaukP0NLnp7k/8b1Rr6Opzz9Uo92sx6DXUei0UOC0UORKfC3PtpJpNZLnMJPnMJNjM8uFjhNIOBqnZzBEtzdEvz9M60CATneQDk+QDk8o+X2AQCQ+9IEIQKdBgcPMtFw7xZkWyrKslGdbKctOfM13mGXkqZgKplZAHanfF2ZtYx/b2zy8uruLbW0eYnGFpsHSimyq8uwsmpbFwvJMSjMtNPUH2NHuZUdbYmDFznYPXd4QyyqzeXd/Hzk2EzMK7MwscDCjwJE43p/vOGp6gHA0TqcnSOtAgHZ3gB5vmAN9PjrcITo8ATrcIXp9IZZWJNZ7qCyrkXyHmVyHmSyriRybiWybmWx74vscmwlnhpEsq4lMq1H2zEZRMBLDHYgw4A/T6wvT74vQ5wvRl/xq1OvY1uamZzBMtzeEOxAZeu78UhebWxJ707l2M4XJDyGFTgtl2RnkOyyUZCUOBxc4zHI9kRBTPaCO5A9H2dQ8wLrGftoGAjyztX1oZF+Ry8LsImdyUEUm88sycWUY6RkMsbfTy852L3s6vezu9LKnw4sveX4qz25G0zjkfMB75wbyHSZ0umP/IYrE4nR7g3R7w3R5E5/Cu70hugeDePwR2txBen1hegdDQ+fQACpyrDT2vrd3ZjboyLKaKM+xEo8rnBlGHBYDDosBp8VIgcuCDrCaDNjMBuxmA1azHrvZQIZRj8WoJ8OkJ8OoRz8BP7XH4opgJIY/HCMYiREIxxgMR/GFovhCscTXcJRYXNHlDeENRvAGo3gCEXLsJra0uBnwR3AHIkPnehaVZ7G+qf+w13FlGFlWmU2/P0yew0yu3UyePfFBIs9uJt9hJs+Z+F7CR4hhkYD6IPG4or57kA3J0X+v7u6ivntw6HxAVZ6N06pzqc6zU1fqYnaRE4tRj1KKNneQPR1emvp8bG31JEYEdg3iDSXCxGrUoWkalXk2KnPtVObaqMq1UZmbGFzhyhh+Z+dILE6/L/mp3h9mwB9J3AJh3P4I/f4wwUjicJPnkD/AnmCUuSUuNjcPHL3dpmWx7sDhf4RNeh0Wo44FZZk09PgwG3SYDHpMBh1mvY5pOVa6B0MYdBp6nYZBr8Og0zDodOQ7zQz4w4CGTgNNA52mkW0z0ecLD21ThUIpyLaa6B4MEYurxE0lvpoNetyBCJFYfOhm0Gm4A1FC0RihaJxQNE5JpoUd7d7DDp0BQ3u7RzqlPJOtLe5keCdCfG6xC3cggivDSKbViMtqxJWR2IO1m43k2E1De6oyAEGIUScBdaI8wQhbW9xsah5gY9MAkVic1/Z0A2DQacwsdFBX6qKuNJN5JU5qCpxD542UUnR7Q9R3DXKgz8/uDi/7e3zs7/HR0u8nntys80qctPQHKM+xUZFjZVq2NTkq0EZpVgZ59tE5B6GUIhCO4Qsn9iQGQ1H8ye/DsTieQCSx1xGJEQjHCUQSeyEWo45OT4hwMgzCsTjhaIw8h5mmXj/RZKhEYnFicUU0rphV4GBLqzv5OyriKvH680pdbGlxo8HQyX0NqCt1sb3Ng0GnoUsGnl7TyHea8QSiGA06jDoNo16H1aRH0zTMBl3iZkzsNcaUIsOox5rcA7QY9TgzjJj0OmxmA7bknqLVZMBq0g+tRwiRFiSgTpZSinZ3kC0tA2xpSQyo2NLiBhTuQBSjPhFa80pc1Ba7mFfsYmaR46hzQ6FojOY+Pw3dPrq8Iba3eWjq89HY46fdHSCuoKbAzp7OQcwGHaVZGZQlRwVW5dkocFqGhrTn2k3yh1YIMdFJQI0FpRQH+nxsbfGwrS3R02xrqxt3IML8Uhfb2jxMz7MlrtwvdjGn2MnsIgfZNvMx1xeKxmjpD9Dclxh23NwfoKk3MSqwuc9PeY6Nba3uoeXNBh0lmRmUZGUwu8iJzWSgKDMxOrDIlRjaLh0xhBBpTgJqvCilaOkPsKfTy6bmAXa0JdortbuD6LXE4a1cu5lZRY6h1jOzCh1U5dqOe1LdHQjT2p8YFdja7098HQjQ2h/AlWHk9b09Rz3HaTFw+ow8BkNRCpxmCpJD3A+OLCtwmsmxmyfkwAghxKQgAZVqvYMhdnV42N7mYVe7lx3tiQEV0bii0GWhzxemOs/OrMLEsPXaYhfV+XYKnOZhH8YLRWN0ukO0uQN0uIO0uQO0DwTRabCxeYAOd5CewdDQObCDlldl09DtI8+RHIWWvJVlW3EkBwkcHK3mzDDIYUUhxGiSgEpH4WicvZ1ednckhq3v6vCyq8NDpyfEomlZrD/Qj9NiGOqrNjc511VNwXsXC5+oaCxOry9MhztIpydIpzdEIBRlX7ePLm+Q7sEQXZ4QPYOJGtY2Hj7Cz6jXyLGZqSt1EYrGk9dnmci2m8i1mcm2Gcm2J67dyraacFgMcrGpEOKDSEBNJP2+MHu6ksHVkbzuqsPLzELHUGBk20zMyLezvCqHHLuJ6nw7M/IdozZwIhZX9A2G6fWH6PGG6RkMJW+J7416Hdvb3PQOhun1hQhGEsO8C51mOjyhofXodRpZViO5djMOi4FMq4nMDCNZNhOugxcaZyTud2Ykhne7rEbsJgk2IaYICaiJTqlEW6W9nYPs6fQmvnZ5sZkMvFn/3rmnkszEAImDFwpX59upzrNTnJkxpn/w/eEovYOJ67MS3RfC9PkS12r1+cP4Q1E6PMHDrt06GGr5DjNd3tBh69Np4MwwUlfiot8fwZmRuODYaTHiTAaazaQfupbp4FeXxYjNkrgQWdpFCTEhSEBNVkopOj2Ja672dnnpcAfZ2DxAfdcgfb4wkGjntKV1gKpcO9OTgTU938aMfDvTcmwpa5MUjMQY8Efo94VwB6O4AxHcyW4OB29GvUZDj2/ogmNPIIInECEUi3Ost+fSiizWJPcyTQYd2VYjZmPiOiib2YAj+bXIZSGu1HtdNUyJ66WcFiNmgw6r2YDNpB/6mmFMXKgs59+EGHUSUFNRny9MfdcgTX0+drZ72dedmEOopT+AUjCr0MHuTi8lmRlMz7NTlWdjbrGTouTP+Y7hD9AYb+FofKhdUeKWCLBwNE6fL4Q3mLggORiJ0e+PDF2gPBhKtD9yWIzUdw0SiMQOW29d8oLiIxW7LHR6Q1iNeqxmPVaTgSyrEZ2mDbWIykheBFzkshCMxIcuGraY9GRmGNE0sBgSy1mMusRjRj1mg+6wr9KtQkwxElDiPcFIjIZuH409PvYmJ2vc152Y+2pWkYONTYmWSHazgcpcG8urchLf5yXaNFXlTZ4Zh2NxhT+c7NcXfq/DxsH7Dn4NxxKB6Asl+vz5IzGUUvT7w4mfkz0A/eEYswodrG7oHZoAE47deqkk00LrQPComhaUZbKn03tYcJkNemYWOmh3BzAfbDuVvNnNBtAYut+k12Ey6Mi0GonF1dDPB59j0CdaVhkNOox6DZNej8mQ6NZh1Osw6DXMej3G5H0GnZa2H1TEpDDxAkrTtPOBXwF64B6l1I8+aHkJqJMXjyva3QEakm2ZGrp97OsexKDTeHVP92GH1AqdFr83G4kAABwhSURBVJZVZQ+FWFWejapcO6VZGdIkNSkSiw+1kApGEt8Hk98HIjEi0UT7qVAkTjCa/BqJoddp9PnCyX6Dib6DwUiMbJuJxh4/4Vjy/kii/ZQGeIJRQpFEkB4MxkMPdx7kyjAe1n29KtdGQ4/vsGWObEQMsKDMxZ7OxHvBZEgE2ewiJ409Pgz6RB/GkswM+v1hDPpEsJVlJXo2GpOPF7ksycff692YaTUQjMSTfR11yXVp2MwGorE4+uSyep2G1agnphQGnS7RA1KXaHuFBnrtYF/IxHo1LTFA52DrLF1yeZ2WWObQ+zTAoE+sU5dcj05DQnn8TKyA0jRND+wBPgy0AGuBq5VSO97vORJQYysYidHY62N/t4+GZHh1eoNsbXEf9gdvRVUOnd7gUEPcqrxEg9zKXCv5DplkbzzE4yoRYpE4kXhiDqrwUC/FxNdI8muiC/x7DXkPNt31hWNED94XU1gMOvr9YSIxNbSs02KkwxMkGlNE44f/HInHKXJl0NCduNYvGotTnW9nW6sn8XM8TiymqCm0s6XFfdjeJsDMQge7O7yH3be0Ips1jYfvhR7rkOyMfDt7uwYPuy/DqCMQObyh8EGVuTb2HxHSZoNGJKbQJ/ceDwbgKeWZbGt1H3W/QZ94X+u0RLjpNA1XhhFfOIYuGZZa8jG9ppFjTzRVPhiIOTYTA4HI0HM1TaMk00KXN5T8OXF/vsNMny88tK6DgRtXHPLcxAeRwVAULdm0WadL3J+ZYcQbjA4tp2mJgNYd+n1yMJXNpCcYjScfP3x5i1FPJBZP3AeQrCfDqCcUjSfWR+LBxAcADaVILv/euipyrNSVZo7oj0Iqj+EsBeqVUg0AmqY9BHwceN+AEmPLYtQzq9DJrELnUY/1+cLs70kcIuwZDLGlxc3+Hh9v7O0ZmrqiIsdKtzd0dFf3vESQOS3Go9YrRkan07Do9BNqjjClEo2Fo/F4IuBi8aGf43GGAi2a7HYfjb3X/T4aixONK+LJhsUKRST6Xmf82CGPxdR738dVYj16HQQj8aHHYnEANdQM+eD9cZXozFKRYyOuEs+PJe+PJ9cXUxBXCqUUZr2OwXBs6LF48rG4UhwcUBuNxwnHIBCJMRCIoA5Zp8mgsa97cOh5SsH0PDu7Oz3E4+9tM6NBIxCOc7A5c1wpZuTb2dXhRR3y3LhSnFKWyfqm/sPW+X5qi51sT07eeqRZhQ52HfEBAo49OwIce6QuwNVLy6krzRzGO+RoqdyDuhw4Xyn12eTP1wPLlFJfPGK5lcBKAKPRuKiurm7caz1Z3d3d5OXlpbqMEzbcug9+ag9F44lDUId0RofEJ7O4UkOHa0wG/SHdyvWY9InDNeNZc7qRusfXVK1bHfGNQh39uDpi2SMWOPQ56hgLHvk8vaaxedOG7UqpuSdabyr3oI71J+mobaKUWgWsArDZbBPyEN/ixYuZinWHojEO9Ppp7EkcMtzfnTjv5Q1F2NnuJQyESRx6KM2ysqwyG7vFQFWenenJc14FzhM7ZDhVt3WqSN3ja6LWrWna0SOFhiGVAdUClB3ycynQlqJaxBgwG/TUFCRaOB3JHYjQODRQY5CGHh/N/X42N7uHhoXXFidO0FcmB2ccHCZflWejMseG1Tw5RhkKIY4tlf/D1wIzNE2rBFqBq4BrUliPGEeuDCPzyzKZX3b4semD3TQaun20DQTY0e5hX7ePDU39PLWlbeh4+vxSFz2DYabn25melxioMT3PRtxkRyklAzWEmARSFlBKqaimaV8EXiAxzPxepdT2D3pObm7uuNQ22lauXJnqEkYkFXVrmkaRK4MiV8ZRjwUjsaGh8S39fnYmw+vhxj784eTFuGf8K6f/5BVy7ebD2kBNz7dRnm1L2ylH5D0yvqTucbdqJE+SC3XFhHdwr6v+4AXJXT7quwap7x6kOzmqaFllNhubB5ieZ2dGfvJWkAiwipzjz9MlhDgpE+s6qJGQgBInyh2IUN81SHOfjx3tXvZ2etnblWgFZTZohKKJLgxVeTZmFDioybdTU+gY6mGYrntcQkwwElBCDJc/HGVf1yB7kp3jE7dBWgcCAEzLyaDDHWJGgZ2ZBYkZkWsKHcwqsJN/giMLhRAjC6i0PK6hadr5mqbt1jStXtO0bx75uFKKL3/5y1RXV1NXV8eGDRtSUeZhnn/+eWbOnEl1dTU/+tHRHZteffVVXC4XCxYsYMGCBXzve99LQZVHu+mmm8jPz2fu3GNfopCO2xqOX/fxtrfVZGBeaSaXLSrl3y6YzR8/s5S3vnkO2777Uf7+hdP42nk1XL98GllWE2/s7ea/nt3Jjfeu4dr/fZeF3/8nV696h+8+tZ1H1jWztcVN8IiGtMfS3NzM2WefzezZs6mtreVXv/rVUcuk4/YeTt3p+P4OBoMsXbqU+fPnU1tby7e//e2jlknH7T2cutNxewPEYjFOOeUULrrooqMe0xJ+nfy7vkXTtIXHXaFKXhGdLjcSAyb2AVWACdgMzFFKsWjRIqWUUs8884w6//zzVTweV6tXr1ZLly5VqRSNRlVVVZXat2+fCoVCqq6uTm3fvv2wZV555RV14YUXpqjC9/faa6+p9evXq9ra2mM+nm7b+qDj1T3a27tvMKRW7+tRD605oL75+GZ1yW/fVLP+4zk17RtPq2nfeFpVfvNp9am73lZffnCD+v1r9eqNPd2qdzB02Dra2trU+vXrlVJKeTweNWPGjKPeJ+m4vYdTdzq+v+PxuPJ6vUoppcLhsFq6dKlavXr1Ycuk4/YeTt3puL2VUurnP/+5uvrqq49VG8AFwHMk9qaWA++q4+RBOl5IctwWSE888QQ33HADmqaxfPlyBgYGaG9vp6ioKCUFr1mzhurqaqqqqgC46qqreOKJJ5gzZ05K6jkRZ5xxBo2Nje/7eLpt64OOV/doy7KZWF6Vw/KqHK5cUg4kuqQ39SVGEx68rd3fxxOb3rucr8hlYU6Rk6WV2VTm2pg3fTZKKRwOB7Nnz6a1tfWw90k6bu+ioqKh13+/utORpmnY7XYAIpEIkUjkqEOz6bi9h1N3OmppaeGZZ57hzjvv5Be/+MWxFvk4cH8yrd7RNC1T07QipVT7+60zHQOqBGg+5OcWYNmhC7S2tlJW9t41vqWlpbS2tqbsTXWset59992jllu9ejXz58+nuLiYn/3sZ9TW1o5nmSOSbtv6RIz19tbrtGSTXBsXzHtve/T5wuxs97C9zc2ONg/b2zxsaXUPjSjMtZuYnmViU6wUX+Z02gYCFLkS57XSfXs3NjayceNGli1bdtRj6fj+jsViLFq0iPr6er7whS8cVXe6bu/j1Q3pt72/+tWv8pOf/ASv9+j+fUnH+tteAkyogDpuCyR1jIEdqfyEMZx6Fi5cyIEDB7Db7Tz77LNceuml7N27d7xKHLF029bDlcrtnW0zcVp1LqdVv3fdXiAcY0e7h22tbjY29vDUW5vRzf0YX318J7CTBWWZ5NhMdOYtZEtXhLn+MJlWE5A+23twcJDLLruMX/7ylzidhzcUTtf3t16vZ9OmTQwMDPCJT3yCbdu2HXbeMl3f38erO92299NPP01+fj6LFi3i1Vdffb/FhtXe7lDpOEjiuC2QSktLaW5+L4hbWlooLi4en+qOYTj1OJ3Ood32Cy64gEgkQk9Pz7jWORLptq2HK922d4ZJz6JpWVyzpISd993JF6a72fG9j/H450/lOxfPYUa+naY+P935S/nJGj8LvvdPzvrpK+wvPJPVPUY2NQ8MTZORCpFIhMsuu4xrr72WT37yk0c9nm7b+0iZmZmcddZZPP/884fdn+7v7/erO92291tvvcWTTz5JRUUFV111FS+//DLXXXfdkYudcHu7dAyooRZImqaZSLRAevLQBS655BLuv/9+lFK88847uFyulO6SL1myhL1797J//37C4TAPPfQQl1xyyWHLdHR0DH1aW7NmDfF4nJycnFSUe0LSbVsPVzpub6UUN998M7Nnz+b2228fCq1Pn1bJTz81n3/efiY/O1WjvP5vfOP8meSbopBVxs9eaebS/3mLed95gcvvepv/emYHz21tp8szov6bJ133saTj9u7u7mZgIDFzdCAQ4MUXX2TWrFmHLZOO7+/h1J1u2/uHP/whLS0tNDY28tBDD3HOOefw5z//+cjFngRuSI7mWw64P+j8E6ThIT71Pi2QNE27tbw8cXL6ggsu4Nlnn6W6uhqr1cof//jHVJaMwWDgt7/9LR/96EeJxWLcdNNN1NbWcvfddwNw66238thjj3HXXXdhMBjIyMjgoYceSotDCVdffTWvvvoqPT09lJaW8t3vfpdIJDE54a233pp22/qg49Wdjtv7rbfe4oEHHmDevHksWLAAgB/84Ac0NTUBibovu+QCXvnHs/zs5o9itVp55I9/pGRGLRsODLCxqZ8NTf3c9/YB/vDGfjKtRpwWI4srslhSkc2Siiym59lH/fccTt3puL3b29u58cYbicVixONxrrjiCi666KLD/l+m4/t7OHWn4/Y+lkNrBp4lMZKvHvADnzne8+VCXSEmmFA0xo42D+sP9LO2sY91jf30+sIAZFmNLJqWzZk1udSVZlJb7JQ2TiIdSCcJIaYipRT7e3ysa0wE1trGPjyBCH3+CHazgUXTslhelcOyqmzmFjsxGSbOLLxi0pCAEkIkdHmDvNvQx7v7e3m3oY+9XYMA1JW6yEmOMjx1ei6zCh3opN+gGHsSUEKIY+sZDPFuQy/vNPTxVn0PDT0+IDEkfsX0HM6dlc+SimzKsq0prlRMUhJQQojhaRsI8Pa+Xt6u7+GtfT3k2s1sb/NQlWvjjJo8zqzJY1lVNlZT2o2jEhOTBJQQ4sQppajvHOT1+h5e39PNu/t7CUbimPQ6llRm8bG5RSyvymF6ni0tR4qJCUECSghx8oKRGGsb+3h9Tzev7e5mIBChyxtiWo6Vc2blc+6sApZWZmMyyOhA8cEGBgb461//ym233SYBJYQYfa0DAV7e1cXLOzt5a18v4Wgcu9nA6TNyOWdWPmfPzCfXYU51mSINNTY2ctFFF7Ft2zYJKCHE2PKHo7xd38tLu7p4eVcnnZ4Q2TYT1Xl2PlJbwEdrC2WghRhycGaHYDC4GfinUupfT+T5ElBCiBFRSrG9zcM/dnTyj+0d7OpIdLGeU+TksoUlnDkzj+p8R4qrFKl0sntQMkRHCDEimqYxt8TF3BIXt3+4hsYeH//Y0cEL2zt5YXsn339mJzMLHFxYV8QF84qozrenumQxwcgelBBi1HW4Azy/rYNnt3aw9kAfSsGsQgcfX1DMBfOKmJZjS3WJYhzIOSghRFrrcAd5bls7z2xpJ64UG5oGWFCWyaULirlofjG5dhlgMVn19vYenLtKAkoIkd5a+/08vaWdv29qY2e7B71O40PVuXzylBI+XFsgFwZPQtdccw0PPvjgduA5GSQhhJgQ9nR6+fvGVp7Y1EaRy8LOdg8X1hXxqcVlLJ6WJRcFTy6yByWEmHjiccW6A/08tr6ZZ7a04wvHqMixcvmiUj5xSiklWRmpLlGcPAkoIcTE5gtFeX5bB4+ub+adhj4WlmfizDByzdJyzpmVL3NbTVwSUEKIyaO5z8+j65t5eG0znZ4QBU4zVywu48olZZRmycXAE4wElBBi8onG4ry8q4sH1zTx6p5uAK5cUsZHaws5c0aezGc1MUhACSEmt5Z+P4+tb+HBNU10ekJU5Fi5fkUFly8qxZVhTHV54v1JQAkhpoZwNM5z29q5f/UB1h/oJ8Oo5xMLS7jptAppr5SeJKCEEFPPtlY3969uZG1jH429fs6syeNzp1dx6vQcGaqePiSghBBTV+9giL++28R9qw/QMxhidpGTz36okovnF8vcVaknASWEEMFIjCc3tfGHNxrY2zXImTW5nFadyzXLpmE3S6eKFJk4AaVp2k+Bi4EwsA/4jFJq4HjPk4ASQgyXUorX9nTz+IYWntrcjivDyI2nVvCZUyvIsplSXd5UM6EC6iPAy0qpqKZpPwZQSn3jeM+TgBJCjMSm5gF+90o9/9jRSYZRzxfOns4Vi8vId1pSXdpUMXEC6rACNO0TwOVKqWuPt6wElBDiZOzp9HL3q/t4a18PA/4I1y2fxi1nVpHvkKAaYxM2oJ4CHlZK/fl9Hl8JrAQoLy9fdODAgfEsTwgxCTX2+PjtK/X8bWMrRr3GdcumsVKCaiylV0BpmvYiUHiMh+5USj2RXOZOYDHwSTWMQmQPSggxmhp7fPzm5Xr+trEFk0HH506v4tOnVpAjc1SNtvQKqOO+sKbdCNwKnKuU8g/nORJQQoixsL/Hx29frqel38/2Ng8rz6ji5g9VYpNRf6Nl4gSUpmnnA78AzlRKdQ/3eRJQQoixVN/l5acv7OaF7Z3k2s18+dxqrlpSLtdRnbwJFVD1gBnoTd71jlLq1uM9TwJKCDEeNjT18+PndvHu/j7Ks618/aMzuWhekTSmHbmJE1AjJQElhBgvSile3dPNj5/bxa4OLx+ZU8BnT69iaWV2qkubiCSghBBitMXjime2tPGD53bR7g5yYV0R3zx/FmXZMifVCZCAEkKIseIPR1n1egN3v7aPuILPnV7JbWdVy0CK4ZGAEkKIsdbuDvDj53bx901tlGZl8PWPzOTjC4qlc/oHk4ASQojxsqGpnwfePsDfNrWyrDKb7186l5oCmYvqfUhACSHEeIrHFQ+va+bHz+9iMBjlpg9V8pVzZ8hhv6ONKKBkcL8QQoyQTqdx9dJyXv6Xs7hsYSmrXm/g3J+/xks7OplIH/7TlQSUEEKcpGybiR9fXsfjnz+V6fk2vvDgBlY+sJ4OdzDVpU1oElBCCDFKFk3L4r7PLOVfPjyTN/Z2c94vXuOBdw4Qj8ve1EhIQAkhxCgy6HV87owq/vHVM1lQlsm3/r6NK37/NvVd3lSXNuFIQAkhxBgoz7HywM1L+fmn5mM1Gbjg12/yh9cbiMne1LBJQAkhxBjRNI3LFpXysyvmc8aMPP7r2Z1c+fvV7O/xpbq0CUECSgghxli+w8IfbljEL66Yz55OLx/71ev86a39cm7qOCSghBBiHGiaxicXlvKPr53J8qocvvPUDu54fAttA4FUl5a2JKCEEGIcFbos/PHTS/jFFfN5bms75//ydZ7e0pbqstKSBJQQQoyzg3tTz3z5dKry7Hzxrxv5l0c24w1GUl1aWpGAEkKIFKnItfHorSv48jnV/G1jCzfeu4bNzQOpLittSEAJIUQKGfU6bv/ITB6+ZQVmo57L736be9/cL62SkIASQoi0sKQim7uuXciZNXl87+kd3PLAetz+qX3ITwJKCCHSRKbVxB9uWMx/XDibl3d1ceFv3mDTFD7kJwElhBBpRNM0Pnt6FY/eugKl4PK73uLhtU1T8pCfBJQQQqShU8qzePbLp3PN0ml84/GtfPXhTQTCsVSXNa5kVi0hhEhTLquR71xSS57DzC9e3MOezkFWXb+IsmxrqksbF7IHJYQQaUyn0/jSuTO499NLaO33c9Fv3uT1Pd2pLmtcSEAJIcQEcPbMfJ760ococlm48Y9ruOuV+kl/XkoCSgghJohpOTb+77ZTuXBeES/v7uKrD28iGJm856UkoIQQYgKxmgz8+qoFnD0znyc2tXHNH96h2xtKdVljQgJKCCEmGJ1Ox21nV3PXtQvZ0e7h0v95i10dnlSXNeokoIQQYoL62LwiHrllBZFYnMt+9zav7OpKdUmjSgJKCCEmsLrSTJ744mlU5Nq4+b613L+6MdUljRoJKCGEmOCKXBk8eusKPjq3kPtXH+D/Pb1jUszWKwElhBCTgNVk4H+uXsiHqnO55839fPXhTYSj8VSXdVKkk4QQQkwSOp3Gty+eQ4HTwo+f30WfL8xd1y3EYTGmurQRkT0oIYSYRDRN4/NnTednn5rP6oZerlr1Dl3eYKrLGhEJKCGEmIQuX1TKPTcupqHbx2V3vc3+Hl+qSzphKQ0oTdO+rmma0jQtN5V1CCHEZHT2zHweXLmcYDjGt/6+jR1tE+taqZQFlKZpZcCHgaZU1SCEEJPdgrJMHlq5gn3dg1z9h3fYPIEmQEzlHtR/A3cAE38spBBCpLHp+XYeuWUFzgwD197zLmsb+1Jd0rCkJKA0TbsEaFVKbR7Gsis1TVunadq67u6p0WJeCCFGW1m2lUduWUG+w8wN/7uGt+p7Ul3ScWlj1a5d07QXgcJjPHQn8O/AR5RSbk3TGoHFSqnjbq3FixerdevWjW6hQggxhXR7Q1x3z7vs7/Vx93ULOWdWwXi8rDaSJ43ZHpRS6jyl1Nwjb0ADUAlsToZTKbBB07RjhZkQQohRlOcw89DK5cwscHDLA+t5bmt7qkt6X+N+iE8ptVUpla+UqlBKVQAtwEKlVMd41yKEEFNRls3EXz63jLrSTO54bDPPpmlIyXVQQggxBTktRu6/aSlnzcznSw9uTMs9qZQHVHJPKv3P1gkhxCRjMxv44WV1zC918aUHN/Lijs5Ul3SYlAeUEEKI1LGbDfzppqXMKXZy21828Pqe9BktLQElhBBT3MHDfdX5dj53/zre3pceB7UkoIQQQpBpNfHAzUuZlmPl5j+tS4uLeSWghBBCAJBjN/Pnzy6jyGXhM39cy6YUt0WSgBJCCDEk32Hhr59bTrbNxM9e2M2ujtQ1mJWAEkIIcZhCl4W/fHYZ9d2D3HjvGlr6/SmpQwJKCCHEUcqyrdz3maUEwjFuuHcNfb7wuNcgASWEEOKYZhY6+N9PL6G1P8BNf1qLPxwd19eXgBJCCPG+llRk85urT2FLywC3/WUDkVh83F5bAkoIIcQH+khtIT/4xDxe3d3NNx7fwljNgnEkw7i8ihBCiAntqqXldHtD/Pyfe8hzmPm3j80e89eUgBJCCDEsXzynmu7BEP/Y3klpVgbXL68Y09eTQ3xCCCGGRdM0vn1xLXOKnHz7ie28vGtsm8tKQAkhhBg2vU7jp5+qY06xky/9dSM728fuQl4JKCGEECfEajJwzw1LcFiM3PyntXR5gmPyOhJQQgghTlihy8I9Ny6m3x/hc/evIxCOjfprSEAJIYQYkbklLn599SlsaXVz+yObiMdHd/i5BJQQQogR+/CcAu68YDbPbevgnjcbRnXdMsxcCCHESbn5Q5X0DIb4wbO7KMm0cmFd0aisV/aghBBCnBRN0/jah2tYWJ7J1x/dPGpTdEhACSGEOGlmg567rluE3WJg5f3rGfCffPdzCSghhBCjosBp4e7rFtLuDvClBzcSO8lBExJQQgghRs2iadl895K5vLG3h5++sPuk1iWDJIQQQoyqa5aVs7XVzd2v7WNuiZOL6opHtB7ZgxJCCDHqvnPJHBaWZ/KdJ3eMeB0Tag8qNzc31SUIIYQYBrNBz93XLaLfHxnxOrTxmnhqlEyoYoUQQgCgjeRJcohPCCFEWpKAEkIIkZYkoIQQQqQlCSghhBBpSQJKCCFEWpKAEkIIkZYkoIQQQqQlCSghhBBpSQJKCCFEWpKAEkIIkZYmVC8+RtguQwghxMQje1BCCCHSkgSUEEKItCQBJYQQIi1JQAkhhEhLElBCCCHSkgSUEEKItPT/Ac4JHpdF7IKcAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "p=sym.plot(solpar.rhs.subs({a:1}), solpar.rhs.subs({a:1.1}),solpar.rhs.subs({a:0.9}),(t,0,4), legend = True,show = False);\n",
    "p[0].label=\"a=1\" \n",
    "p[1].label=\"a=1.1\"\n",
    "p[2].label=\"a=0.9\"\n",
    "p.show();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lorsque $a=1$, le problème est numériquement mal posé, i.e. sensible aux perturbations sur la donnée initiale.  \n",
    "Toute méthode numérique induisant une erreur sur la solution, cette erreur va s'amplifier avec les itérations succéssives. La seule issue pour traiter ce genre de problème est d'utiliser des schémas d'ordre élévé (ici donc le schéma EM qui est d'ordre 2 tandis que les schémas EE et EI sont d'ordre 1) avec un pas $h$ petit."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Circuit LC\n",
    "Un circuit LC est un circuit électrique contenant une bobine (L) et un condensateur (Capacité). C'est ainsi qu'on obtient le phénomène de résonance électrique.\n",
    "![Tuned_circuit_animation_3.gif](Tuned_circuit_animation_3.gif)\n",
    "\n",
    "Dans un circuit électrique de type LC en série, le courant $I$, en ampères (A), et la charge électrique du condensateur $q$, en coulombs (C), évoluent avec le temps selon le système d'équations différentielles\n",
    "$$\n",
    "\\begin{cases}\n",
    "L I'(t)+ \\dfrac{1}{C}q(t)=0\\\\\n",
    "q'(t)=I(t)\\\\\n",
    "\\end{cases}\n",
    "$$\n",
    "avec \n",
    "* $L>0$ l'inductance de la bobine, en henrys (H) ;\n",
    "* $C>0$ la capacité électrique du condensateur, en farads (F) ;\n",
    "* $t$ le temps en secondes (s).\n",
    "\n",
    "> Posons\n",
    "* $L=2$H,\n",
    "* $C=0.5$F,\n",
    "* $q(0)=1$C,\n",
    "* $i(0)=2$A.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ">**Q2 [3 points]**  \n",
    "Écrire le système sous forme matricielle ainsi que sous la forme d'une EDO d'ordre 2 en l'inconnue $t\\mapsto q(t)$. Calculer ensuite la solution exacte.  \n",
    "Afficher $t\\mapsto q(t)$, $t\\mapsto i(t)$ et $q\\mapsto i(q)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Correction Q2**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On a l'écriture matricielle\n",
    "$$\n",
    "\\begin{pmatrix}\n",
    "q\\\\I\n",
    "\\end{pmatrix}'(t)\n",
    "=\n",
    "\\begin{pmatrix}\n",
    "0&1\\\\\n",
    "-\\frac{1}{LC}&0\n",
    "\\end{pmatrix}\n",
    "\\begin{pmatrix}\n",
    "q\\\\I\n",
    "\\end{pmatrix}(t)\n",
    "$$\n",
    "soit encore\n",
    "$$\n",
    "\\begin{cases}\n",
    "q'(t)=\\varphi_1(t,q(t),I(t))\n",
    "\\\\\n",
    "I'(t)=\\varphi_2(t,q(t),I(t))\n",
    "\\end{cases}\n",
    "\\qquad\n",
    "\\text{avec }\n",
    "\\begin{cases}\n",
    "\\varphi_1(t,q(t),I(t))=I(t)\n",
    "\\\\\n",
    "\\varphi_2(t,q(t),I(t))=-\\frac{1}{LC}q(t)\n",
    "\\end{cases}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Transformons le système en une EDO d'ordre 2:\n",
    "$$\n",
    "\\begin{cases}\n",
    "L I'(t)+ \\frac{1}{C}q(t)=E(t)\\\\\n",
    "q'(t)=I(t)\\implies q''(t)=I'(t) \\implies Lq''(t)=LI'(t)\n",
    "\\end{cases}\n",
    "$$\n",
    "donc\n",
    "$$\n",
    "L q''(t)\n",
    "=L I'(t)\n",
    "=-\\frac{1}{C}q(t)+E(t)\n",
    "=-\\frac{1}{C}q(t)+E(t)\n",
    "$$\n",
    "ainsi\n",
    "$$\n",
    "\\boxed{Lq''(t)+\\frac{1}{C}q(t)=0}\n",
    "\\text{ et }\n",
    "\\boxed{I(t)=q'(t)}\n",
    "$$\n",
    "avec\n",
    "$$\n",
    "\\begin{cases}\n",
    "q(0)=q_0,\\\\\n",
    "q'(0)=i_0.\n",
    "\\end{cases}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On peut calculer la solution exacte à la main:\n",
    "$$\n",
    "\\begin{cases}\n",
    "q''(t)+\\frac{1}{LC}q(t)=0\n",
    "\\\\\n",
    "q(0)=q_0,\\\\\n",
    "q'(0)=i_0\n",
    "\\end{cases}\n",
    "$$\n",
    "Il s'agit d'une EDO d'ordre 2 linéaire à coefficients constants et homogène. \n",
    "Le polynôme caractéristique est $\\lambda^2+\\frac{1}{LC}=0$ et a les deux racines complexes conjugués $\\lambda_1=-i\\frac{1}{\\sqrt{LC}}$ et $\\lambda_2=i\\frac{1}{\\sqrt{LC}}$ ainsi la solution générale est\n",
    "$$\n",
    "q(t)= C_1 \\cos\\left( \\frac{1}{\\sqrt{LC}} t\\right) + C_2 \\sin\\left(\\frac{1}{\\sqrt{LC}} t\\right), \\qquad C_1,C_2\\in\\mathbb{R}.\n",
    "$$\n",
    "\n",
    "\n",
    "Sinon, on peut utiliser le module de calcul symbolique mais **attention**, pour ne pas avoir des mauvaises intéractions entre le module `matplotlib` et le module `sympy`, on utilisera l'import de sympy avec un alias:\n",
    "```\n",
    "import sympy as sym\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy as sym\n",
    "sym.init_printing()\n",
    "from IPython.display import display, Math\n",
    "\n",
    "sym.var('t,t_0,q_0,i_0')\n",
    "L=sym.Symbol('L',positive=True)\n",
    "C=sym.Symbol('C',positive=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAKEAAAAuCAYAAABAt3m7AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAJdElEQVR4Ae2c7XEUORCG15QD4CADLgM+IjjIALgIgAyg+Mc/CjIAR8BHBkAENmRwZADlDLj3EZKQNK2ZWXt2ZtZWV42lbbWkVner1dLIc/Dr16/NkvDixYur6v+55+GGTx8Jf7okX63v+SRwOF9X1Z5eyeCehFLl3yj/Vc/fAdfSiy2BKysY3mMZ3t2Ej1fK3xDuZoJr2RVKQDoKK1cvd0N0azBCvOBJ7yha4eokIMN6KqbGOgqcCvQmHMwZE4oRmD7Swwx6r99xGQ7cCYcnvK+0LcdBKCtLpZv7YumO0mdjWRMtRniq9G1ZZ1YjDJ2LEXZDD5R+DDhS/cZIP+i5pXzbmCCUlYH0wkbyi9JbFmu+nJj+nvLfUxr9Bv+P0ky3sy/HYiDEf58LBvGOeMFmgKlg1pdHR2wea4B+0eVPg4B61M9gdiNU7/f0fE9ng/Iw/UwpsweXTQwBrsH6JPBQuuksqQmbHf2GMl+P+njTCEsc0TBTohf0xsYMwQhDoEusODreiKNpmZ1KQPohFsyWWKND9JuFWQUN9R/qiYa8UyP0Bob7peMfejA+DO2lngDECcwM0giq29m0xMKWmVwCia6O1fh1PejosZ7PKvumFMDLRQfiMPqjcgwPfaFHVrCbwhHbHyt9rTQF6tPO7o3QMxY2GW72CPfJcxMHItxfKYctP78EvK5YjeJmQjh0xynFQcLRbeU78aBo0CfGiqe8qxQjq8F/KsBoI1yJuQkzYoIZwSBYYlP3Tf6bcKcTdteaOocEEl3x5irVFa0GDxh6QK/WhiOUY3xlnVAWUupn8f5OjFCdsARf1aCiy/UcMAOiF/S4liwrgZquCJtKXV0Trs+BjNEvho4xR9iVERJ4ZgPwM44ZEJbkyETLLCoBS1foaStdbaHfjiFPboSeGSy9NDYGu1F5ZpzgGiwjgR5duZjN0BVLaebFEs47dXz7CYnLUj9b0ic3wqTHMr6I8YKYI+BlpjVYhwRKXT0QW1ZsB11Nb1G/DAkdV2jxhFl/kxuhOj9VJ3i7yKxniBgjXFSIuzDhFgPxFXnsY2IsXV8bayyr6IpjmVpsh2HeqYwlGpfaxdvxwsEyZF73ZfjDSoPnRTOTjsTEU6WcOb3Tw0xhBwaus80XblbwfDAjs1lZYQKB4r3LM68K+V6hU13BeJBHGU5Rhh459bCAlwtvJCOMeNMjKww8OwPuXGBQZYwEa8edAlgtjH1SWbnbpXzvQONgbJPdAtk7AfQwLNlgRBhTej4YawjPOR+XTzJvFgl6MqrDyoMdZTekOp5QBG62K+WmC+94zdsSPX2tukjjYal4XhuXL+ftTRYyCP9az1c9XEEj5LiokMV2xiA50sGTZd7MoLNQeEvqZ2DGhBIy8RvQ9w7wN8X+/UUIfeEAywUz9qcxNOp1hGjQ7TOK8VdPMGQbrIZbXzBRPWRKvc5qahqhiGEEsOKC3yX7+3fyWyD7K4o/nMs4iNfDe3zi377JRhzZN5H/NPwnB73pPTvLsa/zL6kYqc4IT7dXicZDLBgC7xrvTMC+FYD6nHl2ZnStwX3ASzajby2Jlut2T/Q81TO4WYNOMoDelH3NCFmOtw48lxS2Boi7Z/Ye65n1FsiS416qb29QgwYIf6LtpesYoSrgCYC98YKeZ9x93EwIN9stECct4494gKcgT4PCRHHBg+Xu0kDHCDVydkfAXsSDUhi7XevGDmMovTm01oYDWmBoZwgN9fG6gyDezBhosOIlI7CM0M1cCbDqCVG8nnhMofysn3FQf+kZFkuwdWOHkKKM7TjVj3wbumbsZZ2SjLgGY14MkL86/6JnGz46Z3tz660mMMsIe+NBMe5u26rBqMzCKGp97QrPJiGbMOLHHQcIP9qbqw4Kpd5QnSFDVhO7BfGK7M99fruw3qKQMiMUU84LqvRdpOhmiLtWsTMUHxgOT2k4bhwqz4xTdCyl0FvQqUP7elB4CtTvW9Ijreq2mDBKo57JjFBkIR4sledakFDZapcKr7c+X0m59Sew/2Z0D10tnsviQY01HOeU7eAJy/6MrtyusMWEpmRyZGmE7n2xFFAKfiMcsRdGuNP/CVE/eJrnns1gMI+ELz0SPHFexYQJdODCLRDrWIBx8V7cgmhcagMeON234kOWwY58rAYbbpwEDr3Aj0SOIp0yhUtvSoAjTgQ+qqxjDL+LJvu77Ve68Hrhxg5MBC9leWzCjHRs0Ac48y2Q0EBLzyaBzi2aszVTryWjxYCDkZvfn0lri56dNnGnCwmUMgm4ucGXGQY9kGhmvwWS8n8Z8pIxKwVyJoRJnZK7aeV1xhuSUW9hyuV4chmKEQwHA8K4Mu8kHIMpb6wQR53oOStksZ3RCGEFfZwlXkOo1L+0IJ0RkhEuvdTDsU80QuV55xw2Y6PldGUOaYoxt/NUX+WGBzye7mfgQ7Rv9cSBCY+xcKVs0Av6Nmiz7McXuc3C5LdAYuMXOCP5c1KAE0EffNSIq22pnjb6HWJodFrVQSmmWYxQneKdsu/PeEZqeFesQbGUs1kaPBMT7WK3QPxYLnrC4fhtdCFZ9zkEYm42jCE2H5TLzpdjz0HNO4EPsydjVoNgNuHSR32lS/Sj4g86ES1CmuwWCG1eZJCs0AMOIVt+K2NmVXtfKTPRk29MxHAwHmbCDz24ZeI+BsDuGsPDpRMPBuPErcfvlvg2+HoDdBv/m3T07KJeg/NLwMuejSErWXYt32pdNOiV460+b5lVPcx+nfOHOsaocMd4L2cwSsNmxMUI+k1a/W6JyjFigtv2lS4JYQXgHIH4QCeDIP3hUEYbIA1OZoTqnBmAAY79/kxtFxtu95JGUPtBGBHXMrNIwL3AUE9m2DQFB5MZoZghbmAHVb5XxjtaAzDxqr/TNzJTCO2StXGN8Uovo0Ih0WEDeMPRMOXu2LrN4uIDcROWZMcYjCrDspvhXWH7szYJsNEYBdIrXhO9bgWTGKE3KgyrNCoMc6Py8swIL5jhfRugG6xLAiGWH2Ncd6THreJBhjqJESYyK112jPvEXPr9mYinLmVKxgwy6aplZ5JAePMRYkOzW+kQupdm4QByEiMUA8QAzJhoSN6wOFs68TzE///Qb+KMsHt2S7Pot55Bvt2W7FAC0gt64pIILwN4X5yBcHf1OANUih1sDZOdE4oBjOlIz7Ee/tvtnR6YgkFw8dvHosVY2fKzm97od7mZAd1gRRLw+kWX6C4YG/HiB5WV4dZWnP8PjD+bix9V7vsAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$\\displaystyle L \\frac{d^{2}}{d t^{2}} q{\\left(t \\right)} = - \\frac{q{\\left(t \\right)}}{C}$"
      ],
      "text/plain": [
       "    2               \n",
       "   d          -q(t) \n",
       "L⋅───(q(t)) = ──────\n",
       "    2           C   \n",
       "  dt                "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAAyCAYAAACkuVf4AAAACXBIWXMAAA7EAAAOxAGVKw4bAAARuElEQVR4Ae2d7ZXftBLGNzkpIIQK4HbASwWEDiCpAOgATr7xLQc6CFTASweQCpLQAdwK7pIOcp+f1uP1m2TZlm3t7sw5XtuyNBqNHs2MJP+99969e3fh5BpAA99///0fOn2n81+ukbEGpJeHSv1Nx5e6fjvO4Sm1a8Axnu6hJRi/n2blT++KBgQajOJvOrvjiHS6dIPD+E7Hm2aQRXJ6co0acIzP98oSjLvzmNfnrc/RDKp/dP7p1jd2YwOlI5xrcCAbWXnxAzXgGM9Xdi7G7/myVb5Sb2NOAeVrtesbnT+urX2SiWWiNzo+1/U/NckneX6QPB/q/GVNcrksYw2ojxzjY7XMpsxh3Gcesyq8vRkEjg/VOoxgrQbwsWRDxksdVZF0x+zjo8YwVSWbC3OtAcf4tS6WXs1h3J3HUo3ervxskD8XSKqK6jsq/lzXLKfVujmNA/lB8jFDcqpTA47xbf0Sxbg7j22KvbGlZfCYyj/S+ceKG8HM4/da5ZPukA3Hy+zNqTINOMa3d0gK4w+2s3cON1QDGLzqNsgFVhzGNzqI5lmyYmmIN8Fe6VyjoyMy+0OyMQOpdQYnEe8kOcbLdPskxn3DvIxybxQXGblvJTAD6z1dV7kkJLm+kHy8OnyvduVKRjb1WV6rde+odhUWl0994RgvqNUpjPuyVUEF3yBWRPa/CxBVOo5Gj+x33JTfnDyXrF9In7730XReBSfHeNlOGGHcnUdZBVfPTQbuIwnJchBgqJlYvvqzZgE7spmc7CM5nawBx/guHTDCuDuPXfRcNdNnku6tBli1UX0TwePgeFOmepK8zOAYXES7TudrwDFeuA+mMO7Oo7CSbwC7mxDRI+OFAGvRDte1Lwmxqc+PBnF6TudqwDG+j/57GHfnsY+Sq+Qqw8aSFUa49oi+t98hudk8r90o/9p0enB8zbWfDtaAY3xXhfcwvtl55EZaufl2bbozN8PWRvSVquSR5AqvvQo3ODsi+mqX2dCh5GPpigPH53SeBhzjO+l+iPFNzkPMeB2OaDaHMADkdzpPAxg29jtq/z0C75U/lJxsQD/Rucbfd0z14msl5o6HqfKetl0DjvHtOkxxaDH+IJUr9UwDmqWE93MHtvL9qSN8C0jn6n6cRlslF1EuBisAkLSG+BHYTzpYOuEjghi3otTw5vcC/D+NvfTzifjT+VWT2o9zu4kRPLMji3yr0rF0ehq2D1aEY3xfhbcYXzXzaID4TOdJI8pzHX/rwNi2pHsiSIwvQK6KJBOzov82QvHPftpDaZe6f6EzewV/N3lKn9AVetnl67aNzuFf+6yjtF6P5BewIV1X5UAkz9nYPqQPHOOHqLnF+NqZB79OxpjGiMGDMbycyEA5ylfxWmMDuPAWgWT6TPejtXWl8YM6ImHatMt+gfgzM9vzF9/myPdyflLNnSeb1bF0tQknwkIYQzqvnoWqLMHC6dg+EBWO8f2V3WL8/sq6WIdOgRpDO/k11KYc5WuZfbyUrEx1P5ZMI8fR0Q+DcNf9AtX/tlNf6UsbWD7zKK3Za36m2/9cJ62+YnxsHSPVYHu1FpYVdIwv09ea3C3GFzsPGTj2OoxBrHKipt9jD5vyTxLPD3mktjADIkr8KsNwM4uyV9UOka9wJY8afnN9V7jau8OugyHT9WmNv2PYNj2b3h3jppHC5y7Ge8tWeoDnxqC+0vG+Dj5hwQYySyoWlTOrGE3J9RyHwVIU0RJ8Ul9DpTx8UrMXPd6PmrayFswMKeXoTAgAmVqqs3yTZ9WBk0JHzC4sonyqdGY86MuWF9B1+MBeU+ZnPeM5jou+wXlDn+pY8qVZq5P6nfbVgOl631oi3Bs8HYlt2gs2u9gCm+240rXlsWVTZmd8+LK1JbqOjpFIU4fJpveuHMM8fl9GAw9b56GOw7BhHNt/+ak0DBoffOt+2ZQlnpERVR5AgOHDuD3WGecQIwBEfbMkPtSVlbfD7C+Vm/vCqe25jNrS4dNeih+ANAfapudcqGwYODr3dKL7IIPOOCacCG9btaR76iMdfRFV4ZDDa6s6M9De6Mx+TE6kRTAAMYNy2k8D4IS+OpOOxDaBDbhl9h6chc6MVxwDgRljEaxiS1rbousLpfMWY8C0zskxQv4McoxnKKlAloDxBzBqOo7O5TXRoSEaGkw6OWWAMJDDMkrqEeUB3SxJHhsIs3kXZrAIvo2OlpRvdMbSW87//8bhfkIZHSjeaOi4YnqlT3DIrUPUNYMSPgzUnBkc/ea0vwbow7N1vRrbwhSyP2vUZGM0tayL3RjO3q39bxs+5JkKcnhbkwCIwBOHmzNGlC1KVm80gz8oooGA8eA8xI4pJ4ZtaISIGIbGlU42UOhyRBizYZlhJozh2R1NOy7U5qGzHMoa7pWvNfy6Ri+0E5pth/IzI0Ph/3KtM6/8MpjCLELXOfQ6kmm2/qactTfVd6MqJOO7UeIdTJAeurPvOQ0EXc9l4rn4xmbW1l9TwdPczNrKLsa2ROKfWrV1NvIxsxi9BKBnYO8jHT0cK53xH2yArnmOE3qlo0d6ZgEQy7cErlvHiLU7G+Oq0/GtXpEeluCbfnxkzoMIul175ImY0eEc2d9BUhnAlFOGTs7uYOXdg7IjRLWLSI6B+BeC6J4zwLcIj+Q54vcbRHSUwVkzSPnhYTtQlVYdSb6loKquDTULFOt/pYMTvsrQM8yZbVmNbfH/WnV29yLAKmksLwX8d2RgrEP/uzpN/rU8qfGOg4EOHyNqk+P7SveL/z6Q8jD4HEMnESJrPe85FeVLAXNUBv46hsChPvjMksrGIrNU2bnIjLK0i0HBAJ2L0D5VnrnZVFQe6tDDS52ZpocfVuqaFxFe6MwxHJRRXhseBH2rrqn+2MDWi05oIAvbE+VKJW3BNsFMbJY7lM/GzWhW0sloeRjzMWLZq8QYcYzHNFw+/dJmHrC2TrZqWF+fMmrks2jC8tq5t98hQPQidsukMzOPYX2dx9eX4rFXZB4iKtWEjNHoTvWT7/m1RKuuiKzQWVuP+DLrQF+f6JjSs5KLkjnwGmZ9RRtWGbNsbO8o92psg8uBXIy/sPE9SL9Q3rc6wC4YHpGe8bINy7NgD6z3AjClh2BT6eyJlBgjjnEp8gAKGGfmAQCIVFqHoHuiYjq2NXYdgQALr4lOUTtwxINIg6i+B5imENPTIwzmlIwhTXIxIHCQTNHRQW/Q6J72A/jnPA+Ftv3hcy44jC4vdDSc2ZE2JPR6OUxceG9LC1P8F7LKz672YhRYM99E4rN6eaGR4akE4PXVVbSw/qzAaJUgGYUkaxFsN3ojuGK8xogxxKY3s/h2DOkaB8aqAfSZjpdKG86yyfOj0u0tzdwxAs8pOhzjkh18b8KWNUS8tmJ89ThbUfc/NvMAAD+LgQ0uA/9wKYt2/qKDSGGKWJIBIDifC52nnA+PMMx7zSjgn0WSj6joPWVm/4G2mmHHUONUwhJTFrN0JvgysBhglvN9XVAvAx0Ash8SIjjdo9+vdODQSec5Do50ZkGkm/4YcCyr0Ycpsral8uzxDAMx9+v9Pert8kQGNmRL9WeX9/Aa57zV0Q95Lr5XWzdhW+XBmPVdFDvKB34/IK/OOAv7HQd2INgRnVlGxgHRB8YL/txb8ER6dIzoWQ4Z75y8pfIcia2UzNZXRwTlDyXI5b1378YvG6hDbT1+0hPqOQDhw4GLBVUZQMP73al10pSSqnmmNhCVMWiqb0sjK86Hd+1twO6qS9WD00M/zOBOoSNlUF0Mqn91EE1vclQqD7bWbphv0rXqNsMeApTm/kJnCyo38d+rsORDZ4dhXPWdjm90eaQcqqvFuM08hv3JYE85BrwcwLLod1g+dc+gorzTsRqwgY9hOIqYNZ3d10fKYLq15ZMteiaK5jiUZBxoAzMIZgUYR4hxvskZBi77/zka40diK6W9I+VoMX4/IhHLStHoVKBieklUZIwibPrJTX7KtWuj/Rw344526MAo0mnhWve25FdlIySfBQOp9etisqOjRjdRHBWrLMLoBBlsPJgRi0g2nyzZ2Qc4Y5ywbs7452wHy62HO7J5LfVzSMbDMK66Tsc3rT9BjhbjvZmHBMEgAhymJrwpgXCxiCNsNivfkiUJIpo1sxUVq4ekE4xDTC/1CDqWBLnDvsr4UfEU9JPUkfQIzlgiBUNd41Tqn28lZVD9OPxPdbDcAWF80FGon4SFZAPLjNjC4udnl07YA7zJdBTGk9gyBZ6J8R3wTbNajA+dR3Kwm0I4SzA2cPk0x7c6YhvjbRHy6Yb8dK7TORrAqBEc7ErqY5wCn5qIBgoNHpi5PdfB/lnrPHRN4EKggayrlr1UPkeGgFvlZeOPjd+tszIc0YX4OMZRxDm0O8bVv7PYounKh807DeOqvzS+aVaL8Z7z4MkSknAMklnHAU9ryBL+nre4BnijDMP8UEdrrIvXcmXwcQojom4lsqlJBLPnP9/C6UzKoPSWJI+t60+9Ut7my7zA2Z22TJcp423PdgTGk9iqCeOF8Q12Wozfv+1I8vb1NGC/HAYAu1AzcB7rHDPGL1UxS2dzr+/iYJjdLo7iM2Tott10MfVaejdf8lp14gxxjJv4JCvxhzka2BXjmdiqCeNF8I3ihxh355EDx1uSR53PlJ4Zx9Mdm8Q0fXKpSfWTTqSf+kqriXapi1/tZuE5KsMEn6ALybZ1xmCDNOY0J6r2pNIaOADjSWxViPFS+Karehh/ULrznF/1GsAgP1kjpQYGG8vMBlKGlmWx0e9elEZkzhowews5BpYZxwsdPVLZ1TL0GF3f4MxwqluJTX/atnimtLViLz/SwCqMb8WWyteI8VL4Rsk9jPvMY4S7W5/AchB7HhZFZDVY+TH8P+uYnFXApMkzMvg8E9nmeez5Va7mr3jhpHpGveG/RYZhHaaDlDPslUnc4NSy2pbg4Y/KaGAxxgthqyqMq00l8U3P9DDuzqMMWG8MFwEKQ0l0zKvWWaQyRC/85oCDT3NzP0W8TRd7gQLgQTmzjqucnb+FZOhwDJf2mvmmfQrJFj7HI45n/C5j2KY7f6/+WITxgtiqDeNF8A2gpjDuzuNuDjVmD18LEGzwzpLy8W0i9krs7aVnw0INuFKR9yPKKF/Wso7y9WTTfQkZhmKHyEy8ozOPoRxDBs09r7jb12MjWTz5YA1kY7wgtmrDeCl803UjjLvzOBjRNVSnwUKEjDMYOYGUfCpHGcqyr/HhIC/gSkXebIBnkXgTwQ35h7IbZRjWzwyqtzTWzaC6mFEEg9BN714rDwMUWWm/UyUaUL8sxngBbNWG8c34pjtjGHfnUQnYTxADY8cPPHsRfoYctufRGkvxwNjPRd4hulfeSacwqJevBEeNuvKulaGtRvxtPfiXNnF8wUck52ZK6IG2z+Ubc/eUvTWwBuNbsFUNxgvimz6axLg7j73hWyl/gYvIDAO9dPaBkWSQdJe94GFLWrqcJBuUti48mUlykS/JqzHUa2To1mnrwWHAdx9wrTp4QSC5F6I8OCCOr3Q4VaYB9c9ijKvMWnzT+powvhnfNCiFcXceaOjuEkYPJ5AzG+hqyQYJZXEGr3VmSStKes6gZJOeT7TbBnObX2mPdQTHoXOSV1NosQxtZVcXwYmprtEMp5ED/nO/M2GPp/dplUEdfnu+BtZgfBW2hJuaMF4C3/ReFOOT/8/j/P52CY7SgABPhM3yjEUqWVUrP19cxemEAdMMnNmyyscyGYOTsuYkWCvmn29NzgL0bJKWytDUzau+1M16MNR9+6uXrvzRN9L0LLQhleeKvf89WwPqo8UYV5lV+KatKnsKxpt6i+C7aUcS4+48zkZ2BfULdLwXzw/c2n2MObGUl8iGcqz3R43sHJ8tz8+SQfWyVMUMauuHFLc038su0ID6ahHGz8LWsElnyaF6ZzF+fyis3989DQgoGH+WjcJUN0cDykvEzpJPtsPJ4bskzxkyqE5mLERk/F9upxuiAfXbIoyfga0pVZ4hRy7GfeYx1WN3NE2gYYOY/yA32ge4oyrpNVt6YTmCj97xNWBbcuvl8Zu6NeAYT/fPEoz/H0uFZn7aHAgyAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$\\displaystyle q{\\left(t \\right)} = C_{1} \\sin{\\left(\\frac{t}{\\sqrt{C} \\sqrt{L}} \\right)} + C_{2} \\cos{\\left(\\frac{t}{\\sqrt{C} \\sqrt{L}} \\right)}$"
      ],
      "text/plain": [
       "             ⎛  t  ⎞         ⎛  t  ⎞\n",
       "q(t) = C₁⋅sin⎜─────⎟ + C₂⋅cos⎜─────⎟\n",
       "             ⎝√C⋅√L⎠         ⎝√C⋅√L⎠"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "q=sym.Function('q')\n",
    "edo=sym.Eq(L*sym.diff(q(t),t,t) , -1/C*q(t))\n",
    "display(edo)\n",
    "\n",
    "solgen=sym.dsolve(edo,q(t)).simplify()\n",
    "display(solgen) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "consts = sym.solve( [solgen.rhs.subs(t,t0)-q_0 , solgen.rhs.diff(t).subs(t,t_0)- i_0],  dict=True)[0]\n",
    "display(consts)\n",
    "\n",
    "solpar=solgen.subs(consts)\n",
    "display(solpar)\n",
    "\n",
    "I= sym.Function('I')\n",
    "I= solpar.rhs.diff(t).simplify()\n",
    "display(Math(r\"I(t)=\"+sym.latex( I )))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pb={L:2,C:sym.S(1)/2}\n",
    "solparLC=solpar.subs(pb)\n",
    "display(solparLC)\n",
    "init={t_0:0,  q_0:1 , i_0:2}\n",
    "solparLCinit=solparLC.subs(init)\n",
    "display(solparLCinit)\n",
    "display(Math(r\"I(t)=\"+sym.latex( I.subs(init).subs(pb) )))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On transforme les deux fonctions `sympy` en fonction vectorielle pour l'évaluation et l'affichage:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib.pylab import *\n",
    "\n",
    "\n",
    "Qfunc = sym.lambdify(t,solpar.rhs.subs(init).subs(pb),'numpy')\n",
    "Ifunc = sym.lambdify(t,I.subs(init).subs(pb),'numpy')\n",
    "\n",
    "tt = linspace(0,7,101)\n",
    "yyQ = Qfunc(tt)\n",
    "yyI = Ifunc(tt)\n",
    "\n",
    "figure(1,figsize=(17,12))\n",
    "\n",
    "subplot(2,3,1)\n",
    "plot(tt,yyQ);\n",
    "xlabel('$t$ (s)')\n",
    "ylabel('$q$ (C)')\n",
    "grid(True)\n",
    "\n",
    "subplot(2,3,2)\n",
    "plot(tt,yyI);\n",
    "xlabel('$t$ (s)')\n",
    "ylabel('$I$ (A)')\n",
    "grid(True)\n",
    "\n",
    "subplot(2,3,3)\n",
    "plot(yyQ,yyI);\n",
    "xlabel('$q$ (C)')\n",
    "ylabel('$I$ (A)')\n",
    "grid(True)\n",
    "\n",
    "\n",
    "\n",
    "L=2\n",
    "C=1/2\n",
    "\n",
    "Ee=[0.5*(yyQ[i])**2/C for i in range(len(yyQ))]\n",
    "subplot(2,3,4)\n",
    "plot(tt,Ee);\n",
    "xlabel('$t$ (s)')\n",
    "ylabel('$E_C$ (?)')\n",
    "title(\"Electrique\")\n",
    "grid(True)\n",
    "\n",
    "Em=[0.5*L*(yyI[i])**2 for i in range(len(yyQ))]\n",
    "subplot(2,3,5)\n",
    "plot(tt,Em);\n",
    "xlabel('$t$ (s)')\n",
    "ylabel('$E_I$ (?)')\n",
    "title(\"Magnétique\")\n",
    "grid(True)\n",
    "\n",
    "\n",
    "Et=[Em[i]+Ee[i] for i in range(len(yyQ))]\n",
    "subplot(2,3,6)\n",
    "plot(tt,Et);\n",
    "xlabel('$q$ (C)')\n",
    "ylabel('$E_T$ (A)')\n",
    "title(\"Somme\")\n",
    "grid(True)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ">**Q3 [4 points]**  \n",
    "Afficher $t\\mapsto q(t)$, $t\\mapsto i(t)$ et $q\\mapsto i(q)$ en comparant solution exacte et solution approchée obtenue avec la méthode d'Euler Explicite pour $t\\in[0;7]$ avec $N=101$ points. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib.pylab import *\n",
    "\n",
    "t0, tfinal = 0, 7\n",
    "N = 51\n",
    "q_0,i_0=1,2\n",
    "\n",
    "phi1 = lambda t,q,qp : qp\n",
    "phi2 = lambda t,q,qp : -q/(L*C)\n",
    "\n",
    "tt = linspace(t0,tfinal,N)\n",
    "yy1 = [Qfunc(t) for t in tt]\n",
    "yy2 = [Ifunc(t) for t in tt]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Correction Q3**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{cases}\n",
    "q_{n+1}=q_{n}+h\\varphi_1(t_n,q_n,I_n)\n",
    "\\\\\n",
    "I_{n+1}=I_n+h\\varphi_2(t_n,q_n,I_n)\n",
    "\\end{cases}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def EE(phi1,phi2,tt):\n",
    "    h = tt[1]-tt[0]\n",
    "    uu1 = [q_0]\n",
    "    uu2 = [i_0]\n",
    "    for i in range(len(tt)-1):\n",
    "        uu1.append( uu1[i] + h*phi1(tt[i],uu1[i],uu2[i]) )\n",
    "        uu2.append( uu2[i] + h*phi2(tt[i],uu1[i],uu2[i]) )\n",
    "    return [uu1,uu2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "uu1E,uu2E = EE(phi1,phi2,tt)\n",
    "\n",
    "figure(1,figsize=(15,5))\n",
    "\n",
    "subplot(1,3,1)\n",
    "plot(tt,yy1,'r-',label=('Exacte'))\n",
    "plot(tt,uu1E,'b-o',label=('EE'))\n",
    "xlabel('$t$ (s)')\n",
    "ylabel('$q$ (C)')\n",
    "legend() \n",
    "grid(True)\n",
    "\n",
    "subplot(1,3,2)\n",
    "plot(tt,yy2,'r-',label=('Exacte'))\n",
    "plot(tt,uu2E,'b-o',label=('EE'))\n",
    "xlabel('$t$ (s)')\n",
    "ylabel('$I$ (A)')\n",
    "legend() \n",
    "grid(True)\n",
    "\n",
    "subplot(1,3,3)\n",
    "plot(yy1,yy2,'r-',label=('Exacte'))\n",
    "plot(uu1E,uu2E,'b-o',label=('EE'))\n",
    "xlabel('$q$ (C)')\n",
    "ylabel('$I$ (A)')\n",
    "legend() \n",
    "grid(True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ">**Q4 [4 points]**  \n",
    "Afficher $t\\mapsto q(t)$, $t\\mapsto i(t)$ et $q\\mapsto i(q)$ en comparant solution exacte et solution approchée obtenue avec la méthode de Heun pour $t\\in[0;7]$ avec $N=101$ points."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Correction Q4**  \n",
    "$$\n",
    "\\begin{cases}\n",
    "\\tilde q_{n+1} = q_{n}+h\\varphi_1(t_n,q_n,I_n),\n",
    "\\\\\n",
    "\\tilde I_{n+1} = I_{n}+h\\varphi_2(t_n,q_n,I_n),\n",
    "\\\\\n",
    "q_{n+1}=q_{n}+\\frac{h}{2}\\Big(\\varphi_1\\big(t_n,q_n,I_n\\big)+\\varphi_1\\big(t_{n+1},\\tilde q_{n+1},\\tilde I_{n+1}\\big)\\Big)\n",
    "\\\\\n",
    "I_{n+1}=I_n+\\frac{h}{2}\\Big(\\varphi_2\\big(t_n,q_n,I_n\\big)+\\varphi_2\\big(t_{n+1},\\tilde q_{n+1},\\tilde I_{n+1}\\big)\\Big)\n",
    "\\end{cases}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Heun(phi1,phi2,tt):\n",
    "    h = tt[1]-tt[0]\n",
    "    uu1 = [q_0]\n",
    "    uu2 = [i_0]\n",
    "    for i in range(len(tt)-1):\n",
    "        uu1pred = uu1[i] + h*phi1(tt[i],uu1[i],uu2[i])\n",
    "        uu2pred = uu2[i] + h*phi2(tt[i],uu1[i],uu2[i])\n",
    "        uu1.append( uu1[i] + h/2*( phi1(tt[i],uu1[i],uu2[i])+phi1(tt[i],uu1pred,uu2pred) ) )\n",
    "        uu2.append( uu2[i] + h/2*( phi2(tt[i],uu1[i],uu2[i])+phi2(tt[i],uu1pred,uu2pred) ) )\n",
    "    return [uu1,uu2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "uu1H,uu2H = Heun(phi1,phi2,tt)\n",
    "\n",
    "figure(1,figsize=(15,5))\n",
    "\n",
    "subplot(1,3,1)\n",
    "plot(tt,yy1,'r-',label=('Exacte'))\n",
    "plot(tt,uu1H,'b-o',label=('H'))\n",
    "xlabel('$t$ (s)')\n",
    "ylabel('$q$ (C)')\n",
    "legend() \n",
    "grid(True)\n",
    "\n",
    "subplot(1,3,2)\n",
    "plot(tt,yy2,'r-',label=('Exacte'))\n",
    "plot(tt,uu2H,'b-o',label=('H'))\n",
    "xlabel('$t$ (s)')\n",
    "ylabel('$I$ (A)')\n",
    "legend() \n",
    "grid(True)\n",
    "\n",
    "subplot(1,3,3)\n",
    "plot(yy1,yy2,'r-',label=('Exacte'))\n",
    "plot(uu1H,uu2H,'b-o',label=('H'))\n",
    "xlabel('$q$ (C)')\n",
    "ylabel('$I$ (A)')\n",
    "legend() \n",
    "grid(True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ">**Q5 [4 points]** \n",
    "Afficher $t\\mapsto q(t)$, $t\\mapsto i(t)$ et $q\\mapsto i(q)$ en comparant solution exacte et solution approchée obtenue avec la méthode de Euler-Cromer pour $t\\in[0;7]$ pour $N=101$ points.  \n",
    "\n",
    ">Cette méthode, appellée Euler-Cromer ou Euler symplectique ou Newton-Störmer-Verlet, s'écrit comme suit:\n",
    "$$\n",
    "\\begin{cases}\n",
    "I_{n+1}=I_n+h\\varphi_2\\big(t_{n},q_{n},I_{n}\\big)\n",
    "\\\\\n",
    "q_{n+1}=q_{n}+h\\varphi_1(t_{n+1},q_{n+1},I_{n+1}), \\qquad\\leftarrow\\text{ce calcul est explicite car $\\varphi_1$ ne dépend pas de $q$}\n",
    "\\end{cases}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Correction Q5**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Euler-Cromer, see fdm-book-4print.pdf page 47\n",
    "# The scheme goes under several names: forward-backward scheme, semi-implicit Euler method, \n",
    "# semi-explicit Euler, symplectic Euler, Newton-Störmer-Verlet, and Euler-Cromer.\n",
    "def EC(phi1,phi2,tt):\n",
    "    h = tt[1]-tt[0]\n",
    "    uu1 = [q_0]\n",
    "    uu2 = [i_0]\n",
    "    for i in range(len(tt)-1):\n",
    "        uu2.append( uu2[i] + h*phi2(tt[i],uu1[i],uu2[i]) )\n",
    "        uu1.append( uu1[i] + h*phi1(tt[i+1],uu1[i],uu2[i+1]) ) # en theorie uu1[i+1] mais n'intervient pas dans phi1\n",
    "    return [uu1,uu2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "uu1EC,uu2EC = Heun(phi1,phi2,tt)\n",
    "\n",
    "figure(1,figsize=(15,5))\n",
    "\n",
    "subplot(1,3,1)\n",
    "plot(tt,yy1,'r-',label=('Exacte'))\n",
    "plot(tt,uu1EC,'b-o',label=('EC'))\n",
    "xlabel('$t$ (s)')\n",
    "ylabel('$q$ (C)')\n",
    "legend() \n",
    "grid(True)\n",
    "\n",
    "subplot(1,3,2)\n",
    "plot(tt,yy2,'r-',label=('Exacte'))\n",
    "plot(tt,uu2EC,'b-o',label=('EC'))\n",
    "xlabel('$t$ (s)')\n",
    "ylabel('$I$ (A)')\n",
    "legend() \n",
    "grid(True)\n",
    "\n",
    "subplot(1,3,3)\n",
    "plot(yy1,yy2,'r-',label=('Exacte'))\n",
    "plot(uu1EC,uu2EC,'b-o',label=('EC'))\n",
    "xlabel('$q$ (C)')\n",
    "ylabel('$I$ (A)')\n",
    "legend() \n",
    "grid(True)"
   ]
  }
 ],
 "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.8.5"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": true,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}