Dynamic Systems:

Natural Response

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

# 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$
InÂ [1]:
%%html
width="560" height="315" frameborder="0" allowfullscreen></iframe></center>

InÂ [2]:
% 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*}
InÂ [3]:
%%html
width="560" height="315" frameborder="0" allowfullscreen></iframe></center>

InÂ [4]:
%%html
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}$$
InÂ [5]:
%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

InÂ [6]:
%%html
width="420" height="315" frameborder="0" allowfullscreen></iframe></center>

InÂ [7]:
%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*}
InÂ [8]:
%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')


InÂ [9]:
%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*}
InÂ [10]:
%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')


InÂ [11]:
%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


# 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
InÂ [12]:
%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')


InÂ [13]:
%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


InÂ [14]:
%plot -s 900,350

v = inv(S)*u;

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


InÂ [15]:
%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 (?)}
InÂ [16]:
%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


InÂ [17]:
%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*}
InÂ [18]:
%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


InÂ [19]:
%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*}
InÂ [20]:
%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


InÂ [21]:
%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}$$

InÂ [22]:
%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


InÂ [23]:
%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
InÂ [24]:
%%html
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

InÂ [25]:
%%html
width="560" height="315" frameborder="0" allowfullscreen></iframe></center>

InÂ [26]:
%%html
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