diff --git a/tp/edo-general.ipynb b/tp/edo-general.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..1f8c20ded18e824b6c3433b4443d824716a101c0
--- /dev/null
+++ b/tp/edo-general.ipynb
@@ -0,0 +1,702 @@
+{
+ "cells": [
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "<img src=\"https://gitlab.irit.fr/toc/ens-n7/texCoursN7/-/raw/main/LOGO_INP_N7.png\" alt=\"N7\" height=\"100\"/>\n",
+ "\n",
+ "<img src=\"https://gitlab.irit.fr/toc/ens-n7/texCoursN7/-/raw/main/logo-insa.png\" alt=\"INSA\" height=\"100\"/>\n",
+ "\n",
+ "# Intégration numérique\n",
+ "\n",
+ "- Date : 2023-2024\n",
+ "- Durée approximative : 1h"
+ ]
+ },
+ {
+ "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": [
+ "# activate local project\n",
+ "using Pkg\n",
+ "Pkg.activate(\".\")\n",
+ "\n",
+ "# load packages\n",
+ "using DifferentialEquations\n",
+ "using ForwardDiff\n",
+ "using LinearAlgebra\n",
+ "using Plots"
+ ]
+ },
+ {
+ "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\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=(900, 600))"
+ ]
+ },
+ {
+ "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",
+ " \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, idxs=(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, idxs=(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, idxs=(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, idxs=(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, 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"
+ ]
+ },
+ {
+ "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",
+ " \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, idxs=(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, idxs=(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, idxs=(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=(900, 600))"
+ ]
+ },
+ {
+ "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 (pour les étuidants N7 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=(900, 600))"
+ ]
+ },
+ {
+ "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."
+ ]
+ },
+ {
+ "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",
+ "ω₀ = 1/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=(900, 600))\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 intervient explicitement, le système est non autonome.\n",
+ "\"\"\"\n",
+ "function Lorentz(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",
+ "σ = 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",
+ " # 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",
+ "\n",
+ " # \n",
+ " duration = 10\n",
+ " fps = 10\n",
+ " frames = fps*duration\n",
+ " step = length(sol.t) ÷ frames\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()"
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Julia 1.9.0",
+ "language": "julia",
+ "name": "julia-1.9"
+ },
+ "language_info": {
+ "file_extension": ".jl",
+ "mimetype": "application/julia",
+ "name": "julia",
+ "version": "1.9.0"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}