diff --git a/tp/ode.ipynb b/tp/ode.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..aace7dd676a87a3f52dba7fffe198a48fe4f5137
--- /dev/null
+++ b/tp/ode.ipynb
@@ -0,0 +1,906 @@
+{
+ "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\n",
+ "\n",
+ "- Date : 2023\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\"); # 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\\\\\n",
+ "\\dot{x}_2(t)=-x_1(t)+3x_2(t)\\\\\n",
+ "x_1(0)=-9/25\\\\\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 la fonction $f$ et exécuter le code ci-dessous. Commentaires ?"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Remarque : \n",
+ "# un problème IVP 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 contrôle du pas\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",
+ "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",
+ "\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",
+ "\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=(850, 550))"
+ ]
+ },
+ {
+ "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$ la solution explose en temps fini (à $0.1$ près) inférieur à la valeur 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()\n",
+ "px = plot()\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=(900, 400))"
+ ]
+ },
+ {
+ "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",
+ "theta0 = pi/3\n",
+ "x0 = [theta0, 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",
+ "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 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",
+ " 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",
+ "end\n",
+ "\n",
+ "theta0 = pi-10*eps()\n",
+ "x0 = [theta0 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",
+ "\n",
+ "theta0 = pi+10*eps()\n",
+ "x0 = [theta0 0]\n",
+ "tf = 50\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",
+ "\n",
+ "# circulation case\n",
+ "for thetapoint0 in 0:1.:4 \n",
+ " tf = 10\n",
+ " tspan = (0.,tf)\n",
+ " x0 = [-pi thetapoint0] # thetapoint0 > 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",
+ "end\n",
+ "for thetapoint0 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",
+ " 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",
+ "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"
+ ]
+ },
+ {
+ "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=(850, 550))"
+ ]
+ },
+ {
+ "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=(850, 550))"
+ ]
+ },
+ {
+ "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 : les solutions sont-elles bornées ?\n",
+ "- Remplacer la vitesse angulaire pour avoir $\\omega = \\omega_0$. 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 : lim = 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))$."
+ ]
+ },
+ {
+ "attachments": {},
+ "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",
+ "\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",
+ "\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",
+ "\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 maintenant 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) + \\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
+}