Lecture 23. 미분방정식과 exp(At) — 파이썬 실습 - 선형대수

L23 서술 파트의 출발점은 $\vv{u} = e^{\lambda t}\vv{x}$ 를 대입하니 $A\vv{x} = \lambda\vv{x}$ 가 튀어나오더라는 것이었다.

이 노트북에서는 그 대입을 코드로 직접 해서 잔차가 0이 되는 것을 보고, 앵커 예제를 손으로 푼 답과 수치적분의 답을 맞춰 보며, $e^{At}$ 를 급수로 직접 더해 scipy.linalg.expm 과 대조한다. 마지막으로 미적분에서 외운 특성방정식이 동반행렬의 특성방정식임을 확인한다.

서술 파트의 내용	여기서 확인하는 방법
대입하면 $A\vv{x} = \lambda\vv{x}$	잔차를 찍어 본다. 고유쌍일 때만 0
일반해와 $\vv{c} = S^{-1}\vv{u}(0)$	손으로 푼 식과 수치적분을 대조
열의 합이 0이면 총량 보존	$\vv{1}^{\mathsf T}A$ 와 $u_1 + u_2$
궤적이 직선	위상평면에 그려 본다
$\lvert e^{\lambda t}\rvert = e^{\operatorname{Re}(\lambda)t}$	허수부를 키워도 크기가 그대로
실수부가 운명	슬라이더로 나선에서 중심으로
$e^{At}$ 급수	항을 늘려 가며 `expm` 에 수렴
성분별이 아니다	두 결과를 나란히
$e^Ae^B \neq e^{A+B}$	교환자가 0인지에 달렸다
동반행렬	두 특성방정식이 같은 계수

0. 준비¶

import numpy as np
import plotly.graph_objects as go
from scipy.integrate import solve_ivp
from scipy.linalg import expm

from linalg_viz import COLORS, layout2d, show_matrix, slider_figure

np.set_printoptions(precision=4, suppress=True)
print("numpy", np.__version__)

numpy 2.5.2

1. 대입해 보면 고윳값이 나온다¶

$\vv{u} = e^{\lambda t}\vv{x}$ 를 $\frac{d\vv{u}}{dt} = A\vv{u}$ 에 넣었을 때 남는 것을 잔차라 하자.

\vv{r}(t) = \underbrace{\lambda e^{\lambda t}\vv{x}}_{d\vv{u}/dt} - \underbrace{e^{\lambda t}A\vv{x}}_{A\vv{u}}

(1)

이것이 모든 $t$ 에서 0이어야 해가 된다.

A = np.array([[-1.0, 2.0], [1.0, -2.0]])


def 잔차(A, lam, x, 시각=(0.0, 0.3, 1.0, 2.5)):
    """u = e^(lam t) x 를 대입했을 때 남는 것의 최대 크기."""
    x = np.asarray(x, dtype=float)
    남은것 = [lam * np.exp(lam * t) * x - np.exp(lam * t) * (A @ x) for t in 시각]
    return float(np.abs(남은것).max())

후보 = [
    ("고유쌍  lambda=0,  x=(2,1)", 0.0, [2, 1]),
    ("고유쌍  lambda=-3, x=(1,-1)", -3.0, [1, -1]),
    ("엉뚱한 벡터  lambda=-3, x=(1,0)", -3.0, [1, 0]),
    ("엉뚱한 값    lambda=-1, x=(1,-1)", -1.0, [1, -1]),
]
for 이름, lam, x in 후보:
    print(f"{이름:>34} :  잔차 {잔차(A, lam, x):.4f}")

           고유쌍  lambda=0,  x=(2,1) :  잔차 0.0000
          고유쌍  lambda=-3, x=(1,-1) :  잔차 0.0000
        엉뚱한 벡터  lambda=-3, x=(1,0) :  잔차 2.0000
      엉뚱한 값    lambda=-1, x=(1,-1) :  잔차 2.0000

고유쌍일 때만 정확히 0이다. 대입이 통하는 유일한 경우가 $A\vv{x} = \lambda\vv{x}$ 이다.

2. 앵커 예제를 끝까지 푼다¶

$A = \begin{bmatrix} -1 & 2 \\ 1 & -2\end{bmatrix}$ , $\vv{u}(0) = (1, 0)$ 이다.

print("대각합 :", np.trace(A), "   행렬식 :", np.linalg.det(A))
print("특성다항식 계수 :", np.poly(A), "  ->  lambda^2 + 3 lambda")
print("고윳값 :", np.sort(np.linalg.eigvals(A).real))

S = np.array([[2.0, 1.0], [1.0, -1.0]])          # 열이 고유벡터
L = np.diag([0.0, -3.0])
print()
print(show_matrix(S, "S"))
print("A S = S Lambda 인가 :", np.allclose(A @ S, S @ L))

대각합 : -3.0    행렬식 : 0.0
특성다항식 계수 : [1. 3. 0.]   ->  lambda^2 + 3 lambda
고윳값 : [-3.  0.]

S
[   2    1 ]
[   1   -1 ]
A S = S Lambda 인가 : True

u0 = np.array([1.0, 0.0])
c = np.linalg.solve(S, u0)
print("c = S^-1 u(0) =", c, "   (1/3, 1/3 인가 :", np.allclose(c, [1/3, 1/3]), ")")

c = S^-1 u(0) = [0.3333 0.3333]    (1/3, 1/3 인가 : True )

def 손으로푼해(t):
    """u(t) = c1 x1 + c2 e^(-3t) x2 를 그대로 옮긴 것."""
    t = np.atleast_1d(np.asarray(t, dtype=float))
    return (np.outer(np.ones_like(t), c[0] * S[:, 0])
            + np.outer(c[1] * np.exp(-3 * t), S[:, 1]))

for t in (0.0, 0.5, 1.0, 20.0):
    print(f"t = {t:>5} :  u = {손으로푼해(t)[0]}")
print()
print("t -> 무한대 의 정상상태 :", S[:, 0] / S[:, 0].sum(), "  = (2/3, 1/3)")

t =   0.0 :  u = [1. 0.]
t =   0.5 :  u = [0.741 0.259]
t =   1.0 :  u = [0.6833 0.3167]
t =  20.0 :  u = [0.6667 0.3333]

t -> 무한대 의 정상상태 : [0.6667 0.3333]   = (2/3, 1/3)

수치적분과 대조하자. solve_ivp 는 고윳값을 전혀 모른 채 $\frac{d\vv{u}}{dt} = A\vv{u}$ 를 잘게 쪼개 따라갈 뿐이다. 두 답이 같아야 한다.

시각 = np.linspace(0, 3, 301)
수치 = solve_ivp(lambda t, u: A @ u, (0, 3), u0, t_eval=시각, rtol=1e-10, atol=1e-12)
공식 = 손으로푼해(시각)

print("수치적분 성공 :", 수치.success)
print("두 답의 최대 차이 :", np.abs(수치.y.T - 공식).max())

수치적분 성공 : True
두 답의 최대 차이 : 8.526734873726127e-12

3. 총량이 보존되는 이유¶

$A$ 의 열의 합이 0이면 $\vv{1}^{\mathsf T}A = \vv{0}^{\mathsf T}$ 이고, 그러면 $\frac{d}{dt}(\vv{1}^{\mathsf T}\vv{u}) = \vv{1}^{\mathsf T}A\vv{u} = 0$ 이다.

하나 = np.ones(2)
print("A 의 열의 합 :", A.sum(axis=0))
print("1^T A =", 하나 @ A, "  -> 영벡터인가 :", np.allclose(하나 @ A, 0))
print()
합 = 공식.sum(axis=1)
print(f"u1 + u2 의 최솟값 {합.min():.12f}, 최댓값 {합.max():.12f}")
print("열의 합이 0 이라는 조건이 고윳값 0 을 만든다 :",
      np.isclose(np.linalg.det(A), 0))

A 의 열의 합 : [0. 0.]
1^T A = [0. 0.]   -> 영벡터인가 : True

u1 + u2 의 최솟값 1.000000000000, 최댓값 1.000000000000
열의 합이 0 이라는 조건이 고윳값 0 을 만든다 : True

왜 열의 합이 0이면 고윳값 0이 생기는가. $\vv{1}^{\mathsf T}A = \vv{0}^{\mathsf T}$ 를 전치하면 $A^{\mathsf T}\vv{1} = \vv{0}$ 이므로 $A^{\mathsf T}$ 가 특이하고, $\det A = \det A^{\mathsf T} = 0$ 이기 때문이다. L18의 성질이 여기서 쓰인다.

print("A^T 1 =", A.T @ 하나)
print("det A =", np.linalg.det(A), "   det A^T =", np.linalg.det(A.T))

A^T 1 = [0. 0.]
det A = 0.0    det A^T = 0.0

4. 궤적이 직선인 이유¶

$\vv{u}(t) = c_1\vv{x}_1 + c_2 e^{-3t}\vv{x}_2$ 에서 시간에 따라 변하는 부분이 $\vv{x}_2$ 방향뿐이다. 그러니 궤적이 $\vv{x}_2 = (1,-1)$ 방향의 직선일 수밖에 없다.

자료 = []
격자 = np.linspace(0, 4, 200)
for u초기 in ([1.0, 0.0], [-0.5, 1.0], [1.4, 0.6], [-1.0, -0.3], [0.2, -1.1]):
    길 = solve_ivp(lambda t, u: A @ u, (0, 4), u초기,
                   t_eval=격자, rtol=1e-10).y
    자료.append(go.Scatter(x=길[0], y=길[1], mode="lines",
                          line=dict(color=COLORS["input"], width=3),
                          showlegend=False))
    자료.append(go.Scatter(x=[길[0, -1]], y=[길[1, -1]], mode="markers",
                          marker=dict(color="#2ca02c", size=10),
                          showlegend=False))

축 = np.array([-2.5, 2.5])
자료.append(go.Scatter(x=축 * 2, y=축, mode="lines", name="고유벡터 (2,1), λ=0",
                      line=dict(color="#2ca02c", width=4)))
자료.append(go.Scatter(x=축, y=-축, mode="lines", name="고유벡터 (1,-1), λ=-3",
                      line=dict(color="#d62728", width=2, dash="dash")))

go.Figure(data=자료, layout=layout2d("모든 궤적이 (1,-1) 방향의 직선", extent=2.2))

초록 직선 위의 점은 전부 정지점이다. $\lambda = 0$ 이라 움직일 이유가 없다. 파란 궤적들은 서로 나란하고, 각자 자기 총량 $u_1 + u_2$ 를 지킨 채 초록 직선에 안착한다.

5. 크기는 실수부만 본다¶

$\left|e^{\lambda t}\right| = e^{\operatorname{Re}(\lambda)t}$ 였다. 허수부를 아무리 키워도 크기가 변하지 않는지 확인하자.

t = 2.0
print(f"{'lambda':>16}{'|e^(lambda t)|':>18}{'e^(Re t)':>14}")
for lam in (-0.5 + 0j, -0.5 + 3j, -0.5 + 50j, 0 + 7j, 0.3 - 100j):
    크기 = abs(np.exp(lam * t))
    print(f"{str(lam):>16}{크기:>18.6f}{np.exp(lam.real * t):>14.6f}")

          lambda    |e^(lambda t)|      e^(Re t)
       (-0.5+0j)          0.367879      0.367879
       (-0.5+3j)          0.367879      0.367879
      (-0.5+50j)          0.367879      0.367879
              7j          1.000000      1.000000
      (0.3-100j)          1.822119      1.822119

허수부가 3이든 50이든 크기가 똑같다. 허수부는 돌리기만 한다.

경우 = {
    "안정 노드": np.array([[-1.0, -1.0], [0.0, -2.0]]),
    "안장점":    np.array([[0.0, 1.0], [1.0, 0.0]]),
    "안정 나선": np.array([[-0.4, -2.0], [2.0, -0.4]]),
    "중심":      np.array([[0.0, -1.0], [1.0, 0.0]]),
    "앵커 A":    A,
}
print(f"{'':>10}{'고윳값':>30}{'최대 실수부':>13}{'판정':>10}")
for 이름, M in 경우.items():
    값 = np.linalg.eigvals(M)
    m = 값.real.max()
    판정 = "안정" if m < -1e-12 else ("발산" if m > 1e-12 else "중립")
    보기 = ", ".join(f"{v.real:.1f}{v.imag:+.1f}i" if abs(v.imag) > 1e-9
                    else f"{v.real:.1f}" for v in 값)
    print(f"{이름:>10}{보기:>30}{m:>13.2f}{판정:>10}")

                                     고윳값       최대 실수부        판정
     안정 노드                    -1.0, -2.0        -1.00        안정
       안장점                     1.0, -1.0         1.00        발산
     안정 나선          -0.4+2.0i, -0.4-2.0i        -0.40        안정
        중심            0.0+1.0i, 0.0-1.0i         0.00        중립
      앵커 A                     0.0, -3.0         0.00        중립

앵커 행렬은 최대 실수부가 정확히 0이라 중립이다. 그래서 원점이 아니라 정상상태로 간다.

실수부를 움직여 보자¶

$A(s) = \begin{bmatrix} s & -2 \\ 2 & s\end{bmatrix}$ 의 고윳값은 $s \pm 2i$ 이다. 허수부는 2로 고정한 채 실수부만 $s$ 로 움직인다. 슬라이더를 끌어 보자.

값들 = np.round(np.linspace(-0.7, 0.7, 15), 2)
시간 = np.linspace(0, 6, 400)
프레임 = []
for s in 값들:
    M = np.array([[s, -2.0], [2.0, s]])
    묶음 = []
    for 각 in np.linspace(0, 2 * np.pi, 5)[:-1]:
        시작 = 1.4 * np.array([np.cos(각), np.sin(각)])
        길 = np.array([expm(M * τ) @ 시작 for τ in 시간])
        안쪽 = np.abs(길).max(axis=1) < 3.0
        끝 = int(np.argmax(~안쪽)) if (~안쪽).any() else len(길)
        묶음.append(go.Scatter(x=길[:끝, 0], y=길[:끝, 1], mode="lines",
                              line=dict(color=COLORS["input"], width=2.5),
                              showlegend=False))
    프레임.append(묶음)

배치 = layout2d("실수부 s 를 움직이면", extent=3.0)
배치["height"] = 560
slider_figure(프레임, 값들, 배치, prefix="s = ", initial=7)

$s < 0$ 이면 안으로 감기고, $s = 0$ 이면 닫힌 원이며, $s > 0$ 이면 밖으로 풀린다. 허수축을 넘는 순간 운명이 바뀐다. 도는 속도(허수부 2)는 내내 그대로이다.

6. 행렬 지수함수¶

$e^{At} = I + At + \frac{(At)^2}{2!} + \cdots$ 를 직접 더해 보자.

def 급수(M, 항수):
    """e^M 을 항수 개의 항까지만 더한다."""
    합 = np.eye(M.shape[0])
    항 = np.eye(M.shape[0])
    for k in range(1, 항수):
        항 = 항 @ M / k
        합 = 합 + 항
    return 합

정답 = expm(A)
print(f"{'항수':>6}{'expm 과의 최대 차이':>24}")
for 항수 in (5, 10, 15, 20, 25, 30):
    print(f"{항수:>6}{np.abs(급수(A, 항수) - 정답).max():>24.3e}")
print()
print(show_matrix(정답, "expm(A)  (t=1)"))

    항수           expm 과의 최대 차이
     5               8.835e-01
    10               8.489e-03
    15               6.151e-06
    20               8.354e-10
    25               3.086e-14
    30               1.776e-15

expm(A)  (t=1)
[  0.683   0.633 ]
[  0.317   0.367 ]

급수를 다 더하지 않아도 된다. 대각화가 있으면 $e^{At} = Se^{\Lambda t}S^{-1}$ 이다.

def 대각화지수(t):
    """S e^(Lambda t) S^-1."""
    return S @ np.diag(np.exp(np.diag(L) * t)) @ np.linalg.inv(S)


def 닫힌꼴(t):
    """서술 파트에서 손으로 얻은 e^(At)."""
    e = np.exp(-3 * t)
    return np.array([[2 + e, 2 - 2 * e],
                     [1 - e, 1 + 2 * e]]) / 3

for t in (0.0, 0.7, 1.0, 3.0):
    a, b, c2 = expm(A * t), 대각화지수(t), 닫힌꼴(t)
    print(f"t = {t:>4} :  expm vs S e^L S^-1 {np.abs(a-b).max():.2e},"
          f"   expm vs 손으로 {np.abs(a-c2).max():.2e}")
print()
print(show_matrix(닫힌꼴(0.0), "t = 0  ->  I 인가"))
print(show_matrix(expm(A * 200), "t 가 크면  ->  랭크 1"))
print("그 열 :", expm(A * 200)[:, 0], "  = 정상상태")

t =  0.0 :  expm vs S e^L S^-1 5.55e-17,   expm vs 손으로 0.00e+00
t =  0.7 :  expm vs S e^L S^-1 3.89e-16,   expm vs 손으로 3.33e-16
t =  1.0 :  expm vs S e^L S^-1 1.78e-15,   expm vs 손으로 1.78e-15
t =  3.0 :  expm vs S e^L S^-1 5.55e-17,   expm vs 손으로 1.11e-16

t = 0  ->  I 인가
[  1   0 ]
[  0   1 ]
t 가 크면  ->  랭크 1
[  0.667   0.667 ]
[  0.333   0.333 ]
그 열 : [0.6667 0.3333]   = 정상상태

함정 하나 — 성분별 지수함수가 아니다¶

print(show_matrix(expm(A), "expm(A)  <- 맞는 계산"))
print(show_matrix(np.exp(A), "np.exp(A)  <- 성분별. 전혀 다르다"))
print("두 결과의 최대 차이 :", np.abs(expm(A) - np.exp(A)).max())

expm(A)  <- 맞는 계산
[  0.683   0.633 ]
[  0.317   0.367 ]
np.exp(A)  <- 성분별. 전혀 다르다
[  0.368    7.39 ]
[   2.72   0.135 ]
두 결과의 최대 차이 : 6.755580811175895

함정 둘 — $e^Ae^B = e^{A+B}$ 는 교환할 때만¶

P = np.array([[0.0, 1.0], [0.0, 0.0]])
Q = np.array([[0.0, 0.0], [1.0, 0.0]])
print("교환자 PQ - QP =\n", P @ Q - Q @ P)
print()
print(show_matrix(expm(P) @ expm(Q), "e^P e^Q"))
print(show_matrix(expm(P + Q), "e^(P+Q)"))
print("차이 :", np.abs(expm(P) @ expm(Q) - expm(P + Q)).max())

교환자 PQ - QP =
 [[ 1.  0.]
 [ 0. -1.]]

e^P e^Q
[  2   1 ]
[  1   1 ]
e^(P+Q)
[  1.54   1.18 ]
[  1.18   1.54 ]
차이 : 0.5430806348152437

# 교환하는 짝이면 성립한다
R = np.array([[2.0, 1.0], [0.0, 2.0]])
T = np.array([[5.0, 3.0], [0.0, 5.0]])
print("교환자 RT - TR =\n", R @ T - T @ R)
print("e^R e^T = e^(R+T) 인가 :", np.allclose(expm(R) @ expm(T), expm(R + T)))
print()
print("A 는 자기 자신과 늘 교환한다 :")
print("  e^(A*1) e^(A*2) = e^(A*3) 인가 :",
      np.allclose(expm(A) @ expm(A * 2), expm(A * 3)))
print("  e^(At) 의 역이 e^(-At) 인가 :",
      np.allclose(np.linalg.inv(expm(A * 1.3)), expm(-A * 1.3)))

교환자 RT - TR =
 [[0. 0.]
 [0. 0.]]
e^R e^T = e^(R+T) 인가 : True

A 는 자기 자신과 늘 교환한다 :
  e^(A*1) e^(A*2) = e^(A*3) 인가 : True
  e^(At) 의 역이 e^(-At) 인가 : True

7. 동반행렬 — 두 특성방정식이 같다¶

$y'' + by' + ky = 0$ 을 $\vv{u} = (y', y)$ 로 묶으면 $A = \begin{bmatrix} -b & -k \\ 1 & 0\end{bmatrix}$ 이 된다.

def 동반행렬(b, k):
    """y'' + b y' + k y = 0 의 동반행렬."""
    return np.array([[-b, -k], [1.0, 0.0]])

for b, k in ((5, 6), (2, 5), (0, 4), (2, 1)):
    C = 동반행렬(b, k)
    print(f"y'' + {b}y' + {k}y = 0")
    print(f"   미적분의 특성방정식 계수 : [1, {b}, {k}]")
    print(f"   det(A - lambda I) 의 계수 : {np.round(np.poly(C), 10)}")
    print(f"   근 : {np.sort_complex(np.linalg.eigvals(C))}")

y'' + 5y' + 6y = 0
   미적분의 특성방정식 계수 : [1, 5, 6]
   det(A - lambda I) 의 계수 : [1. 5. 6.]
   근 : [-3.+0.j -2.+0.j]
y'' + 2y' + 5y = 0
   미적분의 특성방정식 계수 : [1, 2, 5]
   det(A - lambda I) 의 계수 : [1. 2. 5.]
   근 : [-1.-2.j -1.+2.j]
y'' + 0y' + 4y = 0
   미적분의 특성방정식 계수 : [1, 0, 4]
   det(A - lambda I) 의 계수 : [1. 0. 4.]
   근 : [0.-2.j 0.+2.j]
y'' + 2y' + 1y = 0
   미적분의 특성방정식 계수 : [1, 2, 1]
   det(A - lambda I) 의 계수 : [1. 2. 1.]
   근 : [-1.+0.j -1.+0.j]

계수가 글자 하나 다르지 않다. 미적분에서 외운 특성방정식이 동반행렬의 특성방정식이었다.

고유벡터도 확인하자. $(A - \lambda I)$ 의 아래 줄이 $x_1 - \lambda x_2 = 0$ 이므로 고유벡터가 $(\lambda, 1)$ 이어야 한다.

C = 동반행렬(5, 6)
for lam in np.sort(np.linalg.eigvals(C).real):
    x = np.array([lam, 1.0])
    print(f"lambda = {lam:>5.1f} :  (A - lambda I) (lambda, 1) = {(C - lam*np.eye(2)) @ x}")

lambda =  -3.0 :  (A - lambda I) (lambda, 1) = [-0.  0.]
lambda =  -2.0 :  (A - lambda I) (lambda, 1) = [-0.  0.]

# y'' + 5y' + 6y = 0,  y(0) = 1,  y'(0) = 0
C = 동반행렬(5, 6)
u초기 = np.array([0.0, 1.0])                      # (y', y)
격자 = np.linspace(0, 3, 300)
수치 = solve_ivp(lambda t, u: C @ u, (0, 3), u초기, t_eval=격자,
                 rtol=1e-12, atol=1e-14).y

Sc = np.array([[-2.0, -3.0], [1.0, 1.0]])         # 열이 (lambda, 1)
cc = np.linalg.solve(Sc, u초기)
y공식 = cc[0] * np.exp(-2 * 격자) + cc[1] * np.exp(-3 * 격자)

print("c =", cc, "  ->  y(t) =", f"{cc[0]:.0f} e^(-2t) - {abs(cc[1]):.0f} e^(-3t)")
print("수치적분과의 최대 차이 :", np.abs(수치[1] - y공식).max())

c = [ 3. -2.]   ->  y(t) = 3 e^(-2t) - 2 e^(-3t)
수치적분과의 최대 차이 : 5.945244296867713e-14

중근이면 동반행렬이 결함 행렬이 된다¶

for b, k, 이름 in ((5, 6, "서로 다른 두 근"), (2, 1, "중근 lambda = -1")):
    C = 동반행렬(b, k)
    값 = np.linalg.eigvals(C)
    _, V = np.linalg.eig(C)
    print(f"{이름:>18} : 고윳값 {np.round(값.real, 3)}, "
          f"고유벡터 행렬의 랭크 {np.linalg.matrix_rank(V)}")

C = 동반행렬(2, 1)
print()
print(show_matrix(C + np.eye(2), "A + I  (lambda = -1)"))
print("랭크 :", np.linalg.matrix_rank(C + np.eye(2)),
      " -> 영공간이 1차원. 고유벡터가 하나뿐이다")

         서로 다른 두 근 : 고윳값 [-3. -2.], 고유벡터 행렬의 랭크 2
    중근 lambda = -1 : 고윳값 [-1. -1.], 고유벡터 행렬의 랭크 1

A + I  (lambda = -1)
[  -1   -1 ]
[   1    1 ]
랭크 : 1  -> 영공간이 1차원. 고유벡터가 하나뿐이다

고유벡터가 하나뿐이라 $c_1e^{-t} + c_2e^{-t}$ 로는 초기조건 두 개를 맞출 수 없다. 그래서 $te^{-t}$ 가 필요해진다. 왜 하필 $t$ 인지는 L28에서 답한다.

C = 동반행렬(2, 1)
격자 = np.linspace(0, 6, 300)
u초기 = np.array([0.0, 1.0])                      # y(0)=1, y'(0)=0
수치 = solve_ivp(lambda t, u: C @ u, (0, 6), u초기, t_eval=격자,
                 rtol=1e-12, atol=1e-14).y
후보 = np.exp(-격자) * (1 + 격자)                  # y = (1 + t) e^(-t)

print("y = (1 + t) e^(-t) 와의 최대 차이 :", np.abs(수치[1] - 후보).max())
go.Figure(
    data=[go.Scatter(x=격자, y=수치[1], mode="lines", name="수치적분",
                     line=dict(color=COLORS["input"], width=5)),
          go.Scatter(x=격자, y=후보, mode="lines", name="(1 + t) e^(-t)",
                     line=dict(color=COLORS["output"], width=2, dash="dash")),
          go.Scatter(x=격자, y=np.exp(-격자), mode="lines", name="e^(-t) 만으로는 부족",
                     line=dict(color="#999999", width=2, dash="dot"))],
    layout=go.Layout(title=dict(text="중근일 때는 t 가 곱해진 항이 필요하다"),
                     xaxis=dict(title=dict(text="t")),
                     yaxis=dict(title=dict(text="y")),
                     height=420, margin=dict(l=70, r=20, t=60, b=50)))

y = (1 + t) e^(-t) 와의 최대 차이 : 6.27831120425526e-14

마치며...¶

서술 파트의 내용	이 노트북의 코드
대입하면 고윳값이 나온다	`잔차` — 고유쌍일 때만 0
일반해와 $\vv{c}$	`손으로푼해` 와 `solve_ivp` 의 차이 10^-10 수준
열의 합이 0	$\vv{1}^{\mathsf T}A = \vv{0}$ , $u_1 + u_2$ 가 상수
궤적이 직선	위상평면. 초록 직선 전체가 정지점
크기는 실수부만	허수부를 50으로 키워도 $\lvert e^{\lambda t}\rvert$ 가 그대로
실수부가 운명	슬라이더 $s$ 가 허수축을 넘는 순간
$e^{At}$ 급수	20항이면 10^-15 까지 맞는다
$Se^{\Lambda t}S^{-1}$	급수와 손으로 푼 닫힌 꼴이 모두 일치
함정 둘	성분별 지수함수, $e^Ae^B \neq e^{A+B}$
동반행렬	`np.poly` 의 계수가 미적분의 특성방정식과 같다
중근	동반행렬의 랭크가 모자라고 $(1+t)e^{-t}$ 가 답

더 해 볼 것¶

5절의 슬라이더에서 허수부 2를 0.5로 바꿔 보자. 감기는 모양이 어떻게 달라지는가? 운명(안정인지 발산인지)은 달라지는가?
잔차 를 복소 고윳값을 가진 행렬에 써 보자. 잔차가 여전히 0인가? 그 복소해 두 개를 더하면 실수가 되는가?
앵커 행렬의 비율 1과 2를 다른 값으로 바꿔 보자. 정상상태가 어떻게 움직이는가? 두 비율의 비와 어떤 관계가 있는가?
3계 미분방정식 $y''' + 6y'' + 11y' + 6y = 0$ 의 동반행렬을 $3 \times 3$ 으로 만들어 보자. np.poly 의 계수가 $[1, 6, 11, 6]$ 인가?

다음 강의에서는 시간을 다시 이산적으로 바꾼다. 이번 강의의 $A$ 는 열의 합이 0이었는데, 마코브 행렬은 열의 합이 1이다. 그 한 글자 차이가 고윳값 0을 고윳값 1로 바꾼다.

Lecture 23. 미분방정식과 exp(At) — 파이썬 실습