"# 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",
"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",
"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 ?"
"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",
"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."
font-size:1em;"> <!-- il faut laisser une ligne vide après celle-ci -->
* Nom :
* Prénom :
</div>
%% Cell type:markdown id: tags:
<div style="width:95%;
margin:10px;
padding:8px;
background-color:#afa;
border:2px solid #bbffbb;
border-radius:10px;
font-weight:bold;
font-size:1.5em;
text-align:center;">
Rendu et consignes
</div>
Une fois le travail terminé, vous enverrez directement le fichier `.ipynb` par email à l'adresse : `olivier.cots@toulouse-inp.fr`.
<div style="width:95%;
margin:10px;
padding:8px;
background-color:#FAA299;
border:2px solid #BF381B;
border-radius:10px;
font-weight:bold;
font-size:1em;
text-align:center;">
Date limite du rendu : vendredi 15/11/2024 à 18h00.
</div>
-**Attention**, à partir de 18h01, 2 points sont enlevés de la note finale toutes les 30 minutes.
-**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é.
-**Documents autorisés :** vous pouvez utiliser votre polycopié et les TPs précédents.
%% Cell type:markdown id: tags:
<div style="width:95%;
margin:10px;
padding:8px;
background-color:#afa;
border:2px solid #bbffbb;
border-radius:10px;
font-weight:bold;
font-size:1.5em;
text-align:center;">
Introduction
</div>
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
$$
\max_{1 \le i \le N}\,\|{x(t_i, t_0, x_0) - x_i}\|\to 0
\quad\text{quand}\quad h \to 0.
$$
Si la convergence est d'ordre $p$ alors il existe une constante $C \ge 0$ telle que l'**erreur globale** $E$ vérifie
$$
E := \max_{1 \le i \le N}\,\|{x(t_i, t_0, x_0) - x_i}\|\le C\, h^p.
$$
Faisons l'hypothèse que $E = M\, h^p$ pour un certain $M \ge 0$. En passant au logarithme, on obtient
$$
\log(E) = \log(M) + p\,\log(h).
$$
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.
### Préambule
%% Cell type:code id: tags:
``` julia
# activate local project
usingPkg
Pkg.activate(".")
# load packages
usingLinearAlgebra
usingPlots
usingPlots.PlotMeasures
usingPolynomials
usingDualNumbers
#
px=PlotMeasures.px;
```
%% Cell type:code id: tags:
``` julia
# Fonctons auxiliaires
# VOUS POUVEZ METTRE VOS FONCTIONS AUXILIAIRES ICI
```
%% Cell type:markdown id: tags:
<div style="width:95%;
margin:10px;
padding:8px;
background-color:#afa;
border:2px solid #bbffbb;
border-radius:10px;
font-weight:bold;
font-size:1.5em;
text-align:center;">
L'exemple d'étude
</div>
On s'intéresse (pour les exercices 1, 2 et 3) au problème de Cauchy
$$
\dot{x}(t) = (1-2t) x(t), \quad x(0) = x_0 = 1
$$
sur l'intervalle $[0, 3]$.
### Modèle et solution
%% Cell type:code id: tags:
``` julia
# Définition du problème de Cauchy
f(t,x)=(1-2t)x# Second membre f(t, x)
x0=1.0# Condition initiale
tspan=(0.0,3.0);# Intervalle de temps
```
%% Cell type:code id: tags:
``` julia
# Solution analytique
function sol(t)
returnexp(t-t^2)
end;
```
%% Cell type:code id: tags:
``` julia
# Estimation de la constante de Lipschitz de f sur [0, 3]
# Voir Théorème 8.2.2 pour l'utilité de cette estimation
La méthode du point milieu implicite est donnée par le tableau de Butcher :
$$
\begin{array}{c | c }
1/2 & 1/2 \\[0.2em]
\hline
& 1 \\
\end{array}
$$
### Exercice 1
1. Implémenter la méthode du point milieu implicite avec le point fixe (penser à voir le polycopié Section 8.2).
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 ?
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.
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 ?
**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 id: tags:
``` julia
# AJOUTER VOTRE CODE ICI
```
%% Cell type:markdown id: tags:
<div style="width:95%;
margin:10px;
padding:8px;
background-color:#afa;
border:2px solid #bbffbb;
border-radius:10px;
font-weight:bold;
font-size:1.5em;
text-align:center;">
La méthode de Hammer & Hollingsworth
</div>
La méthode de Hammer & Hollingsworth est donnée par le tableau de Butcher :
$$
\begin{array}{c | c c}
0 & 0 & 0 \\[0.2em]
2/3 & 1/3 & 1/3 \\[0.2em]
\hline
& 1/4 & 3/4 \\
\end{array}
$$
### Exercice 2
1. Implémenter la méthode de Hammer & Hollingsworth avec le point fixe.
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.
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 id: tags:
``` julia
# AJOUTER VOTRE CODE ICI
```
%% Cell type:markdown id: tags:
<div style="width:95%;
margin:10px;
padding:8px;
background-color:#afa;
border:2px solid #bbffbb;
border-radius:10px;
font-weight:bold;
font-size:1.5em;
text-align:center;">
La méthode de Lobatto IIIA
</div>
La méthode de Lobatto IIIA est donnée par le tableau de Butcher :
$$
\begin{array}{c | c c c}
0 & 0 & 0 & 0 \\[0.2em]
1/2 & 5/24 & 1/3 & -1/24 \\[0.2em]
1 & 1/6 & 2/3 & 1/6 \\[0.2em]
\hline
& 1/6 & 2/3 & 1/6 \\
\end{array}
$$
### Exercice 3
1. Implémenter la méthode de Lobatto IIIA avec le point fixe.
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.
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 id: tags:
``` julia
# AJOUTER VOTRE CODE ICI
```
%% Cell type:markdown id: tags:
<div style="width:95%;
margin:10px;
padding:8px;
background-color:#afa;
border:2px solid #bbffbb;
border-radius:10px;
font-weight:bold;
font-size:1.5em;
text-align:center;">
Calcul de dérivée et équation différentielle
</div>
On considère le problème de Cauchy suivant :
$$
\dot{x}(t) = 1 + x^2(t), \quad x(0) = \tan(x_0).
$$
La solution exacte est donnée par $x(t, x_0) = \tan(t + x_0)$.
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$.
Pour retrouver la taille de la discrétisation, on pourra écrire `N = Int(round(tf/h))`.
### Exercice 4
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.
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 id: tags:
``` julia
# AJOUTER VOTRE CODE ICI
```
%% Cell type:markdown id: tags:
### Exercice 5
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.
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 id: tags:
``` julia
# On fournit ceci pour la gestion des nombres duaux
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.
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.