행렬 응용 - Linear Difference Equation

▶차분방정식: 

차분방정식은, 시간이 지남에 따라 상태가 변하는 문제를 방정식으로 규정한 것이다. 

$ \begin{equation} u_{k+1} = Au_k \end{equation}$

위의 공식은 $u_k$였던 상태가 단위 시간 1이 지났을 경우 행렬 A를 곱하는 것과 같이 변한다는 것을 의미한다. 

위 공식을 $k=0$일때부터 순서대로 전개해보면, 아래와 같이 상태가 변함을 알 수 있다.

$$ \begin{align*} u_1 &= Au_0 \\ u_2 &= Au_1 = A(Au_0) = A^2 u_0 \\ u_3 &= Au_2 = A(A^2u_0) = A^3 u_0 \\ \vdots \\ u_k &= A^k u_0 \end{align*} $$

따라서, 시간 $k$에서의 상태는 행렬 $A$의 거듭제곱을 초기값에 곱한 결과와 동일하다. 

 

▶고유값과 고유벡터를 이용한 해 구하기: 

이런 선형차분방정식(linear difference equation)은 고유값과 고유벡터를 이용하여 아래와 같이 해를 구할 수 있다. 

 

• $x_{k+1} = a x_k + b x_{k-1} (k \geq 1) $라고 한다면, 이에 대응하는 또 다른 방정식을 고려하고, 이를 행렬로 표현한다.
• $$ 
  \begin{align*} x_{k+1} &= a x_k + b x_{k-1}, \\ 
  x_k &= 1 x_k + 0 x_{k-1}, k \geq 1 
  \\ \begin{bmatrix} x_{k+1} \\ x_k \end{bmatrix} &= \begin{bmatrix}  a & b \\ 1 & 0 \end{bmatrix} \begin{bmatrix}  x_{k} \\ x_{k-1} \end{bmatrix} \\
  \boldsymbol{x}_K &= M \boldsymbol{x}_{k-1}, \boldsymbol{x}_k = \begin{bmatrix} x_{k+1} \\ x_k \end{bmatrix}, M = \begin{bmatrix}  a & b \\ 1 & 0 \end{bmatrix}
  \end{align*} 
  $$
• 이때, 행렬차분방정식의 해는 $\boldsymbol{x}_k = \boldsymbol{M}^k \boldsymbol{x}_{0}, (k \geq 1)$로 주어진다. 
• 만약 행렬 $\boldsymbol{M}$이 대각화(diagonalize)가 가능하다면, $M \boldsymbol{v}_1 = \lambda_1 \boldsymbol{v}_1, \quad M \boldsymbol{v}_2 = \lambda_2 \boldsymbol{v}_2$ 를 만족하는 고유값 $(\lambda_1, \lambda_2)$와 고유벡터 $ (\boldsymbol{v}_1, \boldsymbol{v}_2) $를 찾을 수 있다.       
  $$ 
  P^{-1} M P = \begin{bmatrix} \lambda_1 & 0 \\ 0 & \lambda_2 \end{bmatrix}, \quad M^k = P \begin{bmatrix} \lambda_1^k & 0 \\ 0 & \lambda_2^k \end{bmatrix} P^{-1}, \quad P = (\boldsymbol{v}_1, \boldsymbol{v}_2)
  $$
• 행렬의 대각화란, 행렬을 대각요소만 가진 대각행렬로 만들어주는 것이며, 정방행렬의 대각화는 익히 잘 아는 고유값분해(eigen decomposition)으로 불리기도 한다. $n \times n$ 정방행렬 $M$이 $n$개의 선형적으로 독립인 고유벡터들을 가진다면 $P^{-1} M P = \Lambda = \begin{bmatrix} \lambda_1 & 0 \\ 0 & \lambda_2 \end{bmatrix}$가 성립한다. 
  따라서, 위에서 본 것과 같이 차분방정식의 해는 아래와 같이 성립한다.
  $$ \boldsymbol{x}_k = M^k x_0 = (P\Lambda P^{-1})(P\Lambda P^{-1})\cdots(P\Lambda P^{-1}) x_0 = P\Lambda^k P^{-1} x_0 $$
• ex) $x_{k+1} = 5x_k - 6 x_{k-1}, x_0 = x_1 = 1 $ : 
  $$
  \boldsymbol{x}_{k} = \begin{bmatrix} 5 & -6 \\ 1 & 0 \end{bmatrix}^k \boldsymbol{x}_0 \quad \textrm{where} \,\,\, \boldsymbol{x}_0 = \begin{bmatrix} 1 \\ 1 \end{bmatrix}.
  $$ 우선 대각화를 이용해서 풀려면 행렬 $M$이 대각화가 가능한지 확인해 봐야 한다. 따라서 먼저 고유값과 고유벡터들을 구해야 한다. 
  $$
  \begin{align*}
  |M - \lambda I| &= 0 \\
  \begin{vmatrix} 5-\lambda & -6 \\ 1 & 0-\lambda \end{vmatrix} &= (5-\lambda)(-\lambda) + 6 = (\lambda-3)(\lambda-2) \\ \therefore \lambda_1 &= 2, \lambda_2 = 3.
  \end{align*}
  $$ then, $$
  \begin{align*}
  \begin{bmatrix} 3 & -6 \\ 1 & -2 \end{bmatrix} \boldsymbol{v_1} &= \boldsymbol{0} \quad \therefore \boldsymbol{v}_1 = \begin{bmatrix} 2 \\ 1 \end{bmatrix} \\
  \begin{bmatrix} 2 & -6 \\ 1 & -3 \end{bmatrix} \boldsymbol{v_2} &= \boldsymbol{0} \quad \therefore \boldsymbol{v}_2 = \begin{bmatrix} 3 \\ 1 \end{bmatrix}
  \end{align*}
  $$ 고유벡터들이 서로 독립임을 확인할 수 있다. 따라서 행렬 $M$은 대각화 가능하다. 따라서, 다음과 같이 $\boldsymbol{x}_k$를 계산할 수 있다. $$ \begin{align*}
  \boldsymbol{x}_k = M^k x_0 &= P\Lambda^k P^{-1} x_0 \\ &= \begin{bmatrix} 2 & 3 \\ 1 & 1 \end{bmatrix} \begin{bmatrix} 2^k & 0 \\ 0 & 3^k \end{bmatrix} \begin{bmatrix} 2 & 3 \\ 1 & 1 \end{bmatrix}^{-1} \begin{bmatrix} 1 \\ 1\end{bmatrix} = \begin{bmatrix} 2 \times 2^{k+1} - 3^{k-1} \\ 2 \times 2^k - 3^k \end{bmatrix}
  \end{align*} $$

▶위 개념을 코드로 다시 짜보기 : 

import numpy as np
from numpy.linalg import eig
from numpy.linalg import matrix_power, multi_dot, inv

def solve_diff(A, k, x0, x1):
    w, v = eig(A)
    eigenval = w
    eigenvec = v
    diagmat = np.diag(w)
    diagmat = matrix_power(diagmat, k)
    x_vec = np.array([[x0],[x1]])
    x_k = multi_dot([eigenvec, diagmat, inv(eigenvec), x_vec])
    return(x_k[1])
    
a = np.array([[5, -6], [1, 0]])
solve_diff(a, 3, 1, 1)