diff --git a/tp/ode.ipynb b/tp/ode.ipynb
index f0c3a7e244c8a7d306dd902e0c61ba9611892bcc..fe834d689f562688daf4aca1b257d901b6c3d629 100644
--- a/tp/ode.ipynb
+++ b/tp/ode.ipynb
@@ -9,7 +9,7 @@
"\n",
"# Intégration numérique et systèmes dynamiques\n",
"\n",
- "- Date : 2023\n",
+ "- Date : 2025\n",
"- Durée approximative : 1 séance (à finir à la maison)"
]
},
@@ -36,7 +36,7 @@
"using ForwardDiff\n",
"using LinearAlgebra\n",
"using Plots\n",
- "include(\"utils.jl\"); # plot_traj!, plot_flow!"
+ "include(\"utils.jl\"); # importe les fonctions plot_traj!, plot_flow!"
]
},
{
@@ -48,9 +48,9 @@
"\n",
"Nous allons ici résoudre le problème à valeur initiale\n",
"$$(IVP)\\left\\{\\begin{array}{l}\n",
- "\\dot{x}_1(t)=x_1(t)+x_2(t)+\\sin t\\\\\n",
- "\\dot{x}_2(t)=-x_1(t)+3x_2(t)\\\\\n",
- "x_1(0)=-9/25\\\\\n",
+ "\\dot{x}_1(t)=x_1(t)+x_2(t)+\\sin t\\\\[0.5em]\n",
+ "\\dot{x}_2(t)=-x_1(t)+3x_2(t)\\\\[1em]\n",
+ "x_1(0)=-9/25\\\\[0.5em]\n",
"x_2(0)=-4/25,\n",
"\\end{array}\\right.\n",
"$$\n",
@@ -61,7 +61,7 @@
"x_1(t)&= (-1/25)(13\\sin t+9\\cos t)\\\\\n",
"x_2(t)&= (-1/25)(3\\sin t+4\\cos t).\n",
"\\end{align*}$$\n",
- "- Coder la fonction $f$ et exécuter le code ci-dessous. Commentaires ?"
+ "- Coder ci-dessous la fonction $f$ dans `exemple1` et exécuter le code. Retrouvez-vous la solution exacte ?"
]
},
{
@@ -71,7 +71,7 @@
"outputs": [],
"source": [
"# Remarque : \n",
- "# un problème IVP est une équation différentielle ordinaire \n",
+ "# un problème IVP (ou problème de Cauchy) est une équation différentielle ordinaire\n",
"# (EDO en français, ODE en anglais)\n",
"# avec en plus une condition initiale.\n",
"\n",
@@ -109,7 +109,10 @@
"cell_type": "markdown",
"metadata": {},
"source": [
+ "Le calcul précédent est fait à l'aide de la fonction `solve` de DifferentialEquations.jl. L'intégrateur utilisé est `Tsit5` qui est un intégrateur de Runge-Kutta d'ordre 5(4). Il est très efficace, or nous observons que la solution numérique est peu précise. Ceci vient des options par défaut de l'intégrateur, notamment celles concernant le contrôle de la longueur du pas de temps. Examinons ce point.\n",
+ "\n",
"## Le contrôle du pas\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",
@@ -122,10 +125,10 @@
"\n",
"Référence : [Hairer, Nørsett, Wanner (2008) Solving Ordinary Differential Equations I Nonstiff Problems](https://link.springer.com/book/10.1007/978-3-540-78862-1).\n",
"\n",
- "Dans `Julia`, on peut modifier les valeurs par défaut : reltol = 1.e-3 et abstol = 1.e-6\n",
- "(lorsque ces valeurs sont des scalaires, on considère les mêmes quantités pour toutes les composantes).\n",
+ "Lors de l'appel à la fonction `solve`, on peut modifier les valeurs des tolérances relative et absolue : reltol = 1.e-3 et abstol = 1.e-6 (lorsque ces valeurs sont des scalaires, on considère les mêmes quantités pour toutes les composantes).\n",
+ "\n",
+ "**Question**. Réaliser l'intégration numérique pour \n",
"\n",
- "Réaliser l'intégration numérique pour \n",
"- reltol = 1e-3 et abstol = 1e-6\n",
"- reltol = 1e-6 et abstol = 1e-9\n",
"- reltol = 1e-10 et abstol = 1e-15\n",
@@ -157,7 +160,7 @@
"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=(850, 550))"
+ "plot(p1, p2, p3, p4, size=(650, 350))"
]
},
{
@@ -169,7 +172,7 @@
"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) inférieur à la valeur 3."
+ "- Déterminer à partir de quelle valeur de $p\\ge 1$ (à $0.1$ près) la solution explose en temps fini, en un temps inférieur à 3."
]
},
{
@@ -198,8 +201,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",
@@ -227,7 +230,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))"
]
},
{
@@ -277,11 +280,11 @@
"m = 1\n",
"λ = [g, l, k, m] # paramètres constants\n",
"\n",
- "theta0 = pi/3\n",
- "x0 = [theta0, 0] # état initial\n",
+ "α0 = pi/3\n",
+ "x0 = [α0, 0] # état initial\n",
"\n",
"t0 = 0\n",
- "tf = 3*pi*sqrt(l/g)*(1 + theta0^2/16 + theta0^4/3072) # 2*approximation de la période\n",
+ "tf = 3*pi*sqrt(l/g)*(1 + α0^2/16 + α0^4/3072)\n",
"tspan = (t0, tf) # instants initial et terminal\n",
"\n",
"prob = ODEProblem(pendule, x0, tspan, λ, reltol = 1.e-8, abstol = 1.e-8) # définition du problème en Julia\n",
@@ -316,51 +319,50 @@
"\n",
"plt = plot(xlabel = \"x₁\", ylabel = \"x₂\", legend = false) # initialisation du plot\n",
"\n",
- "for theta0 in 0:(2*pi)/10:2*pi\n",
- " theta0_princ = theta0\n",
- " tf = 3*pi*sqrt(l/g)*(1 + theta0_princ^2/16 + theta0_princ^4/3072) # 2*approximation of the period\n",
- " tspan = (0.0,tf)\n",
- " x0 = [theta0 0]\n",
+ "for α0 in 0:(2*pi)/10:2*pi\n",
+ " tf = 3*pi*sqrt(l/g)*(1 + α0^2/16 + α0^4/3072) # 2*approximation of the period\n",
+ " tspan = (0, tf)\n",
+ " x0 = [α0 0]\n",
" prob = ODEProblem(pendule, x0, tspan, λ, reltol = 1.e-8, abstol = 1.e-8)\n",
" sol = solve(prob)\n",
- " plot!(plt, sol, vars=(1,2), color=\"blue\") # lw = linewidth\n",
+ " plot!(plt, sol, idxs=(1,2), color=\"blue\") # lw = linewidth\n",
"end\n",
"\n",
- "theta0 = pi-10*eps()\n",
- "x0 = [theta0 0]\n",
+ "α0 = pi-10*eps()\n",
+ "x0 = [α0 0]\n",
"tf = 50 # problem for tf=50 1/4 of the period!\n",
"tspan = (0.0,tf)\n",
"prob = ODEProblem(pendule, x0, tspan, λ, reltol = 1.e-8, abstol = 1.e-8)\n",
"sol = solve(prob)\n",
- "plot!(plt, sol, vars=(1,2), xlims = (-2*pi,4*pi), color=\"green\") # lw = linewidth\n",
+ "plot!(plt, sol, idxs=(1,2), xlims = (-2*pi,4*pi), color=\"green\") # lw = linewidth\n",
"\n",
- "theta0 = pi+10*eps()\n",
- "x0 = [theta0 0]\n",
+ "α0 = pi+10*eps()\n",
+ "x0 = [α0 0]\n",
"tf = 50\n",
- "tspan = (0.0,tf)\n",
+ "tspan = (0, tf)\n",
"prob = ODEProblem(pendule, x0, tspan, λ, reltol = 1.e-8, abstol = 1.e-8)\n",
"sol = solve(prob)\n",
- "plot!(plt, sol, vars=(1,2), xlims = (-2*pi,4*pi), color=\"green\") # lw = linewidth\n",
+ "plot!(plt, sol, idxs=(1,2), xlims = (-2*pi,4*pi), color=\"green\") # lw = linewidth\n",
"\n",
"# circulation case\n",
- "for thetapoint0 in 0:1.:4 \n",
+ "for α̇0 in 0:1.:4 \n",
" tf = 10\n",
- " tspan = (0.,tf)\n",
- " x0 = [-pi thetapoint0] # thetapoint0 > 0 so theta increases from -pi to ...\n",
+ " tspan = (0, tf)\n",
+ " x0 = [-pi α̇0] # α̇0 > 0 so theta increases from -pi to ...\n",
" prob = ODEProblem(pendule, x0, tspan, λ, reltol = 1.e-8, abstol = 1.e-8)\n",
" sol = solve(prob)\n",
- " plot!(plt, sol, vars=(1,2), color=\"red\") # lw = linewidth\n",
+ " plot!(plt, sol, idxs=(1,2), color=\"red\") # lw = linewidth\n",
"end\n",
- "for thetapoint0 in -4:1.:0\n",
+ "for α̇0 in -4:1.:0\n",
" tf = 10\n",
- " tspan = (0.,tf)\n",
- " x0 = [3*pi thetapoint0] # thetapoint0 < 0 so theta decreases from 3pi to ...\n",
+ " tspan = (0, tf)\n",
+ " x0 = [3*pi α̇0] # α̇0 < 0 so theta decreases from 3pi to ...\n",
" prob = ODEProblem(pendule, x0, tspan, λ, reltol = 1.e-8, abstol = 1.e-8)\n",
" sol = solve(prob)\n",
- " plot!(plt, sol, vars=(1,2), color=\"purple\") # lw = linewidth\n",
+ " 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=(850, 550))\n"
+ "plot!(plt, [-pi 0 pi 2*pi 3*pi], [0 0 0 0 0], seriestype=:scatter, color=:black)\n",
+ "plot!(plt, xlims = (-pi,3*pi), size=(650, 350))\n"
]
},
{
@@ -415,7 +417,7 @@
"\n",
"# Trajectoires intérieures\n",
"for x01 in -2:0.4:0\n",
- " x0 = [x01,0]\n",
+ " x0 = [x01, 0]\n",
" prob = ODEProblem(vdp, x0, tspan, reltol = 1.e-8, abstol = 1.e-8)\n",
" sol = solve(prob)\n",
" plot!(sol, vars=(1,2), color = \"blue\") \n",
@@ -423,7 +425,7 @@
"\n",
"# Trajectoires extérieures\n",
"for x02 in 2.5:0.5:4\n",
- " x0 = [0,x02]\n",
+ " x0 = [0, x02]\n",
" prob = ODEProblem(vdp, x0, tspan, reltol = 1.e-8, abstol = 1.e-8)\n",
" sol = solve(prob)\n",
" plot!(sol, vars=(1,2), color = \"green\")\n",
@@ -437,7 +439,7 @@
"\n",
"#\n",
"plot!([0], [0], seriestype=:scatter) # point d'équilibre\n",
- "plot!(xlims = (-2.5, 2.5), ylims = (-3, 5), size=(850, 550))"
+ "plot!(xlims = (-2.5, 2.5), ylims = (-3, 5), size=(650, 350))"
]
},
{
@@ -534,7 +536,7 @@
"# FIN A COMPLETER/MODIFIER L'AFFICHAGE\n",
"# ------------------------------\n",
"\n",
- "plot!(xlims = (0, 4), ylims = (0, 2.5), size=(850, 550))"
+ "plot!(xlims = (0, 4), ylims = (0, 2.5), size=(650, 350))"
]
},
{
@@ -548,8 +550,8 @@
"$$\\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.\n",
+ "- 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."
]
@@ -600,7 +602,7 @@
" plt = plot(xlabel = \"x₁\", ylabel = \"x₂\", zlabel = \"t\", legend = false)\n",
"\n",
" #\n",
- " ω == ω₀ ? lim = 50 : lim = 5\n",
+ " lim = ω == ω₀ ? 50 : 5\n",
" plot!(xlims = (-lim, lim), ylims = (-lim, lim), zlims = (t0, tf), size=(850, 550))\n",
"\n",
" # get x1, x2, t from sol\n",
@@ -752,7 +754,7 @@
"source": [
"## Flots en dimension 2\n",
"\n",
- "On considère le problème de contrôle optimal suivant:\n",
+ "On considère le problème de contrôle optimal suivant :\n",
"\n",
"$$ \\frac{1}{2}\\int_0^{1} u^2(t)\\,\\mathrm{d}t \\to \\min, $$\n",
"\n",
@@ -762,7 +764,7 @@
"\n",
"$$ x(0)=x_0=-1,\\quad x(1)=0. $$\n",
"\n",
- "Le but est ici de tracer les extrémales et le flot associé pour des points initiaux $z(0) = (x_0, p(0))$."
+ "Le but est ici de tracer les extrémales et le flot associé pour des points initiaux $z(0) = (x_0, p(0))$. On se rappelle que les extrémales sont normales dans ce cas, autrement dit $p^0 = -1$."
]
},
{
@@ -770,17 +772,17 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- "On note $z = (x, p)$ la paire état, état adoint. On note $F(t, z)$ le système hamiltonien donné par le PMP, c-a-d:\n",
+ "On note $z = (x, p)$ la paire état, état adjoint. On note $F(t, z)$ le système hamiltonien donné par les conditions nécessaires d'optimalité :\n",
"\n",
"$$\n",
"F(t, z) = \\left(\\frac{\\partial H}{\\partial p}(t, z, u(z)), -\\frac{\\partial H}{\\partial x}(t, z, u(z))\\right)\n",
"$$\n",
"\n",
- "où $u(z)$ est le contrôle maximisant. \n",
+ "où $u(z)$ est le contrôle solution de $\\partial_u H = 0$. \n",
"\n",
"- Codez ci-dessous la fonction $F$ et visualiser son flot. \n",
"- Visualisez l'extrémale la plus proche de la solution.\n",
- "- Qu'observez-vous sur la forme des courbes vertes, c-a-d des fronts d'ondes ? Pourquoi ?\n",
+ "- Qu'observez-vous sur la forme des courbes vertes (appellées fronts d'ondes) ? Pourquoi ?\n",
"\n",
"**Remarque.** On rappelle que les extrémales s'écrivent $x(t) = p_0\\, \\mathrm{sh}(t) + x_0\\, e^{-t}$, $p(t) = p_0\\, e^t$."
]
@@ -828,7 +830,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- "On maintenant considère le problème de contrôle optimal suivant:\n",
+ "On considère maintenant le problème de contrôle optimal suivant:\n",
"\n",
"$$ \\frac{1}{2}\\int_0^{1} u^2(t)\\,\\mathrm{d}t \\to \\min, $$\n",
"\n",