Chapter 04

변화율로 세상을 쓰다

세상의 많은 법칙은 "지금 값이 얼마인가"가 아니라 "지금 얼마나 빨리 변하고 있는가"로 주어집니다. 열은 온도 차에 비례해서 흐르고, 축전기 전류는 전압 차에 비례하고, 감염은 감염자와 비감염자가 만나는 만큼 늘어납니다. 변화율로 적힌 법칙에서 미래를 읽어 내는 도구, 그것이 미분방정식(Differential Equation)입니다. 이 장에서는 서로 전혀 달라 보이는 물건들이 사실 같은 식 몇 개로 움직인다는 것을 직접 확인합니다.

칩 온도·RC 충전·방사성 붕괴·약물 농도가 모두 $\dot y = -(y-y_\infty)/\tau$ 하나임을 보이고, 시정수 τ로 63%·5τ 규칙을 쓴다.
미분방정식을 "모든 점에 화살표를 그린 지도"(기울기장)로 읽고, 풀이 공식 없이도 해의 모양을 예측한다.
오일러법과 RK4가 어떻게 해를 근사하는지, 단계 크기 h가 오차와 안정성을 어떻게 좌우하는지 안다.
스프링–질량–감쇠기와 RLC 회로가 같은 2차 방정식임을 알고, 감쇠비 ζ와 특성근의 위치로 응답을 읽는다.
공진과 Q, SIR 감염병 모델의 집단면역 임계 $1-1/R_0$, 열확산 방정식과 그 안정 조건을 설명한다.

이 수식이 없었다면스마트폰 AP의 열 관리(쓰로틀링)

게임을 켜면 스마트폰 AP는 수 와트를 쏟아냅니다. 설계자가 알고 싶은 것은 "몇 초 뒤에 접합 온도가 한계에 닿는가, 그때 클록을 얼마나 낮춰야 하는가"입니다. 그런데 손에 쥔 정보는 온도 자체가 아니라 온도가 변하는 속도를 정하는 규칙뿐입니다. 들어오는 전력과 빠져나가는 열의 차이가 온도 상승률을 정한다는 것. 이 규칙을 미분방정식으로 쓰고 풀면 온도 곡선 전체가 나오고, 그 곡선이 '열 시정수'라는 숫자 하나로 요약됩니다. 휴대폰의 쓰로틀링 정책, 서버의 팬 제어, 패키지의 방열판 크기가 모두 이 계산 위에 서 있습니다.

칩 발열 예측→온도는 모르고 변화율만 앎→미분방정식 · 시정수 τ→온도 곡선 전체

하나의 식, 네 개의 세계

칩은 몇 초 만에 뜨거워지는가

AP가 일정 전력 $P$를 소비하기 시작했다. 칩의 열용량은 $C_{th}$(J/K), 칩에서 공기까지 열이 빠져나가는 길의 열저항은 $R_{th}$(K/W)다. 쓰로틀링 한계 80 °C에 닿는 시각을 알고 싶다. 그런데 물리 법칙이 알려 주는 것은 "지금 온도"가 아니라 "지금 온도가 오르는 속도"다.

에너지 보존을 그대로 쓰면 됩니다. 칩에 쌓이는 열(열용량 × 온도 상승률)은 들어오는 전력에서 빠져나가는 열류를 뺀 것입니다. 빠져나가는 열류는 주변과의 온도 차에 비례합니다(뉴턴의 냉각 법칙이라고 불리는 근사).

$$C_{th}\frac{dT}{dt} = P - \frac{T - T_{amb}}{R_{th}} \quad\Longrightarrow\quad \frac{dT}{dt} = -\frac{T - T_\infty}{\tau},\qquad \tau = R_{th}C_{th},\; T_\infty = T_{amb} + P R_{th}$$

$T_\infty$: 충분히 오래 기다렸을 때의 정상 온도. $\tau$: 열 시정수(thermal time constant). 식의 오른쪽은 "목표까지 남은 거리 ÷ τ" — 멀리 있을수록 빨리 다가가고, 가까워질수록 느려진다.

그림 4-1. 칩 열 모델의 등가 회로. 발열 $P$는 전류원, 열용량은 축전기, 방열 경로는 저항이 된다. 이 회로의 식은 RC 회로와 기호만 다를 뿐 똑같다. 실제 패키지는 다이–기판–섀시로 여러 단의 RC가 이어진 '열 회로망'(포스터·카우어 모델)으로 다룬다.

이 식의 해는 한 번만 구해 두면 됩니다. $u = T - T_\infty$로 두면 $\dot u = -u/\tau$, 즉 "자기 크기에 비례해서 줄어드는 양"이고, 그런 함수는 지수함수뿐입니다.

$$y(t) = y_\infty + \left(y_0 - y_\infty\right)e^{-t/\tau}$$

$t=\tau$에서 남은 거리는 $e^{-1}\approx 36.8\%$, 즉 목표의 63.2%까지 왔다. $t = 5\tau$에서는 $e^{-5}\approx 0.67\%$만 남아 엔지니어는 이를 "정착했다"고 본다. 반감기는 $\tau\ln 2 \approx 0.693\,\tau$.

수치 예를 들어 봅시다(설명용 값). $T_{amb}=25$ °C, $P = 5$ W, $R_{th}=12$ K/W, $C_{th}=5$ J/K이면 $T_\infty = 85$ °C, $\tau = 60$ s입니다. 80 °C에 닿는 시각은 $e^{-t/\tau} = 5/60$에서 $t = 60\ln 12 \approx 149$ 초. 정상 온도가 한계를 넘는다는 사실(85 > 80)은 $T_\infty$만 봐도 알 수 있고, 언제 넘는지는 τ가 알려 줍니다. 방열판을 키우면 $R_{th}$가 줄어 $T_\infty$가 내려가고, 금속 덩어리(열 확산판)를 붙이면 $C_{th}$가 커져 같은 $T_\infty$로 더 천천히 갑니다 — 짧은 버스트 부하에는 후자가 효과적입니다.

같은 곡선, 다른 이름표

놀라운 것은 이 식이 칩에만 해당하지 않는다는 점입니다. "목표와의 차이에 비례해서 변한다"는 구조를 가진 현상은 모두 같은 곡선을 그립니다.

현상	y	y∞ (목표)	τ
RC 회로 충전 (카메라 플래시 축전기)	축전기 전압 $V_C$	전원 전압 $V_s$	$RC$
칩·커피 냉각/가열 (뉴턴 냉각)	온도 $T$	$T_{amb}+PR_{th}$	$R_{th}C_{th}$
방사성 붕괴	남은 핵 수 $N$	0	평균 수명 $=t_{1/2}/\ln2$
일정 속도 정맥 주입 (1구획 약동학)	혈중 농도 $C$	주입 속도 / 청소율	분포 용적 / 청소율
소수 캐리어 재결합 (반도체)	과잉 캐리어 $\Delta n$	0 (또는 생성률 × τ)	캐리어 수명

SIMULATOR

해석 렌즈: 같은 곡선에 이름표만 바꾸기

해석 렌즈

시정수 τ 시작점의 접선 보이기

이 렌즈의 식

—

τ (물리 단위)—

반감기 τ·ln2—

5τ (99.3%)—

t = τ에서의 값—

왼쪽 점(시작값 $y_0$)과 오른쪽 점(목표 $y_\infty$)을 위아래로 끌 수 있습니다. 해볼 것: ① 렌즈를 바꿔 보세요 — 축의 단위와 이름만 바뀌고 곡선 모양은 그대로입니다. ② 시작점의 접선(점선)이 언제나 정확히 t = τ에서 목표선과 만나는지 확인. ③ $y_0$를 목표보다 위로 끌면 '충전'이 '방전·냉각'이 되지만 63% 규칙은 그대로. · 모델: 선형 1차 계(열저항·청소율·R·C가 상수). 실제 칩은 다단 RC, 약물은 다구획 모델이 더 정확하다. 단위 환산용 기준값(RC 1 s, 칩 40 s, I-131 반감기 약 8일, 약물 4 h)은 설명용 예시다.

τ만 알면 된다 1차 계의 응답은 $y_0, y_\infty, \tau$ 세 숫자로 완전히 정해진다. 그래서 데이터시트는 열 저항과 열 시정수를, 의약품 설명서는 반감기를, 회로도는 RC 값을 적어 둔다. 오실로스코프로 RC를 잴 때 "최종값의 63%에 닿는 시각"을 읽는 것도 같은 이유다.

기울기장: 모든 점에 화살표를 그린 지도

공식으로 안 풀리는 식

개체 수가 자원 한계 $K$에 가까워지면 성장이 둔해지는 로지스틱 식 $\dot y = r y(1-y/K)$이나, 부하가 주기적으로 변하는 칩의 온도 $\dot T = -(T - T_{amb})/\tau + P(t)/C_{th}$ 같은 식은 공식을 외워서 풀기 어렵거나 아예 닫힌 해가 없다. 그래도 해가 어떤 모양일지는 알고 싶다.

미분방정식 $\dot y = f(t, y)$를 다시 읽어 봅시다. 이 식은 "$(t, y)$ 지점에 있다면 기울기는 $f(t,y)$다"라고 말합니다. 즉 평면의 모든 점에 화살표 하나씩을 지정합니다. 이 화살표 지도가 기울기장(Direction Field, Slope Field)이고, 해곡선은 화살표를 따라 흐르는 물길입니다. 해를 구한다는 것은 "출발점을 정하면 물길이 정해진다"는 사실을 이용하는 일입니다.

SIMULATOR

기울기장 위에 해곡선 흘려 보내기

—

방정식 기울기 0인 곡선(널클라인) 보기 화살표를 기울기 크기로 색칠

드래그 점 (t, y)—

그 점의 기울기 f(t,y)—

그린 해곡선 수—

빈 곳을 클릭(터치)하면 그 점에서 출발하는 해곡선이 앞뒤로 흘러갑니다. 큰 점은 끌 수 있는 '현재 상태'입니다. 해볼 것: ① 로지스틱에서 y = 0과 y = 4 근처에 점을 놓아 보기 — 하나는 밀어내고(불안정) 하나는 끌어당긴다(안정). ② y′ = y − t에서 직선 y = t + 1 바로 위·아래에 놓으면 해가 정반대로 갈라진다. ③ 주기 부하에서 출발점이 어디든 결국 같은 진동으로 모이는지(과도 응답 → 정상 응답). · 모델: 각 해곡선은 RK4로 작은 간격으로 적분했다. 해가 범위 밖으로 발산하면 그 지점에서 멈춘다.

시뮬레이터에서 본 것들을 정리하면, 공식 없이도 해에 대해 많은 것을 말할 수 있습니다.

평형점(equilibrium): $f(y^*) = 0$인 수평선. 화살표가 양쪽에서 모이면 안정, 흩어지면 불안정. 로지스틱의 $y=K$는 안정, $y=0$은 불안정. 1차 선형계에서 $y_\infty$가 언제나 안정 평형인 것은 τ > 0이기 때문이다.
유일성: $f$가 매끄러우면 해곡선은 서로 교차하지 않는다(피카르–린델뢰프 정리). 그래서 두 곡선 사이에 갇힌 해는 영원히 그 사이에 있다 — 공식 없이 해의 범위를 가두는 강력한 논리다.
과도 응답과 정상 응답: 감쇠가 있는 계에 주기 입력이 들어오면, 처음 상태의 기억은 $e^{-t/\tau}$로 사라지고 모든 해가 하나의 주기 해로 모인다. 칩 온도가 부하 주기를 따라 출렁이는 것이 이것이다. 진폭과 위상 지연은 2장의 페이저로 바로 계산된다.

로지스틱 식의 해 $\dot y = ry(1-y/K)$는 변수분리로 풀린다: $y(t) = K/\big(1 + (K/y_0 - 1)e^{-rt}\big)$. 처음에는 지수 성장($y \ll K$이면 $\dot y\approx ry$), 중간 $y = K/2$에서 기울기 최대, 끝에서는 $K$로 1차 계처럼 다가간다. 기술 보급 곡선(S-커브), 반도체 결함 밀도 학습 곡선의 경험식이 이 모양을 빌려 쓴다. 6절의 감염병 모델에서 다시 나온다.

컴퓨터는 어떻게 푸는가: 오일러와 RK4

회로 시뮬레이터 안에서 일어나는 일

SPICE, 열 해석, 차량 동역학 시뮬레이터는 손으로 풀 수 없는 수천 개짜리 연립 미분방정식을 매 순간 푼다. 컴퓨터가 할 수 있는 것은 "이 점에서 기울기 $f(t,y)$를 계산해 줘"뿐이다. 기울기만으로 어떻게 곡선을 따라가는가? 한 걸음을 얼마나 크게 내디뎌도 되는가?

가장 정직한 방법은 기울기장의 화살표를 따라 한 걸음씩 걷는 것입니다. 지금 점의 기울기로 $h$만큼 직선으로 나아가고, 도착한 점에서 다시 기울기를 묻습니다. 이것이 오일러법(Euler Method)으로, 오일러가 1768년의 적분학 교재에서 제시했습니다.

$$y_{n+1} = y_n + h\, f(t_n, y_n)$$

한 걸음의 오차는 $O(h^2)$(곡선의 휘어짐을 무시한 만큼), 걸음 수가 $T/h$이므로 끝점 오차는 $O(h)$ — 걸음을 절반으로 줄이면 오차도 절반. 1차 방법이다.

그림 4-2. 왼쪽: 오일러법은 출발점의 접선만 믿고 걷기 때문에 곡선이 휘는 쪽 바깥으로 오차가 쌓인다. 오른쪽: 룽게–쿠타 4차(RK4)는 한 걸음 안에서 기울기를 네 번 물어($k_1$~$k_4$) 심프슨 공식처럼 1:2:2:1로 평균한다.

더 나은 방법은 한 걸음 안에서 기울기를 여러 번 물어보는 것입니다. 룽게(1895)와 쿠타(1901)가 발전시킨 방법 중 가장 널리 쓰이는 RK4(classical Runge–Kutta)는 이렇습니다.

$$\begin{aligned}k_1 &= f(t_n, y_n), & k_2 &= f\!\left(t_n+\tfrac h2,\; y_n + \tfrac h2 k_1\right),\\ k_3 &= f\!\left(t_n+\tfrac h2,\; y_n + \tfrac h2 k_2\right), & k_4 &= f(t_n + h,\; y_n + h k_3),\end{aligned}\qquad y_{n+1} = y_n + \tfrac h6\left(k_1 + 2k_2 + 2k_3 + k_4\right)$$

끝점 오차 $O(h^4)$: 걸음을 절반으로 줄이면 오차가 1/16. 걸음당 계산은 4배지만, 같은 정확도를 훨씬 큰 걸음으로 얻는다. 이 책의 모든 시뮬레이터(MB.rk4)도 이것을 쓴다.

SIMULATOR

오일러 vs RK4: 단계 크기 h와 오차·불안정

해 비교

끝점 오차 vs h (로그–로그)

문제

단계 크기 h RK4도 함께 보기

걸음 수 / f 호출(오일러·RK4)—

오일러 끝점 오차—

RK4 끝점 오차—

오일러 한 걸음 증폭률—

—

해볼 것: ① 감쇠 문제에서 h를 1 → 2 → 2.5로 올려 보세요. h > 1이면 오일러 해가 부호를 바꾸며 진동하고, h > 2이면 0으로 가야 할 해가 폭발합니다. ② 진동 문제에서는 h가 아무리 작아도 오일러 해가 바깥으로 나선형으로 커집니다(에너지가 생김). ③ 오른쪽 그래프에서 두 직선의 기울기 — 오일러는 1, RK4는 4 — 를 확인. · 모델: 감쇠는 $y(0)=1$, 진동은 $x(0)=1, \dot x(0)=0$. 오차는 감쇠 $t\approx 2$, 진동 $t\approx 20$(정수 걸음으로 가장 가까운 시각)에서 정확해와의 차이.

두 번째 관찰이 더 중요합니다. 정확도와 별개로 안정성(Stability)의 벽이 있습니다. $\dot y = -\lambda y$에 오일러를 쓰면 $y_{n+1} = (1 - h\lambda)\,y_n$이므로, 증폭률 $|1-h\lambda|$가 1을 넘는 순간($h\lambda > 2$) 해가 폭발합니다. 시정수가 짧은 성분(큰 λ) 하나가 전체의 걸음 크기를 묶어 버리는 이 현상을 강성(stiffness)이라 합니다. 1 ps짜리 기생 RC와 1 ms짜리 열 시정수가 한 회로에 섞여 있으면, 명시적 방법은 1 ps 걸음으로 1 ms를 걸어야 합니다. 그래서 SPICE 같은 회로 시뮬레이터는 후진 오일러·사다리꼴·기어(BDF) 같은 암시적 방법을 기본으로 씁니다. 이 이야기는 10장에서 이어집니다.

2차 방정식: 서스펜션과 RLC는 같은 식이다

바늘이 흔들리지 않고 빨리 멈추게

아날로그 계기의 바늘, 자동차 서스펜션, 하드디스크 헤드, 카메라 손떨림 보정 액추에이터. 모두 "목표 위치로 빨리 가되, 지나쳐서 출렁이지 말라"는 요구를 받는다. 스프링이 너무 세면 출렁이고, 댐퍼가 너무 세면 굼뜨다. 그 사이 어디가 최적인가?

질량 $m$이 스프링(강성 $k$)과 댐퍼(감쇠 계수 $c$)에 매달려 있으면 뉴턴의 제2법칙은 2차 미분방정식이 됩니다. 놀랍게도 저항·인덕터·축전기를 직렬로 이은 회로에서 전하 $q$가 따르는 식도 글자만 다르고 똑같습니다.

그림 4-3. 기계–전기 상사(相似). 질량은 관성(인덕턴스), 댐퍼는 에너지를 열로 버리는 저항, 스프링은 에너지를 저장했다 돌려주는 축전기(강성은 1/C)에 대응한다. 그래서 서스펜션 설계와 필터 설계가 같은 언어를 쓴다.

$$\ddot x + 2\zeta\omega_0\,\dot x + \omega_0^2\,x = 0,\qquad \omega_0 = \sqrt{\frac km} = \frac{1}{\sqrt{LC}},\qquad \zeta = \frac{c}{2\sqrt{mk}} = \frac R2\sqrt{\frac CL}$$

$\omega_0$: 감쇠가 없을 때의 고유 각진동수. $\zeta$: 감쇠비(damping ratio). 계의 성격은 이 두 숫자로 전부 정해진다. 필터 설계에서는 $Q = 1/(2\zeta)$를 더 자주 쓴다.

이 식을 푸는 요령은 $x = e^{st}$를 넣어 보는 것입니다. 지수함수는 미분해도 모양이 그대로이므로 미분방정식이 대수 방정식으로 바뀝니다.

$$s^2 + 2\zeta\omega_0 s + \omega_0^2 = 0 \quad\Longrightarrow\quad s_{1,2} = -\zeta\omega_0 \pm \omega_0\sqrt{\zeta^2 - 1}$$

이 특성방정식의 근이 어디 있느냐가 전부다. $\zeta>1$: 서로 다른 음의 실근 둘 → 진동 없이 느리게(과감쇠). $\zeta=1$: 중근 → 진동 없이 가장 빠르게(임계 감쇠). $\zeta<1$: 켤레 복소근 $-\sigma \pm i\omega_d$ → $e^{-\sigma t}$로 줄어드는 진동(부족 감쇠).

복소근이 왜 진동일까요? 2장의 오일러 공식으로 $e^{(-\sigma + i\omega_d)t} = e^{-\sigma t}\,(\cos\omega_d t + i\sin\omega_d t)$. 실수부는 얼마나 빨리 줄어드는가, 허수부는 얼마나 빨리 도는가입니다. 복소평면의 점 하나가 응답의 성격을 다 말해 줍니다.

SIMULATOR

감쇠비 ζ와 특성근: 계단 응답 ↔ 복소평면

계단 응답 x(t)

특성근 s (복소평면)

감쇠비 ζ 고유 각진동수 ω₀

해석

특성근 s₁, s₂—

종류—

오버슈트—

2% 정착 시간—

오른쪽 복소평면의 점(위쪽 근)을 직접 끌 수 있습니다. 해볼 것: ① 근을 허수축 쪽으로 끌면(ζ → 0) 감쇠 없는 진동, 왼쪽으로 끌면 빨리 죽는다 — 실수부 = 감쇠 속도. ② 원(|s| = ω₀)을 따라 돌리면 ω₀는 그대로 ζ만 바뀐다. ③ '임계 감쇠로'를 누른 뒤 ζ를 1보다 키우면 근이 실수축에서 둘로 갈라지고, 느린 근이 원점 쪽으로 다가가 응답이 오히려 굼떠진다. · 모델: 선형 2차 계, 단위 계단 입력(목표 위치 0 → 1). 실제 서스펜션은 타이어 강성·비선형 댐퍼가 더해진 다자유도 계다. 승용차 서스펜션은 대략 ζ ≈ 0.2~0.4, 계측기는 0.6~0.7 근처로 설계하는 경우가 많다.

왜 계측기는 ζ ≈ 0.7인가 임계 감쇠(ζ = 1)는 '오버슈트 없이 가장 빠른' 응답이지만, 약간의 오버슈트(ζ = 0.7에서 약 4.6%)를 허용하면 목표 근처(±5%)에 더 일찍 들어온다. 버터워스 2차 필터가 $Q = 1/\sqrt2$, 즉 $\zeta = 0.707$인 것도 같은 타협이다 — 주파수 응답이 공진 봉우리 없이 가장 평평하다.

강제 진동과 공진: 작은 힘이 큰 진폭을 만드는 이유

엔진 회전수가 특정 값일 때만 떨리는 차

특정 rpm에서만 대시보드가 웅웅 떨리고, 세탁기는 탈수 초반 특정 속도를 지날 때 크게 흔들린다. 라디오는 수많은 방송 중 하나만 골라 듣는다. 모두 외부에서 주기적으로 미는 힘이 계의 고유 진동수와 맞을 때 생기는 일이다. 얼마나 커지고, 무엇이 그 크기를 정하는가?

2차 계를 $F_0\sin\omega t$로 계속 밀면, 과도 응답이 $e^{-\zeta\omega_0 t}$로 사라진 뒤 같은 주파수로 진동하는 정상 응답만 남습니다. 그 진폭을 정적 처짐 $F_0/k$로 나눈 값(동적 증폭률)은 주파수비 $r = \omega/\omega_0$의 함수입니다.

$$\left|\frac{X}{F_0/k}\right| = \frac{1}{\sqrt{(1-r^2)^2 + (2\zeta r)^2}},\qquad \varphi = \operatorname{atan2}\!\left(2\zeta r,\; 1 - r^2\right),\qquad Q = \frac{1}{2\zeta}$$

$r = 1$에서 증폭률은 정확히 $Q$이고 위상 지연은 90°. 봉우리 자체는 $r = \sqrt{1-2\zeta^2}$에서 약간 낮은 쪽에 생긴다(ζ가 작으면 거의 1). 반전력(−3 dB) 대역폭은 $\Delta\omega \approx \omega_0/Q$.

SIMULATOR

주파수 응답 곡선과 공진의 성장

동적 증폭률 |X|/(F₀/k) vs r = ω/ω₀

정지 상태에서 시작한 응답 x(t), f₀ = 1 Hz

감쇠비 ζ 구동 주파수비 r

Q = 1/(2ζ)—

정상 진폭(증폭률)—

위상 지연 φ—

−3 dB 대역폭(f₀ = 1 Hz)—

해볼 것: ① r = 1에 두고 ζ를 0.01까지 낮추면 진폭이 Q배(50배)까지, 그러나 천천히 자란다 — 정상 상태까지 약 Q/π 주기 이상이 걸린다. ② r = 2 이상에서는 증폭률이 1보다 작고 위상이 거의 180° — 질량이 힘과 반대로 움직인다(진동 절연의 원리). ③ ζ를 0.7 이상으로 하면 봉우리가 사라진다. · 모델: 선형 2차 계, $\ddot x + 2\zeta\omega_0\dot x + \omega_0^2 x = \omega_0^2\sin(r\omega_0 t)$, RK4로 적분. 점선은 정상 진폭.

진폭이 커지는 이유는 에너지의 장부로 보면 분명합니다. 힘이 한 주기 동안 계에 넣는 일은 $\oint F\,dx = \oint F\,\dot x\,dt$인데, 공진에서는 속도가 힘과 같은 위상(변위는 90° 지연)이라 매 주기 일이 최대로 들어갑니다. 진폭은 "매 주기 들어오는 에너지 = 댐퍼가 매 주기 버리는 에너지"가 될 때까지 자랍니다. 댐퍼가 버리는 에너지는 진폭의 제곱에 비례하므로, 감쇠가 작을수록 그 균형점은 높아집니다. 그래서 Q는 두 얼굴을 가집니다: 공진에서의 증폭률이자, 자유 진동이 얼마나 오래 울리는가(에너지가 $1/e$로 줄 때까지 약 $Q/2\pi$ 주기)입니다.

계	Q의 대략적 크기	공진을 어떻게 쓰나/피하나
자동차 서스펜션	1~3	바퀴 진동 주파수 대역에서 차체로 전달을 줄이도록 ζ 조정
RLC 동조 회로(라디오)	수십~수백	Q가 클수록 옆 채널을 잘 거른다(선택도)
수정 진동자	10⁴~10⁶	극도로 좁은 공진 → 시계·클록의 기준 주파수
MEMS 자이로·가속도계	수~수만(진공 패키징)	감지 모드 공진으로 감도 증폭, 대역폭과 맞바꿈

타코마 다리는 '단순 공진'이 아니다 1940년 타코마 내로스 다리 붕괴는 흔히 바람과의 공진으로 소개되지만, 바람은 일정한 주기로 밀지 않았다. 현재의 설명은 구조물의 움직임이 공기력을 바꾸고 그 공기력이 다시 움직임을 키우는 공탄성 플러터(자려 진동)다. 식으로 쓰면 유효 감쇠 $\zeta$가 음수가 된 경우로, 특성근이 복소평면의 오른쪽 반평면으로 넘어간 것이다. 앞 절의 시뮬레이터에서 근을 허수축 오른쪽으로 보내면 어떻게 될지 상상해 보라.

감염병 SIR 모델: 연립 미분방정식

백신은 몇 %가 맞아야 하는가

모두가 맞을 수는 없다. 그렇다면 몇 %가 면역을 가져야 유행이 스스로 꺼지는가? 접종률이 그보다 조금 모자라면 무슨 일이 생기는가? 그리고 정점은 언제 오는가 — 병상을 미리 준비하려면 그 날짜가 필요하다.

인구를 세 칸으로 나눕니다: 감염될 수 있는 $S$, 감염된 $I$, 회복(또는 면역)된 $R$. 감염은 $S$와 $I$가 만나는 만큼(질량 작용) 일어나고, 회복은 평균 $D$일 걸린다고 하면 1927년 커맥–맥켄드릭이 정식화한 SIR 모델이 됩니다.

그림 4-4. SIR 구획 모델. 각 화살표가 변화율 하나다. 칸 사이의 흐름만 정하면 연립 미분방정식이 저절로 써진다 — 화학 반응 속도론, 약동학 구획 모델, 반도체의 캐리어 생성·재결합 모델이 모두 같은 방식으로 만들어진다.

$$\dot S = -\beta S I,\qquad \dot I = \beta S I - \gamma I,\qquad \dot R = \gamma I,\qquad R_0 = \frac{\beta}{\gamma}$$

$S+I+R = 1$(비율). $R_0$: 기초 감염 재생산수 — 모두가 감수성일 때 감염자 한 명이 회복 전까지 옮기는 평균 인원. $\dot I = \gamma I\,(R_0 S - 1)$이므로 $R_0 S > 1$일 때만 감염이 늘어난다.

마지막 줄에서 두 가지 결론이 바로 나옵니다. 첫째, 감염자 곡선의 정점은 $S$가 $1/R_0$로 떨어지는 순간입니다(그때 $\dot I = 0$). 둘째, 처음부터 $S < 1/R_0$이면 감염은 처음부터 줄어듭니다. 접종률 $p$이면 $S_0 \approx 1-p$이므로

$$p_c = 1 - \frac{1}{R_0}$$

집단면역 임계. $R_0 = 2.5$면 60%, $R_0 = 4$면 75%. 홍역처럼 $R_0$가 10을 훌쩍 넘는 것으로 알려진 병은 90%를 넘어야 한다. 임계에 못 미친 접종도 정점을 낮추고 늦춘다.

SIMULATOR

SIR 모델: R₀, 회복 기간, 접종률

S · I · R (인구 비율) vs 일

위상 평면 (S, I)

기초 재생산수 R₀ 평균 회복 기간 D 접종률 p

세로축

실효 재생산수 R₀(1−p)—

집단면역 임계 1−1/R₀—

감염자 정점 / 날짜—

최종 누적 감염—

해볼 것: ① 접종률을 집단면역 임계 바로 아래·위로 움직여 보세요 — 임계를 넘으면 곡선이 일어서지 못한다. ② 오른쪽 위상 평면에서 I의 정점이 언제나 점선 S = 1/R₀ 위에 있음을 확인. ③ 로그 축으로 바꾸면 초기 감염자 곡선이 직선 — 지수 성장 $e^{(\beta-\gamma)t}$이다. 같은 R₀에서 D만 바꾸면 최종 감염 비율은 같고 시간 축만 늘어난다. · 모델: 잘 섞인 단일 집단, 출생·사망 없음, 영구 면역, 접종은 100% 효과, 초기 감염자 0.01%. 행동 변화·연령 구조·잠복기(SEIR)는 무시한 교육용 모델이다.

SIR과 로지스틱 성장의 관계도 보입니다. 회복이 없으면($\gamma = 0$) $S = 1 - I$이므로 $\dot I = \beta I (1-I)$ — 바로 2절의 로지스틱 식입니다. 회복이 있으면 감염자 수는 로지스틱처럼 포화하지 않고, 감수성 인구라는 '연료'가 $1/R_0$ 아래로 떨어지는 순간 꺾여서 줄어듭니다. 또 하나 눈여겨볼 점: 유행이 끝나도 $S$는 0이 되지 않습니다. 최종 감염 비율 $z$는 $z = 1 - e^{-R_0 z}$를 만족하며, $R_0=2.5$면 약 89%입니다. 집단면역 임계(60%)를 한참 넘어선 '관성 감염'이 생기는 셈입니다.

편미분방정식 맛보기: 열확산 방정식

칩 위의 핫스팟은 어떻게 퍼지는가

1절의 모델은 칩 전체를 온도 하나로 보았다. 그러나 실제 다이에는 CPU 코어 하나가 몰래 뜨거워지는 핫스팟이 있고, 온도 센서는 그 옆 몇 mm에 있다. 온도가 위치 $x$와 시간 $t$ 모두의 함수가 되면, 변화율도 두 방향으로 생긴다.

막대를 잘게 나눈 칸 하나를 생각합시다. 칸으로 들어오는 열은 왼쪽 이웃과의 온도 차, 나가는 열은 오른쪽 이웃과의 온도 차에 비례합니다(푸리에의 열전도 법칙). 둘을 합치면 칸의 온도 상승률은 $(T_{i-1} - T_i) + (T_{i+1} - T_i)$에 비례합니다. 정리하면 놀랍도록 단순한 문장이 됩니다.

그림 4-5. 2차 미분의 직관. $\frac{T_{i-1} - 2T_i + T_{i+1}}{\Delta x^2} = \frac{2}{\Delta x^2}\left(\frac{T_{i-1}+T_{i+1}}{2} - T_i\right)$. 2차 미분은 "이웃 평균에서 얼마나 튀어나와 있는가"를 잰다. 봉우리는 깎이고 골짜기는 메워진다 — 확산은 세상을 평평하게 만드는 연산이다.

$$\frac{\partial T}{\partial t} = \alpha\,\frac{\partial^2 T}{\partial x^2}\qquad\xrightarrow{\text{유한 차분}}\qquad T_i^{n+1} = T_i^n + r\left(T_{i-1}^n - 2T_i^n + T_{i+1}^n\right),\quad r = \frac{\alpha\,\Delta t}{\Delta x^2}$$

$\alpha = k/(\rho c_p)$: 열확산도(m²/s). 실리콘은 상온에서 약 0.8~0.9 cm²/s. 확산 거리는 $\sqrt{\alpha t}$로 자라므로, 실리콘에서 1 mm 퍼지는 데는 대략 $t \sim L^2/\alpha \approx 10$ ms 정도가 걸린다. 1822년 푸리에가 열의 해석적 이론에서 이 식을 풀려고 푸리에 급수를 꺼냈다(5장).

오른쪽 식은 1절의 오일러법을 칸마다 적용한 것뿐입니다. 그래서 같은 안정성 문제가 따라옵니다. 가장 들쭉날쭉한 모양(칸마다 +, − 번갈아)은 한 걸음에 $1 - 4r$배가 되므로, $|1-4r| \le 1$, 즉

$$r = \frac{\alpha\,\Delta t}{\Delta x^2} \le \frac12$$

격자를 두 배 촘촘히($\Delta x/2$) 하면 시간 걸음은 네 배 짧아야 한다. 명시적 확산 계산이 비싼 이유이고, TCAD·열 해석 도구가 크랭크–니콜슨 같은 암시적 방법을 쓰는 이유다. 2차원이면 한계가 1/4로 더 엄격해진다.

SIMULATOR

1D 열확산: 핫스팟과 도펀트가 퍼지는 모습

—

초기 상태

양 끝 경계

r = αΔt/Δx² 프레임당 걸음 수 안정 조건 깨 보기 (r = 0.6) 이웃 평균 − 자기 (변화 방향) 보기

경과 시간—

최댓값—

총량 (∑, 초기 대비)—

걸음 수—

캔버스를 클릭하거나 끌면 그 위치에 열을 더합니다.

해볼 것: ① 핫스팟에서 재생 — 봉우리가 낮아지며 넓어지고, 폭은 $\sqrt{t}$로 자란다. 단열 경계에서는 총량이 보존되고, 방열판 경계에서는 끝으로 빠져나간다. ② '도펀트 확산'은 표면 농도를 고정한 채 실리콘 속으로 불순물이 들어가는 공정이다. 점선(해석해 $\mathrm{erfc}$)과 계산이 겹치는지 보라. ③ '안정 조건 깨 보기'를 켜면 칸마다 번갈아 튀는 톱니가 자라 계산이 폭발한다. · 모델: 80칸 명시적 유한 차분(FTCS). 열 모드는 길이 10 mm 실리콘 막대(α ≈ 88 mm²/s), 도펀트 모드는 깊이 0.5 µm, 확산계수 D ≈ 10⁻¹⁴ cm²/s(1000 °C 안팎의 붕소 정도 크기)로 시간을 환산했다. 표면에서 열이 공기로 빠지는 효과는 무시했다.

도펀트 확산 모드가 이 장의 마지막 의외성입니다. 반도체 공정에서 고온로에 웨이퍼를 넣어 붕소나 인을 실리콘 속으로 밀어 넣는 일은, 칩 위에서 열이 퍼지는 것과 완전히 같은 방정식을 따릅니다(픽의 제2법칙, $\partial C/\partial t = D\,\partial^2 C/\partial x^2$). 표면 농도를 고정하면 해는 $C = C_s\,\mathrm{erfc}\big(x/2\sqrt{Dt}\big)$이고, 확산 깊이는 $\sqrt{Dt}$로 자랍니다. D가 10⁻¹⁴ cm²/s 정도라면 한 시간에 $\sqrt{Dt}\approx 60$ nm — 접합 깊이를 수십 nm 단위로 제어할 수 있는 근거입니다. 실제 공정에서는 D가 온도에 지수적으로(아레니우스) 의존하고 농도·결함에 따라 달라지므로, TCAD가 이 식을 비선형으로 확장해 수치로 풉니다. 같은 식은 확률 쪽에서도 나타납니다: 무작위 걸음의 분포가 퍼지는 모습이 바로 이것이며, 그 폭이 $\sqrt t$로 자라는 이유는 7장의 중심극한정리가 설명합니다.

이 도구가 쓰이는 곳

미분방정식은 시리즈 거의 모든 책의 바닥에 깔려 있습니다. 아래 책들에서 "더 깊이" 링크를 따라왔다면, 이 장의 어느 식이 그 현상과 대응하는지 표로 확인해 보세요.

ANALOGBOOKRC·RLC 과도 응답
필터의 Q, 증폭기 정착 시간, 위상 여유 — 모두 특성근의 위치 문제. CARBOOK서스펜션과 승차감
감쇠비 ζ, 바퀴·차체의 2자유도 진동, 공진 회피. PACKAGINGBOOK열저항·열 시정수
다이–기판–방열판의 RC 열 회로망, 과도 열 임피던스. PROCESSBOOK도펀트 확산·열처리
픽의 법칙 $\partial C/\partial t = D\,\partial^2 C/\partial x^2$, 확산 깊이 $\sqrt{Dt}$. TCADBOOKPDE의 수치 해
확산–드리프트, 유한 차분·유한 요소, 안정성과 암시적 방법. DEVICEBOOK캐리어 연속 방정식
소수 캐리어 수명 τ, 생성–재결합 = 1차 감쇠, 확산 길이 $\sqrt{D\tau}$. BODYBOOK약동학·감염병
약물 반감기, 구획 모델, 체온 조절 피드백. PHONEBOOK스마트폰 열 관리
AP 쓰로틀링, 배터리 충전 곡선, 손떨림 보정 액추에이터.

핵심 정리

미분방정식은 "지금 값"이 아니라 "변화의 규칙"을 적는다. 출발점 하나가 주어지면 미래 전체가 정해진다.
"목표와의 차이에 비례해 변한다"는 모든 현상(RC 충전, 칩 온도, 방사성 붕괴, 약물 주입)은 $\dot y = -(y-y_\infty)/\tau$이고, 해는 $y_\infty + (y_0-y_\infty)e^{-t/\tau}$. τ에서 63%, 5τ에서 99.3%.
기울기장은 방정식을 그림으로 바꾼다. 평형점의 안정성, 해의 비교차, 정상 응답으로의 수렴을 공식 없이 읽을 수 있다.
오일러는 1차($O(h)$), RK4는 4차($O(h^4)$) 정확도다. 정확도와 별개로 $h\lambda$ 같은 안정성 한계가 있고, 강성 문제는 암시적 방법을 쓴다.
스프링–질량–감쇠기와 RLC는 같은 식 $\ddot x + 2\zeta\omega_0\dot x + \omega_0^2 x = 0$. 특성근이 실수 둘(과감쇠)·중근(임계)·켤레 복소수(진동)인지가 응답을 정한다. 실수부 = 감쇠, 허수부 = 진동.
공진에서 증폭률은 $Q = 1/2\zeta$, 대역폭은 $\omega_0/Q$. SIR 모델에서 정점은 $S = 1/R_0$, 집단면역 임계는 $1-1/R_0$. 열확산 $\partial_t T = \alpha\,\partial_x^2 T$의 2차 미분은 "이웃 평균과의 차이"이며, 명시적 계산은 $\alpha\Delta t/\Delta x^2 \le 1/2$을 지켜야 한다.

확인 퀴즈

열 시정수 τ = 60 s인 칩이 25 °C에서 정상 온도 85 °C를 향해 가열된다. 60초 뒤 온도는 대략?

t = τ에서 목표까지 거리의 63.2%를 간다. 25 + 0.632 × 60 ≈ 62.9 °C. 절반(55 °C)에 닿는 시각은 τ ln2 ≈ 42 s다.

$\dot y = -10\,y$를 오일러법, h = 0.25로 풀면 어떻게 되는가?

증폭률 1 − hλ = 1 − 2.5 = −1.5. |−1.5| > 1이므로 진동하며 발산한다. 안정하려면 hλ ≤ 2, 즉 h ≤ 0.2여야 한다. 정확도가 아니라 안정성의 문제다.

RLC 회로의 특성근이 $s = -2 \pm 5i$ (단위 1/s)일 때 옳은 것은?

켤레 복소근이므로 진동. 허수부 5 rad/s가 감쇠 진동수 $\omega_d$, 실수부 −2가 포락선 $e^{-2t}$(시정수 1/2 s). ω₀ = √(4+25) ≈ 5.39, ζ = 2/5.39 ≈ 0.37.

공진 주파수 1 MHz, Q = 50인 동조 회로의 −3 dB 대역폭은 대략?

대역폭 ≈ f₀/Q = 1 MHz / 50 = 20 kHz. Q가 높을수록 좁은 대역만 통과시켜 선택도가 좋아지지만, 응답이 정착하는 데 더 오래 걸린다.

$R_0 = 4$인 감염병에서 집단면역에 필요한 최소 면역 비율과, SIR 모델에서 감염자 수가 정점에 이르는 조건은?

$p_c = 1 - 1/R_0 = 0.75$. $\dot I = \gamma I(R_0 S - 1) = 0$이 되는 $S = 1/R_0 = 0.25$에서 감염자가 최대다.

명시적 유한 차분으로 열확산을 풀다가 정확도를 위해 격자 간격 Δx를 절반으로 줄였다. 안정성을 유지하려면 시간 걸음 Δt는?

안정 조건 $\alpha\Delta t/\Delta x^2 \le 1/2$에서 Δx²이 1/4이 되므로 Δt도 1/4. 칸 수는 두 배, 걸음 수는 네 배라 계산량은 8배가 된다. 그래서 암시적 방법이 쓰인다.

현상	y	y∞ (목표)	τ
RC 회로 충전 (카메라 플래시 축전기)	축전기 전압 \(V_C\)	전원 전압 \(V_s\)	\(RC\)
칩·커피 냉각/가열 (뉴턴 냉각)	온도 \(T\)	\(T_{amb}+PR_{th}\)	\(R_{th}C_{th}\)
방사성 붕괴	남은 핵 수 \(N\)	0	평균 수명 \(=t_{1/2}/\ln2\)
일정 속도 정맥 주입 (1구획 약동학)	혈중 농도 \(C\)	주입 속도 / 청소율	분포 용적 / 청소율
소수 캐리어 재결합 (반도체)	과잉 캐리어 \(\Delta n\)	0 (또는 생성률 × τ)	캐리어 수명