[수치해석] Partial Differential Equation (2) Matrix stability analysis

 

앞선 글에서는 PDE를 공간에 대해 discretization을 적용해 system of ODEs로 바꾸고 그 결과 얻은 matrix의 eigenvalue를 확인했다.

 

https://normal-engineer.tistory.com/511

 

 

system of ODEs에 대해서도 이미 다뤘고 어떤 scheme을 적분에 적용하느냐에 따라서 numerical solution의 stability가 달라진다. 마찬가지로 PDEs에서 scheme의 stability를 분석한다.

 

 

Stability 분석하는 방법은 세 가지가 있다.

 

 

1) Matrix stability analysis
앞서 했던 것처럼 eigenvalue를 구하지만 TDMA가 아니면 Computational cost가 크다.
2) Von Neumann stability analysis
3) Modified wavenumber analysis

 

 

 

이 글에서는 Matrix stability analysis를 다룬다.

Matrix stability analysis

1) Wave equation

다음과 같은 wave equation이 있다고 한다.

$$\frac{\partial u}{\partial t}+c\frac{\partial u}{\partial x}=0$$

 

$$
\frac{\partial u_i}{\partial t} + \frac{c}{2\Delta x}(u_{i+1} - u_{i-1}) = 0
$$

과 같이 중심차분으로 대체한다.

 

* $j = 1$

$$
\frac{\partial u_1}{\partial t} = \frac{c}{2\Delta x} u_0 - \frac{c}{2\Delta x} u_2
$$

* $j = 2$

$$
\frac{\partial u_2}{\partial t} = \frac{c}{2\Delta x} u_1 - \frac{c}{2\Delta x} u_3
$$

* $j = N-1$

$$
\frac{\partial u_{N-1}}{\partial t} = \frac{c}{2\Delta x} u_{N-2} - \frac{c}{2\Delta x} u_N
$$

 

 

이를 행렬 형태로 만들면 다음과 같다.

 

$$
\frac{d}{dt}
\begin{bmatrix}
u_1 \\
u_2 \\
\vdots \\
u_{N-1}
\end{bmatrix}
=
\frac{c}{2\Delta x}
\begin{bmatrix}
0 & -1 & 0 & \cdots & 0 \\
1 & 0 & -1 & \cdots & 0 \\
\vdots & \ddots & \ddots & \ddots & \vdots \\
0 & \cdots & 1 & 0
\end{bmatrix}
\begin{bmatrix}
u_1 \\
u_2 \\
\vdots \\
u_{N-1}
\end{bmatrix}
+
\begin{bmatrix}
\frac{c}{2\Delta x} u_0 \\
0 \\
\vdots \\
-\frac{c}{2\Delta x} u_N
\end{bmatrix}
$$

이전 글을 통해서 우리는 이 matrix의 eigenvalue를 구하는 공식을 알고 있다.

(참고 : $A = B[a\ b\ c]^{n \times n}_{\text{대칭 삼각행렬}}, d_j = \frac{j\pi}{n+1}, \quad j=1,\ldots,n $
$\lambda_j = b + 2 \sqrt{ac} \cos d_j$)

 

따라서 다음 matrix가 주어질 때 eigenvalue는 다음과 같다.

$$
\frac{dU}{dt} = A U, \quad
A = \frac{c}{2\Delta x} B, \quad
B = 
\begin{bmatrix}
1 & 0 & -1 \end{bmatrix}
\Rightarrow
\lambda_j = \pm i \frac{c}{2\Delta x} \cos \frac{j\pi}{N}
$$

 

여기서 알 수 있는 것은 eigenvalue가 허수라는 것이다. 즉, 사인파 해를 가지는 scheme이다. 즉, solution이 stable하게 나올 수 없으므로 high order Runge-Kutta를 사용하거나 implict Euler method를 쓴다.

 

 

 

2) Heat equation

 

 

$$\frac{\partial \phi}{\partial t} = \alpha \frac{\partial^2 \phi}{\partial x^2}$$

 

다음 Heat equation을 semi-discretization을 위해 2CD를 적용하고 Explicit Euler를 적용한다.

 

 

$$
\lambda_i = -2 + 2 \cos \frac{j\pi}{N} \quad (j = 1, 2, \dots, N)
$$

 

 

$$
\begin{aligned}
\Phi^{n+1} &= \Phi^n + \left[ \frac{\alpha}{\Delta x^2} \Delta t \cdot B  \begin{bmatrix} 1 & -2 & 1 \end{bmatrix} \right] \Phi^n \\
&= \left[ I + \frac{\alpha \Delta t}{\Delta x^2} B  \begin{bmatrix} 1 & -2 & 1 \end{bmatrix} \right] \Phi^n \\
\Rightarrow \Phi^n &= \left[ I + \frac{\alpha \Delta t}{\Delta x^2} B \begin{bmatrix}1 & -2 & 1\end{bmatrix} \right]^n \Phi^0
\end{aligned}
$$

이 방정식에서 $$\left[ I+\frac{\alpha \Delta t}{\Delta x^2}\lambda_i \right]^n \Phi^0$$이기 때문에 앞의 matrix eigenvalue를 분석해야한다.

 

$\lambda_{i}=(-2+2\cos\frac{j\pi}{N})\frac{\alpha}{\Delta x^2}$이라는 것을 $B[1,-2,1]$에서 알고 있기 때문에

$|\lambda_{max}|=\left| -\frac{4 \alpha }{\Delta x^2}\right|$

 

따라서 Explicit Euler 방법을 이용할 때 timestep $\Delta t$는 다음을 만족해야한다.

 

$\Delta t \leq \frac{2}{|\lambda_{max}|} =\frac{(\Delta x)^2}{2\alpha}$.

 

여기서 알 수 있는 사실은 $\alpha$가 클수록 확산계수가 크다는 것이고, 시간에 따라 변화가 빠르기 때문에 그만큼 timestep의 크기가 작아져야한다.