From 44ceb1fae0f552334a8646d7dd7b993415a78cf8 Mon Sep 17 00:00:00 2001
From: Olivier Cots <olivier.cots@irit.fr>
Date: Mon, 2 Sep 2024 12:43:03 +0200
Subject: [PATCH] up tp edo
---
tp/edo-general.ipynb | 176 ++++++++++++++++++++++++++++++++-----------
1 file changed, 134 insertions(+), 42 deletions(-)
diff --git a/tp/edo-general.ipynb b/tp/edo-general.ipynb
index d125db7..108f464 100644
--- a/tp/edo-general.ipynb
+++ b/tp/edo-general.ipynb
@@ -10,7 +10,7 @@
"\n",
"# Intégration numérique\n",
"\n",
- "- Date : 2023-2024\n",
+ "- Date : 2024-2025\n",
"- Durée approximative : 1h15"
]
},
@@ -19,12 +19,23 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- "## Introduction\n",
- "\n",
- "Pour la documentation du package `Julia` sur l'intégration numérique, voir\n",
- "[documentation de DifferentialEquations.jl](https://diffeq.sciml.ai/stable/).\n",
- "\n",
- "Il nous faut tout d'abord importer les bons packages de Julia."
+ "<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",
+ "Pour calculer numériquement des solutions d'équations différentielles ordinaires, nous allons utiliser le \n",
+ "package `OrdinaryDiffEq.jl`. Celui-ci est inclus dans un package plus volumineux `DifferentialEquations.jl`. \n",
+ "\n",
+ "Pour la documentation de ces packages `Julia` sur l'intégration numérique, voir\n",
+ "[la documentation de DifferentialEquations.jl](https://diffeq.sciml.ai/stable/)."
]
},
{
@@ -38,7 +49,7 @@
"Pkg.activate(\".\")\n",
"\n",
"# load packages\n",
- "using DifferentialEquations\n",
+ "using OrdinaryDiffEq\n",
"using ForwardDiff\n",
"using LinearAlgebra\n",
"using Plots"
@@ -49,7 +60,17 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- "## Exemple 1\n",
+ "<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",
+ "Exemple 1\n",
+ "</div>\n",
"\n",
"Nous allons ici résoudre le problème à valeur initiale\n",
"$$(IVP)\\left\\{\\begin{array}{l}\n",
@@ -84,7 +105,7 @@
" Second membre de l'IVP\n",
" x : vecteur d'état\n",
" λ : vecteur de paramètres\n",
- " t : variable de temps\n",
+ " t : variable de temps. Ici le temps intervient explicitement, le système est non autonome.\n",
"\"\"\"\n",
"function exemple1(x, λ, t)\n",
"\n",
@@ -96,7 +117,7 @@
" return xpoint\n",
"end\n",
"\n",
- "λ = [] # pas de paramètres\n",
+ "λ = Float64[] # pas de paramètres\n",
"\n",
"t0 = 0\n",
"tf = 10\n",
@@ -114,7 +135,18 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- "## Le contrôle du pas\n",
+ "<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",
+ "Le contrôle du pas\n",
+ "</div>\n",
+ "\n",
"Les bons intégrateurs numériques ont une mise à jour automatique du pas. Pour les \n",
"[méthodes de Runge-Kutta](https://fr.wikipedia.org/wiki/Méthodes_de_Runge-Kutta) \n",
"par exemple, sur un pas $[t_0, t_{1}]$ on calcule par 2 schémas différents deux solutions approchées à l'instant $t_{1}$ : $x_{1}$ et $\\hat{x}_{1}$ dont l'erreur locale est respectivement en $O(h^p)$ et $O(h^{p+1})$ (les erreurs globales sont elles en $O(h^{p-1})$ et $O(h^p)$). Ainsi on peut estimer l'erreur locale du schéma le moins précis à l'aide de la différence $x_{1}-\\hat{x}_{1}$. On peut alors estimer le pas de façon à avoir la norme de cette différence inférieure à une tolérance donnée `Tol` ou encore à avoir \n",
@@ -162,19 +194,29 @@
"p4 = plot(T,[(-1/25)*(13*sinT+9*cosT) (-1/25)*(3*sinT+4*cosT)], label = [\"x₁(t)\" \"x₂(t)\"], title = \"Solution exacte\")\n",
"\n",
"# affichage\n",
- "plot(p1, p2, p3, p4, size=(900, 600))"
+ "plot(p1, p2, p3, p4, size=(650, 350))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "## Sur l'explosion des solutions en temps fini\n",
+ "<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",
+ "Sur l'explosion des solutions en temps fini\n",
+ "</div>\n",
"\n",
"Nous allons ici résoudre le problème à valeur initiale $\\dot{x}(t) = 1 + x^p(t)$ avec $p \\ge 1$ et $x(0)=0$. \n",
"\n",
"- Coder la fonction $f$ qui permet d'écrire l'équation différentielle sous la forme $\\dot{x}(t)=f(t,x(t))$.\n",
- "- Déterminer à partir de quelle valeur de $p\\ge 1$ la solution explose en temps fini (à $0.1$ près)."
+ "- Déterminer à partir de quelle valeur de $p\\ge 1$ la solution explose en temps fini (à $0.1$ près) inférieur à la valeur 3."
]
},
{
@@ -190,7 +232,7 @@
" t : variable de temps. Ici le temps n'intervient pas explicitement, le système est autonome.\n",
"\"\"\"\n",
"function f(x, p, t)\n",
- " return 1+x^p\n",
+ " return 0.0\n",
"end\n",
"\n",
"#\n",
@@ -203,8 +245,8 @@
"xspan = 0:0.01:2\n",
"\n",
"#\n",
- "pf = plot()\n",
- "px = plot()\n",
+ "pf = plot(title=\"xᵖ\")\n",
+ "px = plot(title=\"x(t, p)\")\n",
"\n",
"#\n",
"p = 1\n",
@@ -230,7 +272,7 @@
"#\n",
"plot!(pf, xlims=(0, 2), ylims=(0, 4))\n",
"plot!(px, xlims=(t0, tf), ylims=( 0, 100))\n",
- "plot(pf, px, size=(900, 400))"
+ "plot(pf, px, size=(650, 350))"
]
},
{
@@ -238,7 +280,17 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- "## Pendule simple\n",
+ "<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",
+ "Pendule simple\n",
+ "</div>\n",
"\n",
"On s'intéresse ici au [pendule simple](https://fr.wikipedia.org/wiki/Pendule_simple). Les principes physiques de la mécanique classique donnent comme équation qui régit l'évolution du mouvement\n",
"\n",
@@ -364,7 +416,7 @@
" plot!(plt, sol, idxs=(1,2), color=\"purple\") # lw = linewidth\n",
"end\n",
"plot!(plt, [-pi 0 pi 2*pi 3*pi], [0 0 0 0 0], seriestype=:scatter)\n",
- "plot!(plt, xlims = (-pi,3*pi), size=(900, 600))\n"
+ "plot!(plt, xlims = (-pi,3*pi), size=(650, 350))\n"
]
},
{
@@ -372,7 +424,17 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- "## Modèle de Van der Pol\n",
+ "<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",
+ "Modèle de Van der Pol\n",
+ "</div>\n",
"\n",
"L'équation différentielle considérée est l'[équation de Van der Pol](https://fr.wikipedia.org/wiki/Oscillateur_de_Van_der_Pol)\n",
"$$(IVP)\\left\\{\\begin{array}{l}\n",
@@ -442,7 +504,7 @@
"\n",
"#\n",
"plot!([0], [0], seriestype=:scatter) # point d'équilibre\n",
- "plot!(xlims = (-2.5, 2.5), ylims = (-3, 5), size=(900, 600))"
+ "plot!(xlims = (-2.5, 2.5), ylims = (-3, 5), size=(650, 350))"
]
},
{
@@ -477,7 +539,17 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- "## Modèle proie-prédateur (modèle de Lotka-Volterra)\n",
+ "<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",
+ "Modèle proie-prédateur (modèle de Lotka-Volterra)\n",
+ "</div>\n",
"\n",
"L'équation différentielle considérée est donnée par les [équations de prédation de Lotka-Volterra](https://fr.wikipedia.org/wiki/Équations_de_prédation_de_Lotka-Volterra)\n",
"\n",
@@ -503,7 +575,7 @@
"**Questions.**\n",
"\n",
"- Le point $(0, 0)$ est clairement un point d'équilibre stable. Il existe un autre point d'équilibre, lequel ?\n",
- "- Afficher l'autre point d'équilibre ainsi que quelques trajectoires du système dans le plan de phase, qui rentrent dans les limites du plot ci-dessous. Qu'observez-vous sur les trajectoires."
+ "- Afficher l'autre point d'équilibre ainsi que quelques trajectoires du système dans le plan de phase, qui rentrent dans les limites du plot ci-dessous. Qu'observez-vous sur les trajectoires ?"
]
},
{
@@ -539,22 +611,34 @@
"# FIN A COMPLETER/MODIFIER L'AFFICHAGE\n",
"# ------------------------------\n",
"\n",
- "plot!(xlims = (0, 4), ylims = (0, 2.5), size=(900, 600))"
+ "plot!(xlims = (0, 4), ylims = (0, 2.5), size=(650, 350))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "## Phénomène de résonance\n",
+ "<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",
+ "Phénomène de résonance\n",
+ "</div>\n",
"\n",
"Soit $\\omega$ et $\\omega_0$ deux réels strictement positifs. On considère l'équation différentielle\n",
"\n",
"$$\\ddot{y}(t) + \\omega_0^2 y(t) = \\cos(\\omega t).$$\n",
"\n",
"- En prenant comme variable d'état $x=(x_1,x_2)=(y, \\dot{y})$, écrire la fonction $f$ qui permet d'écrire l'équation différentielle sous la forme $\\dot{x}(t) = f(t,x(t))$.\n",
- "- Coder cette fonction ci-dessous et exécuter le code : les solutions sont-elles bornées ?\n",
- "- Remplacer la vitesse angulaire pour avoir $\\omega = \\omega_0$. Commentaires."
+ "- Coder cette fonction ci-dessous et exécuter le code avec $\\omega_0 = \\omega/2$ : les solutions sont-elles bornées ?\n",
+ "- Remplacer la vitesse angulaire pour avoir $\\omega_0 = \\omega$. Commentaires.\n",
+ "\n",
+ "Remarque : voir la page Wikipedia sur le [phénomène de résonance](https://fr.wikipedia.org/wiki/Résonance#:~:text=Le%20phénomène%20de%20résonance%20n,en%20phase%20»%20avec%20cette%20dernière.) pour plus d'informations."
]
},
{
@@ -582,7 +666,7 @@
"\n",
"# parameters\n",
"ω = 1\n",
- "ω₀ = 1/2\n",
+ "ω₀ = ω/2\n",
"λ = [ω, ω₀]\n",
"\n",
"# integration times\n",
@@ -604,7 +688,7 @@
"\n",
" #\n",
" ω == ω₀ ? lim = 50 : lim = 5\n",
- " plot!(xlims = (-lim, lim), ylims = (-lim, lim), zlims = (t0, tf), size=(900, 600))\n",
+ " plot!(xlims = (-lim, lim), ylims = (-lim, lim), zlims = (t0, tf), size=(850, 550))\n",
"\n",
" # get x1, x2, t from sol\n",
" x1 = [sol.u[i][1] for i ∈ 1:length(sol.u)]\n",
@@ -629,7 +713,17 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- "## Attracteur de Lorenz\n",
+ "<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",
+ "Attracteur de Lorenz\n",
+ "</div>\n",
"\n",
"L'équation différentielle considérée est l'[équation de Lorenz](https://fr.wikipedia.org/wiki/Attracteur_de_Lorenz) donnée par\n",
"\n",
@@ -665,11 +759,12 @@
" Second membre de l'IVP\n",
" x : vecteur d'état\n",
" λ : vecteur de paramètres\n",
- " t : variable de temps. Ici le temps intervient explicitement, le système est non autonome.\n",
+ " t : variable de temps. Ici le temps n'intervient pas explicitement, le système est autonome.\n",
"\"\"\"\n",
"function Lorentz(x, λ, t)\n",
"\n",
" xpoint = similar(x)\n",
+ " σ, ρ, β = λ\n",
" \n",
" # A COMPLETER/MODIFIER\n",
" xpoint[:] = zeros(size(x))\n",
@@ -717,16 +812,8 @@
" x2 = [sol.u[i][2] for i ∈ 1:length(sol.u)]\n",
" x3 = [sol.u[i][3] for i ∈ 1:length(sol.u)]\n",
"\n",
- " # print min and max of each coordinate\n",
- " #println(\"x1 : \", minimum(x1), \" \", maximum(x1))\n",
- " #println(\"x2 : \", minimum(x2), \" \", maximum(x2))\n",
- " #println(\"x3 : \", minimum(x3), \" \", maximum(x3))\n",
- "\n",
" #\n",
- " plot!(xlims = (-20, 20), ylims = (-25, 30), zlims = (0, 50), size=(900, 600))\n",
- "\n",
- " # print length of the time grid\n",
- " #println(\"t : \", length(sol.t))\n",
+ " plot!(xlims = (-20, 20), ylims = (-25, 30), zlims = (0, 50), size=(850, 550))\n",
"\n",
" # \n",
" duration = 10\n",
@@ -734,6 +821,11 @@
" frames = fps*duration\n",
" step = length(sol.t) ÷ frames\n",
"\n",
+ " # return if there is nothing to compute\n",
+ " if (step == 0)\n",
+ " return;\n",
+ " end\n",
+ "\n",
" # plot the trajectory step by step and make a gif\n",
" i = 1\n",
" animation = @animate for j ∈ 2:step:length(sol.t)\n",
--
GitLab