[수치해석] Runge Kutta 4th order 증명

Introduction

 

설명을 시작하기 앞서 우리가 구하고자 하는 Ordinary Differential Equation(ODE)는 다음과 같습니다.

$y' = f(y,t)$

 

저번 글에서 작성했던 RK4 식은 아래와 같습니다.

 

$y^{n+1}=y^{n}+\frac{1}{6} k_{1}+\frac{1}{3}\left(k_{2}+k_{3}\right)+\frac{1}{6} k_{4}$
$k_{1}=\Delta t f\left(y^{n}, t^{n}\right)$
$k_{2}=\Delta t f\left(y^{n}+\frac{1}{2} k_{1}, t^{n}+\frac{\Delta t}{2}\right)$
$k_{3}=\Delta t f\left(y^{n}+\frac{1}{2} k_{2}, t^{n}+\frac{\Delta t}{2}\right)$
$k_{4}=\Delta t f\left(y^{n}+k_{3}, t^{n}+\Delta t\right)$

 

RK2에 대해서는 증명을 했지만, RK4가 더 많이 쓰이기 때문에 RK4에 대해서도 동일하게 증명을 하도록 하겠습니다. 

(되도록이면 이전 글을 봐주시기 바랍니다.)

 

RK2를 다시 떠올리면 다음과 같습니다.

$t^{n}$에서 $t^{n+1}$로 $\Delta t$만큼 적분을 할 때 다음 그림처럼 중간에 한 스텝을 더 추가합니다.

이 때 그 스텝의 시간은 $t^{n}+\alpha \Delta t$이고, 그 시간에 대응하는 상태값은 $y^{n}+\beta \Delta t f(y^{n},t^{n})$을 사용합니다. 시간이 $\alpha \Delta t$만큼 지났으므로 상태도 변했을 것이기 때문입니다. $\beta=0$이면 원래 시작했던 상태인 $y^{n}$을 쓰는 것이고 $\beta = 1$을 쓰면 $y^{n}$에서 구한 기울기를 가지고 예측한 다음 스텝의 상태인 $y^{n}+\Delta t f(y)$를 쓰는 것입니다.

 

$y^{n+1}=y^{n}+\gamma_{1} k_{1}+\gamma_{2} k_{2}$

$k_{1}=\Delta t f\left(y^{n}, t^{n}\right)$
$k_{2}=\Delta t f\left(y^{n}+\beta k_{1}, t^{n}+\alpha \Delta t\right)$

 

이렇게 $y=y^{n}, \; t=t^{n}$일 때 한 번 기울기를 구하고 $y=y^{n}+\beta \Delta t f(y^{n},t^{n}), \; t=t^{n}+\alpha \Delta t$일 때 한 번 기울기를 구한 다음에 이 기울기에 가중치를 곱해서 조합한 결과로 예측하는 것이 RK2입니다.

 

RK2 방법에서는 중간 스텝이 $t^{n}$과 $t^{n+1}$의 정가운데라고 생각하면 $\alpha=\beta = \frac{1}{2} \quad \gamma_{1}=\gamma_{2}=\frac{1}{2}$

 

 

마찬가지로 RK4는 기울기를 구하는 과정이 4번이 됩니다. 위에서는 시작 위치에서 한 번, 중간 위치에서 한 번 기울기를 구하고 이 기울기를 weighted sum을 했었습니다.

마찬가지로 RK4 역시 첫 번째는 현재 위치에서 기울기를 구하고($k_{1}$), 그 기울기를 바탕으로 $\Delta t$의 절반만큼 state $y$와 $t$가 이동했을 때 기울기를 또 구합니다.($k_{2}$)

이 때 구한 기울기만큼의 값을 가지고 또 중간 지점의 상태와 시간을 전제로 기울기를 구합니다. ($k_{3}$) 마지막으로 구한 기울기로 $\Delta t$만큼 시간이 지났을 때 그 상태와 시간에서 기울기를 구합니다. ($k_{4}$)

 

이렇게 구했을 때 최대한 정확도(accuracy)를 높일 수 있는 $\gamma_{1},\gamma_{2},\gamma_{3},\gamma_{4}$를 구하도록 합니다.

 

 

(주의 : 여기 증명을 다시 검토하는 과정에서 증명 과정은 똑같은데 globally 4th order, locally 5th order가 나오도록 해야하는데 증명에서 부족한 부분이 있어서 거기까지 나오지 않았습니다. 나중에 다시 수정하도록 하겠습니다.

 

이후 수정 : 항이 너무 많아서 MATLAB으로 돌려가면서 해보았는데 아래의 계수로 선정했을 때 $\Delta t^4$항까지 사라지면서 locally $\Delta t^5$이 성립하였습니다. 수식을 보충하고 싶으나 너무 길어져서.. 나중에야 가능할 것 같습니다.)

---

증명과정

 

 

\begin{align} y^{n+1}&=y^{n}+\gamma_{1} k_{1}+\gamma_{2} k_{2}+ \gamma_{3} k_{3}+\gamma_{4} k_{4}&\cdots(1)\\
k_{1}&=\Delta t f\left(y^{n}, t^{n}\right)&\cdots (2)\\
k_{2}&=\Delta t f\left(y^{n}+\frac{1}{2} k_{1}, t^{n}+\frac{\Delta t}{2}\right)& \cdots (3)\\
k_{3}&=\Delta t f\left(y^{n}+\frac{1}{2} k_{2}, t^{n}+\frac{\Delta t}{2}\right)& \cdots (4)\\
k_{4}&=\Delta t f\left(y^{n}+k_{3}, t^{n}+\Delta t\right)&\cdots (5)\end{align}

 

 

(1)번 식의 목적은 $k_{1,2,3,4}$의 조합으로 최대한 accuracy를 높이는 것, 즉 $\theta(\Delta t^{n})$에서 n을 높이는 것입니다.

 

원래 taylor expansion으로 얻을 수 있는 $y^{n+1}$(exact solution)은 다음과 같습니다.

 

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

 

여기서

$y' = f(y,t)=\frac{dy}{dt}$입니다.

 

또한 f는 y와 t에 대한 함수이고, y 역시 t에 dependent하기 때문에 미분할 때 더 복잡해집니다. 만약 y와 t에 independent했으면 여러 번 미분하는 것이 비교적 덜 복잡할 것입니다.

 

위 식의 $y^{\prime\prime}$과 $y^{\prime \prime \prime}$을 구해보겠습니다. 이렇게 taylor series expansion을 자세히 구하려는 이유는 나중에 위의 (1)번 식이 taylor expansion에서 나온 식(exact solution)과의 차이가 적도록 $\gamma$를 선정하기 위해서입니다.

 

$$y^{\prime \prime}=\frac{d f\left(y, t\right)}{d t}=\frac{\partial f}{\partial y} \frac{d y}{d t}+\frac{\partial f}{\partial t}=f_{y} y^{\prime}+f_{t}=f_{y} f+f_{t} \quad \cdots (6)$$

$$\begin{aligned} y^{\prime \prime \prime}=\frac{d\left(f_{y} f+f_{t}\right)}{d t} &=\left(f_{y} f+f_{t}\right) f_{y}+f\left(f_{y y} f+f_{y t}\right)+f_{y t} f+f_{t t} \\ &=f_{y}^{2} f+f_{y} f_{t}+f_{y y} f^{2}+2 f f_{y t}+f_{t t} \quad \cdots (7)\end{aligned}$$

 

위 식을 쓸 때 보충설명을 한 식은 아래에 있습니다.

 

$$\begin{aligned}
\frac{d f_{y} f}{d t} &=\frac{d f}{d t} f_{y}+f \frac{d f_{y}}{d t} \\
&=\left(f_{y} f+f_{t}\right) f_{y}+f\left(f_{yy} f+f_{y t}\right)
\end{aligned}$$

$$\frac{d f_{t}}{d t} =\frac{\partial f_{t}}{\partial y} \frac{d y}{d t}+\frac{\partial f_{t}}{\partial t}$$

 

이제 (6)번식과 (7)번식을 염두에 두고 있어야 합니다.

 

다시 RK4 식으로 돌아오겠습니다.

각 $k_{1},k_{2},k_{3},k_{4}$가 어떻게 전개되는지 해보겠습니다.

 

여기서 시작하기 앞서 미리 알아야하는게 있는데 변수가 2개인 함수가 있을 때 개별 변수의 변화에 대해 다음과 같이 나타낼 수 있습니다.

 

$$f(y+\Delta y, t+\Delta t)=f(y,t) + \Delta y \frac{\partial f}{\partial y} + \Delta t \frac{\partial f}{\partial t} + \frac{1}{2}\left( (\Delta y)^{2}\frac{\partial^{2}f}{\partial y^{2}}+2\Delta y \Delta t \frac{\partial^{2} f}{\partial y \partial t} + (\Delta t)^{2}\frac{\partial^2 f}{\partial t^2} \right)+\cdots$$

 

---

참고 : 여러 개의 변수일 때 일반화된 방법은 아래와 같이 위키피디아 문서에서 찾을 수 있었다.

 

---

이 식을 기억해두고 $k_{1},k_{2},k_{3},k_{4}$를 구하겠습니다.

또한 앞에서도 이야기했던 것처럼 $f^{n}$은 $f(y^{n},t^{n})$을 의미합니다.

 

$k_{1}=\Delta t f\left(y^{n}, t^{n}\right) $

 

$$\begin{align} k_{2}&=\Delta t f\left(y^{n}+\frac{1}{2} k_{1}, t^{n}+\frac{\Delta t}{2}\right)\\ &=\Delta t f\left(y^{n}, t^{n}\right)+\Delta t \frac{1}{2} k_{1} f_{y}^{n}+\Delta t\left(\frac{1}{2} \Delta t\right) f_{t}^{n}+\Delta t \frac{1}{2}\left(\frac{k_{1}^{2}}{4} \frac{\partial^{2} f}{\partial y^{2}}+2 \frac{\Delta t}{2} \frac{k_{1}}{2} \frac{\partial f}{\partial y \partial t}+\frac{\Delta t^{2}}{4} \frac{\partial^{2} f}{\partial t^{2}}\right) + \cdots \\ &=\Delta t f\left(y^{n}, t^{n}\right)+\frac{\Delta t^{2}}{2} f^{n} f_{y}^{n}+\frac{\Delta t^{2}}{2} f_{t}^{n}+\frac{\Delta t^{3}}{8} F_{3}+\theta\left(\Delta t^{4}\right) \end{align}$$

 

여기서 $F_{3}$은 제가 임의로 정한 문자입니다.

구체적인 식은 다음과 같습니다.

 

$\begin{align}&\frac{1}{2}\left(\frac{k_{1}^{2}}{4} \frac{\partial^{2} f}{\partial y^{2}}+2 \frac{\Delta t k_{1}}{4} \frac{\partial f}{\partial y \partial t}+\frac{\Delta t^{2}}{4} \frac{\partial^{2} f}{\partial t^{2}}\right)\\
&\approx \frac{1}{8}\left(\Delta t^{2} f^{2} f_{y y}+2 \Delta t^{2} f f_{y t}+\Delta t^{2} f_{t t}\right) (\because k_{1}=\Delta t f)\\
&\approx \frac{\Delta t^{2}}{8}\left(f^{2} f_{y y}+2 f f_{y t}+f_{t t}\right)\end{align}$

$F_{3}:= f^{2} f_{y y}+2 f f_{y t}+f_{t t}$

 

이제 $k_{3}$를 구하겠습니다.

 

$\begin{aligned} k_{3}&=\Delta t f\left(y^{n}+\frac{1}{2} k_{2}, t^{n}+\frac{\Delta t}{2}\right)\\ &=\Delta t f\left(y^{n}, t^{n}\right)+\Delta t \frac{1}{2} k_{2} f_{y}^{n}+\Delta t\left(\frac{1}{2} \Delta t\right) f_{t}^{n}+\Delta t \frac{1}{2}\left(\frac{k_{2}^{2}}{4} \frac{\partial^{2} f}{\partial y^{2}}+2 \frac{\Delta t}{2} \frac{k_{2}}{2} \frac{\partial f}{\partial y \partial t}+\frac{\Delta t^{2}}{4} \frac{\partial^{2} f}{\partial t^{2}}\right)\cdots+\theta\left(\Delta t^{4}\right) \\ &=\Delta t f^{n}+\frac{\Delta t^{2}}{2} f^{n} f_{y}^{n}+\frac{\Delta t^{3}}{4} f^{n}\left(f_{y}^{n}\right)^{2}+\frac{\Delta t^{3}}{4} f_{t}^{n} f_{y}^{n}+\frac{\Delta t^{2}}{2} f_{t}^{n} +\frac{\Delta t^{3}}{8} F_{3}+\theta\left(\Delta t^{4}\right) \end{aligned}$

 

이 식을 보면 $k_{2}$를 직접 대입하되 $\Delta t^{4}$ 이상의 항은 $\theta(\Delta t^{4})$에 넣어서 표시하지 않았습니다. 또한 여기서도 위와 같이 식을 간편하게 만들기 위해 $F_{3}$으로 표현했습니다.

 

 

$\frac{1}{8}\left(k_{2}^{2} f_{y y}+2 \Delta t k_{2} f_{y t}+\Delta t^{2} f_{t t}\right)$

$\approx \frac{\Delta t^{2}}{8}\left(f^{2} f_{yy}+2 f_{yt}+f_{t t}\right)\; (\because k_{3}\approx \Delta t f)$

$\approx \frac{\Delta t^{2}}{8}F_{3}$

 

마지막으로 $k_{4}$를 구합니다.

 

 

$\begin{aligned} k_{4}&=\Delta t f\left(y^{n}+k_{3}, t^{n}+\Delta t\right)\\ &=\Delta t f\left(y^{n}, t^{n}\right)+\Delta t k_{3} f_{y}^{n}+\Delta t(\Delta t) f_{t}^{n}+\frac{\Delta t^{3}}{2}F_{3}+\cdots+\theta\left(\Delta t^{4}\right) \\ &= \Delta t f^{n}+\Delta t^{2} f^{n} f_{y}^{n}+\frac{\Delta t^{3}}{2}\left(f^{n}\left(f_{y}^{n}\right)^{2}+f_{t}^{n} f_{y}^{n}\right)+\Delta t^{2} f_{t}^{n} +\frac{\Delta t^{3}}{2} F_{3}+\theta\left(\Delta t^{4}\right) \end{aligned}$

 

$k_{4}$에서는 $\frac{\Delta t^{3}}{8}F_{3}$ 대신 $\frac{\Delta t^{3}}{2}F_{3}$를 씁니다.

 

$\frac{1}{2}\left(k_{3}^{2} f_{4 y}+2 \Delta t k_{3} f_{y t}+\Delta t^{2} f_{t t}\right)$
$\approx \frac{\Delta t^{2}}{2}\left(f_{4 y}+2 f_{4 t}+f_{t t}\right)\; (\because k_{3}\approx \Delta t f)$

 

이렇게 구한 $k$를 가중치 $\gamma$를 곱해서 정리합니다.

 

$\begin{aligned}y^{n+1}&=\gamma_{1} k_{1}+\gamma_{2} k_{2}+\gamma_{3} k_{3}+\gamma_{4} k_{4}\\&=\left(\gamma_{1}+\gamma_{2}+\gamma_{3}+\gamma_{4}\right) \Delta t f^{n}+\left(\gamma_{2}+\gamma_{3}+2 \gamma_{4}\right) \frac{\Delta t^{2}}{2}\left(f^{n} f_{y}^{n}+f_{t}^{n}\right)+\\&\Delta t^3\left[\frac{\gamma_{2}}{8} F_{3}+\frac{\gamma_{3}}{8} F_{3}+\frac{\gamma_{4}}{2} F_{3}+ \left(\frac{\gamma_{3}}{4}+\frac{\gamma_{4}}{2}\right)\left\{ f^{n}\left(f_{n}^{n}\right)^{2}+f_{t}^{n} f_{y}^{n} \right\} \right]+\theta\left(\Delta t^{4}\right)\end{aligned}$

 

이 식을 보고 다시 위에서 exact solution을 구한 것과 비교합니다.

 

$\text{Exact solution : }$

$\begin{aligned}y^{n+1}&=y^{n}+y^{\prime} \Delta t+y^{\prime \prime} \frac{\Delta t^{2}}{2 !}+y^{\prime \prime \prime} \frac{\Delta t^{3}}{3 !}+\cdots \\&=y^{n}+f^{n} \Delta t +\left(f_{y}^{n} f^{n}+f_{t}^{n}\right)\frac{\Delta t^{2}}{2 !} + \left((f_{y}^{n})^{2} f^{n}+f_{y}^{n} f_{t}^{n}+f_{y y}^{n} (f^{n})^{2}+2 f^{n} f_{y t}^{n}+f_{t t}^{n} \right)\frac{\Delta t^{3}}{3 !}\cdots\end{aligned}$

 

이 식과 위 식에서 $\Delta t^{3}$항까지 같게 만들려면 다음 식을 만족해야합니다.

 

또한 주목할 점은 exact solution의 $(f_{y}^{n})^{2} f^{n}+f_{y}^{n} f_{t}^{n}+f_{y y}^{n} (f^{n})^{2}+2 f^{n} f_{y t}^{n}+f_{t t}^{n}$와

$F_{3}+f^{n}\left(f_{n}^{n}\right)^{2}+f_{t}^{n} f_{y}^{n}$이 같은 식이라는 것입니다.

 

따라서 $\gamma$에 대한 조건은 아래와 같습니다. (각 항을 일대일 대응해서 같도록 만들면 됩니다.)

 

$\left\{\begin{array}{l}\gamma_{1}+\gamma_{2}+\gamma_{3}+\gamma_{4}=1 \\ \gamma_{2}+\gamma_{3}+2 \gamma_{4}=1 \\ \frac{\gamma_{3}}{4}+\frac{\gamma_{4}}{2}=\frac{1}{6} \\ \frac{1}{8}\left(\gamma_{1}+\gamma_{3}+4 \gamma_{4}\right)=\frac{1}{6}\end{array}\right.$

 

각각을 풀면

 

$\gamma_{1}=\gamma_{4}=\frac{1}{6}, \quad \gamma_{2}=\gamma_{3}=\frac{2}{6}$

 

다음과 같은 일반적인 RK4가 나옵니다.

또한 이렇게 $\gamma$를 설정할 시에는 accuracy가 globally 4th order가 나온다고 합니다.

 

 

 

 

---

Conclusion

 

이 글에서는 일반적인 RK4를 RK2를 증명한 방법과 유사하게 증명하였습니다.

 

그러나 위의 방법은 일반적으로 사용되는 Runge Kutta 4th order 방법임을 알립니다. 원래는 더 일반화해서 각 계수를 구하므로, 무수히 많은 RK4 방법이 존재합니다.

 

Explicit RK4 방법의 일반적인 공식은 다음과 같습니다.

 

$y_{n+1}=y_{n}+h \sum_{i=1}^{N} \gamma_{i} k_{i}$

 

$\begin{aligned} k_{1} &=f\left(t_{n}, y_{n}\right) \\ k_{2} &=f\left(t_{n}+c_{2} h, y_{n}+h\left(a_{21} k_{1}\right)\right) \\ k_{3} &=f\left(t_{n}+c_{3} h, y_{n}+h\left(a_{31} k_{1}+a_{32} k_{2}\right)\right) \\ & \cdots \\ k_{N} &=f\left(t_{n}+c_{N} h, y_{n}+h\left(a_{N,1} k_{1}+a_{N,2} k_{2}+\cdots+a_{N, N-1} k_{N-1}\right)\right) \end{aligned}$

 

$\gamma_{i}$는 가중치이고, $k_{i}$는 위에서 구한 $k_{1},k_{2}...$와 동일한 의미입니다.

다만 $k_{n}$일 때 $f$ 안에 $k_{n-1}$일 뿐 아니라 $k_{n-2},...$ 등 이전에 구한 $k$을 조합해서 쓸 수 있습니다.

따라서 엄청나게 많은 계수에 자유도가 존재합니다. 그래서 RK4에는 unique한 방법이 있는게 아니라 무수히 많은 것입니다.

 

자세한 내용은 위키피디아를 참고해주세요.

https://ko.wikipedia.org/wiki/%EB%A3%BD%EA%B2%8C-%EC%BF%A0%ED%83%80_%EB%B0%A9%EB%B2%95#%EC%9D%BC%EB%B0%98%EC%A0%81%EC%9D%B8_4%EC%B0%A8_%EB%A3%BD%EA%B2%8C-%EC%BF%A0%ED%83%80_%EB%B0%A9%EB%B2%95

룽게-쿠타 방법 - 위키백과, 우리 모두의 백과사전