📂수치해석

디리클레 경계 조건이 주어진 열방정식에 대한 초기값 문제의 수치해석적 풀이

예제 ¹

$$\begin{cases} u_{t} = \gamma u_{xx} \\ u(t,0) = u(t,l) = 0 \\ u(0,x) = f(x) \end{cases}$$

주어진 문제는 대수적 풀이가 있을 정도로 쉽고 간단하지만, 미분방정식을 푸는 방법으로써의 수치해석을 왜 배우는지 명쾌하게 알려주는 예시가 되기도 한다. 단순히 $y ' = f(x,y)$ 꼴의 미분방정식을 푸는 게 편미분방정식의 풀이로도 이어지는 것이다.

풀이

대수적 풀이에서와 달리 $\gamma$ 의 값이 중요한 것은 아니므로 편의 상 $\gamma = 1$ 이라고 하자.

$$x_{j} = \delta j$$ 선공간^{line space} $x$ 의 스텝사이즈를 $\delta$ 로 두면 $j= 0, 1, \cdots , m$ 에 대해 위와 같이 나타낼 수 있다.

$x$ 에 대한 이계미분의 근사로써 $$U_{xx} \approx {{ U(t, x_{j+1}) - 2 U(t, x_{j}) + U(t, x_{j-1}) } \over { \delta^2 }}$$ 라고 하면 오직 $t$ 에 대한 식 $$u_{t} (t, x_{j}) = {{ u(t, x_{j+1}) - 2 u(t, x_{j}) + u(t, x_{j-1}) } \over { \delta^2 }}$$ 을 $m-1$ 개만큼 얻을 수 있다. 이것을 행렬로 나타내보면 $$\begin{bmatrix} u_{t} (t, x_{1}) \\ u_{t} (t, x_{2}) \\ u_{t} (t, x_{3}) \\ \vdots \\ u_{t} (t, x_{m-2}) \\ u_{t} (t, x_{m-1}) \end{bmatrix} = {{1} \over {\delta^2}} \begin{bmatrix} -2 & 1 & 0 & \cdots & 0 & 0 \\ 1 & -2 & 1 & \cdots & 0 & 0 \\ 0 & 1 & -2 & \cdots & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & -2 & 1 \\ 0 & 0 & 0 & \cdots & 1 & -2 \end{bmatrix} \begin{bmatrix} u (t, x_{1}) \\ u (t, x_{2}) \\ u (t, x_{3}) \\ \vdots \\ u (t, x_{m-2}) \\ u (t, x_{m-1}) \end{bmatrix}$$ 이고, 지면을 아끼기 위해 $$\mathbf{u} ' (t) :=\begin{bmatrix} u_{t} (t, x_{1}) \\ u_{t} (t, x_{2}) \\ u_{t} (t, x_{3}) \\ \vdots \\ u_{t} (t, x_{m-2}) \\ u_{t} (t, x_{m-1}) \end{bmatrix}$$

$$\Lambda := {{1} \over {\delta^2}} \begin{bmatrix} -2 & 1 & 0 & \cdots & 0 & 0 \\ 1 & -2 & 1 & \cdots & 0 & 0 \\ 0 & 1 & -2 & \cdots & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & -2 & 1 \\ 0 & 0 & 0 & \cdots & 1 & -2 \end{bmatrix}$$

$$\mathbf{u} (t) := \begin{bmatrix} u (t, x_{1}) \\ u (t, x_{2}) \\ u (t, x_{3}) \\ \vdots \\ u (t, x_{m-2}) \\ u (t, x_{m-1}) \end{bmatrix}$$ 라고 두면 $\mathbf{u} ’ (t) = \Lambda \mathbf{u} (t)$ 이다. 이제 시간 $t$ 의 스텝사이즈를 $h$ 로 두어 $t_{k} = hk$ 라 하고 여러가지 메소드를 적용해서 수치적인 솔루션을 구하면 된다.

■

교훈: 스티프한 미분방정식

그러나 이것은 어디까지나 이론적인 이야기고, 실제로 수치적인 해를 구하는데에는 상당한 어려움이 있다.

$\Lambda$ 의 고유값을 알기 위해 $\Lambda - \lambda I_{m-1}$ 의 행렬식을 구해보자.

삼중대각행렬의 행렬식: 삼중대각행렬 $X_{n} := \begin{bmatrix} x & 1 & 0 & \cdots & 0 & 0 \\ 1 & x & 1 & \cdots & 0 & 0 \\ 0 & 1 & x & \cdots & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & x & 1 \\ 0 & 0 & 0 & \cdots & 1 & x \end{bmatrix}$ 에 대해 $\displaystyle | x_{n}| = U_{n} \left( {{x} \over {2}} \right)$

$$\Lambda - \lambda I_{m-1} = {{1} \over { \delta^2 }} \begin{bmatrix} -2 - \delta^2 \lambda & 1 & 0 & \cdots & 0 & 0 \\ 1 & -2 - \delta^2 & 1 & \cdots & 0 & 0 \\ 0 & 1 & -2 - \delta^2 & \cdots & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & -2 - \delta^2 & 1 \\ 0 & 0 & 0 & \cdots & 1 & -2 - \delta^2 \end{bmatrix}$$ 은 $(m - 1) \times (m-1)$ 삼중대각행렬이므로, 행렬식은 $$\det \left( \Lambda - \lambda I_{m-1} \right) = {{1} \over { \delta^2 }} U_{m-1} \left( {{-2 - \delta^2 \lambda } \over {2}} \right)$$ 와 같이 구해진다. 여기서 $U_{m-1}$ 은 $m-1$ 차 제2종 체비셰프 다항함수를 의미한다.

제2종 체비셰프 다항함수의 근과 대칭성:

• [1]: $\displaystyle U_{n} (X)$ 의 근은 $\displaystyle x_{k} = \cos \left( {{k} \over {n+1}} \pi \right)$, $k=1, \cdots , n$
• [2]: $U_{n} (-x) = (-1)^{n} U_{n} (x)$

괄호에서 $-1$ 을 빼내면 $${{(-1)^{m-1}} \over { \delta^2 }} U_{m-1} \left( {{2 + \delta^2 \lambda } \over {2}} \right)$$ 이고, 고유값을 구하기 위해선 $$U_{m-1} \left( {{2 + \delta^2 \lambda } \over {2}} \right)$$ 만 신경쓰면 충분하다. $\displaystyle U_{m-1} \left( {{2 + \delta^2 \lambda } \over {2}} \right)$ 의 $j$ 번째 근은 $j=1, \cdots , (m-1)$ 에 대해 $${{2 + \delta^2 \lambda_{j}} \over {2}} = \cos \left( {{j} \over {m}} \pi \right)$$ 이고, 양변에 $2$ 을 곱하면 $$\begin{align*} \displaystyle 2+\delta^2 \lambda_{j} =& 2 \cos \left( {{j} \over {m}} \pi \right) \\ =& 2 \cos \left( 2 \cdot {{j} \over {2m}} \pi \right) \\ =& 2 \left[ 1 - 2 \sin^2 \left( {{ j \pi } \over { 2m }} \right) \right] \end{align*}$$ 따라서 고유값은 $$\lambda_{j} = - {{4} \over { \delta^2 }} \sin^2 \left( {{j \pi } \over {2m}} \right)$$ 과 같이 구해진다. 그런데 $$\begin{align*} \left| \lambda_{m-1} \right| =& \left| {{4} \over { \delta^2 }} \sin^2 \left( {{(m-1) \pi } \over {2m}} \right) \right| \\ \approx & {{4} \over {\delta^2}} \\ =& 4 m^2 \end{align*}$$

이므로, $\lambda_{j}$ 는 $\displaystyle m = {{1} \over {\delta}}$ 이 크면 클수록 고유값의 크기도 엄청나게 커진다. 수치해를 가능한 한 정확히 알기 위해 $\delta$ 를 작게 잡아 $[0,1]$ 을 많이 쪼갰다면 계산이 무척 어려워지는 것이다. 가령 $m=100$ 으로만 잡아도 $\mathbf{u} ’ (t) = \Lambda \mathbf{u} (t)$ 은 스티프^stiff 한 문제가 되어서 수치적으로 풀기엔 꽤 어렵다.

같이보기

• 대수적 풀이