diff --git a/6.1 Equations Differentielles Ordinaires.ipynb b/6.1 Equations Differentielles Ordinaires.ipynb index d496469..27fe099 100644 --- a/6.1 Equations Differentielles Ordinaires.ipynb +++ b/6.1 Equations Differentielles Ordinaires.ipynb @@ -1,416 +1,410 @@ { "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Equations Différentielles Ordinaires\n", "\n", "## Problème de Cauchy\n", "\n", "$\\newcommand{\\rr}{\\mathbb R}\n", "f:\\rr_+\\times\\rr \\rightarrow\\rr$ continue, $y_0\\in\\rr$ donné.\n", "On cherche $y:t\\in I \\subset \\rr_+\\rightarrow y(t)\\in\\rr$\n", "qui satisfait le problème suivant\n", "$$\\left\\{\n", "\\begin{array}{l}\n", "y'(t)=f(t,y(t))\\qquad \\,\\forall t \\in I \\\\\n", "y(t_0)=y_0\n", "\\end{array}\n", "\\right.$$\n", "où $y'(t)=\\displaystyle\\frac{d y(t)}{d t}$.\n", "\n", "#### Exemple\n", "\n", "Écrivez la discretisation par la méthode d'Euler progressive et rétrograde du problème de Cauchy\n", "$$\\left\\{\n", "\\begin{array}{l}\n", "y'(t)=-t \\,y(t)^2\\qquad \\,\\forall t \\in [0,4] \\\\\n", "y(t_0)=2\n", "\\end{array}\n", "\\right.$$\n", "La solution de ce problème est $y(t) = \\frac{2}{1+t^2}$\n", "\n", "![Approximation par le méthodes de Euler Retrograde et Heun](EulerRetroHeun.png)\n", "\n", "Figure: Approximation par le méthodes de Euler Retrograde et Heun\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Euler Progressif\n", + "## Euler Progressif (Série 10, Ex 3)\n", "\n", "Ecrivez une fonction `forwardEuler` \n", "qui approche la solution du problème\n", "$$y'(t) = f(t,y(t)), \\; t \\in (T_0,T_f), \\quad y(0) = y_0,$$\n", "en utilisant la méthode d'Euler progressif. L'entête de la fonction doit être la suivante:\n", "```python\n", "def forwardEuler( fun, interval, y0, N ) :\n", " # FORWARDEULER Solve differential equations using the forward Euler method.\n", " # [T, U] = FORWARDEULER( FUN, INTERVAL, Y0, N ), with INTERVAL = [T0, TF],\n", " # integrates the system of differential equations y'=f(t, y) from time T0 \n", " # to time TF, with initial condition Y0, using the forward Euler method on \n", " # an equispaced grid of N intervals. Function FUN(T, Y) must return\n", " # a column vector corresponding to f(t, y). Each row in the solution \n", " # array U corresponds to a time returned in the column vector T.\n", "```\n", "\n", "Une fois écrite la fonction `forwardEuler`, utilisez les commandes suivantes pour calculer la solution:\n", "```python\n", "f = lambda t,x : ( 2/15*x*(1-x/1000) ) # C=2/15 et B = 1000\n", "tsp = [0,100]\n", "y0 = 100; Nh = 25;\n", "t25,u25 = forwardEuler(f,tsp,y0,Nh)\n", "\n", "plt.plot(t25,u25,'o')\n", "```\n", "Approchez la solution de l'équation différentielle \n", "$$ y'(t) = \\frac{2\\,y}{15}(1-y/1000), \\; t \\in (0,100), \\quad y(0) = 100,$$\n", "avec $N_{h}=25$.\n", "\n", "Que se passe-t-il avec $N_{h}=7$? Discutez les résultats obtenus.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# importing libraries used in this book\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def forwardEuler( fun, interval, y0, N ) :\n", " # FORWARDEULER Solve differential equations using the forward Euler method.\n", " # [T, U] = FORWARDEULER( FUN, INTERVAL, Y0, N ), with INTERVAL = [T0, TF],\n", " # integrates the system of differential equations y'=f(t, y) from time T0 \n", " # to time TF, with initial condition Y0, using the forward Euler method on \n", " # an equispaced grid of N intervals. Function FUN(T, Y) must return\n", " # a column vector corresponding to f(t, y). Each row in the solution \n", " # array U corresponds to a time returned in the column vector T.\n", " \n", " # time step\n", - " h = ( interval[1] - interval[0] ) / N\n", + " h = 0 # Complete here\n", " \n", " # time snapshots \n", " t = np.linspace( interval[0], interval[1], N+1 )\n", "\n", " # initialize the solution vector\n", " u = np.zeros(N+1)\n", " u[0] = y0\n", "\n", " # time loop (n=0,...,n, but array indeces in Matlab start at 1)\n", " for n in range(N) :\n", - " u[n+1] = u[n] + h * fun( t[n], u[n] )\n", + " u[n+1] = 0 # Complete here\n", " \n", " return t, u\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Nous pouvons alors résoudre l'exercice avec les commandes:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "f = lambda t,x : ( 2/15*x*(1-x/1000) ) # C=2/15 et B = 1000\n", "tsp = [0,100]\n", "y0 = 100; Nh = 25;\n", "t25,u25 = forwardEuler(f,tsp,y0,Nh)\n", "\n", "plt.plot(t25,u25,'o')\n", "\n", "# labels, title, legend\n", "plt.xlabel('$t_n$'); plt.ylabel('$u_n$'); #plt.title('data')\n", "plt.legend(['$u_n$'])\n", "plt.title('$u_n\\\\approx u(t_n)$')\n", "plt.grid(True)\n", "plt.show()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Pour résoudre le problème avec $N_{h} = 7$ et tracer le graphe de la solution, on utilise les commandes suivantes: " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "Nh = 7\n", "t7,u7 = forwardEuler(f,tsp,y0,Nh);\n", "\n", "plt.plot(t25,u25,'o-')\n", "plt.plot(t7,u7,'o-')\n", "\n", "# labels, title, legend\n", "plt.xlabel('$t_n$'); plt.ylabel('$u_n$'); #plt.title('data')\n", "plt.legend(['$N_h=25$','$N_h=7$'])\n", "plt.title('$u_n\\\\approx u(t_n)$')\n", "plt.grid(True)\n", "plt.show()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Solutions de l'équation différentielle* $y'(t) = \\frac{2\\,y}{15}(1-y/1000), \\; t \\in (0,100), \\quad y(0) = 100$ *pour diverses valeurs de* $N_h$.\n", "\n", "On voit que la solution avec $N_{h}=7$ est oscillante et ne tend pas vers 1000 lorsque $t \\to \\infty$: en effet,\n", "la condition de stabilité sur le pas de discretisation $N_{h}$ n'est pas satisfaite.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Euler Retrograde\n", + "## Euler Retrograde (Serie 10, Ex 4)\n", "\n", "Ecrivez une fonction Matlab `backwardEuler` \n", "qui approche la solution du problème\n", "$$y'(t) = f(t,y(t)), \\; t \\in (T_0,T_f), \\quad y(0) = y_0,$$\n", "en utilisant la méthode d'Euler progressif. L'entête de la fonction doit être la suivante:\n", "```python\n", "def backwardEuler( fun, interval, y0, N ) :\n", " # BACKWARDEULER Solve differential equations using the backward Euler method.\n", " # [T, U] = BACKWARDEULER( FUN, INTERVAL, Y0, N ), with INTERVAL = [T0, TF],\n", " # integrates the system of differential equations y'=f(t, y) from time T0 \n", " # to time TF, with initial condition Y0, using the backward Euler method on \n", " # an equispaced grid of N intervals. Function FUN(T, Y) must return\n", " # a column vector corresponding to f(t, y). Each row in the solution \n", " # array U corresponds to a time returned in the column vector T.\n", " from scipy.optimize import fsolve\n", "```\n", "\n", "Vous disposez de la fonction `scipy.optimize.fsolve` pour résoudre une équation non-linéaire\n", "(`fsolve` cache plusieurs méthodes de resolution, comme Newton ou le point fixe).\n", "\n", "Exemple d'utilisation de `fsolve`:\n", "on peut chercher le zéro de `F` près de `x0` en utilisant le code suivant:\n", "```python\n", "from scipy.optimize import fsolve\n", "a = 5; h = 0.1; \n", "F = lambda x : a + np.sin(x) - h*x; \n", "\n", "x0 = a+h\n", "zero = fsolve(F, x0)\n", "\n", "print(f'F({zero[0]:5.2f}) = {F(zero)[0]:4.1e}')\n", "```\n", "\n", "Approchez la solution de l'équation différentielle \n", "$$ y'(t) = \\frac{2\\,y}{15}(1-y/1000), \\; t \\in (0,100), \\quad y(0) = 100,$$\n", "avec $N_{h}=25$. Que se passe-t-il avec $N_{h}=7$? Discuter les résultats obtenus.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# importing libraries used in this book\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from scipy.optimize import fsolve\n", "a = 5; h = 0.1; \n", "F = lambda x : a + np.sin(x) - h*x; \n", "\n", "x0 = a+h\n", "zero = fsolve(F, x0)\n", "\n", "print(f'F({zero[0]:5.2f}) = {F(zero)[0]:4.1e}')\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def backwardEuler( fun, interval, y0, N ) :\n", " # BACKWARDEULER Solve differential equations using the backward Euler method.\n", " # [T, U] = BACKWARDEULER( FUN, INTERVAL, Y0, N ), with INTERVAL = [T0, TF],\n", " # integrates the system of differential equations y'=f(t, y) from time T0 \n", " # to time TF, with initial condition Y0, using the backward Euler method on \n", " # an equispaced grid of N intervals. Function FUN(T, Y) must return\n", " # a column vector corresponding to f(t, y). Each row in the solution \n", " # array U corresponds to a time returned in the column vector T.\n", " from scipy.optimize import fsolve\n", " \n", " # time step\n", - " h = ( interval[1] - interval[0] ) / N\n", - " \n", + " h = # Complete here\n", + " \n", " # time snapshots \n", " t = np.linspace( interval[0], interval[1], N+1 )\n", "\n", " # initialize the solution vector\n", " u = np.zeros(N+1)\n", " u[0] = y0\n", "\n", " # time loop (n=0,...,n, but array indeces in Matlab start at 1)\n", " for n in range(N) :\n", " # non-linear function\n", " F = lambda x : u[n] + h * fun(t[n+1],x) - x \n", " # solve the non-linear equation using the built-in matlab function \"fsolve\"\n", " # to compute u[n+1] \n", - " u[n+1] = fsolve( F, u[n]);\n", - " # u[n+1] = fsolve( F, u[n]+ h * fun( t[n], u[n] ));\n", - "\n", - " # NOTE:\n", - " # in the call of fsolve, a more accurate initial guess is obtained \n", - " # by replacing u[n] with the forward euler method:\n", - " # u[n+1] = fsolve( F, u(n) + h * fun( t(n), u(n) ), options ); \n", + " u[n+1] = 0 # Complete here, fsolve ( ?, ?);\n", "\n", " return t, u" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Nous pouvons alors résoudre l'exercice avec les commandes:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "f = lambda t,x : ( 2/15*x*(1-x/1000) ) # C=2/15 et B = 1000\n", "tsp = [0,100]\n", "y0 = 100; Nh = 25;\n", "t25,u25 = backwardEuler(f,tsp,y0,Nh)\n", "\n", "plt.plot(t25,u25,'o')\n", "\n", "# labels, title, legend\n", "plt.xlabel('$t_n$'); plt.ylabel('$u_n$'); #plt.title('data')\n", "plt.legend(['$u_n$'])\n", "plt.title('$u_n\\\\approx u(t_n)$')\n", "plt.grid(True)\n", "plt.show()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Pour résoudre le problème avec $N_{h} = 7$ et tracer le graphe de la solution, on utilise les commandes suivantes: " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "Nh = 7\n", "t7,u7 = backwardEuler(f,tsp,y0,Nh);\n", "\n", "plt.plot(t25,u25,'o-')\n", "plt.plot(t7,u7,'o-')\n", "\n", "# labels, title, legend\n", "plt.xlabel('$t_n$'); plt.ylabel('$u_n$'); #plt.title('data')\n", "plt.legend(['$N_h=25$','$N_h=7$'])\n", "plt.title('$u_n\\\\approx u(t_n)$')\n", "plt.grid(True)\n", "plt.show()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Solutions de l'équation différentielle $y'(t) = \\frac{2\\,y}{15}(1-y/1000), \\; t \\in (0,100), \\quad y(0) = 100$ pour diverses valeurs de $N_h$.*\n", "\n", "On voit que `fsolve` n'arrive pas à resoudre l'équation non-linéaire pour $N_{h}=7$.\n", "Il faut chercher un meilleur `x0`. \n", "\n", "Dans `backwardEuler`, replacez le `u[n]`\n", "```python\n", " u[n+1] = fsolve( F, u[n]);\n", "```\n", "par la solution correspondante à la méthode d'Euler progressive, i.e.\n", "```python\n", " u[n+1] = fsolve( F, u[n]+ h * fun( t[n], u[n] ));\n", "```\n", "\n", "Maintenant, pas seulement `fsolve` trouve une solution, mais en plus l'approximation de l'EDO ne présente pas d'oscillations.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.7.16" } }, "nbformat": 4, "nbformat_minor": 4 }