INFO-F-305 - TP 4

Dessin qualitiatif - Portrait de phase III

Jacopo De Stefani - jdestefa@ulb.ac.be

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

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

Project

• Analyse des systèmes ODE (ordre 1-2) avec Octave
• Simulation des systèmes et comparaison avec des données réelles
• 2 x Séances TP en salle machine + Travail personnel
• Durée: Environ 4 semaines
• A remettre: Jupyter Notebook/Rapport LaTeX sur UV

A = [-5 -3; 1 -1]
trace_det_condition(A)
portrait_phase(A)

A =

  -5  -3
   1  -1

Tr(A) = -6
det(A) = 8
Tr(A) > 4*det(A)? = 0
sys_str = [-5*x(1)+-3*x(2);1*x(1)+-1*x(2)]
sys_ode = f(t, x) = [-5*x(1)+-3*x(2);1*x(1)+-1*x(2)]

A = [2 1; 5 -2]
trace_det_condition(A)
portrait_phase(A)

A =

   2   1
   5  -2

Tr(A) = 0
det(A) = -9
Tr(A) > 4*det(A)? = 1
sys_str = [2*x(1)+1*x(2);5*x(1)+-2*x(2)]
sys_ode = f(t, x) = [2*x(1)+1*x(2);5*x(1)+-2*x(2)]

A = [1/2 2; 1 1/2]
trace_det_condition(A)
portrait_phase(A)

A =

   0.50000   2.00000
   1.00000   0.50000

Tr(A) = 1
det(A) = -1.75
Tr(A) > 4*det(A)? = 1
sys_str = [0.5*x(1)+2*x(2);1*x(1)+0.5*x(2)]
sys_ode = f(t, x) = [0.5*x(1)+2*x(2);1*x(1)+0.5*x(2)]

A = [3/4 11; -1/8 1/2]
trace_det_condition(A)
portrait_phase(A)

A =

    0.75000   11.00000
   -0.12500    0.50000

Tr(A) = 1.25
det(A) = 1.75
Tr(A) > 4*det(A)? = 0
sys_str = [0.75*x(1)+11*x(2);-0.125*x(1)+0.5*x(2)]
sys_ode = f(t, x) = [0.75*x(1)+11*x(2);-0.125*x(1)+0.5*x(2)]

A = [0 1; -7 3]
trace_det_condition(A)
portrait_phase(A)

A =

   0   1
  -7   3

Tr(A) = 3
det(A) = 7
Tr(A) > 4*det(A)? = 0
sys_str = [0*x(1)+1*x(2);-7*x(1)+3*x(2)]
sys_ode = f(t, x) = [0*x(1)+1*x(2);-7*x(1)+3*x(2)]

Système linéaire, regulier, à dimensions finies - Th 5.1

$\begin{cases} \dot{\mathbf{X}}(t) & = \mathbf{A}(t)\mathbf{X}(t) + \mathbf{B}(t)\mathbf{U}(t) \\ \mathbf{Y}(t) & = \mathbf{C}(t)\mathbf{X}(t) \\ \end{cases}$

où $\mathbf{A}(\cdot),\mathbf{B}(\cdot),\mathbf{C}(\cdot)$ sont des matrices continues en $T$.

Si on considère le mouvement libre (c-à-d $\mathbf{U}(t) = 0$ $\forall t$), et on se concentre sur l'état, avec $\mathbf{A}(t)$ constante:

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

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é et valeurs propres

Stabilité asymptotique - Th 5.3

Un système $\dot{\mathbf{X}}(t) = \mathbf{A}\mathbf{X}(t)$ est asymptotiquement stable si la partie réelle de tous les valeurs propres de A est negative.

Stabilité simple - Th 5.4

Un système $\dot{\mathbf{X}}(t) = \mathbf{A}\mathbf{X}(t)$ est (simplement) stable si la partie réelle de tous les valeurs propres de A est negative ou nulle et celles ayant partie réelle nulle ont une multiplicité égale à 1.

Instabilité - Th 5.5

Un système $\dot{\mathbf{X}}(t) = \mathbf{A}\mathbf{X}(t)$ est instable s'il existe soit une valeur propre de A avec partie réelle positive soit une valeur propre ayant partie réelle nulle et une multiplicité supérieure à 1.

Points d'équilibre

Definition

Un point d’équilibre est un point où la trajectoire coincide avec la condition initiale.

• Si $\det(\mathbf{A}) \neq 0$ alors le seul état d’équilibre est (0, 0) et le système corréspondant est dit simple
• Si $\det(\mathbf{A}) = 0$ alors il existe plusieurs état d'équilibre et le système corréspondant est dit non simple

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

Solution générale - Théoreme 5.11

$\mathbf{X}(t) = C_1 \mathbf{X}^{(1)}(t) + C_2 \mathbf{X}^{(2)}(t) $

Solution particulières

$\mathbf{X}^{(1)}(t) = e^{\lambda_1 t} \mathbf{v_1} $

$\mathbf{X}^{(2)}(t) = e^{\lambda_2 t} \mathbf{v_2} $

Exemple

$\mathbf{\dot{x}}=\begin{bmatrix}1 & 3\\ 3 & 1\end{bmatrix}\mathbf{x}$

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 positive.

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.

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

Solution générale - Théoreme 5.11

$\mathbf{X}(t) = C_1 \mathbf{X}^{(1)}(t) + C_2 \mathbf{X}^{(2)}(t) $

Solution particulières

$\mathbf{X}^{(1)}(t) = e^{\lambda t} \mathbf{v_1} $

$\mathbf{X}^{(2)}(t) = e^{\lambda t} (\mathbf{v_2} + t\mathbf{v_1}) $

Exemple

$\mathbf{\dot{x}}=\begin{bmatrix}3 & -4\\ 1 & -1\end{bmatrix}\mathbf{x}$

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.

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

Solution générale - Théoreme 5.11

$\mathbf{X}(t) = C_1 \mathbf{X}^{(1)}(t) + C_2 \mathbf{X}^{(2)}(t) $

Solution particulières - $\mathbf{v_{1,2}} = \mathbf{u} + i \mathbf{v}$

$\mathbf{X}^{(1)}(t) = e^{\alpha t} (\cos(\beta t)\mathbf{u} - \sin(\beta t)\mathbf{v} ) $

$\mathbf{X}^{(2)}(t) = e^{\alpha t} (\sin(\beta t)\mathbf{u} + \cos(\beta t)\mathbf{v}) $

Exemple

$\mathbf{\dot{x}}=\begin{bmatrix}0 & 4\\ -1 & 0\end{bmatrix}\mathbf{x}$

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.

Plan Trace-Determinant

$tr(A)= \sum_i \lambda_i$

$det(A)= \prod_i \lambda_i$

$det(A) = \frac{1}{4} tr(A)^2$

x = -10:0.1:10;
plot (x, 1/4*x.^2);
hold on;
A_vec = [[-2 1; 1 -2];[2 1; 2 3];[5 9; 6 2];[0 1; 0 1];[0 -1; 1 0];[1/3 -2; 3 -1]]; #cbind like
i = 0;
while i<(rows(A_vec)/2)# Column-wise iteration
    A = A_vec(2*i+1:2*(i+1),1:2); # Get the proper submatrix
    plot(trace(A),det(A),sprintf(";[%d %d];o",trace(A),det(A)))
    i++;
endwhile
set(gca, "xaxislocation", "origin");
set(gca, "yaxislocation", "origin");
h = legend; legend(h,"location","southwest");
xlabel("tr(A)");
ylabel("det(A)");
box off;
hold off;

Exercice 1

$ \mathbf{A} = \begin{bmatrix} 2 & 1 \\ 1 & 2 \end{bmatrix} $

1. Calculer les valeurs propres et vecteurs propres et en déduire le type du système
2. Calculer les isoclines et dessiner des vecteurs vitesse sur ces isoclines
3. Associer les signes des portions du plan délimitées par les isoclines.
4. Tracer des trajectoires qualitativement.

# 1. Calcul valeur propres et vecteurs propres
A = [2 1; 1 2]
[V,L] = eig(A)

A =

   2   1
   1   2

V =

  -0.70711   0.70711
   0.70711   0.70711

L =

Diagonal Matrix

   1   0
   0   3

printf("Tr(A) = %d\n",trace(A));
printf("det(A) = %d\n",det(A));
printf("Tr(A) > 4*det(A)? = %d\n",(trace(A) > 4*det(A)));

Tr(A) = 4
det(A) = 3
Tr(A) > 4*det(A)? = 0

# Dessiner les droites correspondants aux vecteurs propres et le sense des trajectoires associés
line_range = -1.5:.1:1.5;
eigenvector_1 = (V(2,1)/V(1,1)) * line_range;
eigenvector_2 = (V(2,2)/V(1,2)) * line_range;
hold on;
plot(line_range,eigenvector_1,"linewidth",10);
plot(line_range,eigenvector_2,"linewidth",10);
quiver([0;0],[0;0],V(1,:),V(2,:),"linewidth",10,"color","k");
legend("v_1","v_2","location","south");

# 2.Calculer les isoclines et dessiner des vecteurs vitesse sur ces isoclines
line_range = -1.5:.1:1.5;
isocline_1 = -(A(1,1)/A(1,2)) * line_range;
isocline_2 = -(A(2,1)/A(2,2)) * line_range;
hold on;
plot(line_range,isocline_1,"linewidth",10);
plot(line_range,isocline_2,"linewidth",10);
legend("isocline_1","isocline_2","location","southwest");

point_sequence = [line_range;isocline_1];
point_sequence = horzcat(point_sequence,[line_range;isocline_2]);
velocity_sequence = [];
for point_vec = point_sequence # Column-wise iteration
    velocity_vec = A*point_vec;
    velocity_sequence = [velocity_sequence,velocity_vec];
endfor

quiver(point_sequence(1,:),point_sequence(2,:),velocity_sequence(1,:),velocity_sequence(2,:),0.2,"linewidth",5,"color","k");
hold on;
plot(point_sequence(1,:),point_sequence(2,:),"ok")
plot(line_range,isocline_1,"linestyle",":","color","k");
plot(line_range,isocline_2,"linestyle",":","color","k");

# Define grid for plotting
x1range=-1.5:.1:1.5;
x2range=-1.5:.1:1.5;
[x1,x2] = meshgrid(x1range, x2range);

# Define the system to plot (based on matrix A)
x1p = A(1,1)*x1+A(1,2)*x2;
x2p = A(2,1)*x1+A(2,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);
# Eigenvectors plot
plot(line_range,eigenvector_1,"linewidth",10);
plot(line_range,eigenvector_2,"linewidth",10);
# Isoclines plot
plot(line_range,isocline_1,"linestyle",":","color","k");
plot(line_range,isocline_2,"linestyle",":","color","k");
grid on;
axis tight;

pkg load symbolic
syms x1(t) x2(t);
ode_sys = [diff(x1(t),t) == A(1,1)*x1(t) + A(1,2)*x2(t);  diff(x2(t),t) == A(2,1)*x1(t) + A(2,2)*x2(t)]
solutions = dsolve(ode_sys);
solutions{1}
solutions{2}

Symbolic pkg v2.7.1: Python communication link active, SymPy v1.3.
ode_sys = (sym 2×1 matrix)

  ⎡d                          ⎤
  ⎢──(x₁(t)) = 2⋅x₁(t) + x₂(t)⎥
  ⎢dt                         ⎥
  ⎢                           ⎥
  ⎢d                          ⎥
  ⎢──(x₂(t)) = x₁(t) + 2⋅x₂(t)⎥
  ⎣dt                         ⎦

ans = (sym)

              3⋅t       t
  x₁(t) = C₁⋅ℯ    + C₂⋅ℯ 

ans = (sym)

              3⋅t       t
  x₂(t) = C₁⋅ℯ    - C₂⋅ℯ 

Exercice 2

$ \mathbf{A} = \begin{bmatrix} -1 & -1 \\ 1 & -1 \end{bmatrix} $

1. Calculer les valeurs propres et vecteurs propres et en déduire le type du système
2. Calculer les isoclines et dessiner des vecteurs vitesse sur ces isoclines
3. Associer les signes des portions du plan délimitées par les isoclines.
4. Tracer des trajectoires qualitativement.

# 1. - Calcul valeur propres et vecteurs propres
A = [-1 -1; 1 -1]
[V,L] = eig(A)

A =

  -1  -1
   1  -1

V =

   0.70711 + 0.00000i   0.70711 - 0.00000i
   0.00000 - 0.70711i   0.00000 + 0.70711i

L =

Diagonal Matrix

  -1 + 1i        0
        0  -1 - 1i

printf("Tr(A) = %d\n",trace(A));
printf("det(A) = %d\n",det(A));
printf("Tr(A) > 4*det(A)? = %d\n",(trace(A) > 4*det(A)));

Tr(A) = -2
det(A) = 2
Tr(A) > 4*det(A)? = 0

# 2.Calculer les isoclines et dessiner des vecteurs vitesse sur ces isoclines
line_range = -1.5:.1:1.5;
isocline_1 = -(A(1,1)/A(1,2)) * line_range;
isocline_2 = -(A(2,1)/A(2,2)) * line_range;
hold on;
plot(line_range,isocline_1,"linewidth",5);
plot(line_range,isocline_2,"linewidth",5);
legend("isocline_1","isocline_2","location","south");

point_sequence = [line_range;isocline_1];
point_sequence = horzcat(point_sequence,[line_range;isocline_2]);
velocity_sequence = [];
for point_vec = point_sequence # Column-wise iteration
    velocity_vec = A*point_vec;
    velocity_sequence = [velocity_sequence,velocity_vec];
endfor

hold on;
quiver(point_sequence(1,:),point_sequence(2,:),velocity_sequence(1,:),velocity_sequence(2,:),0.2,"linewidth",5,"color","k");
plot(point_sequence(1,:),point_sequence(2,:),"ok");
plot(line_range,isocline_1,"linestyle",":","color","k");
plot(line_range,isocline_2,"linestyle",":","color","k");

# Define grid for plotting
x1range=-3:.1:3;
x2range=-3:.1:3;
[x1,x2] = meshgrid(x1range, x2range);

# Define the system to plot (based on matrix A)
x1p = A(1,1)*x1+A(1,2)*x2;
x2p = A(2,1)*x1+A(2,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);
# Trajectory plot
sys_str = sprintf('[%d*x(1)+%d*x(2);%d*x(1)+%d*x(2)]',A(1,1),A(1,2),A(2,1),A(2,2));
sys_ode = inline(sys_str,'t','x');
for x_0=0:0.1:1
 [ts,ys] = ode45(sys_ode,[0,10],[0;x_0]);
 plot(ys(:,1),ys(:,2))
end
# Isoclines plot
plot(line_range,isocline_1,"linestyle",":","color","k","linewidth",5);
plot(line_range,isocline_2,"linestyle",":","color","k","linewidth",5);
grid on;
axis tight;

pkg load symbolic
syms x1(t) x2(t);
ode_sys = [diff(x1(t),t) == A(1,1)*x1(t) + A(1,2)*x2(t);  diff(x2(t),t) == A(2,1)*x1(t) + A(2,2)*x2(t)]
solutions = dsolve(ode_sys);
solutions{1}
solutions{2}

ode_sys = (sym 2×1 matrix)

  ⎡d                         ⎤
  ⎢──(x₁(t)) = -x₁(t) - x₂(t)⎥
  ⎢dt                        ⎥
  ⎢                          ⎥
  ⎢d                         ⎥
  ⎢──(x₂(t)) = x₁(t) - x₂(t) ⎥
  ⎣dt                        ⎦

ans = (sym)

                                    -t
  x₁(t) = (-C₁⋅cos(t) - C₂⋅sin(t))⋅ℯ  

ans = (sym)

                                    -t
  x₂(t) = (-C₁⋅sin(t) + C₂⋅cos(t))⋅ℯ  

Exercice 3

$ \mathbf{A} = \begin{bmatrix} -2 & 0 \\ -4 & -2 \end{bmatrix} $

1. Calculer les valeurs propres et vecteurs propres et en déduire le type du système
2. Calculer les isoclines et dessiner des vecteurs vitesse sur ces isoclines
3. Associer les signes des portions du plan délimitées par les isoclines.
4. Tracer des trajectoires qualitativement.

# 1. - Calcul valeur propres et vecteurs propres
A = [-2 0; -4 -2]
[V,L] = eig(A)

A =

  -2   0
  -4  -2

V =

   0.00000   0.00000
   1.00000   1.00000

L =

Diagonal Matrix

  -2   0
   0  -2

printf("Tr(A) = %d\n",trace(A));
printf("det(A) = %d\n",det(A));
printf("Tr(A) > 4*det(A)? = %d\n",(trace(A) > 4*det(A)));

Tr(A) = -4
det(A) = 4
Tr(A) > 4*det(A)? = 0

# Dessiner les droites correspondants aux vecteurs propres et le sense des trajectoires associés
line_range = -1.5:.1:1.5;
eigenvector_1 = (V(2,1)/(V(1,1)+1e-170)) * line_range;
eigenvector_2 = (V(2,2)/(V(1,2)+1e-170)) * line_range;
hold on;
plot(line_range,eigenvector_1,"linewidth",10);
plot(line_range,eigenvector_2,"linewidth",10);
quiver([0;0],[0;0],V(1,:),V(2,:),"linewidth",10,"color","k")

# 2.Calculer les isoclines et dessiner des vecteurs vitesse sur ces isoclines
line_range = -1.5:.1:1.5;
isocline_1 = -(A(1,1)/(A(1,2)+1e-170)) * line_range;
isocline_2 = -(A(2,1)/A(2,2)) * line_range;
hold on;
plot(line_range,isocline_1,"linewidth",10);
plot(line_range,isocline_2,"linewidth",10);

point_sequence = [line_range;isocline_1];
point_sequence = horzcat(point_sequence,[line_range;isocline_2]);
velocity_sequence = [];
for point_vec = point_sequence # Column-wise iteration
    velocity_vec = A*point_vec;
    velocity_sequence = [velocity_sequence,velocity_vec];
endfor

hold on;
quiver(point_sequence(1,:),point_sequence(2,:),velocity_sequence(1,:),velocity_sequence(2,:),0.2,"linewidth",5,"color","k");
plot(point_sequence(1,:),point_sequence(2,:),"ok");
plot(line_range,isocline_1,"linestyle",":","color","k");
plot(line_range,isocline_2,"linestyle",":","color","k");

# Define grid for plotting
x1range=-1.5:.1:1.5;
x2range=-1.5:.1:1.5;
[x1,x2] = meshgrid(x1range, x2range);

# Define the system to plot (based on matrix A)
x1p = A(1,1)*x1+A(1,2)*x2;
x2p = A(2,1)*x1+A(2,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);
# Eigenvectors plot
plot(line_range,eigenvector_1,"linewidth",10);
plot(line_range,eigenvector_2,"linewidth",10);
# Isoclines plot
plot(line_range,isocline_1,"linestyle",":","color","k");
plot(line_range,isocline_2,"linestyle",":","color","k");
grid on;
axis tight;

pkg load symbolic
syms x1(t) x2(t);
ode_sys = [diff(x1(t),t) == A(1,1)*x1(t) + A(1,2)*x2(t);  diff(x2(t),t) == A(2,1)*x1(t) + A(2,2)*x2(t)]
solutions = dsolve(ode_sys);
solutions{1}
solutions{2}

ode_sys = (sym 2×1 matrix)

  ⎡     d                       ⎤
  ⎢     ──(x₁(t)) = 2⋅x₁(t)     ⎥
  ⎢     dt                      ⎥
  ⎢                             ⎥
  ⎢d                            ⎥
  ⎢──(x₂(t)) = 4⋅x₁(t) + 2⋅x₂(t)⎥
  ⎣dt                           ⎦

ans = (sym)

              2⋅t
  x₁(t) = C₂⋅ℯ   

ans = (sym)

                           2⋅t
  x₂(t) = (4⋅C₁ + 4⋅C₂⋅t)⋅ℯ