미분방정식을 연립방정식으로, 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계 상미분방정식)

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