diff --git a/tutoriels/julia/ode-exemple-sol.ipynb b/tutoriels/julia/ode-exemple-sol.ipynb
deleted file mode 100644
index 1af496d28d3756a2373309193b4cdcba56222386..0000000000000000000000000000000000000000
--- a/tutoriels/julia/ode-exemple-sol.ipynb
+++ /dev/null
@@ -1,536 +0,0 @@
-{
- "cells": [
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "# Intégration numérique avec Julia, une introduction"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Introduction\n",
- "Il nous faut tout d'abord importer les bons packages de Julia"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "#import Pkg; \n",
- "#Pkg.add(\"DifferentialEquations\")\n",
- "#Pkg.add(\"Plots\")\n",
- "using DifferentialEquations\n",
- "using Plots\n",
- "#ENV[\"JULIA_GR_PROVIDER\"] = \"GR\""
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Pendule simple\n",
- "### Introduction\n",
- "On s'intéresse ici au pendule simple. Les principes physiques de la mécanique classique donnent comme équation qui régit l'évolution du mouvement\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",
- "On prend ici comme variable d'état qui décrit le système $x(t)=(x_1(t),x_2(t))=(\\alpha(t), \\dot{\\alpha}(t))$. Le système différentiel du premier ordre que l'on obtient s'écrit alors\n",
- "$$\n",
- "\\left\\{\\begin{array}{l}\n",
- "\\dot{x}_1(t) = x_2(t)\\\\\n",
- "\\dot{x}_2(t) = -\\frac{g}{l}\\sin(x_1(t))-kx_2(t)\\\\\n",
- "x_1(0) = x_{0,1}=\\alpha_0\\\\\n",
- "x_2(0) = x_{0,2}=\\dot{\\alpha}_0\n",
- "\\end{array}\\right.\n",
- "$$\n",
- "### Cas où la fonction second membre renvoie xpoint"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "function pendule(x,p,t)\n",
- "# second member of the IVP\n",
- "# x : state\n",
- "# real(2)\n",
- "# p : parameter vector\n",
- "# t : time, variable not use here\n",
- "# real\n",
- "# Output\n",
- "# xpoint : vector of velocity\n",
- "# same as x\n",
- " g = p[1]; l = p[2]; k = p[3]\n",
- " xpoint = similar(x)\n",
- " xpoint[1] = x[2]\n",
- " xpoint[2] = -(g/l)*sin(x[1]) - k*x[2]\n",
- " return xpoint\n",
- "end\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# Main\n",
- "#\n",
- "g = 9.81; l = 10; k = 0;\n",
- "p = [g,l,k] # constantes\n",
- "theta0 = pi/3\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) # instant initial et terminal\n",
- "x0 = [theta0,0] # état initial\n",
- "prob = ODEProblem(pendule,x0,tspan,p) # défini le problème en Julia\n",
- "sol = solve(prob) # réalise l'intégration numérique\n",
- "plot(sol, label = [\"x_1(t)\" \"x_2(t)\"]) "
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "### Cas où xpoint est le premier argument modifié de la fonction second membre"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "function pendule1!(xpoint,x,p,t)\n",
- "# second member of the IVP\n",
- "# x : state\n",
- "# real(2)\n",
- "# p : parameter vector\n",
- "# t : time, variable not use here\n",
- "# real\n",
- "# Output\n",
- "# xpoint : vector of velocity\n",
- "# same as x\n",
- " g = p[1]; l = p[2]; k = p[3]\n",
- " xpoint = similar(x)\n",
- " xpoint[1] = x[2]\n",
- " xpoint[2] = -(g/l)*sin(x[1]) - k*x[2]\n",
- " return nothing\n",
- "end"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# Main\n",
- "#\n",
- "g = 9.81; l = 10; k = 0;\n",
- "p = [g,l,k] # constantes\n",
- "theta0 = pi/3\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) # instant initial et terminal\n",
- "x0 = [theta0,0] # état initial\n",
- "prob = ODEProblem(pendule1!,x0,tspan,p) # défini le problème en Julia\n",
- "sol = solve(prob) # réalise l'intégration numérique\n",
- "plot(sol, label = [\"x_1(t)\" \"x_2(t)\"]) "
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "L'instruction `xpoint = similar(x)` **crée une variable locale xpoint** et on a perdu le paramètre formel `xpoint`."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "function pendule2!(xpoint,x,p,t)\n",
- "# second member of the IVP\n",
- "# x : state\n",
- "# real(2)\n",
- "# p : parameter vector\n",
- "# t : time, variable not use here\n",
- "# real\n",
- "# Output\n",
- "# xpoint : vector of velocity\n",
- "# same as x\n",
- " g = p[1]; l = p[2]; k = p[3]\n",
- " xpoint[1] = x[2]\n",
- " xpoint[2] = -(g/l)*sin(x[1]) - k*x[2]\n",
- " return nothing\n",
- "end"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# Main\n",
- "#\n",
- "g = 9.81; l = 10; k = 0.;\n",
- "p = [g,l,k] # constantes\n",
- "theta0 = pi/3\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) # instant initial et terminal\n",
- "x0 = [theta0,0] # état initial\n",
- "prob = ODEProblem(pendule2!,x0,tspan,p) # défini le problème en Julia\n",
- "sol = solve(prob) # réalise l'intégration numérique\n",
- "plot(sol, label = [\"x_1(t)\" \"x_2(t)\"]) "
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Diagrammes de phases"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# Main\n",
- "#\n",
- "g = 9.81; l = 10; k = 0.; # si k = 0.15 donc on a un amortissement\n",
- "p = [g l k]\n",
- "plot()\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,p)\n",
- " sol = solve(prob)\n",
- " plot!(sol,vars=(1,2), xlabel = \"x_1\", ylabel = \"x_2\", legend = false) # lw = linewidth\n",
- "end\n",
- "theta0 = pi-10*eps()\n",
- "x0 = [theta0 0]\n",
- "tf = 70 # problem for tf=50 1/4 of the period!\n",
- "tspan = (0.0,tf)\n",
- "prob = ODEProblem(pendule,x0,tspan,p)\n",
- "sol = solve(prob)\n",
- "plot!(sol,vars=(1,2), xlims = (-2*pi,4*pi), xlabel = \"x_1\", ylabel = \"x_2\", legend = false) # lw = linewidth\n",
- "\n",
- "theta0 = pi+10*eps()\n",
- "x0 = [theta0 0]\n",
- "tf = 70\n",
- "tspan = (0.0,tf)\n",
- "prob = ODEProblem(pendule,x0,tspan,p)\n",
- "sol = solve(prob)\n",
- "plot!(sol,vars=(1,2), xlims = (-2*pi,4*pi), xlabel = \"x_1\", ylabel = \"x_2\", legend = false) # lw = linewidth\n",
- "\n",
- "# circulation case\n",
- "for thetapoint0 in 0:0.2:2 \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,p)\n",
- " sol = solve(prob)\n",
- " plot!(sol,vars=(1,2), xlabel = \"x_1\", ylabel = \"x_2\", legend = false) # lw = linewidth\n",
- "end\n",
- "for thetapoint0 in -2:0.2: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,p)\n",
- " sol = solve(prob)\n",
- " plot!(sol,vars=(1,2), xlabel = \"x_1\", ylabel = \"x_2\", legend = false) # lw = linewidth\n",
- "end\n",
- "plot!([-pi 0 pi 2*pi 3*pi], [0 0 0 0 0], seriestype=:scatter)\n",
- "plot!(xlims = (-pi,3*pi))\n"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Modèle de Van der Pol\n",
- "### Exemple de cycle limite\n",
- "L'équation différentielle considérée est l'équation de Van der Pol\n",
- "$$(IVP)\\left\\{\\begin{array}{l}\n",
- "\\dot{y}_1(t)=y_2(t)\\\\\n",
- "\\dot{y}_2(t)=(1-y_1^2(t))y_2(t)-y_1(t)\\\\\n",
- "y_1(0)=2.00861986087484313650940188\\\\\n",
- "y_2(0)=0\n",
- "\\end{array}\\right.\n",
- "$$\n",
- "$t_f=T=6.6632868593231301896996820305$\n",
- "\n",
- "\n",
- "La solution de ce problème de Cauchy est périodique de période $T$.\n",
- "\n",
- "Résoudre et de visualiser la solution pour les points de départ :\n",
- "- x0 = [2.00861986087484313650940188,0]\n",
- "- [x01 , 0] pour x01 in -2:0.4:0\n",
- "- [0 , x02] pour x02 in 2:1:4\n",
- "\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "function vdp(x,p,t)\n",
- "# Van der Pol model\n",
- "# second member of the IVP\n",
- "# x : state\n",
- "# real(2)\n",
- "# p : parameter vector\n",
- "# t : time, variable not use here\n",
- "# real\n",
- "# Output\n",
- "# xpoint : vector of velocity\n",
- "# same as x\n",
- "# To complete\n",
- "# xpoint = similar(x)\n",
- "# xpoint[1] = x[2]\n",
- "# xpoint[2] = (1-x[1]^2)*x[2] - x[1]\n",
- " xpoint = [x[2] , (1-x[1]^2)*x[2] - x[1]]\n",
- " return xpoint\n",
- "end"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# Main\n",
- "#\n",
- "tf = 6.6632868593231301896996820305\n",
- "tspan = (0.0, tf)\n",
- "x0 = [2.00861986087484313650940188,0]\n",
- "mu = 1\n",
- "p = [mu]\n",
- "t0 = 0.;\n",
- "tspan = (t0,tf)\n",
- "println(\"x0 = $x0, p = $p, tspan = $tspan\")\n",
- "prob = ODEProblem(vdp,x0,tspan,p)\n",
- "sol = solve(prob)\n",
- "plot(sol,vars=(1,2), xlabel = \"x_1\", ylabel = \"x_2\", legend = false)\n",
- "for x01 in -2:0.4:0 #4.5:4.5\n",
- " x0 = [x01,0]\n",
- " prob = ODEProblem(vdp,x0,tspan,p)\n",
- " sol = solve(prob)#,force_dtmin=true)# \n",
- " plot!(sol,vars=(1,2), xlabel = \"x_1\", ylabel = \"x_2\", legend = false) \n",
- "end\n",
- "for x02 in 2:1:4\n",
- " x0 = [0.,x02]\n",
- " prob = ODEProblem(vdp,x0,tspan,p)\n",
- " sol = solve(prob)#solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)\n",
- " plot!(sol,vars=(1,2), xlabel = \"x_1\", ylabel = \"x_2\", legend = false)\n",
- "end\n",
- "plot!([0], [0], seriestype=:scatter) # point visualisation\n",
- "plot!(xlims = (-2.5,5))\n",
- " \n",
- "\n"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "La solution périodique est ici ce qu'on appelle un **cycle limite**"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "# Modèle de Lorentz\n",
- "## Chaos\n",
- "$$(IVP)\\left\\{\\begin{array}{l}\n",
- "\\dot{x}_1(t)=-\\sigma x_1(t)+\\sigma x_2(t)\\\\\n",
- "\\dot{x}_2(t)=-x_1(t)x_3(t)+rx_1(t)-x_2(t)\\\\\n",
- "\\dot{x}_3(t)=x_1(t)x_2(t)-bx_3(t)\\\\\n",
- "x_1(0)=-8\\\\\n",
- "x_2(0)=8\\\\\n",
- "x_3(0)=r-1\n",
- "\\end{array}\\right.\n",
- "$$\n",
- "avec $\\sigma=10, r=28, b=8/3$.\n",
- "\n",
- "Calculer et visualisé la solution de ce problème\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "function lorentz(x,p,t)\n",
- "# Lorentz model\n",
- "# second member of the IVP\n",
- "# x : state\n",
- "# real(3)\n",
- "# p : parameter vector\n",
- "# t : time, variable not use here\n",
- "# real\n",
- "# Output\n",
- "# xpoint : vector of velocity\n",
- "# same as x\n",
- "# To complete\n",
- " sigma = p[1]; rho = p[2]; beta = p[3];\n",
- " xpoint = similar(x)\n",
- " xpoint[1] = sigma*(x[2]-x[1])\n",
- " xpoint[2] = x[1]*(rho-x[3]) - x[2]\n",
- " xpoint[3] = x[1]*x[2] - beta*x[3]\n",
- " return xpoint \n",
- "end"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# test de la fonction lorent!\n",
- "sigma = 10.; rho = 28.; beta = 8/3;\n",
- "p = [sigma rho beta]\n",
- "#x0 = [1.1,0.0,0.0]\n",
- "x0 = [-8, 8, rho-1]\n",
- "xpoint = x0\n",
- "xpoint = lorentz(x0,p,0)\n",
- "println(\"xpoint = $xpoint\")"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# Main\n",
- "#\n",
- "x0 = [-8, 8, rho-1]\n",
- "tspan = (0.0,100.0)\n",
- "prob = ODEProblem(lorentz,x0,tspan,p)\n",
- "sol = solve(prob)\n",
- "plot(sol,vars=(1,2,3), xlabel = \"x_1\", ylabel = \"x_2\", zlabel = \"x_3\", legend = false)\n",
- "plot!([x0[1]], [x0[2]], [x0[3]], seriestype=:scatter) # point visualisation\n",
- "#annotate!([x0[1],x0[2],x0[3],text(\". Point de départ\", 10,:left)])"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "p1 = plot(sol,vars=(0,1),ylabel = \"x_1(t)\", legend = false)\n",
- "p2 = plot(sol,vars=(0,2),ylabel = \"x_2(t)\", legend = false)\n",
- "p3 = plot(sol,vars=(0,3),xlabel = \"t\", ylabel = \"x_3(t)\", legend = false)\n",
- "plot(p1,p2,p3,layout=(3,1))#,legend=false)"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Attention à bien utiliser les codes\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",
- "dont la solution 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*}"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "function exemple1(x,p,t)\n",
- "# second member of the IVP\n",
- "# Input\n",
- "# x : state\n",
- "# real(2)\n",
- "# p : parameter vector\n",
- "# t : time\n",
- "# real\n",
- "# Output\n",
- "# xpoint : vector of velocity\n",
- "# same as x\n",
- " xpoint = similar(x) # xpoint est un objet de même type que x\n",
- " xpoint[1] = x[1] + x[2] + sin(t)\n",
- " xpoint[2] = -x[1] + 3*x[2] \n",
- "# xpoint = [x[1]+x[2]+sin(t), -x[1] + 3*x[2]]\n",
- " return xpoint \n",
- "end\n",
- "p = []\n",
- "t0 = 0.; tf = 8\n",
- "tspan = (t0,tf)\n",
- "x0 = [-9/25 , -4/25]\n",
- "prob = ODEProblem(exemple1,x0,tspan,p) # défini le problème en Julia\n",
- "sol = solve(prob, DP5()) # réalise l'intégration numérique\n",
- " # avec ode45 de matlab\n",
- "p1 = plot(sol, label = [\"x_1(t)\" \"x_2(t)\"], title = \"ode45 de Matlab\") #, lw = 0.5) # lw = linewidth\n",
- "sol = solve(prob, DP5(), reltol = 1.e-6, abstol = 1.e-9)\n",
- "p2 = plot(sol, label = [\"x_1(t)\" \"x_2(t)\"], title = \"reltol=1.e-6, abstol=1.e-9\") #, lw = 0.5)\n",
- "sol = solve(prob, DP5(), reltol = 1.e-10, abstol = 1.e-15)\n",
- "p3 = plot(sol, label = [\"x_1(t)\" \"x_2(t)\"], title = \"reltol=1.e-10, abstol=1.e-16\") #, lw = 0.5)\n",
- "T = t0:(tf-t0)/100:tf\n",
- "sinT = map(sin,T) # opération vectorielle\n",
- "cosT = map(cos,T)\n",
- "p4 = plot(T,[(-1/25)*(13*sinT+9*cosT) (-1/25)*(3*sinT+4*cosT)], label = [\"x_1(t)\" \"x_2(t)\"], xlabel = \"t\", title = \"Solution exacte\")\n",
- "plot(p1,p2,p3,p4)\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": []
- }
- ],
- "metadata": {
- "kernelspec": {
- "display_name": "Julia 1.6.0",
- "language": "julia",
- "name": "julia-1.6"
- },
- "language_info": {
- "file_extension": ".jl",
- "mimetype": "application/julia",
- "name": "julia",
- "version": "1.6.0"
- }
- },
- "nbformat": 4,
- "nbformat_minor": 4
-}