열방정식

수치해석 계산물리 마지막 수정: 2026-08-27 04:28:19

1. 개요[편집]

열방정식 Heat Equation
형태	u_t = α∇²u
유형	2계 선형 포물형 PDE
계수	열확산도 α = k/(ρc_p) [m²/s]
기본해	가우스 핵 (4παt)^−n/2 exp(−\|x\|²/4αt)
기원	푸리에, 열의 해석적 이론 (1822)
대표 성질	최대원리 · 즉시 평활화 · 무한 전파속도
수치적 특징	명시적 도식은 Δt ∝ Δx² 로 묶임

쌍곡형은 정보를 실어 나르고, 타원형은 균형을 잡고, 포물형은 잊는다.

열방정식(heat equation)은 어떤 스칼라장이 자기 자신의 곡률에 비례해 시간에 따라 변하는 2계 선형 포물형 편미분방정식 $u_t=\alpha\nabla^2 u$ 이다. 여기서 $\alpha$ 는 열확산도(thermal diffusivity)로 차원이 $\text{m}^2/\text{s}$ — 확산 현상의 시간·공간 스케일이 $L^2\sim\alpha t$ 로 묶인다는 사실이 이 상수 하나에 다 들어 있다.

이름은 열이지만 무대는 훨씬 넓다. 농도 확산(픽의 제2법칙), 중성자 확산, 소용돌이의 점성 감쇠, 브라운 운동의 전이확률 밀도, 영상 처리의 가우시안 블러, 금융의 블랙-숄즈(변수 변환하면 정확히 열방정식이다)까지 전부 같은 방정식이다. 그래서 열방정식은 “열 계산용 식”이라기보다 포물형 PDE의 원형(prototype) 으로 취급된다. 수치해석 교과서가 예제로 이것부터 꺼내는 이유도, 이 한 방정식에서 안정조건·강성·암시법의 필요성이 전부 유도되기 때문이다.

이 문서는 방정식 자체의 구조와 그 구조가 수치해석에 강제하는 것들을 다룬다. 열전도·대류·복사를 물리적으로 구분해 쓰는 공학적 관점은 열전달 해석, 대류항이 붙어 방향성이 생기는 순간의 이야기는 대류-확산 방정식, 시간 이산화 도식 자체의 세부는 크랭크-니콜슨법이 맡는다.

2. 어디서 오는가[편집]

유도는 두 줄이다. 첫째, 푸리에 법칙 — 열속은 온도 기울기의 반대 방향으로 흐른다.

\mathbf{q}=-k\nabla T

둘째, 에너지 보존 — 검사체적 안의 내부에너지 변화율은 들어오고 나간 열속의 차이다.

\rho c_p\frac{\partial T}{\partial t}=-\nabla\cdot\mathbf{q}+\dot q_v

둘을 합치고 $k$ 가 상수라 가정하면

\frac{\partial T}{\partial t}=\alpha\nabla^2 T+\frac{\dot q_v}{\rho c_p},\qquad \alpha=\frac{k}{\rho c_p}

주목할 점은 구성 법칙(푸리에 법칙)이 근사라는 것이다. 열속이 기울기에 즉각 비례한다는 가정은 곧 정보가 무한히 빨리 전달된다는 뜻이고, 여기서 유명한 역설이 나온다.

3. 기본해와 무한 전파속도[편집]

$n$ 차원 전 공간에서 초기조건이 원점의 델타 함수일 때의 해가 열핵(heat kernel)이다.

\Phi(\mathbf{x},t)=\frac{1}{(4\pi\alpha t)^{n/2}}\exp\!\left(-\frac{|\mathbf{x}|^2}{4\alpha t}\right),\qquad t>0

임의의 초기조건은 이것과의 합성곱으로 풀린다: $u(\mathbf{x},t)=\int\Phi(\mathbf{x}-\mathbf{y},t)\,u_0(\mathbf{y})\,d\mathbf{y}$ . 확산 폭이 $\sqrt{4\alpha t}$ 로 자란다는 것, 즉 거리는 시간의 제곱근이라는 확산의 핵심 스케일링이 지수부에 그대로 박혀 있다.

그런데 이 식은 $t>0$ 인 모든 $\mathbf{x}$ 에서 엄밀히 양수다. 원점에 성냥불을 켠 순간 $10^{9}$ 광년 밖 온도가 $10^{-10^{40}}$ K 만큼 오른다는 뜻이다. 값이 터무니없이 작아 실무에선 아무 문제가 없지만, 상대론적으로는 명백한 결함이다. 이를 고치려면 푸리에 법칙에 완화 시간 $\tau$ 를 넣어

\tau\frac{\partial\mathbf{q}}{\partial t}+\mathbf{q}=-k\nabla T

로 바꾸고(카타네오-베르노트), 그러면 방정식이 $\tau u_{tt}+u_t=\alpha\nabla^2u$ 라는 쌍곡형(전신 방정식)이 되어 전파속도가 $\sqrt{\alpha/\tau}$ 로 유한해진다. 금속의 $\tau$ 는 피코초 수준이라 일상 스케일에서는 무의미하지만, 초단펄스 레이저 가공이나 극저온 결정에서는 실제로 두 번째 소리(second sound)로 관측된다.¹

4. 분리변수와 푸리에급수[편집]

유한 구간 $[0,L]$ 에서 양끝 온도를 0으로 고정하면, $u=X(x)T(t)$ 로 놓는 순간 문제가 고유값 문제로 바뀐다.

u(x,t)=\sum_{n=1}^{\infty}b_n\sin\frac{n\pi x}{L}\;\exp\!\left(-\alpha\frac{n^2\pi^2}{L^2}t\right)

각 모드의 감쇠율이 파수의 제곱에 비례한다는 것이 이 식의 전부이자 열방정식 이해의 핵심이다. 파장이 절반이 되면 사라지는 속도가 네 배다. 그래서 초기조건에 아무리 험한 톱니가 있어도 고주파 성분은 순식간에 지워지고, 남는 것은 $n=1$ 모드뿐이다. 긴 시간 거동이 항상 “가장 부드러운 모드 하나”로 수렴한다는 이 성질이 뒤에서 강성(stiffness)의 원인으로 되돌아온다.

역사적으로 이 급수가 푸리에 변환이라는 분야 전체를 낳았다. 푸리에는 1807년 파리 과학원에 논문을 냈다가 라그랑주가 “임의의 함수를 삼각급수로 쓸 수 있다는 주장은 근거가 없다”고 반대하는 바람에 게재가 막혔고, 1811년 현상 공모에서 상은 받았지만 엄밀성 지적을 함께 받았으며, 단행본 Théorie analytique de la chaleur 는 1822년에야 나왔다.² 열 문제를 풀려다 해석학의 기초를 다시 짜게 만든 셈이다.

5. 경계조건과 무차원화[편집]

문제를 닫으려면 초기조건 하나와 경계조건이 필요하다. 세 가지 형태가 표준이고, 이름은 다르지만 결국 “경계에서 값을 주느냐, 기울기를 주느냐, 둘의 조합을 주느냐”의 차이다.

종류	형태	물리적 의미	이산화 시 주의
디리클레	$u=g$	벽 온도 고정	가장 쉬움. 무한 열용량 가정임을 잊지 말 것
노이만	$\partial u/\partial n=q$	열속 지정( $q=0$ 이면 단열)	유령 셀로 2차 정확도 유지, 전부 노이만이면 해가 상수만큼 부정
로빈	$-k\,\partial u/\partial n=h(u-u_\infty)$	대류 열전달	$h$ 가 크면 디리클레로, 0이면 노이만으로 퇴화

로빈 조건에서 나오는 무차원수가 비오 수 $Bi=hL/k$ 다. $Bi\ll0.1$ 이면 물체 내부 온도가 사실상 균일해져서 PDE를 아예 버리고 상미분방정식 하나(럼프드 모델)로 갈아탈 수 있다. 풀기 전에 $Bi$ 부터 계산하라는 조언이 열전달 수업에서 반복되는 이유다.

시간축을 무차원화하면 푸리에 수 $Fo=\alpha t/L^2$ 가 나온다. 확산이 특징 길이 $L$ 을 가로지르는 데 걸리는 시간이 $Fo\sim1$ 이라는 뜻이고, 실무적으로는 계산 종료 시각을 정하는 기준이다. $Fo=0.5$ 면 아직 과도기, $Fo>2$ 면 거의 정상 상태다. 눈치챘겠지만 FTCS의 안정 파라미터 $r=\alpha\Delta t/\Delta x^2$ 는 격자 셀 하나에 대한 푸리에 수 그 자체다. 격자 셀 안에서 확산이 한 셀을 넘어가면 명시적 도식이 터진다 — 안정조건을 외우지 말고 이렇게 기억하는 편이 낫다.

6. 최대원리와 평활화, 그리고 비가역성[편집]

포물형 방정식의 성격을 규정하는 정리가 최대원리다. 소스항이 없다면, 시공간 영역에서 $u$ 의 최댓값과 최솟값은 반드시 포물형 경계(초기 시각 + 옆면 경계)에서 달성된다. 내부에서 갑자기 새로운 최댓값이 생겨날 수 없다.

실무적 함의가 크다.

오버슈트가 없다. 초기 온도가 300~500 K 사이였다면 이후 어떤 시각에도 그 범위를 벗어나지 않는다. 계산 결과에 505 K 가 찍혔다면 물리가 아니라 도식을 의심해야 한다.
즉시 평활화. 초기조건이 불연속이어도 $t>0$ 이면 해는 $C^\infty$ 다. 계단 함수를 넣으면 다음 순간 오차함수(erf)가 되어 있다. 쌍곡형이 불연속을 그대로 실어 나르는 것과 정반대다.
비가역성. 시간을 거꾸로 돌리는 후방 열방정식은 비적정(ill-posed)이다. 고주파가 $e^{+\alpha k^2 t}$ 로 폭발하므로 입력의 미세한 잡음이 지수적으로 증폭된다. 열화상에서 내부 열원 분포를 되짚는 문제가 어려운 근본 이유이고, 역문제에서 정칙화가 반드시 등장하는 교과서적 사례다.

7. 수치해석 (1) — FTCS와 Δx² 장벽[편집]

1차원에서 시간은 전진 오일러, 공간은 중심차분으로 이산화한 것이 FTCS(Forward-Time Central-Space)다.

u_j^{n+1}=u_j^n+r\left(u_{j+1}^n-2u_j^n+u_{j-1}^n\right),\qquad r=\frac{\alpha\,\Delta t}{\Delta x^2}

폰 노이만 안정성 해석을 하면 증폭인자가

G=1-4r\sin^2\!\frac{k\Delta x}{2}

이고, $|G|\le1$ 은 곧 $r\le 1/2$ 다. 즉

\Delta t\le\frac{\Delta x^2}{2\alpha}

2차원 등간격에서는 $\Delta t\le\Delta x^2/(4\alpha)$ , 3차원에서는 $\Delta x^2/(6\alpha)$ 로 더 빡빡해진다. 이 조건의 악랄함은 지수에 있다. 격자를 두 배로 조밀하게 만들면 스텝 수가 네 배로 늘고, 스텝당 비용도 두 배(1D)라 총 비용이 여덟 배가 된다. CFL 조건( $\Delta t\propto\Delta x$ )이 선형인 것과 대비하면 포물형이 왜 명시적 도식과 궁합이 나쁜지 분명해진다.

$r=1/2$ 정확히 잡으면 재밌는 일이 생긴다. 위 식이 $u_j^{n+1}=(u_{j+1}^n+u_{j-1}^n)/2$ 가 되어 자기 자신이 사라진다 — 격자가 짝·홀로 완전히 분리되는 이른바 체스판 분리. 형식적으로는 안정하지만 실전에서는 $r\le0.4$ 근처로 여유를 두는 것이 관행이다.

8. 수치해석 (2) — 강성의 원천[편집]

FTCS 안정조건은 사실 행렬 고유값의 문제다. 디리클레 경계의 1차원 라플라시안을 이산화한 행렬의 고유값은

\lambda_m=-\frac{4\alpha}{\Delta x^2}\sin^2\!\frac{m\pi\Delta x}{2L},\qquad m=1,\dots,N

이고, 최소 크기는 $|\lambda_1|\approx\alpha\pi^2/L^2$ , 최대 크기는 $|\lambda_N|\approx 4\alpha/\Delta x^2$ 다. 둘의 비 — 강성비(stiffness ratio) — 는

\frac{|\lambda_{\max}|}{|\lambda_{\min}|}\approx\frac{4L^2}{\pi^2\Delta x^2}=O(N^2)

격자점 1000개짜리 1차원 문제만 해도 강성비가 $10^5$ 를 넘는다. 명시적 도식은 가장 빠른 모드에 시간 스텝을 맞춰야 하는데, 정작 그 모드는 시작하자마자 죽어서 해에 기여하지도 않는다. 관심 있는 것은 $\lambda_1$ 스케일의 느린 완화인데 $\lambda_N$ 이 시간 스텝의 목줄을 쥐고 있는 상황 — 이것이 강성 방정식의 정의 그 자체다. 열방정식이 강성 문제의 교과서 예제인 이유가 이 계산 한 줄이다.

9. 수치해석 (3) — 암시법과 ADI[편집]

해법은 암시법이다. 확산항을 새 시각에서 평가하면( $\theta$ -방법의 $\theta$ 를 1 또는 1/2 로) 안정조건이 사라진다.

도식	$\theta$	정확도	안정성	특징
FTCS (전진 오일러)	0	$O(\Delta t,\Delta x^2)$	$r\le1/2$	행렬 풀이 없음
후진 오일러	1	$O(\Delta t,\Delta x^2)$	무조건 안정	강한 소산, 매우 견고
크랭크-니콜슨법	1/2	$O(\Delta t^2,\Delta x^2)$	무조건 안정	고파수에서 $G\to-1$

1차원에서 암시법의 대가는 삼중대각 연립방정식 하나이고, 토머스 알고리즘으로 $O(N)$ 에 풀리므로 사실상 공짜다. 1차원에서 명시적 도식을 쓸 이유는 거의 없다.

문제는 2·3차원이다. 암시 이산화가 만드는 계수 행렬은 더 이상 삼중대각이 아니라 밴드 폭이 $N$ 인 희소행렬이 되고, 직접 분해하면 채움(fill-in) 때문에 비용이 폭발한다. 여기서 나온 고전적 우회로가 ADI(Alternating Direction Implicit, 피스먼-래치퍼드 1955)다. 한 반스텝은 $x$ 방향만 암시로, 다음 반스텝은 $y$ 방향만 암시로 처리하면, 매 반스텝이 독립적인 삼중대각 문제의 묶음으로 분해된다. 무조건 안정성과 2차 정확도를 유지하면서 비용은 $O(N)$ — 자세한 것은 ADI 법 참고. 현대적으로는 ADI 대신 다중격자법이나 크릴로프 부분공간법으로 암시 계를 푸는 쪽이 일반적이고, 주기 경계에서는 스펙트럴 방법으로 각 푸리에 모드를 정확히 지수 적분해 버리는 선택지도 있다.

10. 선형을 벗어나면[편집]

실제 해석에서 만나는 것은 순수한 상수계수 열방정식이 아니다. 어디서 선형성이 깨지는지, 그리고 그때 무엇이 살아남는지가 실무 감각이다.

변수 계수. 물성이 온도에 의존하면 $\rho c_p\,T_t=\nabla\cdot(k(T)\nabla T)$ 가 된다. 이때 $k$ 를 셀 중심값의 산술평균으로 면에 보간하면 열속 보존이 깨진다. 직렬 저항 논리에 맞는 것은 조화평균이고, 물성 대비가 큰 다층 구조에서는 이 한 줄이 결과를 몇 십 %씩 바꾼다.
상변화. 융해·응고가 끼면 잠열 때문에 엔탈피가 온도의 다가 함수가 되고, 이동 경계를 가진 스테판 문제가 된다. 온도를 미지수로 두는 대신 엔탈피를 미지수로 두는 엔탈피법이 열처리 시뮬레이션이나 용접·적층 해석의 표준 처방이다.
비선형 확산. $u_t=\nabla\cdot(u^m\nabla u)$ 같은 다공질 매체 방정식은 확산계수가 $u=0$ 에서 사라져서, 놀랍게도 해가 유한 전파속도를 갖는다. 선형 열방정식의 무한 전파속도 문제가 비선형성으로 자연스럽게 해소되는 사례다.
복사가 끼면. 표면 복사는 $T^4$ 항을 경계조건에 넣어 문제를 비선형으로 만든다. 고온 해석에서 뉴턴법 반복이 붙는 지점이 대개 여기다. 자세한 것은 열전달 해석 참고.

검증은 오히려 쉬운 편이다. 정확해가 널려 있기 때문이다. 반무한 고체의 계단 응답이 오차함수, 유한 판의 과도 응답이 푸리에급수, 점열원 응답이 가우스 핵 — 어느 쪽이든 닫힌 형태와 직접 비교할 수 있다. 정확해가 없는 형상에서는 원하는 해를 먼저 정해 놓고 그것이 만족하도록 소스항을 역산하는 제조해법(MMS)이 표준이며, 확산 문제는 소스항이 그냥 라플라시안이라 MMS를 붙이기 가장 쉬운 방정식 축에 든다. 여기서 관측 수렴 차수가 2가 안 나오면 대개 경계조건 이산화가 1차로 떨어져 있다.³

11. 여담: 확산은 어디에나 있다[편집]

점성 소실 극한. 버거스 방정식의 점성항이 하는 일이 정확히 이 평활화이고, 콜-호프 변환은 비선형 버거스를 정확히 열방정식으로 바꿔 버린다. 비선형 보존법칙의 엔트로피 조건이 “점성 소실 극한”으로 정의되는 것도 같은 맥락 — 물리적 해를 고르는 심판이 결국 열방정식이다.
수치 확산. 대류항을 1차 상류차분으로 이산화하면 수정 방정식에 $\frac{U\Delta x}{2}u_{xx}$ 라는 항이 튀어나온다. 즉 의도하지 않은 열방정식이 코드에 끼어 있는 것이고, 이것이 대류-확산 방정식에서 말하는 거짓 확산이다. 열방정식을 아는 사람만 자기 코드가 얼마나 뭉개고 있는지 정량적으로 안다.
영상과 기하. 이미지에 가우시안 블러를 씌우는 것은 열방정식을 짧게 푸는 것과 같다. 여기서 확산계수를 기울기에 따라 바꾸면 경계는 살리고 잡음만 지우는 비등방 확산(페로나-말릭)이 되고, 곡면을 자기 평균곡률로 흘리면 평균곡률 흐름이 된다. 레벨셋 방법과 그래프 라플라시안 기반 다양체 학습이 같은 언어를 쓰는 이유다.
무작위 행보와의 쌍대성. 격자 위 무작위 행보의 확률분포가 만족하는 방정식이 정확히 이산 열방정식이다. FTCS 계수 $r\le1/2$ 는 “왼쪽·오른쪽·제자리 확률이 전부 음이 아니어야 한다”는 요구와 같고, 안정조건이 곧 확률 해석의 정합성이 된다. 안정조건을 어긴 도식이 음의 온도를 뱉는 이유도 여기서 설명된다.⁴

12. 관련 문서[편집]

“두 번째 소리”는 초유체 헬륨이나 고순도 결정에서 열이 파동처럼 전파되는 현상이다. 열이 파동이라는 말은 처음 들으면 헛소리 같지만, 포논이 서로 충돌하며 유체처럼 굴 때 실제로 그렇게 된다. 다만 이걸 보려면 온도를 몇 켈빈까지 내려야 해서, 상온 CFD 엔지니어의 현생과는 관계가 없다. ↩
라그랑주의 반대는 심술이 아니라 정당한 지적이었다. 불연속 함수를 연속함수의 무한합으로 쓴다는 발상은 당시 해석학의 기초를 흔드는 것이었고, 이게 제대로 정리되기까지 디리클레·리만·칸토어를 거쳐 한 세기가 걸렸다. 실용주의자가 먼저 답을 맞히고 이론가가 뒤늦게 왜 맞았는지 밝히는 패턴은 수치해석 바닥에서 지금도 매일 반복된다. ↩
관측 차수가 2.0 대신 1.0 근처로 나오면 십중팔구 노이만 경계에서 유령 셀을 안 쓰고 한쪽차분으로 때웠기 때문이다. 내부는 2차인데 경계 한 줄이 1차면 전체가 1차가 된다 — 수치해석의 가장 흔한 자해 중 하나. ↩
그래서 확산 문제에서 “안정성”과 “유계성”과 “최대원리”는 사실상 같은 이야기의 세 얼굴이다. 도식이 음수 계수를 갖는 순간 셋이 동시에 무너진다. 계산 결과에 −3 K 가 찍혀 있다면 방금 절대영도를 뚫은 게 아니라 $r$ 이 0.5를 넘었을 뿐이다. ↩