Dynamical Systems:

Natural Response

Table of Contents





1. Natural Response to Non-zero Initial Conditions

1.1. First-Order Ordinary Differential Equation

Consider the first-order ordinary differential equation (ODE):


$$ \frac{dx(t)}{dt} = kx(t), \quad x(0) = x_0 $$


The solution describes how the system responds over time from an initial condition $x_0$.


General solution:


$$ x(t) = x_0 e^{kt} $$


The solution is an exponential function, where the constant $x_0$ is determined by the initial condition.


Stability:

  • If $k > 0$, the system is unstable (exponential growth).
  • If $k < 0$, the system is stable (exponential decay).

In [ ]:
from IPython.display import YouTubeVideo
YouTubeVideo('4X0SGGrXDiI?start=399&end=459?rel=0', width = "560", height = "315")
Out[ ]:
In [ ]:
% plot an analytic solution

k = -0.9;
x0 = 2;

t = linspace(0,5,100);
x = x0*exp(k*t);

plot(t,x);
xlabel('time [sec]')
ylabel('x(t)')
ylim([0,2])

% but, all we did is just plotting (not computing)












No description has been provided for this image



















Time constant $\tau$:




$$
\dot{x} + \frac{1}{\tau}x = 0 \implies \dot{x} = -\frac{1}{\tau}x = kx
$$




$$
x(t) = x_0 e^{-t/\tau}
$$




    

  • Large $\tau$: slow response.
    • 

  • Small $\tau$: fast response.
    • 





After time $t = \tau$:




$$
\frac{x(\tau)}{x_0} = e^{-1} \approx 0.368
$$


















1.2. Two Independent First-Order ODEs
Consider two decoupled systems:




$$
\begin{aligned}
\dot{u}_1 &= \lambda_1 u_1, \quad u_1(0) = u_{10} \\\\
\dot{u}_2 &= \lambda_2 u_2, \quad u_2(0) = u_{20}
\end{aligned}
$$




Solutions:




$$
\begin{aligned}
u_1(t) &= u_{10} e^{\lambda_1 t} \\\\
u_2(t) &= u_{20} e^{\lambda_2 t}
\end{aligned}
$$




In vector form:




$$
u(t) =
\begin{bmatrix} u_1(t) \\ u_2(t) \end{bmatrix}, \quad
\dot{u}(t) =
\begin{bmatrix} \lambda_1 & 0 \\\\ 0 & \lambda_2 \end{bmatrix}
\begin{bmatrix} u_1(t) \\ u_2(t) \end{bmatrix} =
\Lambda u(t)
$$


















1.3. Coupled First-Order ODEs
For coupled variables:




$$
\begin{aligned}
\dot{u}_1 &= a_{11} u_1 + a_{12} u_2 \\\\
\dot{u}_2 &= a_{21} u_1 + a_{22} u_2
\end{aligned}
$$




In matrix form:




$$
\dot{u}(t) = A u(t), \quad A = \begin{bmatrix} a_{11} & a_{12} \\\\ a_{21} & a_{22} \end{bmatrix}
$$




General solution:




$$
u(t) = c_1 v_1 e^{\lambda_1 t} + c_2 v_2 e^{\lambda_2 t}
$$




where $\lambda_{1,2}$ and $v_{1,2}$ are eigenvalues and eigenvectors of $A$, and $c_1, c_2$ are determined by initial conditions.




Key ideas:


    

  • Natural response describes the system's behavior in the absence of external forcing.
    • 

  • Eigenvalues determine growth/decay rates and stability.
    • 

  • Time constant characterizes the speed of decay in first-order systems.
    • 






















2. Systems of Differential Equations
This section follows the matrix approach to systems of ODEs, inspired by MIT 18.03SC (Differential Equations) by Gilbert Strang.




2.1. General System and Superposition Principle



Consider the system of linear differential equations:




$$
\dot{\vec{u}} = A \vec{u}, \quad \vec{u}(0) = \vec{u}_0
$$




Superposition principle:

If $\vec{u}_1$ and $\vec{u}_2$ are solutions, any linear combination is also a solution:




$$
\begin{aligned}
\dot{\vec{u}} &= c_1 \dot{\vec{u}}_1 + c_2 \dot{\vec{u}}_2 \\\\
&= c_1 A \vec{u}_1 + c_2 A \vec{u}_2 \\\\
&= A(c_1 \vec{u}_1 + c_2 \vec{u}_2)
\end{aligned}
$$




For a scalar case:




$$
\dot{u} = a u \implies u(t) = c e^{a t}
$$




We look for a vector solution of the form:




$$
\vec{u}(t) = \vec{x} e^{\lambda t}
$$




Then:




$$
\dot{\vec{u}} = \lambda \vec{x} e^{\lambda t}, \quad A \vec{u} = A \vec{x} e^{\lambda t}
$$




This leads to the eigenvalue problem:




$$
A \vec{x} = \lambda \vec{x}
$$




Thus, solving linear ODEs is equivalent to finding eigenvalues and eigenvectors of $A$.
















2.2. Eigenanalysis and General Solution
Given:




$$
\dot{\vec{u}} = A \vec{u}, \quad \vec{u}(0) = \vec{u}_0
$$




Find eigenvalues $\lambda_i$ and eigenvectors $\vec{x}_i$:




$$
A \vec{x}_i = \lambda_i \vec{x}_i
$$




General solution:




$$
\vec{u}(t) = c_1 e^{\lambda_1 t} \vec{x}_1 + c_2 e^{\lambda_2 t} \vec{x}_2
$$




where $c_1$ and $c_2$ are determined from the initial condition:




$$
\vec{u}(0) = c_1 \vec{x}_1 + c_2 \vec{x}_2
$$




In matrix form:




$$
\vec{u}(0) =
\begin{bmatrix} \vec{x}_1 & \vec{x}_2 \end{bmatrix}
\begin{bmatrix} c_1 \\\\ c_2 \end{bmatrix}
$$
































2.3. Generalization to $n \times n$ Systems
For an $n \times n$ matrix $A$, the solution involves:




    

  • Eigenvalues $\lambda_i$, $i=1,\cdots,n$
    • 

  • Corresponding eigenvectors $\vec{x}_i$
    • 





Each $\vec{u}_i(t) = \vec{x}_i e^{\lambda_i t}$ solves the system.




Then, the general solution is:




$$
\vec{u}(t) = c_1 \vec{u}_1 + c_2 \vec{u}_2 + \cdots + c_n \vec{u}_n
$$




where $c_1, c_2, \ldots, c_n$ are determined by the initial condition $\vec{u}(0)$.


















2.4. ODE and Linear Algebra (2D Case)
Consider the system of linear ODEs:




$$
\dot{x}(t) = A x(t), \quad A \in \mathbb{R}^{2 \times 2}, \quad x(t) \in \mathbb{R}^2
$$




Find eigenvalues $\lambda_1, \lambda_2$ and corresponding eigenvectors $x_1, x_2$:




$$
A x_1 = \lambda_1 x_1, \quad A x_2 = \lambda_2 x_2
$$




Step 1: Build Eigenvector and Eigenvalue Matrices


Define:




$$
S = \begin{bmatrix} x_1 & x_2 \end{bmatrix}, \quad
\Lambda = \begin{bmatrix} \lambda_1 & 0 \\\\ 0 & \lambda_2 \end{bmatrix}
$$




These satisfy:




$$
A S = S \Lambda, \quad S^{-1} A S = \Lambda, \quad A = S \Lambda S^{-1}
$$




Step 2: Express Solution Using Decoupling


Consider the ansatz:




$$
x(t) = c_1(t) x_1
$$




Then:




$$
\dot{x}(t) = \dot{c}_1(t) x_1, \quad A x(t) = c_1(t) A x_1 = c_1(t) \lambda_1 x_1
$$




Thus:




$$
\dot{c}_1(t) x_1 = c_1(t) \lambda_1 x_1 \implies \dot{c}_1(t) = \lambda_1 c_1(t)
$$




which gives:




$$
c_1(t) = c_1(0) e^{\lambda_1 t}
$$




Similarly for $x_2$:




$$
c_2(t) = c_2(0) e^{\lambda_2 t}
$$




Step 3: General Solution in Matrix Form




$$
x(t) = c_1(0) e^{\lambda_1 t} x_1 + c_2(0) e^{\lambda_2 t} x_2
$$




Expressed using $S$ and $\Lambda$:




$$
A S = S \Lambda \implies A = S \Lambda S^{-1}
$$




Introduce:




$$
v(t) = S^{-1} x(t), \quad x(t) = S v(t)
$$




Then:




$$
\dot{x}(t) = S \dot{v}(t), \quad A x(t) = A S v(t) = S \Lambda v(t)
$$




This leads to:




$$
S \dot{v}(t) = S \Lambda v(t) \implies \dot{v}(t) = \Lambda v(t)
$$




Step 4: Decoupled Equations in $v$-Frame




$$
v(t) = \begin{bmatrix} v_1(0) e^{\lambda_1 t} \\\\ v_2(0) e^{\lambda_2 t} \end{bmatrix}
$$




where:




$$
v(0) = S^{-1} x(0)
$$




Thus:




$$
x(t) = S v(t) = S e^{\Lambda t} v(0)
$$




where:




$$
e^{\Lambda t} = \begin{bmatrix} e^{\lambda_1 t} & 0 \\\\ 0 & e^{\lambda_2 t} \end{bmatrix}
$$




Step 5: Final General Solution




$$
x(t) = S e^{\Lambda t} S^{-1} x(0)
$$




which gives a closed-form expression for the system’s evolution.




Summary


    

  • Eigenvalues $\lambda_1,\lambda_2$ determine growth/decay rates and stability.
    • 

  • Eigenvectors $x_1,x_2$ define the system's natural modes.
    • 

  • The system decouples into independent ODEs in the eigenvector frame ($v$-frame).
    • 

  • The original solution is reconstructed via $x(t) = S e^{\Lambda t} S^{-1} x(0)$.
    • 















In [ ]:






from IPython.display import YouTubeVideo
YouTubeVideo('YUjdyKhWt6E', width = "560", height = "315")















Out[ ]:



















In [ ]:






from IPython.display import YouTubeVideo
YouTubeVideo('lXNXrLcoerU', width = "560", height = "315")















Out[ ]:


























3. Real Eigenvalues
The general solution for a system with real eigenvalues is:




$$
\vec{u}(t) = c_1 \vec{x}_1 e^{\lambda_1 t} + c_2 \vec{x}_2 e^{\lambda_2 t}
$$




where:




$$
\begin{bmatrix}
c_1 \\ c_2
\end{bmatrix}
=
\begin{bmatrix}
\vec{x}_1 & \vec{x}_2
\end{bmatrix}^{-1}
\begin{bmatrix}
u_1(0) \\ u_2(0)
\end{bmatrix}
$$
















3.1. Examples
3.1.1. Example 1
Consider:




$$
\begin{align*}
\dot{u_1} &= -3u_1 + u_2 \\
\dot{u_2} &= u_1 - 3u_2
\end{align*}
$$




$$
\dot{\vec{u}} =
\begin{bmatrix}
-3 & 1\\
1 & -3
\end{bmatrix}
\vec{u} \qquad \text{where} \quad
A =
\begin{bmatrix}
-3 & 1\\
1 & -3
\end{bmatrix}
$$




(1) Eigenvalues




$$
\begin{align*}
\begin{vmatrix}
-3-\lambda & 1\\
1 & -3-\lambda
\end{vmatrix} &=
(3 + \lambda)^2 - 1 = \lambda^2 + 6\lambda + 8 \\
&= (\lambda + 2)(\lambda + 4) = 0 \quad\quad \\ \\
&\therefore \lambda_1 = -2,\; \lambda_2 = -4
\end{align*}
$$




(2) Eigenvectors




$$
\begin{align*}
\begin{bmatrix}
-1 & 1\\
1 & -1
\end{bmatrix}
\begin{bmatrix}
x\\
y
\end{bmatrix} =
\begin{bmatrix}
0\\
0
\end{bmatrix}
&\implies
\vec{x}_1 =
\begin{bmatrix}
1\\
1
\end{bmatrix} \qquad \text{for } \lambda_1\\\\
\begin{bmatrix}
1 & 1\\
1 & 1
\end{bmatrix}
\begin{bmatrix}
x\\
y
\end{bmatrix} =
\begin{bmatrix}
0\\
0
\end{bmatrix}
&\implies
\vec{x}_2 =
\begin{bmatrix}
1\\
-1
\end{bmatrix}  \qquad \text{for } \lambda_2
\end{align*}
$$




(3) General solution $\vec{u}(t)$




$$
\begin{align*}
\vec{u}(t) &= c_1\vec{x}_1e^{\lambda_1 t} + c_2\vec{x}_2e^{\lambda_2 t} \\ \\
&= c_1
\begin{bmatrix}
1 \\
1
\end{bmatrix} e^{-2t} + c_2
\begin{bmatrix}
1 \\
-1
\end{bmatrix} e^{-4t} =
\begin{bmatrix}
c_1e^{-2t} + c_2e^{-4t} \\
c_1e^{-2t} - c_2e^{-4t} \\
\end{bmatrix}
\end{align*}
$$




(4) Initial Condition




$$u(0) = c_1
\begin{bmatrix}
1 \\
1
\end{bmatrix} + c_2
\begin{bmatrix}
1 \\
-1
\end{bmatrix}$$


















In the eigenvector basis ($v$-coordinates)




$$
\begin{align*}
\dot{\vec{v}} =
\begin{bmatrix}
-2 & 0 \\
0 & -4
\end{bmatrix} \vec{v} \qquad \text{where} \quad
&\Lambda =
\begin{bmatrix}
\lambda_1 & 0 \\
0 & \lambda_2
\end{bmatrix} \\\\
&\vec{x}_1' =
\begin{bmatrix}
1 \\
0
\end{bmatrix}, \quad \vec{x}_2' =
\begin{bmatrix}
0 \\
1
\end{bmatrix}  
\end{align*}
$$




$$
\begin{align*}
v_1(t) &= v_1(0)e^{\lambda_1 t} = v_1(0)e^{-2t} \\
v_2(t) &= v_2(0)e^{\lambda_2 t} = v_2(0)e^{-4t}
\end{align*}
$$




    

  • Relationship between $\vec{u}$ and $\vec{v}$
    • 





$$
\begin{align*}
u &= Sv\\
AS &= S\Lambda \quad \text{where} \quad
S =
\begin{bmatrix}
\vec{x}_1 & \vec{x}_2
\end{bmatrix} =
\begin{bmatrix}
1 & 1 \\
1 & -1
\end{bmatrix}
\end{align*}
$$




$$
\vec{u}(t) = S v(t) = S e^{\Lambda t} S^{-1} \vec{u}(0)
$$
















    

  • Simulation
    • 



$$
\dot{\vec{u}} =
\begin{bmatrix}
-3 & 1\\
1 & -3
\end{bmatrix}
\vec{u}, \qquad u(0) =
\begin{bmatrix}
4\\
8
\end{bmatrix}
$$












In [ ]:






%plot -s 900,350

A = [-3  1;
     1 -3];

%% eigen-analysis
[S,D] = eig(A);

[lamb,idx] = sort(diag(D),'descend');
S = S(:,idx);

x1 = S(:,1);
x2 = S(:,2);

u0 = [4;8];
C = inv(S)*u0;

t = 0:0.01:3;
u = C(1)*x1*exp(lamb(1)*t) + C(2)*x2*exp(lamb(2)*t);

% plot u1 and u2 as a function of time
plot(t,u(1,:),t,u(2,:))
xlabel('time')
legend('u_1', 'u_2')






























No description has been provided for this image

















Phase portrait













In [ ]:






from IPython.display import YouTubeVideo
YouTubeVideo('e3FfmXtkppM', width = "560", height = "315")















Out[ ]:



















In [ ]:






%plot -s 560,420

% plot eigenvectors (X1 and X2)
k = -20:0.1:20;
y1 = S(:,1)*k;
y2 = S(:,2)*k;

plot(y1(1,:), y1(2,:)); hold on
plot(y2(1,:), y2(2,:));
xlabel('u_1', 'fontsize', 12)
ylabel('u_2', 'fontsize', 12)
title('Trajectory', 'fontsize', 12)
axis equal
axis([-10 10 -10 10])

% plot a trajectory of u1 and u2
for i = 1:length(t)
    plot(u(1,i), u(2,i), 'k.');
end
hold off






























No description has been provided for this image

















3.1.2. Example 2
Consider the decoupled system:




$$\dot{\vec{u}} =
\begin{bmatrix}
-2 & 0 \\
0 & -4
\end{bmatrix}\vec{u}, \qquad u(0) =
\begin{bmatrix}
4\\
8
\end{bmatrix}
$$




    

  • Eigenvalues:
    • 





$$
\begin{vmatrix}
-2 -\lambda & 0 \\
0 & -4 - \lambda
\end{vmatrix} = (\lambda + 2)(\lambda + 4) = 0 \quad \implies \quad \lambda_1 = -2, \lambda_2 = -4
$$




    

  • Eigenvectors:
    • 





$$
\begin{align*}
\begin{bmatrix}
0 & 0 \\
0 & -2
\end{bmatrix}
\begin{bmatrix}
x\\
y
\end{bmatrix} =
\begin{bmatrix}
0\\
0
\end{bmatrix} \quad &\implies \quad \vec{x}_1 =
\begin{bmatrix}
1\\
0
\end{bmatrix}\\
\begin{bmatrix}
2 & 0 \\
0 & 0
\end{bmatrix}
\begin{bmatrix}
x\\
y
\end{bmatrix} =
\begin{bmatrix}
0\\
0
\end{bmatrix} \quad &\implies \quad \vec{x}_2 =
\begin{bmatrix}
0\\
1
\end{bmatrix}
\end{align*}
$$




The system is already diagonalized, and solutions are:




$$
\vec{u}(t) =
\begin{bmatrix}
4 e^{-2t} \\\\ 8 e^{-4t}
\end{bmatrix}
$$














In [ ]:






%plot -s 900,350

A = [-2  0;
     0 -4];

%% eigen-analysis
[S,D] = eig(A);

[lamb,idx] = sort(diag(D),'descend');
S = S(:,idx);

x1 = S(:,1);
x2 = S(:,2);

u0 = [4;8];
C = inv(S)*u0;

t = 0:0.01:3;
u = C(1)*x1*exp(lamb(1)*t) + C(2)*x2*exp(lamb(2)*t);

% plot u1 and u2 as a function of time
plot(t,u(1,:),t,u(2,:))
xlabel('time')






























No description has been provided for this image















In [ ]:






%plot -s 560,420

% plot eigenvectors (X1 and X2)
k = -20:0.1:20;
y1 = S(:,1)*k;
y2 = S(:,2)*k;

plot(y1(1,:), y1(2,:)); hold on
plot(y2(1,:), y2(2,:));
xlabel('u_1', 'fontsize', 12)
ylabel('u_2', 'fontsize', 12)
title('Trajectory', 'fontsize', 12)
axis equal
axis([-10 10 -10 10])

% plot a trajectory of u1 and u2
for i = 1:length(t)
    plot(u(1,i), u(2,i), 'k.');
end
hold off






























No description has been provided for this image

















3.1.3. Example 3
Consider:




$$\dot{\vec{u}} =
\begin{bmatrix}
-2 & 0 \\
1 & -4
\end{bmatrix}\vec{u}, \qquad u(0) =
\begin{bmatrix}
4\\
8
\end{bmatrix}
$$




    

  • Eigenvalues:
    • 





$$
\begin{vmatrix}
-2 -\lambda & 0 \\
1 & -4 - \lambda
\end{vmatrix} = (\lambda + 2)(\lambda + 4) = 0 \quad \implies \quad \lambda_1 = -2, \lambda_2 = -4
$$




    

  • Eigenvectors:
    • 





$$
\begin{align*}
\begin{bmatrix}
0 & 0 \\
1 & -2
\end{bmatrix}
\begin{bmatrix}
x\\
y
\end{bmatrix} =
\begin{bmatrix}
0\\
0
\end{bmatrix} \quad &\implies \quad \vec{x}_1 =
\begin{bmatrix}
2\\
1
\end{bmatrix}\\
\begin{bmatrix}
2 & 0 \\
1 & 0
\end{bmatrix}
\begin{bmatrix}
x\\
y
\end{bmatrix} =
\begin{bmatrix}
0\\
0
\end{bmatrix} \quad &\implies \quad \vec{x}_2 =
\begin{bmatrix}
0\\
1
\end{bmatrix}
\end{align*}
$$














In [ ]:






%plot -s 900,350

A = [-2  0;
     1 -4];

%% eigen-analysis
[S,D] = eig(A);

[lamb,idx] = sort(diag(D),'descend');
S = S(:,idx);

x1 = S(:,1);
x2 = S(:,2);

u0 = [4;8];
C = inv(S)*u0;

t = 0:0.01:3;
u = C(1)*x1*exp(lamb(1)*t) + C(2)*x2*exp(lamb(2)*t);

% plot u1 and u2 as a function of time
plot(t,u(1,:),t,u(2,:))
xlabel('time')






























No description has been provided for this image















In [ ]:






%plot -s 560,420

% plot eigenvectors (X1 and X2)
k = -20:0.1:20;
y1 = S(:,1)*k;
y2 = S(:,2)*k;

plot(y1(1,:), y1(2,:)); hold on
plot(y2(1,:), y2(2,:));
xlabel('u_1', 'fontsize', 12)
ylabel('u_2', 'fontsize', 12)
title('Trajectory', 'fontsize', 12)
axis equal
axis([-10 10 -10 10])

% plot a trajectory of u1 and u2
for i = 1:length(t)
    plot(u(1,i), u(2,i), 'k.');
end
hold off






























No description has been provided for this image



















Summary






























3.2. Different Eigenvectors with the Same Eigenvalues
In linear systems, a matrix $A$ may have a repeated eigenvalue $\lambda$. Depending on its structure, $A$ can have:


    

  • Multiple linearly independent eigenvectors (diagonalizable)
    • 

  • Only one eigenvector (defective)
    • 





Case 1: Repeated Eigenvalues with Multiple Independent Eigenvectors


If $A$ has a repeated eigenvalue $\lambda$ with two linearly independent eigenvectors $x_1, x_2$, then:




$$
A x_1 = \lambda x_1, \quad A x_2 = \lambda x_2
$$




The general solution is:




$$
u(t) = c_1 x_1 e^{\lambda t} + c_2 x_2 e^{\lambda t}
$$




Here, $A$ is diagonalizable, and the system behaves like a standard real-eigenvalue system with repeated $\lambda$.




Case 2: Repeated Eigenvalues with Only One Eigenvector (Defective Case)


If $A$ is defective (non-diagonalizable), it has a repeated eigenvalue $\lambda$ but only one eigenvector $x_1$.




To complete the solution, we introduce a generalized eigenvector $x_g$ satisfying:




$$
(A - \lambda I) x_g = x_1
$$




The general solution becomes:




$$
u(t) = c_1 x_1 e^{\lambda t} + c_2 (x_g + t x_1) e^{\lambda t}
$$




This introduces a $t e^{\lambda t}$ term, indicating a defective mode with linear growth along the eigenspace.




3.2.1. Example: Defective Matrix with Repeated Eigenvalue
Consider:




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




    

  • Eigenvalues:
    • 





$$
\det(A - \lambda I) = (2 - \lambda)^2 = 0 \implies \lambda = 2
$$




    

  • Eigenvectors:
    • 





$$
(A - 2I) x = 0 \implies \begin{bmatrix} 0 & 1 \\\\ 0 & 0 \end{bmatrix} \begin{bmatrix} x_1 \\\\ x_2 \end{bmatrix} = 0
$$




$$
x_2 = 0, \quad x_1 \text{ free} \implies x_1 = \begin{bmatrix} 1 \\\\ 0 \end{bmatrix}
$$




    

  • Generalized Eigenvector:
    • 





$$
(A - 2I) x_g = x_1
$$




$$
\begin{bmatrix} 0 & 1 \\\\ 0 & 0 \end{bmatrix} \begin{bmatrix} a \\\\ b \end{bmatrix} = \begin{bmatrix} 1 \\\\ 0 \end{bmatrix}
$$




$$
b = 1, \quad a \text{ arbitrary} \implies x_g = \begin{bmatrix} 0 \\\\ 1 \end{bmatrix}
$$




    

  • General Solution:
    • 





$$
u(t) = c_1 \begin{bmatrix} 1 \\\\ 0 \end{bmatrix} e^{2t} + c_2 \left( \begin{bmatrix} 0 \\\\ 1 \end{bmatrix} + t \begin{bmatrix} 1 \\\\ 0 \end{bmatrix} \right) e^{2t}
$$





















4. De-coupling through Linear Transformation
We have studied how to solve linear ODEs. In this section, we will examine the physical meaning of decoupling in detail using the following example.
































Consider the system:




$$
\begin{aligned}
\dot{u}_1 &= -2(u_1 - u_2) \\
\dot{u}_2 &= (u_1 - u_2)
\end{aligned}
$$




In matrix form:




$$
\dot{\vec{u}} =
\begin{bmatrix}
-2 & 2 \\\\
1 & -1
\end{bmatrix} \vec{u}
$$




Step 1: Eigenvalues and Eigenvectors


Compute eigenvalues:




$$
\begin{aligned}
\det(A - \lambda I) &= (-2-\lambda)(-1-\lambda) - 2 = \lambda^2 + 3\lambda \\\\
&= \lambda(\lambda + 3) \implies \lambda_1=0,\quad \lambda_2=-3
\end{aligned}
$$




For $\lambda_1=0$:




$$
\begin{bmatrix}
-2 & 2 \\\\
1 & -1
\end{bmatrix}
\begin{bmatrix}
x \\\\ y
\end{bmatrix} = 0 \implies x = y \implies \vec{x}_1 = \begin{bmatrix} 1 \\\\ 1 \end{bmatrix}
$$




For $\lambda_2=-3$:




$$
\begin{bmatrix}
1 & 2 \\\\
1 & 2
\end{bmatrix}
\begin{bmatrix}
x \\\\ y
\end{bmatrix} = 0 \implies x = -2y \implies \vec{x}_2 = \begin{bmatrix} 2 \\\\ -1 \end{bmatrix}
$$




Step 2: Interpretation of Modes


    

  • $\lambda_1=0$ corresponds to a constant mode ($e^{0t}=1$), representing a conserved quantity.
    

    • 

  • $\lambda_2=-3$ corresponds to a decaying mode ($e^{-3t}$), representing a dissipative effect.
    

    • 





Step 3: Change of Variables to Decouple


Total amount of water:




$$
v_1 = u_1 + 2u_2
$$




Difference in height:




$$
v_2 = u_1 - u_2
$$




Compute derivatives:




$$
\begin{aligned}
\dot{v}_1 &= \dot{u}_1 + 2\dot{u}_2 = -2(u_1 - u_2) + 2(u_1 - u_2) = 0 \\\\
\dot{v}_2 &= \dot{u}_1 - \dot{u}_2 = -2(u_1 - u_2) - (u_1 - u_2) = -3(u_1 - u_2) = -3v_2
\end{aligned}
$$




Thus:




$$
\begin{aligned}
\dot{v}_1 &= 0 \implies v_1(t) = v_1(0) \quad \text{(constant)} \\\\
\dot{v}_2 &= -3v_2 \implies v_2(t) = v_2(0) e^{-3t} \quad \text{(decay)}
\end{aligned}
$$




Step 4: Interpretation


    

  • $v_1$ is invariant over time, representing the total water, which is conserved.
    

    • 

  • $v_2$ decays exponentially to zero, representing the difference in water levels, which equalizes over time.
    

    • 





Key Takeaways


    

  • Eigenvalue $\lambda_1=0$ indicates a conservation law in the system.
    

    • 

  • Eigenvalue $\lambda_2=-3$ reveals an exponentially decaying mode.
    

    • 

  • Changing variables to align with eigenvectors decouples the system, simplifying analysis.
    

    • 















In [ ]:






%plot -s 900,350

A = [-2  2;
     1 -1];

%% eigen-analysis
[S,D] = eig(A);

[lamb,idx] = sort(diag(D),'descend');
S = S(:,idx);

x1 = S(:,1);
x2 = S(:,2);

u0 = [0;6];
C = inv(S)*u0;

t = 0:0.01:3;
u = C(1)*x1*exp(lamb(1)*t) + C(2)*x2*exp(lamb(2)*t);

% plot u1 and u2 as a function of time
plot(t,u(1,:),t,u(2,:))
xlabel('time')
legend('u_1', 'u_2')






























No description has been provided for this image















In [ ]:






%plot -s 560,420

% plot eigenvectors (X1 and X2)
k = -20:0.1:20;
y1 = S(:,1)*k;
y2 = S(:,2)*k;

plot(y1(1,:), y1(2,:)); hold on
plot(y2(1,:), y2(2,:));
xlabel('u_1', 'fontsize', 12)
ylabel('u_2', 'fontsize', 12)
title('Trajectory', 'fontsize', 12)
axis equal
axis([-10 10 -10 10])

% plot a trajectory of u1 and u2
for i = 1:length(t)
    plot(u(1,i), u(2,i), 'k.');
end
hold off






























No description has been provided for this image















In [ ]:






%plot -s 900,350

v = inv(S)*u;

plot(t,v(1,:),t,v(2,:))
xlabel('time')
legend('v_1', 'v_2')






























No description has been provided for this image















In [ ]:






%plot -s 560,420

q1 = inv(S)*y1;
q2 = inv(S)*y2;

plot(q1(1,:), q1(2,:)); hold on
plot(q2(1,:), q2(2,:));

for i = 1:length(t)
    plot(v(1,i), v(2,i), 'k.'), hold on
end
hold off
axis equal
axis([-10 10 -10 10])
xlabel('v_1', 'fontsize', 12)
ylabel('v_2', 'fontsize', 12)
title('Trajectory', 'fontsize', 12)






























No description has been provided for this image






















5. Complex Eigenvalues and Oscillatory Behavior
In dynamic systems, the eigenvalues $\lambda$ of the system matrix can be complex:




$$
\lambda = \sigma + i \omega,
$$




where $\sigma = \operatorname{Re}(\lambda)$ determines exponential growth or decay, and $\omega = \operatorname{Im}(\lambda)$ determines oscillation frequency.




The exponential of a complex eigenvalue appears as:




$$
\begin{align*}
e^{\lambda t} &= e^{(\sigma + i \omega)t} = e^{\sigma t} e^{i \omega t} = e^{\sigma t}(\cos \omega t + i \sin \omega t).
\end{align*}
$$




    

  • The factor $e^{\sigma t}$ describes exponential decay or growth.
    • 

  • The factor $e^{i \omega t}$ describes oscillations with angular frequency $\omega$.
    • 





The behavior depends on the real part of the eigenvalue:




$$
\begin{aligned}
\operatorname{Re}(\lambda) = \sigma < 0 &\implies \text{stable} \ (\text{decaying oscillations}), \\
\operatorname{Re}(\lambda) = \sigma > 0 &\implies \text{unstable} \ (\text{growing oscillations}).
\end{aligned}
$$




Thus, complex eigenvalues correspond to oscillatory solutions, with the decay (or growth) rate determined by $\sigma$ and oscillation frequency by $\omega$.




















5.1. Example 1
Consider the system:




$$
\dot{\vec{u}} =
\begin{bmatrix}
0 & 1 \\
-1 & 0
\end{bmatrix}\vec{u}, \quad \vec{u}(0) =
\begin{bmatrix}
1 \\ 0
\end{bmatrix}.
$$




The characteristic equation is:




$$
\begin{vmatrix}
-\lambda & 1 \\
-1 & -\lambda
\end{vmatrix} = \lambda^2 + 1 = 0,
$$




which yields complex conjugate eigenvalues:




$$
\lambda_1 = i, \quad \lambda_2 = -i.
$$




Find the corresponding eigenvectors:




$$
\begin{aligned}
\begin{bmatrix}
-i & 1 \\
-1 & -i
\end{bmatrix}
\begin{bmatrix}
x \\ y
\end{bmatrix} = 0 \implies \vec{x}_1 =
\begin{bmatrix}
1 \\ i
\end{bmatrix}, \quad
\begin{bmatrix}
-i & 1 \\
-1 & i
\end{bmatrix}
\begin{bmatrix}
x \\ y
\end{bmatrix} = 0 \implies \vec{x}_2 =
\begin{bmatrix}
1 \\ -i
\end{bmatrix}.
\end{aligned}
$$




Note that the system matrix represents a rotation:




$$
\begin{bmatrix}
0 & 1 \\
-1 & 0
\end{bmatrix} =
\begin{bmatrix}
\cos(-\frac{\pi}{2}) & -\sin(-\frac{\pi}{2}) \\
\sin(-\frac{\pi}{2}) & \cos(-\frac{\pi}{2})
\end{bmatrix} = R\left(-\frac{\pi}{2}\right).
$$




The general solution is:




$$
\begin{aligned}
\vec{u}(t) &= c_1 \vec{x}_1 e^{i t} + c_2 \vec{x}_2 e^{-i t} \\
&= c_1
\begin{bmatrix}
1 \\ i
\end{bmatrix} e^{i t} + c_2
\begin{bmatrix}
1 \\ -i
\end{bmatrix} e^{-i t}.
\end{aligned}
$$




Expanding the terms:




$$
\vec{u}(t) =
\begin{bmatrix}
c_1(\cos t + i \sin t) + c_2(\cos t - i \sin t) \\
c_1(i \cos t - \sin t) + c_2(-i \cos t - \sin t)
\end{bmatrix}.
$$




Apply the initial condition $\vec{u}(0) = [1,0]^T$:




$$
\vec{u}(0) =
\begin{bmatrix}
c_1 + c_2 \\
i c_1 - i c_2
\end{bmatrix} =
\begin{bmatrix}
1 \\ 0
\end{bmatrix}
\implies c_1 = c_2 = \frac{1}{2}.
$$




Thus,




$$
\vec{u}(t) =
\begin{bmatrix}
\frac{1}{2}(e^{i t} + e^{-i t}) \\
-\frac{1}{2i}(e^{i t} - e^{-i t})
\end{bmatrix} =
\begin{bmatrix}
\cos t \\
-\sin t
\end{bmatrix}.
$$




This describes a circular trajectory in the $u_1$-$u_2$ plane.






Physical Interpretation


This system models simple harmonic motion:




$$
\ddot{u} + \omega_n^2 u = 0, \quad \omega_n = 1,
$$




where:




$$
u_1 = u \quad (\text{position}), \quad u_2 = \dot{u} \quad (\text{velocity}).
$$




The system is:




$$
\begin{aligned}
\begin{bmatrix}
\dot{u}_1 \\ \dot{u}_2
\end{bmatrix} =
\begin{bmatrix}
u_2 \\ -u_1
\end{bmatrix}
\implies \dot{\vec{u}} =
\begin{bmatrix}
0 & 1 \\
-1 & 0
\end{bmatrix} \vec{u}.
\end{aligned}
$$




Thus, $\lambda_{1,2} = \pm i$ correspond to purely imaginary eigenvalues, indicating undamped oscillations with angular frequency $\omega_n=1$.














In [ ]:






%plot -s 900,350

A = [0  1;
     -1 0];

%% eigen-analysis
[S,D] = eig(A);

lamb = diag(D)

x1 = S(:,1);
x2 = S(:,2);

u0 = [0;6];
C = inv(S)*u0;

t = 0:0.01:10;
u = C(1)*x1*exp(lamb(1)*t) + C(2)*x2*exp(lamb(2)*t);

% plot u1 and u2 as a function of time
plot(t,u(1,:),t,u(2,:),t,zeros(size(t)),'k--')
xlabel('time')
legend('u_1', 'u_2')



















lamb =

   0.0000 + 1.0000i
   0.0000 - 1.0000i














No description has been provided for this image















In [ ]:






%plot -s 900,420

subplot(1,2,1)
% plot a trajectory of u1 and u2
for i = 1:length(t)
    plot(u(1,i), u(2,i), 'k.'), hold on
end
hold off

axis equal
axis([-12 12 -12 12])
grid on
xlabel('u_1', 'fontsize', 12)
ylabel('u_2', 'fontsize', 12)
title('Trajectory', 'fontsize', 12)

subplot(1,2,2)
plot(real(lamb),imag(lamb),'o')
axis equal
grid on
axis([-4,4,-6,6])
xlabel('Re', 'fontsize', 12)
ylabel('Im', 'fontsize', 12)
title('Eigenvalues','fontsize',12)






























No description has been provided for this image





















5.2. Example 2



Consider the system:




$$
\dot{\vec{u}} =
\begin{bmatrix}
0 & 1 \\
-4 & 0
\end{bmatrix} \vec{u}, \quad \vec{u}(0) =
\begin{bmatrix}
1 \\ 0
\end{bmatrix}.
$$




The characteristic equation is:




$$
\begin{vmatrix}
-\lambda & 1 \\
-4 & -\lambda
\end{vmatrix} = \lambda^2 + 4 = 0,
$$




which gives purely imaginary eigenvalues:




$$
\lambda_1 = 2i, \quad \lambda_2 = -2i.
$$




These indicate oscillations with angular velocity $\omega = 2$.




Find the corresponding eigenvectors:




$$
\begin{aligned}
\begin{bmatrix}
-2i & 1 \\
-4 & -2i
\end{bmatrix}
\begin{bmatrix}
x \\ y
\end{bmatrix} = 0 \implies \vec{x}_1 =
\begin{bmatrix}
1 \\ 2i
\end{bmatrix}, \quad
\begin{bmatrix}
2i & 1 \\
-4 & 2i
\end{bmatrix}
\begin{bmatrix}
x \\ y
\end{bmatrix} = 0 \implies \vec{x}_2 =
\begin{bmatrix}
1 \\ -2i
\end{bmatrix}.
\end{aligned}
$$




The general solution is:




$$
\begin{aligned}
\vec{u}(t) &= c_1 \vec{x}_1 e^{i 2t} + c_2 \vec{x}_2 e^{-i 2t} \\
&= c_1
\begin{bmatrix}
1 \\ 2i
\end{bmatrix} e^{i 2t} + c_2
\begin{bmatrix}
1 \\ -2i
\end{bmatrix} e^{-i 2t}.
\end{aligned}
$$




Applying the initial condition $\vec{u}(0) = [1, 0]^T$:




$$
\vec{u}(0) =
\begin{bmatrix}
c_1 + c_2 \\
2i(c_1 - c_2)
\end{bmatrix} =
\begin{bmatrix}
1 \\ 0
\end{bmatrix} \implies c_1 = c_2 = \frac{1}{2}.
$$




Thus,




$$
\vec{u}(t) =
\begin{bmatrix}
\frac{1}{2}(e^{i 2t} + e^{-i 2t}) \\
-2 \frac{1}{2i}(e^{i 2t} - e^{-i 2t})
\end{bmatrix} =
\begin{bmatrix}
\cos 2t \\
-2 \sin 2t
\end{bmatrix}.
$$




This represents a faster oscillation than in Example 1, with angular frequency $\omega = 2$ and amplitude ratio $2:1$ between $u_2$ and $u_1$.














In [ ]:






%plot -s 900,350

A = [0  1;
     -4 0];

%% eigen-analysis
[S,D] = eig(A);

lamb = diag(D)

x1 = S(:,1);
x2 = S(:,2);

u0 = [0;6];
C = inv(S)*u0;

t = 0:0.01:10;
u = C(1)*x1*exp(lamb(1)*t) + C(2)*x2*exp(lamb(2)*t);

% plot u1 and u2 as a function of time
plot(t,u(1,:),t,u(2,:),t,zeros(size(t)),'k--')
xlabel('time')
legend('u_1', 'u_2')



















lamb =

   0.0000 + 2.0000i
   0.0000 - 2.0000i














No description has been provided for this image















In [ ]:






%plot -s 900,420

subplot(1,2,1)
% plot a trajectory of u1 and u2
for i = 1:length(t)
    plot(u(1,i), u(2,i), 'k.'), hold on
end
hold off

axis equal
axis([-12 12 -12 12])
grid on
xlabel('u_1', 'fontsize', 12)
ylabel('u_2', 'fontsize', 12)
title('Trajectory', 'fontsize', 12)

subplot(1,2,2)
plot(real(lamb),imag(lamb),'o')
axis equal
grid on
axis([-4,4,-6,6])
xlabel('Re', 'fontsize', 12)
ylabel('Im', 'fontsize', 12)
title('Eigenvalues','fontsize',12)






























No description has been provided for this image





















5.3. Example 3
Consider the system:




$$
\dot{\vec{u}} =
\begin{bmatrix}
3 & -9 \\
4 & -3
\end{bmatrix} \vec{u}.
$$




The characteristic equation is:




$$
\begin{vmatrix}
3 - \lambda & -9 \\
4 & -3 - \lambda
\end{vmatrix} = (3 - \lambda)(-3 - \lambda) - (-9)(4)
= (\lambda^2 - 9) + 36 = \lambda^2 + 27.
$$




Solving for $\lambda$:




$$
\lambda^2 + 27 = 0 \implies \lambda = \pm 3\sqrt{3} i.
$$




Thus, the system has purely imaginary eigenvalues:




$$
\lambda_1 = 3\sqrt{3} i, \quad \lambda_2 = -3\sqrt{3} i,
$$




which indicate oscillatory solutions with angular frequency $\omega = 3\sqrt{3}$.














In [ ]:






%plot -s 900,350

A = [3  9;
     -4 -3];

%% eigen-analysis
[S,D] = eig(A);

lamb = diag(D)

x1 = S(:,1);
x2 = S(:,2);

u0 = [0;6];
C = inv(S)*u0;

t = 0:0.01:10;
u = C(1)*x1*exp(lamb(1)*t) + C(2)*x2*exp(lamb(2)*t);

% plot u1 and u2 as a function of time
plot(t,u(1,:),t,u(2,:),t,zeros(size(t)),'k--')
xlabel('time')
legend('u_1', 'u_2')



















lamb =

   0.0000 + 5.1962i
   0.0000 - 5.1962i














No description has been provided for this image















In [ ]:






%plot -s 900,420

subplot(1,2,1)
% plot a trajectory of u1 and u2
for i = 1:length(t)
    plot(u(1,i), u(2,i), 'k.'), hold on
end
hold off

axis equal
axis([-12 12 -12 12])
grid on
xlabel('u_1', 'fontsize', 12)
ylabel('u_2', 'fontsize', 12)
title('Trajectory', 'fontsize', 12)

subplot(1,2,2)
plot(real(lamb),imag(lamb),'o')
axis equal
grid on
axis([-4,4,-6,6])
xlabel('Re', 'fontsize', 12)
ylabel('Im', 'fontsize', 12)
title('Eigenvalues','fontsize',12)






























No description has been provided for this image






















6. Complex Eigenvalues with Damping
6.1. Example
Consider the system:




$$
\dot{\vec{u}} =
\begin{bmatrix}
-1 & 1 \\
-1 & -1
\end{bmatrix} \vec{u}.
$$




The characteristic equation is:




$$
\begin{vmatrix}
-1 - \lambda & 1 \\
-1 & -1 - \lambda
\end{vmatrix} = (-1 - \lambda)^2 - (-1)(1) = (1 + \lambda)^2 - 1 = \lambda^2 + 2\lambda + 2 = 0,
$$




which yields complex eigenvalues:




$$
\lambda_1 = -1 + i, \quad \lambda_2 = -1 - i.
$$




Find the eigenvectors:




For $\lambda_1 = -1 + i$:




$$
\left(A - \lambda_1 I\right) =
\begin{bmatrix}
-1 - (-1 + i) & 1 \\
-1 & -1 - (-1 + i)
\end{bmatrix} =
\begin{bmatrix}
- i & 1 \\
-1 & - i
\end{bmatrix}.
$$




Solving:




$$
\begin{aligned}
-ix + y = 0 \implies y = ix,
\end{aligned}
$$




gives:




$$
\vec{x}_1 = \begin{bmatrix} 1 \\ i \end{bmatrix}.
$$




Similarly, for $\lambda_2 = -1 - i$:




$$
\left(A - \lambda_2 I\right) =
\begin{bmatrix}
-1 - (-1 - i) & 1 \\
-1 & -1 - (-1 - i)
\end{bmatrix} =
\begin{bmatrix}
i & 1 \\
-1 & i
\end{bmatrix},
$$




which yields:




$$
\vec{x}_2 = \begin{bmatrix} 1 \\ -i \end{bmatrix}.
$$




Thus, the general solution is:




$$
\vec{u}(t) = e^{-t} \left( c_1 \vec{x}_1 e^{i t} + c_2 \vec{x}_2 e^{-i t} \right),
$$




which represents a damped oscillation with angular frequency $\omega = 1$ and exponential decay rate $\sigma = -1$.




Note that the system matrix can be decomposed as:




$$
\begin{bmatrix}
-1 & 1 \\
-1 & -1
\end{bmatrix} =
\begin{bmatrix}
-1 & 0 \\
0 & -1
\end{bmatrix} +
\begin{bmatrix}
0 & 1 \\
-1 & 0
\end{bmatrix},
$$




where the first matrix represents exponential decay and the second represents rotation.














In [ ]:






%plot -s 900,350

A = [-1 1;
     -1 -1];

%% eigen-analysis
[S,D] = eig(A);

lamb = diag(D)

x1 = S(:,1);
x2 = S(:,2);

u0 = [0;6];
C = inv(S)*u0;

t = 0:0.01:10;
u = C(1)*x1*exp(lamb(1)*t) + C(2)*x2*exp(lamb(2)*t);

% plot u1 and u2 as a function of time
plot(t,u(1,:),t,u(2,:),t,zeros(size(t)),'k--')
xlabel('time')
legend('u_1', 'u_2')



















lamb =

  -1.0000 + 1.0000i
  -1.0000 - 1.0000i














No description has been provided for this image















In [ ]:






%plot -s 900,420

subplot(1,2,1)
% plot a trajectory of u1 and u2
for i = 1:length(t)
    plot(u(1,i), u(2,i), 'k.'), hold on
end
hold off

axis equal
axis([-12 12 -12 12])
grid on
xlabel('u_1', 'fontsize', 12)
ylabel('u_2', 'fontsize', 12)
title('Trajectory', 'fontsize', 12)

subplot(1,2,2)
plot(real(lamb),imag(lamb),'o')
axis equal
grid on
axis([-4,4,-6,6])
xlabel('Re', 'fontsize', 12)
ylabel('Im', 'fontsize', 12)
title('Eigenvalues','fontsize',12)






























No description has been provided for this image





















6.2. The Second-Order ODE



This section follows the approach from MIT 2.087 Engineering Mathematics: Linear Algebra and ODEs, Fall 2014, by Gilbert Strang.














In [ ]:






from IPython.display import YouTubeVideo
YouTubeVideo('xvTYUnqn2wY', width = "560", height = "315")















Out[ ]:























6.2.1. Spring-Mass-Damper System
Consider the generic second-order linear differential equation:




$$
a \frac{d^2x(t)}{dt^2} + b \frac{dx(t)}{dt} + c x(t) = 0, \quad x(0) = x_0, \quad \dot{x}(0) = v_0.
$$




A classic example of this is the mass-spring-damper system.
































This equation describes free vibration of a simple mass-spring-damper system, particularly for non-zero initial conditions. The equation of motion is:




$$
\begin{aligned}
-kx - c\dot{x} &= m\ddot{x}, \\
m\ddot{x} + c\dot{x} + kx &= 0, \\
\ddot{x} + \frac{c}{m}\dot{x} + \frac{k}{m}x &= 0.
\end{aligned}
$$




Rewriting this in terms of the damping ratio $\zeta$ and natural frequency $\omega_n$:




$$
\ddot{x} + 2\zeta \omega_n \dot{x} + \omega_n^2 x = 0,
$$




where:




$$
\begin{aligned}
\omega_n^2 &= \frac{k}{m}, \\
2\zeta \omega_n &= \frac{c}{m}, \\
\zeta &= \frac{c}{2\sqrt{km}} = \frac{1}{2}\frac{c}{\sqrt{km}}.
\end{aligned}
$$




The underdamped solution is:




$$
x(t) = e^{-\zeta \omega_n t} \left(a_1 e^{i \omega_d t} + a_2 e^{-i \omega_d t} \right),
$$




or equivalently:




$$
x(t) = e^{-\zeta \omega_n t} \left(b_1 \cos \omega_d t + b_2 \sin \omega_d t \right),
$$




where $\omega_d = \omega_n \sqrt{1 - \zeta^2}$.




Using initial conditions $x(0) = x_0$, $\dot{x}(0) = v_0$, we find:




$$
x(t) = e^{-\zeta \omega_n t} \left(x_0 \cos \omega_d t + \frac{\zeta \omega_n x_0 + v_0}{\omega_d} \sin \omega_d t \right).
$$
















Experimental Demo for Dampled Oscillating












In [ ]:






from IPython.display import YouTubeVideo
YouTubeVideo('ZqedDWEAUN4?start=80&end=114', width = "560", height = "315")















Out[ ]:



















In [ ]:






from IPython.display import YouTubeVideo
YouTubeVideo('pjzggPk8DO4', width = "560", height = "315")















Out[ ]:























6.2.2. State-Space Representation
Define:




$$
\vec{x} =
\begin{bmatrix}
x_1 \\
x_2
\end{bmatrix}
=
\begin{bmatrix}
x \\
\dot{x}
\end{bmatrix}.
$$




Then the system becomes:




$$
\begin{align*}
\dot{\vec{x}} &= \left[ {\begin{matrix}
   \dot{x}_1 \\
   \dot{x}_2 \\   
  \end{matrix} } \right]
  =
\left[ {\begin{matrix}
   0 & 1 \\
   -\omega_n^2 & -2\zeta \omega_n \\   
  \end{matrix} } \right]
\left[ {\begin{array}{cc}
   x_1 \\
   x_2 \\   
\end{array} } \right] = \left[ {\begin{matrix}
   0 & 1 \\
   -\omega_n^2 & -2\zeta \omega_n \\   
  \end{matrix} } \right] \vec{x}\\ \\y & =
\left[ {\begin{matrix}
   1 & 0    
\end{matrix} } \right]  
\left[ {\begin{array}{cc}
   x_1 \\
   x_2 \\   
\end{array} } \right]  
\end{align*}
$$




Eigenanalysis


Compute the eigenvalues of:




$$
A =
\begin{bmatrix}
0 & 1 \\
-\omega_n^2 & -2\zeta \omega_n
\end{bmatrix}.
$$




The characteristic equation is:




$$
\begin{aligned}
\det(A - \lambda I) &= \lambda^2 + 2\zeta \omega_n \lambda + \omega_n^2 = 0, \\
\lambda &= -\zeta \omega_n \pm \omega_n \sqrt{\zeta^2 - 1}.
\end{aligned}
$$




The behavior depends on $\zeta$:




$$
\begin{aligned}
\zeta > 1 &\implies \lambda \text{ real, negative (overdamped)}, \\
0 < \zeta < 1 &\implies \lambda \text{ complex conjugate} \implies \text{oscillatory decay (underdamped)}.
\end{aligned}
$$






Interpretation of $0 < \zeta < 1$




$$
\begin{align*}
\lambda &= -\zeta\omega_n \pm j\omega_n\sqrt{1-\zeta^2} \\\\
e^{\lambda t}&= \underbrace{e^{-\zeta\omega_n t}}_{\text{decaying}} \cdot \underbrace{e^{j\omega_n\sqrt{1-\zeta^2}t}}_{\text{oscillating}}
\end{align*}
$$




showing the system response decays exponentially while oscillating.






Eigenvalues in the $s$-plane




$$
\lambda = -\zeta \omega_n \pm j \omega_n \sqrt{1 - \zeta^2}.
$$




Define:




$$
\begin{align*}
\zeta \omega_n &= \sigma\\\\
\omega_n \sqrt{1-\zeta^2} &= \omega_d \\\\
\zeta &= \frac{\sigma}{\omega_n} = \sin \theta \\\\
\frac{\omega_d}{\omega_n} &= \cos \theta = \sqrt{1-\sin^2 \theta} = \sqrt{1-\zeta^2}
\end{align*}
$$




This captures the physical interpretation: oscillations with damping in the underdamped regime.
































Pure oscillating


    

  • Corresponds to the case when the damping ratio $\zeta=0$. This is the undamped case, representing pure oscillation with eigenvalues located purely on the imaginary axis.
    • 

































Critical damping


    

  • Illustrates the special case when $\zeta = 1$, leading to real and repeated eigenvalues.
    

    • 

  • No oscillation: the response is non-oscillatory.
    

    • 

  • Fastest settling: the system returns to equilibrium more quickly than any overdamped system but without overshooting.
    

    • 

































Over damping


    

  • No oscillations: the system response is purely exponential decay.
    

    • 

  • Slow return: slower than the critically damped case; the system takes longer to return to equilibrium.
    

    • 

  • Distinct real poles: both eigenvalues are on the real axis, both negative.
    

    • 



































Let's simulate the response of a second-order system with natural frequency $\omega_n = 2$ and damping ratio $\zeta = 0.1$, which corresponds to the underdamped case.












In [ ]:






%plot -s 900,350

wn = 2; zeta = 0.1;
A = [0 1;-wn^2 -2*zeta*wn];

[S,D] = eig(A);
lamb = diag(D)

t = 0:0.01:10;
U0 = [2 10]';
C = inv(S)*U0;

u = C(1)*S(:,1)*exp(lamb(1)*t) + C(2)*S(:,2)*exp(lamb(2)*t);

% plot u1 and u2 as a function of time
plot(t,u(1,:),t,u(2,:),t,zeros(size(t)),'k--')
axis tight
legend('u_1','u_2')



















lamb =

  -0.2000 + 1.9900i
  -0.2000 - 1.9900i














No description has been provided for this image















In [ ]:






%plot -s 900,420

subplot(1,2,1)
% plot a trajectory of u1 and u2
for i = 1:length(t)
    plot(u(1,i), u(2,i), 'k.'), hold on
end
hold off

axis equal
axis([-12 12 -12 12])
grid on
xlabel('u_1', 'fontsize', 12)
ylabel('u_2', 'fontsize', 12)
title('Trajectory', 'fontsize', 12)

subplot(1,2,2)
plot(real(lamb),imag(lamb),'o')
axis equal
grid on
axis([-4,4,-6,6])
xlabel('Re', 'fontsize', 12)
ylabel('Im', 'fontsize', 12)
title('Eigenvalues','fontsize',12)






























No description has been provided for this image






















7. Stability



When analyzing dynamic systems, one of the fundamental questions is: Will the system “blow up,” or will it settle down?

In control design, stability is a primary objective, alongside:


    

  • Accurate tracking of desired signals
    • 

  • Robustness to model errors and disturbances
    • 





It’s often instructive to build intuition starting from scalar systems before moving to higher-dimensional systems.




7.1. Scalar System Example
Consider a simple linear scalar differential equation:




$$
\dot{x} = a x \quad \implies \quad x(t) = e^{a t} x(0)
$$




Here, the behavior of the system depends entirely on the sign of the parameter $a$:




$$\begin{cases}
a>0 & : \quad \text{unstable} \ (\text{exponential growth}) \\
a<0 & : \quad \text{asymptotically stable} \ (\text{exponential decay}) \\
a=0 & : \quad \text{critically stable} \ (\text{constant, marginally stable})
\end{cases}$$




7.2. From Scalars to Matrices: General Linear Systems
For vector systems:




$$
\dot{x} = A x
$$




we can no longer simply look at the entries of $A$ to determine stability. Unlike the scalar case where we can say $a>0$, matrices require a different tool: eigenvalues.




An eigenvalue-eigenvector pair $(\lambda, v)$ of $A$ satisfies:




$$
A v = \lambda v
$$




The eigenvalues $\lambda$ reveal how the matrix “acts” along different directions (given by eigenvectors $v$). Specifically:




    

  • If $\operatorname{Re}(\lambda) > 0$, solutions grow exponentially in the direction of $v$ (instability).
    • 

  • If $\operatorname{Re}(\lambda) < 0$, solutions decay exponentially (asymptotic stability).
    • 

  • If $\operatorname{Re}(\lambda) = 0$, solutions neither grow nor decay, which may indicate marginal stability.
    • 





$$\begin{cases}
\operatorname{Re}(\lambda) > 0 & : \quad \text{unstable} \\
\operatorname{Re}(\lambda) < 0 & : \quad \text{asymptotically stable} \\
\operatorname{Re}(\lambda) = 0 & : \quad \text{critically stable}
\end{cases}$$




7.3. Summary
    

  • Scalar systems provide a clear picture: stability depends on the sign of $a$.
    • 

  • Matrix systems generalize this: stability depends on the real parts of the eigenvalues of $A$.
    • 

  • Stability is a cornerstone for designing controllers that ensure performance and robustness.
    • 















In [1]:






%%javascript
$.getScript('https://kmahelona.github.io/ipython_notebook_goodies/ipython_notebook_toc.js')