Kwangmin Kim - Ch.6 §6.3 — 미분방정식에의 응용 (Applications to Differential Equations)

1 이산에서 연속으로 — 같은 이야기, 다른 시간

앞 절에서는 이산 동역학 \(\mathbf{u}_{k+1} = \mathbf{A}\mathbf{u}_k\) 와 그 해 \(\mathbf{u}_k = \mathbf{A}^k \mathbf{u}_0\) 를 다루었다. 고유값을 찾고 대각화하면 \(\mathbf{A}^k = \mathbf{S}\mathbf{\Lambda}^k\mathbf{S}^{-1}\) 로 거듭제곱이 쉽게 풀린다는 이야기였다.

이제 시간을 연속으로 바꾼다:

\[\frac{d\mathbf{u}}{dt} = \mathbf{A}\mathbf{u}, \quad \mathbf{u}(0) = \mathbf{u}_0\]

이 식은 “매 순간 \(\mathbf{u}\) 의 변화율이 \(\mathbf{A}\mathbf{u}\) 이다”라는 뜻이다. 이산 버전에서 한 스텝에 \(\mathbf{A}\) 를 곱했다면, 연속 버전에서는 무한히 작은 스텝마다 \(\mathbf{A}\, dt\) 를 곱한다. 그 결과 이산의 \(\lambda^k\) 가 연속에서는 \(e^{\lambda t}\) 로 바뀐다.

핵심 직관 — 이산과 연속의 대응표

이산 (\(\mathbf{u}_{k+1} = \mathbf{A}\mathbf{u}_k\))	연속 (\(d\mathbf{u}/dt = \mathbf{A}\mathbf{u}\))
해: \(\mathbf{u}_k = \mathbf{A}^k \mathbf{u}_0\)	해: \(\mathbf{u}(t) = e^{\mathbf{A}t}\mathbf{u}_0\)
고유 성분: \(c_i \lambda_i^k \mathbf{x}_i\)	고유 성분: \(c_i e^{\lambda_i t} \mathbf{x}_i\)
안정: \(\lvert\lambda_i\rvert < 1\)	안정: \(\text{Re}(\lambda_i) < 0\)
지배 모드: 가장 큰 \(\lvert\lambda_i\rvert\)	지배 모드: 가장 큰 \(\text{Re}(\lambda_i)\)

두 세계는 “\(\lambda^k\) vs \(e^{\lambda t}\)” 라는 함수 하나 차이이다. 나머지 논리(고유벡터 분해, 지배 모드, 안정성)는 완전히 똑같이 작동한다.

이 대응은 우연이 아니다. \(\lambda^k\) 와 \(e^{\lambda t}\) 모두 “단위 스텝을 곱셈으로 반복하는 함수”이며, 미분 \(\frac{d}{dt}\) 와 \(e^{\lambda t}\) 의 관계(\(\frac{d}{dt}e^{\lambda t} = \lambda e^{\lambda t}\))가 \(\mathbf{A}\) 와 \(\mathbf{x}\) 의 관계(\(\mathbf{A}\mathbf{x} = \lambda\mathbf{x}\))와 구조적으로 같기 때문이다.

2 스칼라 준비 운동 — \(du/dt = \lambda u\)

벡터 버전으로 가기 전에 가장 단순한 1차원 문제를 확실히 해두자.

\[\frac{du}{dt} = \lambda u, \quad u(0) = u_0\]

해는 \(u(t) = u_0 e^{\lambda t}\) 이다. 이것은 미적분학의 가장 기본 결과이다. 검증은 간단하다: 양변을 미분하면 좌변 \(\frac{du}{dt} = \lambda u_0 e^{\lambda t}\) 와 우변 \(\lambda u = \lambda u_0 e^{\lambda t}\) 가 일치한다.

\(\lambda\) 의 의미

\(\lambda\)	해의 거동
\(\lambda > 0\)	지수적 성장 (\(u \to \infty\))
\(\lambda = 0\)	상수 유지 (\(u = u_0\))
\(\lambda < 0\)	지수적 감쇠 (\(u \to 0\))
\(\lambda = i\omega\) (순허수)	순수 진동 (\(u = u_0 e^{i\omega t}\), 크기 불변)
\(\lambda = r + i\omega\)	나선형 (크기는 \(e^{rt}\), 위상은 \(\omega t\))

이 스칼라 분류가 벡터 버전에도 그대로 적용된다. 벡터 문제는 이런 스칼라 문제를 “고유벡터 방향마다 독립적으로” 푸는 것일 뿐이다.

3 벡터 방정식의 해 — 왜 \(e^{\lambda t}\mathbf{x}\) 인가

\(n\) 차원 연립 미분방정식을 푸는 핵심 가설은:

\[\mathbf{u}(t) = e^{\lambda t}\mathbf{x}\]

라는 “순수 지수 해(pure exponential solution)” 를 찾는 것이다. 이것을 원래 방정식에 대입해 보자:

좌변: \(\dfrac{d\mathbf{u}}{dt} = \dfrac{d}{dt}(e^{\lambda t}\mathbf{x}) = \lambda e^{\lambda t}\mathbf{x}\)
우변: \(\mathbf{A}\mathbf{u} = \mathbf{A}(e^{\lambda t}\mathbf{x}) = e^{\lambda t}\mathbf{A}\mathbf{x}\)

양변이 같으려면 \(\mathbf{A}\mathbf{x} = \lambda\mathbf{x}\). 즉, \(\lambda\) 는 고유값이고 \(\mathbf{x}\) 는 고유벡터여야 한다.

\[\boxed{\,\mathbf{A}\mathbf{x} = \lambda\mathbf{x} \iff \mathbf{u}(t) = e^{\lambda t}\mathbf{x}\ \text{가}\ \frac{d\mathbf{u}}{dt} = \mathbf{A}\mathbf{u}\ \text{의 해}\,}\]

왜 이것이 기적인가? 원래 방정식 \(\frac{d\mathbf{u}}{dt} = \mathbf{A}\mathbf{u}\) 는 \(n\) 개의 얽힌 1계 방정식이다. 각 성분 \(u_i\) 의 변화율이 다른 성분들에 의존한다. 그런데 고유벡터 방향으로 움직이는 해는 자기 자신에만 의존하여 진화한다. 벡터 전체가 “스칼라 방정식 하나처럼” 행동하는 것이다.

일반 해 (General Solution)

\(\mathbf{A}\) 가 \(n\) 개의 독립 고유벡터 \(\mathbf{x}_1, \ldots, \mathbf{x}_n\) 을 가지면, \(d\mathbf{u}/dt = \mathbf{A}\mathbf{u}\) 의 일반 해는:

\[\mathbf{u}(t) = c_1 e^{\lambda_1 t}\mathbf{x}_1 + c_2 e^{\lambda_2 t}\mathbf{x}_2 + \cdots + c_n e^{\lambda_n t}\mathbf{x}_n\]

상수 \(c_1, \ldots, c_n\) 은 초기 조건 \(\mathbf{u}(0) = c_1\mathbf{x}_1 + \cdots + c_n\mathbf{x}_n\) 에서 결정된다.

구조적 해석: 초기 벡터 \(\mathbf{u}(0)\) 를 고유벡터 기저로 분해하면 각 성분이 자기만의 지수 속도 \(e^{\lambda_i t}\) 로 독립 진화한다. 이것은 이전 절의 \(\mathbf{A}^k\) 이야기와 완벽하게 평행하다. 차이는 단 하나, “\(\lambda^k\)” 가 “\(e^{\lambda t}\)” 로 바뀌었다는 것뿐이다.

4 풀이 절차 (3단계)

이 절차는 이산의 \(\mathbf{A}^k\mathbf{u}_0\) 를 풀 때와 구조적으로 동일하다.

Step 1 — 고유값·고유벡터 구하기: \(\det(\mathbf{A} - \lambda\mathbf{I}) = 0\) 로 \(\lambda_1, \ldots, \lambda_n\) 을 구하고 각각의 고유벡터 \(\mathbf{x}_i\) 를 찾는다.

Step 2 — 초기 조건 분해: 초기 벡터를 고유벡터의 선형결합으로 표현한다.

\[\mathbf{u}(0) = c_1\mathbf{x}_1 + c_2\mathbf{x}_2 + \cdots + c_n\mathbf{x}_n\]

행렬 형태로는 \(\mathbf{S}\mathbf{c} = \mathbf{u}(0)\) 를 풀어 \(\mathbf{c} = \mathbf{S}^{-1}\mathbf{u}(0)\).

Step 3 — 시간 진화: 각 고유 성분에 \(e^{\lambda_i t}\) 를 곱해 더한다.

\[\mathbf{u}(t) = \sum_{i=1}^n c_i e^{\lambda_i t}\mathbf{x}_i\]

4.1 예시 1 — \(2 \times 2\) 정상 상태 문제

\[\mathbf{A} = \begin{bmatrix} -1 & 2 \\ 1 & -2 \end{bmatrix}, \quad \mathbf{u}(0) = \begin{bmatrix} 1 \\ 0 \end{bmatrix}\]

Step 1: \(\det(\mathbf{A} - \lambda\mathbf{I}) = (-1-\lambda)(-2-\lambda) - 2 = \lambda^2 + 3\lambda = \lambda(\lambda + 3) = 0\) 에서 \(\lambda_1 = 0\), \(\lambda_2 = -3\).

고유벡터:

\(\lambda_1 = 0\): \(\mathbf{A}\mathbf{x} = \mathbf{0}\) 를 풀면 \(\mathbf{x}_1 = \begin{bmatrix}2\\1\end{bmatrix}\)
\(\lambda_2 = -3\): \((\mathbf{A} + 3\mathbf{I})\mathbf{x} = \begin{bmatrix}2 & 2\\1 & 1\end{bmatrix}\mathbf{x} = \mathbf{0}\) 를 풀면 \(\mathbf{x}_2 = \begin{bmatrix}1\\-1\end{bmatrix}\)

Step 2: \(c_1\begin{bmatrix}2\\1\end{bmatrix} + c_2\begin{bmatrix}1\\-1\end{bmatrix} = \begin{bmatrix}1\\0\end{bmatrix}\)

첫 식 \(2c_1 + c_2 = 1\), 둘째 식 \(c_1 - c_2 = 0\) → \(c_1 = c_2 = \frac{1}{3}\).

Step 3:

\[\mathbf{u}(t) = \frac{1}{3}e^{0}\begin{bmatrix}2\\1\end{bmatrix} + \frac{1}{3}e^{-3t}\begin{bmatrix}1\\-1\end{bmatrix} = \frac{1}{3}\begin{bmatrix}2\\1\end{bmatrix} + \frac{1}{3}e^{-3t}\begin{bmatrix}1\\-1\end{bmatrix}\]

거동 분석: \(t \to \infty\) 일 때 \(e^{-3t} \to 0\) 이므로 \(\mathbf{u}(t) \to \frac{1}{3}(2, 1)^\top\) 라는 정상 상태(steady state)로 수렴한다. 이것은 \(\lambda_1 = 0\) 에 대응하는 고유벡터 방향이다. 고유값 0 = “변하지 않는 방향” = 정상 상태라는 원리가 드러난다.

4.2 예시 2 — 상삼각 \(3 \times 3\)

\[\frac{d\mathbf{u}}{dt} = \begin{bmatrix} 1 & 1 & 1 \\ 0 & 2 & 1 \\ 0 & 0 & 3 \end{bmatrix}\mathbf{u}, \quad \mathbf{u}(0) = \begin{bmatrix}9\\7\\4\end{bmatrix}\]

삼각 행렬이므로 고유값은 대각 원소 그대로: \(\lambda_1 = 1, \lambda_2 = 2, \lambda_3 = 3\).

고유벡터는 \(\mathbf{x}_1 = (1, 0, 0)^\top\), \(\mathbf{x}_2 = (1, 1, 0)^\top\), \(\mathbf{x}_3 = (1, 1, 1)^\top\).

\((9, 7, 4) = 2\mathbf{x}_1 + 3\mathbf{x}_2 + 4\mathbf{x}_3\) 이므로:

\[\mathbf{u}(t) = 2e^t\mathbf{x}_1 + 3e^{2t}\mathbf{x}_2 + 4e^{3t}\mathbf{x}_3\]

모든 고유값이 양수이므로 해는 발산한다. 장기적으로는 가장 큰 고유값 \(\lambda_3 = 3\) 에 대응하는 \(\mathbf{x}_3\) 방향이 지배한다:

\[\mathbf{u}(t) \sim 4e^{3t}\mathbf{x}_3 \quad (t \to \infty)\]

5 안정성 분석

\(t \to \infty\) 에서 \(\mathbf{u}(t) \to \mathbf{0}\) 인지 여부를 안정성(stability)이라 한다. 해가 \(\sum c_i e^{\lambda_i t}\mathbf{x}_i\) 이므로, 안정성은 지수 \(e^{\lambda_i t}\) 의 거동에만 의존한다.

5.1 실수 고유값

\(\lambda \in \mathbb{R}\) 일 때 \(e^{\lambda t} \to 0\) 이려면 \(\lambda < 0\) 이어야 한다.

5.2 복소 고유값

실수 행렬의 복소 고유값은 반드시 켤레쌍 \(\lambda = r \pm i\omega\) 로 나타난다. Euler 공식으로 분해하면:

\[e^{\lambda t} = e^{(r + i\omega)t} = e^{rt}(\cos\omega t + i\sin\omega t)\]

크기는 \(\lvert e^{\lambda t}\rvert = e^{rt}\) 이므로:

실수부 \(r\) 이 성장/감쇠를 결정한다 (\(r < 0\) 이면 안정).
허수부 \(\omega\) 는 진동의 각주파수이다 (안정성과 무관).

직관: 복소 고유값은 “회전 + 스케일” 조합이다. 실수부가 음이면 나선형으로 원점에 수렴, 양이면 바깥으로 나선, 0이면 원 궤도로 순수 진동한다.

연속 시스템 안정성 정리

\(d\mathbf{u}/dt = \mathbf{A}\mathbf{u}\) 는 \(\mathbf{A}\) 의 모든 고유값 실수부가 음수일 때 점근적으로 안정하다:

\[\text{Re}(\lambda_i) < 0 \quad \forall i \iff \mathbf{u}(t) \to \mathbf{0}\ \text{as}\ t \to \infty\]

적어도 하나의 고유값이 \(\text{Re}(\lambda) = 0\) 이고 나머지는 \(< 0\) 이면 정상 상태(수렴은 하되 0이 아님)에 도달한다.

하나라도 \(\text{Re}(\lambda) > 0\) 이면 불안정(해가 발산).

5.3 \(2 \times 2\) 행렬의 안정성 지도

\(2 \times 2\) 행렬의 특성 다항식은 \(\lambda^2 - T\lambda + D = 0\) 이다. 여기서 \(T = \text{tr}(\mathbf{A}) = \lambda_1 + \lambda_2\) 이고 \(D = \det(\mathbf{A}) = \lambda_1\lambda_2\) 이다.

안정성은 \((T, D)\) 평면 위의 영역으로 깔끔하게 분류된다.

\(2 \times 2\) 안정성 판별

\(d\mathbf{u}/dt = \mathbf{A}\mathbf{u}\) 가 안정할 필요충분조건:

\[\text{tr}(\mathbf{A}) < 0 \quad \text{and} \quad \det(\mathbf{A}) > 0\]

이유를 풀어쓰면:

\(\lambda_1, \lambda_2\) 가 모두 실수 음수라면 \(\lambda_1 + \lambda_2 = T < 0\), \(\lambda_1\lambda_2 = D > 0\). 역으로 \(D > 0\) 이면 두 고유값이 같은 부호이고, \(T < 0\) 이면 그 부호가 음이다.
복소 켤레쌍 \(\lambda = r \pm i\omega\) 이면 \(D = r^2 + \omega^2 > 0\) 는 항상 양수이고 \(T = 2r\). \(T < 0 \iff r < 0 \iff\) 안정.

\((T, D)\) 평면의 분류:

영역	고유값 상태	동역학
\(T < 0, D > 0, T^2 > 4D\)	음의 실수 두 개	안정 노드 (직선 수렴)
\(T < 0, D > 0, T^2 < 4D\)	음의 실수부 복소쌍	안정 나선 (수렴 나선)
\(T < 0, D > 0, T^2 = 4D\)	중복 음의 실수	임계 수렴
\(T = 0, D > 0\)	순허수쌍	중립 (원형 궤도)
\(T > 0, D > 0\)	양의 실수부	불안정 (발산)
\(D < 0\)	부호 다른 실수	안장점(saddle) — 방향에 따라 발산

포물선 \(T^2 = 4D\) 가 실수와 복소 영역을 가른다. 아래쪽(판별식 양)은 실수 고유값, 위쪽은 복소쌍.

실용 관점: 2×2 시스템이라면 고유값을 직접 풀 필요 없이 trace와 determinant만 보고 안정성을 즉시 판정할 수 있다. 제어공학·생태학·전기회로 해석에서 광범위하게 쓰이는 빠른 판단 도구이다.

6 2계 미분방정식을 1계 시스템으로

역학에서 가장 중요한 방정식은 감쇠 조화 진동자이다:

\[my'' + by' + ky = 0\]

\(m\): 질량 — 관성을 나타낸다
\(b\): 감쇠 계수 — 속도에 비례하는 저항
\(k\): 강성 — 변위에 비례하는 복원력

이 2계 방정식은 변수 하나를 추가하는 고전적 기법으로 1계 연립 방정식이 된다.

6.1 변환 기법

상태 벡터를 \(\mathbf{u} = \begin{bmatrix}y\\y'\end{bmatrix}\) 로 두면:

첫 성분의 미분: \(\dfrac{d}{dt}y = y'\) (정의상 자명)
둘째 성분의 미분: \(\dfrac{d}{dt}y' = y'' = -\dfrac{k}{m}y - \dfrac{b}{m}y'\) (원 방정식에서 \(y''\) 를 풀어냄)

두 식을 한 번에 쓰면:

\[\frac{d}{dt}\begin{bmatrix}y\\y'\end{bmatrix} = \begin{bmatrix}0 & 1\\ -k/m & -b/m\end{bmatrix}\begin{bmatrix}y\\y'\end{bmatrix}\]

이제 \(d\mathbf{u}/dt = \mathbf{A}\mathbf{u}\) 형태이므로 위의 3단계 절차로 풀 수 있다. 이 변환은 \(n\) 계 스칼라 방정식을 항상 \(n\) 차원 1계 시스템으로 바꾼다 — 상태 변수로 \(y, y', y'', \ldots, y^{(n-1)}\) 를 쌓으면 된다.

왜 이 기법이 우아한가? 2계 방정식의 특성 방정식 \(m\lambda^2 + b\lambda + k = 0\) 을 그냥 푸는 것과 1계 시스템의 특성 방정식 \(\det(\mathbf{A} - \lambda\mathbf{I}) = 0\) 이 완전히 같은 방정식이 된다. 선형대수가 미분방정식의 계산을 흡수하는 순간이다.

6.2 감쇠의 3가지 경우

\(m\lambda^2 + b\lambda + k = 0\) 의 판별식 \(b^2 - 4mk\) 의 부호에 따라:

조건	고유값	물리적 의미
\(b^2 > 4mk\)	서로 다른 음의 실수	과감쇠 (진동 없이 수렴)
\(b^2 = 4mk\)	중복 음의 실수	임계감쇠 (가장 빠른 수렴)
\(b^2 < 4mk\)	음의 실수부 복소쌍	저감쇠 (진동하며 수렴)
\(b = 0\)	순허수쌍 \(\pm i\sqrt{k/m}\)	무감쇠 (순수 진동)

자동차 서스펜션 설계는 “임계감쇠에 가깝지만 약간 저감쇠”를 목표로 한다 — 빠르게 수렴하면서 약간의 부드러움을 유지하기 위함이다.

7 행렬 지수 함수 \(e^{\mathbf{A}t}\)

7.1 정의 — 테일러 급수

스칼라 지수 함수 \(e^x = 1 + x + \frac{x^2}{2!} + \frac{x^3}{3!} + \cdots\) 에서 \(x\) 자리에 행렬 \(\mathbf{A}t\) 를 넣는다:

행렬 지수 함수 (Matrix Exponential)

\[e^{\mathbf{A}t} = \mathbf{I} + \mathbf{A}t + \frac{(\mathbf{A}t)^2}{2!} + \frac{(\mathbf{A}t)^3}{3!} + \cdots = \sum_{k=0}^\infty \frac{(\mathbf{A}t)^k}{k!}\]

이 급수는 모든 정방행렬 \(\mathbf{A}\) 와 실수 \(t\) 에 대해 수렴한다.

미분하면 \(\frac{d}{dt}e^{\mathbf{A}t} = \mathbf{A}e^{\mathbf{A}t} = e^{\mathbf{A}t}\mathbf{A}\) 가 성립한다 (급수의 항별 미분으로 확인). 따라서:

\[\mathbf{u}(t) = e^{\mathbf{A}t}\mathbf{u}(0)\]

는 항상 \(d\mathbf{u}/dt = \mathbf{A}\mathbf{u}\) 의 해가 된다. 이 공식의 위력은 결함 행렬(독립 고유벡터가 부족한 경우)에도 작동한다는 점이다.

7.2 대각화 경로 — 실제 계산

\(\mathbf{A}\) 가 대각화 가능하면 \(\mathbf{A} = \mathbf{S}\mathbf{\Lambda}\mathbf{S}^{-1}\). 거듭제곱 공식 \(\mathbf{A}^k = \mathbf{S}\mathbf{\Lambda}^k\mathbf{S}^{-1}\) 를 테일러 급수에 대입하면:

\[e^{\mathbf{A}t} = \sum_{k=0}^\infty \frac{t^k}{k!}\mathbf{S}\mathbf{\Lambda}^k\mathbf{S}^{-1} = \mathbf{S}\left(\sum_{k=0}^\infty \frac{(\mathbf{\Lambda}t)^k}{k!}\right)\mathbf{S}^{-1} = \mathbf{S}\, e^{\mathbf{\Lambda}t}\, \mathbf{S}^{-1}\]

대각 행렬의 지수는 각 대각 원소의 지수이므로:

\[e^{\mathbf{\Lambda}t} = \begin{bmatrix} e^{\lambda_1 t} & & \\ & \ddots & \\ & & e^{\lambda_n t} \end{bmatrix}\]

최종 공식

\[\boxed{\,\mathbf{u}(t) = e^{\mathbf{A}t}\mathbf{u}(0) = \mathbf{S}\, e^{\mathbf{\Lambda}t}\, \mathbf{S}^{-1}\mathbf{u}(0)\,}\]

이 공식을 풀어쓰면 바로 3단계 절차가 된다:

\(\mathbf{S}^{-1}\mathbf{u}(0) = \mathbf{c}\): 초기 조건을 고유벡터 좌표로 변환
\(e^{\mathbf{\Lambda}t}\mathbf{c} = (c_1 e^{\lambda_1 t}, \ldots, c_n e^{\lambda_n t})^\top\): 각 성분 독립 진화
\(\mathbf{S}(\cdots)\): 표준 좌표로 복귀

7.3 주의 — 행렬 곱의 비가환성

스칼라에서는 \(e^{a+b} = e^a e^b\) 가 항상 성립하지만, 행렬은 일반적으로 \(e^{\mathbf{A}+\mathbf{B}} \neq e^\mathbf{A} e^\mathbf{B}\) 이다. 이 항등식은 \(\mathbf{A}\mathbf{B} = \mathbf{B}\mathbf{A}\) (교환 가능)일 때에만 성립한다. 이로 인해 시간에 따라 변하는 \(\mathbf{A}(t)\) 의 미분방정식은 본 절의 방법으로 풀리지 않고 시간 순서 지수(time-ordered exponential)라는 더 복잡한 도구가 필요하다.

7.4 결함 행렬의 경우

\(\mathbf{A}\) 가 독립 고유벡터 부족으로 대각화 불가능해도 \(e^{\mathbf{A}t}\) 는 여전히 급수로 정의된다. Jordan 표준형 \(\mathbf{A} = \mathbf{P}\mathbf{J}\mathbf{P}^{-1}\) 을 사용하면:

\[e^{\mathbf{A}t} = \mathbf{P}\, e^{\mathbf{J}t}\, \mathbf{P}^{-1}\]

이때 Jordan 블록의 지수에는 \(t e^{\lambda t}\), \(t^2 e^{\lambda t}\) 같은 다항식 × 지수 항이 나타난다. 예를 들어 \(\mathbf{J} = \begin{bmatrix}\lambda & 1\\0 & \lambda\end{bmatrix}\) 의 경우:

\[e^{\mathbf{J}t} = e^{\lambda t}\begin{bmatrix}1 & t\\0 & 1\end{bmatrix}\]

물리적으로는 임계감쇠에서 \(te^{-\alpha t}\) 형태로 자주 나타난다.

8 응용 분야

분야	모델	고유값의 역할
제어공학	선형 시스템 \(\dot{\mathbf{x}} = \mathbf{A}\mathbf{x} + \mathbf{B}\mathbf{u}\)	극점 = 고유값, 안정성·응답 특성 결정
전기회로	RLC 회로 ODE	저감쇠/과감쇠 분류, 공진 주파수
화학 반응	1차 반응 네트워크	반응 속도 상수 = 고유값
약동학	다중 구획 모델	흡수·분포·소실 반감기
인구 모델	Lotka-Volterra 선형화	평형점 주변 안정성
역학	질량-스프링 시스템	정규 모드 주파수
경제학	선형 거시 모델	장기 균형 수렴
수치적분	ODE 해법기 안정성	수치 스킴의 안정 영역

8.1 약동학 2-구획 모델

약물이 혈장(central compartment)과 조직(peripheral compartment) 사이를 오가며 배설되는 시스템을 생각해 보자.

\[\frac{d}{dt}\begin{bmatrix}C_1\\C_2\end{bmatrix} = \begin{bmatrix}-(k_{10}+k_{12}) & k_{21}\\ k_{12} & -k_{21}\end{bmatrix}\begin{bmatrix}C_1\\C_2\end{bmatrix}\]

\(k_{12}\): 혈장 → 조직 전달 속도
\(k_{21}\): 조직 → 혈장 전달 속도
\(k_{10}\): 혈장에서 배설 속도

\(\mathbf{A}\) 의 두 고유값은 모두 음수이고 해는 \(C_1(t) = A e^{\lambda_1 t} + B e^{\lambda_2 t}\) 형태의 두 지수의 합으로 나타난다. 작은 \(\lvert\lambda\rvert\) 는 “말단 반감기”, 큰 \(\lvert\lambda\rvert\) 는 “분포 반감기”에 해당한다. 이것이 임상 약동학에서 약물 혈중농도 곡선을 bi-exponential fit 으로 분석하는 수학적 근거이다.

9 코드 예시

import numpy as np

def solve_ode_eigen(A, u0, t):
    """
    du/dt = A u, u(0) = u0 를 고유값·고유벡터로 풀이
    반환: u(t) 벡터
    """
    # Step 1: 고유값·고유벡터
    eigenvalues, eigenvectors = np.linalg.eig(A)
    S = eigenvectors

    # Step 2: 초기 조건 분해 (c = S^-1 u0)
    c = np.linalg.solve(S, u0)

    # Step 3: 시간 진화
    exp_lambda_t = np.exp(eigenvalues * t)
    u_t = S @ (c * exp_lambda_t)
    return u_t.real  # 실수 행렬이면 허수부 수치 잔차 제거

# 예시 1: 정상 상태로 수렴
A = np.array([[-1.0, 2.0],
              [ 1.0, -2.0]])
u0 = np.array([1.0, 0.0])

print("예시 1: 정상 상태 수렴")
print(f"{'t':>6} {'u1(t)':>10} {'u2(t)':>10}")
for t in [0, 0.5, 1, 2, 5, 10]:
    u = solve_ode_eigen(A, u0, t)
    print(f"{t:>6.1f} {u[0]:>10.6f} {u[1]:>10.6f}")

# 정상 상태 이론값 = (2/3, 1/3)
print(f"\n이론적 정상 상태: (0.666667, 0.333333)")

import numpy as np
from scipy.linalg import expm

# scipy.linalg.expm 이 정확한 e^(At) 계산 제공
A = np.array([[-1.0, 2.0],
              [ 1.0, -2.0]])
u0 = np.array([1.0, 0.0])

def solve_via_matrix_exp(A, u0, t):
    return expm(A * t) @ u0

print("scipy.linalg.expm 기반 해")
for t in [0, 0.5, 1, 2, 5, 10]:
    u = solve_via_matrix_exp(A, u0, t)
    print(f"t={t:>4.1f}: u = {u}")

# 고유값 직접 확인
w, V = np.linalg.eig(A)
print(f"\n고유값: {w}")
print(f"안정성: 모두 Re(λ) ≤ 0 → {np.all(w.real <= 0)}")

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import expm

# my'' + by' + ky = 0 을 상태 공간으로 변환
m, k = 1.0, 10.0
cases = {'과감쇠': 8.0, '임계감쇠': 2*np.sqrt(m*k), '저감쇠': 1.0, '무감쇠': 0.0}

fig, ax = plt.subplots(figsize=(9, 5))
times = np.linspace(0, 10, 500)
u0 = np.array([1.0, 0.0])  # y(0)=1, y'(0)=0

for name, b in cases.items():
    A = np.array([[0.0, 1.0],
                  [-k/m, -b/m]])
    eigenvalues = np.linalg.eigvals(A)
    y_series = [(expm(A*t) @ u0)[0] for t in times]
    ax.plot(times, y_series, label=f'{name} (b={b:.2f}, λ={eigenvalues})')

ax.axhline(0, color='k', lw=0.5)
ax.set_xlabel('시간 t')
ax.set_ylabel('변위 y(t)')
ax.set_title('감쇠 조화 진동자: 감쇠 계수에 따른 거동')
ax.legend(fontsize=8)
ax.grid(alpha=0.3)
plt.tight_layout()

import numpy as np
import matplotlib.pyplot as plt

# (trace, det) 평면에서 무작위 2x2 행렬의 안정성 분류
rng = np.random.default_rng(42)
n = 500
As = rng.uniform(-3, 3, size=(n, 2, 2))

traces = np.array([np.trace(A) for A in As])
dets = np.array([np.linalg.det(A) for A in As])
eigenvalues = np.array([np.linalg.eigvals(A) for A in As])
stable = np.all(eigenvalues.real < 0, axis=1)

fig, ax = plt.subplots(figsize=(8, 6))
ax.scatter(traces[stable], dets[stable], c='tab:blue', s=15, label='안정', alpha=0.7)
ax.scatter(traces[~stable], dets[~stable], c='tab:red', s=15, label='불안정', alpha=0.7)

# 판별 곡선
T = np.linspace(-3, 3, 400)
ax.plot(T, T**2/4, 'k--', alpha=0.5, label='$T^2=4D$')
ax.axhline(0, color='k', lw=0.5)
ax.axvline(0, color='k', lw=0.5)
ax.set_xlabel('Trace $T$')
ax.set_ylabel('Determinant $D$')
ax.set_title('2x2 안정성 지도: $T < 0$ & $D > 0$ 영역만 안정')
ax.legend()
ax.grid(alpha=0.3)
plt.tight_layout()

library(expm)  # expm::expm 이 행렬 지수

A <- matrix(c(-1, 2, 1, -2), nrow = 2, byrow = TRUE)
u0 <- c(1, 0)

# 고유값 분해
ev <- eigen(A)
cat("고유값:", ev$values, "\n")

# 해의 시간 경로
times <- c(0, 0.5, 1, 2, 5, 10)
solve_ode <- function(t) as.vector(expm(A * t) %*% u0)
for (t in times) {
  u <- solve_ode(t)
  cat(sprintf("t=%.1f: u = (%.6f, %.6f)\n", t, u[1], u[2]))
}

# 안정성 판정 (2x2)
stable_2x2 <- function(A) {
  T_val <- sum(diag(A))
  D_val <- det(A)
  T_val < 0 && D_val > 0
}
cat("\n안정성:", stable_2x2(A), "\n")
cat("Trace:", sum(diag(A)), "Determinant:", det(A), "\n")

10 정리

개념	핵심 공식	직관
순수 지수 해	\(\mathbf{u}(t) = e^{\lambda t}\mathbf{x}\)	고유벡터 방향은 스칼라처럼 진화
일반 해	\(\mathbf{u}(t) = \sum c_i e^{\lambda_i t}\mathbf{x}_i\)	초기 조건을 고유벡터로 분해
행렬 지수	\(e^{\mathbf{A}t} = \mathbf{S}\, e^{\mathbf{\Lambda}t}\, \mathbf{S}^{-1}\)	좌표 변환 → 독립 진화 → 복귀
안정성	\(\text{Re}(\lambda_i) < 0\ \forall i\)	모든 모드가 감쇠
2×2 판별	\(T < 0\) & \(D > 0\)	Trace·Determinant 즉시 확인
복소 고유값	\(\lambda = r \pm i\omega\)	실수부=성장률, 허수부=진동 주파수
2계 변환	\(\mathbf{u} = (y, y')^\top\)	상태 변수 추가로 1계화

이 절의 본질은 한 줄이다: 미분을 곱셈으로 바꾼다. 미분 연산자 \(\frac{d}{dt}\) 와 행렬 \(\mathbf{A}\) 는 일반적으로 다루기 어렵지만, \(e^{\lambda t}\) 와 고유벡터 축을 도입하면 두 연산이 같은 스칼라 \(\lambda\) 의 곱으로 변한다. 이것이 선형대수가 미분방정식의 모든 연산을 흡수할 수 있는 이유이며, 공학과 과학의 수많은 동역학이 결국 \(\mathbf{A}\mathbf{x} = \lambda\mathbf{x}\) 한 식으로 수렴하는 이유이다 (Strang, 2009, Ch.6.3).

11 관련 주제

선행 지식

Ch.6 Overview — 고유값과 고유벡터 — 정의, 특성 방정식, 기하 직관
Ch.6 §6.2 — 행렬의 대각화 — \(\mathbf{A} = \mathbf{S}\mathbf{\Lambda}\mathbf{S}^{-1}\)

Ch.6 세부 주제

후속 주제