[수치해석] 수치해석 과정 예시

이번에는 수치해석 과정을 열전달 예시를 통해서 볼 수 있습니다.

열전달 방정식은 이미 solution이 잘 나와있지만 이 예시를 통해서 수치해석을 하는 과정에 대한 이해를 높일 수 있습니다.

 

 

벽에 붙어있는 전도체가 있을 때, 두께는 길이에 비해 매우 작다고 가정하기 때문에 1-Dimensional heat transfer 문제입니다.

 

수치해석은 크게 세 단계로 나눌 수 있습니다.

modeling - transform to algebraic equation - solve algebraic equation

 

이 문제도 같은 과정을 거치도록 하겠습니다.

 

1) modeling

 

위의 시스템을 differential equation(DE)로 만드는 과정입니다.

이 시스템에서 알아내고자 하는 것(solution)은 x에 따른 전도체의 온도(T) 분포입니다.

 

$\frac{\partial T}{\partial t}=k\frac{\partial^2 T}{\partial x^2}-q$

(T : 온도, x : 위치, t : 시간, q : heat source)

이미 열전달(Heat Transfer)을 아는 분들이라면 친숙한 식을 도출해낼 수 있습니다. 열전달을 모르더라도 시스템의 DE가 이렇게 되어있다, 정도로 생각해주시면 됩니다.

 

여기서 T는 종속변수(dependent variable)이고, t와 x는 독립변수(independent variable)입니다.

그러나 위 식은 linear unsteady partial DE입니다.

(독립변수의 제곱이나 세제곱 등이 존재하지 않으므로 linear, 시간에 의존하므로 unsteady, 종속변수가 2개 이상이므로 partial)

 

풀기 어려운 형태이므로 일단은 multiple ordinary DE로 transform한 다음에 algebraic equation으로 변형해주도록 합니다. algebraic equation으로 바꾸기 전에 한 단계를 더 추가하는 셈입니다.

---

2) Transform to multiple ordinary DE

 

온도 T가 시간에 대해서 cyclic하다고 생각하면, 온도 T는 fourier tranform을 적용하여 harmonic form의 조합으로 나타낼 수 있습니다.

 

이 경우 $\hat{T}(\omega)=\int_{-\infty}^{\infty}T(t)e^{-i\omega t}dt$

$T(t)=\frac{1}{2\pi} \int_{-\infty}^{\infty}\hat{T}(\omega) e^{i\omega t}d\omega$

 

fourier transform으로 얻은 것은 특정 주파수 성분에 대한 amplitude 값입니다. 또한 이것을 역으로 transform하면 온도를 주파수 성분들의 합으로 표현할 수 있습니다.

 

$\frac{1}{2\pi}\int_{-\infty}^{\infty}[i\omega \hat{T}(\omega)=k\frac{\partial^2 \hat{T}}{\partial x^2}-\hat{q} ]d\omega$ (이해부족으로 보충 필요)

 

$k\frac{\partial^2 \hat{T}}{\partial x^2}-i\omega \hat{T}(\omega)-\hat{q}=0 \text{ for }\forall \omega$

다음과 같은 ODE로 변형합니다.

---

3) Transform to algebraic equation

 

algebraic equation으로 바로 바꾸기 전에 배경이 되는 이야기를 하고자 합니다.

 

$f'(x)$를 배울 때 보통은 극한 개념을 이용해서 배우게 됩니다. 두 점 사이의 기울기를 구하고 이 두 점 사이의 간격을 극한으로 작게 보냄으로써 미분값을 정의할 수 있습니다.

 

그래서 $f'(x)$는 다양한 방식으로 표현이 될 수 있는데요,

$f'(x)=\underset{\Delta x\rightarrow 0}{\lim}\frac{f(x+\Delta x)-f(x)}{\Delta x}$

$f'(x)=\underset{\Delta x\rightarrow 0}{\lim}\frac{f(x+\Delta x)-f(x-\Delta x)}{2\Delta x}$

$f'(x)=\underset{\Delta x\rightarrow 0}{\lim}\frac{f(x)-f(x-\Delta x)}{\Delta x}$

 

만약 극한을 취하지 않는다면 세 식의 값은 다를 것입니다.

또한 극한을 취한다하더라도, 위의 세 식이 모두 같은 값을 가지려면 전제가 필요합니다. $f(x)$가 continuous하면서 smooth해야합니다.

 

$f(x)$가 continuous하고 differential할 때

현실에서 $\Delta x$를 무한하게 0으로 만들 수는 없기 때문에 finite method를 사용합니다.

 

분야마다 조금씩 다르지만 

finite difference method(FDM)

finite volume method(FVM)

finite element method(FEM)

 

대표적으로 세 가지 방법이 finite method입니다.

 

위에서 하던 예시에 이어서는, FDM 방식을 사용합니다.

 

$\frac{\partial^2 \hat{T}}{\partial x^2}$를 finite method를 통해서 구해볼 것입니다.

 

미분가능성과 연속성이 보장된 $f(x)$는 talyer expansion으로 표현할 수 있습니다.

 

$f(x+\Delta x)=f(x)+\Delta x f'(x)+\frac{\Delta x^{2}}{2!}f''(x)+\frac{\Delta x^{3}}{3!} f'''(x)+\cdots (1)$

$f(x-\Delta x)=f(x)-\Delta x f'(x)+\frac{\Delta x^{2}}{2!}f''(x)-\frac{\Delta x^{3}}{3!} f'''(x)+\cdots (2)$

 

두 식을 합하면((1)+(2))

 

$\frac{f(x-\Delta x)-2f(x)+f(x+\Delta x)}{\Delta x^{2}}=f''(x)+\frac{\Delta x^{2}}{12}f^{(4)}(x)+\frac{\Delta x^{4}}{360}f^{(6)}(x)+\cdots$

 

$\Delta x$가 충분히 작다면 등식의 우변에서 $f''(x)$를 제외한 나머지는 Leading-order error로 볼 수 있습니다.

 

$\frac{f(x-\Delta x)-2f(x)+f(x+\Delta x)}{\Delta x^{2}}=f''(x)+\theta(\Delta x^{2})$

만약 $\Delta x$에 따라 $\frac{1}{2}$배 줄어들면 error는 $\frac{1}{4}$배 줄어듭니다.

 

---

 

이 예시에서 궁금했던 것은 $\hat{T}(x)$이므로 위의 $k\frac{\partial^2 \hat{T}}{\partial x^2}-i\omega \hat{T}(\omega)-\hat{q}=0$ 식을 풀어보겠습니다.

 

위와 같이 전도체를 $x=0$부터 $x=L$까지 grid로 finite하게 나눈 다음에 위 ODE를 적용합니다.

 

$k[\frac{\hat{T}(x_{i}+\Delta x)+\hat{T}(x_{i})+\hat{T}(x_{i}-\Delta x)}{\Delta x^{2}}]-i\omega \hat{T}(x_{i})-\hat{q}(x_{i})=0$

이 식을 각 항에 대해서 정리하면

$\frac{k}{\Delta x^{2}}\hat{T}(x_{i}-\Delta x)+(\frac{k}{\Delta x^{2}}-i\omega)\hat{T}(x_{i})+\frac{k}{\Delta x^{2}}\hat{T}(x_{i}+\Delta x)=\hat{q}(x_{i})$

 

그리고 $x_{i}$를 $x_{L}$, L번째 위치로 바꿔주면서 $x-\Delta x$와 $x+\Delta x$를 $x_{L-1},x_{L+1}$로 바꿔줍니다.

 

$\frac{k}{\Delta x^{2}}\hat{T}(x_{L-1})+(\frac{k}{\Delta x^{2}}-i\omega)\hat{T}(x_{L})+\frac{k}{\Delta x^{2}}\hat{T}(x_{L+1})=\hat{q}(x_{L})$

 

식에서 계수가 상당히 복잡하므로 $\frac{k}{\Delta x^{2}}=\alpha$, $(\frac{k}{\Delta x^{2}}-i\omega)=\beta$로 치환하겠습니다.

 

$\left\{\begin{matrix} \alpha \hat{T}(x_{0})+\beta \hat{T}(x_{1})+\alpha \hat{T}(x_{0})=\hat{q}(x_{1})\\ \alpha \hat{T}(x_{1})+\beta \hat{T}(x_{2})+\alpha \hat{T}(x_{3})=\hat{q}(x_{2})\\ \vdots\\ \alpha \hat{T}(x_{N-2})+\beta \hat{T}(x_{N-1})+\alpha \hat{T}(x_{N})=\hat{q}(x_{N-1})\end{matrix}\right.$

 

위 식은 $i=1$부터 $i=N-1$까지의 식을 순차적으로 배열한 것입니다. 이런 구조를 structured method라고 합니다.

 

참고로 ODE를 풀기 위해서는 문제에서 Boundary condition을 알고 있어야 합니다. 즉 $\hat{T}(x_{0})$과 $\hat{T}(x_{N})$은 알고 있습니다.

 

위 식들은 algebraic equation이므로 matrix form으로 나타낼 수 있습니다.

 

$\begin{bmatrix}\beta & \alpha & 0 & \cdots & 0\\ \alpha & \beta & \alpha & \cdots & 0\\ 0 & \alpha & \beta & \cdots & 0\\  &  &  &  & \end{bmatrix}\begin{bmatrix}\hat{T}(x_{1})\\ \hat{T}(x_{2})\\ \vdots\\ \hat{T}(x_{N-1})\end{bmatrix}=\begin{bmatrix}\hat{q}(x_{1})-\alpha\hat{T}(x_{0})\\ \hat{q}(x_{2})\\ \vdots\\ \hat{q}(x_{N-1})-\alpha\hat{T}(x_{N})\end{bmatrix}$

 

 

위와 같이 대각선을 중심으로 band를 형성한 matrix를 Tri-diagonal matrix(TDMA)라고 합니다. 이와 같은 banded matrix($\leftrightarrow$ sparse matrix)는 계산 cost가 적습니다.

 

$Ax=b \Rightarrow x=A^{-1}b \sim \theta (n^{1})$

같은 matrix라고 해도 위의 matrix의 row 순서를 바꿔서 sparse matrix로 바꾸면 $\sim \theta (n^{3})$입니다.

즉 같은 선형 방정식이어도 구조에 따라서 cost가 다를 수 있습니다.

 

이 식에서 $A^{-1}$이 어떻게 이뤄지는지(마지막 단계)는 다음으로 넘어가겠습니다.

 

---

 

(참고) 수치해석은 interpolation, (수치해석적으로) 미분, 적분 이 세 가지 operation을 통해 이루어집니다.