logo

룽게-쿠타 메소드에서 계수 결정하는 방법 📂수치해석

룽게-쿠타 메소드에서 계수 결정하는 방법

개요1

다음과 같이 주어진 상미분 방정식을 생각하자. yytt에 대한 함수이며, 프라임()(^{\prime})tt에 대한 미분을 의미한다.

y=f(t,y),tt0,y(t0)=y0,tn+1=tn+h y^{\prime} = f(t,y),\quad t \ge t_{0},\quad y(t_{0}) = y_{0}, \quad t_{n+1} = t_{n} + h

명시적 룽게-쿠타 메소드란, 주어진 yn=y(tn)y_{n} = y(t_{n})에 대해서 yn+1=y(tn+1)y_{n+1} = y(t_{n+1})을 다음과 같이 근사하는 방법이다.

yn+1=yn+hj=1νbjf(tn+cjh,ξj),n=0,1,2, y_{n+1} = y_{n} + h\sum_{j=1}^{\nu} b_{j} f(t_{n} + c_{j}h, \xi_{j}),\quad n=0,1,2,\cdots

여기서 각각의 ξj\xi_{j}는,

ξ1=ynξ2=yn+ha2,1f(tn,ξ1)ξ3=yn+ha3,1f(tn,ξ1)+ha3,2f(tn+c2h,ξ2)ξν=yn+hi=1ν1aν,if(tn+cih,ξi) \begin{align*} \xi_{1} &= y_{n} \\ \xi_{2} &= y_{n} + ha_{2,1}f(t_{n}, \xi_{1}) \\ \xi_{3} &= y_{n} + ha_{3,1}f(t_{n}, \xi_{1}) + ha_{3,2}f(t_{n} + c_{2}h, \xi_{2}) \\ \vdots & \\ \xi_{\nu} &= y_{n} + h\sum_{i=1}^{\nu-1}a_{\nu, i}f(t_{n} + c_{i}h, \xi_{i}) \end{align*}

RK-메소드에서 아래의 계수들은 선택하여 사용하는 것이다. 적절한 계수를 선택하는 방법에 대해서 알아본다. A=[aj,i]=[a2,1a3,1a3,2aν,1aν,2aν,ν1]andb=[b1b2bν]andc=[c1c2cν] A = [a_{j,i}] = \begin{bmatrix} \\ a_{2,1} & \\ a_{3,1} & a_{3,2} \\ \vdots & \\ a_{\nu,1} & a_{\nu,2} & \cdots & a_{\nu,\nu-1} & \end{bmatrix} \quad \text{and} \quad \mathbf{b} = \begin{bmatrix} b_{1} \\ b_{2} \\ \vdots \\ b_{\nu} \end{bmatrix} \quad \text{and} \quad \mathbf{c} = \begin{bmatrix} c_{1} \\ c_{2} \\ \vdots \\ c_{\nu} \end{bmatrix}

ν=2\nu = 2인 경우

우선 가장 간단한 ν=2\nu = 2인 경우에 대해서 살펴보자. 그러면 ξ1=yn\xi_{1} = y_{n}이고 ξ2=yn+ha2,1f(tn,yn)\xi_{2} = y_{n} + ha_{2,1}f(t_{n}, y_{n})이다. 따라서 yn+1y_{n+1}은 다음과 같다.

yn+1=yn+h[b1f(tn+c1h,ξ1)+b2f(tn+c2h,ξ2)]=yn+h[b1f(tn,yn)+b2f(tn+c2h,yn+ha2,1f(tn,yn))] \begin{align} y_{n+1} &= y_{n} + h \Big[ b_{1}f(t_{n} + c_{1}h, \xi_{1}) + b_{2}f(t_{n} + c_{2}h, \xi_{2}) \Big] \nonumber \\ &= y_{n} + h \Big[ b_{1}f(t_{n}, y_{n}) + b_{2}f(t_{n} + c_{2}h, y_{n} + ha_{2,1}f(t_{n}, y_{n}) )\Big] \end{align}

여기서 c1=0c_{1} = 0이다. 벡터값 함수 ff가 충분히 스무스하다고 가정하자. 그러면 f(tn+c2h,yn+ha2,1f(tn,yn))f(t_{n} + c_{2}h, y_{n} + ha_{2,1}f(t_{n}, y_{n}) )의 점 (tn,yn)(t_{n}, y_{n})에서 테일러 전개는 다음과 같다.

f(tn+c2h,yn+ha2,1f(tn,yn))=f(tn,yn)+(c2h,ha2,1f(tn,yn))Tf(tn,yn)+O(h2)=f(tn,yn)+h[c2f(tn,yn)t+a2,1f(tn,yn)yf(tn,yn)]+O(h2) \begin{align*} & f(t_{n} + c_{2}h, y_{n} + ha_{2,1}f(t_{n}, y_{n}) ) \\ &= f(t_{n}, y_{n}) + (c_{2}h, ha_{2,1}f(t_{n}, y_{n}))^{T} \nabla f(t_{n}, y_{n}) + \mathcal{O}(h^{2}) \\ &= f(t_{n}, y_{n}) + h\left[c_{2} \dfrac{\partial f(t_{n}, y_{n})}{\partial t} + a_{2,1}\dfrac{\partial f(t_{n}, y_{n})}{\partial y}f(t_{n}, y_{n}) \right] + \mathcal{O}(h^{2}) \\ \end{align*}

f\nabla fff그래디언트를 의미한다. 이를 (1)(1)에 대입하면,

yn+1=yn+h[b1f(tn,yn)+b2f(tn,yn)+hb2c2f(tn,yn)t+hb2a2,1f(tn,yn)yf(tn,yn)+O(h2)]=yn+h(b1+b2)f(tn,yn)+h2b2[c2f(tn,yn)t+a2,1f(tn,yn)yf(tn,yn)]+O(h3) \begin{align*} y_{n+1} &= y_{n} + h \bigg[ b_{1}f(t_{n}, y_{n}) + b_{2}f(t_{n}, y_{n}) \\ & \qquad + hb_{2}c_{2} \dfrac{\partial f(t_{n}, y_{n})}{\partial t} + hb_{2}a_{2,1}\dfrac{\partial f(t_{n}, y_{n})}{\partial y}f(t_{n}, y_{n}) + \mathcal{O}(h^{2}) \bigg] \\ &= y_{n} + h (b_{1} + b_{2})f(t_{n}, y_{n}) \\ & \qquad + h^{2}b_{2} \left[ c_{2} \dfrac{\partial f(t_{n}, y_{n})}{\partial t} + a_{2,1}\dfrac{\partial f(t_{n}, y_{n})}{\partial y}f(t_{n}, y_{n}) \right] + \mathcal{O}(h^{3}) \tag{2} \end{align*}

그런데 y=f(t,y)y^{\prime} = f(t, y)전미분d(y)=f(t,y)tdt+f(t,y)ydyd (y^{\prime}) = \dfrac{\partial f(t, y)}{\partial t}dt + \dfrac{\partial f(t, y)}{\partial y}dy이고, =ddt^{\prime} = \dfrac{d}{dt}이므로, yy^{\prime \prime}는 다음과 같다.

y=f(t,y)t+f(t,y)ydydt=f(t,y)t+f(t,y)yf(t,y) y^{\prime \prime} = \dfrac{\partial f(t, y)}{\partial t} + \dfrac{\partial f(t, y)}{\partial y}\dfrac{dy}{dt} = \dfrac{\partial f(t, y)}{\partial t} + \dfrac{\partial f(t, y)}{\partial y}f(t, y)

테일러 정리에 의해 다음을 얻는다.

yn+1=yn+hy+h22y+O(h3)=yn+hf(tn,yn)+h22(f(t,y)t+f(t,y)yf(t,y))+O(h3) \begin{align*} y_{n+1} &= y_{n} + hy^{\prime} + \dfrac{h^{2}}{2}y^{\prime \prime} + \mathcal{O}(h^{3}) \\ &= y_{n} + hf(t_{n},y_{n}) + \dfrac{h^{2}}{2}\left( \dfrac{\partial f(t, y)}{\partial t} + \dfrac{\partial f(t, y)}{\partial y}f(t, y) \right) + \mathcal{O}(h^{3}) \tag{3} \\ \end{align*}

(2)(2)(3)(3)을 비교하면 다음의 관계식을 얻는다.

b1+b2=1,b2c2=12,a2,1=c2 b_{1} + b_{2} = 1, \qquad b_{2}c_{2} = \dfrac{1}{2}, \qquad a_{2,1} = c_{2}

이를 만족하는 RK 테이블로는 아래의 것들이 있다.

cAbt0121201,02323121434,011121212, \begin{array}{c|c} \mathbf{c} & A \\ \hline & \mathbf{b}^{t} \end{array} \qquad \begin{array}{c|c} 0 \\ \frac{1}{2} & \frac{1}{2} \\[0.5em] \hline & 0 & 1 \end{array}, \qquad \begin{array}{c|c} 0 \\ \frac{2}{3} & \frac{2}{3} \\[0.5em] \hline \phantom{\dfrac{1}{2}} & \frac{1}{4} & \frac{3}{4} \end{array}, \qquad \begin{array}{c|c} 0 \\ 1 & 1 \\[0.5em] \hline \phantom{\dfrac{1}{2}} & \frac{1}{2} & \frac{1}{2} \end{array}, \qquad

ν=3\nu = 3인 경우

우선 테일러 정리에 의해 다음을 얻는다.

yn+1=yn+hy+h22y+h36y+O(h4)=yn+hfn+h22(ft+fyf)+h36(2ft2+2fy2f2+ftfy)+O(h4)(4) \begin{aligned} y_{n+1} &= y_{n} + hy^{\prime} + \dfrac{h^{2}}{2}y^{\prime \prime} + \dfrac{h^{3}}{6}y^{\prime \prime \prime} + \mathcal{O}(h^{4}) \\ &= y_{n} + hf_{n} + \dfrac{h^{2}}{2}\left( \dfrac{\partial f}{\partial t} + \dfrac{\partial f}{\partial y}f \right) + \dfrac{h^{3}}{6}\left( \dfrac{\partial^{2} f}{\partial t^{2}} + \dfrac{\partial^{2} f}{\partial y^{2}}f^{2} + \dfrac{\partial f}{\partial t}\dfrac{\partial f}{\partial y} \right) + \mathcal{O}(h^{4}) \tag{4} \end{aligned}

그리고 i=1j1aj,i=cj\sum\limits_{i=1}^{j-1} a_{j,i} = c_{j}라고 가정하자. 그러면 ν=3\nu = 3인 경우에는,

c1=0,c2=a2,1,c3=a3,1+a3,2 c_{1} = 0,\qquad c_{2} = a_{2,1},\qquad c_{3} = a_{3,1} + a_{3,2}

fn=f(tn,yn)f_{n} = f(t_{n}, y_{n})이라 표기하면,

ξ1=yn,ξ2=yn+hc2fn,ξ3=yn+ha3,1fn+ha3,2f(tn+c2h,yn+hc2fn) \xi_{1} = y_{n}, \quad \xi_{2} = y_{n} + hc_{2} f_{n}, \\[1em] \xi_{3} = y_{n} + ha_{3,1}f_{n} + ha_{3,2}f(t_{n} + c_{2}h, y_{n} + hc_{2} f_{n})

따라서 yn+1y_{n+1}은 다음과 같다.

yn+1=yn+hb1fn+hb2f(tn+c2h,ξ2)+hb3f(tn+c3h,ξ3)=yn+hb1fn+hb2f(tn+c2h,yn+hc2fn)+hb3f(tn+c3h,yn+ha3,1fn+ha3,2f(tn+c2h,yn+hc2fn)) \begin{align*} y_{n+1} &= y_{n} + hb_{1}f_{n} + hb_{2}f(t_{n} + c_{2}h, \xi_{2}) + hb_{3}f(t_{n} + c_{3}h, \xi_{3}) \\ &= y_{n} + hb_{1}f_{n} + hb_{2}f(t_{n} + c_{2}h, y_{n} + hc_{2} f_{n}) \\ & \qquad + hb_{3}f(t_{n} + c_{3}h, y_{n} + ha_{3,1}f_{n} + ha_{3,2}f(t_{n} + c_{2}h, y_{n} + hc_{2} f_{n})) \tag{5} \end{align*}

b2b_{2}가 포함된 항은 다음과 같다.

f(tn+c2h,yn+hc2f(tn,yn))=fn+(hc2)ft+(hc2fn)fy+(c2h)222ft2+(c2hfn)222fy2+O(h3) \begin{align*} & f(t_{n} + c_{2}h, y_{n} + hc_{2}f(t_{n}, y_{n}) ) \\ &= f_{n} + (hc_{2}) \dfrac{\partial f}{\partial t} + (hc_{2}f_{n})\dfrac{\partial f}{\partial y} + \dfrac{(c_{2}h)^{2}}{2} \dfrac{\partial^{2} f}{\partial t^{2}} + \dfrac{(c_{2}h f_{n})^{2}}{2} \dfrac{\partial^{2} f}{\partial y^{2}} + \mathcal{O}(h^{3}) \\ \end{align*}

이를 (5)(5)에 대입하면,

yn+1=yn+h(b1+b2)fn+h2b2c2ft+h2b2c2fyfn+h3c22b222ft2+h3c22b222fy2fn2+hb3f(tn+c3h,yn+ha3,1fn+ha3,2f(tn+c2h,yn+hc2fn))+O(h4)(6) \begin{aligned} y_{n+1} &= y_{n} + h(b_{1} + b_{2})f_{n} + h^{2}b_{2}c_{2} \dfrac{\partial f}{\partial t} + h^{2}b_{2}c_{2} \dfrac{\partial f}{\partial y}f_{n} \\ & \qquad + \dfrac{h^{3}c_{2}^{2}b_{2}}{2} \dfrac{\partial^{2} f}{\partial t^{2}} + \dfrac{h^{3}c_{2}^{2}b_{2}}{2} \dfrac{\partial^{2} f}{\partial y^{2}}f_{n}^{2} \\ & \qquad + hb_{3}f(t_{n} + c_{3}h, y_{n} + ha_{3,1}f_{n} + ha_{3,2}f(t_{n} + c_{2}h, y_{n} + hc_{2} f_{n})) + \mathcal{O}(h^{4}) \tag{6} \end{aligned}

이제 b3b_{3}가 포함된 항을 마저 계산해보면, (계수 비교에 필요하지 않은 항은 생략한다)

f(tn+c3h,yn+ha3,1fn+ha3,2f(tn+c2h,yn+hc2fn))=fn+c3hft+(ha3,1fn+ha3,2f(tn+c2h,yn+hc2fn))fy+h2c3222ft2+(ha3,1fn+ha3,2f(tn+c2h,yn+hc2fn))222fy2+O(h3)=fn+c3hft+ha3,1fnfy+ha3,2[fn+hc2ft+O(h2)]fy+h2c3222ft2+(ha3,1fn+ha3,2fn)222fy2+O(h3)=fn+c3hft+hc3fnfy+h2c2a3,2ftfy+h2c3222ft2+(hc3fn)222fy2+O(h3) \begin{align*} & f(t_{n} + c_{3}h, y_{n} + ha_{3,1}f_{n} + ha_{3,2}f(t_{n} + c_{2}h, y_{n} + hc_{2} f_{n})) \\ &= f_{n} + c_{3}h \dfrac{\partial f}{\partial t} + \left( ha_{3,1}f_{n} + ha_{3,2}f(t_{n} + c_{2}h, y_{n} + hc_{2} f_{n}) \right) \dfrac{\partial f}{\partial y} \\ & \qquad + \dfrac{h^{2}c_{3}^{2}}{2} \dfrac{\partial^{2} f}{\partial t^{2}} + \dfrac{\left( ha_{3,1}f_{n} + ha_{3,2}f(t_{n} + c_{2}h, y_{n} + hc_{2} f_{n}) \right)^{2}}{2} \dfrac{\partial^{2} f}{\partial y^{2}} + \mathcal{O}(h^{3}) \\ &= f_{n} + c_{3}h \dfrac{\partial f}{\partial t} + ha_{3,1}f_{n}\dfrac{\partial f}{\partial y} + ha_{3,2} \left[ f_{n} + hc_{2}\dfrac{\partial f}{\partial t} + \mathcal{O}(h^{2}) \right] \dfrac{\partial f}{\partial y} \\ & \qquad + \dfrac{h^{2}c_{3}^{2}}{2} \dfrac{\partial^{2} f}{\partial t^{2}} + \dfrac{(ha_{3,1}f_{n} + ha_{3,2}f_{n})^{2}}{2} \dfrac{\partial^{2} f}{\partial y^{2}} + \mathcal{O}(h^{3}) \\ &= f_{n} + c_{3}h\dfrac{\partial f}{\partial t} + hc_{3}f_{n}\dfrac{\partial f}{\partial y} + h^{2}c_{2}a_{3,2}\dfrac{\partial f}{\partial t} \dfrac{\partial f}{\partial y} \\ & \qquad + \dfrac{h^{2}c_{3}^{2}}{2} \dfrac{\partial^{2} f}{\partial t^{2}} + \dfrac{(hc_{3}f_{n})^{2}}{2} \dfrac{\partial^{2} f}{\partial y^{2}} + \mathcal{O}(h^{3}) \\ \end{align*}

이를 (6)(6)에 대입하고 정리하면 다음과 같다.

yn+1=yn+h(b1+b2)fn+h2b2c2ft+h2b2c2fyfn+h3c22b222ft2+h3c22b222fy2fn2+hb3[fn+c3hft+hc3fnfy+h2c2a3,2ftfy+h2c3222ft2+(hc3fn)222fy2]+O(h4)=yn+h(b1+b2+b3)fn+h2(b2c2+b3c3)ft+h2(b2c2+b3c3)fyfn+h3(b2c22+b3c32)22ft2+h3b3c2a3,2ftfy+h3(b2c22+b3c32)2fn22fy2+O(h4) \begin{align*} y_{n+1} &= y_{n} + h(b_{1} + b_{2})f_{n} + h^{2}b_{2}c_{2} \dfrac{\partial f}{\partial t} + h^{2}b_{2}c_{2} \dfrac{\partial f}{\partial y}f_{n} \\ & \qquad + \dfrac{h^{3}c_{2}^{2}b_{2}}{2} \dfrac{\partial^{2} f}{\partial t^{2}} + \dfrac{h^{3}c_{2}^{2}b_{2}}{2} \dfrac{\partial^{2} f}{\partial y^{2}}f_{n}^{2} \\ & \qquad + hb_{3} \left[ f_{n} + c_{3}h\dfrac{\partial f}{\partial t} + hc_{3}f_{n}\dfrac{\partial f}{\partial y} + h^{2}c_{2}a_{3,2}\dfrac{\partial f}{\partial t} \dfrac{\partial f}{\partial y} \right. \\ & \qquad + \left. \dfrac{h^{2}c_{3}^{2}}{2} \dfrac{\partial^{2} f}{\partial t^{2}} + \dfrac{(hc_{3}f_{n})^{2}}{2} \dfrac{\partial^{2} f}{\partial y^{2}} \right] + \mathcal{O}(h^{4}) \\ &= y_{n} + h(b_{1} + b_{2} + b_{3})f_{n} + h^{2}(b_{2}c_{2} + b_{3}c_{3}) \dfrac{\partial f}{\partial t} + h^{2}(b_{2}c_{2} + b_{3}c_{3}) \dfrac{\partial f}{\partial y}f_{n} \\ & \qquad + \dfrac{h^{3}(b_{2}c_{2}^{2} + b_{3}c_{3}^{2})}{2} \dfrac{\partial^{2} f}{\partial t^{2}} + h^{3}b_{3}c_{2}a_{3,2}\dfrac{\partial f}{\partial t} \dfrac{\partial f}{\partial y} + \dfrac{h^{3}(b_{2}c_{2}^{2} + b_{3}c_{3}^{2})}{2} f_{n}^{2}\dfrac{\partial^{2} f}{\partial y^{2}} + \mathcal{O}(h^{4}) \end{align*}

위 식과 (4)(4)의 계수를 비교하면, 다음의 관계식을 얻는다.

b1+b2+b3=1,b2c2+b3c3=12,b2c22+b3c32=13,b3c2a3,2=16 b_{1} + b_{2} + b_{3} = 1, \quad b_{2}c_{2} + b_{3}c_{3} = \dfrac{1}{2}, \quad b_{2}c_{2}^{2} + b_{3}c_{3}^{2} = \dfrac{1}{3}, \quad b_{3}c_{2}a_{3,2} = \dfrac{1}{6}

이를 만족하는 RK 테이블로는 아래의 것들이 있다.

0121211212162316(classical RK method) \begin{array}{c|c} 0 \\ \frac{1}{2} & \frac{1}{2} \\[0.5em] 1 & -1 & 2 \\[0.5em] \hline \phantom{\dfrac{1}{2}} & \frac{1}{6} & \frac{2}{3} & \frac{1}{6} \end{array} \quad (\text{classical RK method}) 023232302312143838(Nystrom scheme) \begin{array}{c|c} 0 \\ \frac{2}{3} & \frac{2}{3} \\[0.5em] \frac{2}{3} & 0 & \frac{2}{3} \\[0.5em] \hline \phantom{\dfrac{1}{2}} & \frac{1}{4} & \frac{3}{8} & \frac{3}{8} \end{array} \quad (\text{Nystrom scheme})

ν=4\nu = 4인 경우

마찬가지의 방법으로 계수에 대한 조건을 구할 수 있으며, 위에서 본 규칙대로라면 다음과 같음을 짐작할 수 있다.

b1+b2+b3+b4=1,b2c2+b3c3+b4c4=12,b2c22+b3c32+b4c42=13,b2c23+b3c33+b4c43=14 b_{1} + b_{2} + b_{3} + b_{4} = 1, \quad b_{2}c_{2} + b_{3}c_{3} + b_{4}c_{4} = \dfrac{1}{2}, \\[1em] b_{2}c_{2}^{2} + b_{3}c_{3}^{2} + b_{4}c_{4}^{2} = \dfrac{1}{3}, \quad b_{2}c_{2}^{3} + b_{3}c_{3}^{3} + b_{4}c_{4}^{3} = \dfrac{1}{4}

RK4에서 가장 유명한 계수는 다음과 같다. 아래의 테이블은 실제로 위의 조건을 만족한다.

012121201210010016131316 \begin{array}{c|cccc} 0 & \\[0.5em] \frac{1}{2} & \frac{1}{2} \\[0.5em] \frac{1}{2} & 0 & \frac{1}{2} \\[0.5em] 1 & 0 & 0 & 1 \\[0.5em] \hline {\displaystyle {\color{white}\dfrac{0}{0}}} & \frac{1}{6} & \frac{1}{3} & \frac{1}{3} & \frac{1}{6} \end{array}

같이보기


  1. Arieh Iserles, A First Course in the Numerical Analysis of Differential Equations (2nd Edition, 2009), p38-41 ↩︎