3단계: 수치 해석, 혼돈의 끝자락, 그리고 무작위성의 방정식

서론: 지금까지 온 길을 먼저 돌아보자

1단계에서 너는 **상미분방정식(ODE)**을 해석적으로 풀었다. $y^{''} + p (x) y^{'} + q (x) y = 0$ 같은 형태, 라플라스 변환으로 미분 연산을 대수 연산으로 바꾸는 기술, 그리고 Frobenius 급수 해법까지. 2단계에서는 공간과 시간이 동시에 얽히는 **편미분방정식(PDE)**으로 무대를 넓혔다. 열전도 방정식 $\frac{\partial u}{\partial t} = α \frac{\partial ^{2} u}{\partial x ^{2}}$ , 변수 분리법, 푸리에 급수로 경계 조건에 맞는 해를 찾는 법. 아름답고 정교했다.

그런데 지금 솔직하게 질문하자. 그 해석적 해가 항상 존재하는가? 잠깐 스스로 생각해 봐라. 만약 열전도 계수가 온도에 따라 달라진다면? $\frac{\partial u}{\partial t} = \frac{\partial}{\partial x} [α (u) \frac{\partial u}{\partial x}]$ . 변수 분리법은 이 순간 무너진다. 도시 전체의 공기 흐름, 비행기 날개처럼 복잡한 기하학, CPU 칩의 비균등 열 분포 — 자연계의 95%는 해석적으로 풀리지 않는다. 그래서 수학자와 공학자들은 전혀 다른 전략을 개발했다. 컴퓨터에게 시키는 것이다. 이것이 3단계의 출발점이다.

PART I. 유한차분법 (FDM, Finite Difference Method)

연속을 이산으로 쪼개다

일곱 살짜리 아이에게 미분을 설명한다면 이렇게 할 수 있다. "산길을 걸을 때 한 발짝 앞이 얼마나 높아지는지 재는 거야. 발짝을 무한히 작게 줄이면 그게 순간 기울기야." 수치 해법은 이것을 뒤집는다. 발짝을 유한하게 유지하는 것이다.

수학적으로 $\frac{d u}{d x} = lim_{Δ x \to 0} \frac{u ( x + Δ x ) - u ( x )}{Δ x}$ 인데, 극한을 취하지 않고 $Δ x$ 를 작지만 유한한 값으로 유지하면 **차분(difference)**이 된다. 세 가지 방식이 있다:

$\frac{d u}{d x} \approx \frac{u ( x + Δ x ) - u ( x )}{Δ x} (전진 차분, Forward Difference)$ $\frac{d u}{d x} \approx \frac{u ( x ) - u ( x - Δ x )}{Δ x} (후진 차분, Backward Difference)$ $\frac{d u}{d x} \approx \frac{u ( x + Δ x ) - u ( x - Δ x )}{2Δ x} (중앙 차분, Central Difference)$

[노트 기록] 위 세 공식을 써라. 그리고 테일러 전개 $u (x + Δ x) = u (x) + u^{'} Δ x + \frac{u ^{''}}{2} (Δ x)^{2} + \dots$ 를 이용해서 전진 차분의 오차가 $O (Δ x)$ 이고, 중앙 차분의 오차가 $O ((Δ x)^{2})$ 임을 직접 유도해 봐라. 왜 중앙 차분이 더 정확한지 스스로 납득이 가는가?

2계 도함수는 중앙 차분을 한 번 더 적용하면 나온다:

$\frac{d ^{2} u}{d x ^{2}} \approx \frac{u ( x + Δ x ) - 2 u ( x ) + u ( x - Δ x )}{( Δ x ) ^{2}}$

이 공식의 구조를 잘 봐라. $u (x)$ 를 기준으로 양쪽 이웃의 값을 "잡아당기는" 것처럼 보이지 않는가? 2단계에서 열전도 방정식을 물리적으로 해석할 때 "이웃과의 온도 차이가 열의 흐름을 결정한다"고 했다. 그 직관과 정확히 일치한다.

열방정식의 이산화와 명시적 스킴

2단계에서 해석적으로 풀었던 1차원 열방정식을 다시 꺼내자. 공간 $x$ 를 $N + 1$ 개 균등 격자 $x_{i} = i \cdot Δ x$ 로, 시간 $t$ 를 $t^{n} = n \cdot Δ t$ 로 이산화하고, $u (x_{i}, t^{n})$ 을 $u_{i}^{n}$ 으로 표기한다. 시간 도함수를 전진 차분으로, 공간 2계 도함수를 중앙 차분으로 근사하면:

$\frac{u _{i}^{n + 1} - u _{i}^{n}}{Δ t} = α \cdot \frac{u _{i + 1}^{n} - 2 u _{i}^{n} + u _{i - 1}^{n}}{( Δ x ) ^{2}}$

정리하면 FTCS 명시적 스킴(Forward Time, Centered Space):

$u_{i}^{n + 1} = u_{i}^{n} + r (u_{i + 1}^{n} - 2 u_{i}^{n} + u_{i - 1}^{n}), r = \frac{α Δ t}{( Δ x ) ^{2}}$

이 식의 의미가 보이는가? 다음 시각의 온도 $u_{i}^{n + 1}$ 를 현재 시각의 값들만으로 명시적으로 계산할 수 있다. 행렬 방정식이 필요 없다. 코드로 직접 반복 계산하면 된다.

폭발하는 계산: 안정성 조건

[노트 기록] $r$ 을 **푸리에 수(Fourier number)**라고 부른다. 물리적으로는 시간 보폭 $Δ t$ 동안 열이 공간 보폭 $Δ x$ 를 얼마나 "통과"하는지를 나타내는 무차원수다.

이제 핵심적인 질문이다. $Δ t$ 를 키우면 $r$ 이 커지고 계산이 빨라진다. 그런데 실제로 그렇게 하면 어떤 일이 벌어질까? $r = 1$ 이라고 해보자. $u_{i}^{n + 1} = - u_{i}^{n} + u_{i + 1}^{n} + u_{i - 1}^{n}$ 이 된다. 초기에 모든 격자값이 0인데 $u_{0}^{0} = 1$ 인 작은 펄스 하나가 있다고 상상해 봐라. 한 스텝, 두 스텝 후 어떻게 될지 손으로 직접 계산해 봐라. 뭔가 이상하지 않은가?

폰 노이만 안정성 분석(Von Neumann stability analysis)은 수치 해의 오차를 푸리에 모드 $e^{i ξ x}$ 로 분해하고 각 모드의 증폭 인자(amplification factor)를 분석하는 기법이다. 이를 통해 FTCS 스킴은 반드시 $r \leq \frac{1}{2}$ 이어야 안정적임을 증명할 수 있다. 이것이 열방정식에서의 **CFL 조건(Courant-Friedrichs-Lewy condition)**이다 (Strikwerda, J.C., Finite Difference Schemes and Partial Differential Equations, 2004). 공간 해상도를 두 배 높이면( $Δ x$ 를 반으로 줄이면), $Δ t$ 는 네 배 줄여야 한다. 계산 비용은 2D 문제에서 8배, 3D에서 16배로 폭발한다. 이것이 기상 예측이나 유체 시뮬레이션에 슈퍼컴퓨터가 필요한 실질적인 이유다.

이 제약을 우회하는 방법이 **암시적 스킴(Implicit Scheme)**이다. 공간 차분을 미래 시각 $n + 1$ 에서 평가하면, 미지수 $u_{i}^{n + 1}$ 에 대한 연립방정식이 생기지만 어떤 $Δ t$ 에서도 안정적이다. Crank-Nicolson 스킴은 현재와 미래의 평균을 사용해 2차 정확도와 무조건 안정성을 동시에 달성한다. "공짜 점심은 없다" — 연립방정식 풀이라는 계산 비용을 지불하고 안정성을 얻는다.

PART II. 유한요소법 기초 (FEM, Finite Element Method)

FDM과 FEM: 무엇이 근본적으로 다른가?

FDM은 격자 위의 **점(point)**에서 해를 근사한다. FEM은 영역을 작은 **요소(element)**로 쪼개고 각 요소 위에서 해를 함수의 조합으로 근사한다. 비행기 날개, 인체 뼈, 이형 단면을 가진 구조물처럼 불규칙한 기하학에서 FEM이 압도적으로 유리한 이유가 여기에 있다.

약형식의 탄생

FEM의 핵심은 **약형식(Weak Formulation)**이다. Poisson 방정식 $- \nabla^{2} u = f$ 를 모든 점에서 정확히 만족해야 하는 형태를 **강형식(Strong Form)**이라 한다. 약형식에서는 임의의 시험함수(test function) $v$ 를 곱하고 영역 전체에 적분한다:

$- \int_{Ω} (\nabla^{2} u) v d Ω = \int_{Ω} f v d Ω$

좌변에 Green의 정리(부분 적분의 고차원 버전)를 적용하면:

$\int_{Ω} \nabla u \cdot \nabla v d Ω = \int_{Ω} f v d Ω$

[노트 기록] 약형식의 공식을 써라. 강형식에는 2계 도함수가 있었는데, 약형식에는 1계 도함수만 남았다. 이것이 왜 중요한가? 해 $u$ 가 2번 미분 가능하지 않더라도 — 예를 들어 재료 경계에서 불연속적인 도함수를 가지더라도 — 약형식의 해는 존재할 수 있다. 즉, 부드러움 요구 조건이 낮아진 것이다.

이후 해를 기저함수(basis function, 1차원에서는 hat function)의 선형 조합 $u \approx \sum_{j} u_{j} ϕ_{j} (x)$ 로 근사하면 약형식은 **강성 행렬(stiffness matrix) 방정식 $Ku = f$ **로 귀결된다. FEM은 미분방정식 문제를 선형대수 문제로 바꾸는 변환이다. 컴퓨터는 미분을 직접 다루지 않아도 된다.

PART III. 비선형 동역학과 안정성 분석

선형 vs. 비선형: 왜 차이가 극적인가?

1단계에서 선형 ODE의 해는 지수함수와 삼각함수의 조합이었다. 단조롭고 예측 가능했다. 비선형 시스템은 다르다. 매개변수 하나가 바뀌는 것만으로 조용한 평형이 주기적 진동으로, 진동이 혼돈으로 변할 수 있다. 물을 천천히 데울 때 처음엔 고요히 가열되다가 어느 순간 대류 소용돌이가 생기고, 그다음엔 완전히 불규칙한 끓음이 되는 것처럼.

평형점과 야코비안 선형화

비선형 2차원 자율 시스템을 고려하자:

$\overset{x}{˙} = f (x, y), \overset{y}{˙} = g (x, y)$

평형점(equilibrium point) $(\overset{x}{ˉ}, \overset{y}{ˉ})$ 는 $f (\overset{x}{ˉ}, \overset{y}{ˉ}) = 0$ 이고 $g (\overset{x}{ˉ}, \overset{y}{ˉ}) = 0$ 을 동시에 만족하는 점으로, 시스템이 도달하면 영원히 머무른다. 평형점 근방에서 $x = \overset{x}{ˉ} + ξ$ , $y = \overset{y}{ˉ} + η$ 로 미소 교란을 도입하면 **선형화(linearization)**가 가능하다:

$(\dot{ξ} \overset{η}{˙}) \approx J (ξ η), J = (\partial f / \partial x \partial g / \partial x \partial f / \partial y \partial g / \partial y)_{(\overset{x}{ˉ}, \overset{y}{ˉ})}$

$J$ 가 **야코비안 행렬(Jacobian matrix)**이다. 비선형 세계를 평형점 근방에서 1단계에서 배운 선형 시스템으로 근사하는 것이다.

[노트 기록] 야코비안의 고유값 $λ_{1}, λ_{2}$ 와 평형점 성질의 관계를 표로 만들어라. 핵심만 적으면: (a) 둘 다 실수 음수 → 안정 절점(stable node), (b) 둘 다 실수 양수 → 불안정 절점, (c) 하나 양수 하나 음수 → 안장점(saddle point), (d) 복소수이고 실수부 음수 → 안정 나선(stable spiral), (e) 복소수이고 실수부 양수 → 불안정 나선, (f) 순허수 → 중심(center) — 이 경우 선형화만으로는 실제 안정성을 판단할 수 없다.

리아푸노프 직접법

선형화로 판단이 불가능한 경우 — 특히 고유값이 순허수일 때 — 에 알렉산드르 리아푸노프가 1892년 개발한 **리아푸노프 직접법(Lyapunov's Direct Method)**이 있다 (The General Problem of the Stability of Motion, 1892). 아이디어는 물리적 에너지와 유사하다. "에너지 같은 함수" $V (x)$ 를 찾아서 두 조건을 확인한다. 첫째 $V (x) > 0$ (원점에서만 0, 나머지에서는 양수), 둘째 $\dot{V} = \frac{d V}{d t} \leq 0$ (시간에 따라 감소하거나 유지). 이 두 조건을 만족하는 리아푸노프 함수가 존재하면 평형점은 안정적이다. $\dot{V} < 0$ 이면 **점근 안정(asymptotically stable)**이다. 마찰이 있는 진자가 에너지를 잃고 평형으로 수렴하는 것을 생각해 봐라 — 그 에너지 함수가 리아푸노프 함수다.

분기 이론 (Bifurcation Theory)

**분기(bifurcation)**는 매개변수 $μ$ 가 임계값 $μ_{c}$ 를 넘을 때 시스템의 위상학적(topological) 행동 자체가 질적으로 변하는 현상이다. 단순히 값이 커지거나 작아지는 것이 아니라, 평형점의 개수가 바뀌거나 안정성이 뒤집어진다.

가장 단순한 예로 안장-절점 분기(Saddle-Node Bifurcation): $\overset{x}{˙} = μ - x^{2}$ . $μ < 0$ 이면 평형점이 없고, $μ = 0$ 이면 $x = 0$ 하나, $μ > 0$ 이면 $x = \pm μ$ 두 개가 생긴다. 분기점 $μ_{c} = 0$ 에서 "없음 → 둘"로 점프한다.

더 흥미로운 것이 **호프 분기(Hopf Bifurcation)**다. 안정적인 나선형 평형점이 $μ = μ_{c}$ 에서 불안정해지면서 자발적으로 극한 순환(limit cycle), 즉 주기적 진동이 탄생한다. 심장 박동, 신경 세포의 발화(firing), Belousov-Zhabotinsky 화학 진동자 모두 호프 분기의 실제 사례다. 매개변수 하나의 변화가 "정지 → 진동"이라는 근본적 행동 변화를 만든다는 것이 얼마나 놀라운가.

**로렌츠 시스템(Lorenz System)**은 이 모든 것의 정점이다:

$\overset{x}{˙} = σ (y - x), \overset{y}{˙} = x (ρ - z) - y, \overset{z}{˙} = x y - β z$

에드워드 로렌츠는 1963년 대기 대류 모델인 이 방정식이 초기 조건에 극도로 민감함을 발견했다. 초기값 0.506127과 0.506128에서 출발한 두 궤적은 단기적으로는 거의 같지만 장기적으로는 완전히 다른 경로를 간다. **나비 효과(Butterfly Effect)**의 수학적 실체다. $ρ = 28, σ = 10, β = 8/3$ 에서 로렌츠 끌개는 카오스적이다 (Strogatz, S.H., Nonlinear Dynamics and Chaos, 2nd ed., 2015). 이 시스템은 결정론적이지만 — 초기값이 정해지면 미래가 완전히 결정된다 — 실용적으로는 예측 불가능하다. 이것이 카오스의 역설이다.

[노트 기록] 로렌츠 시스템의 평형점을 직접 구하는 절차를 써라. $\overset{x}{˙} = \overset{y}{˙} = \overset{z}{˙} = 0$ 을 풀면 $ρ > 1$ 일 때 세 평형점이 나온다. 각 평형점의 좌표를 구하고 (계산은 프로젝트에서), 야코비안을 구성하는 절차를 기술하라.

PART IV. 확률 미분방정식 (SDE)과 금융 모델링

무작위성을 방정식에 집어넣다

지금까지 모든 방정식은 **결정론적(deterministic)**이었다. 초기 조건을 알면 미래를 이론적으로 완벽히 예측할 수 있다. 그런데 세상에는 본질적으로 불확실성이 내재된 현상들이 있다. 주가, 열 잡음, 이온 채널의 개폐, 양자 입자의 운동.

1827년 식물학자 로버트 브라운(Robert Brown)은 물 위 꽃가루 입자의 불규칙한 운동을 관찰했다. 1905년 아인슈타인이 이것을 수백만 개 물 분자의 무작위 충돌 결과임을 보였고, 1923년 노버트 위너(Norbert Wiener)가 위너 과정(Wiener Process) $W (t)$ 를 수학적으로 엄밀히 정의했다. 핵심 성질은 네 가지다. $W (0) = 0$ 이고, 겹치지 않는 구간의 변화량은 서로 독립적이며, $W (t) - W (s) \sim N (0, t - s)$ (정규 분포이고 분산이 시간 차이에 비례)이고, 경로는 연속이지만 어디서도 미분 불가능하다. 마지막 성질을 깊이 생각해 봐라. 위너 과정의 경로는 무한히 들쭉날쭉해서 기울기가 존재하지 않는다. 그렇다면 $d W$ 를 어떻게 다루어야 하는가?

이토 미적분의 핵심

SDE의 일반 형태는 **이토 확률 미분방정식(Itô SDE)**이다:

$d X = μ (X, t) d t + σ (X, t) d W$

$μ$ 는 드리프트(drift) 계수 — 평균적인 방향, $σ$ 는 확산(diffusion) 계수 — 무작위성의 세기. 여기서 일반 미적분과 근본적으로 다른 점이 있다. 일반 미적분에서 $(d x)^{2} \approx 0$ 이지만, 확률 과정에서 (dW)² = dt가 된다 (이토 등거리성, Itô Isometry). 이 작은 차이가 모든 것을 바꾼다.

이로부터 **이토의 보조정리(Itô's Lemma)**가 나온다. $V = V (X, t)$ 라면:

$d V = 일반 연쇄법칙 + 이토 수정항 (\frac{\partial V}{\partial t} + μ \frac{\partial V}{\partial X} + \frac{σ ^{2}}{2} \frac{\partial ^{2} V}{\partial X ^{2}}) d t + σ \frac{\partial V}{\partial X} d W$

일반 연쇄법칙과 비교해 봐라. 추가 항 $\frac{σ ^{2}}{2} \frac{\partial ^{2} V}{\partial X ^{2}}$ 이 새로 나타났다. 이것이 무작위성이 만드는 **이토 수정항(Itô correction term)**이다. 이 항은 2계 도함수를 포함하기 때문에, 볼록한(convex) 함수에는 비선형적으로 작용한다 — 이것이 금융에서 **감마 효과(gamma effect)**로 알려진 현상의 수학적 뿌리다.

[노트 기록] 이토의 보조정리를 써라. 일반 연쇄법칙 $d V = \frac{\partial V}{\partial t} d t + \frac{\partial V}{\partial X} d X$ 와 어떻게 다른지, 추가 항이 어디서 비롯되는지 명확히 표시하라.

블랙-숄즈: 열방정식의 변장

1973년 피셔 블랙(Fischer Black)과 마이런 숄즈(Myron Scholes)는 유럽형 콜 옵션 가격을 결정하는 PDE를 도출했다. 주가 $S$ 가 기하 브라운 운동(Geometric Brownian Motion) $d S = μ S d t + σ S d W$ 를 따른다고 가정하고, 옵션 $V (S, t)$ 에 이토의 보조정리를 적용한 뒤, 델타 헤징(delta hedging)으로 무위험 포트폴리오를 구성하면 $μ$ 가 소거되면서:

$\frac{\partial V}{\partial t} + \frac{1}{2} σ^{2} S^{2} \frac{\partial ^{2} V}{\partial S ^{2}} + r S \frac{\partial V}{\partial S} - r V = 0$

이것이 블랙-숄즈 PDE다. 놀랍게도, 변수 변환 $S = K e^{x}$ , $V = K e^{α x + β τ} u (x, τ)$ 를 적용하면 이 방정식이 정확히 2단계에서 배운 열전도 방정식으로 변환된다. 금융과 열역학이 동일한 수학적 구조를 공유한다는 것이다. 우주의 심오한 농담처럼 느껴지지 않는가?

수치적으로는 FDM으로 직접 풀거나, 오일러-마루야마(Euler-Maruyama) 방법으로 SDE를 이산화하는 몬테카를로 시뮬레이션으로 풀 수 있다:

$S_{n + 1} = S_{n} + r S_{n} Δ t + σ S_{n} Δ t \cdot Z_{n}, Z_{n} \sim N (0, 1)$

수천 개의 주가 경로를 시뮬레이션하고 만기 페이오프를 평균하면 현재 옵션 가격이 나온다. 숄즈와 머튼은 이 연구로 1997년 노벨 경제학상을 받았다. 그런데 이 모델을 실제로 사용한 헤지펀드 LTCM(Long-Term Capital Management)은 1998년 러시아 채무 불이행으로 40억 달러 손실을 내고 붕괴했다 (Hull, J.C., Options, Futures, and Other Derivatives, 11th ed., 2021). 모델의 가정(연속 헤징, 정규 분포 수익률)이 실제 시장 위기에서 깨진 것이다. 수학 모델의 힘과 한계를 동시에 보여주는 역사다.

프로젝트 (각자 40분, 정답 없음)

프로젝트 1: 1차원 열방정식 FDM 솔버

설정. 막대 길이 $L = 1$ , 열확산 계수 $α = 0.01$ , 초기 온도 $u (x, 0) = sin (π x)$ , 경계 조건 $u (0, t) = u (1, t) = 0$ . 2단계에서 구한 해석적 해는 $u (x, t) = e^{- α π^{2} t} sin (π x)$ 이다.

문제 1-1. FTCS 명시적 스킴으로 $t = 0$ 에서 $t = 2$ 까지 시뮬레이션하는 알고리즘을 Python 의사 코드 또는 실제 코드로 작성하라. $N = 20$ (공간 격자), $Δ t = 0.1$ 을 사용하라. 먼저 CFL 수 $r = \frac{α Δ t}{( Δ x ) ^{2}}$ 를 계산하고, 이 설정이 안정 조건 $r \leq 1/2$ 을 만족하는지 확인하라.

문제 1-2. $r > 1/2$ 가 되도록 $Δ t$ 를 증가시키면 어떤 일이 일어나는지 예측하고 실제로 확인하라. 불안정이 발생하는 $t$ 를 수치로 특정하라 (예: 수치 해의 최솟값이 $- 10$ 이하가 되는 시각).

문제 1-3. 같은 문제를 Crank-Nicolson 스킴으로 구현하라. $t = 0.5$ 시점에서 FTCS, Crank-Nicolson, 해석적 해의 $L^{2}$ 오차 $∥ u_{numeric} - u_{exact} ∥_{2}$ 를 계산하고 비교하라. 정확도 차이가 이론적 예측(1차 vs. 2차 정확도)과 일치하는가?

문제 1-4. 열방정식에 원천항을 추가한 $\frac{\partial u}{\partial t} = α \frac{\partial ^{2} u}{\partial x ^{2}} + sin (2 π x) cos (t)$ 를 고려하라. FTCS 스킴을 어떻게 수정해야 하는가? 이 경우 해석적 특수해를 가정할 수 있는가? (힌트: $u_{p} = A (x) sin (t) + B (x) cos (t)$ 형태를 시도해 봐라.)

프로젝트 2: 비선형 안정성 분석 — 포식자-피식자 모델과 분기

설정. 로트카-볼테라 방정식(Lotka-Volterra):

$\overset{x}{˙} = α x - β x y, \overset{y}{˙} = δ x y - γ y$

$α = 1.0$ , $β = 0.1$ , $δ = 0.075$ , $γ = 1.5$ . $x$ : 토끼(피식자), $y$ : 여우(포식자).

문제 2-1. 시스템의 모든 평형점을 구하라. $(0, 0)$ 을 포함해 몇 개인가? 각 평형점에서 야코비안을 계산하고 고유값을 구해 유형을 분류하라. Part III에서 정리한 표를 활용하라.

문제 2-2. 리아푸노프 함수 후보 $V (x, y) = δ x - γ ln x + β y - α ln y$ 에 대해 $\dot{V} = \frac{d V}{d t}$ 를 계산하라. 결과가 예상과 다르다면 그 물리적 의미가 무엇인지 해석하라. (힌트: 로트카-볼테라 시스템은 보존 시스템이다.)

문제 2-3. 사냥 제어를 추가한다: $\overset{x}{˙} = α x - β x y - h x$ , $\overset{y}{˙} = δ x y - γ y - h y$ , $h \geq 0$ . $h$ 가 증가함에 따라 비자명 평형점의 위치가 어떻게 이동하는가? 분기가 발생하는 임계값 $h_{c}$ 를 찾고 분기의 종류를 판별하라.

문제 2-4. 반 데르 폴 진동자 $\overset{x}{¨} - μ (1 - x^{2}) \overset{x}{˙} + x = 0$ 을 $x_{1} = x$ , $x_{2} = \overset{x}{˙}$ 로 2차원 시스템으로 변환하라. $μ = 0$ 일 때와 $μ > 0$ 일 때 원점의 안정성이 어떻게 달라지는가? $μ$ 를 0에서 증가시킬 때 어떤 유형의 분기가 발생하는가? (힌트: 야코비안 고유값의 실수부가 $μ$ 에 따라 어떻게 변하는지 추적하라.)

프로젝트 3: 블랙-숄즈 구현 — FDM vs. 몬테카를로

설정. 유럽형 콜 옵션. $S_{0} = 100$ , $K = 100$ , $r = 0.05$ , $σ = 0.2$ , $T = 1$ . 해석적 해:

$C = S_{0} N (d_{1}) - K e^{- r T} N (d_{2}), d_{1} = \frac{ln ( S _{0} / K ) + ( r + σ ^{2} /2 ) T}{σ T}, d_{2} = d_{1} - σ T$

문제 3-1. 오일러-마루야마 방법으로 기하 브라운 운동을 이산화하고 $N = 1000$ 개의 주가 경로를 시뮬레이션하라. 각 경로에서 만기 페이오프 $max (S_{T} - K, 0)$ 를 계산하고 $e^{- r T}$ 를 곱해 현재 가치로 환산한 뒤 평균하여 몬테카를로 가격을 추정하라. 해석적 해와 비교하라.

문제 3-2. 몬테카를로 표준 오차는 $\frac{s}{N}$ 이다. $N = 100, 1000, 10000, 100000$ 으로 변화시키며 오차를 측정하라. 오차 감소율이 $O (N^{- 1/2})$ 와 일치하는가? 정확도를 10배 높이려면 계산 비용이 몇 배 필요한가?

문제 3-3. FDM으로 블랙-숄즈 PDE를 직접 풀어라. 공간 격자 $S \in [0, 2 K]$ , $N_{S} = 100$ , $N_{t} = 200$ . 경계 조건: $V (0, t) = 0$ , $V (2 K, t) = 2 K - K e^{- r (T - t)}$ . 만기 조건: $V (S, T) = max (S - K, 0)$ . 시간을 만기에서 현재로 역방향으로 진행해야 한다는 점에 주의하라. 왜 역방향인가? 수치해와 해석해의 $L^{2}$ 오차를 계산하라.

문제 3-4. (심화) 블랙-숄즈의 핵심 가정은 변동성 $σ$ 가 상수라는 것이다. 실제 시장에서는 행사 가격별 내재 변동성(implied volatility)이 "스마일(smile)" 곡선을 그린다 — 이 현상이 왜 블랙-숄즈 가정의 위반을 나타내는가? 변동성을 $σ (S, t)$ 로 일반화한 **국소 변동성 모델(Local Volatility Model)**에서 블랙-숄즈 PDE는 어떻게 바뀌는가? FDM 코드에서 어느 부분을 수정해야 하는가?

마치며

오늘 우리는 세 개의 세계를 연속적으로 연결했다. FDM은 2단계의 해석적 해를 근사하는 이산화이고, 블랙-숄즈는 열방정식의 변장이며, 리아푸노프 함수는 1단계 에너지 보존 개념의 추상화다. 수치 해법, 비선형 동역학, 그리고 확률론은 분리된 기술이 아니라 하나의 언어로 자연과 금융을 기술하는 방법들이다. 이 연결이 보이기 시작할 때, 수식을 외우는 것이 아니라 이해하고 있는 것이다.

미분방정식 및 동역학 해석