INFO-F-305 - TP 7

Dessin qualitiatif - Portrait de phase IV

Jacopo De Stefani - jdestefa@ulb.ac.be

Github: https://github.com/jdestefani/ulb-infof305

Slides: https://jdestefani.github.io/ulb-infof305

Système du second ordre non-linéaire

$ \begin{cases} \dot{x_{1}} &= f_1(x_{1},x_{2}) \\ \dot{x_{2}} &= f_2(x_{1},x_{2}) \\ \end{cases} $

Matrice de Jacobi

$J=\begin{bmatrix} \frac{\partial f_1}{\partial x_1} & \frac{\partial f_1}{\partial x_2} \\ \frac{\partial f_2}{\partial x_1} & \frac{\partial f_2}{\partial x_2} \\ \end{bmatrix}$

N.B. La matrice jacobienne, évaluée en un point stationnaire $(\bar{x_1}, \bar{x_2})$ donne une approximation linéaire du comportement du système au voisinage de ce point.

Stabilité d'un système du second ordre non-linéaire

Stabilité asymptotique - Th 6.1

Si le système linéarisé obtenu par linéarisation du système autour de $\mathbf{\bar{x}}$ est stable asymptotiquement (c-à-d. toutes les valeurs propres de A sont négatives), alors l’état d’équilibre $\mathbf{\bar{x}}$ est stable asymptotiquement.

Instabilité asymptotique - Th 6.2

Si la matrice A du système linéarisé a une ou plusieurs valeurs propres avec partie réelle positive, l’état d’équilibre $\mathbf{\bar{x}}$ du système non-lineaire est instable.

N.B. Quand la matrice a toutes ses valeurs propres avec partie réelle négative mis à part quelques valeurs propres avec partie réelle nulle, nous ne pouvons rien déduire sur la stabilité de $\mathbf{\bar{x}}$ à partir de l’analyse du système linéarisé.

Rappels Théoriques - Système du second ordre

$\dot{\mathbf{X}}(t) = \mathbf{A}\mathbf{X}(t)$

$\begin{bmatrix}\dot{x_{1}}\\ \dot{x_{2}}\end{bmatrix}=\begin{bmatrix}a_{11} & a_{12}\\ a_{21} & a_{22}\end{bmatrix}\begin{bmatrix}x_{1}\\ x_{2}\end{bmatrix}$

Equation caractéristique

$\det(A-\lambda I) = \lambda^{2}-\tau\lambda+\Delta=0$

$\lambda_{1,2} = \frac{\tau \pm \sqrt{\tau^2 - 4 \Delta}}{2}$

où:

• $\tau=a_{11}+a_{22}=tr(A)$
• $\Delta=a_{11}a_{22}-a_{12}a_{21}=\det(A)$

Stabilité - Racines réelles et distinctes - $\lambda_1,\lambda_2$

Noeud stable

Si valeurs propres de $\mathbf{A}$ sont réelles et négatives.

Noeud instable

Si valeurs propres de $\mathbf{A}$ sont réelles et positives.

Selle

Si valeurs propres de $\mathbf{A}$ sont réelles et ayant signe discorde.

Système non simple

Considérons le cas où $\det(A) = \lambda_1\lambda_2 = 0$ et la valeur propre $\lambda_2 \neq 0$.

Tous les états qui appartiennent à la droite $a_{11}x_1 + a_{12}x_2 = 0$ sont des états d’équilibre. Aussi, toutes les trajectoires sont des droites parallèles à la droite $v_2$.

Deux configurations sont possibles:

• $0 = \lambda_1 > \lambda_2$ : ceci implique que tous les états d’équilibre sont stables.
• $ 0 = \lambda_1 < \lambda_2$: ceci implique que tous les états d’équilibre sont instables.

Stabilité - Racines réelles et multiples - $\lambda=\lambda_1=\lambda_2$

Systéme simple - $\lambda_1 = \lambda_2 = \lambda \neq 0$

• A diagonalisable - : Tous les vecteurs sont vecteurs propres. Chaque droite qui passe par l’origine est une trajectoire. Si $\lambda < 0(> 0)$ nous avons la (in)stabilité asymptotique. L’origine est dénommée noeud singulier.
• A non diagonalisable: il existe un seul vecteur propre et donc seul une droite qui contient une trajectoire. L’origine est dénommée noeud dégénéré.

Systéme non simple - $\lambda_1 = \lambda_2 = \lambda = 0$

• Il y a une infinité de points d’équilibres instables (multiplicité plus grande que 1) et toutes les trajectoires se trouvent sur des droites parallèles.

Stabilité - Racines complexes et conjuguées - $\lambda_{1,2} = \alpha \pm i\beta$

Centre

• $\alpha = 0$: Les trajectoires sont des ellipses fermées avec période $T = \frac{2\pi}{b}$. L’origine est dit un centre.

Foyer stable

• $\alpha < 0$: Le système est asymptotiquement stable et les trajectoires convergent vers l’origine en suivant des spirales.

Foyer instable

• $\alpha > 0$: Le système est instable et les trajectoires s’éloignent de l’origine en suivant des spirales.

Demarche pour systèmes non linéaires

1. Calcul isoclines
2. Calcul point stationaires (d'équilibre)
3. Calcul matrice de Jacobi
4. Etude des systèmes linéarisés aux alentours des point stationaires
5. Dessin qualitatif

Exercice 1

$ \begin{cases} \dot{x_1} &= x_1 \\ \dot{x_2} &= x_1^2 + x_2^2 -1 \\ \end{cases} $

1. Calcul des isoclines et partition du plan en régions.
2. Calcul des points stationnaires (par annulation des dérivées).
3. Calcul de la matrice jacobienne.
4. Pour chaque point stationnaire : a. Appliquer la matrice stationnaire au point. b. Etudier le système linéaire qui en découle, et qui donnera une approximation du comportement du système au voisinage de ce point stationnaire.
5. Tracer des trajectoires qualitativement.

# 3. Calcul de la matrice jacobienne.
pkg load symbolic
syms x1 x2;
sys = [x1 ; x1^2 + x2^2 - 1]
jacobian(sys)

Symbolic pkg v2.7.1: Python communication link active, SymPy v1.4.
sys = (sym 2x1 matrix)

  [     x1      ]
  [             ]
  [  2     2    ]
  [x1  + x2  - 1]

ans = (sym 2x2 matrix)

  [ 1     0  ]
  [          ]
  [2*x1  2*x2]

# 4. Linearisation dans les points stationaires
A = [1 0; 0 2]; # Linearisation pour [0,1]
[V,L] = eig(A)

V =

   1   0
   0   1

L =

Diagonal Matrix

   1   0
   0   2

portrait_phase(A,[0 1])

A = [1 0; 0 -2] # Linearisation pour [0,-1]
[V,L] = eig(A)

A =

   1   0
   0  -2

V =

   0   1
   1   0

L =

Diagonal Matrix

  -2   0
   0   1

portrait_phase(A,[0 -1])

#5. Tracer des trajectoires qualitativement.
# Define grid for plotting
x1range=-3:.25:3;
x2range=-3:.25:3;
[x1,x2] = meshgrid(x1range, x2range);

# Define the system to plot (based on matrix A)
x1p = x1;
x2p = x1^2+x2^2-1;

# Plot circle
r = 1;
t = linspace(0,2*pi,100)'; 
circsx = r.*cos(t); 
circsy = r.*sin(t); 

#Normalize values for plotting
arrow=sqrt(x1p.^2+x2p.^2);

# Vector field plot
hold on;
quiver(x1,x2,x1p./arrow,x2p./arrow,0.7);

plot(circsx,circsy,"linewidth",6,"linestyle","--");
plot([0,0],[-2,2],"linewidth",6,"linestyle","--")

# Plot stationary points
plot(0,1, "*" ,"linewidth",4);
plot(0,-1,"*","linewidth",4);

grid on;
axis tight;

Exercice 2

$ \begin{cases} \dot{x_1} &= x_2 \\ \dot{x_2} &= x_1(1-x_1^2) + x_2 \\ \end{cases} $

1. Calcul des isoclines et partition du plan en régions.
2. Calcul des points stationnaires (par annulation des dérivées).
3. Calcul de la matrice jacobienne.
4. Pour chaque point stationnaire : a. Appliquer la matrice stationnaire au point. b. Etudier le système linéaire qui en découle, et qui donnera une approximation du comportement du système au voisinage de ce point stationnaire.
5. Tracer des trajectoires qualitativement.

# 3. Calcul de la matrice jacobienne.
pkg load symbolic
syms x1 x2;
sys = [x2 ; x1*(1-x1^2) + x2]
jacobian(sys)

sys = (sym 2x1 matrix)

  [       x2        ]
  [                 ]
  [   /      2\     ]
  [x1*\1 - x1 / + x2]

ans = (sym 2x2 matrix)

  [    0      1]
  [            ]
  [        2   ]
  [1 - 3*x1   1]

# 4. Linearisation dans les points stationaires
A = [0 1; 1 1] # Linearisation pour [0,0]
[V,L] = eig(A)

A =

   0   1
   1   1

V =

  -0.85065   0.52573
   0.52573   0.85065

L =

Diagonal Matrix

  -0.61803         0
         0   1.61803

portrait_phase(A,[0 0])

# 4. Linearisation dans les points stationaires
A = [0 1; -2 1] # Linearisation pour [0,1] et [0,-1]
[V,L] = eig(A)

A =

   0   1
  -2   1

V =

   0.20412 - 0.54006i   0.20412 + 0.54006i
   0.81650 + 0.00000i   0.81650 - 0.00000i

L =

Diagonal Matrix

   0.5000 + 1.3229i                  0
                  0   0.5000 - 1.3229i

portrait_phase(A,[0 1])

#5. Tracer des trajectoires qualitativement.
# Define grid for plotting
x1range=-3:.25:3;
x2range=-3:.25:3;
[x1,x2] = meshgrid(x1range, x2range);

# Define the system to plot (based on matrix A)
x1p = x2;
x2p = -(x1.^3)+x1+x2;

#Normalize values for plotting
arrow=sqrt(x1p.^2+x2p.^2);

# Vector field plot
hold on;
quiver(x1,x2,x1p./arrow,x2p./arrow,0.7);

# Plot isocline
x1_red = -1.65:.1:1.65;
plot(x1_red,x1_red.^3-x1_red,"linewidth",6,"linestyle","--")
plot([-3,3],[0,0],"linewidth",6,"linestyle","--")


# Plot stationary points
plot(0,0, "*" ,"linewidth",4);
plot(1,0, "*" ,"linewidth",4);
plot(-1,0,"*","linewidth",4);

grid on;
axis tight;

Exercice 3 - Lotka-Volterra

$ \begin{cases} \dot{x_1} &= 0.2x_1 - 0.08x_1x_2 \\ \dot{x_2} &= 0.1x_1x_2 - 0.2x_2 \\ \end{cases} $

1. Calcul des isoclines et partition du plan en régions.
2. Calcul des points stationnaires (par annulation des dérivées).
3. Calcul de la matrice jacobienne.
4. Pour chaque point stationnaire : a. Appliquer la matrice stationnaire au point. b. Etudier le système linéaire qui en découle, et qui donnera une approximation du comportement du système au voisinage de ce point stationnaire.
5. Tracer des trajectoires qualitativement.

warning("off","all")
# 3. Calcul de la matrice jacobienne.
pkg load symbolic
syms x1 x2;
sys = [0.2*x1 - 0.08*x1*x2 ; 0.1*x1*x2 - 0.2*x2]
jacobian(sys)

sys = (sym 2x1 matrix)

  [  2*x1*x2   x1]
  [- ------- + --]
  [     25     5 ]
  [              ]
  [  x1*x2   x2  ]
  [  ----- - --  ]
  [    10    5   ]

ans = (sym 2x2 matrix)

  [1   2*x2  -2*x1 ]
  [- - ----  ------]
  [5    25     25  ]
  [                ]
  [   x2     x1   1]
  [   --     -- - -]
  [   10     10   5]

# 4. Linearisation dans les points stationaires
A = [0.2 0; 0 0.2] # Linearisation pour [0,0]
[V,L] = eig(A)

A =

   0.20000   0.00000
   0.00000   0.20000

V =

   1   0
   0   1

L =

Diagonal Matrix

   0.20000         0
         0   0.20000

portrait_phase(A,[0 0])

# 4. Linearisation dans les points stationaires
A = [0 -0.16; 0.25 0] # Linearisation pour [2,2.5]
[V,L] = eig(A)

A =

   0.00000  -0.16000
   0.25000   0.00000

V =

   0.00000 + 0.62470i   0.00000 - 0.62470i
   0.78087 + 0.00000i   0.78087 - 0.00000i

L =

Diagonal Matrix

   0.00000 + 0.20000i                    0
                    0   0.00000 - 0.20000i

portrait_phase(A,[2 2.5])

#5. Tracer des trajectoires qualitativement.
# Define grid for plotting
x1range=-1:.2:5;
x2range=-1:.2:5;
[x1,x2] = meshgrid(x1range, x2range);

# Define the system to plot (based on matrix A)
x1p = 0.2*x1 - 0.08*(x1.*x2);
x2p = 0.1*(x1.*x2) - 0.2*x2;

#Normalize values for plotting
arrow=sqrt(x1p.^2+x2p.^2);

# Vector field plot
hold on;
quiver(x1,x2,x1p./arrow,x2p./arrow,0.5);

# Plot isoclines
plot([-1,5],[2.5,2.5],"linewidth",6,"linestyle","--")
plot([2,2],[-1,5],"linewidth",6,"linestyle","--")

# Plot stationary points
plot(0,0, "*" ,"linewidth",4);
plot(2,2.5, "*" ,"linewidth",4);

grid on;
axis tight;