미분방정식은 응용수학, 물리학, 공학에서 가장 많이 나오는 수학적 문제입니다. 특히나 물리학과 공학에서는 미분 방정식이 거의 전부라고 해도 좋을 정도로 미분 방정식은 지배적인 문제입니다. 자연을 기술하는 거의 모든 방정식이 미분 방정식의 형태로 쓰여질 수 있기 때문에, 물리학과 공학에서 미분 방정식을 자주 볼 수 있습니다. 예를들면, 고전 역학을 기술하는 뉴턴의 운동 방정식, 전자기 현상을 기술하는 맥스웰 방정식, 열과 관련된 현상을 설명하는 열 방정식 등이 모두 미분 방정식의 형태로 주어집니다. 파동 방정식, 나비에-스토크스 방정식은 뉴턴의 운동 방정식으로 부터 유도된 방정식이기에 당연히 미분 방정식입니다.
따라서 자연 현상을 해석하고 예측하기 위해서는 상황에 맞는 미분 방정식을 선택하고, 경계 조건을 정하고, 그 다음엔 문제를 풀면 됩니다. 주로 학부 과정의 교과서에서는 종이와 펜으로 풀 수 있는 해석적인 방법으로 문제를 풀지만, 실제 자연 현상을 최대한 비슷하게 기술하는 문제를 풀거나 대학원 과정 수준에서는 해석적인 방법으로 문제를 풀 수 없는 경우가 대부분입니다. 이런 경우에는 어쩔 수 없이 컴퓨터의 도움을 받아서 수치적인 방법으로 문제를 풀어야 합니다.
유한요소법(Finite Element Method)는 컴퓨터를 이용한 수치적인 방법으로 미분 방정식을 풀기 위해 고안된 방법론입니다. (아래에서는 미분 방정식의 입력 변수가 시간과 공간좌표라고 가정 하겠습니다) 미분 방정식에서는 함수를 시간이나 공간에 따라서 미분을 하게 되는데, 이는 시간과 공간이 연속적이라는 것을 기본적으로 가정합니다. 어떤 값이 유한한 구간에서 연속이면 그 값의 갯수는 무한개가 되는데, 컴퓨터로는 무한개의 값에 대한 계산이 불가능 하기 때문에, 유한요소법에서는 시간과 공간 변수를 유한한 크기의 작은 조각으로 불연속화 합니다. 이렇게 하면 시간과 공간의 값은 유한한 갯수를 갖게 되며, 이 유한한 갯수의 값들을 이용하여 다양한 계산을 수행하여 답을 얻게 됩니다.
유한요소법에는 다양한 형태의 하위 방법이 존재합니다. 이 중에서 개념적으로 가장 쉬운 방법론은 유한분차법(Finite Difference Method)입니다. 유한분차법은 함수의 미분을 테일러 급수를 이용하여 함수값의 차이로 근사합니다. 이 경우, 미분 방정식은 연립 방정식의 형태로 바뀌고, 이 연립 방정식을 대수적인 방법을 이용하여 풀게 됩니다. 즉, 시간과 위치 $t \in [t_{min}, t_{max}]$과 $x \in [x_{min}, x_{max}]$에 대하여 $(x, t)$를 $(x_n, t_m), n, m =1, 2, 3, 4, ...$와 같이 유한개의 값을 갖도록 불연속화 시키고 이때의 함수값 $f(x_n, t_m) = f_{n, m}$ 을 계산합니다. 그 다음엔 풀려고 하는 미분 방정식으로 부터 $\{ x_n, t_m, f_{n,m}\}$ 이 만족해야하는 조건을 찾고, 이 조건을 만족하는 $f_{n,m}$을 찾는 과정을 거칩니다.
역시, 설명만해서는 이해하기 어렵고, 가장 간단한 형태의 미분 방정식을 직접 유한분차법으로 풀어 보도록 하겠습니다. 매우 간단한 예시이지만, 유한분차법 더 나아가 좀 더 상위 개념인 유한요소법의 기본적인 철학과 방법론을 이해할 수 있습니다.
유한분차법으로 1차원 헬름홀츠 방정식 풀기
헬름홀츠 방정식은 물리학과 공학 전반에서 볼 수 있는 방정식으로 특히 파동, 진동을 다루는 분야에서 자주 접할 수 있습니다. 일반적인 공간 차원에서 헬름홀츠 방정식은
$$(\nabla^2 + k^2) f(\vec{r}) = 0$$
으로 주어지고, 이번에 풀어 볼 가장 간단한 공간 차원이 1차원인 경우에는
$$\Big(\frac{d^2}{dx^2} + k^2 \Big)f(x) = 0$$
와 같이 간단히 주어집니다. 미분 방정식에는 경계 조건이 있어야 하는데, 역시나 가장 간단한 경계 조건인 디리클레 경계 조건인 $x \in [0, L]$ 일 때, $f(0) = f(L) = 0$으로 두도록 하겠습니다.
과정-1 : 입력 변수, 즉 공간 변수를 불연속화 하기
이 문제의 입력 변수는 1차원 공간 변수 $x$ 입니다. $x$은 당연히 $0$과 $L$ 사이의 값을 갖는 연속적인 실수 입니다. 앞서 설명한 대로, 유한분차법에서는 입력 변수의 공간을 유한개로 불연속화 시킵니다. $[0, L]$의 공간을 $N$개의 조각으로 나누면 $n$번째 조각의 점의 $x$값은 $x_n = \frac{L}{N}n$이 됩니다. $x=0$과 $x=L$에서의 값을 모두 포함 할 수 있어야 하기 때문에 $n$값은 $0, 1, 2, 3, ..., N$ 이렇게 $N+1$개의 값을 갖습니다.
위 예시는 $N=7$ 인 경우 입니다. $N$ 등분된 공간의 각 점에서의 $x$값은 $\{ x_0, x_2, ..., x_6, x_7\}$의 값을 갖습니다.
과정-2 : 미분 방정식을 분차 방정식으로 바꾸기
테일러 전개
$$f(x+h) = f(x) + f'(x)h + ...$$
$$f(x-h) = f(x) - f'(x)h + ...$$
로 부터, 위 식에서 아래식을 빼고 정리하면
$$f'(x) = \frac{f(x+h) - f(x-h)}{2h}$$
를 얻습니다. 여기서 $x$값은 불연속화된 점 $x_n$, $h = \frac{L}{N}$ 라고 하면,
$$f'(x_n) = \frac{f(x_{n+1}) - f(x_{n-1})}{\frac{2L}{N}}$$
이 됩니다. 즉 $x_n$에서의 함수의 미분값을 인접한 점인 $x_{n+1}, x_{n-1}$에서의 함수값 $f(x_{n+1}), f(x_{n-1})$으로 쓸 있는 것 입니다. 미분값은 함수의 순간 기울기인데, 이를 $[x_{n-1}, x_{n+1}]$ 구간의 평균 기울기로 근사한 것 과 같습니다. 위 식은 $\frac{L}{N}$이 작아질 수록, 즉 $N$이 커질 수록 정확해 집니다.
헬름홀츠 방정식에서는 1계 미분이 아니라, 2계 미분이 나오기 때문에, 위와 같은 방식으로 2계 미분이 어떻게 근사 될 수 있는지를 구해야 합니다. 테일러 전개에서 2계 미분이 나올 때 까지 근사를 하면 되는데,
$$f(x+h) = f(x) + f'(x)h + \frac{1}{2}f''(x)h^2$$
$$f(x-h) = f(x) - f'(x)h + \frac{1}{2}f''(x)h^2$$
로 부터, 위 식과 아래식을 더하고 정리하면,
$$f''(x) = \frac{f(x+h) + f(x-h) - 2f(x)}{h^2}$$
를 얻습니다. 여기서 $x$값은 불연속화된 점 $x_n$, $h = \frac{L}{N}$ 라고 하면,
$$f''(x_n) =\frac{f(x_{n+1}) + f(x_{n-1}) - 2f(x)}{\frac{L^2}{N^2}}$$
이 됩니다. 역시 위와 같이, $x_n$에서의 함수의 2계 미분값을 인접한 점인 $x_{n+1}, x_{n-1}$에서의 함수값 $f(x_{n+1}), f(x_{n-1})$으로 쓸 있는 것 입니다. 위 식의 정확도는 $\frac{L}{N}$이 작아질 수록, 즉 $N$이 커질 수록 정확해 집니다.
헬름홀츠 방정식 전체를 쓰면,
$$-f(x_{n+1}) - f(x_{n-1}) + 2f(x_n) = \frac{L^2}{N^2} k^2 f(x_n)$$
$$\text{경계 조건 : } f(x_0) = f(x_N) = 0$$
이 됩니다. 편의상 미분 부분과 $k^2 f(x)$ 부분을 좌변과 우변으로 분리하였고, 좌변을 최대한 간단히 하기 위해서 $f(x)$값을 제외한 모든 상수들은 우변으로 넘겼습니다. 위 식에서 $f(x_n)$을 $f_n$로 간단히 표현하면,
$$-f_{n+1} - f_{n-1} + 2f_n = \frac{L^2}{N^2} k^2 f_n$$
$$f_0 = f_N = 0$$
이 되고, 위 식을 만족하는 $N+1$개의 값 $f_0, f_1, ..., f_N$를 찾으면 됩니다. 함수의 미분으로 주어졌던 문제를 $N+1$개의 변수 $f_n$의 연립 방정식으로 바꾸었습니다.
과정-3 : 연립 방정식의 구조를 파악하기
위에서 식을 하나 밖에 쓰지 않았지만, 위 식은 모든 $n=0, 1, 2, ... N$에서 만족해야 하는 식 입니다. 따라서 $f_0, f_1, ..., f_N$에 대한 연립 방정식이 됩니다. 혹시 아직까지 이해를 하지 못 하신 분을 위해서, 한 번 더 explicit 하게 순서대로 쓰면 아래와 같습니다.
$$2f_1 - f_{2} = \frac{L^2}{N^2} k^2 f_1$$
$$-f_{1} + 2f_2 - f_{3} = \frac{L^2}{N^2} k^2 f_2$$
$$-f_{2} + 2f_3 - f_{4} = \frac{L^2}{N^2} k^2 f_3$$
$$......$$
$$-f_{N_2} + 2f_{N-1} = \frac{L^2}{N^2} k^2 f_{N-1}$$
$f_{n-1}, f_n, f_{n+1}$의 순서를 눈에 보이게 하기 위해서 좌변의 배열을 조금 바꾸었습니다. 또한 $f_0$과 $f_N$은 값이 $0$으로 제한 돼 있기 때문에, 위 연립 방정식에서 완전히 제거 했습니다. 이렇게 쓰면 연립 방정식이라는 것이 눈에 바로 들어 옵니다. 그리고 위 식의 대칭적인 구조 역시 확인할 수 있습니다. 위 식을 한 번 더 정리 하면 더 좋은데요, 편의상 복잡하게 보이는 $\frac{L^2}{N^2} k^2$ 를 $\lambda$로 간단히 하고 우변을 좌변으로 이항하여 우변은 $0$으로 만들겠습니다. 이렇게 하면,
$$(2-\lambda)f_1 - f_{2} = 0$$
$$-f_{1} + (2-\lambda)f_2 - f_{3} = 0$$
$$-f_{2} + (2-\lambda)f_3 - f_{4} = 0$$
$$......$$
$$-f_{N-2} + (2-\lambda)f_{N-1}= 0$$
와 같이 더 깔끔하게 정리 됩니다. 선형 연립 방정식은 행렬을 이용하여 표현하면 더 깔끔하게 정리 할 수 있는데요, $\vec{f} = (f_1, f_2, ..., f_{N-2}, f_{N-1})$로 정의하면 위 연립 방정식은
$$(\textbf{A} - \lambda \textbf{I})\vec{f}=
\begin{pmatrix}
2-\lambda&-1\\
-1&2-\lambda &-1\\
&-1&2-\lambda &-1&...\\
&&...\\
&&&&-1&2-\lambda &-1\\
&&&&&-1&2-\lambda\\
\end{pmatrix}
\begin{pmatrix}
f_1\\
f_2\\
f_3\\
...\\
f_{N-2}\\
f_{N-1}\\
\end{pmatrix}
=
\begin{pmatrix}
0\\
0\\
0\\
...\\
0\\
0\\
\end{pmatrix}$$
가 됩니다. 이 선형 연립 방정식이 자명하지 않은, 즉 $f_n=0$이 아닌, 해를 갖기 위해서는 $\det|\textbf{A} - \lambda \textbf{I}|$=0를 만족해야 합니다. 선형 대수학을 공부한 적이 있으신 분이라면 이 같은 형태의 행렬 방정식 및 $\det|\textbf{A} - \lambda \textbf{I}|$=0 조건에 대해서 매우 익숙하실 것 입니다. 바로 행렬의 고유 벡터와 고유 값을 구하는 것 인데요. $\lambda$는 행렬의 고유값, 그리고 $\vec{f}$는 고유 벡터가 됩니다. $\frac{L^2}{N^2} k^2$ 를 $\lambda$라고 정의 했었으니, $k = \frac{N}{L}\sqrt{\lambda}$가 됩니다. $\lambda$가 특정한 값(고유값)만을 갖기 때문에 $k$ 역시 특정한 값만 갖게 됩니다. 헬름홀츠 방정식에 대해 소개했던 이전 포스팅에서 미분 방정식의 경계 조건에 의해서 $k$가 특정한 값만 갖게 된다고 설명하였습니다. 미분 방정식을 해석적인 방법으로 풀었을 때와, 이 포스팅에서 처럼 유한분차법을 이용하여 대수적인 방식으로 풀었을 때 모두 $k$값은 경계조건에 의해서 특정한 값을 갖게 되는 동일한 수학적 성질을 보여줍니다.
과정-4 : 연립 방정식을 풀기. 즉 행렬 $\textbf{A}$의 고유값과 고유벡터를 구하기
행렬 대각화를 하는 프로그래밍 라이브러리를 이용해서 위 행렬을 대각화 할 수 있지만, 위 행렬 $\textbf{A}$는 대각성분을 중심으로 대칭성이 있는 간단한 행렬이기 때문에, 해석적인 방법으로 고유값과 고유벡터를 구할 수 있습니다. 고유값은 $n=1, 2, 3, ... N-1$에 대해서, $\lambda_n = 2-2\cos(\frac{\pi n}{N})$으로, 이에 대응되는 고유 벡터 $\vec{f}_n$의 $i$번째 성분($i=1, 2, 3, ..., N-1$)은
$$f_n(i) = \sin \Big(\frac{\pi n}{N} i\Big)$$
로 주어집니다. 위 식은 $f(0) = f(N) = 0$을 만족시키기 때문에 $i$의 범위를 $i=0, 1, ..., N$으로 확장할 수 있습니다.
$N=20$ 일때, $f_1(i), f_2(i)$의 그래프 입니다. $i$는 $i=0, 1, 2, ..., 20$과 같이 정수값에서만 정의되는 값이기 때문에 위 점과 점 사이에는 아무것도 없어야 하지만, 그래프를 보는 편의상 인접한 점과 점 사이를 선으로 연결 하였습니다.
$\lambda$를 원래값인 $k$에 대한 값으로 환원하면,
$$k = \frac{N}{L} \sqrt{2-2 \cos \Big(\frac{\pi n }{N}\Big)} = \frac{N}{L} 2 \sin \Big(\frac{\pi n}{2N} \Big)$$
이 됩니다.
과정-5 : 해석적인 방법의 결과와 비교하기
지난 포스팅에서 이미 풀어본 바와 같이 1차원 헬름홀츠 방정식은 해석적인 방법으로 간단하게 풀 수 있고, 그 해는
$$f_n(x) = \sin \Big( \frac{\pi n }{L} x\Big)$$
$$k = \frac{\pi n}{L}$$
이 됩니다. 고유 함수를 비교하면, 해석적입 방법의 해에서 $x = \frac{L}{N}i$를 대입하면 정확히 유한분차법 방법의 해와 일치함을 확인할 수 있습니다. 고유 값$k$를 비교 하면, 함수의 개형이 좀 많이 달라 보입니다. 유한분차법은 $N$이 매우 커지면 커질 수록, 결과가 정확해 집니다. 즉, $N \rightarrow \infty$ ($n \ll N$)를 대입하여 유한분차법에서 얻어진 고유값을 전개하면,
$$k = \frac{N}{L} 2 \sin \Big(\frac{\pi n}{2N} \Big) \approx \frac{N}{L} 2 \frac{\pi n}{2N} = \frac{\pi n}{L}$$
이 되어, 유한분차법으로 얻을 수 있는 $k$값과 해석적입 방법으로 얻은 $k$값이 일치함을 확인할 수 있습니다.
해석적인 해 $k$와 유한분차법의 해 $k$를 조금 더 비교해 보도록 하겠습니다. $\sin(x) = x - \frac{1}{6}x^3$ 이므로, 이 식을 이용하여 유한분차법의 해 $k$를 전개하고, 보기 좋게 정리하면,
$$k = k_{\text{analytic}} \Big( 1 - \frac{\pi^2}{24} \frac{n^2}{N^2} \Big) $$
이 됩니다. 즉, 비율로 본다면 $\frac{\pi^2}{24} \frac{n^2}{N^2}$ 만큼의 에러가 생기는 것인데, 이 값은 $\frac{1}{N^2}$에 비례합니다. 유한분차법은 공간을 $N$조각으로 나누고, 함수의 미분값을 평균기울기로 근사한 것이기 때문에, 항상 $\frac{1}{N}$과 관련된 에러를 갖게 됩니다.
이번 포스팅에서도 왠지 닭 잡는데 소 잡는 칼을 쓴 것 같다.
해석적으로 정확히 풀 수 있는 방정식을 유한요소법(유한분차법)이라는 다소 귀찮고도 부적확한(에러가 있는) 방법으로 풀었습니다. 하지만 소 잡는 칼을 만들어 놓아야, 나중에 소를 잡을 수 있습니다. 해석적인 방법으로는 도저히 풀 수 없는 복잡한 형태의 미분방정식을 풀어야 할 때, 혹은 방정식 자체는 간단하지만 경계 조건이 매우 복잡한 경우에는 유한요소법과 같이 수치적인 방법이 꼭 필요하게 됩니다. 다음 포스팅 부터는 더 흥미로운 문제를 풀어 보도록 하겠습니다.
