diff --git a/be/rk-implicites.ipynb b/be/rk-implicites.ipynb new file mode 100644 index 0000000000000000000000000000000000000000..483f16b34118e4d4d5bbc63ccc1e81211e4d5333 --- /dev/null +++ b/be/rk-implicites.ipynb @@ -0,0 +1,460 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "[<img src=\"https://gitlab.irit.fr/toc/etu-n7/controle-optimal/-/raw/master/ressources/Logo-toulouse-inp-N7.png\" alt=\"N7\" height=\"80\"/>](https://gitlab.irit.fr/toc/etu-n7/controle-optimal)\n", + "<img src=\"https://gitlab.irit.fr/toc/ens-n7/texCoursN7/-/raw/main/logo-insa.png\" alt=\"INSA\" height=\"80\" style=\"margin-left:50px\"/>\n", + "\n", + "# Méthodes de Runge-Kutta implicites - TP Projet\n", + "\n", + "- Date : 2024-2025\n", + "- Durée approximative : inconnue" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "<div style=\"width:95%;\n", + " margin:10px;\n", + " padding:8px;\n", + " background-color:#FAA299;\n", + " border:2px solid #BF381B;\n", + " border-radius:10px;\n", + " font-weight:bold;\n", + " font-size:1em;\"> <!-- il faut laisser une ligne vide après celle-ci -->\n", + "\n", + "* Nom :\n", + "* Prénom :\n", + "</div>" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "<div style=\"width:95%;\n", + " margin:10px;\n", + " padding:8px;\n", + " background-color:#afa;\n", + " border:2px solid #bbffbb;\n", + " border-radius:10px;\n", + " font-weight:bold;\n", + " font-size:1.5em;\n", + " text-align:center;\">\n", + "Rendu et consignes\n", + "</div>\n", + "\n", + "Une fois le travail terminé, vous enverrez directement le fichier `.ipynb` par email à l'adresse : `olivier.cots@toulouse-inp.fr`.\n", + "\n", + "<div style=\"width:95%;\n", + " margin:10px;\n", + " padding:8px;\n", + " background-color:#FAA299;\n", + " border:2px solid #BF381B;\n", + " border-radius:10px;\n", + " font-weight:bold;\n", + " font-size:1em;\n", + " text-align:center;\">\n", + "Date limite du rendu : vendredi 15/11/2024 à 18h00.\n", + "</div>\n", + "\n", + "- **Attention**, à partir de 18h01, 2 points sont enlevés 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": "markdown", + "metadata": {}, + "source": [ + "<div style=\"width:95%;\n", + " margin:10px;\n", + " padding:8px;\n", + " background-color:#afa;\n", + " border:2px solid #bbffbb;\n", + " border-radius:10px;\n", + " font-weight:bold;\n", + " font-size:1.5em;\n", + " text-align:center;\">\n", + "Introduction\n", + "</div>\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.\n", + "\n", + "### Préambule" + ] + }, + { + "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", + "using DualNumbers\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": [ + "<div style=\"width:95%;\n", + " margin:10px;\n", + " padding:8px;\n", + " background-color:#afa;\n", + " border:2px solid #bbffbb;\n", + " border-radius:10px;\n", + " font-weight:bold;\n", + " font-size:1.5em;\n", + " text-align:center;\">\n", + "L'exemple d'étude\n", + "</div>\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]$.\n", + "\n", + "### Modèle et solution" + ] + }, + { + "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": [ + "<div style=\"width:95%;\n", + " margin:10px;\n", + " padding:8px;\n", + " background-color:#afa;\n", + " border:2px solid #bbffbb;\n", + " border-radius:10px;\n", + " font-weight:bold;\n", + " font-size:1.5em;\n", + " text-align:center;\">\n", + "La méthode du point milieu implicite\n", + "</div>\n", + "\n", + "La méthode du point milieu implicite est donnée par le tableau de Butcher :\n", + "\n", + "$$\n", + " \\begin{array}{c | c }\n", + " 1/2 & 1/2 \\\\[0.2em]\n", + " \\hline\n", + " & 1 \\\\\n", + " \\end{array}\n", + "$$\n", + "\n", + "### Exercice 1\n", + "\n", + "1. Implémenter la méthode du point milieu implicite avec le point fixe (penser à voir le polycopié Section 8.2).\n", + "2. Pourquoi si $h \\ge 0.4$, 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 du point milieu 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 du point milieu implicite en fonction de $h$. Quel est l'ordre de convergence de la méthode du point milieu implicite ?\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": null, + "metadata": {}, + "outputs": [], + "source": [ + "# AJOUTER VOTRE CODE ICI" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "<div style=\"width:95%;\n", + " margin:10px;\n", + " padding:8px;\n", + " background-color:#afa;\n", + " border:2px solid #bbffbb;\n", + " border-radius:10px;\n", + " font-weight:bold;\n", + " font-size:1.5em;\n", + " text-align:center;\">\n", + "La méthode de Hammer & Hollingsworth\n", + "</div>\n", + "\n", + "La méthode de Hammer & Hollingsworth est donnée par le tableau de Butcher :\n", + "\n", + "$$\n", + " \\begin{array}{c | c c}\n", + " 0 & 0 & 0 \\\\[0.2em]\n", + " 2/3 & 1/3 & 1/3 \\\\[0.2em]\n", + " \\hline\n", + " & 1/4 & 3/4 \\\\\n", + " \\end{array}\n", + "$$\n", + "\n", + "### Exercice 2\n", + "\n", + "1. Implémenter la méthode de Hammer & Hollingsworth 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 Hammer & Hollingsworth. Quel est l'ordre de convergence de la méthode de Hammer & Hollingsworth ?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# AJOUTER VOTRE CODE ICI" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "<div style=\"width:95%;\n", + " margin:10px;\n", + " padding:8px;\n", + " background-color:#afa;\n", + " border:2px solid #bbffbb;\n", + " border-radius:10px;\n", + " font-weight:bold;\n", + " font-size:1.5em;\n", + " text-align:center;\">\n", + "La méthode de Lobatto IIIA\n", + "</div>\n", + "\n", + "La méthode de Lobatto IIIA est donnée par le tableau de Butcher :\n", + "\n", + "$$\n", + " \\begin{array}{c | c c c}\n", + " 0 & 0 & 0 & 0 \\\\[0.2em]\n", + " 1/2 & 5/24 & 1/3 & -1/24 \\\\[0.2em]\n", + " 1 & 1/6 & 2/3 & 1/6 \\\\[0.2em]\n", + " \\hline\n", + " & 1/6 & 2/3 & 1/6 \\\\\n", + " \\end{array}\n", + "$$\n", + "\n", + "### Exercice 3\n", + "\n", + "1. Implémenter la méthode de Lobatto IIIA 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 Lobatto IIIA. Quel est l'ordre de convergence de la méthode de Lobatto IIIA ?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# AJOUTER VOTRE CODE ICI" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "<div style=\"width:95%;\n", + " margin:10px;\n", + " padding:8px;\n", + " background-color:#afa;\n", + " border:2px solid #bbffbb;\n", + " border-radius:10px;\n", + " font-weight:bold;\n", + " font-size:1.5em;\n", + " text-align:center;\">\n", + "Calcul de dérivée et équation différentielle\n", + "</div>\n", + "\n", + "On considère le problème de Cauchy suivant :\n", + "\n", + "$$\n", + " \\dot{x}(t) = 1 + x^2(t), \\quad x(0) = \\tan(x_0).\n", + "$$\n", + "\n", + "La solution exacte est donnée par $x(t, x_0) = \\tan(t + x_0)$. \n", + "\n", + "On considère la méthode du point milieu implicite et on note $x_{h}(x_0)$ l'approximation de la solution à l'instant finale $t_f = \\pi/4$, pour un pas de temps $h$.\n", + "Pour retrouver la taille de la discrétisation, on pourra écrire `N = Int(round(tf/h))`.\n", + "\n", + "### Exercice 4\n", + "\n", + "1. Ecrire une fonction Julia `dx_by_finite_diff(x0, h; η)` qui fournit la dérivée de $x_h(x_0)$ par rapport à $x_0$, par différences finies avant, où `η` est le pas de la méthode des différences finies. On pourra écrire une fonction auxiliaire pour calculer la solution approchée.\n", + "2. Fixons $x_0 = 0.5$. Pour $h \\in \\{ 1e^{-1}, 1e^{-2}, 1e^{-3}, 1e^{-4}, 1e^{-5} \\}$, tracer sur le même graphique (en échelle log) l'erreur de la dérivée de $x_h(x_0)$ obtenue par différences finies par rapport à la dérivée exacte en fonction de $\\eta$. Commentez." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# AJOUTER VOTRE CODE ICI" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Exercice 5\n", + "\n", + "1. Ecrire une fonction Julia `dx_by_dual(x0, h)` qui fournit la dérivée de $x_h(x_0)$ par rapport à $x_0$, par différentiation automatique par les nombres duaux.\n", + "2. Fixons $x_0 = 0.5$. Tracer sur le même graphique (en échelle log), en fonction de $h$, l'erreur de la dérivée de $x_h(x_0)$ obtenue par la différentiation automatique par rapport à la dérivée exacte et l'erreur de la solution approchée $x_h(x_0)$ par rapport à la solution exacte. Commentez." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# On fournit ceci pour la gestion des nombres duaux\n", + "\n", + "# creation of dual number ε\n", + "import Base.*\n", + "*(e::Function, x::Union{Number, Vector{<:Number}}) = e(x)\n", + "*(x::Union{Number, Vector{<:Number}}, e::Function) = e(x)\n", + "ε(x=1) = begin \n", + " if x isa Number\n", + " return Dual.(0.0, x)\n", + " else\n", + " return Dual.(zeros(length(x)), x)\n", + " end\n", + "end\n", + "em = ε\n", + "dual(x::Union{Dual, Vector{<:Dual}}) = dualpart.(x)\n", + "real(x::Union{Dual, Vector{<:Dual}}) = realpart.(x);\n", + "\n", + "# example of a dual number:\n", + "1+2ε" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# AJOUTER VOTRE CODE ICI" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Exercice 6\n", + "\n", + "1. Ecrire une fonction Julia `dx_by_eq_var(x0, h)` qui fournit la dérivée de $x_h(x_0)$ par rapport à $x_0$, par les équations variationnelles.\n", + "2. Fixons $x_0 = 0.5$. Tracer sur le même graphique (en échelle log), en fonction de $h$, l'erreur de la dérivée de $x_h(x_0)$ obtenue par les équations variationnelles par rapport à la dérivée exacte et l'erreur de la solution approchée $x_h(x_0)$ par rapport à la solution exacte. Commentez." + ] + }, + { + "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 +}