"# 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",
"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",
"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",
"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",
Une fois le travail terminé, vous enverrez directement le fichier `.ipynb` par email à l'adresse : `olivier.cots@toulouse-inp.fr`.
-**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.
-**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:code id: tags:
``` julia
a=(1,2)
b=(3,4)
a.-b
```
%% Output
%% Cell type:code id: tags:
``` julia
f(x,y)=x+y
```
%% Output
%% Cell type:code id: tags:
``` julia
f(a...)
```
%% Output
%% Cell type:code id: tags:
``` julia
v=[a...]
```
%% Output
%% Cell type:code id: tags:
``` julia
usingLinearAlgebra
k1=[1,2]
k2=[3,4]
ks=(k1,k2)
G(k1,k2)=(k1,k2)
norm(ks.-G(ks...))
```
%% Output
%% Cell type:markdown id: tags:
## Introduction
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.
%% Cell type:code id: tags:
``` julia
# activate local project
usingPkg
Pkg.activate(".")
# load packages
usingLinearAlgebra
usingPlots
usingPlots.PlotMeasures
usingPolynomials
#
px=PlotMeasures.px;
```
%% Cell type:code id: tags:
``` julia
# Fonctons auxiliaires
# VOUS POUVEZ METTRE VOS FONCTIONS AUXILIAIRES ICI
```
%% Cell type:markdown id: tags:
## L'exemple d'étude
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]$.
%% 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
1. Implémenter la méthode d'Euler implicite avec le point fixe (penser à voir le polycopié Section 8.2).
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 ?
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 d'Euler implicite en fonction de $h$ et vérifier que l'erreur est bien en $O(h)$.
**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
# F(y) = y - G(y)
G(y)=(y+sin(y))/3
y=1
y=G(y)
y=G(y)
y=G(y)
y=G(y)
y=G(y)
y=G(y)
y=G(y)
y=G(y)
y=G(y)
y=G(y)
y=G(y)
y=G(y)
y=G(y)
```
%% Output
%% Cell type:code id: tags:
``` julia
F(y)=y-G(y)
F(y)
```
%% Output
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
## La méthode des trapèzes
La méthode des trapèzes est donnée par le tableau de Butcher :
$$
\begin{array}{c | c c}
0 & 0 & 0 \\[0.2em]
1 & 1/2 & 1/2 \\[0.2em]
\hline
& 1/2 & 1/2 \\
\end{array}
$$
### Exercice 2
1. Implémenter la méthode des trapèzes 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 des trapèzes. Quel est l'ordre de convergence de la méthode des trapèzes ?
%% Cell type:code id: tags:
``` julia
# AJOUTER VOTRE CODE ICI
```
%% Cell type:markdown id: tags:
## La méthode de Gauss à 2 étages
La méthode de Gauss à 2 étages est donnée par le tableau de Butcher :
1. Implémenter la méthode de Gauss à 2 étages 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 Gauss à 2 étages. Quel est l'ordre de convergence de la méthode de Gauss à 2 étages ?
%% Cell type:code id: tags:
``` julia
# AJOUTER VOTRE CODE ICI
```
%% Cell type:markdown id: tags:
## Un autre exemple
On considère à partir de maintenant l'équation différentielle en dimension 2 :