{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# L3 Analyse Numérique – TP7\n", "\n", "[Dequay A](mailto:antoine.dequay@ens-rennes.fr) &\n", "[Le Barbenchon P](mailto:pierre.le-barbenchon@ens-rennes.fr). TP ENS Rennes\n", "\n", "[Boutin B](mailto:benjamin.boutin@univ-rennes1.fr). Cours et TP Université de Rennes 1 - UFR Mathématiques\n", "\n", "Dans ce TP, on étudiera des méthodes numériques de calcul de solutions d'équations différentielles.\n", "\n", "- Exercice 1 : *Schéma d'Euler explicite*\n", "- Exercice 2 : *Méthodes symplectiques*\n", "\n", "Contrairement au TP précédent, **on n'utilisera pas la bibliothèque `scipy.integrate` et sa fonction `odeint`**. Le but est de programmer des méthodes qui la remplacent à la main.\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from matplotlib import pyplot as plt\n", "#%pip install ffmpeg # --> optionnel" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercice 1 : Schéma d'Euler explicite\n", "\n", "On considère un problème de Cauchy :\n", "\n", "$$ \\left \\{\n", "\\begin{aligned}\n", "\ty' &= f(t,y)\\\\\n", "\ty(t_0) &= y^0.\n", "\\end{aligned}\n", "\\right . $$\n", "\n", "où $f$ est une fonction définie sur un ouvert $U$ de $\\mathbb{R} \\times \\mathbb{R}^{d}$, à valeurs dans $\\mathbb{R}^d$, $(t_0,y^0) \\in U$ est la condition initiale.\n", "\n", "Sous l'hypothèse que $f$ est continue sur $U$ et *localement lipschitzienne* par rapport à sa seconde variable $y$, le **théorème de Cauchy-Lipschitz** garantit l'existence d'une unique solution maximale $y$, définie et continûment différentiable sur un intervalle $I$ ouvert contenant $t_0$.\n", "\n", "Pour simplifier, on se limitera au cas des équations autonomes (la fonction $f$ ne dépendant alors pas de $t$) vérifiant de plus que $f\\in\\mathcal{C}^1(\\mathbb{R}^d;\\mathbb{R}^d)$ est de différentielle bornée, si bien que toutes les solutions sont globales et définies de manière unique par leur valeur en $t_0=0$.\n", "\n", "En intégrant le problème de Cauchy entre deux instants $t$ et $t+h$ où $h>0$, on obtient la formulation:\n", "\n", "$$\n", "y(t+h)=y(t)+\\int_t^{t+h} f(s,y(s))\\mathrm{d} s.\n", "$$\n", "\n", "Les méthodes numériques permettant de résoudre le problème de Cauchy sont obtenues sur la base d'une formule d'approximation de cette intégrale par une méthode de quadrature.\n", "\n", "La méthode d'**Euler explicite** définit une suite d'instants $(t_i)$ et des approximations $y_i\\in\\mathbb{R}^d$ de la solution exacte $y(t_i)$. Elle se base sur une approximation de l'intégrale précédente par l'approximation des rectangles à gauche. La relation de récurrence prend alors la forme:\n", "$$\n", "\\begin{aligned}\n", "&y_{0} = y^0,\\\\\n", "&\\forall i\\in\\mathbb{N},\\quad y_{i+1} = y_{i}+h_{i}f(t_i,y_{i}),\\\\\n", "&\\forall i\\in\\mathbb{N},\\quad t_{i+1} = t_{i}+h_{i}.\n", "\\end{aligned}\n", "$$\n", "\n", "Typiquement, considérant un intervalle $[0,t_f]$, on fixe un entier $N$. Le pas $h_i$ est alors choisi égal à $h=t_f/N$ pour tout $i$, $0\\leq i\\leq N$, de sorte que l'approximation se fait aux points de la subdivision uniforme $0=t_0" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "from Animation import creer_animation\n", "\n", "def exemple():\n", " def f_exemple(t):#renvoie le graphe (X,Y) à tracer à la frame numéro t\n", " X = np.linspace(0, 2*np.pi, 100)\n", " Y = np.sin(X + t * 2 * np.pi / 100)\n", " return (X, Y)\n", "\n", " nb_frames = 100 #le nombre de frames à afficher\n", "\n", " creer_animation(f_exemple,\n", " nb_frames,\n", " \"exemple.gif\",#nom du fichier dans lequel on enregistre l'animation\n", " frame_time=20,#temps entre l'affichage de deux frames\n", " xmin = 0, xmax = 2*np.pi, ymin = -1, ymax = 1 #cadre dans lequel est dessiné l'animation\n", " )\n", "\n", "exemple()\n", "# dans le même dossier que votre fichier .ipynb, vous trouverez un gif de votre animation !" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Question 4\n", "\n", "Résoudre ce même problème avec la méthode d'Euler implicite :\n", "\n", "$$\\forall i \\in \\mathbb{N},\\quad y_{i+1} = y_{i}+h_{i}f(t_{i+1},y_{i+1}).$$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercice 2 : Méthodes symplectiques\n", "\n", "Étant donné une fonction $H\\in\\mathcal{C}^1(\\mathbb{R}^d\\times \\mathbb{R}^d;\\mathbb{R})$, on définit le système hamiltonien :\n", "$$ \\left\\{\n", "\\begin{aligned}\n", "&p'(t)=-\\nabla_q H(p,q),\\\\\n", "&q'(t)=\\nabla_p H(p,q).\n", "\\end{aligned}\n", "\\right. $$\n", "De nombreux systèmes de la physique prenne la forme précédente (cf. cours de Philippe Chartier).\n", "À titre d'exemple, citons le cas du mouvement d'une particule ponctuelle de masse $m$, de position $q\\in\\mathbb{R}^d$ et de moment $p\\in\\mathbb{R}^d$, soumise à une force $F$ dérivant d'un potentiel $V:\\mathbb{R}^d\\longrightarrow \\mathbb{R}$, avec donc $F=-\\nabla V(q)$. Le principe fondamental de la dynamique s'écrit alors $mq'' = -\\nabla V$, soit encore $mq'=p$ et $p'=-\\nabla V(q)$. Le Hamiltonien concerné est alors l'énergie totale du système physique : $H(p,q)=\\frac{1}{2m}\\vert p\\vert^2+V(q)$ et le système différentiel prend la forme ci-dessus.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Question 1\n", "\n", "Vérifier brièvement que pour toute solution maximale $(p,q)\\in\\mathcal{C}^1(\\mathbb{R},\\mathbb{R}^d\\times\\mathbb{R}^d)$, le hamiltonien $H(p,q)$ est constant au cours du temps. En déduire que si les ensembles de niveaux de $H$ sont compacts, alors les trajectoires sont toutes globales." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Question 2\n", "\n", "On se place dans le cas de l'oscillateur harmonique : $H(p,q)=p^2+q^2$ (avec $d=1$).\n", "Comment l'énergie du système évolue-t-elle numériquement avec les méthodes d'Euler explicite et d'Euler implicite ?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Question 3\n", "\n", "Que se passe-t-il avec la méthode d'\\textbf{Euler symplectique} suivante ?\n", "$$\n", "\\left\\{\n", "\\begin{aligned} \n", "\tp_{i+1}=p_i-h_i \\nabla_q H(p_{i},q_{i+1})\\\\\n", "\tq_{i+1}=q_i+h_i \\nabla_p H(p_{i},q_{i+1})\n", "\\end{aligned}\n", "\\right.\n", "$$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### **Application : le problème des ressorts**\n", "\n", "Un système de cinq masses assimilées à des points matériels $P\\in\\{A,B,C,D,E\\}$ tous de même masse $m=1$ sont reliés par sept ressorts comme l'indique le schéma suivant. Les points sont soumis aux seules forces issues de l'effet de ces ressorts, en particulier la gravité et les frottements éventuels ne sont pas considérés.\n", "\n", "![](TP07_schema.svg)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "Chaque ressort $r$ est caractérisé par une constante de raideur $k_r$ et une longueur d'équilibre $L^0_r$ ; sa longueur effective dans la configuration courant est notée $L_r$. Les valeurs indiquées sur le schéma correspondent aux données $(k_r,L_r^0)$. Les points de base $A=(0,0)$ et $B=(3,0)$ sont fixés dans le plan $\\mathbb{R}^2$. L'énergie potentielle et l'énergie cinétique du système sont respectivement :\n", "$$\n", "E_{\\mathrm{pot}} = \\sum_{r} \\dfrac{1}{2} k_r (L^0_r-L_r)^2,\\qquad\n", "E_{\\mathrm{kin}} = \\sum_{P} \\dfrac{1}{2} m \\left\\vert \\dfrac{dP}{dt}\\right\\vert^2.\n", "$$\n", "Le système est hamiltonien avec :\n", "$$\n", "\\begin{aligned}\n", "&q={}^t(x_C,y_C,x_D,y_D,x_E,y_E),\\\\\n", "&p={}^t(\\dot x_C,\\dot y_C,\\dot x_D,\\dot y_D,\\dot x_E,\\dot y_E),\\\\\n", "&H(p,q) = E_{\\mathrm{pot}}+E_{\\mathrm{kin}}.\n", "\\end{aligned}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Nous vous fournissons des fonctions pré-programmées qui pourront être librement utilisées dans `RessortLibrary.py`." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Question 3\n", "\n", "Déterminer par une méthode d'Euler explicite, la configuration du système à l'instant $t_f=20$ avec un pas $h=0.1$, pour la donnée de Cauchy $p^0=0_{\\mathbb{R}^6}$ et $q^0={}^t(0.5,1.5,3.5,1.5,1,-0.5)$." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "from RessortLibrary import *\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Question 4\n", "Tracer l'évolution de l'énergie potentielle, de l'énergie cinétique et de l'énergie totale du système au cours du temps, pour la solution numérique précédemment obtenue." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Question 5\n", "\n", "Reprendre ces deux questions en utilisant cette fois-ci une méthode d'Euler symplectique." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Question 6\n", "\n", "On ajoute une force de frottement de la forme $-\\alpha p$ dans la seconde équation, avec $\\alpha>0$. Que donnent les deux méthodes précédentes (configuration géométrique et énergie totale) ? On pourra résoudre sur $[0,100]$ avec un pas $h=0.1$ et pour un coefficient de frottement $\\alpha=0.1$.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "interpreter": { "hash": "e7370f93d1d0cde622a1f8e1c04877d8463912d04d973331ad4851f04de6915a" }, "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.6.5" } }, "nbformat": 4, "nbformat_minor": 2 }