[수치해석] System of ODEs - Boundary value problem

 

Initial value problem과 달리 boundary value problem은 boundary condition이 initial에 있지 않은 경우다.

이럴 경우에는 솔루션을 구할 때 기존의 수치해석 방법을 적용할 수 없다.

 

2nd order ODE가 있다고 하자.

 

$$y^{''}=f(y,y^{'},x)$$

$$y(0)=y_0, y(L)=y_L$$

이 방정식을 풀기 위해서는 $y(0)=0,y^{'}(0)=0$이 주어져야 한다. 그러나 $y(L)=y_L$이 주어졌기 때문에 다른 방법을 적용해야한다.

 

두 가지 방법이 있는데

첫 번째는 shooting method,

두 번째는 direct method

 

첫 번째 방법은 생각나는 솔루션을 적용하는 방법이고

두 번쨰 방법은 finite method를 적용하는 방법이다.

 

1. Shooting method

 

1-1)  BVP를 1st-order ODEs

$u=y, \; v=y^{'}$

$$\begin{aligned}u^{'}=v \\ v^{'}=f(u,v,x) \end{aligned}$$

BCs $$u(0)=y_0, \; u(L)=y_L$$

 

1-2) $v(0)=E$로 Guess한다.

1-3) ODE를 적분한다 until $x=L$.

1-4) $u(L)=y_{L}$인지 확인한다.

1-5) If not, 1-2)로 돌아간다.

 

Linear 방정식이면 2개의 initial guess로 바로 솔루션을 구할 수 있다.

1. $y_1'(0) \rightarrow y_1(x) \rightarrow y_1(L) \neq y_L$
2. $y_2'(0) \rightarrow y_2(x) \rightarrow y_2(L) \neq y_L$

 

인 경우라고 하자.

 

미지수 2개, boundary condition 2개이고 솔루션을 다음과 같이 작성할 수 있다.
$$
y(x) = c_1 y_1(x) + c_2 y_2(x)
$$

경계조건 적용:
$$
\begin{aligned}
y_L &= c_1 y_1(L) + c_2 y_2(L), \\
y_0 &= c_1 y_1(0) + c_2 y_2(0).
\end{aligned}
$$

 

Shooting 시 두 초기해 모두 $y(0)=y_0$ 로 설정하므로

$$
y_1(0) = y_2(0) = y_0 \;\Rightarrow\; c_1 + c_2 = 1.
$$

따라서 $c_2 = 1 - c_1$.

$$
\begin{aligned}
y_L &= c_1 y_1(L) + (1-c_1) y_2(L) \\
    &= (y_1(L) - y_2(L)) c_1 + y_2(L).
\end{aligned}
$$

이 방정식을 풀면 

 

$$
c_1 = \frac{\,y_L - y_2(L)\,}{\,y_1(L) - y_2(L)\,}, \qquad
c_2 = \frac{\,y_1(L) - y_L\,}{\,y_1(L) - y_2(L)\,}.
$$

 

바로 solution을 구할 수 있다.

 

비선형일 경우에는 2번의 initial guess로 해결할 수 없다. 대신 다음 guess를 어떻게 하는지에 대한 방법을 제안할 수 있다.

secant method로 새로운 $y^{'}(0)$를 구할 수 있다.

 

 

$$
M = \frac{y_1(L) - y_2(L)}{y_1'(0) - y_2'(0)}
$$

라고 할 때, 다음에 쓸 새로운 $y^{'}(0)$를 다음과 같이 계산할 수 있다.

$$
y'(0) = y_2'(0) + \frac{y(L) - y_2(L)}{M}
$$

 

이 값을 써서 맞지 않으면 다시 적용한다. 

$$
y_{\alpha+1}'(0) = y_{\alpha}'(0) + \frac{1}{m_{\alpha-1}} \left[ y_L - y_{\alpha}(L) \right]
$$

 

이 경우 $y(L)$이 $y^{'}(0)$에 민감하면 converge하기 어렵다.

 

 

2. Direct method

$$
\left\{\begin{array}{l}
y^{\prime \prime}(x)+A(x) y^{\prime}(x)+B(x) y=f(x) \\
y(0)=y_0, y(L)=y_L
\end{array}\right.
$$

 

다음과 같이 boundary value problem이 주어질 때 

Discretization을 적용해서 식을 전개하면 다음과 같다.

 

  $$
  \frac{y_{i+1} - 2y_i + y_{i-1}}{\Delta x^2} + \frac{A_i (y_{i+1} - y_{i-1})}{2\Delta x} + B_i y_i + \theta(\Delta x^2) = f_i
  $$

 

위 식은 $y_i$에 대한 algebraic equation이 된다.
  $i = 1, \ldots, N-1$ → 총 $(N-1)$개의 방정식을 matrix로 정리할 수 있다.

$$
y_{i-1} \left( \frac{1}{\Delta x^2} - \frac{A_i}{2\Delta x} \right)
+ y_i \left( \frac{2}{\Delta x^2} + B_i \right)
+ y_{i+1} \left( \frac{1}{\Delta x^2} + \frac{A_i}{2\Delta x} \right)
= f_i
$$

 

이 식을 보면

$$
\alpha_i y_{i-1}+\beta_i y_i+\gamma_i y_{i+1}=f_i
$$

으로 정리할 수 있고 총 ODEs를 matrix로 표현하면 다음과 같다.

 

$$
\begin{bmatrix}
\beta_1 & \gamma_1 & 0 & \cdots & 0 \\
\alpha_2 & \beta_2 & \gamma_2 & \cdots & 0 \\
0 & \alpha_3 & \beta_3 & \gamma_3 & \cdots \\
\vdots & & \ddots & \ddots & \ddots \\
0 & \cdots & \alpha_{N-1} & \beta_{N-1}
\end{bmatrix}
\begin{bmatrix}
y_1 \\
y_2 \\
\vdots \\
y_{N-1}
\end{bmatrix}
=
\begin{bmatrix}
f_1 - \alpha_1 y_0 \\
f_2 \\
\vdots \\
f_{N-1} - \gamma_{N-1} y_N
\end{bmatrix}
$$

 

이 매트릭스는 TDMA (Tri-Diagonal Matrix Algorithm)이기 때문에 inversion이 쉽다. $\mathcal{O}(n)$

이 방법의 단점은 order of accuracy를 올리고 싶을 때 문제가 된다. Pade scheme을 적용해서도 가능하다.