이번 포스팅은 지난 강의 Lecture 23-(1)에 이은 미분방정식과 선형대수의 두 번째 이야기이다. 앞선 강의와 이어지는 내용이니만큼 앞 강의 내용을 잘 이해하고 보는 것을 추천한다. 

 

 

4. 행렬 지수함수(Matrix exponential)

 

- Diagonalization for differential equation

 

지난 강의에서 계속 다루어왔던 아래의 식을 다시 생각해보자. 

 

 

식 (18)의 행렬 A는 벡터 u와 그의 미분값 u'을 연결시키는(couple)역할을 한다. 그런데 우리는 앞선 강의에서 A의 고유값과 고유벡터를 구하여 지수함수(exponential function)로 이루어진 일반해(general solution)를 도출한바 있다. 일반해는 각각의 상수와 고유값, 그리고 고유벡터들로 이루어진 식이며, 선형조합(Linear combination)의 형태로 분리가 되어 있다. 이는 기존의 행렬 A가 고유값/고유벡터에 의해서 분리(uncouple)가 되는 것이다. 이것이 무슨 말인지 이해하기 위해 식 (18)을 고유벡터의 관점에서 다시 정리해보자. 

 

 

식 (19)는 (18)의 u를 고유벡터와 가중치(weight)곱의 선형조합에 대한 식으로 정리한 것이다. 이를 행렬에 대한 식으로 정의한 것이 (19.1)이다. 여기서 v는 고유값에 대한 지수함수를, S는 고유벡터들로 이루어진 고유벡터행렬을 의미한다. u=Sv를 풀어서 써보면 식 (19.2)와 같이 일반해로 정리되기 때문에 식에 무리는 없다. 이제 식 (19)의 내용을 (18)에 대입하여 A를 분해하여 정리해보자. 

 

 

(20.1)은 (18)에서 u대신 Sv를 대입하여 정리한 것이다. 이때 S는 상수이기 때문에 좌변의 미분식에서 S는 앞으로 뺄 수 있다. 양변의 좌측에 S의 역함수를 곱해주면 (20.2)와 같이 정리할 수 있는데, 우변의 S^-1AS가 어디서 많이 본 형태이다. 바로 Lecture 22의 행렬의 대각화(diagonalization)에서 배웠던 형태인데, S^-1AS는 Λ(capital lambda)와 같음을 배웠다. 이를 통해 (20.2)를 다시 (20.3)과 같이 정리할 수 있다. 

 

식 (20.3)을 보면 dv/dt=Λv가 되는데, 식 자체의 의미를 놓고 보면 벡터 변수 v를 시간으로 미분한 것은(시간의 변화에 따른 변수 변화량의 비율) 변수에 어떤 상수값을 곱한 것과 같다이다. 이 의미가 성립하려면 v는 식 (19.1)과 같이 지수함수로 정의되어야 한다. 지수함수는 미분한 값이 바로 자기 자신이 되기 때문이다(지수함수의 지수부의 상수가 미분하면 지수함수 앞에 곱해지고 이것이 람다가 되는 꼴). v는 벡터이기 때문에 여러 단일 변수를 동시에 정의할 수 있으며 아래 식은 이에 대한 내용을 정리한 것이다. 

 

 

식 (21)은 dv/dt=Λv를 2x2크기로 가정하고 식을 풀어 정리한 것이다. 실제로 위의 식과 같이 v에 지수함수를 대입하여 t로 미분해보면 식이 성립함을 알 수 있다. 여기서 중요한 것은 위의 식을 활용하여 미분방정식을 푸는것이다. 우리가 결국 구하고자하는 것은 u이고, v는 식 (19)와 같이 u=Sv의 꼴로 나타냈을 때 S의 각 고유벡터들에 곱해지는 계수(coefficient)에 해당한다. 그러므로 먼저 v의 시간 t에 대한 함수를 정의해보자. 지난 강의 Lecture 23-(1) 에서 배웠듯이 미분방정식의 해(solution)는 지수함수(exponential function)를 도입하여 정의할 수 있다. 단일 변수 대신 다변수에 대해 정의한 행렬에 지수함수를 도입한다는 것이 다를 뿐 기본적인 개념은 같다. 

 

 

(22.1)은 식 (21)에서 도출한 식이고, 여기서 양변에 dt를 곱해주고 1/v를 곱해준 뒤 적분(integral)을 해준 것이 (22.2)이다. 적분을 하면 (22.3)이 되고, 좌변의 C1을 우변으로 이항하고 자연로그를 없애주기 위해 양변에 지수함수(exponential function)를 취하면 (22.4)가 된다. 절대값은 +, 혹은 -값이 되지만 임의의 값이 될 수 있으므로 +값을 선택하여 e^(c2-c1)자체를 상수 C로 치환하여 (22.5)를 정의할 수 있다. 

여기서 우리는 식을 미분방정식의 해를 구하기 위한 함수의 형태로 정리해야 한다. 즉 특정 시간 t에서 해당 미분방정식의의 값이 무엇인지 구할 수 있도록 v의 t에 대한 함수 v(t)의 꼴로 정의해야 한다. v(t)의 꼴로 정의하는 것은 반드시 초기값이 있어야 한다. v에서 t=0인 경우는 지수함수부분이 1이 되므로 상수 c1, c2만 남게 된다. 이 상수 벡터가 v의 초기값이 되고 (22.5)의 C를 대체하여 생각할 수 있다. C를 v의 초기값으로 대체하여 (22.6)과 같이 v(t)를 정의하였다. 

 

v(t)를 구했으니 이제 u(t)를 구해보자. 어떻게 구할 수 있을까? 바로 식 (19)의 u=Sv로부터 구할 수 있다. 아래의 식을 보자. 

 

 

u=Sv를 u(t)의 꼴로 정리하면 시간의 변수 t를 포함하고 있는 v도 역시 v(t)의 형태로 변해서 (23.2)와 같이 정리할 수 있다. 그 다음 (23.2)의 v(t)를 (23.1)의 v(t)의 우변으로 식을 치환하면 (23.3)과 같이 정리가 된다. u(t)도 마찬가지로 초기값 u(0)가 포함되어야 하는데, u(0)=Sv(0)이고 이로부터 (23.4)와 같이 v(0)를 정의할 수 있다. 여기서 정리한 v(0)를 (23.3)에 대입하여 정리하면 최종적으로 (23.5)와 같이 정리할 수 있다. 

 

정리하자면 우리가 최초에 풀고자 했던 문제는 식 (18)에서 정의한 것과 같이 다변수(multi-variable) 함수에 대한 미분방정식을 푸는 것이었다. 여기서 A는 u와 u의 미분값 사이를 정의하며, 이들을 연결해주는(couple) 역할을 한다. 특정 시간 t에서의 미분방정식의 해(solution)를 구하기 위해선 u를 시간 t에 대한 함수 u(t)의 꼴로 정의를 해줘야 하는데, 이때 사용한 방법이 이 행렬 A의 고유값과 고유벡터를 이용하여 정리하는 것이다. A의 고유값과 고유벡터를 이용하여 대각화를 하고, 여기에 지수함수를 적용하여 식 (23.5)와 같이 미분방정식의 해를 구할 수 있는 식을 만든 것이다. 결국 A의 고유값과 고유벡터를 이용하여 최초의 식 A를 분해하고 식 (19)와 같이 지수함수를 도입하여 다변수 미분방정식의 해(solution)를 구하기 위한 식을 (23.5)와 같이 도출하였다

 

그런데 식 (23.5)의 가운데 부분은 다시 아래와같이 정리할 수 있다. 

 

 

u(t)의 가운데 대각화된 식을 e^At라고 할 수 있는데, 이와 같이 지수함수를 밑으로 하는 형태의 행렬을 행렬지수함수(matrix exponential)라 한다. 뒤이어 알아볼 내용은 식 (24.1)의 행렬지수함수이다

 

행렬지수함수를 공부하기 전에 먼저 알아야 할 내용이 있다. 바로 테일러 급수(Taylor series)이다. 인터넷에 테일러 급수에 대한 수많은 자료가 있지만, 매우 중요한 개념이니만큼 간략히 정리하고 넘어가도록 하자. 

 

 

- Taylor series

 

테일러 급수(Taylor series)란 어떤 함수를 특정 위치 x=a에서 근사(approximation)하는 방법이다. 테일러 급수의 정의는 아래의 식과 같다. 

 

 

식 (25)를 풀어서 설명해보자면 n번 미분이 가능한 어떤 함수 f(x)가 있을 때, 이 함수의 특정 위치에서의 값, 즉 x=a에서의 값은 f(x)와 그의 미분값들로 계산한 항들의 무한대의 합으로 표현할 수 있다는 것이다. 이때 각 항은 미분차수값의 factorial로 나누어주어야 한다. 즉 1차 미분이면 1 factorial(1!), 2차 미분이면 2 factorial(2!)로 나누어줘야 한다. 테일러 급수의 식이 왜 이와 같이 표현되었는지를 이해하려면 멱급수(power series)의 전개부터 이해를 해야하는데, 이 부분은 본 포스팅의 범위를 넘어가기 때문에 다음의 동영상을 보는 것을 추천한다(Power series/Euler's great formula). 

 

위에서 설명한 테일러 급수에서 한 가지 중요한 것은  테일러 급수가 근사하는 것은 함수 전체가 아니라 함수의 특정 위치(x=a)에서의 근사만을 의미한다. 물론 어떤 함수의 경우엔 근사하는 위치 주변에서의 값도 잘 맞아떨어질 순 있겠지만, 기본적으로는 x=a에서의 근사를 하기 위함이다. 이와 같은 설명만으론 이해하기 어려울 수 있으므로 예를 한 번 들어보자. 

 

 

식 (26)의 함수 f(x)는 3차 함수이며 3번 미분이 가능하기 때문에 3차 미분까지만 고려하면 된다. 그 이상은 0이 되기 때문에 의미가 없다. 위의 함수를 x=0일때와 x=-2일때 각각 근사를 해보자. 참고로 x=0일때의 테일러 급수를 매클로린 급수(Maclaurin series)라 한다. 

 

 

식 (27)은 (26)의 함수 f(x)를 x=0에서 근사한 것이다. (27.1)은 f(x=0)일때의 값과 3차 미분까지의 값들을 나열하였다. (27.2)는 원래의 함수 f(x)를 테일러 급수로 근사한 결과이다. 각 미분항에 곱해지는 (x-a)^n에서 x=0에서의 근사이기 때문에 a에는 0이 대입된다. 정리한 결과 원래의 함수 f(x)가 그대로 복원되고 0을 대입하면 마찬가지로 2가 나오는 것을 볼 수 있다. 

 

이제 x=-2일 때의 근사(approximation)를 수행해보자. 

 

 

식 (28)은 x=-2에서 테일러 급수를 통해 f(x)를 근사한 모습이다. 마찬가지로 원래의 함수를 잘 복원하는 것을 볼 수 있다. 

 

그런데 여기서 어떤 사람들은 이렇게 생각할 수 있다. "이렇게 그대로 다시 복원할거면 도대체 테일러 급수를 왜 사용하는거지? 그냥 그대로 사용하면 될텐데..". 사실 위에서 예로든 함수 f(x)는 다항식의 유한번의 사칙연산으로 구할 수 있는 대수함수(algebraic function)이고, 따라서 테일러 급수를 통해 원래의 함수를 그대로 복원할 수 있다. 테일러 급수가 진정 빛을 발하는 경우는 대수적 연산만으로는 정확한 값을 구할 수 없고 극한의 과정을 거쳐야만 정확한 값을 구할 수 있는 함수, 이를 테면 지수 함수(exponential function), sin, cos, log등 초월 함수(transcendental function)를 구할 때이다. 초월 함수의 대표적인 테일러 급수의 예로 지수함수(exponential series)와 기하급수(geometric series)에 대한 식을 살펴보자. 

 

 

식 (29.1)과 (29.2)는 각각 지수함수와 기하 급수를 나타내며, (29.2)의 기하급수는 x의 크기가 1보다 작을 때에만 좌변과 우변의 식이 성립(수렴)한다. 이 외에도 sin, cos, log등의 다양한 초월 함수들이 존재하며 자세한 사항은 테일러 급수 위키(Taylor series)를 참고하기 바란다. 이번 강의를 위해선 일단 위의 두 함수에 대한 테일러 전개(Taylor expansion)를 알면 된다. 참고로 위의 지수함수와 기하급수 두 개의 테일러 급수는 각각 x=0일 때의 매클로린 급수이다. 

 

그렇다면 테일러 급수가 이러한 초월함수에서 힘을 발휘하는 이유는 무엇일까? 그것은 바로 근사(approximation)하는 능력에 있다. 제아무리 성능 좋은 수퍼컴퓨터라해도 무한번의 연산을 수행할 수는 없을 것이다. 따라서 우리는 해당 함수의 값을 원하는 정도의 정밀도가 나올 때까지만 근사하여 구하면 된다. 이때 이러한 근사의 과정을 테일러 급수를 통해서 수행하게 된다. 실제로 공학용 계산기에서 sin함수의 값은 테일러 급수를 통해 어느 정도로 근사된 값이 계산되어 나오는 것이다. 테일러 급수라는 개념이 없었다면 이러한 초월함수들의 계산 과정에 큰 어려움이 있었을 것이다. 근사(approximation)에 대한 이해를 돕기 위해 (29.2)의 기하 급수의 수렴에 대한 예를 살펴보자. 

 

 

식 (30.1)은 x=1/2일 때의 기하급수의 계산 결과를 나타낸다. 기하급수의 특성상 |x|<1일 때 더해가는 항, 즉 n을 늘려나갈수록 특정 값에 수렴하게 되는데, (30.2)는 두 번째 항까지, (30.3), (30.4)는 각각 세 번째, 네 번째 항까지 더한 결과를 나타낸다. 항이 늘어날 수록 수렴값 2에 가깝게 되는 것을 볼 수 있다. 사실 무한대로 더해야지만 2에 완전히 수렴하겠지만, 무한대로 더하는 것은 불가능하기 때문에 만족할만한 수준의 정밀도를 설정해놓고 그에 상응하는 항까지 더하여 근사를 하면 된다. 

 

이 외에도 테일러 급수는 어떤 복잡한 모델을 단순화 시켜서 풀거나, 미분방정식의 해를 구할 때, 정적분을 구할 때 등 과학과 공학에 있어서 다양하게 활용되는 굉장히 중요한 개념이므로 잘 이해하도록 하자. 

 

 

- Matrix exponential

 

앞서 테일러 급수를 간략히 공부했으니 이제 행렬지수함수(Matrix exponential)를 공부해보도록 하자. 미리 말하지만 행렬지수함수를 공부하는 이유는 식 (24)에서와 같이 미분방정식의 해(solution)에 대한 식을 도출하는 과정에서 행렬지수함수가 나오기 때문에 이를 정의하는 법을 알기 위함이다. 우리가 알아볼 행렬지수함수의 형태는 아래의 식과 같다. 

 

 

우리가 학창시절에 공부하면서 봤던 지수함수의 형태는 대부분 지수함수 e의 지수부에 상수, 혹은 변수가 있는 형태였다. 그런데 식 (31)에 정의된 행렬지수함수는 지수함수의 지수부에 행렬(matrix)이 있는 경우이다. 여기서 행렬 A는 고정된 값이고 t가 기존의 x와 같은 변수이다. 이러한 행렬지수함수는 어떻게 정의할 수 있을까? 그 힌트는 바로 앞서 공부했던 테일러 급수에 있다. 식 (29.1)의 지수함수의 테일러 전개가 (31)의 행렬지수함수에도 그대로 적용되는 것이다. 다만 단일 변수 대신 행렬의 형태로 전개가 되는 것이다. 아래의 식을 살펴보자.  

 

 

비교해보기 쉽도록 원래의 식 (29.1)과 함께 정리하였다. (29.1)의 1이 (32)에선 단위 행렬 I가 되고, x는 At가 되며 나머지 항들도 같은 패턴을 보이는 것을 알 수 있다. 이때 (29.1)의 1+x 부분이나 (32)의 I + At등의 부분은 사실은 n^0 !, n^1 !로 각각 나누어 준 것이다. 지수함수와 같은 패턴으로 행렬지수함수(Matrix exponential)이 계산된다는건 알겠는데, 정확히 어떤 형태로 나온다는건지 감이 오지 않을 수 있다. 아래의 행렬을 이용하여 행렬지수함수를 계산해 보도록 하자. 

 

 

식 (33.1)의 행렬 A와 t=1임을 가정하여 행렬지수함수를 계산하였다. 계산 항은 총 10차 까지 계산하였다. 계산 결과값은 (33.3)과 같은 결과를 보인다. 이 값이 과연 맞는 값인지 MATLAB을 통해 확인해보자. MATLAB에서 행렬지수함수를 계산하기 위해선 expm()함수를 이용하면 된다. 아래 그림은 계산 결과와 소스코드이다. 

 

  

 

Fig. 1 Matrix Exponential계산 결과와 소스코드. MATLAB 내장함수결과(ans =), 테일러급수 n=10차 계산결과(M =)

 

Fig. 1의 좌측은 MATLAB을 이용한 계산결과를, 우측은 소스코드를 나타낸다. 좌측 결과에서 ans = 는 MATLAB의 내장함수인 expm()함수를 이용한 계산결과이고 M = 은 테일러급수를 이용하여 10차까지 근사한 결과이다. 근사값이 내장함수값과 거의 일치함을 보이며, 차수를 더해갈 수록 오차는 줄어들게 되고 보다 정확한 근사를 할 수 있다. 이와 같이 행렬지수함수가 테일러 급수로 표현될 수 있음을 증명하였다. 

 

 

위에서 행렬지수함수를 테일러 급수를 통해 정의하는 법을 배웠는데, 우리는 앞서 지수함수와 더불어 식 (29.2)의 기하급수(geometric series)에 대해서도 배웠다. 기하급수 역시 식 (29.2)의 기존의 형태에 행렬을 대입하여 정리할 수 있다. 아래 식을 보자. 

 

 

비교를 위해 역시 기존의 식 (29.2)와 함께 정리했다. (29.2)의 1/1-x를 다시 쓰면 (1-x)^-1와 같이 쓸 수 있는데, 식 (34)와 같이 행렬에 대해 정리했을 때 이는 역함수로 해석할 수 있다. 즉 (34)의 (I-At)의 역함수로 생각할 수 있다. 1/(I-At)가 아니다. (29.2)에서 x의 크기가 1보다 작은 경우에만 식이 성립한다고 했다. 마찬가지로 행렬의 형태로 정리한 (34)에서도 이와 같은 조건이 있는데 바로 At의 모든 고유값이 1보다 작은 경우에만 식이 성립한다는 것이다. 이때 중요한것은 행렬 A에 대한 고유값이 아니라 상수 t가 곱해진 At의 모든 고유값이 1보다 작아야 한다는 것이다. 실제 값을 대입하여 계산 결과를 살펴보자. 

 

 

 

위의 계산 결과가 실제로 (I-At)의 역함수를 계산했을 때와 같은지 MATLAB을 통해 확인해보자. 

 

 

Fig. 2  기하급수(Geometric series)의 Matrix Exponential계산 결과와 소스코드. MATLAB 계산결과(ans =), 기하급수 n=100차 계산결과(M =)

 

Fig. 2의 왼쪽은 MATLAB 계산결과이고 오른쪽은 소스코드이다. V와 D는 각각 At의 고유벡터와 고유값을, mag1과 mag2는 At의 λ1과 λ2의 크기를 각각 의미한다. 

 

다음으로 ans = 는 (I-At)의 역함수를 계산한 결과이며, M은 t=2.2일 때 기하급수를 테일러 급수를 이용하여 계산한 결과이다. 이때 차수는 n=100차 까지 계산하였다. 결과적으로 역함수를 계산한 결과와 급수를 통해 근사한 결과인 M을 비교해보면 거의 같음을 알 수 있다. 물론 차수를 높일 수록 값은 더욱 같아진다. 이로써 기하급수를 행렬지수함수의 형태로 정리한 식 (34)가 조건이 맞을 경우 식이 성립함을 증명하였다. 위의 MATLAB코드에서 t와 A의 값을 바꾸어가며 테스트 해보면 mag1과 mag2가 1이 넘는 경우엔 식이 성립하지 않음을 볼 수 있으니 테스트해보기 바란다. 

 

급수(series)의 관점으로 봤을 땐 식 (32)가 (34)보다는 훨씬 뛰어난 특성을 보인다. (32)의 경우엔 A나 t가 얼마나 크던 간에 항을 더해갈수록 큰 수로 나누어주기 때문에 나중으로 갈 수록 0에 가까워지고 결국 수렴하게 된다. 반면에 (34)의 경우엔 At의 고유값이 1보다 큰 경우 발산할 수 있다. 기하급수의 행렬지수함수를 알아봤지만, 이번 강의에서 우리가 눈여겨봐야할 급수는 바로 식 (32)의 행렬지수함수이다. 이 행렬지수함수를 이용하여 앞에서 다루었던 식 (24)의 행렬지수함수를 증명해보도록 하자. 그 사이에 많은 내용을 다루었기 때문에 편의를 위해서 식 (24)를 다시 써보도록 하자. 

 

 

우리가 알고자 하는 것은 식 (24.2)가 어떻게 성립하는지, 즉 $Se^{\Lambda t}S^{-1}$를 어떻게 $e^{At}$로 정의할 수 있는지이다. 식 (32)의 행렬지수함수를 이용하여 이를 정의할 수 있다. 아래의 식을 보자. 

 

 

식 (36.2)는 (36.1)의 행렬지수함수에서 A를 $S \Lambda S^{-1}$로 치환한 것이다. $A=S \Lambda S^{-1}$는 Lecture 22의 행렬의 대각화 부분에서 배웠던 내용이다. (36.2)의 첫 번째 항은 t의 0제곱과 zero factorial이기 때문에 지수함수 e가 사라져서 고유벡터행렬 S만 남은 것이다. (36.2)에서 고유벡터행렬 $S$와 $S^{-1}$는 상수이기 때문에 양쪽으로 뺄 수 있다. 이렇게 S들을 빼서 정리한것이 (36.3)이고, 괄호안에는 고유값 행렬(Λ)에 대한 식이 남는데, 이는 $e^{\Lambda t}$로 볼 수 있다. 따라서 (36.4)와 같이 정의할 수 있고 식 (24.1)을 증명한 것이다. 

 

한 가지 기억해야 할 것은 식 (36.1)은 언제나 성립하는 반면 식 (36.4)는 항상 성립하는 것은 아니라는 것이다. (36.4)가 성립하기위해선 A가 대각화가 가능해야 한다는 것이다. 다시말하면 행렬 A는 n개의 독립(independent)인 고유벡터(eigenvector)를 가져야 한다. 이는 행렬의 대각화를 배울 때 강조했던 내용이며 (36.4)가 성립하기 위해선 반드시 충족시켜야 할 조건이다. 중요한 내용이므로 다시 한 번 강조하니 꼭 기억하도록 하자. 

 

 

이즈음에서 다시 한 번 생각해보자. 우리가 이렇게 행렬지수함수를 계산하는 법을 공부하는 이유는 무엇인가? 바로 선형연립미분방정식(simultaneous linear differential equation)의 해(solution)를 구하기 위함이며, 그 해는 식 (24)와 같이 u(t)=의 형태이다. u(t)의 식에서 결국 구해야 할 것은 주어진 행렬 A를 기반으로 식 (24.1)의 행렬지수함수를 구하는 것이다. 사실 (36.1)과 같이 테일러급수의 형태로 근사하여 해를 구할 수도 있다. 그러나 계산과정에 있어 비용이 많이 드는 것이 사실이다. 따라서 고유값과 고유벡터를 기반으로 한 (36.4)와 같은 형태로 해를 구하는 것이 훨씬 효율적이다(n개의 독립인 고유벡터가 존재한다는 가정하에서). 

 

행렬 A의 고유값과 고유벡터를 구하는 것은 어렵지 않지만, 문제는고유값 행렬(Λ)에 대한 행렬지수함수이다. 결국 람다에 대한 행렬지수함수를 얻기 위해 (36.1)의 과정을 거쳐야 한다면 (36.4)의 형태가 큰 장점은 없어보인다. 하지만 고유값 행렬(Λ)은 대각행렬이고, 이에 대한 행렬지수함수는 그 계산방법이 훨씬 간단하다는 이점이 있다. 아래의 식을 보자. 

 

 

식 (37)은 고유값 행렬(Λ)의 행렬지수함수를 나타낸 것이다. 보다시피 고유값만 알면 쉽게 계산할 수 있는 형태이다. 과연 위의 형태대로 계산하는 것이 맞는지 MATLAB을 통해 확인해보자. 

 

 

Fig. 3 고유값 행렬(Λ)의 행렬지수함수 계산결과(Left)와 MATLAB소스코드(Right)

 

 

Fig. 3은 고유값 행렬의 행렬지수함수의 계산 검증 결과를 나타낸다. 왼쪽의 계산결과에서 V와 D는 각각 행렬 A의 고유벡터행렬과 고유값 행렬이고, ans는 MATLAB 내장함수인 expm()으로 고유값 행렬의 행렬지수함수를 계산한결과를, 그리고 expD는 식 (37)의 형태로 직접 계산한 결과이다. 비교해보면 정확히 일치하는 것을 볼 수 있다. 이로써 계산이 훨씬 간편해졌다. 물론 코드가 더 추가되긴 했지만 행렬이 커짐에 따라 이 계산방식이 빛을 발할 것이다. (내장함수인 expm()의 계산 알고리즘이 이미 이와 같은 방식을 사용한다거나, 혹은 더 좋은 방법을 사용한다면 이야기는 달리진다)

 

 

정리하자면 최초에 식 (18)과 같이 행렬의 형태로 선형연립미분방정식이 있고 초기값(initial condition)이 주어졌을 때 우리는 식 (24)와 같은 형태로 해를 구할 수 있다. 해를 구하기 위해선 행렬 A의 고유값과 고유벡터를 기반으로 식 (24.2)와 같이 해를 정의할 수 있으며, 이때 고유값 행렬의 행렬지수함수는 (37)의 형태로 계산할 수 있다. 

 

마지막으로 중요한 포인트 한 가지를 다시 짚어보자. 바로 지난 강좌에서 일부 다루었던 내용이지만 중요한 부분이므로 다시 언급하도록 하겠다.

일반적으로 어떤 시스템의 해를 구한다는 것은 시간(t)이 지날수록 그 값이 수렴하는 형태일 것이다. 그렇다면 언제 수렴하는 것일까? 그 해답은 식 (37)의 고유값 행렬에 있다. 식 (24.2)에서 u(t)는 고유벡터행렬 S와 고유값 행렬로 정의되어 있는데, 여기서 고유벡터행렬은 상수이기 때문에 고정되어 있다. 움직이는 것은 t가 곱해져서 영향을 받는 고유값 행렬의 행렬지수함수인데, 식 (37)을 보면 대각원소가 각 고유값에 대한 지수함수로 정의되어 있다. 이것이 의미하는 것은 t가 양의 방향으로 나아간다고 가정했을 경우 해가 수렴하기 위해선 행렬 A의 모든 고유값들이 0보다 작아야 한다는 것이다. 보다 더 정확히 얘기하자면 행렬 A의 모든 실수부(real part)의 고유값이 0보다 작아야지만 시간 t가 지날 수록 해가 안정상태(stable)가 된다. 이는 어떤 시스템 행렬 A를 설계하는데에 있어서 중요한 포인트이므로 잘 기억하도록 하자. 

 

 

- Test

 

지금까지 배운 행렬지수함수(Matrix exponential)를 이용한 미분방정식의 해를 구하는 방법이 과연 맞는지 확인해보도록 하자. 지난 강의 Lecture 23-(1)  에서 직접 구했던 해를 이번 강좌에서 배운 행렬지수함수를 이용하여 구해서 결과가 같은지 확인해보자. 아래는 MATLAB을 이용해 두 방식을 테스트 한 결과이다. 

 

 

Fig. 4 직접 구한 해와 행렬지수함수를 이용한 해의 비교

 

 

Fig. 4의 왼쪽이 계산 결과, 오른쪽이 소스코드이다. 왼쪽의 결과에서 A는 행렬, t는 임의의 시간 4.75를 가정하였고 초기값인 v0는 동일하게 [1 0]T으로 설정하였다. Sol=로 표시된 값이 행렬지수함수를 이용하여 계산한 값이고 u=로 표시된 값이 기존에 직접 계산했던 값의 결과이다. 계산 결과가 동일한 것을 볼 수 있다. 

 

 

 

5. 마치며(continue)

 

이번 강좌에선 지난 시간에 이어 미분방정식과 선형대수에 대한 내용을 다루었다. 어떤 시스템의 미분방정식을 세우고 그 해를 u(t)의 꼴로 구하는 과정을 배웠다. 핵심은 결국 시스템 행렬 A의 고유값과 고유벡터에 있으며, 이들을 기반으로한 행렬지수함수를 통해 그 해를 정의할 수 있다. 또한 행렬지수함수를 테일러 급수를 이용하여 정의하였으며, 행렬의 대각화(Diagonalization)를 통해 행렬지수함수를 정의하는 과정을 배웠다. 다음 포스팅에서 미분방정식과 선형대수의 마지막 강좌를 하도록 하겠습니다. 

 

이번에 포스팅할 내용은 미분방정식(Differential equation)을 선형대수를 이용하여 푸는 방법이다. 이번 포스팅의 내용을 이해하기 위해선 고유값/고유벡터 강의와 대각화(Diagonalization)내용을 먼저 학습하고 오는 것을 추천한다. 특히 지난 강의 Lecture 22에서 다루었던 대각화와 행렬의 거듭제곱, 그리고 차분방정식(Difference equation)과 관련이 깊기 때문에 이를 확실하게 이해한 뒤 이번 강좌를 보는 것이 좋다. 본격적인 내용을 학습하기전에 미분방정식이란 무엇이며, 이번 강의에 필요한 관련 내용들에 대해서 간략히 살펴보고 가도록 하자. 

 

 

1. 미분방정식의 개요(Introduction to Differential equation)

 

- What is a differential equation?

 

미분방정식(Differential equation)이란 어떤 방정식(equation)에 도함수, 즉 미분(derivative)이 포함된 것이다. 이는 물리학에서의 운동방정식, 가령 스프링에 매달린 추의 운동 분석에 대한 방정식 등을 세우는데에 사용되기도 하고 그밖에 경제학, 생물학 등 다양한 분야에 걸쳐 사용될 수 있는 수학적 방법이다. 미분방정식은 일반적인 방정식보다 식을 세우고 이해하거나 해(solution)를 구하기가 굉장히 까다롭다. 미분방정식에 대해서 본격적으로 알아보기 전에 다른 형태의 방정식들은 어떤 것들이 있는지 살펴보자. 아래는 미분방정식이 아닌 일반적인 형태의 방정식들을 나열해논 것이다. 

 

 

식 (1.1)부터 차례로 선형(일차) 방정식(Linear equation), 이차 방정식(Quadratic equation), 삼차 방정식(Cubic equation)을 각각 나타내며, 최고차항에 따라 구분되는 형태이다. 이들 방정식은 미지수(unknown)와 그들의 앞에 곱해져있는 상수인 계수(coefficient)로 구성되어 있으며, 덧셈, 뺄셈 연산을 통한 조합으로 구성되어 있다. 이와 같이 어떤 미지수의 거듭제곱과 계수들의 조합으로 구성되어있는 식을 대수방정식(algebraic equation)이라고 한다. 미지수의 최고차항에 따라 식을 구분할 때 차수가 일차인 식을 특히 선형방정식이라고 하는데, 이는 식 자체가 선형성(Linearity)을 가지기 때문이다. 2차식부터는 비선형성(non-linearity)을 가진다. 

 

사실 엄밀히 따지면 식 (1.1)은 대수방정식은 아닌데, 미지수가 x와 y 두 개로 구성되어있기 때문이다. (1.1)처럼 미지수가 한 개 이상으로 구성되어 있는 경우 이들을 다항방정식(polynomial equation)이라고 하며 종종 다항식(polynomial)이라고 쓰인다. 다항방정식은 대수방정식을 포괄하기 때문에 보통 대수방정식이라는 말 보다는 다항식이라는 용어를 사용한다. 

 

 

이제 미분방정식이 어떤 형태인지 살펴보자. 미분방정식의 형태는 식 (1)의 다항식(polynomial)형태의 식에 미분(derivative)의 개념이 더해진 식이기 때문에 훨씬 복잡해진다. 아래의 식을 보자. 

 

 

 

식 (2)는 하나의 미분방정식을 여러 가지 형태로 표현한 것이다. 식 (2.1)은 derivative form이고, y의 미분식을 prime을 이용하여 표현하였다. prime하나당 한 번의 미분을 나타내며, 이차미분은 y의 double prime으로 표현한다. 

 

식 (2.2)는 differential form이고 미분을 dy/dx(delta x, delta y로 읽음)로 표현한다. 여기서 분모에 존재하는 미지수 y는 종속 변수(dependent variable)이고, 분자인 x는 독립 변수(independent variable)라고 한다. 이는 (2.2)에서 y는 x에 대한 함수이기 때문이다. 즉 y는 x가 변함에 따라 그 변화량이 얼마나 되는지를 봐야 하는 x에 종속되는 변수이고, x는 어디에도 종속되지 않는 독립적인 변수이다. 예를 들어 어떤 공장의 온도에 따른 가스의 압력에 대한 관계를 미분방정식으로 만들었다고 생각해보자. 여기서 온도를 x, 가스의 압력을 y라고 했을 때 온도의 변화에 따른 가스의 압력을 dy/dx와 같은 미분식을 넣어서 만들 수 있을 것이다. x로 표기된 온도는 독립적으로 우리가 조정하는 것이고, y로 표시된 가스의 압력은 x를 임의로 변화시킴에 따라 변화된 결과를 관찰하는 종속변수라고 이해하면 된다. 

 

마지막 (2.3)은 differential operator이고 미분을 D문자로 표기한다. 많이 사용되지는 않는 표현법이지만 이러한 표현법이 있다는 것 정도만 알아두자. 

 

이러한 미분방정식은 크게 두 가지 형태로 분류할 수 있다. 첫 번째는 ODE(Ordinary differential equation), 두 번째는 PDE(Partial differential equation)이다. 각각 한글로 상미분방정식, 편미분방정식으로 부른다. ODE는 1개 혹은 여러 개의 종속 변수(dependent variable)를 가질 수 있으나, 오직 단 한개의 독립 변수(independent variable)만을 가진다. 반면에 PDE는 ODE와 달리 여러 개의 독립 변수(independent variable)를 가질 수 있다. 단 한 개의 종속 변수를 가지느냐, 아니면 여러 개를 가질 수 있냐, 이것이 가장 큰 차이점이다. 아래 식을 보자. 

 

 

식 (3.1)은 (2.2)와 같은 식이고, 여러 개의 dependent variable을 가질 수 있으나, 오직 하나의 독립 변수(independent variable)인 dx만을 가지는 형태의 식을 ODE라고 한다. 일반적으로 미분방정식이라고 하면 상미분방정식(ODE)을 의미한다. 반면 식 (3.2)와 같이 여러 개의 독립변수를 가지는 형태를 편미분방정식이라 하며, 종속 변수를 서로 다른 독립변수로 각각 편미분하기 때문에 delta인 d대신에 round인 를 이용하여 표기한다. 식 (3.2)에선 하나의 종속 변수 u를 각각 x와 y로 2차 편미분 하였다. 이번 강의에서는 편미분방정식은 따로 다루지 않고 상미분방정식만을 다룬다. 

 

 

- Relationship between Differential equation and Exponential function

 

미분방정식(Differential equation)은 지수함수(Exponential function)와 밀접한 관련이 있다. 이번 강좌를 학습하기 위해선 이 둘의 관계를 아는 것이 좋다. 아래의 식을 보자. 

 

 

식 (4.1)은 t의 변화에 대한 y의 변화율, 즉 y를 t로 미분한것을 나타낸다. 이 변화율은 y의 값에 어떤 상수 k를 곱한 것으로 나타낼 수 있다. 이 미분방정식에서 y의 값을 구하기위해 우리는 지수함수(exponential function)를 이용할 수 있다. 식 (4.2)는 지수함수를 이용한 y의 일반해(general solution)이고, (4.2)에서 C는 y의 초기값(initial value)이며 k는 비례 상수(proportionality constant)이다. 

 

그렇다면 식 (4.2)는 어떻게 (4.1)의 일반해가 될 수 있는 것일까? (4.1)의 식을 약간만 변형시킨다면 (4.2)와 같이 만들 수 있다. 아래 식을 보도록 하자. 

 

 

 

위의 식 (5.1)은 (4.1)과 같은 식이다. 양변을 y로 나누고 dt를 곱해준 뒤, 적분(integral)을 해주면 식 (5.2)와 같이 되고, 적분을 수행한 결과가 (5.3)이다. 좌변의 적분 상수 C1을 우변으로 넘겨주고 C1+C2=C를 만든 다음 양변에 지수함수(exponential function)를 취해주면 (5.4)와 같이 된다. 지수법칙에 의해 e^c를 분리한 뒤 C로 치환해주면 식 (5.5)와 같이 정리할 수 있다. 이때 (5.4)의 y의 절대값에 의해 C는 +C와 -C의 두 가지 경우가 발생하게 되는데, C자체가 임의의 값이 될 수 있으므로 +C를 선택하여 (5.5)와 같이 정리할 수 있다. 이렇게 하여 (5.5)와 같이 미분방정식의 일반해(general solution)를 도출하였다. 

 

이해를 돕기위해 아래의 연습문제를 풀어보도록 하자. 

 

Ex) 당신은 연이율 2.4% 복리의 거치식 적금에 가입하였고 100만원을 예금하였다. 이 적금이 10년후, 50년 후엔 각각 얼마가 되어있을까? 

 

식 (5.5)에 이 문제를 적용하면 각 변수는 다음과 같이 설정할 수 있다. 우선 예금금 100만원은 초기값으로 C=100만원이 되고, 연이율 2.4%는 비례 상수(proportionality constant)가 되어 k=0.024가 된다. 문제에서 나타난 수치들을 식 (5.5)의 수식으로 만들어보면 다음과 같다. 

 

 

식 (6.1)은 연습 문제의 내용을 식 (4.1)의 미분방정식의 형태로 나타낸 것이다. 해(t)가 지날 때마다 내가 저축한 원금(y)이 얼마만큼씩 변하는지가 (6.1)의 좌변이 나타내는 것이고, 이 변화율은 현재의 시간 t에서의 원금 y에 이율 2.4%를 곱한 것과 같다는 것이 우변이 의미하는 것이다. 문제의 내용을 미분방정식의 형태로 정리하긴 했는데, 이를 어떻게 풀지가 막막할 것이다. 이때 활용할 수 있는 것이 식 (5.5)의 일반해(general solution)이다. 식 (6.2)는 (6.1)의 미분방정식을 일반해로 나타내어 푼 것이다. 

 

계산 결과 10년 후에는 약 127만원, 50년 후에는 약 332만원이 통장에 남게된다. 이 값은 근사값(approximate value)으로 실제 값과는 약간 차이가 있으나 거의 근접한 값을 구할 수 있다. 그러나 t가 커질 수록 오차는 벌어지게 된다. 이와 같이 미분방정식(differential equation)은 지수함수(exponential function)가 포함된 일반해(general solution)의 형태로 표현할 수 있으며 이를 통해 미분방정식의 해(solution)를 구할 수 있다. 이 관계를 잘 숙지하도록 하자. 

 

 

 

2. 미분방정식과 선형대수

 

- Differential equations and Linear algebra

 

미분방정식은 식 자체가 굉장히 복잡한 만큼 여러 가지 풀이 방법이 존재한다. 이번 챕터에서 배울 내용은 이러한 미분방정식을 선형대수의 행렬을 이용해서 푸는 방법이다. 이곳에서 다룰 내용은 상미분방정식(ODE)만 해당한다. 우선 아래의 두 개의 선형미분방정식(Linear Differential Equation)을 살펴보자. 

 

 

식 (7)은 u1과 u2에 대한 선형미분방정식을 differential form으로 각각 나타낸 것이다. 이 표기법이 익숙치 않다면 u1, u2를 x, y로 각각 생각하면 된다. u1과 u2의 미분값은 각각 서로를 포함하고 있기 때문에 이 둘은 서로 관계가 있고 따라서 하나의 벡터로 만들어서 행렬을 이용해 함께 풀어낼 수 있다. 행렬을 이용하여 풀기 위해 이들을 벡터로 정의하면 아래와 같다.  

 

 

여기서 우리가 구하고자 하는 것은 식 (7)의 두 개의 미분방정식을 행렬을 이용하여 1차 식으로 만들어 푸는 것이다. 이를 위해 u1과 u2를 원소(component)로 하는 벡터 u(t)를 식 (8.1)과 같이 정의할 수 있다. (8.2)는 u를 미분한 du/dt를 나타낸 것이고, (8.3)은  t=0일때의 초기 조건(initial condition)에서의 값을 의미한다. 

 

u와 u'을 정의했으니 이제 남은 행렬 A를 정의하면 우리가 구하고자 하는 미분방정식의 행렬 형태를 정의할 수 있다. A는 어떻게 정할 수 있을까? 이를 위해 먼저 아래와 같이 식을 나열해보자. 

 

 

식 (9)의 u1'은 (7.1)의 좌변에 표현된 u1의 미분값에 해당한다. 이 u1'이 되기 위해선 우변의 u1과 u2에 각각 얼마씩이 곱해져야하는지를 생각해보면 A의 row1을 알 수 있다. -1과 2를 각각 u1과 u2에 곱하면 u1'이 된다. 따라서 A의 row1은 -1, 2가 되고, 마찬가지로 식 (7.2)를 참고하여 row2를 구해보면 1, -2가 되는 것을 알 수 있다. 

이렇게 하여 식 (9)를 다시 써보면 아래와 같다. 

 

 

위와 같이 식 (7)의 두 개의 미분방정식을 행렬에 대한 식으로 표현하였다. 이 식을 통해 두 미분방정식에 대하여 무엇을 알 수 있을까? 우리는 식 (7)에 표현된 u1과 u2에 대한 미분방정식이 시간이 지남에 따라 어떻게 변하는지를 식 (9)를 통해 알 수 있다. A의 고유값(eigenvalue)과 고유벡터(eigenvector)를 통해서 말이다. A의 고유값과 고유벡터를 구해보자. 

 

 

식 (10.1)을 통해 행렬 A의 고유값이 각각 0과 -3임을 알았다. A가 2x2이면서 특이행렬(singular matrix)이기 때문에 두 개의 고유값중 하나가 0이 나오는 것이다. 이 특성을 알고 2x2행렬에선 trace가 모든 고유값들의 합과 같다는 성질을 이용하여 계산과정 없이 고유값을 빠르게 구할 수도 있다. 다음으로 고유벡터를 구하기 위해 (A-λI)x=0 의 해인 null space를 구해야 한다. 이때의 행렬이 특이 행렬이므로 x2가 free variable이다. 이 free variable을 1로 설정한 뒤 나머지 해를 구하면 (10.2)와 (10.3)과 같이 고유벡터를 구할 수 있다. 

 

 

- General Solution

 

식 (9)의 행렬 A에 대한 고유값(eigenvalue)과 고유벡터(eigenvector)를 구했으니 이제 이들을 활용하여 u(t)에 대한 해를 구해보도록 하자. A는 2x2크기의 행렬이므로 두 개의 고유값과 고유벡터를 가진다. 따라서 u(t)에 대한 일반해(general solution)는 두 개의 special solution, 즉 두 개의 exponential solution의 조합으로 표현할 수 있다. 아래의 식을 보자. 

 

 

식 (11)은 (9)에 나타난 u(t)의 일반해(general solution)를 나타낸 것이다. 식 (5)에서 보인것과 같이 미분방정식의 일반해를 exponential solution의 형태로 보이긴 했는데, 과연 이들이 미분방정식의 해의 형태가 될 수 있는가에 대한 의문이 생길 수 있다. 우리가 최종적으로 구하고자 하는 du/dt = Au에 exponential solution을 대입해서 이를 확인해보도록 하자. 

 

 

(12.1)에서와 같이 초기값을 고려하지 않은 순수 지수함수(pure exponential)를 u로 여겨 원래의 식에 대입하여 정리하면 (12.2)와 같이 정리할 수 있다. 좌변은 u의 pure exponential을 t로 미분한 결과이며, 따라서 exponential의 지수부인 람다(λ)가 앞으로 나와 계수로 곱해진 것이다. (12.2)의 양변을 지수함수로 나누어주면 결과적으로 (12.3)과 같이 정리할 수 있으며, 이는 앞서 공부했던 고유값/고유벡터의 정의와 같은 것을 볼 수 있다. 마찬가지로 λ2에 대한 pure exponential을 대입 해도 해를 구할 수 있다. 이를 통해 우리가 구하고자 하는 (12.1)의 행렬로 정의된 미분방정식의 해를 지수함수(exponential function)를 통해 구할 수 있음을 증명한 셈이다. 

 

식 (11)에 정의된 지수함수(exponential function)로 정의된 미분방정식의 일반해는 지난 강의(Lecture 22)에서 배웠던 대각화(diagonalization)에서의 일반해와 같은 꼴로 생각할 수 있다. 단지 t로 정의된 연속적인(continuous) 개념인지, k로 정의된 이산적인(discrete) 개념인지가 다르다. 아래의 두 식은 2x2크기 행렬의 미분방정식(differential equation)과 차분방정식(difference equation)의 일반해를 각각 나타낸다.  

 

 

이제 지수함수와 벡터는 알겠는데, 각 항의 앞에 곱해진 상수 c1과 c2는 어떻게 알 수 있을까? 이들 상수는 초기 조건(initial condition)으로부터 계산할 수 있다. 이는 앞부분에 미분방정식의 개념을 설명할 때 배웠던 식 (4)의 일반해의 초기값(initial value) C에 해당하는 값이다. 식 (8.3)의 초기값으로부터 c1과 c2를 구해보도록 하자.  

 

 

우선 (14.1)은 앞서 구한 고유값과 고유벡터를 대입하여 일반해의 식을 정리한 것이다. 고유값이 각각 λ1=0, λ2=-3이므로 첫 번째 지수 함수는 1이 되고, 두 번째 지수함수의 지수부는 -3t가 되어 (14.1)과 같이 정리할 수 있다. 이때 t=0일 때의 초기값 u가 u(0)=[1 0]T이므로 t에 0을 대입하여 정리하면 (14.2)와 같이 정리할 수 있다. c1과 c2는 각각 column vector에 곱해진 계수가 되기 때문에 이 column vector(고유 벡터들)들을 행렬로 정리하면 (14.3)과 같이 나타낼 수 있으며, 이때의 행렬을 고유벡터 행렬(eigenvector matrix) S라고 한다. S를 가우스 소거(Gauss elimination)를 통해 c1과 c2를 계산하면 각각 c1=1/3과 c2=-1/3이 된다. (가우스 소거법은 Lecture 2참조)

 

고유벡터 행렬 S를 통해 우리는 해(solution)의 계수(coefficient) c1, c2, ...를 구할 수 있으며, 이 계수들은 해에서 각 순수 지수함수(pure exponential)들이 얼마만큼인지를 나타낸다. 이 계수들은 t=0일 때의 u인 u(0)로부터 구할 수 있으며, 이렇게 구한 계수들은 t가 무한대로 갔을 때 경우에 따라 극한값이 되기도 한다. 

 

c1과 c2를 구했으니 이들을 일반해 식에 대입하여 정리하면 아래와 같다. 

 

 

식 (15)가 (9)의 미분방정식의 일반해(general solution)이다. 일반해를 도출하기까지의 과정을 살펴보면 먼저 미분방정식을 행렬의 형태로 만든 다음, 초기 조건(initial condition)에 따른 초기 값을 구하고 행렬의 고유값과 고유벡터를 구한 뒤, 고유벡터행렬을 통해 각 term의 계수값들(c1, c2, ...)을 구한다. 행렬 A의 고유값, 고유벡터, 계수 등 모든 값들을 구하여 일반해 형태의 식에 대입하고 정리하여 최종적인 해를 구할 수 있다. 

 

식 (15)와 같이 일반해를 구했기 때문에 시간이 지남에 따라 미분방정식의 값이 어떻게 변하는지를 쉽게 알아낼 수 있다. 시간에 대한 변수인 t값만 대입해주면 해당 시간에서의 값을 알 수 있다. 그런데 우리가 세운 미분방정식이 만약 어떤 기계의 진동에 따른 기계 수명 예측과 같은 문제라고 한다면, 경영진 입장에서는 당연히 아주 오랜 시간이 지났을 경우 기계의 기대 수명이 궁금할 것이다. 즉 우리가 여기서 알고자 하는 것은 시간이 무한대로 지났을 때 미분방정식의 값이다. 이는 앞서 구했던 미분방정식 행렬의 고유값과 관련이 깊다. 다음 파트에서 고유값과 미분방정식의 무한대 시간에서의 해(infinite solution)의 관계에 대해 알아보자. 

 

 

- Stability

 

우리는 지난 대각화 강의에서 고유값을 통해 해당 차분방정식의 특성을 파악할 수 있음을 언급했다. 마찬가지로 이번 강의에서 다뤘던 미분방정식 행렬의 경우에도 고유값을 통해서 관련 특성을 파악할 수 있다. 즉 식 (9)의 A의 고유값을 통해 우리는 행렬 A에 내포된 미분방정식들이 시간이 지남에 따라 어떻게 변하는지를 예측할 수 있다. 그리고 궁극적으로 시간이 무한대로 지났을 때 그 극값이 어떻게 되는 지를 알 수 있다. 

 

식 (10)에서 계산한 첫 번째 고유값은 0이고, 이 고유값은 식 (11)에서 일반해의 첫 번째 term인 c1*exp(0*t)*x1의 지수 함수의 지수부가 된다. 따라서 지수함수는 t가 아무리 증가해도 0의 거듭제곱 꼴이 되기때문에 항상 1이 되므로 식 (15)에서 첫 번째 term에서 지수함수 부분이 사라진 것이다. 이는 결국 미분방정식이 시간이 지날수록 정상상태(steady state)로 수렴함을 의미한다. 

 

두 번째 고유값은 -3이고 식 (11)에서 c2*exp(-3*t)*x2의 꼴이 되기 때문에 t가 커질수록 지수함수의 지수부는 작아지고 점점 작은 값이 된다. 따라서 두 번째 term의 식 자체가 시간이 증가할수록 굉장히 작은 값이 된다. 결과적으로 우리가 구한 일반해 (15)의 극한값(t가 무한대로 갔을 때의 값)은 아래와 같다. 

 

 

 

식 (16.2)가 극한 값이며, 실제 값은 [0.666... , 0.333...]T이다. 식 (7.1)의 미분방정식의 t가 무한대일때의 극한값은 0.666..., 식 (7.2)의 극한값은 0.333...이다. 이들은 2개의 미분방정식 (7.1), (7.2)가 연립된 형태인 시스템 (9)에 대한 해 벡터(solution vector)의 극한값 u(t=)=[0.666..., 0.333...]T이다. 

 

미분방정식은 대수방정식(Algebraic equation)처럼 어떤 값이나 값들의 집합의 형태로 해가 나오는 것이 아니다. 미분방정식의 해는 (15)와 같이 식의 형태로 나오며, MATLAB과 같은 프로그램을 통해 시간이 지남에 따라, 혹은 종속 변수(dependent variable)가 변함에 따라 해가 어떤 식으로 변하는지를 알 수 있다. 아래 그래프는 (15)의 해가 시간이 지남에 따라 변하는 과정을 나타낸 것이다. 

 

 

Fig. 1 미분방정식의 해가 시간이 지남에 따라 변화하는 과정

 

 

Fig. 1에서 빨간색 선과 파란색 선이 각각 일반해의 u1과 u2를 의미한다. t=0에서부터 시간이 지남에 따라 특정 값에 수렴하는 것을 볼 수 있는데, 최종적으로 수렴하는 값이 u1=0.666..., u2=0.333... 임을 알 수 있다. 이는 (16)에서 구한 극값과 일치한다. uf(first)라고 표시된 벡터는 일반해에서 첫 번째 고유값인 0이 속해있는 식을 의미한다. 녹색과 자홍색 그래프가 uf의 u1과 u2에 각각 해당하는데, 시간이 지나도 어떤 특정값에서 변화가 없음을 확인할 수 있다. 이는 고유값이 0이기 때문에 지수함수 부분이 항상 1이 되기 때문이다. 

 

us(second)라고 표기된 벡터는 일반해에서 두 번째 고유값인 -3이 속해있는 식을 의미한다. 지수부의 고유값이 음수이기 때문에 시간이 지날수록 값이 작아지는 것을 볼 수 있다. 하늘색과 검은색 그래프가 초기 값에서 시간이 지날수록 0으로 수렴하는 것을 볼 수 있다. 

 

이와 같이 미분방정식을 행렬을 이용하여 풀 때 행렬의 고유값에 의해서 우리는 연립미분방정식의 해가 시간이 지남에 따라 정상상태(steady state)수렴할 수도, 값이 0이 되어 사라질 수도 혹은 값이 무한대로 커져서 발산할수도 있다. 그렇다면 고유값이 어떤 상태일 때 발산 혹은 수렴이 될까? 아래의 정리를 살펴보자. 

 

(1) Stability  (2) Steady State  (3) Divergence 



 

 

미분방정식의 해는 크게 세 가지 상태가 될 수 있는데, 위의 표에 나타난 것과 같이 0으로 수렴하는 안정화 상태인 stability, 어떤 특정 값으로 수렴하는 정상 상태(steady state) 그리고 아예 값 자체가 무한정 커지는 divergence상태가 있다. 

 

만약 모든 고유값이 0보다 작다면 stability 상태가 된다. 이때는 해(solution) u(t)가 시간이 지날수록 0으로 수렴하게 되는데, 모든 고유값이 0보다 작은 음수인 경우 모든 지수함수의 값들이 점점 작아지기 때문이다. 여기서 고유값을 실수부(real part)라고 지정한 것은 해가 허수부(imaginary part)를 포함할 수도 있기 때문이다. 아래의 예를 보자. 

 

 

식 (17.1)은 허수부를 포함한 지수함수이다. 지수법칙에 의해 실수부와 허수부를 나눌 수 있는데, 이때 허수부의 크기는 오일러 공식에 의해 항상 1이 되기 때문에 t가 커진다고 해도 해에 실질적인 영향은 주지 않는다. 허수부는 오일러 공식 위키(https://en.wikipedia.org/wiki/Euler%27s_formula)에서 볼 수 있듯이 t가 증가해도 그저 복소 평면에서 단위 원(unit circle)을 반복해서 돌 뿐이다. 결국 Stability는 고유값의 실수부(real part)에 의해서만 결정된다

Stability를 만족시키는 조건은 2x2 크기의 행렬에서는 trace와 행렬식(determinant)로 쉽게 구할 수 있다. 2x2행렬에선 trace(A)=λ1+λ2, det(A)=λ1*λ2 임을 이용하여 trace < 0와 det > 0를 만족한다면 두 개의 고유값이 모두 음수이므로 Stability상태가 됨을 쉽게 알 수 있다. 

 

두 번째로 정상 상태(Steady state)의 경우는 해가 어떤 특정 값으로 수렴하는 경우를 의미하는데, (15), (16)에서 구한 해가 바로 이 경우 이다. 즉 어떤 하나의 고유값이 0이고 나머지 고유값이 0보다 작은 경우에 해당하는데, (16)의 경우에 첫 번째 고유값이 0이고, 두 번째 고유값이 0보다 작은 -3이므로 시간이 지날수록 두 번째 term은 0으로 수렴하여 사라지고 첫 번째 term의 상수인 c1=1/3과 첫 번째 고유값 x1=[2 1]T의 곱의 값으로 수렴하게 된다. Fig. 1의 빨간색과 파란색 그래프가 바로 이 특정 값으로 수렴하는 모습을 보여준다. 

 

마지막의 경우는 어떤 고유값이 하나라도 0보다 큰 경우 시간이 지날 수록 무한대로 커져서 발산하게 된다. 

 

종합해보면 복소수까지 고려했을 때 오직 실수부의 고유값만이 미분방정식의 해에 영향을 미치고, 고유값(eigenvalue)이 0보다 크고 작음, 혹은 같음에 따라 해의 상태가 결정된다. 

 

아래 그림은 Fig. 1을 plot하기 위한 MATLAB 소스코드이다. 

 

 

 

3. 마치며(continue)

 

이번 강의에선 미분방정식의 간단한 소개와 함께 선형대수를 이용하여 미분방정식을 푸는 방법을 배웠다. 기본적으로 미분방정식에 대한 식을 세운 뒤 고유값과 고유벡터를 구하여 일반해에 대한 식을 세운 다음 초기 조건으로부터 일반해의 계수를 구하는 순서로 해를 구하게 된다. 해는 지수함수의 고유값의 거듭제곱꼴로 나타나는데, 고유값이 0보다 크거나 같음 또는 작음에 따라 해의 상태가 결정됨을 배웠다. 다음 포스팅에서는 이번 강의에 이어서 미분방정식과 선형대수의 나머지 이야기를 하도록 하겠다. 

 

이번 포스팅에서 다룰 내용은 바로 행렬의 대각화(Diagonalization)이다. 행렬의 대각화는 지난 시간에 배운 고유값(eigenvalue)과 고유벡터(eigenvector)를 활용하기 위한 하나의 방법이라고 할 수 있으며, 다른 말로는 고유값분해(Eigendecomposition)라고도 불린다. 또한 행렬의 대각화를 통해 LU 분해, QR분해와 같이 행렬을 고유값과 고유벡터로 구성된 부분 행렬들로 분해할 수 있으며, 이는 어떤 반복적인 선형방정식을 풀 때 굉장히 유용한 특성을 가지고 있다. 대각화에 대해 공부해보자. 

 

 

1. 행렬의 대각화(Diagonalization)

 

- Diagonalizing a matrix

 

지난 시간에 우리는 고유값(eigenvalue)과 고유벡터(eigenvector)에 대한 내용을 배웠다. 일단 관련 식을 다시 써보자. 

 

 

식 (1)을 다시 한 번 설명해보면 람다(lambda, λ)는 고유값을, x는 고유벡터를 각각 나타내며, 행렬 A에 의해 선형 변환(Linear transformation)을 시켜도 변환 전과 후가 평행한(parallel) 벡터를 고유벡터, 그리고 길이 변화량의 정도를 나타낸 것이 고유값이다. 이 고유값/고유벡터를 알면 행렬 A의 중요한 특성과 정보를 알 수 있으며 다양한 곳에 응용할 수 있다. 그 중 하나가 바로 행렬의 대각화(diagonalization)이다. 그렇다면 고유값과 고유벡터를 이용하여 행렬 A를 어떻게 대각화 할 수 있는 걸까? 

 

우선 식 (1)은 하나의 고유벡터와 고유값에 대한 식이다. 그러나 일부 경우를 제외하곤 대부분의 nxn크기의 행렬은 n개의 고유값과 고유벡터를 갖는다. 어떤 행렬 A의 n개의 고유값과 고유벡터를 찾은 뒤엔 아래의 식에 따라 대각화를 수행하면 된다. 

 

 

식 (2)는 대각화에 대한 식이며, 이번 포스팅의 핵심적인 식이라고 보면 된다. 식 (2)에서 A는 원래의 행렬을 의미하고, 행렬 S는 A의 고유벡터들을 column vector형태로 차례로 끼워 넣은 nxn크기의 고유벡터 행렬이다. 이때 A의 앞에는 S의 역행렬이, A의 뒤에는 S가 각각 곱해진다. 이렇게 좌변을 계산하면 우변의 Λ (대문자 람다)행렬이 만들어지는데, 이 람다 행렬(Λ)은 대각 행렬(diagonal matrix)이고 각각의 대각 원소들은 고유값들로 차례로 채워져있다. 이제 식 (2)가 대충 뭘 나타내는지는 알겠는데, 도대체 저게 어떻게 만들어졌고, 무슨 의미가 있다는 걸까? 이제부터 차근차근 알아보도록 하자. 

 

그 전에 먼저 식 (2)에서 눈여겨봐야할 것이 있다. 앞에서 S는 A의 고유벡터들로 이루어진 고유벡터행렬(eigenvector matrix) 이라고 배웠다. 그런데 식 (2)의 A앞에 S의 역행렬(Inverse matrix)이 곱해져 있는 것을 볼 수 있다. 이것이 의미하는 것은 무엇일까? 바로 S가 역행렬을 가질 수 있어야 하고, 이는 S가 특이 행렬(singular matrix)이 아니어야 함을 의미하며, 결국 A의 고유벡터들이 n개의 독립(independent)인 벡터를 가져야 한다는 것이다

 

일단 A가 n개의 독립인 고유벡터를 가진다고 가정할 때, 이 고유벡터들을 column vector의 형태로 차례로 붙여서 S를 만들었다고 해보자. 이때 A와 S를 곱하면 어떤 일이 벌어질까? 아래의 수식을 통해 확인해보자. 

 

 

식 (3)은 행렬 A와 그의 고유벡터(eigenvector)들로 만들어진 행렬 S와의 곱을 나타낸다. S에서 x1, x2, ... xn은 각 고유벡터들을 column 형태로 삽입한 것이다. A와 S를 곱할 때 지난 강의(행렬 곱셈) Lecture 3에서와 같이 column-wise로 생각할 수 있다. 즉 A와 S의 첫 번째 column vector x1과 곱해져서 Ax1이라는 하나의 column vector가 만들어지고, 차례로 Ax2, ... Axn이 계산되어 결과적으로 이들 Ax1, Ax2, 등의 column vector들로 이루어진 동일한 크기의 행렬이 만들어진다. 

 

그런데 가만히 보면 Ax1은 원래의 행렬 A와 그의 첫 번째 고유벡터와의 곱이다. 여기서 우리는 식 (1)을 떠올릴 수 있다. Ax1은  λ1x1이다. 따라서 (3.2)와 같이 쓸 수 있다. 결국 AS의 곱을 (3.2)와 같이 각 column vector들을 고유값과 고유벡터들의 곱으로 표현할 수 있으며, 이는 식 (1)의 개념으로부터 나온 것이다. 

 

최초에 AS의 곱을 (3.2)와 같이 표현하였다. 여기서 고유값들을 고유벡터로부터 한 번 더 분리하여 나타낼 수 있다. 이때 중요한 것은 (3.2)를 S와 고유값들과의 곱으로 분리할 때 행렬끼리의 곱으로 나타내고 싶다는 것이다. 어떻게 할 수 있을까? 다시 한 번 column-wise의 행렬곱을 생각해보자. S를 고유벡터행렬, 여기에 곱해질 어떤 행렬을 가령 V라고 생각해 봤을 때, Sv1이 (3.2)의 첫 번째 column vector인 λ1x1이 되어야한다. 이렇게 되기 위해선 v1의 첫 번째 component만 λ1이고 나머지 component들은 모두 0이어야 한다.  λ2x2가 되기 위해선 v2의 두 번째 component만 λ2이고 나머지는 0이어야 한다. 이런식으로 만든 행렬이 바로 (3.3)의 오른쪽 고유값(eigenvalue)들로 이루어진 행렬이고, 대각 행렬(diagonal matrix)의 형태이다. 우린 이 고유값들로 이루어진 행렬을 대문자 람다(capital lambda)를 써서 Λ라고 표현한다. (※ 일반적으로 대문자는 행렬을, 소문자 볼드체는 벡터를, 소문자는 벡터의 원소나 스칼라 값을 나타냄)

 

결과적으로 (3.1)에 표현된 최초의 AS의 행렬곱셈을 (3.3)과 같이 SΛ로 표현할 수 있으며 AS=SΛ이다. 이렇게 표현할 수 있는 근간에는 우리가 지금까지 공부해왔던 고유값/고유벡터에 대한 식 (1)이 있다. 즉 식 (3)은 식 (1)의 확장판이라고 생각하면 된다. 식 (3)에서 얻은 결론을 다시 써보면 아래의 식과 같다. 

 

 

식 (4)가 의미하는 것은 어떤 행렬 A에 고유벡터행렬 S를 곱하면, 고유벡터행렬에 고유값행렬을 곱한 것과 같다, 즉 식 (1)을 모든 고유값/고유벡터에 대해서 한번에 정리한 식이라고 할 수 있다. 결국 식 (1)로부터 (4)를 유도한 셈이다. 식 (4)로부터 우리는 다양한 식을 유도할 수 있다. 양변의 왼쪽에 S의 역행렬(inverse matrix)을 곱해보자. 

 

 

식 (4)의 양변의 왼쪽에 S의 역행렬을 곱했더니 우리가 처음에 핵심적인 식이라고 배웠던 식 (2)가 만들어졌다. 이 과정이 임의의 정방행렬 A를 대각화(diagonalization)하는 과정이다. 여기서 중요한 것은 양변에 S 역행렬을 곱했는데, 이것은 S가 역행렬을 가질 수 있어야 한다는 조건을 만족해야 하며, A가 n개의 독립인 고유벡터를 가져야 함을 의미한다. 지난 강의 Lecture 21-(2)의 마지막에 n개의 독립인 고유벡터(eigenvector)를 가지지 못하는 행렬(triangular matrix)에 대해서 알아봤다. 그러나 이러한 경우는 일부이며, 대부분의 행렬은 n개의 독립인 벡터를 가지기 때문에 여기서는 고유벡터행렬(eigenvector matrix)이 역행렬을 가질 수 있음을 가정하겠다. 

 

 

- Factorization of a matrix

 

우리는 식 (4)로부터 행렬 A를 대각화(diagonalization)하여 고유값들로만 이루어진 식 (2)와 같은 대각 행렬을 만들었다. 그런데 이런 고유값행렬(eigenvalue matrix)말고도 다른 형태의 행렬을 만들 수 있다. 식 (4)에서 양변의 좌측에 S의 역행렬을 곱했다면, 이번에는 양변의 우측에 역행렬을 곱하는 것이다. 아래 식을 보자. 

 

 

식 (5.1)로부터 양변의 우측에 S의 역행렬을 곱했더니 식 (5.2)와 같은 식을 만들어냈다. 즉 A를 고유벡터행렬(eigenvector matrix)인 S와 고유값행렬(eigenvalue matrix)인 Λ, 그리고 고유벡터행렬의 역행렬의 곱으로 정의한 것이다. 이는 결국 A를 이전에 공부했던 LU factorization(가우스 소거로부터 행렬을 분해하는 방법, Lecture 4), QR decomposition(그람 슈미트 정규 직교화 방법으로 행렬을 분해하는 방법, Lecture 17-(2))과 같이 행렬을 분해하는 하나의 방법이며, 고유값과 고유벡터들의 행렬의 조합으로 행렬을 인수 분해(factorization)하는 방법이다. 조합되는 방법이 S와 대각 행렬, 그리고 다시 S의 역행렬이 반복됨을 주목하자. 

 

 

- Useful attribute

 

이제 이 행렬이 얼마나 유용한 특성을 가지는지 알아보도록 하자. 우리가 어떤 문제를 풀때 선형대수를 이용하여 선형연립방정식을 세우고, 이를 반복적으로 곱해서 계산해야 할 때가 있다. 이런 경우 A를 제곱, 3제곱, ... n제곱 과 같은 식으로 반복해서 곱해야할 때가 있는데, 이렇게 반복하여 자기 자신을 곱해야할 때 대각화에 의한 행렬 분해가 굉장히 유용한 특성을 가진다. 먼저 이 분해가 고유값과 고유벡터 식으로부터 유도된 만큼 A를 제곱했을 때 고유값/고유벡터 식 (1)이 어떻게 변하는지 알아보자. 

 

 

식 (6.1)은 고유값/고유벡터 식이다. 여기서 A를 제곱했다는 것은 식 (6.1)의 양변의 좌측에 각각 행렬 A를 곱한 것과 같다. 이렇게 양변에 A를 곱하면 식 (6.2)가 만들어지는데, (6.2)에서 우변의 람다는 상수이기 때문에 앞으로 빼서 정리했다. 이때 (6.2)의 우변에서 Ax는 λx와 같기 때문에 바꿔서 써주면 λλx가 되고, 최종적으로 (6.3)과 같이 정리할 수 있다. 결국 고유값/고유벡터 식에서 A를 제곱하면 우변의 고유값(eigenvalue)인 람다(λ)가 제곱이 된다. 이때 고유벡터는 A의 제곱과는 무관하게 변함없이 그대로이다. 

 

식 (6)에서 단일 고유값/고유벡터 식의 A제곱에 대한 내용을 살펴봤으니 이번에는 고유값행렬과 고유벡터행렬에 대한 식 (5)를 가지고 A의 제곱을 하면 어떤 결과가 나오는지 살펴보자. 식 (5.2)의 A를 제곱한 결과는 아래 식과 같다. 

 

 

식 (7.1)의 A를 제곱하면 식 (7.2)의 밑줄 친 부분이 단위행렬(identity matrix)이 되어 소거(cancel)되고, 결국 (7.3)과 같이 고유값행렬인 Λ의 제곱이 된다. 이는 단일 고유값/고유벡터에 대한 A제곱식 (6)과 같은 꼴이며 다만 그 형태가 n개의 고유값/고유벡터에 대한 행렬의 형식인 것이다. (6)과 마찬가지로 A를 제곱해도 고유벡터행렬(eigenvector matrix)은 그대로 유지되며 고유값행렬(eigenvalue)만 제곱의 형태로 나타난다. 

 

그런데 제곱이 아니라 세제곱, 네제곱, k제곱을 한다면 어떻게 될까? 실제로 계산을 해보면 세제곱을 하는 경우 역시 중간의 고유벡터행렬들은 단위행렬로써 소거되고, 우변의 고유값행렬만 세제곱이 된다. 결과적으로 A가 k번의 거듭제곱(power)를 하는 경우, 우변에서 고유벡터행렬(eigenvector matrix)은 그대로 유지되고 고유값행렬(eigenvalue matrix)만이 k번의 거듭제곱이 일어난다는 일반적인 식을 도출할 수 있다. 아래의 식 (8)은 이러한 일반식을 나타낸다. 

 

 

식 (8)을 잘 곱씹어보면 다음과 같은 통찰력을 얻을 수 있다. 어떤 행렬 A의 고유값(eigenvalue)과 고유벡터(eigenvector)는 A의 k번의 거듭제곱을 이해하는데에 아주 좋은 방법을 제공해준다. 예를 들어 LU분해나 QR분해 등의 방법을 이용하여 행렬을 인수분해 했을 경우, A를 1000제곱하면 LU*LU*...를 1000번 수행한 행렬이나 QR*QR*...를 1000번 수행한 행렬을 가지고 A의 1000제곱이 어떤 결과를 보일지를 예측하거나 A가 어떤 행렬인지를 분석해야 한다. 이들 분해 방법으로는 거의 불가능하다. 그러나 고유값행렬과 고유벡터행렬로 분해한 식 (8)의 방법은 A의 1000제곱에 대해서 단지 고유값행렬(eigenvalue matrix)의 1000제곱만 분석하면 된다. 이마저도 고유값행렬의 1000제곱을 할 필요도 없이 각 고유값들의 제곱을 계산해도 된다. 

 

 

식 (9)는 고유값행렬의 제곱을 나타내며, 그 결과는 대각 원소들인 각 고유값들이 자기자신을 제곱한 것과 같은 것을 볼 수 있다. 따라서 A의 거듭제곱에 대한 계산이 훨씬 쉽고 효율적으로 진행된다. 

 

또한 행렬 A를 고유값행렬과 고유벡터행렬로 대각화(diagonolization)하면 A가 어떤 특성을 가지는지, 혹은 k번의 거듭제곱을 했을 때, 어떤 결과가 나올지를 대략 유추해 볼수도 있다. 가령 아래와 같은 조건을 갖출 경우, 우리는 행렬 A를 안정 행렬(stable matrix)라고 할 수 있다. 

 

 

식 (10)의 조건에 따르면 어떤 행렬 A의 모든 고유값의 절대값이 1보다 작을 때 A의 거듭제곱(power)을 지속해 나갈 수록 A는 0으로 수렴한다. 이때 A는 stable하다 라고 할 수 있다. 그 아래 예를 보면 행렬 A의 고유값은 각각 λ1=0.4172, λ2=0.9828인데 둘 다 절대값이 1보다 작다. 이때 A를 거듭제곱한 결과가 그 아래 나와있는데, A의 10제곱, 50제곱, 100제곱을 거듭 할수록 A의 원소의 값들이 점점 작아지는 것을 볼 수 있다. 결국 어떤 행렬의 고유값(eigenvalues)들의 상태를 보고 해당 행렬의 특성이나 거듭제곱을 무한대로 반복했을 때의 결과를 끝까지 해보지 않아도 예측할 수 있는 것이다. 

 

 

2. 대각화 가능한 행렬(Diagonalizable matrices)

 

- Which matrices are diagonalizable?

 

지금까지 우리는 어떤 행렬 A를 대각화(diagonolization)과정을 통해 고유값행렬(eigenvalue matrix)과 고유벡터행렬(eigenvalue matrix)로 분해한다면 A의 거듭제곱을 계산하는데에 있어 굉장히 유용한 특성이 있음을 배웠다. 다시 한 번 강조하지만 대각화(diagonalization)를 위해 가장 중요한것은 A가 n개의 독립(independent)인 고유값과 고유벡터를 가지고 있어야 이 모든 것이 성립한다는 것이다. 그렇다면 어떤 행렬들이 대각화가 가능할까? 바로 아래의 정의를 만족시키는 행렬들이 대각화 가능한(diagonalizable)행렬들이다. 

 

 



대각화 가능한(diagonalizable) 행렬의 조건:
    • 어떤 행렬 A의 고유값(lambda, λ)들이 전부 서로 다른 값을 가진다면, 즉 반복되는 고유값이 없다면, A는 반드시 n개의 독립인 고유벡터(eigenvectors)를 가지며 대각화가 가능하다(diagonalizable)

 

위의 노란 박스안에 정의된 조건을 만족하는 행렬은 대각화가 가능한 행렬이다. 위의 정의에 대한 증명이 궁금하다면 여기를 참조하기 바란다. 혹시나 위의 정의가 맞는지 의심이 들 수도 있기에 MATLAB을 이용하여 실험을 해보았다. MATLAB의 랜덤 함수 rand()를 이용하여 20x20크기의 랜덤값을 가진 행렬을 만들고, 고유값행렬(eigenvalue matrix)과 고유벡터행렬(eigenvector matrix)을 eig()함수를 이용하여 구한 다음 고유값행렬에서 같은 고유값이 존재하는지 확인한다. 그 다음으로 n개의 독립인 고유벡터를 가지는지 확인하기 위해 고유벡터행렬의 rank를 계산하여 n보다 작은지를 검사하여 그 횟수를 세는 간단한 프로그램이다. for문을 이용하여 1000번 반복하였다. 아래 그림은 MATLAB코드와 그 결과를 캡쳐한 화면이다. 

 

 

 

 

Fig. 1 MATLAB코드와 n개의 서로 다른 고유값 실험

 

실험 결과 1000번의 반복에도 단 한번의 카운트가 되지 않았음을 알 수 있다. 이는 결국 서로 다른 n개의 고유값을 가지면 n개의 독립(independent)인 고유벡터가 존재함을 실험적으로 보인 것이다. 물론 완벽한 증명법은 아니겠지만, 이미 수학적으로 증명된 정의를 실험을 통해 한 번 더 확인하는 과정이라고 생각해도 좋을것이다. 

 

 

- When A has repeated eigenvalues (positive case)

 

A가 만약 반복되는 고유값을 가진다면 무조건 n개의 독립인 벡터를 가질 수 없는 걸까? 반드시 그렇지는 않다. 다만 반복되는 고유값을 가진 경우, 좀 더 면밀히 살펴봐야한다. 즉 아래의 정의로 설명할 수 있다. 

 



    • A가 반복되는 고유값(eigenvalues, lambdas)을 가진 경우, A는 n개의 독립(independent)인 고유벡터(eigenvectors)를 가질 수도, 혹은 가지지 않을 수도 있다.  

 

 

그렇다면 어떤 행렬이 반복되는 고유값을 가져도 n개의 독립인 고유벡터를 가질까? 한 예가 바로 단위 행렬(identity matrix)이다. 아래 단위 행렬의 고유값/고유벡터 예를 보자. 

 

 

 

 

식 (11)에서 2x2크기의 단위 행렬의 2개의 고유값은 1로 반복되었다. 그런데 고유값이 반복되었음에도 불구하고 고유벡터는 2개의 독립인 고유벡터가 나왔다. 이는 단위행렬이 어떤 벡터를 곱하여 변형시켜도 자기 자신이 나오게 만드는 특성을 고려해보면 당연한 것이다. 이 단위 행렬을 식 (2)와 같이 대각화(diagonalization)를 하면 어떻게 되는지 살펴보자. 

 

 

 

식 (12)와 같이 대각화에 대한 식에도 성립하는 것을 알 수 있다. 사실 식 (11)에서 계산한 고유벡터 말고도 단위행렬에 대한 2차원 공간의 무수히 많은 고유벡터가 존재한다. 2x2 이외에도 nxn크기의 단위행렬의 경우에도 마찬가지로 모두 동일한 고유값을 가지지만, n개의 독립인 고유벡터가 존재한다. 결과적으로 단위 행렬의 경우 동일한 고유값을 가졌음에도 n개의 독립인 고유벡터가 존재하는 case이다. 

 

또 한 가지 예를 살펴보자. 두 번째로 살펴볼 예는 지난 강의 Lecture 21-(2)에서 살펴봤던 180도에 대한 회전행렬(rotation matrix)이다. 자세한 관련 내용은 전 강의를 참고하도록 하고, 식을 다시 써보면 다음과 같다. 

 

 

단위 행렬과 마찬가지로 고유값은 -1로 동일하게 두 번 반복됐지만, 고유벡터는 2개의 독립인 벡터가 존재함을 볼 수 있다. 지난 강의에서 봤겠지만 180도를 회전시키는 행렬의 경우 변환 전과 후의 벡터가 비록 방향은 완전히 반대로 바뀌지만 평행(parallel)하기 때문에 역시 무수히 많은 고유벡터들이 존재하게 된다. 

 

 

- When A has repeated eigenvalues (negative case)

 

이번엔 고유값이 반복되었을 때 n개의 독립인 고유벡터가 존재하지 않는 경우를 살펴보자. 사실 이 역시 지난 강의 Lecture 21-(2)에서 이미 본 적이 있는 행렬이다. 바로 삼각행렬(triangular matrix)이다. 아래 식을 보자. 

 

 

 

식 (14)의 행렬 A는 삼각행렬의 형태를 띄고 있으며, 고유값은 전부 2로 같다. 지난 강의에서 배웠듯이 삼각행렬은 대각 원소들이 곧 고유값이 됨을 기억하자. det(A-λI)를 계산할 때 아래쪽 원소들 값이 모두 0이기 때문에 반대편(cross) 대각 원소들은 방정식에 포함되지 않는다. 따라서 대각 원소들만이 방정식 계산에 포함되기 때문에 대각 원소들 값이 곧 고유값이 된다. 

 

고유벡터를 구해보면 식 (14)의 아래와 같이 되는데, 보다시피 한 개의 고유벡터만이 존재한다. (A-λI)의 null space인 고유벡터를 구할 때, 식을 만족시키기 위해서는 x2에는 어떤 값도 들어가서는 안된다. 따라서 무조건 0이 되어야 하고 x1에만 임의의 값을 넣을 수 있다. 그런데 이렇게 되면 x는 1차원 공간에 존재하게 되고, 결국 두 고유벡터가 같은 부분 공간에 존재하게 된다. 따라서 2x2행렬에서 단 1개의 고유벡터만 존재하게 된다. 이렇게 되면 n개의 독립인 고유벡터가 존재하지 않기 때문에 삼각행렬의 경우 대각화(diagonolization)가 불가능한 경우다. 

 

 

정리해보면 어떤 행렬 A가 n개의 서로 다른 고유값을 가지는 경우엔 A는 반드시 서로 다른 독립인 고유 벡터를 가지며, A는 대각화가 가능하다. 

반면 A가 어떤 반복되는 고유값을 가지는 경우엔 독립인 고유벡터를 가질 수도 있지만(단위 행렬, 회전 행렬), 그렇지 않을 경우도 존재한다(삼각 행렬)

 

 

 

3. 대각화와 차분방정식(Diagonalization and Difference Equation)

 

- Difference Equation

 

앞서 배운 대각화를 이용하면 방정식(equation)을 손쉽게 풀 수 있다. 지금껏 공부한 선형대수(Linear Algebra)가 선형연립방정식을 손쉽게 풀기 위함인데 갑자기 방정식이라니? 그냥 방정식이 아니라 우리나라말로 차분 방정식, 혹은 계차방정식이라 불리는 Difference equation이다. 계차방정식(difference equation)은 시간이 지남에 따라 상태가 변화하는 문제를 방정식으로 만들어놓은 것이다. 이를테면 바이러스가 1초마다 자기 자신을 둘로 분열시킨다고 했을 때, 처음 1마리가 1000초 후엔 몇 마리가 되어 있는가? 와 같은 문제들 말이다. 즉 방정식(equation)시간(time)과 그에 따른 미분(derivative)의 개념이 들어가 있는 것이다. 

 

여기서 어떤 사람은 "그렇다면 계차방정식이 결국 미분방정식(differential equation) 아닌가?"라고 생각할지도 모르겠다. 이 둘은 비슷해보이지만 약간 다르다. 둘 다 미분(derivative)이라는 개념이 들어가 있지만, 시간을 어떻게 보느냐에 따라 달라진다계차방정식(difference equation)은 시간을 정수단위로 끊어서 생각하는, 즉 이산적인(discrete)개념이고, 미분방정식(differential equation)은 시간을 끊어지지 않고 쭈욱 이어지는 개념, 즉 연속적인(continuous) 개념으로 생각하는 것이다. 아래 식은 똑같은 개념을 각각 계차방정식과 미분방정식의 방법으로 표현한 것이다. 

 

 

식 (15)에서 계차방정식(Difference equation)의 경우 a의 시간에 따른 변화를 n으로 표현하였다. n에 들어갈 수 있는 숫자는 오직 정수뿐이며 소수 등의 다른 숫자는 들어갈 수 없다. 계차방정식이 말하고자 하는 것은 a가 n일때를 기준으로 a의 이전 상태(n-1)는 a의 바로 앞의 상태(n+1)에서 현재 상태(n)를 뺀 값이 되며, 이것이 계속 반복된다는 것이다. n이 +1씩 변화할 때 그에 따라 a가 어떻게 변화는지를 보는 것, 즉 미분(derivative)을 이산적(discrete)으로 나타낸 것이다. n는 보통 시간이나 시퀀스(sequence)를 나타내는데 결국에는 시간(time)의 개념이 들어갈 수밖에 없다. 

 

반면 식 (15)의 우측의 미분방정식(differential equation)을 보면, 여기서 dx는 delta x를 의미하고 x의 변화를 의미한다. dy도 마찬가지로 delta y를 의미하는데, 이때 변화량을 의미하는 delta가 나타내는 그 "변화량"은 무한대로 작은 값을 의미한다. 즉 상상하기 힘들 정도로 아주 찰나의 순간을 의미하며, 결국 이 변화량은 끊어지지 않고 연속적인(continuous) 개념이 된다. dx는 시간일수도 있고 혹은 다른 물리량을 의미할 수도 있다. 여기서는 일단 시간으로 생각해보자. dy는 dx, 즉 x가 아주 찰나의 순간의 변화가 발생했을 때 y가 얼마만큼 변하는지를 나타내며, dy/dx의 나눗셈을 했기 때문에 결국 x가 변화했을 때 y가 얼만큼 변하는지에 대한 비율(ratio)을 의미하게 된다. 계차방정식(difference equation)과의 관계를 살펴보면 x는 n과 관련이 있고, y는 a와 관련이 있다. 식 (15)는 미분에 대한 똑같은 내용을 계차방정식과 미분방정식으로 각각 다르게 표현한 것이다. 

 

사실 현실세계의 문제는 미분방정식이지만, 이를 인간이 보다 이해하기 쉽고 무엇보다도 컴퓨터가 계산하기 좋게 만든 것이 계차방정식이다. 계차방정식도 n과 n+1사이의 간격(or 주기 등)이 짧으면 짧을수록 미분방정식에 가깝게 된다. 

 

 

다시 선형대수 문제로 돌아와서, 우리가 어떤 순서나 시간의 개념이 들어가있는 계차방정식(difference equation)을 선형대수를 이용하여 푼다고 가정해보자. 시간의 개념이 있기때문에 분명 초기값(initial value or initial condition)이라는 것이 존재할 것이다. 이 방정식을 u에 대한 계차방정식이라고 했을 때, 식은 아래와 같다. 

 

 

식 (16)은 u에 대한 계차방정식이며 k번째 u를 시스템 A를 통해 계산했을 때 k+1번째 u의 값을 계산하는 방정식이다. 여기서 u가 순차적으로 변화한다고 했을 때 어느 타이밍 혹은 어느 순서에서 어떤 값이 나올지를 계산하기 위함이므로 u는 초기값 u0상태에서 시작한다. 이를 실제 문제에 빗대어 설명하자면 어느 도시의 인구 증가율을 A라는 방정식으로 만들고, u0가 도시의 초기 인구, k를 초기 년도로부터 몇년 이후의 년도인지를 나타내는 인덱스라고 하자. uk는 도시의 초기 인구 u0로 부터 k년도가 지났을 시점의 인구의 값이다. 이와 같은 문제를 계차방정식(difference equation)으로 표현하면 아래의 식과 같다. 

 

 

식 (17.1)을 보면 초기값 u0에서 A라는 방정식을 통해 u의 다음 시점의 변화된 값 u1을 계산할 수 있다. u1은 다시 A를 통해 그 다음 시점의 값 u2가 되고, k번째 u의 값 uk는 A를 통해 그 다음 차례의 u값인 uk+1이 된다. 그런데 u2를 계산할 때를 보면 u1을 A에 곱하여 계산하는데, u1은 Au0이다. 따라서 식 (17.1)의 u1을 Au0로 치환시켜주면 u2=AAu0가 되고, 결국 A가 2번 거듭제곱 된 형태이다. 이 식을 일반적으로 쓰면 u의 k번째 값 uk를 계산하기 위해선 u0를 A의 k번째 거듭제곱식에 곱해주면 된다. 결국 식 (17.2)의 상자속 식과 같은 일반적인 식이 도출된다. 이 식이 바로 계차방정식(difference equation)을 나타내며, 1차 시스템(first order system)이다. 위의 식이 1차(first order)인 이유는 식 (17.1)과 같이 k에서 k+1로 딱 한 단계만을 연결시키기 때문이고, 시스템(system)이라고 부를 수 있는 것은 우리가 구하고자 하는 미지수(unknown) 단일 숫자가 아닌 벡터(vector)이기 때문이다. 

 

 

- Real solution of difference equation

 

결국 우리는 식 (17.2)의 상자속 식을 통해 u의 100번째, 혹은 1000번째 값을 u의 초기값 u0에 A의 거듭제곱을 곱해서 쉽게 구할 수 있다. 사실 A행렬을 직접 100번 곱해서 u100을 구할 수도 있지만, 앞서 배운 행렬의 대각화(diagonolization)를 이용하면 이 계차방정식을 진짜로 풀 수 있다. 여기서 진짜로 풀 수 있다는 의미는 이 계차방정식에서 k가 커질수록 미지수의 값이 얼마만큼 증가하고 감소하는지 그 증가폭을 알 수 있다는 것이다. 즉 이 계차방정식의 Dynamics를 아는 것이 우리가 구하고자 하는 진정한 해답이고 이를 대각화를 통해 알 수 있다. 단순히 A만 거듭해서 곱해서는 이 시스템에 대한 자세한 정보를 알 수 없다. 이제 계차방정식에 대한 진정한 해를 어떻게 구하는지 알아보자. 

 

우선 주어진 방정식은 초기값 u0가 주어진다. 계차방정식을 진짜로 푸는 과정은 이 초기값인 u0 벡터를 행렬 A의 고유벡터들의 선형 결합(Linear combination)으로 표현하는 것으로부터 시작한다. 아래 식을 보자. 

 

 

식 (18)은 초기값 벡터 u0를 행렬 A의 고유벡터(eigenvectors)들의 선형 결합으로 표현한 것이다(x는 고유벡터, c는 상수). 이를 행렬의 곱 형태로 나타내면 고유벡터행렬(eigenvector matrix)인 S와 상수값 벡터인 c와의 곱인 u0=Sc로 나타낼 수 있다. 그런데 u0는 A의 고유벡터들의 선형 결합으로 표현이 가능한가? 정답은 Yes이다. 단 A가 n개의 독립(independent)인 고유벡터를 가지고 있다는 전제가 있어야 한다. n개의 독립인 고유벡터는 n차원 공간인 Rn의 기저(basis)를 형성하기 때문이다. 따라서 u0는 어떤 값이던 Rn공간내에 존재하기 때문에 A의 고유벡터들의 선형결합으로 표현이 가능한다. 

 

이제 식 (18)의 u0에 행렬 A를 곱하면 어떻게 정리가 되는지 살펴보자. 

 

 

식 (18)의 u0에 행렬 A를 곱하면 (19)와 같이 각각의 분리된 파트에 A가 곱해지는 것과 같다. c는 상수이므로 앞으로 빼서 정리하면 (19.1)처럼 정리할 수 있다. 여기서 각 파트는 cAx의 형태가 되는데, Ax는 고유값/고유벡터 식 (1)에 의해 (19.2)와 같이 λx로 바꿔서 정리할 수 있다. 결국 초기값 u0를 A를 통해 변환시킨 결과를 식 (19)와 같이 고유값과 고유벡터의 선형결합으로 나타냈다. 그러나 우리가 알고자하는 것은 A를 지속적으로 곱했을 때의 결과, 즉 예를 들면 A의 100제곱을 하면 초기값 u0가 어떻게 변화하는지, 그 변화의 정도는 얼마만큼인지 등이다. 식 (19)에서 A의 100제곱을 했을 때 식이 어떻게 되는지 살펴보자. 

 

 

우선 식 (19)에서 A를 한 번 더 곱하면 상수인 c와 λ를 앞으로 빼서 정리하여 (20.1)과 같이 만들 수 있다. 다시 고유값/고유벡터 식을 통해 치환하여 정리하면 식 (20.2)와 같이 되고, 이는 A를 제곱한 수 만큼 각 파트의 람다값이 제곱되는 꼴이 된다. A의 100제곱은 식 (20.3)과 같이 각 파트의 고유값인 람다(lambda)가 100제곱이 되고 이들의 조합으로 정리할 수 있다. 

 

다음으로 식 (20.3)을 앞서 배웠던 행렬의 대각화(diagonalization)방법을 통해 정리해보자. 행렬 A를 대각화하고 u0를 Sc로 나타내면 가운데 S행렬이 단위행렬이 되어 소거되고 식 (20.4)와 같이 정리가 된다. 이를 행렬(matrix)의 형태로 나타낸 것이 식 (20.5)이며, 이 식을 정리하면 (20.3)과 같은 결과가 나오는 것을 알 수 있다. 이를 통해 우리는 초기값 u0가 어떤 방정식 시스템 A에 의해 100번 변환된 결과인 u100을 고유값과 고유벡터로 풀 수 있다. 결과적으로 초기값 u0를 A의 k거듭제곱번 변화시킨 계차방정식(difference equation)은 고유벡터행렬(eigenvector matrix)과 고유값 행렬(eigenvalue matrix)의 k거듭제곱, 그리고 상수벡터 c의 곱으로 식 (20.6)과 같이 정리할 수 있다. 

 

정리해보면 어떤 계차방정식을 제대로 푼다는 것은 이를 고유값과 고유벡터로 풀어내는 것을 의미하며, 여기서 중요한 것은 고유값을 통해 우리는 이 계차방정식의 변화의 추이를 가늠해볼 수 있다는 것이다. 이는 단순히 A를 k번 거듭제곱하여 결과값만을 도출하는 것과는 달리 해당 시스템의 특성을 정확하게 파악하는데에 그 의의가 있다고 볼 수 있다. 다음 section에서 실제 예를 통해 이에 대한 이해를 높여보자. 

 

 

4. 대각화의 응용(Application of the Diagonalization)

 

- Fibonacci number

 

대각화(diagonalization)를 피보나치 수열(Fibonacci sequence)에 응용해보도록 하자. 피보나치 수열은 맨 처음 0과 1로 시작해서 다음 숫자는 앞의 숫자 두 개의 합이 되는 규칙을 가지고 있다. 

 

 

식 (21)은 피보나치 수열을 나타내며 F1, F2는 피보나치 수열에서 각 원소의 인덱스를 나타낸다. 여기서 우리가 알고자하는 것은 이런식으로 수열이 지속적으로 증가할 때 F100번째 값은 무엇인지, 그리고 얼마나 빠르게 값이 증가하는지 이다. 피보나치 수열은 보다시피 F1, F2, ... 와 같이 순서가 더해질 수록 값이 증가하는 체제이다. 이와 같이 시간이나 순서에 따라 특정 규칙으로 값이 변화하는 체제를 우리는 앞서 배웠던 계차방정식(Difference equation)으로 만들 수 있고, 대각화를 통해 고유값과 고유벡터들의 조합으로 표현할 수 있다. 피보나치 수열을 계차방정식으로 만들면 특정 인덱스에서의 값이 무엇인지 알 수 있고, 고유값을 통해 이 계차방정식의 증가/감소 폭이 얼마나 되는지 알 수 있다. 이제 피보나치 수열의 계차방정식을 만들어보자. 

 

피보나치 수열의 규칙은 다음과 같이 정리할 수 있다. 

 

 

식 (22)와 같이 Fk+2는 이전의 두 개의 값 Fk+1과 Fk의 합으로 나타낼 수 있다. 이제 이 수식을 식 (17.1)과 같이 선형계차방정식(Linear Difference Equation)으로 만들어보자. 그런데 두 가지 문제가 있다. 일단 현재로썬 식이 (22) 하나뿐이다. 어떤 시스템이 되기 위해선 식이 두 개 이상 존재해야 하고, 해(solution)도 벡터로써 존재해야 한다. 또 하나의 문제는 식 (22)가 2차미분방정식(second-order differential equation)과 같다는 것이다(여기서 식 자체는 계차방정식으로 표현함). 식을 보면 Fk가 두 단계 이후의 Fk+2와 연결되어 있고 영향을 미치게 된다. 따라서 k가 두 번에 걸쳐 변화할 때 그 값에 대한 변화가 어떻게 이루어지는지 보기 때문에 이차미분과 같다. 식 (17.1)과 같이 1차식의 시스템으로 만들기 위해선 약간의 트릭이 필요하다. 먼저 시스템의 입력값인 uk가 필요한데, 벡터형태여야 한다. 이 문제에서 uk벡터는 어떤식으로 만들어야 할까? 

 

피보나치 수열을 정의한 식 (22)를 보면 Fk+2라는 결과값은 Fk+1 + Fk의 이전 두 개의 입력값으로 정의되어 있다. 즉 입력값은 Fk+1과 Fk가 될 것이고 이들을 uk벡터로 정의하면 될 것이다. 또한 식 (17.1)과 같이 1차 시스템이 되기 위해선 어떤 선형시스템 A에 uk를 곱했을 때 uk+1이 나와야한다. 즉 uk+1은 uk의 원소들인 Fk+1과 Fk가 각각 한 단계씩 더 나아간 상태, 즉 Fk+2와 Fk+1이 될 것이다. 이를 식으로 정리하면 아래와 같다. 

 

 

uk와 uk+1를 만들었으니 다음 단계는 시스템 행렬 A를 만드는 것이다. 우리가 현재 하고 있는 작업은 단일 방정식으로 구성된 2차 미분방정식 (22)를 1차 선형계차방정식으로 만드는 것이다. 따라서 식 (22)의 내용이 그대로 반영이 된 1차 시스템식이 나와야 한다. 그렇다면 식 (23)의 uk를 곱했을 때 uk+1이 되면서 식 (22)의 내용이 반영되기 위해선 시스템 A가 어떻게 구성되어야 할까? 아래의 식을 보자. 

 

 

식 (24.1)이 피보나치 수열의 시스템 방정식(system equation)이다. 첫 번째 줄은 식 (22)와 동일하고 두 번째 줄은 1차식으로 만들기 위한 트릭이다. (24.1)의 좌변이 결국 uk+1이 되는 것을 볼 수 있다. 우변의 식들은 미지수(unknown) uk가 어떤 시스템 행렬 A에 곱해져서 만들어지는 행렬인데, 이와 같은 식이 되기 위해서는 식 (24.2)의 A와 같은 식이 만들어지면 된다. 

 

이렇게 하여 우리는 단일 식으로 구성된 2차 미분방정식의 꼴인 피보나치 수열 문제를 1차 시스템으로 바꾸었다. 이제 다음으로 할 일은 이 1차 선형방정식을 풀기 위해 식 (20.6)과 같이 고유값과 고유벡터로 표현하는 것이다. 먼저 행렬 A의 고유값과 고유벡터를 구해보자. 

 

 

우선 식 (25.1)의 행렬 A를 보면 대각 행렬(diagonal matrix)형태인 것을 알 수 있다. Lecture 21-(2)에서 배웠듯이 대각 행렬의 고유값은 허수(imaginary number)가 아닌 실수(real number)가 나온다. 또한 trace A가 고유값들의 합과 같은 규칙을 통해 고유값들의 합이 1이 된다는 것을 알 수 있고, det A가 고유값들의 곱과 같다는 규칙을 통해 고유값들의 곱이 -1이 됨을 알 수 있다. 실제로 식 (25.2)와 같이 고유값을 계산하여 더하고 곱해보면 trace A=1, det A=-1와 각각 같음을 알 수 있다. 

 

여기서 한 가지 주목할 점은 (25.1)의 고유값에 대한 2차식이다. 이 고유값의 방정식은 피보나치 수열의 계차방정식 행렬 A를 통해 만든 것이고, 행렬 A는 피보나치 수열의 규칙을 나타내는 식 (22)에 약간의 트릭을 써서 1차식으로 만든 것이다. 결국 최초의 2차 미분 방정식 (22)로부터 만들어온 것인데, 실제로 고유값에 대한 식을 구해보니 이 최초의 식과 같은 양상을 보인다. 

 

 

결국 행렬 A의 고유값에 대한 det(A-λI)는 A의 특성방정식(characteristic equation)을 만들어내고, 이 특성방정식은 해당 행렬이 수행하는 어떤 작업에 대한 특성을 나타내는 식이 된다. 최초의 피보나치 수열에 대한 식 (22)의 형태가 행렬 A의 특성방정식으로 고스란히 드러남을 통해 이러한 사실을 알 수 있다

 

이렇게 하여 식 (25.2)와 같이 2x2행렬 A에 대한 2개의 고유값을 구했다. 첫 번재 고유값 λ1은 1보다 크고, λ2는 1보다 작다. 이와 같이 2개의 고유값이 서로 다르기 때문에 A는 2개의 독립인 고유벡터가 존재하며, 대각화(diagonalization)가 가능하다

 

 

이제 처음에 가졌던 질문을 다시 생각해보자. 피보나치 수열은 계속 증가하는 함수이다. 그렇다면 얼마나 빠르게 증가하는가? 그리고 임의의 순서의 수열의 값, 예를 들어 F100의 값은 얼마가 되는가? 첫 번째 질문에 대한 답은 고유값에 있다. 얼마나 빠르게 증가하는지, 즉 피보나치 수열의 Dynamics는 고유값(eigenvalue)에 있다. 두 번째 질문에 대한 답은 앞서 정리했던 식 (17.2)를 활용하면 된다. 식 (17.2)를 다시 고유값과 고유벡터에 대한 식 (20.3)으로 분리하면 아래의 식과 같이 정리할 수 있다. 

 

 

F100의 값은 초기값 u0에 A의 100제곱을 곱한 u100을 구하는 것과 같다. 이 식은 (27.2)와 같이 고유값과 고유벡터들의 선형 조합으로 정리할 수 있고, 실제 고유값을 넣어서 (27.3)과 같이 정리할 수 있다. 이때 두 번째 고유값인 λ2는 1보다 작은 수이기 때문에 거듭제곱을 할 수록 0으로 수렴하여 실제로 F100에 미치는 영향은 0으로 봐도 무방하다. 따라서 (27.3)의 첫 번째 term만 고려해도 된다. 결국 위의 식에서 피보나치 수열의 값의 증가를 제어하는 요소는 고유값이고, 그 중에서도 1보다 큰 첫 번째 고유값이다. 이제 식을 완성시키기 위해 고유벡터(eigenvector)를 구해보자. 

 

 

기존에 고유벡터를 구하던 방법과 같이 각 고유값을 대입하고 가우스소거를 한 다음, free variable에 1을 설정하여 고유벡터를 구할 수도 있으나, 고유값을 뺀 행렬이 특이 행렬(singular matrix)이고, 곱한 결과가 영벡터가 되야 함을 염두해 봤을 때 고유벡터는 x=[λ 1]T임을 유추할 수 있다. 이때 식 (28.2)의 row1의 식은 식 (25.1)의 특성방정식이다. 고유벡터의 식에 고유값인 λ1, λ2를 각각 대입하면 (28.3)과 같이 고유벡터를 구할 수 있다. 

 

마지막으로 구할 파라미터는 계수값 c1과 c2이다. 이 둘은 아래의 식으로부터 구할 수 있다. 

 

 

먼저 우리가 이미 알고있는 초기값 u0를 구하는 계차방정식을 세운 다음, 이를 고유값과 고유벡터에 관한 식으로 정리한다. 초기값을 구하는 계차방정식은 행렬 A가 0제곱이기 때문에 식 (29.2)와 같이 고유값도 역시 0제곱이 되어 사라지고 u0=c1x1+c2x2만 남게 된다. 이를 (29.3)과 같이 고유벡터행렬(eigenvector matrix) S와 계수 벡터(coefficient vector) c의 곱으로 나타낼 수 있다. 이제 식 (29.3)을 가우스소거법(Gauss Elimination)으로 풀면 계수벡터 c를 구할 수 있다. S행렬의 소거 과정과 c의 값을 구하는 과정은 아래 식과 같다. 

 

 

식 (30.1)에서 S의 첫 번째 pivot원소가 있는 row1에 역수를 곱해 1을 만들어 주고, row2에서 row1빼주면 (30.2)와 같이 정리가 된다. 이때 식에 루트가 있기 때문에 계산 과정에서 분모의 유리화(Rationalization)를 해준다. 소거된 식으로부터 후방대입법(back-substitution)을 통해 c1과 c2를 각각 계산할 수 있다. 

 

필요한 파라미터들인 c, xλ들을 다 구했기 때문에 식 (27.3)에 이 값들을 대입하여 피보나치 수열의 100번째 값인 F100에 대한 값을 행렬 A를 100번 제곱하지 않아도 고유값의 100제곱을 통해 계산할 수 있다. 또한 고유값을 통해 해당 시스템(피보나치 수열)이 얼마나 빠르게 증가하는지를 파악할 수 있다.  

 

 

5. 마치며

 

이번 강좌에선 행렬의 대각화(diagonalization)에 대해 공부하였다. 대각화의 핵심엔 고유값(eigenvalue)과 고유벡터(eigenvector)가 있으며, 어떤 행렬을 고유값행렬(eigenvalue matrix)과 고유벡터행렬(eigenvalue matrix)로 분해하여 표현하는 방법이다. 대각화를 통해 어떤 반복적인 행렬곱셈(Matrix multiplication)을 해야 하는 문제를 훨씬 적은 계산량으로 효과적으로 풀 수 있으며, 이를 활용해 계차방정식(difference equation)의 형태로 정리한 피보나치 수열 문제를 해결하였다. 중요한것은 대각화를 통해 고유값/고유벡터로 분해하여 문제를 해결한다면 그 시스템 A에 대한 Dynamics를 파악할 수 있다는 것이다. 행렬 A의 거듭제곱(power)이 지속될 수록 얼마만큼 빠르게, 혹은 느리게 값이 증가 혹은 감소하는지에 대해서 분석할 수 있다. 

또한 행렬의 대각화는 다음 포스팅에서 다룰 미분방정식(differential equation)을 푸는 중요한 방법이 되므로 잘 공부하도록 하자. 

 

+ Recent posts