{
"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": 31,
"metadata": {},
"outputs": [],
"source": [
"# import Pkg; Pkg.add(\"DifferentialEquations\")\n",
"# Pkg.add(\"Plots\")\n",
"# import Pkg; Pkg.add(\"ODEInterfaceDiffEq\")\n",
"using DifferentialEquations\n",
"using Plots"
]
},
{
"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": 32,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"pendule (generic function with 1 method)"
]
},
"execution_count": 32,
"metadata": {},
"output_type": "execute_result"
}
],
"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": 33,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 33,
"metadata": {},
"output_type": "execute_result"
}
],
"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": 34,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"pendule1! (generic function with 1 method)"
]
},
"execution_count": 34,
"metadata": {},
"output_type": "execute_result"
}
],
"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": 35,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 35,
"metadata": {},
"output_type": "execute_result"
}
],
"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": 36,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"pendule2! (generic function with 1 method)"
]
},
"execution_count": 36,
"metadata": {},
"output_type": "execute_result"
}
],
"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": 37,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 37,
"metadata": {},
"output_type": "execute_result"
}
],
"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": 38,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 38,
"metadata": {},
"output_type": "execute_result"
}
],
"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": 39,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"vdp (generic function with 1 method)"
]
},
"execution_count": 39,
"metadata": {},
"output_type": "execute_result"
}
],
"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": 40,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"x0 = [2.0086198608748433, 0.0], p = [1], tspan = (0.0, 6.66328685932313)\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 40,
"metadata": {},
"output_type": "execute_result"
}
],
"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": 41,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"lorentz (generic function with 1 method)"
]
},
"execution_count": 41,
"metadata": {},
"output_type": "execute_result"
}
],
"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": 42,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"xpoint = [160.0, -16.0, -136.0]\n"
]
}
],
"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": 43,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 43,
"metadata": {},
"output_type": "execute_result"
}
],
"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": 44,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 44,
"metadata": {},
"output_type": "execute_result"
}
],
"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": 45,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 45,
"metadata": {},
"output_type": "execute_result"
}
],
"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.4.1",
"language": "julia",
"name": "julia-1.4"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.4.1"
}
},
"nbformat": 4,
"nbformat_minor": 2
}