diff --git a/6.2 EDO - Stabilite et convergence.ipynb b/6.2 EDO - Stabilite et convergence.ipynb index 49ce24f..1eea3b0 100644 --- a/6.2 EDO - Stabilite et convergence.ipynb +++ b/6.2 EDO - Stabilite et convergence.ipynb @@ -1,208 +1,223 @@ { "cells": [ { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "tags": [] + }, "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 problè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": [ + "# 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": [ "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", "\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 }