미분방정식을 연립방정식으로, FDM #3 (2계 상미분방정식)

 

이 글은

2022.07.22 - [수학의 재미/아름다운 이론] - 미분방정식을 연립방정식으로, FDM #2

미분방정식을 연립방정식으로, FDM #2

에서 이어집니다.

 

앞선 글에서 1계 도함수가 포함된 미분 방정식을 FDM을 이용하여 해결하는 방법을 알아보았습니다.

 

이 글에서는 2계도함수가 포함된 미분방정식을 풀어보겠습니다. 

 

$x$의 범위를 $0\leq x\leq1$로 하는 미분방정식

$$f''(x) + p f'(x) + q f(x)= 0~,~ f(0)=y_0 ,  f(1) = y_1 \tag{1}$$

인 함수를 생각해 봅시다. $p,q, y_0, y_1$은 모두 실수입니다.

 

미분방정식의 풀이

식(1)을 풀기 위하여 $f(x) = e^{\alpha x}$라고 두어 봅시다. 그러면 식(1)은

$$(\alpha^2 +p \alpha +q ) e^{\alpha x} =0$$

이 되므로 $\alpha$는 2차 방정식

$$t^2 + pt +q=0 \tag{2}$$

의 근입니다. 따라서 다음의 경우들로 나뉩니다.

 

$p^2 - 4q>0$ 일 때,

이 때는 식(2)이 서로 다른 두 근 $\alpha_1 , \alpha_2$를 가집니다. 따라서 식(1)의 해는

$$f(x) = c_1 e^{\alpha_1 x} + c_2 e^{\alpha_2 x}$$

로 쓸 수 있습니다. $f(0)=y_0 , f(1)= y_1$을 이용하면

$$c_1 +c_2 = y_0 ~,~ e^{\alpha_1} c_1 + e^{\alpha_2}c_2 =y_1$$

을 만족하죠. 이를 풀면

$$c_1 = \frac{e^{\alpha_2}y_0 - y_1}{e^{\alpha_2} - e^{\alpha_1}} ~,~ c_2 = \frac{-e^{\alpha_1}y_0 + y_1}{e^{\alpha_2} - e^{\alpha_1}}$$

입니다.

 

 

$p^2 - 4q=0$ 일 때,

이 때는 식(2)이 중근 $\alpha$를 갖습니다. 따라서

$$\alpha = -\frac{p}2$$

입니다. 이때 식(1)의 해는 

$$f(x) = c_1 e^{\alpha x} + c_2 x e^{\alpha x}$$  

라 쓸 수 있습니다.

 

$xe^{\alpha x}$ 가 근이 되는 이유

$$(x e^{\alpha x} )' = (\alpha x+1) e^{\alpha x}$$ 이고 $$(x e^{\alpha x} )'' = (\alpha^2 x+ 2\alpha)e^{\alpha x}$$ 이므로 식(1)의 좌변은 $$(\alpha^2 x+ 2\alpha+p(\alpha x+1) +qx)e^{\alpha x} = (\alpha^2 +p\alpha+q)x +(2\alpha+p)=0$$
이므로 근이 됩니다.

 

$f(0)=y_0 , f(1)= y_1$ 을 이용하면 $f(0) = c_1, f(1)= e^\alpha (c_1+c_2)$이므로

$$c_1 = y_0 ~.~ c_2 = y_1 e^{-\alpha} - y_0$$

가 됩니다.

 

$p^2 - 4q<0$ 일 때,

이 때는 식(2)이 허근 $\alpha = z_1 \pm z_2 i$를 갖습니다.

$$e^{\alpha x} = e^{ (z_1 \pm z_2 i)x} = e^{z_1 x} (\cos (z_2 x) \pm i \sin(z_2 x))$$

이므로 

$$e^{z_1 x} \cos(z_2 x) ~,~ e^{z_1 x} \sin(z_2 x)$$

가 식(1)의 두 근이 되고, 따라서 일반 해는

$$f(x) = c_1 e^{z_1 x} \cos(z_2 x) + c_2 e^{z_1 x} \sin(z_2 x)$$

의 형태로 쓰임을 알 수 있습니다.

$f(0)=y_0 , f(1)= y_1$ 을 이용하면 $f(0) = c_1, f(1)= c_1 e^{z_1}\cos(z_2)+c_2 e^{z_1}\sin(z_2) =y_1$이므로

$$c_1= y_0~,~ c_2 = \frac{y_1}{e^{z_1}\sin(z_2)} - \frac{y_0}{\tan(z_2)}$$

입니다.

 

 

 

FDM을 이용한 풀이

 

$x$의 범위를 $0\leq x\leq1$로 하는 미분방정식 식(1)인

$$f''(x) + p f'(x) + q f(x)= 0~,~ f(0)=y_0 ,  f(1) = y_1 \tag{1}$$

을  $f'(x)$에 대한 forward difference와 backward differece 각각을 사용하여 풀어보겠습니다.

이계도함수 $f''(x)$는 central difference로 정의합니다. 

 

 

풀이를 위해 $x$의 최댓값을 $x_{\max}$라 합시다. 그리고 $[0, x_{\max}]$을 $N$등분을 하여

$$0=x_0 < x_1 <x_2 <\cdots < x_{N-1} < x_N = x_{\max}$$

로 나눕니다. 총 $N+1$개의 $x$ 수열이 생기겠죠.

 

$x_{i+1}$과 $x_i$사이의 간격을 $h$라 합시다. 그리고 

$$f(x_i)= u_i$$

라 정의합니다(FDM을 풀 땐 $u$라는 기호를 많이 쓰더군요) 

그러면 점 $x_i$에서의 이계도함수는

$$f''(x_i ) = \frac{f(x_{i}+h)-2f(x_i) +f(x_{i}-h)}{h^2} = \frac{f(x_{i+1})-2f(x_i) +f(x_{i-1})}{h^2} =\frac{u_{i+1}-2u_i +u_{i-1}}{h^2}$$

입니다.

 

 

Forward Difference

$x_i$에서 $f'$을 forward difference로 구하면

$$f'(x_i ) = \frac{f(x_i +h) -f(x_i)}h = \frac{f(x_{i+1})-f(x_i)}h = \frac{u_{i+1}-u_{i}}h$$

입니다.

따라서 식(1) 은

 

$$\frac{u_{i+1}-2u_i +u_{i-1}}{h^2} + p \frac{u_{i+1}-u_{i}}h + q u_i =0$$ 이고

$u_0 = y_0~,~ u_N =y_1$입니다.

 

따라서 정리하면

$$\frac1{h^2} u_{i-1} -\left( \frac2{h^2} + \frac{p}h -q \right) u_i + \left( \frac1{h^2} +\frac{p}h\right) u_{i+1}=0~,~ i\geq 1  \tag{2}$$

입니다. 간략하게

$$\begin{align} a_{i} &=  \frac1{h^2} \\ b_i  & = -\left( \frac2{h^2} + \frac{p}h -q \right) \\ c_i &=  \frac1{h^2} +\frac{p}h \end{align}$$

라 하면,  식(2)는

$$a_i u_{i-1} + b_i u_i + c_i u_{i+1} =0$$

입니다.

 

 

특별히 $i=1$일 때는
$$b_1 u_1 + c+1 u_2 = -a_1 y_0$$  
이고, $i=N-1$일 때는
$$a_{N-1} u_{N-2} + b_{N-1} u_{N-1} = - c_{N-1} y_1$$
입니다. 따라서 행렬식으로 표현하면


$$\pmatrix { b_1 & c_1 &   &  &   &\\ a_2 & b_2 & c_2 &  & &\\  & a_3 & b_3 &\ddots &  &\\  &  & &\ddots & \ddots & \\  & &  & a_{N-2} & b_{N-2} & c_{N-2} \\  &  &  & & a_{N-1} & b_{N-1} } \pmatrix { u_1 \\ u_2 \\ u_3 \\ \vdots \\ u_{N-2} \\ u_{N-1} } = \pmatrix { -a_1 y_0 \\ 0\\ 0\\ \vdots\\ 0\\ - c_{N-1} y_1 }$$

 

행렬을 차례대로 $\mathbf{T, u, k}$ 라 하면

$$\mathbf{Tu} = \mathbf{k} \tag{3}$$

입니다.

 

여기서 우선 짚고 넘어갈 개념이 있습니다. 바로 $\mathbf{T}$의 생김새인데요. 이 행렬은 대각원소, 대각선 바로 위, 아래 대각선의 원소를 제외하고는 모두 0 입니다. 이를 tridiagonal matrix라 합니다.  이에 대해서는 다른 글에서 심도 있게 다뤄 보도록 하겠습니다. 어쨌든 식(3)에서

$$\mathbf{u} = \mathbf{T}^{-1} \mathbf{k}$$

입니다. 이렇게 $\mathbf{u}$를 구하면 함수 $f$를 이산적이지만 촘촘히 구한 결과가 나오겠지요.

 

 

Backward Difference

이번에는 $x_{i+1}$에서의 $f'$의 backward difference를 사용하여 구해볼까요?

 

$$f'(x_{i+1}) = \frac{f(x_{i+1}-h)-f(x_{i+1})}{-h} = \frac{u_{i+1}-u_{i}}h$$

이고 식(1)은

 

$$\frac{u_{i+1}-2u_i +u_{i-1}}{h^2} + p \frac{u_{i+1}-u_{i}}h + q u_{i+1} =0$$ 이고

$u_0 = y_0~,~ u_N =y_1$ 입니다.

 

따라서 정리하면

$$\frac1{h^2} u_{i-1} -\left( \frac2{h^2} + \frac{p}h \right) u_i + \left( \frac1{h^2} +\frac{p}h +q \right) u_{i+1}=0~,~ i\geq 1  \tag{2'}$$

입니다. 간략하게

$$\begin{align} a_{i} &=  \frac1{h^2} \\ b_i  & = -\left( \frac2{h^2} + \frac{p}h  \right) \\ c_i &=  \frac1{h^2} +\frac{p}h + q \end{align}$$

라 하면,  식(2')는

$$a_i u_{i-1} + b_i u_i + c_i u_{i+1} =0$$

입니다.

나머지 부분은 forward difference때와 동일합니다. 계수만 달라졌을 뿐이지 수식(3)이 그대로 나오게 됩니다.

 

한 가지 숙제.

 

$x_{i+1}$에서의 difference 라면

$$f''(u_{i+1}) = \frac{u_{i+2}-2u_{i+1}+u_i}{h^2}$$

라 써야 할 것 같은데, 마치 $u_i$에서의 central difference처럼 썼습니다. 그래도 되는 걸까요? 이유는 뭘까요~

 

 

Python으로 풀어보기

 

글의 양이 많아지는 관계로 다음 글에 바로 이어서 설명하도록 하겠습니다. 

2022.07.25 - [수학의 재미/아름다운 이론] - 미분방정식을 연립방정식으로, FDM #4 (2계 상미분방정식)

미분방정식을 연립방정식으로, FDM #4 (2계 상미분방정식)

위 글을 참고하시기 바랍니다.