[고등자동제어] Similarity Transformation for matrix exponential (2)

 

(1)에서 이어지는 내용이기 때문에 (1)먼저 보는 것을 추천드리지만, (2)만 보는 분들을 위해

matrix A에 대해 eigenvalues, eigenvectors, 그리고 transformation matrix T까지 구하는 과정들을 첨부합니다. (읽으셔야 아래 내용이 더 이해가 잘 됩니다.)

 

 

4) Meaning of eigenvalues and eigenvectors

 

이렇게 구한 eigenvalues, eigenvectors 그리고 matrix exponential까지..

이에 대한 물리적인 의미에 대해서 찾아보려고 합니다.

 

어떤 dynamic (LTI) system에 대해서 우리는 equation을 만들었고, equation의 solution도 구했습니다.

(normal-engineer.tistory.com/32 참고)

 

이 솔루션은 matrix exponential이 포함된 식이었습니다. 아래 그림을 보시면 DE에 대해 free response를 구한 답입니다.

 

이와 같이 dynamic system에 대해 $x_{1}(t)$와 $x_{2}(t)$에 대해 결과를 구할 수 있습니다.

 

그런데 이제 다른 방식으로 생각을 해보려고 합니다. initial condtion을 eigenvector의 선형조합으로 나타냈을 때, solution이 어떻게 수렴하는지에 대해서 보겠습니다.

 

즉 식에서 나온 대로, x(t)는 초기 조건에서 $e^{-2t}$의 속도로, $(1,-2)^{T}$ 방향으로 수렴해가는 성분과 $e^{-t}$의 속도로, $(1,-1)^{T}$ 방향으로 수렴해가는 성분이 존재합니다.

 

 

그림으로 나타내면 다음과 같습니다. 점차 0으로 수렴해가는 것은 맞지만, 어디에서 출발하느냐에 따라 $x_{1}(t), x_{2}(t)$의 경로가 달라집니다. 또한 $t_{1}$ 방향은 $e^{-2t}$라는 더 빠른 속도로 수렴하기 때문에 $t_{1}$방향으로 먼저 따라가다가 나중에 $t_{2}$방향으로 수렴하는 것을 알 수 있습니다.

 

이것을 확인해보기 위해 MATLAB 코드로 그래프를 그려보았습니다.

 

예를 들어서 (1,1)에서 시작하는 이 Dynamic system은,

$t_{1}$ 방향이 더 빨리 수렴하므로 그 방향으로 따라가다 나중에는 $t_{2}$방향으로도 수렴하면서 원점에 도달하게 됩니다. (MATLAB코드는 아래에)

 

 

---

지금까지는 모든 eigenvalues가 distinct하다고 가정하고 설명했습니다.

하지만 diagonal matrix 외에도 jordan form이거나 complex eigenvalue matrix인 경우도 있습니다.

 

1) Complex eigenvalue인 경우 

transformation matrix T를 어떻게 설정하느냐에 따라 두 가지 방법으로 나뉩니다. (결과는 동일)

 

첫번째 방법은, 기존에 했던 eigenvalue를 구해서 diagonal matrix를 구하는 방법과 동일하게 하는 것입니다.

 

 

두 번째 방법은, 기존에 알고 있었던 complex eigenvalue matrix를 활용하는 방법입니다.

 

 

 

2) Repeated eigenvalue인 경우

 

$det(\lambda I-A)=(\lambda-\lambda_{m})^{3}$이라고 가정하겠습니다. 즉, 3중근일 때

matrix exponential은 $nullity(\lambda_{m}I-A)$에 따라서 달라집니다. 

 

 

 

따라서 중근이 나오더라도 각 case에 대해서 구하면 됩니다.

 

---

MATLAB code

f = @(t,X) [X(2); -2*X(1)-3*X(2)]; 
x1 = linspace(-10,10,20); 
x2 = linspace(-10,10,20); 
[x, dx] = meshgrid(x1,x2); 

u = zeros(size(x)); 
v = zeros(size(x)); 

% we can use a single loop over each element to compute the derivatives at 
% each point (y1, y2) 
t=0; % we want the derivatives at each point at t=0, i.e. the starting time 
for i = 1:numel(x) 
    Xprime = f(t,[x(i); dx(i)]); 
    u(i) = Xprime(1); 
    v(i) = Xprime(2); 
end 

for i = 1:numel(x1) 
    a1(i) = -x1(i); 
    a2(i) = -(x1(i)-1)+1; 
    b1(i) = -2*x1(i); 
    b2(i) = -2*(x1(i)-1)+1; 
end 

quiver(x,dx,u,v,'r'); figure(gcf) 
xlabel('x1') 
ylabel('x2') 
grid; 
axis tight equal; 
hold on 
[ts,ys] = ode45(f,[0,50],[1;1]); 
plot(ys(:,1),ys(:,2)) 
plot(ys(1,1),ys(1,2),'bo') % starting point 
plot(ys(end,1),ys(end,2),'ks') % ending point 
plot(x1,a1); 
plot(x1,a2); 
plot(x1,b1); 
plot(x1,b2); 
axis([-10 10 -10 10]); 
hold off