From 14d015734c3366cd982fb80e6640fcfb091ea243 Mon Sep 17 00:00:00 2001
From: Olivier Cots <olivier.cots@irit.fr>
Date: Tue, 21 May 2024 14:47:11 +0200
Subject: [PATCH] foo
---
tp/rk-implicites.ipynb | 467 +++++++++++++++++++++++++++++++++++++++++
1 file changed, 467 insertions(+)
create mode 100644 tp/rk-implicites.ipynb
diff --git a/tp/rk-implicites.ipynb b/tp/rk-implicites.ipynb
new file mode 100644
index 0000000..565a768
--- /dev/null
+++ b/tp/rk-implicites.ipynb
@@ -0,0 +1,467 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "<img src=\"https://gitlab.irit.fr/toc/ens-n7/texCoursN7/-/raw/main/LOGO_INP_N7.png\" alt=\"N7\" height=\"80\"/>\n",
+ "\n",
+ "<img src=\"https://gitlab.irit.fr/toc/ens-n7/texCoursN7/-/raw/main/logo-insa.png\" alt=\"INSA\" height=\"80\"/>\n",
+ "\n",
+ "# Méthodes de Runge-Kutta implicites - TP Projet\n",
+ "\n",
+ "- Date : 2023-2024\n",
+ "- Durée approximative : inconnue"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "**Nom** : \"Entrez votre nom ici\"\n",
+ "\n",
+ "**Prénom** : \"Entrez votre prénom ici\""
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Rendu et consignes\n",
+ "\n",
+ "Une fois le travail terminé, vous enverrez directement le fichier `.ipynb` par email à l'adresse : `olivier.cots@toulouse-inp.fr`.\n",
+ "\n",
+ "- **Date limite du rendu : mercredi 15/11/2023 à 23h59.** Attention, à partir de 24h, 2 points est enlevé de la note finale toutes les 30 minutes.\n",
+ "- **Attention :** Le fichier doit être renommé de la façon suivante : `rk_implicites_NOM_Prenom.ipynb`. 4 points enlevés si le nom du fichier n'est pas respecté.\n",
+ "- **Documents autorisés :** vous pouvez utiliser votre polycopié et les TPs précédents."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 13,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "(-2, -2)"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "a = (1, 2)\n",
+ "b = (3, 4)\n",
+ "a.-b"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 17,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "f (generic function with 1 method)"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "f(x, y) = x+y "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 19,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "3"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "f(a...)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 16,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "2-element Vector{Int64}:\n",
+ " 1\n",
+ " 2"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "v = [a...]"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 23,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "0.0"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "using LinearAlgebra\n",
+ "k1 = [1, 2]\n",
+ "k2 = [3, 4]\n",
+ "ks = (k1, k2)\n",
+ "G(k1, k2) = (k1, k2)\n",
+ "norm(ks.-G(ks...))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Introduction\n",
+ "\n",
+ "Nous allons dans ce TP, implémenter quelques méthodes de Runge-Kutta **implicites** (voir **polycopié Section 8.2**) et étudier leur convergence. On considère un pas de temps $h$ uniforme. Une méthode à un pas implicite est convergente si pour toute solution $x(\\cdot, t_0, x_0)$ d'un problème de Cauchy, la suite approximante ${(x_i)}_i$ donnée par la méthode à un pas implicite vérifie \n",
+ "$$\n",
+ " \\max_{1 \\le i \\le N}\\, \\|{x(t_i, t_0, x_0) - x_i}\\| \\to 0 \n",
+ " \\quad\\text{quand}\\quad h \\to 0.\n",
+ "$$\n",
+ "\n",
+ "Si la convergence est d'ordre $p$ alors il existe une constante $C \\ge 0$ telle que l'**erreur globale** $E$ vérifie\n",
+ "\n",
+ "$$\n",
+ " E := \\max_{1 \\le i \\le N}\\, \\|{x(t_i, t_0, x_0) - x_i}\\| \\le C\\, h^p.\n",
+ "$$\n",
+ "\n",
+ "Faisons l'hypothèse que $E = M\\, h^p$ pour un certain $M \\ge 0$. En passant au logarithme, on obtient\n",
+ "\n",
+ "$$\n",
+ " \\log(E) = \\log(M) + p\\, \\log(h).\n",
+ "$$\n",
+ "\n",
+ "Nous en déduisons que si on trace $\\log(E)$ en fonction de $\\log(h)$, on doit obtenir une droite de pente $p$. C'est ce que nous allons vérifier dans ce TP."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# activate local project\n",
+ "using Pkg\n",
+ "Pkg.activate(\".\")\n",
+ "\n",
+ "# load packages\n",
+ "using LinearAlgebra\n",
+ "using Plots\n",
+ "using Plots.PlotMeasures\n",
+ "using Polynomials\n",
+ "\n",
+ "#\n",
+ "px = PlotMeasures.px;"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Fonctons auxiliaires \n",
+ "\n",
+ "# VOUS POUVEZ METTRE VOS FONCTIONS AUXILIAIRES ICI"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## L'exemple d'étude\n",
+ "\n",
+ "On s'intéresse (pour les exercices 1, 2 et 3) au problème de Cauchy\n",
+ "\n",
+ "$$\n",
+ " \\dot{x}(t) = (1-2t) x(t), \\quad x(0) = x_0 = 1\n",
+ "$$\n",
+ "\n",
+ "sur l'intervalle $[0, 3]$."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Définition du problème de Cauchy\n",
+ "f(t, x) = (1-2t)x # Second membre f(t, x)\n",
+ "x0 = 1.0 # Condition initiale\n",
+ "tspan = (0.0, 3.0); # Intervalle de temps"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Solution analytique\n",
+ "function sol(t)\n",
+ " return exp(t-t^2)\n",
+ "end;"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Estimation de la constante de Lipschitz de f sur [0, 3]\n",
+ "# Voir Théorème 8.2.2 pour l'utilité de cette estimation\n",
+ "function dfx(t)\n",
+ " return 1-2t\n",
+ "end \n",
+ "L = maximum(abs.(dfx.(range(0, stop=3, length=1000))))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## La méthode d'Euler implicite\n",
+ "\n",
+ "La méthode d'Euler implicite est donnée par :\n",
+ "\n",
+ "$$\n",
+ "\\left\\{\n",
+ "\\begin{array}{l}\n",
+ "x_{n+1} = x_n + h f(t_n + h, x_{n+1}) = G(x_{n+1}) \\\\\n",
+ "x_0 = x(t_0)\n",
+ "\\end{array}\n",
+ "\\right.\n",
+ "$$\n",
+ "\n",
+ "### Exercice 1\n",
+ "\n",
+ "1. Implémenter la méthode d'Euler implicite avec le point fixe (penser à voir le polycopié Section 8.2).\n",
+ "2. Pourquoi si $h \\ge 0.2$, l'algorithme du point fixe peut s'arrêter au bout du nombre d'itérations max et donc ne pas converger pour la méthode d'Euler implicite ?\n",
+ "3. Tracer la solution approchée et la solution exacte sur le même graphique pour différentes valeurs de $h$ que vous choisirez pour illustrer la convergence de la méthode.\n",
+ "4. Tracer l'erreur globale de la méthode d'Euler implicite en fonction de $h$ et vérifier que l'erreur est bien en $O(h)$.\n",
+ "\n",
+ "**Attention** : pour l'algorithme du point fixe, faites attention aux critères d'arrêts (il y en a 2) ! Voir votre polycopié Section 8.2. Vous fixerez la valeur de la tolérance à $10^{-6}$ et le nombre maximum d'itérations à $1000$."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 10,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "0.00447987015488227"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "# AJOUTER VOTRE CODE ICI\n",
+ "\n",
+ "# F(y) = y - G(y)\n",
+ "\n",
+ "G(y) = (y + sin(y))/3\n",
+ "\n",
+ "y = 1\n",
+ "y = G(y)\n",
+ "y = G(y)\n",
+ "y = G(y)\n",
+ "y = G(y)\n",
+ "y = G(y)\n",
+ "y = G(y)\n",
+ "y = G(y)\n",
+ "y = G(y)\n",
+ "y = G(y)\n",
+ "y = G(y)\n",
+ "y = G(y)\n",
+ "y = G(y)\n",
+ "y = G(y)\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 11,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "0.0014932950464876392"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "F(y) = y - G(y)\n",
+ "F(y)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": []
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## La méthode des trapèzes\n",
+ "\n",
+ "La méthode des trapèzes est donnée par le tableau de Butcher :\n",
+ "\n",
+ "$$\n",
+ " \\begin{array}{c | c c}\n",
+ " 0 & 0 & 0 \\\\[0.2em]\n",
+ " 1 & 1/2 & 1/2 \\\\[0.2em]\n",
+ " \\hline\n",
+ " & 1/2 & 1/2 \\\\\n",
+ " \\end{array}\n",
+ "$$\n",
+ "\n",
+ "### Exercice 2\n",
+ "\n",
+ "1. Implémenter la méthode des trapèzes avec le point fixe.\n",
+ "2. Tracer la solution approchée et la solution exacte sur le même graphique pour différentes valeurs de $h$ que vous choisirez pour illustrer la convergence de la méthode.\n",
+ "3. Tracer l'erreur globale de la méthode des trapèzes. Quel est l'ordre de convergence de la méthode des trapèzes ?"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# AJOUTER VOTRE CODE ICI"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## La méthode de Gauss à 2 étages\n",
+ "\n",
+ "La méthode de Gauss à 2 étages est donnée par le tableau de Butcher :\n",
+ "\n",
+ "$$\n",
+ " \\begin{array}{c | c c}\n",
+ " 1/2 - \\sqrt{3}/6 & 1/4 & 1/4 - \\sqrt{3}/6 \\\\[0.2em]\n",
+ " 1/2 + \\sqrt{3}/6 & 1/4 + \\sqrt{3}/6 & 1/4 \\\\[0.2em]\n",
+ " \\hline\n",
+ " & 1/2 & 1/2 \\\\\n",
+ " \\end{array}\n",
+ "$$\n",
+ "\n",
+ "### Exercice 3\n",
+ "\n",
+ "1. Implémenter la méthode de Gauss à 2 étages avec le point fixe.\n",
+ "2. Tracer la solution approchée et la solution exacte sur le même graphique pour différentes valeurs de $h$ que vous choisirez pour illustrer la convergence de la méthode.\n",
+ "3. Tracer l'erreur globale de la méthode de Gauss à 2 étages. Quel est l'ordre de convergence de la méthode de Gauss à 2 étages ?"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# AJOUTER VOTRE CODE ICI"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Un autre exemple\n",
+ "\n",
+ "On considère à partir de maintenant l'équation différentielle en dimension 2 :\n",
+ "\n",
+ "$$\n",
+ " \\dot{x}_1(t) = x_2(t), \\quad \\dot{x}_2(t) = - x_1(t).\n",
+ "$$\n",
+ "\n",
+ "On peut montrer facilement que la norme de $x(t) = (x_1(t), x_2(t))$ est constante le long des solutions :\n",
+ "\n",
+ "$$\n",
+ " \\frac{\\mathrm{d}}{\\mathrm{d} t} \\|x(t)\\|^2 = 2\\, \\left( x(t) \\,|\\, \\dot{x}(t) \\right) = 2 \\left( x_1(t) x_2(t) - x_2(t) x_1(t) \\right) = 0.\n",
+ "$$\n",
+ "\n",
+ "### Exercice 4\n",
+ "\n",
+ "On considère le problème de Cauchy associé de condition initiale $x_0 = (1, 0)$.\n",
+ "\n",
+ "1. Afficher l'approximation de la solution sur $[0, 10]$ pour les méthodes :\n",
+ "- Euler explicite ;\n",
+ "- Euler implicite ;\n",
+ "- Trapèzes ;\n",
+ "- Gauss à 2 étages.\n",
+ "2. Commentaires.\n",
+ "\n",
+ "**Attention :** vous ferez un affichage dans le plan $(x_1, x_2)$. Vous fixerez le nombre de pas à $N=100$."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# AJOUTER VOTRE CODE ICI"
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Julia 1.9.0",
+ "language": "julia",
+ "name": "julia-1.9"
+ },
+ "language_info": {
+ "file_extension": ".jl",
+ "mimetype": "application/julia",
+ "name": "julia",
+ "version": "1.9.0"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
--
GitLab