{
"cells": [
{
"cell_type": "markdown",
"id": "b8b9cbb6-7d46-4378-9df9-9a6287861d89",
"metadata": {},
"source": [
"# MT94/P26/TD - Systèmes dynamiques"
]
},
{
"cell_type": "markdown",
"id": "abd43863-06ac-4f17-923a-697695ff4fdd",
"metadata": {},
"source": [
"## Equation de Duffing\n",
"### Système autonome\n",
"On considère l'équation différentielle\n",
"$$\n",
"x'' + k x' -x + x^3 = 0,\n",
"$$\n",
"avec $k=0.15$.\n",
"1. Mettre cette équation différentielle sous forme d'un système d'équations différentielles du premier ordre en temps"
]
},
{
"cell_type": "markdown",
"id": "dfb42dae-31e9-4bcd-a8ee-fddc56136bc5",
"metadata": {},
"source": [
"$$\n",
"X' = f(t,X),\n",
"$$\n",
"puis déterminer les points d'équilibre de ce système. Ces points d'équilibre sont-ils stables ou instables ? "
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "9f19ec41-1c36-415a-9db5-2c862cd1261f",
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"id": "2a8ad920-4a75-447f-ae24-31d4b2e23594",
"metadata": {},
"source": [
"Utiliser la fonction `cvode` de Scilab pour calculer la solution de l'équation différentielle avec des conditions initiales $X(0)$ choisies de manière à illustrer les résultats établis à la question précédente. On pourra par exemple représenter le champ de vecteurs $f(X)$ et les solutions obtenues dans le plan de phase $(O,X_1,X_2)$."
]
},
{
"cell_type": "code",
"execution_count": 15,
"id": "a48c883f-e7dd-48a3-814b-67018aae239a",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"function dXdt=f_duffing(t,X)\n",
" k=0.15;\n",
" dXdt = [X(2); -k*X(2)+X(1)-X(1)^3];\n",
"endfunction\n",
"t = linspace(0,5,1000);\n",
"for X10=[-0.5:0.5:0.5]\n",
" for X20=[-0.5:0.5:0.5]\n",
" [_,X] = cvode(f_duffing,t,[X10;X20]);\n",
" plot(X(1,:),X(2,:),X10,X20,\"or\")\n",
" end\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "2f39fe00-c169-4238-997e-f40b5e72c505",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\u001b[4l\n",
" fchamp: direction field of a 2D first order ODE\n",
"\n",
" Syntax:\n",
" \n",
" fchamp(f, t, xr, yr, [arfact, rect, strf])\n",
" fchamp(f, t, xr, yr, )\n",
" h = fchamp(...)\n",
"\n",
" See also: champ, champ_properties, xarrows, polyline_properties\n",
"\n",
" Online reference page : https://help.scilab.org/fchamp.html \n",
"\n"
]
}
],
"source": [
"help fchamp"
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "e8aeb9f2-2fff-4ef0-a95e-d3dad785881d",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fchamp(f_duffing,0, 1+linspace(-0.1,0.1,10), linspace(-0.1,0.1,10))\n",
"plot(1,0,'xr')\n",
"[_,X] = cvode(f_duffing,t,[1.05;0]);\n",
"plot(X(1,:),X(2,:),1.05,0,\"or\")\n",
"gca().data_bounds=[0.9 1.1 -0.1 0.1];\n",
"\n",
"scf;\n",
"\n",
"fchamp(f_duffing,0, linspace(-0.1,0.1,10), linspace(-0.1,0.1,10))\n",
"plot(0,0,'xr')\n",
"[_,X] = cvode(f_duffing,t,[0.05;0]);\n",
"plot(X(1,:),X(2,:),0.05,0,\"or\")\n",
"gca().data_bounds=[-0.1 0.1 -0.1 0.1];\n"
]
},
{
"cell_type": "markdown",
"id": "c972b65b-8543-4395-b631-d54285d1ff58",
"metadata": {},
"source": [
"## Système non autonome\n",
"On considère maintenant l'équation différentielle\n",
"$$\n",
"x'' + k x' -x + x^3 = b\\cos t,\n",
"$$\n",
"avec $b=0.3$ et pour conditions initiales $x(0)=\\dot x(0)=0$.\n",
"1. Calculez avec `cvode` la solution de l'équation différentielle pour $t\\in[0,40\\pi]$.\n",
" "
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "3b57c209-7852-40a3-a02e-74fc749925c1",
"metadata": {},
"outputs": [],
"source": [
"function dXdt=f_duffing(t,X)\n",
" k=0.15;\n",
" b=0.3;\n",
" dXdt = [X(2); -k*X(2)+X(1)-X(1)^3+b*cos(t)];\n",
"endfunction\n",
"t = linspace(0,40*%pi,1000);\n",
"[_,Xa] = cvode(f_duffing,t,[0;0]);"
]
},
{
"cell_type": "markdown",
"id": "adb6fb5f-9ec9-43e5-a756-2feeff8143cd",
"metadata": {},
"source": [
"2. Calculez la solution en prenant les conditions initiales $x(0)=10^{-6},\\dot x(0)=0$. Superposer les deux solutions obtenues et commentez le phénomène observé."
]
},
{
"cell_type": "code",
"execution_count": 16,
"id": "6d47302e-6b5d-46b3-af91-df3f3c40975d",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"[_,Xb] = cvode(f_duffing,t,[1e-6;0]);\n",
"plot(t,Xa(1,:),t,Xb(1,:))"
]
},
{
"cell_type": "markdown",
"id": "b6cdd6aa-05e0-49d8-8257-8b0daed1a510",
"metadata": {},
"source": [
"On constate un comportement chaotique (dépendance à des perturbations infinitésimales des conditions initiales)"
]
},
{
"cell_type": "markdown",
"id": "4122297d-05ac-493b-9623-2e5daad97677",
"metadata": {},
"source": [
"## Equation de Van der Pol\n",
"On considère l'équation différentielle\n",
"$$\n",
"x'' - \\mu (1-x^2)x' + x = 0.\n",
"$$\n",
"1. Calculer les points d'équilibre et déterminer leur stabilité."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "9ce0e954-096a-4dfb-8939-f320e5bfdb4d",
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"id": "29981870-125e-4352-89d0-6afda01c4269",
"metadata": {},
"source": [
"2. Illustrer le résultat obtenu à la question précédente en approchant la solution de l'équation différentielle pour un choix adéquat de condition initiale. On pourra prendre des valeurs de $\\mu$ comprises entre 0.01 et 5. Commenter les résultats obtenus."
]
},
{
"cell_type": "code",
"execution_count": 36,
"id": "a1bfc426-728f-432d-9a16-104e890ab7e3",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"function dXdt=f_vdp(t,X)\n",
" mu=5;\n",
" dXdt = [X(2); mu*(1-X(1)^2)*X(2)-X(1)];\n",
"endfunction\n",
"t = linspace(0,20,1000);\n",
"for X10=[-0.5:0.5:0.5]\n",
" for X20=[-1:1:1]\n",
" [_,X] = cvode(f_vdp,t,[X10;X20]);\n",
" plot(X(1,:),X(2,:),X10,X20,'or')\n",
" end\n",
"end\n",
"h = plot(X(1,t>5),X(2,t>5),'r','linewidth',2);\n",
"legend(h,\"Cycle limite\",-1)\n"
]
},
{
"cell_type": "markdown",
"id": "ccb6b668-2a01-4f1f-a33e-610a23944cad",
"metadata": {},
"source": [
"## Système proie-prédateur (Lokta-Volterra)\n",
"\n",
"$$\\begin{align}\n",
"x'&=x(\\alpha-\\beta y)\\\\\n",
"y'&=-y(\\gamma-\\delta x)\n",
"\\end{align}\n",
"$$\n",
"Etudiez les points fixes et simulez le système pour\n",
"$\\alpha=2$, $\\beta=10^{-3}$, $\\gamma=10$, $\\delta=2.10^{-3}$."
]
},
{
"cell_type": "code",
"execution_count": 59,
"id": "39f60d5c-be63-4f64-8423-e82c97ceb307",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"function dXdt=f_lv(t,X)\n",
" alpha=2;\n",
" beta=1e-3;\n",
" gamma=10;\n",
" delta=2e-3;\n",
" x=X(1);\n",
" y=X(2);\n",
" dXdt = [x*(alpha-beta*y); -y*(gamma-delta*x)];\n",
"endfunction\n",
"\n",
"[t,X] = cvode(f_lv,[0 10],[100;100]);\n",
"plot(t,X)\n",
"legend(\"x (prédateurs)\",\"y (proies)\",-1)\n",
"scf;\n",
"plot(X(1,:),X(2,:))\n",
"title(\"Portrait de phase\")"
]
},
{
"cell_type": "markdown",
"id": "0fb02a9b-e519-49de-8664-987ab9e2f2e4",
"metadata": {},
"source": [
"## Système de Lorenz\n",
"\n",
"Etudiez les éventuels points fixes du système de Lorenz\n",
"\\begin{align*}\n",
"x'&=\\sigma(y-x),\\\\\n",
"y'&=\\rho x-y-xz,\\\\\n",
"z'&=xy-\\beta z,\n",
"\\end{align*}\n",
"où $\\sigma$, $\\rho$ et $\\beta$ sont des nombres positifs. Représentez les trajectoires du système de Lorenz pour $\\sigma=10$, $\\rho=28$, $\\beta=\\frac{8}{3}$. "
]
},
{
"cell_type": "code",
"execution_count": 75,
"id": "8659229d-e2c6-4a1f-a522-2b415f0fc813",
"metadata": {},
"outputs": [],
"source": [
"function dXdt=f_lorenz(t,X)\n",
" beta=8/3;\n",
" rho=28;\n",
" sigma=10;\n",
" x=X(1);\n",
" y=X(2);\n",
" z=X(3);\n",
" dXdt = [sigma*(y-x); rho*x-y-x*z; x*y-beta*z]\n",
"endfunction\n",
"\n",
"[t,X] = cvode(f_lorenz,[0 100],[1;0;0], refine=2);"
]
},
{
"cell_type": "code",
"execution_count": 77,
"id": "12022410-29a4-44ac-a897-9e19dba8a681",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"h=param3d(X(1,:),X(2,:),X(3,:));\n",
"h.foreground=color(\"blue\");\n",
"isoview on\n",
"gca().rotation_angles=[65 -57];\n",
"gcf().figure_size(2)=800;"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d1e8e254-4471-451f-bbff-91c995aca968",
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Scilab",
"language": "scilab",
"name": "scilab"
},
"language_info": {
"file_extension": ".sci",
"help_links": [
{
"text": "MetaKernel Magics",
"url": "https://metakernel.readthedocs.io/en/latest/source/README.html"
}
],
"mimetype": "text/x-scilab",
"name": "scilab",
"version": "0.10.2"
}
},
"nbformat": 4,
"nbformat_minor": 5
}