{
"cells": [
{
"attachments": {},
"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=\"100\"/>](https://gitlab.irit.fr/toc/etu-n7/controle-optimal)\n",
"\n",
"# Intégration numérique et systèmes dynamiques\n",
"\n",
"- Date : 2025\n",
"- Durée approximative : 1 séance (à finir à la maison)"
]
},
{
"attachments": {},
"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."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"using DifferentialEquations\n",
"using ForwardDiff\n",
"using LinearAlgebra\n",
"using Plots\n",
"include(\"utils.jl\"); # importe les fonctions plot_traj!, plot_flow!"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"## Exemple 1\n",
"\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\\\\[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",
"\n",
"- Ecrire la fonction $f$ qui permet d'écrire l'équation différentielle sous la forme $\\dot{x}(t)=f(t,x(t))$.\n",
"- Vérifier que la solution de $(IVP)$ est\n",
"$$\\begin{align*}\n",
"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 ci-dessous la fonction $f$ dans `exemple1` et exécuter le code. Retrouvez-vous la solution exacte ?"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Remarque : \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",
"\"\"\"\n",
" 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",
"\"\"\"\n",
"function exemple1(x, λ, t)\n",
"\n",
" xpoint = similar(x)\n",
" \n",
" # A COMPLETER/MODIFIER\n",
" xpoint[:] = zeros(size(x))\n",
" \n",
" return xpoint\n",
"end\n",
"\n",
"λ = [] # pas de paramètres\n",
"\n",
"t0 = 0\n",
"tf = 10\n",
"tspan = (t0, tf) # intervalle d'intégration\n",
"\n",
"x0 = [-9/25, -4/25] # condition initiale\n",
"\n",
"prob = ODEProblem(exemple1, x0, tspan, λ) # définition du problème IVP\n",
"sol = solve(prob) # résolution du problème IVP\n",
"p1 = plot(sol, label = [\"x₁(t)\" \"x₂(t)\"], title = \"default options\") # affichage de la solution"
]
},
{
"attachments": {},
"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",
"$$\\frac{\\|x_{1}-\\hat{x}_{1}\\|}{\\mathrm{Tol}}<1.$$\n",
"\n",
"En pratique on considère des tolérences absolue et relative composante par composante et on définit : $\\mathrm{sc}_i = \\mathrm{abstol}_i + \\mathrm{reltol}_i \\times \\max(|x_{0i}|,|x_{1i}|)$. On calcule alors le pas de façon à avoir \n",
"\n",
"$$\\mathrm{err} = \\sqrt{\\frac{1}{n}\\sum_{i=1}^n\\left(\\frac{x_{1i}-\\hat{x}_{1i}}{\\mathrm{sc}_i}\\right)^2}<1.$$\n",
"\n",
"\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",
"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",
"- reltol = 1e-3 et abstol = 1e-6\n",
"- reltol = 1e-6 et abstol = 1e-9\n",
"- reltol = 1e-10 et abstol = 1e-15\n",
"\n",
"et afficher les résultats en ajoutant la solution."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# calcul avec les options par défaut\n",
"sol = solve(prob, reltol = 1e-3, abstol = 1e-6)\n",
"p1 = plot(sol, label = [\"x₁(t)\" \"x₂(t)\"], title = \"reltol=1e-3, abstol=1e-6\")\n",
"\n",
"# A COMPLETER/MODIFIER\n",
"#\n",
"p2 = plot()\n",
"\n",
"#\n",
"p3 = plot()\n",
"\n",
"# solution exacte\n",
"T = t0:(tf-t0)/100:tf\n",
"sinT = sin.(T) # opération vectorielle\n",
"cosT = cos.(T)\n",
"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=(650, 350))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Sur l'explosion des solutions en temps fini\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$ (à $0.1$ près) la solution explose en temps fini, en un temps inférieur à 3."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"\"\"\"\n",
" Second membre de l'IVP\n",
" x : vecteur d'état\n",
" p : scalaire\n",
" 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 0.0\n",
"end\n",
"\n",
"#\n",
"x0 = 0.0\n",
"t0 = 0.0\n",
"tf = 3\n",
"tspan = (t0, tf) # instants initial et terminal\n",
"\n",
"#\n",
"xspan = 0:0.01:2\n",
"\n",
"#\n",
"pf = plot(title=\"xᵖ\")\n",
"px = plot(title=\"x(t, p)\")\n",
"\n",
"#\n",
"p = 1\n",
"prob = ODEProblem(f, x0, tspan, p, reltol = 1.e-8, abstol = 1.e-8) # définition du problème en Julia\n",
"sol = solve(prob) # intégration numérique\n",
"pf = plot!(pf, xspan, x -> (f(x, p, t0)-1), label = \"p=$p\") # affichage du second membre\n",
"px = plot!(px, sol, label = \"p=$p\") # affichage de la solution\n",
"\n",
"#\n",
"####################################\n",
"p = 1 ######### A MODIFIER #########\n",
"####################################\n",
"prob = ODEProblem(f, x0, tspan, p, reltol = 1.e-8, abstol = 1.e-8) # définition du problème en Julia\n",
"sol = solve(prob) # intégration numérique\n",
"pf = plot!(pf, xspan, x -> (f(x, p, t0)-1), label = \"p=$p\") # affichage du second membre\n",
"px = plot!(px, sol, label = \"p=$p\") # affichage de la solution\n",
"\n",
"#\n",
"p = 2\n",
"prob = ODEProblem(f, x0, tspan, p, reltol = 1.e-8, abstol = 1.e-8) # définition du problème en Julia\n",
"sol = solve(prob) # intégration numérique\n",
"pf = plot!(pf, xspan, x -> (f(x, p, t0)-1), label = \"p=$p\") # affichage du second membre\n",
"px = plot!(px, sol, label = \"p=$p\") # affichage de la solution\n",
"\n",
"#\n",
"plot!(pf, xlims=(0, 2), ylims=(0, 4))\n",
"plot!(px, xlims=(t0, tf), ylims=( 0, 100))\n",
"plot(pf, px, size=(650, 350))"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"## Pendule simple\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",
"$$ ml^2\\ddot{\\alpha}(t)+mlg\\sin(\\alpha(t)) +k\\dot{\\alpha}(t)=0,$$\n",
"où $\\ddot{\\alpha}(t)$ désigne la dérivée seconde de l'angle $\\alpha$ par rapport au temps $t$. \n",
"\n",
"- En prenant comme variable d'état $x=(x_1,x_2)=(\\alpha, \\dot{\\alpha})$, é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. D'après-vous, observez-vous une oscillation ou une rotation du pendule ? \n",
"- Remplacer la vitesse angulaire initiale $\\dot{\\alpha}(0)$ par $2$. Commentaires."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"\"\"\"\n",
" Second membre de l'IVP\n",
" x : vecteur d'état\n",
" λ : vecteur de paramètres\n",
" t : variable de temps. Ici le temps n'intervient pas explicitement, le système est autonome.\n",
"\"\"\"\n",
"function pendule(x, λ, t)\n",
"\n",
" xpoint = similar(x)\n",
" g, l, k, m = λ\n",
" \n",
" # A COMPLETER/MODIFIER\n",
" xpoint[:] = zeros(size(x))\n",
" \n",
" return xpoint \n",
"end\n",
"\n",
"#\n",
"g = 9.81\n",
"l = 10\n",
"k = 0\n",
"m = 1\n",
"λ = [g, l, k, m] # paramètres constants\n",
"\n",
"α0 = pi/3\n",
"x0 = [α0, 0] # état initial\n",
"\n",
"t0 = 0\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",
"sol = solve(prob) # intégration numérique\n",
"\n",
"plot(sol, label = [\"x₁(t)\" \"x₂(t)\"]) # affichage de la solution"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"**Diagramme de phases.**\n",
"\n",
"- Exécutez le code ci-dessous et interprétez le graphique: où sont les oscillations, les rotations et les points d'équilibre stables et instables ?\n",
"- On considère le cas où $k=0.15$ (on introduit un frottement). Que se passe-t-il ?"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"#\n",
"g = 9.81\n",
"l = 1.5\n",
"k = 0\n",
"m = 1\n",
"λ = [g, l, k, m] # paramètres constants\n",
"\n",
"plt = plot(xlabel = \"x₁\", ylabel = \"x₂\", legend = false) # initialisation du plot\n",
"\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, idxs=(1,2), color=\"blue\") # lw = linewidth\n",
"end\n",
"\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, idxs=(1,2), xlims = (-2*pi,4*pi), color=\"green\") # lw = linewidth\n",
"\n",
"α0 = pi+10*eps()\n",
"x0 = [α0 0]\n",
"tf = 50\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, idxs=(1,2), xlims = (-2*pi,4*pi), color=\"green\") # lw = linewidth\n",
"\n",
"# circulation case\n",
"for α̇0 in 0:1.:4 \n",
" tf = 10\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, idxs=(1,2), color=\"red\") # lw = linewidth\n",
"end\n",
"for α̇0 in -4:1.:0\n",
" tf = 10\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, 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, color=:black)\n",
"plot!(plt, xlims = (-pi,3*pi), size=(650, 350))\n"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"## Modèle de Van der Pol\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",
"\\dot{x}_1(t)=x_2(t)\\\\\n",
"\\dot{x}_2(t)=(1-x_1^2(t))x_2(t)-x_1(t)\\\\\n",
"x_1(0)=2.00861986087484313650940188\\\\\n",
"x_2(0)=0\n",
"\\end{array}\\right.\n",
"$$\n",
"jusqu'au temps $t_f=T=6.6632868593231301896996820305$, où $T$ est la période de la solution.\n",
"\n",
"Codez le second membre de l'IVP ci-dessous et exécutez le code. Commentaires."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"\"\"\"\n",
" Second membre de l'IVP: modèle de Van der Pol\n",
" x : vecteur d'état\n",
" λ : vecteur de paramètres\n",
" t : variable de temps. Ici le temps n'intervient pas explicitement, le système est autonome.\n",
"\"\"\"\n",
"function vdp(x, λ, t)\n",
"\n",
" xpoint = similar(x)\n",
" \n",
" # A COMPLETER/MODIFIER\n",
" xpoint[:] = zeros(size(x))\n",
" \n",
" return xpoint\n",
"end\n",
"\n",
"#\n",
"t0 = 0\n",
"tf = 6.6632868593231301896996820305\n",
"tspan = (t0, tf)\n",
"\n",
"#\n",
"plot(xlabel = \"x₁\", ylabel = \"x₂\", legend = false)\n",
"\n",
"# Trajectoires intérieures\n",
"for x01 in -2:0.4: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",
"end\n",
"\n",
"# Trajectoires extérieures\n",
"for x02 in 2.5:0.5:4\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",
"end\n",
"\n",
"# Cycle limite périodique\n",
"x0 = [2.00861986087484313650940188,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 = \"red\", lw = 2.)\n",
"\n",
"#\n",
"plot!([0], [0], seriestype=:scatter) # point d'équilibre\n",
"plot!(xlims = (-2.5, 2.5), ylims = (-3, 5), size=(650, 350))"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"**Stabilité.**\n",
"\n",
"On calcule ci-dessous les valeurs propres de la matrice jacobienne du second membre de l'IVP au point d'équilibre $x_e = (0, 0)$. Commentaires (rappelez-vous votre cours d'[automatique](https://gitlab.irit.fr/toc/etu-n7/controle-optimal/-/raw/master/ressources/notes-autom-2021.pdf))."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Jacobienne de la fonction vdp \n",
"dvdp(x) = ForwardDiff.jacobian(x -> vdp(x, [], 0), x)\n",
"\n",
"# Point d'équilibre et la matrice Jacobienne en ce point\n",
"xe = [0, 0]\n",
"A = dvdp(xe)\n",
"\n",
"# Valeurs propres de la matrice Jacobienne\n",
"eigvals(A)"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"## Modèle proie-prédateur (modèle de Lotka-Volterra)\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",
"$$\\left\\{\\begin{array}{l}\n",
"\\dot{x}_1(t)= \\phantom{-}x_1(t) ( \\alpha - \\beta x_2(t)) \\\\[0.5em]\n",
"\\dot{x}_2(t)= -x_2(t) ( \\gamma - \\delta x_1(t))\n",
"\\end{array}\\right.\n",
"$$\n",
"\n",
"où\n",
"\n",
"- $t$ est le temps ;\n",
"- $x(t)$ est l'effectif des proies en fonction du temps ;\n",
"- $y(t)$ est l'effectif des prédateurs en fonction du temps ;\n",
"\n",
"Les paramètres suivants caractérisent les interactions entre les deux espèces :\n",
"\n",
"- $\\alpha$, taux de reproduction intrinsèque des proies (constant, indépendant du nombre de prédateurs) ;\n",
"- $\\beta$, taux de mortalité des proies dû aux prédateurs rencontrés ;\n",
"- $\\delta$, taux de reproduction des prédateurs en fonction des proies rencontrées et mangées ;\n",
"- $\\gamma$, taux de mortalité intrinsèque des prédateurs (constant, indépendant du nombre de proies) ;\n",
"\n",
"**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 ?"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# ------------------------------\n",
"# DEBUT A COMPLETER/MODIFIER LE SECOND MEMBRE DE L'IVP\n",
"# ...\n",
"\n",
"# FIN A COMPLETER/MODIFIER LE SECOND MEMBRE DE L'IVP\n",
"# ------------------------------\n",
"\n",
"# paramètres\n",
"α = 2/3\n",
"β = 4/3\n",
"γ = 1\n",
"δ = 1\n",
"λ = [α, β, γ, δ]\n",
"\n",
"t0 = 0\n",
"tf = 20\n",
"tspan = (t0, tf)\n",
"\n",
"plt = plot(xlabel = \"x₁\", ylabel = \"x₂\", legend = false)\n",
"\n",
"# ------------------------------\n",
"# DEBUT A COMPLETER/MODIFIER L'AFFICHAGE\n",
"# ...\n",
"\n",
"# FIN A COMPLETER/MODIFIER L'AFFICHAGE\n",
"# ------------------------------\n",
"\n",
"plot!(xlims = (0, 4), ylims = (0, 2.5), size=(650, 350))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Phénomène de résonance\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 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."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"\"\"\"\n",
" 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",
"\"\"\"\n",
"function resonance(x, λ, t)\n",
"\n",
" xpoint = similar(x)\n",
" \n",
" # A COMPLETER/MODIFIER\n",
" xpoint[:] = zeros(size(x))\n",
" \n",
" return xpoint\n",
"\n",
"end\n",
"\n",
"# parameters\n",
"ω = 1\n",
"ω₀ = ω/2\n",
"λ = [ω, ω₀]\n",
"\n",
"# integration times\n",
"t0 = 0\n",
"tf = 100\n",
"tspan = (t0, tf)\n",
"\n",
"# initial condition\n",
"x0 = [0, 0]\n",
"\n",
"# define the problem and solve it\n",
"prob = ODEProblem(resonance, x0, tspan, λ, reltol = 1.e-8, abstol = 1.e-8)\n",
"sol = solve(prob)\n",
"\n",
"function makegif()\n",
"\n",
" #\n",
" plt = plot(xlabel = \"x₁\", ylabel = \"x₂\", zlabel = \"t\", legend = false)\n",
"\n",
" #\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",
" x1 = [sol.u[i][1] for i ∈ 1:length(sol.u)]\n",
" x2 = [sol.u[i][2] for i ∈ 1:length(sol.u)]\n",
" t = sol.t\n",
"\n",
" # plot the trajectory step by step and make a gif\n",
" i = 1\n",
" @gif for j ∈ 2:2:length(sol.t)\n",
" # plot the trajectory part from i to j\n",
" plot!(plt, x1[i:j], x2[i:j], t[i:j], color=:blue)\n",
" # update i \n",
" i = j\n",
" end\n",
"\n",
"end\n",
"\n",
"makegif()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Attracteur de Lorenz\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",
"$$\\left\\{\\begin{array}{l}\n",
"\\dot{x}_1(t)= \\phantom{-}\\sigma (x_2(t) - x_1(t)) \\\\[0.5em]\n",
"\\dot{x}_2(t)= \\rho\\, x_1(t) - x_2(t) - x_1(t)\\, x_3(t) \\\\[0.5em]\n",
"\\dot{x}_3(t)= x_1(t)\\, x_2(t) - \\beta x_3(t)\n",
"\\end{array}\\right.\n",
"$$\n",
"\n",
"où \n",
"\n",
"- $t$ est le temps ;\n",
"- $x(t)$ est proportionnel à l'intensité du mouvement de convection ;\n",
"- $y(t)$ est proportionnel à la différence de température entre les courants ascendants et descendants ;\n",
"- $z(t)$ est proportionnel à la déviation du profil vertical de température par rapport à la valeur linéaire de référence.\n",
"- $\\sigma$, $\\beta$ et $r$ sont des paramètres réels positifs.\n",
"- $\\sigma$ est le nombre de Prandtl ;\n",
"- $\\rho$ est le nombre de \"Rayleigh\".\n",
"\n",
"On fixe souvent $\\sigma = 10$, $\\beta = 8/3$ et $\\rho = 28$ pour illustrer le phénomène de chaos.\n",
"\n",
"Compléter le code ci-dessous et exécuter le. Observer le phénomène de chaos. Vous pouvez modifer la condition initiale."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"\"\"\"\n",
" Second membre de l'IVP\n",
" x : vecteur d'état\n",
" λ : vecteur de paramètres\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",
" \n",
" return xpoint\n",
" \n",
"end\n",
"\n",
"# parameters\n",
"σ = 10\n",
"ρ = 28\n",
"β = 8/3\n",
"λ = [σ, ρ, β]\n",
"\n",
"# integration times\n",
"t0 = 0\n",
"tf = 100\n",
"tspan = (t0, tf)\n",
"\n",
"# initial condition\n",
"x0 = [0, 1, 0]\n",
"\n",
"# define the problem and solve it\n",
"prob = ODEProblem(Lorentz, x0, tspan, λ, reltol = 1.e-8, abstol = 1.e-8)\n",
"sol = solve(prob)\n",
"\n",
"function makegif()\n",
"\n",
" #\n",
" plt = plot(xlabel = \"x₁\", ylabel = \"x₂\", zlabel = \"x₃\", legend = false)\n",
"\n",
" # \n",
" plot!(plt, background_color = :royalblue4)\n",
"\n",
" # update the color of axis of the figure to white\n",
" plot!(plt, axiscolor = :white, tickfontcolor = :white, bordercolor = :white, fg_color_guide = :white)\n",
"\n",
" # plot the three equilibrium points\n",
" plot!(plt, [0], [0], [0], seriestype=:scatter, color=:red)\n",
" plot!(plt, [ sqrt(β*(ρ-1))], [ sqrt(β*(ρ-1))], [ρ-1], seriestype=:scatter, color=:red)\n",
" plot!(plt, [-sqrt(β*(ρ-1))], [-sqrt(β*(ρ-1))], [ρ-1], seriestype=:scatter, color=:red)\n",
"\n",
" # get x1, x2, t from sol\n",
" x1 = [sol.u[i][1] for i ∈ 1:length(sol.u)]\n",
" 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",
" #\n",
" plot!(xlims = (-20, 20), ylims = (-25, 30), zlims = (0, 50), size=(850, 550))\n",
"\n",
" # \n",
" duration = 10\n",
" fps = 10\n",
" 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",
" # plot the trajectory part from i to j\n",
" plot!(plt, x1[i:j], x2[i:j], x3[i:j], color=:gold2)\n",
" # update i \n",
" i = j\n",
" end\n",
"\n",
" gif(animation, \"lorentz.gif\", fps=fps)\n",
"\n",
"end\n",
"\n",
"makegif()"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"## Flots en dimension 2\n",
"\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",
"sous les contraintes\n",
"\n",
"$$ \\dot{x}(t) = -x(t) + u(t)$$\n",
"\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))$. On se rappelle que les extrémales sont normales dans ce cas, autrement dit $p^0 = -1$."
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"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 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 (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$."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# système hamiltonien\n",
"function F(z)\n",
" zpoint = similar(z)\n",
" \n",
" # A COMPLETER/MODIFIER\n",
" zpoint = zeros(size(z))\n",
"\n",
" return zpoint\n",
"end\n",
"\n",
"t0 = 0\n",
"tf = 1\n",
"\n",
"# initialisation du plot\n",
"plt = plot()\n",
"ymin = -0.2\n",
"ymax = 2.5\n",
"\n",
"# plot de trajectoires\n",
"z0 = [[-1, p0] for p0 ∈ range(ymin, ymax, length=50)]\n",
"plot_traj!(plt, F, z0, tf)\n",
"\n",
"# plot du flot\n",
"z0 = [[-1, p0] for p0 ∈ range(ymin, ymax, length=50)]\n",
"plot_flow!(plt, F, z0, tf)\n",
"\n",
"# affichages additionels\n",
"plot!(plt, xlims = (-1.2, 2), ylims = (ymin, ymax))\n",
"plot!(plt, [-1, -1], [ymin, ymax], linestyle = :dot, color = :black)\n",
"plot!(plt, [0, 0], [ymin, ymax], linestyle = :dot, color = :black, size=(850, 550), xlabel = \"x\", ylabel = \"p\")"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"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",
"sous les contraintes\n",
"\n",
"$$ \\dot{x}(t) = -x(t) + \\alpha x^2(t) + u(t)$$\n",
"\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))$. Visualisez le flot et le tp est terminé."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# paramètre\n",
"α = 1.5\n",
"\n",
"# système hamiltonien\n",
"function F(z)\n",
" zpoint = similar(z)\n",
" zpoint[1] = -z[1] + α*z[1]^2 + z[2]\n",
" zpoint[2] = (1 - 2*α*z[1])*z[2]\n",
" return zpoint\n",
"end\n",
"\n",
"t0 = 0\n",
"tf = 1\n",
"\n",
"# initialisation du plot\n",
"plt = plot()\n",
"ymin = -0.2\n",
"ymax = 2.5\n",
"\n",
"# plot de trajectoires\n",
"z0 = [[-1, p0] for p0 ∈ range(ymin, ymax, length=50)]\n",
"plot_traj!(plt, F, z0, tf)\n",
"\n",
"# plot du flot\n",
"z0 = [[-1, p0] for p0 ∈ range(ymin, ymax, length=100)]\n",
"plot_flow!(plt, F, z0, tf)\n",
"\n",
"# affichages additionels\n",
"plot!(plt, xlims = (-1.2, 2), ylims = (ymin, ymax))\n",
"plot!(plt, [-1, -1], [ymin, ymax], linestyle = :dot, color = :black)\n",
"plot!(plt, [0, 0], [ymin, ymax], linestyle = :dot, color = :black, size=(850, 550), xlabel = \"x\", ylabel = \"p\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 1.10.2",
"language": "julia",
"name": "julia-1.10"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.10.2"
}
},
"nbformat": 4,
"nbformat_minor": 4
}