[수치해석] Numerical solution of ODE (1) introduction

이번에는 ODE를 수치해석적인 방법으로 푸는 것에 대해 다루겠습니다.

 

1. Introduction

 

예를 들어 $y^{''}+wy=f(x)$

$\frac{d^{2}y}{dx^{2}}+wy=f(x)$

 

적분 상수가 2개 나오므로 조건이 2개가 필요합니다.

 

initial value problem(IVP)이면 초기조건이 주어지고

$y(0)=y_{0}$

$\left. \frac{dy}{dx} \right |_{x=0}=v_{0}$

 

boundary value problem(BVP)이면

$y(0)=y_{0}$

$y(L)=y_{L}$

 

일단 처음에는 IVP를 먼저 풀고, 그 다음에 BVP로 확장하겠습니다.

 

모든 방법은 $0\leq t \leq t^{n}$에서의 solution을 알고 있다고 가정하고 그 solution을 $t^{n+1}-t^{n}+\Delta t$에서의 solution을 알아내는데 사용합니다.

 

# Notation 참고 
$x_{B}^{A}$ : $A$는 시간에서의 위치(Timesteps), $B$는 공간에서의 위치를 의미합니다. 

 

예를 들면,

$y^{'}=f(y,t)$이라고 할 때,

$y^{n+1}=y(t^{n}+\Delta t)=y^{n}+\Delta t (y^{n})^{'}+\frac{\Delta t^{2}}{2!}(y^{n})^{''}+\frac{\Delta t^{3}}{3!}(y^{n})^{'''}+\cdots$

 

이 식을 $y^{'}$에 대해서 정리합니다.

 

$(y^{n})^{'}=\frac{y^{n+1}-y^{n}}{\Delta t}-\frac{\Delta t}{2!}(y^{n})^{''}+\cdots$

$\frac{y^{n+1}-y^{n}}{\Delta t}=f(y^{n},t^{n})+\theta(\Delta t^{1})$

 

$y^{n+1}=y^{n}+\Delta t \cdot f(y^{n},t^{n})+\theta(\Delta t^{2})$

 

여기에서 $y^{n}, f(y^{n},t^{n})$은 과거의 solution에서 알 수 있는 항들이고 $y^{n+1}$은 explicit euler라고 합니다.

 

$\left\{\begin{matrix}
n=0, & y^{1}=y^{0}+\Delta f(y^{0},t^{0})\\ 
n=1, & y^{2}=y^{1}+\Delta f(y^{1},t^{1})
\end{matrix}\right.$

$t_{1}=t^{0}+\Delta t$

 

따라서 이 내용을 정리하면 explicit euler는

$y^{n+1}=y^{n}+\Delta f(y^{n},t^{n})+\theta(\Delta t^{2})$

이 방법은 locally 2nd order in time

그리고 전체에 대한 적분으로 생각하면 order가 줄어들어들기 때문에(by mean value theorem)

globally 1st-order in time입니다.

 

 

2. Explicit Euler / Implicit Euler

 

1st order 미분방정식을 푸는 방법을 배우고 고차 미분방정식을 일차 미분방정식으로 바꿔서 풀 예정입니다.

따라서 처음에는 1st order 미분방정식을 푸는 방법에 대해서 배울 것입니다.

 

위에서 explicit euler에 대해 배웠는데 이 식에 대해 더 알아보겠습니다.

 

$y^{n+1}=y^{n}+\Delta f(y^{n},t^{n})+\theta(\Delta t^{2})$

 

이 식을 알고 나면, 궁금한 것은 long-term에서 이 시스템의 solution이 어떻게 되느냐는 것입니다.

즉, n이 증가할 때 solution이 어떻게 되는지 알고자 합니다.

 

 

$\text{ex) }\frac{dy}{dt}+y^{2}=0$

$y(0)=1=y^{0}$ 다음과 같은 initial condition이 주어졌을 때,

 

$\text{by explicit euler, }y^{n+1}=y^{n}+\Delta t (-y^{n})^{2} (\because \frac{dy}{dt}=-y^{2})$

$\begin{cases}y^{1}=y^{0}-\Delta t (y^{0})^{2}=1-\Delta t & \text{when }n=0\\ y^{2}=y^{1}-\Delta t (y^{1})^{2}=(1-\Delta t)-\Delta t (1-\Delta t)^{2} & \text{when }n=1\\ \vdots \end{cases}$

 

 

이와 같이 n을 늘려가면서 계속 구하면 됩니다.

 

이렇게 explicit euler 방법이 아니더라도 다른 방법들도 존재합니다.

 

 

* Classification

 

ODE를 푸는 방법에는 explicit과 implicit 방법으로 나눌 수 있습니다.

 

$\begin{cases}\text{Explicit : }& y^{n+1}=f(y^{n},y^{n-1},y^{n-2},\ldots)\\ \text{Implicit : }& y^{n+1}=f(y^{n+1},y^{n},y^{n-1},\ldots) \end{cases}$

 

explicit 방법은 이미 알고 있는 solution으로부터 다음 solution을 알아내는 것이고, implicit은 미래의 값을 알고 있는 것처럼 식에 넣어서 풀어내는 것입니다.

이 분류가 ODE의 stability에서 차이를 줍니다. explicit보다 implicit이 더 stable합니다.

 

또한 multi-step 방법과 substep 방법이 있는데, substep 방법은 하나의 시간간격을 더 잘게 쪼개서 작은 substep을 만들어 계산을 진행하는 것입니다. Runge-kutta(룽게-쿠타) 방법이 그 예시입니다.

multi-step 방법은 $t^{n+1}$ 시점을 알아내기 위해 그 전 시간들을 사용하는 방법입니다.

 

$\begin{cases}\text{Substep methods : }& t^{n}\leq t^{n_{1}},t^{n_{2}},t^{n_{3}},\ldots \leq t^{n+1}\\ \text{Multi-step methods : }& t^{n-2},t^{n-1},t^{n}\rightarrow t^{n+1} \end{cases}$

 

 

3. issues of solving ODE

 

ODE를 풀 때 고려해야 하는 요소에 대해 이야기하겠습니다.

 

1) Accuracy

2) Cost

3) Stability

 

이 세 가지입니다.

 

정확도와 cost는 지금까지 계속해서 얘기했던 것으로, 정확도는 $\theta(\Delta t^{n})$과 같이 시간 간격을 얼마나 줄이느냐에 따라서 몇 제곱으로 줄어드는지 나타낸 것이고, cost는 알고리즘의 계산복잡도를 말합니다. 예를 들면 $O(n^{3})$

그리고 지금까지 언급하지 않은 나머지 하나는 stability입니다.

수치해석에서 말하는 stability는 solution이 bounded되느냐 아니냐로 따집니다. solution이 oscillate하더라도 그 진폭이 무한으로 가지 않는다면 stable합니다.

 

기본적으로 numerical method는 unstable한 시스템을 해석할 수 없습니다.

 

수치해석 방법을 stability 관점에서 나누면 세 가지로 분류할 수 있습니다.

 

1) A-stable method (Absolutely stable) : 항상 stable한 결과를 내는 방법입니다. 이 방법은 발산하지는 않지만 항상 좋은 것이 아닙니다. unstable한 시스템을 넣으면 그만큼 unstable하게 나와야 하는데 그렇지 않아서 stable하게 해석하면, 잘못된 해석이기 때문입니다.

 

2) Conditionally stable

 

3) A-unstable : 항상 unstable한 결과를 내놓는 방법은 잘못된 것이므로 사용하지 않습니다.

 

implicit method는 일반적으로 A-stable하고, explicit method는 일반적으로 conditionally stable합니다.

방금 전에 살펴보았던 explicit euler를 다시 가져오겠습니다.

 

$y^{n+1}=y^{n}+\Delta f(y^{n},t^{n})+\theta(\Delta t^{2})$

 

이 방법에서 우리가 바꿀 수 있는 변수는 $\Delta t$입니다.

이 $\Delta t$가 작아질 수록 error는 작아지므로 정확도가 올라갑니다. 또한 $\Delta t$가 너무 크면 결과가 unstable하며, 또 너무 작으면 계산 cost가 증가합니다.

 

따라서 적절한 $\Delta t$를 고르는 것이 중요합니다.

앞으로 다룰 것은 식에 대한 이해를 바탕으로 stable한 결과를 낼 수 있는지에 대한 내용입니다.

 

 

4. Model problem

 

많은 ODE가 존재하지만 그 모든 ODE에 어떤 method를 적용해보면서 stability를 테스트해보는 것은 불가능합니다.

따라서 대표성을 가지는 ODE를 하나 정해서 풀어보면 좋을 것입니다.

 

결론적으로는 $y^{'}=\lambda y$를 model problem으로 정해서 푸는데, 이 방정식으로 결정하게 된 과정을 설명하겠습니다.

 

$y^{'}=\frac{dy}{dt}=f(y,t)$

$\begin{align*} f(y,t)&=f(y^{0},t^{0})+(t-t^{0})\left. \frac{\partial f}{\partial t} \right|_{(y^{0},t^{0})}+\frac{(t-t_{0})^{2}}{2!}\left. \frac{\partial^2 f}{\partial t^2}\right|_{(y^{0},t^{0})}+\cdots \\ & +(y-y^{0})\left. \frac{\partial f}{\partial y} \right|_{(y^{0},t^{0})}+\frac{(y-y_{0})^{2}}{2!}\left. \frac{\partial^2 f}{\partial y^2}\right|_{(y^{0},t^{0})}+\cdots \\ & = \lambda y + \alpha_{1}+\alpha_{2}t +\theta(\Delta t^{2},\Delta y^{2})+\cdots \end{align*}$

 

taylor expansion을 통해 $\lambda y + \alpha_{1}+\alpha_{2}t$와 error term으로 바꿀 수 있습니다.

 

$\text{Homogeneous solution : }y^{'}-\lambda y=0 \Rightarrow y_{h}(t)=y^{0}e^{\lambda t}$

$\text{Inhomogeneous solution : }y^{'}-\lambda y=\alpha_{1}+\alpha_{2}t\Rightarrow y_{p}(t) = a+bt $라고 가정한 다음에 대입해서 $a,b$를 구한다.

$y^{'}-\lambda y = \alpha_{1}+\alpha_{2}t$

$b-\lambda(a+bt)=\alpha_{1}+\alpha_{2}t$

$a=-\frac{\lambda \alpha_{1}-\alpha_{2}}{\lambda^{2}}$

$b=-\frac{\alpha_{2}}{\lambda}$

 

따라서 total solution은

$y(t)=y^{0}e^{\lambda t}-\frac{(\lambda \alpha_{1}-\alpha_{2})}{\lambda^{2}}-\frac{\alpha_{2}}{\lambda}t$

 

long term에 대해서 해를 구하면 inhomogeneous solution 보다는 homogeneous solution이 더 지배적이게 됩니다

 

즉 $y(t)\approx y^{0}e^{\lambda t}\Rightarrow y^{'}=\lambda y \text{ : model problem}$

 

이렇게 해서 $y^{'}=\lambda y$를 model problem으로 정하고, 여러가지 수치해석 기법을 이 문제에 적용해보면서 수치해석 기법의 accuracy, stability, phase error를 살펴보려고 합니다.

 

여기서 $\lambda$는 앞서 본 식에서 $\left. \frac{\partial f}{\partial y} \right |_{(y^{0},t^{0})}$으로 구한 값입니다.

시스템의 eigenvalue와 같습니다.

---

다음부터는 이번 내용을 바탕으로

수치해석 기법의 특징을 분석하도록 하겠습니다.