Dynamic Systems:

Natural Response

By Prof. Seungchul Lee
http://iai.postech.ac.kr/
Industrial AI Lab at POSTECH

Table of Contents

1. Natural Response to Non-zero Initial Conditions

1.1. The First Order ODE

• MIT 2.087 Engineering Mathematics: Linear Algebra and ODEs, Fall 2014, by Gilbert Strang

$$ \begin{align*} \frac{dx(t)}{dt} &= kx(t), \qquad x(0)=x_0 \\ \\ \rightarrow \; x(t) &= x_0e^{\,kt} \end{align*} $$

• Solution will be exponential functions
  • unknown coefficient determined by I.C.

• Stability
  • unstable if $k > 0$
  • stable if $k < 0$

%%html
<center><iframe src="https://www.youtube.com/embed/4X0SGGrXDiI?start=399&end=459?rel=0" 
width="560" height="315" frameborder="0" allowfullscreen></iframe></center>

% 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)

• $\tau$: time constant
  • large $\tau$: slow response
  • small $\tau$: fast response

$$ \begin{align*} \dot{x} + \frac{1}{\tau}x &= 0 \quad \implies \quad \dot{x} = - \frac{1}{\tau} x = kx\\ \\ x(t) &= x(0) e ^{-\frac{1}{\tau}t}\\ \\ \frac{x(\tau)}{x(0)} &= e^{-1} = \frac{1}{e} = 0.368\cdots \end{align*} $$

1.2. Two First Order ODEs (Independent)

$$ \begin{align*} \dot{u_1} &= \lambda_1 u_1 \implies u_1(t) = u_1(0) e^{\lambda_1 t}\\ \\ \dot{u_2} &= \lambda_2 u_2 \implies u_2(t) = u_2(0) e^{\lambda_2 t} \end{align*} $$

• Suppose $u_1$ and $u_2$ are independent

$$ \begin{align*} u &= \begin{bmatrix}u_1 \\ u_2 \end{bmatrix} \\ \\ \dot{u} &= \begin{bmatrix}\dot{u}_1 \\ \dot{u}_2 \end{bmatrix} = \begin{bmatrix}\lambda_1 u_1 \\ \lambda_2 u_2 \end{bmatrix} = \begin{bmatrix}\lambda_1 & 0 \\ 0 & \lambda_2 \end{bmatrix} \begin{bmatrix}u_1 \\ u_2 \end{bmatrix} = \Lambda u \end{align*} $$

1.3. ODE in Vector Form (dependent)

• Suppose $u_1$ and $u_2$ are dependent

$$ \begin{align*} \dot{u_1} &= a_{11} u_1 + a_{12} u_2\\ \dot{u_2} &= a_{21} u_1 + a_{22} u_2 \\\\ \dot{u} &= \begin{bmatrix}a_{11} & a_{12} \\ a_{21} & a_{22} \end{bmatrix} u = A u \\\\ u(t) &= \begin{bmatrix}u_1(t) \\ u_2(t) \end{bmatrix} = ? \end{align*} $$

2. Systems of Differential Equations

• Matrix Methods | MIT 18.03SC Differential Equations, Fall 2011

• Given

$\quad\;$ (I will use mixed notations, $\vec{u}$ or $u$)

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

• (Supperposition) If $\vec{u}_1$ and $\vec{u}_2$ are solutions of $\dot{\vec{u}} = A\vec{u}$, so is any linear combination $\vec{u} = c_1\vec{u}_1 + c_2\vec{u}_2$

$$ \begin{align*} \dot{\vec{u}}_1 &= A\vec{u}_1 \\ \dot{\vec{u}}_2 &= A\vec{u}_2 \\ \\ \dot{\vec{u}} &= c_1\dot{\vec{u}}_1 + c_2\dot{\vec{u}}_2 = c_1A\vec{u}_1 + c_2A\vec{u}_2 \\ &= A(c_1\vec{u}_1 + c_2\vec{u}_2) \\ &= A\vec{u} \end{align*} $$

• For a single and scalar ODE

$$\dot{u} = au \quad \implies \quad u(t) = ce^{at}$$

• Let us try

$$\vec{u}(t) = \vec{x}e^{\lambda t} \quad \text{(in a vector form)} $$

$$ \begin{align*} \dot{\vec{u}} &= \vec{x}\lambda e^{\lambda t} =\underline{\lambda\vec{x}} e^{\lambda t} \\ A\vec{u} &= \underline{A\vec{x}} e^{\lambda t} \qquad \qquad \Longleftrightarrow \quad A\vec{x} = \lambda\vec{x} \quad \text{(eigenvalue problem)} \end{align*}$$

• Linear ODE $=$ Eigenvalue problem

2.1. Eigenanalysis

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

• Eigenanalysis

$$ A\vec{x}_1 = \lambda_1 \vec{x}_1$$$$ A\vec{x}_2 = \lambda_2 \vec{x}_2$$

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

• $\vec{u}(0)$ can be decomposed to $c_1$ and $c_2$ with respect to $\vec{x}_1$ and $\vec{x}_2$

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

$$ \begin{align*} c_1\vec{x}_1 \quad &\Longrightarrow \quad c_1\vec{x}_1 e^{\lambda_1 t} \\ c_1\vec{x}_2 \quad &\Longrightarrow \quad c_2\vec{x}_2 e^{\lambda_2 t} \\\\ \vec{u}(0) = c_1\vec{x}_1 + c_2\vec{x}_2 \quad &\Longrightarrow \quad \vec{u}(t) = c_1\vec{x}_1 e^{\lambda_1 t} + c_2\vec{x}_2 e^{\lambda_2 t} \\ \end{align*} $$

2.2. General Solutions for ODE

• For $n \times n$ matrix $A$

$$ \begin{align*} \lambda_i &: \;\; i^{\text{th}} \text{ eigenvalues of } A \\ \vec{x}_i &: \;\; i^{\text{th}} \text{ eigenvectors of } A \\ \\ \vec{u}_i &= \vec{x}_i e^{\lambda_1 t} \; \text{ solution of } \; \dot{\vec{u}} = A\vec{u} \end{align*} $$

• Then, a general solution

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

• $c_1, c_2, \cdots, c_n$ are determined by initial conditions

• Systems of differential equations $\quad \leftrightarrow \quad $ Eigenanalysis (Gilbert Strang Lecture 21, video below)

• System stability $\quad \leftrightarrow \quad$ eigenvalues

2.3. ODE and Linear Algebra

$$\dot{x}(t) = Ax(t)$$

$$ \begin{align*} Ax_1 &= \lambda_1x_1, \qquad x_1, x_2: \quad \text{eigenvectors} \\ Ax_2 &= \lambda_2x_2, \qquad \lambda_1, \lambda_2: \quad \text{eigenvalues} \end{align*} $$

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

$$ \begin{align*} \dot{c}_1(t)x_1 &= Ac_1(t)x_1 \\ &= c_1(t)Ax_1 = c_1(t)\lambda_1x_1 \\\\ &\Downarrow\\\\ \dot{c}_1(t) &= \lambda_1c_1(t) \\ \therefore \; c_1(t) &= c_1(0)e^{\lambda_1t} \end{align*} $$

• Suppose $x(t) = c_2(t)x_2$

$$\therefore \quad c_2(t) = c_2(0)e^{\lambda_2t}$$

$$ \begin{align*} A \begin{bmatrix} x_1 & x_2 \end{bmatrix} = \begin{bmatrix} \lambda_1 x_1 & \lambda_2 x_2 \end{bmatrix} &= \begin{bmatrix} x_1 & x_2 \end{bmatrix} \begin{bmatrix} \lambda_1 & 0 \\ 0 & \lambda_2 \end{bmatrix}\\ \\ \\ S &\triangleq \begin{bmatrix} x_1 & x_2 \end{bmatrix} \quad \text{eigenvector matrix}\\\\ \Lambda &= \begin{bmatrix} \lambda_1 & 0 \\ 0 & \lambda_2 \end{bmatrix} \quad \text{diagnonal matrix} \end{align*} $$

$$ \begin{align*} AS &= S\Lambda \\ S^{-1}AS& = \Lambda \\ A &= S\Lambda S^{-1} \end{align*} $$

• Linear Transformation

$$ \begin{align*} u &= Sv \\ \\ \dot{u} &= Au \\ S\dot{v} &= ASv \\ \\ \dot{v} &= S^{-1}ASv \\ &= \Lambda v, \quad \qquad v(0) = S^{-1}u(0) \end{align*} $$

• Decoupled in $v$-frame

• Solution

$$u(t) = v_1(0)e^{\lambda_1t}x_1 + v_2(0)e^{\lambda_2t}x_2$$

$\quad\;$ I.C.

$$ \begin{align*} u(0) &= v_1(0)x_1 + v_2(0)x_2 \\ &= \begin{bmatrix} x_1 & x_2 \end{bmatrix} \begin{bmatrix} v_1(0) \\ v_2(0) \end{bmatrix} = Sv(0) \end{align*} $$

• $v$-frame is decoupled in the coordinates of $\vec{x}_1$ and $\vec{x}_2$

$$ v(t) = \begin{bmatrix} v_1(0)e^{\lambda_1t} \\ v_2(0)e^{\lambda_2t} \end{bmatrix} = \begin{bmatrix} e^{\lambda_1t} & 0 \\ 0 & e^{\lambda_2t} \end{bmatrix} \begin{bmatrix} v_1(0) \\ v_2(0) \end{bmatrix}, \quad \quad v(0) = \begin{bmatrix} v_1(0) \\ v_2(0) \end{bmatrix} $$

$$ \begin{align*} u(t) &= \begin{bmatrix} x_1 & x_2 \end{bmatrix} \begin{bmatrix} e^{\lambda_1t} & 0 \\ 0 & e^{\lambda_2t} \end{bmatrix} v(0) \\ &= \begin{bmatrix} x_1 & x_2 \end{bmatrix} \begin{bmatrix} e^{\lambda_1t} & 0 \\ 0 & e^{\lambda_2t} \end{bmatrix} \begin{bmatrix} x_1 & x_2 \end{bmatrix}^{-1} u(0) \\ \\ &= S e^{\Lambda t} S^{-1} u(0) \end{align*} $$

%%html
<center><iframe src="https://www.youtube.com/embed/YUjdyKhWt6E?rel=0" 
width="560" height="315" frameborder="0" allowfullscreen></iframe></center>

%%html
<center><iframe src="https://www.youtube.com/embed/lXNXrLcoerU?rel=0" 
width="420" height="315" frameborder="0" allowfullscreen></iframe></center>

3. Real Eigenvalues

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

$\quad \;$ where $$ \left[ {\begin{matrix} c_1 \\ c_2 \\ \end{matrix} } \right] = \left[ {\begin{array}{cc} \vec{x}_1 \, \, \vec{x}_2 \\ \end{array} } \right]^{-1} \left[ {\begin{array}{cc} u_1(0) \\ u_2(0) \\ \end{array} } \right]. $$

3.1. Example 1

$$ \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} $$

$\quad\;$ 1) Find the 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*} $$

$\quad\;$ 2) Find the corresponding 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*} $$

$\quad\;$ 3) Find $\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*} $$

$\quad\;$ 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 (mode, coordinate, axsis)

$$ \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*} $$

• Simulation

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

%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')

Phase portrait

• https://en.wikipedia.org/wiki/Phase_portrait
• http://tutorial.math.lamar.edu/Classes/DE/PhasePlane.aspx
• Lec 27 | MIT 18.03 Differential Equations, Spring 2006

%%html
<center><iframe src="https://www.youtube.com/embed/e3FfmXtkppM?rel=0" 
width="420" height="315" frameborder="0" allowfullscreen></iframe></center>

%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

3.2. Example 2

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

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

$$ \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*} $$

%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')

%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

3.3. Example 3

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

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

$$ \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*} $$

%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')

%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

3.4. Different Eigenvectors with the Same Eigenvalues

4. De-coupling via Linear Transformation

- ODE $$ \begin{align*} \,\dot{u}_1 &= -2(u_1 - u_2) \\ \,\dot{u}_2 &=\;\;\;\,\,(u_1 - u_2) \end{align*} $$
$$ \dot{\vec{u}} = \begin{bmatrix} -2 & 2 \\ 1 & -1 \end{bmatrix} \vec{u} $$
$$ \begin{align*} \begin{vmatrix} -2 -\lambda & 2 \\ 1 & -1 - \lambda \end{vmatrix} &= (2 + \lambda)(1 + \lambda) - 2 = \lambda^2 + 3\lambda = \lambda (\lambda + 3) = 0 \end{align*} $$
$$ \begin{align*} \lambda_1 = 0 \quad \implies \quad \begin{bmatrix} -2 & 2 \\ 1 & -1 \end{bmatrix} \begin{bmatrix} x \\ y \end{bmatrix} &= \begin{bmatrix} 0 \\ 0 \end{bmatrix} \quad \implies \quad \vec{x}_1 = \begin{bmatrix} 1 \\ 1 \end{bmatrix} \\ \lambda_2 = -3 \quad \implies \quad \begin{bmatrix} 1 & 2 \\ 1 & 2 \end{bmatrix} \begin{bmatrix} x \\ y \end{bmatrix} &= \begin{bmatrix} 0 \\ 0 \end{bmatrix} \quad \implies \quad \vec{x}_2 = \begin{bmatrix} 2 \\ -1 \end{bmatrix} \end{align*} $$
- Change variables $\lambda_1 = 0 \quad \implies \quad e^{0t} = 1 \quad \implies \quad$ does not change over time (= invariant)

$\quad\;$ 1) Total amount of water

$$v_1 = 1u_1 + 2u_2$$

$\quad \;$ 2) Difference in height

$$v_2 = u_1 - u_2$$

• De-coupled

$$ \begin{align*} \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{align*} $$

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

• Trajectory Comparison

%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')

%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

%plot -s 900,350

v = inv(S)*u;

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

%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)

5. Complex Eigenvalues (Starting Oscillation)

• $\lambda$ can be a complex number $\lambda = \sigma + i \omega$

$$ \begin{align*} e^{\lambda t} = e^{(\sigma + i\omega)t} &= \underbrace{e^{\sigma t}}_{\text{decay}}\; \underbrace{e^{i\omega t}}_{\text{oscillates}} \\\\ &= e^{\sigma t}(\cos \omega t + i\sin \omega t) \end{align*} $$

$$ \begin{align*} \text{If Re}(\lambda) = \sigma < 0 &\quad \implies \quad \text{stable} \\ &\quad \implies \quad \text{decay by } e^{\sigma t} \end{align*} $$

5.1. Example 1

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

$$ \begin{vmatrix} -\lambda & 1 \\ -1 & -\lambda \end{vmatrix} = \lambda^2 + 1 = 0, \qquad \begin{align*} &\lambda_1 = i \\ &\lambda_2 = -i \end{align*} \qquad \text{complex conjugate} $$

$$ \begin{align*} \begin{bmatrix} -i & 1 \\ -1 & -i \end{bmatrix} \begin{bmatrix} x \\ y \end{bmatrix} \begin{bmatrix} 0 \\ 0 \end{bmatrix} \quad &\implies \quad \vec{x}_1 = \begin{bmatrix} 1 \\ i \end{bmatrix} \\ \begin{bmatrix} -i & 1 \\ -1 & i \end{bmatrix} \begin{bmatrix} x \\ y \end{bmatrix} \begin{bmatrix} 0 \\ 0 \end{bmatrix} \quad &\implies \quad \vec{x}_2 = \begin{bmatrix} 1 \\ -i \end{bmatrix} \end{align*} \qquad \text{ vector complex conjugate} $$

• Note

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

• Solution

$$ \begin{align*} \vec{u}(t) &= c_1\vec{x}_1e^{\lambda_1t} + c_2\vec{x}_2e^{\lambda_2t} \\\\ &= c_1 \begin{bmatrix} 1 \\ i \end{bmatrix} e^{it} + c_2 \begin{bmatrix} 1 \\ -i \end{bmatrix} e^{-it} \\ \\ &= \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} \\ \\ \vec{u}(0) &= \begin{bmatrix} c_1 + c_2 \\ ic_1 - ic_2 \end{bmatrix} = \begin{bmatrix} 1 \\ 0 \end{bmatrix} \quad \implies \quad c_1 = c_2 = {1 \over 2} \\ \\ \vec{u}(t) &= \begin{bmatrix} {1 \over 2} (e^{it} + e^{-it}) \\ -{1 \over 2i} (e^{it} - e^{-it}) \end{bmatrix} = \begin{bmatrix} \cos t \\ -\sin t \end{bmatrix} \end{align*} $$

• What is the corresponding physical system?

  • Simple harmonic motion Revisited

$$ \ddot{u} + \omega^2_nu = 0, \qquad \text{assume } \omega_n = 1 $$

$$ \begin{align*} \ddot{u} + u = 0, \qquad &u_1 = u \text{ : displacement or position} \\ &u_2 = \dot{u} \text{ : velocity} \end{align*} $$

$$ \begin{align*} \begin{bmatrix} \dot{u}_1 \\ \dot{u}_2 \end{bmatrix} &= \begin{bmatrix} \dot{u} \\ \ddot{u} \end{bmatrix} = \begin{bmatrix} u_2 \\ -u_1 \end{bmatrix} \\ \\ \therefore \; \dot{\vec{u}} &= \begin{bmatrix} 0 & 1 \\ -1 & 0 \end{bmatrix} \vec{u} \end{align*} $$

$$ \begin{vmatrix} -\lambda & 1 \\ -1 & -\lambda \end{vmatrix} = \lambda^2 + 1 = 0, \qquad \begin{align*} &\lambda_1 = i \\ &\lambda_2 = -i \end{align*} \qquad \text{two angular velocities (?)} $$

%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

%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)

5.2. Example 2

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

$$ \begin{vmatrix} -\lambda & 1 \\ -4 & \lambda \end{vmatrix} = \lambda^2 + 4 = 0, \quad \begin{matrix} \begin{align*} \lambda_1 &= 2i \\ \lambda_2 &= -2i \end{align*} \end{matrix} \quad \implies \quad \text{angular velocity} $$

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

$$ \begin{align*} \vec{u}_t &= c_1\vec{x}_1e^{\lambda_1t} + c_2\vec{x}_2e^{\lambda_2t} \\ \\ &= c_1 \begin{bmatrix} 1 \\ 2i \end{bmatrix} e^{i2t} + c_2 \begin{bmatrix} 1 \\ -2i \end{bmatrix} e^{-i2t} \\ \\ \vec{u}(t) &= \begin{bmatrix} {1 \over 2}\left(e^{i2t} + e^{-i2t}\right) \\ -2{1 \over 2i}\left(e^{i2t} - e^{-i2t}\right) \end{bmatrix} = \begin{bmatrix} \cos t \\ -2\sin t \end{bmatrix} \\ \\ \vec{u}(0) &= \begin{bmatrix} c_1 + c_2 \\ 2i(c_1 - c_2) \end{bmatrix} = \begin{bmatrix} 1 \\ 0 \end{bmatrix} \quad \implies \quad c_1 = c_2 = {1 \over 2} \end{align*} $$

%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

%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)

5.3. Example 3

$$ \begin{align*} \dot{u} = \begin{bmatrix} 3 & -9 \\ 4 & -3 \end{bmatrix} u \end{align*} $$

$$ \begin{vmatrix} 3 - \lambda & -9 \\ 4 & -3 - \lambda \end{vmatrix} = \lambda^2 - 9 + 36 = \lambda^2 + 27 = 0, \qquad \begin{align*} \lambda_1 &= 3\sqrt{3}i \\ \lambda_2 &= -3\sqrt{3}i \end{align*} $$

%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

%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)

6. Complex Eigenvalues with Damping

6.1. Example 1

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

$$ \begin{vmatrix} -1-\lambda & 1 \\ -1 & -1 - \lambda \end{vmatrix} = (1+\lambda)^2 - 1 = \lambda^2 + 2\lambda + 2 = 0, \qquad \begin{align*} \lambda_1 = -1 + i \\ \lambda_2 = -1 - i \end{align*} $$

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

$$ u(t) = e^{-t}(c_1\vec{x}_1e^{it} + c_2\vec{x}_2e^{-it}) $$

• Note

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

%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

%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)

6.2. The Second Order ODE

• MIT 2.087 Engineering Mathematics: Linear Algebra and ODEs, Fall 2014, by Gilbert Strang

%%html
<center><iframe src="https://www.youtube.com/embed/xvTYUnqn2wY?rel=0" 
width="560" height="315" frameborder="0" allowfullscreen></iframe></center>

Generic form

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

One of examples is the mass-spring-damper system

This notebook simluates the free vibration of a simple mass-spring-damper system. More specifically, we'll look at how system response to non-zero initial conditions.

The equation of motion for the system is:

$$ \begin{align*} -kx - c\dot{x} &= m\ddot{x} \\ \\ m\ddot{x} + c\dot{x} + kx& = 0\\\\ \ddot{x} + {c \over m}\dot{x} + {k \over m}x &= 0 \end{align*} $$

We could also write this equation in terms of the damping ratio, $\zeta$, and natural frequency $\omega_n$.

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

where

$$ \begin{align*} \omega^2_n &= {k \over m}\\ 2\zeta\sqrt{k \over m} = {c \over m} \quad \implies \quad \zeta &= {1 \over 2} {c \over m} \sqrt{m \over k} = {1 \over 2}{c \over \sqrt{mk}} \end{align*} $$

We'll use the solution to the differential equation that we developed in class to plot the response. The solution for the underdamped case is:

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

or

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

To use this equation, we need to solve for $a_1$ and $a_2$ or $b_1$ and $b_2$ using the initial conditions. Here, let's use the sin/cosine form. Solving the equation for generic intial velocity, $\dot{x} = v_0$, and a generic initial displacement, $x = x_0$, we find:

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

Experiment

%%html
<center><iframe src="https://www.youtube.com/embed/ZqedDWEAUN4?start=80&end=114?rel=0" 
width="560" height="315" frameborder="0" allowfullscreen></iframe></center>

%%html
<center><iframe src="https://www.youtube.com/embed/pjzggPk8DO4?rel=0" 
width="560" height="315" frameborder="0" allowfullscreen></iframe></center>

• Define states

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

• State space representation

$$ \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*} $$

$$\dot{\vec{x}} = A\vec{x}, \qquad A = \begin{bmatrix} 0 & 1 \\ -\omega^2_n & -2\zeta\omega_n \end{bmatrix}$$

• Eigenanalysis

$$ \begin{align*} \begin{vmatrix} -\lambda & 1 \\ -\omega^2_n & -2\zeta\omega_n - \lambda \end{vmatrix} &= \lambda(2\zeta\omega_n + \lambda) + \omega^2_n \\ &= \lambda^2 + 2\zeta\omega_n\lambda + \omega^2_n = 0 \end{align*} $$

$$ \begin{align*} \therefore \lambda &= -\zeta\omega_n \pm \sqrt{\zeta^2\omega^2_n - \omega^2_n} \\ &= -\zeta \omega_n \pm \omega_n\sqrt{\zeta^2-1} \end{align*} $$

$$ \begin{align*} &\text{if} \quad \zeta > 1, \qquad \;\; \lambda < 0 \text{ and real} \\ &\text{if} \quad 0 < \zeta < 1, \quad \lambda \;\text{ complex number} \implies \text{start oscilating} \end{align*} $$

• Physical 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*} $$

• Eigenvalues in s-plane

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

$$ \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*} $$

• Oscillating with damping (under damping)

• Pure oscillating

• Critical damping

• Over damping

• Simulation with $\omega_n = 2$, and $\zeta = 0.1$

%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

%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)

7. Stability

Always trying to figure out if the system "blows up" or not

Recall the control design objectives:

• stability
• tracking
• robustness

It is useful to start with scalar systems to get some intuition about what is going on

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

$\begin{cases} a>0: \text{unstable}\\ a<0: \text{asymptotically stable}\\ a=0: \text{critically stable} \end{cases}$

From scalars to mattices?

We cannot say that $A>0$, but we can do the next best thing - eigenvalues !

$$Av = \lambda v$$

The eigenvalues tell us how the matrix $A$ 'acts' in different directions (eigenvectors)

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

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