diff --git a/6.1 Equations Differentielles Ordinaires.ipynb b/6.1 Equations Differentielles Ordinaires.ipynb new file mode 100644 index 0000000..75478d9 --- /dev/null +++ b/6.1 Equations Differentielles Ordinaires.ipynb @@ -0,0 +1,281 @@ +{ + "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", + "Approximer la solution 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", + "Avec les méthodes de Euler Proressive et Rétrograde, Heun, Crank-Nicolson, et Euler modifié.\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", + "\n", + "from OrdinaryDifferentialEquationsLib import forwardEuler,backwardEuler,Heun,CrankNicolson,modifiedEuler\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "f = lambda t,x : -t*x**2\n", + "y0 = 2; tspan=[0, 4]\n", + "Nh = 20\n", + "\n", + "for method in [forwardEuler,backwardEuler,Heun,CrankNicolson,modifiedEuler] :\n", + "\n", + " t, y = method(f, tspan, y0, Nh)\n", + " plt.plot(t, y,'o-')\n", + " plt.plot(t, y,'o-')\n", + "\n", + "y = lambda t : 2/(1+t**2)\n", + "t = np.linspace(tspan[0],tspan[1],100)\n", + "plt.plot(t, y(t),'-')\n", + "# labels, title, legend\n", + "plt.xlabel('$t$'); plt.ylabel('$y$')\n", + "plt.legend(['forwardEuler','backwardEuler','Heun','CrankNicolson','modifiedEuler','$y(t)$'])\n", + "plt.grid(True)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Stabilité\n", + "\n", + "On considère le problème de Cauchy\n", + "$$\\begin{cases}\n", + " y'(t) = -2y(t), \\quad t>0 \\\\\n", + " y(0) = 1\n", + "\\end{cases}$$\n", + "La solution exacte de ce prolème est $y(t)= e^{-2t}$\n", + "\n", + "Résolvez ce problème par les méthodes d'Euler Progressive et Rétrograde\n", + "sur l'intervalle $[0,10]$ avec\n", + "un pas de temps $h=0.9$ et $1.1$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "l = -2\n", + "f = lambda t,x : l*x\n", + "y0 = 1; tspan=[0, 10]\n", + "h = 1.1; Nh = np.ceil((tspan[1] - tspan[0])/h).astype(int)\n", + "t_EP, y_EP = forwardEuler(f, tspan, y0, Nh)\n", + "t_ER, y_ER = backwardEuler(f, tspan, y0, Nh)\n", + "\n", + "plt.plot(t_EP, y_EP,'o-')\n", + "plt.plot(t_ER, y_ER,'o-')\n", + "\n", + "y = lambda t : np.exp(l*t)\n", + "t = np.linspace(tspan[0],tspan[1],100)\n", + "plt.plot(t, y(t),'-')\n", + "# labels, title, legend\n", + "plt.xlabel('$t$'); plt.ylabel('$y$')\n", + "plt.legend(['EP','ER','$y(t)$'])\n", + "plt.grid(True)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Convergence\n", + "\n", + "On considère le problème de Cauchy\n", + "$$\\begin{cases}\n", + " y'(t) = - y(0.1-\\cos(t)), \\quad t>0 \\\\\n", + " y(0) = 1\n", + "\\end{cases}$$\n", + "\n", + "Resolvez ce problème par les méthodes d'Euler progressive et de\n", + "Heun sur l'intervalle $[0,12]$ avec\n", + "un pas de temps $h=0.4$.\n", + "\n", + "La solution exacte est $y(t) = e^{-0.1t +\\sin(t)}$. \n", + "On remarque que la solution obtenue par la méthode de\n", + "Heun est beaucoup plus précise que celle d'Euler progressive.\n", + "Par ailleurs, on peut voir que si on réduit le pas de temps, la solution obtenue par la\n", + "méthode d'Euler progressive s'approche de la solution exacte." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "f = lambda t,y : (np.cos(t) - 0.1)*y\n", + "\n", + "tspan = [0,12]; y0 = 1;\n", + "h = 0.4; Nh = np.ceil((tspan[1] - tspan[0])/h).astype(int)\n", + "Nh = np.ceil(12/h).astype(int)\n", + "\n", + "t_EP, y_EP = forwardEuler(f, tspan, y0, Nh)\n", + "t_H, y_H = Heun(f, tspan, y0, Nh)\n", + "\n", + "plt.plot(t_EP, y_EP,'o-')\n", + "plt.plot(t_H, y_H,'o-')\n", + "\n", + "y = lambda t : np.exp(-0.1*t+np.sin(t))\n", + "t = np.linspace(tspan[0],tspan[1],100)\n", + "plt.plot(t, y(t),'-')\n", + "# labels, title, legend\n", + "plt.xlabel('$t$'); plt.ylabel('$y$')\n", + "plt.legend(['EP','Heun','$y(t)$'])\n", + "plt.grid(True)\n", + "plt.show()\n", + "print(\"Figure: Approximation par le méthodes de Euler Retrograde et Heun.\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "tspan = [0,12]; y0 = 1;\n", + "NhRange = [30, 50, 100, 500]\n", + "for Nh in NhRange :\n", + " t, y = forwardEuler(f,tspan,y0,Nh)\n", + " plt.plot(t, y,'-')\n", + "\n", + "yt = lambda t : np.exp(-0.1*t+np.sin(t))\n", + "t = np.linspace(tspan[0],tspan[1],100)\n", + "plt.plot(t, yt(t),':')\n", + "# labels, title, legend\n", + "plt.xlabel('$t_n$'); plt.ylabel('$u_n$')\n", + "plt.legend(NhRange+['$y(t)$'])\n", + "plt.title('$u_n\\\\approx u(t_n)$')\n", + "plt.grid(True)\n", + "plt.show()\n", + "print(\"Figure: Solutions obtenues par la méthode d'Euler progressive pour différents pas de temps.\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "On veut, maintenant, estimer l'ordre de convergence de ces deux\n", + "méthodes. Pour cela, on résout le problème avec\n", + "différents pas de temps et on compare les résultats obtenus à\n", + "l'instant $t=6$ avec la solution exacte." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "tspan=[0,6]\n", + "NhRange = [30, 50, 100, 500]\n", + "errEP = []\n", + "errH = []\n", + "# Solution at end time\n", + "y6 = yt(tspan[1])\n", + "for Nh in NhRange :\n", + " # Forward Euler\n", + " t, y = backwardEuler(f,tspan,y0,Nh)\n", + " # Error at the end of the simulation\n", + " errEP.append( np.abs(y6 - y[-1] ) )\n", + " \n", + " # Heun\n", + " [t, y] = Heun(f, tspan, y0, Nh);\n", + " # Error at the end of the simulation\n", + " errH.append( np.abs(y6 - y[-1] ) )\n", + "\n", + "h = (tspan[1] - tspan[0])/np.array(NhRange)\n", + "plt.loglog(h,errEP,'o-b',h,errH,'o-r')\n", + "plt.loglog(h,h*(errEP[0]/h[0]),':',h,(h**2*(errH[0]/h[0]**2)),':')\n", + "plt.xlabel('$h$'); plt.ylabel('$|y(6)-u_{N_h}|$')\n", + "plt.legend(['EP','Heun','$h$','$h^2$'])\n", + "plt.title('Decay of the error')\n", + "plt.grid(True)\n", + "plt.show()\n", + "\n", + "print(\"Figure: Erreurs en échelle logarithmique commises par les méthodes\" +\n", + " \" d'Euler progressive et de Heun dans le calcul de y(6).\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "La figure montre, en échelle logarithmique, les erreurs\n", + "commises par les deux méthodes en fonction de $h$.\n", + "On voit bien que la\n", + "méthode d'Euler progressive converge à l'ordre 1 tandis que celle de\n", + "Heun à l'ordre 2.\n" + ] + }, + { + "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 +}