기타 등등

다된밥통과 함께하는 Adaptive RK4

다된밥통 2020. 11. 16. 10:57

안녕하세요 다된밥통입니다! 오늘은 만약 우리가 어떤 운동의 미분방정식만을 알고 있고 이를 exact하게 풀 수 없을 때, numerical solution을 구하는 수치해석적 방법인 RK4에 대해 설명하려 하는데요! 그중에서도 RK4의 step size를 adaptive하게 바꾸는 알고리듬인 adaptive RK4에 대한 얘기입니다! (이는 한양대학교 손승우 교수님 강의 슬라이드를 기초로 하고 있음을 알려드립니다!)

 

 

우리 주변 (특히 물리학과!)에서 보는 다양한 시스템에서 우리는 운동방정식을 기술할 수 있고, 특히 고전물리적 시각으로는 이러한 시스템의 운동방정식을 풀어냄으로써 운동상태를 기술할 수 있어요!

우리가 고등교육부터 배워온 대표적인 미분방정식으로 쓰이는 시스템은 행성의 궤도운동이 있겠는데요! 뉴턴의 만유인력법칙 (가속도는 변위의 이차미분!)으로 이 궤도운동들을 기술할 수 있겠죠?!

 

물론 이러한 시스템은 여러 가정 밑에서 exact solution을 구할 수 있지만요! 오늘 우리는 이 궤도운동에 대한 numerical solution을 구하는 방법을 살펴보려 해요! 어디 한번 볼까요!?

 

 

일반적인 미분방정식에 대해서는 우리는 크게 2가지 문제로 나누어볼 수 있는데요! 하나는 초기값 문제, 다른 하나는 경계값 문제입니다! 오늘 이 포스팅에서는 초기값 문제 (initial value problem)에 대해서만 살펴보려 해요!

 

numerical solution을 구하고자하는 아이디어는! 우리는 미분방정식을 알고 있기 때문에 이를 이용하자인데요! 임의의 함수 y에 대한 1차 미분방정식은 함수 y 위에서의 접선방향 기울기로 생각할 수 있죠! 만약 우리가 exact한 form을 모른다면 이 미분방정식을 이용해서 (즉, 기울기!), 아주 작은 미소 step size (작은 걸음이라 표현할게요!)를 나아가보자라는 것이죠!

 

정리하면, 함수 y 위에서 초기값을 알 때!, 그 초기값으로부터의 기울기 방향으로 아주 작은 걸음을 가는 과정을 반복하면 이 함수의 form을 나타낼 수 있다는 것이죠!

 

 

함수 y의 미분은 사실 정의상 작은 걸음동안의 함수 y의 변화량으로 쓰이기 때문인데요! 이를 이용해서 업데이트해 나가는 방법을 Euler's method라고 합니다!

 

 

뿐만 아니라, 우리가 그 미분에서의 Taylor expansion을 이용하면 시간 t에 대한 higher order term들도 고려해볼 수 있는데요! 몇차항까지 고려하느냐에 따라 Taylor method of order n (n차항까지 고려)로 불리게 됩니다!

 

물론, 이러한 numerical solution들은 그 걸음폭에 대한 error가 생기게 되는데요! local하게는 h (걸음폭)의 n+1승만큼, global하게는 h의 n승 order로 쌓이게 됩니다! (그렇다면, h를 굉장히 작게, 걸음을 아주 작게 가져가면 error가 무시할 정도로 되겠죠?)

 

 

이러한 Taylor method들 중 order 4가 그 유명한 RK4 이구요! 자세한 유도과정은 skip하겠지만... 그 의미를 살펴보면! RK4는 초기값이 주어졌을 때!, 그때의 현재 위치의 기울기, 절반 걸음 후의 기울기, k2를 이용한 절반 걸음 후의 기울기, 그리고 k3를 이용한 한 스텝 후의 기울기를 적절한 가중치를 둔 평균 기울기로 다음 위치를 업데이트해나가겠다라고 생각하면 좋을거 같아요! 즉, 여러 기울기를 바라보고 그 평균으로 나아가겠다 라는 것이지요!

 

RK4뿐만아니라 더 higher order를 고려하면 더 정확해지지 않을까? 네! 물론 당연히 더 정확해지지만, 아래 표에서처럼 그 효율이 살짝 떨어지게 되고, 이 때문에 많은 곳에서 RK4를 자주 쓰는 것 같아요! (물론 이 모든 내용은 개인적인 제 생각이.. 다분하지만...)

 

 

다시 행성의 궤도운동으로 돌아와보면! 우리는 미분방정식으로 이 시스템을 기술할 수 있고 (뉴턴의 만유인력 법칙을 이용!), 이때 우리가 초기값을 알면, RK4와 같은 numerical method를 이용해 풀 수 있는데요! 특히, 궤도운동은 근일점에서 매우 큰 힘을 느끼게 되고 이로인해, 충분히 작은 걸음으로 업데이트하지 않으면, 아주 큰 미분값으로 인해 그 궤도가 아래처럼 휘어져버리거나 밖으로 튕겨버리게 되는 결과를 나타내게 될 수 있어요!

 

물론, 정말정말 작은 걸음으로 가게 된다면 그 문제는 해결되겠죠? 그치만, 이 근일점을 대비해 아주 작은 걸음을 하게되면 나머지 궤도부분에서 불필요하게 계산을 많이 하게 되고, 이로인해 전체 계산 속도가 오래 걸릴 수 밖에 없는데요!

 

여기서 생각한 아이디어는, 걸음걸이를 adaptive하게 해서 매 순간 적절한 걸음으로 나아갈 수 없을까!? 인 것이며, 오늘의 주제가 되는 adaptive RK4의 메인 아이디어가 되겠습니다~!

 

 

정리해서 말하면, 아래 그림처럼! 우리가 구하고자 하는 함수가 특정부분에서 굉장히 stiff하게되면, 거기에 알맞는 걸음폭을 생각해야하는데요! 그렇게 되면 다른 부분에서 불필요한 계산이 많이 발생하게 되니! 그때그때마다 적절한 걸음폭을 생각해서 나아가자가 바로 adaptive RK4가 되겠습니다!

 

이러한 방법을 어떻게 해나가자하면요! 현재 위치에서 h라는 걸음폭으로 두걸음 가본 위치와 2h로 한걸음 가본 위치를 계산한 다음, 그 두 위치의 차이가 무시할 정도로 작다면! 다음부턴 2h로 가자!이고, 아닌 경우에는 좀더 작은 걸음으로 다시 계산해서 보정을 해나가는 것이지요!

 

두배라는 효과덕에 아주 기하급수적으로 빠르게 충분히 넓은 h로도 가게 되며, 반대의 경우도 (아래 설명하겠지만) 기하급수적으로 줄어들게 되므로 빠르게 원하는 accuracy를 유지하며 계산해나갈 수 있게 되는 것이지요!

 

 

좀 더 수식적으로 살펴보면,

RK4는 local하게 h의 n+1승 즉 5승 order로 error가 쌓이게 되고, 이를 표현하면 다음 스텝에서의 함수 값은 정확한 값보다 어떤 비례상수 c를 곱해진 ch^5 만큼의 error가 쌓이게 되겠죠? 이와 같은 방법으로 두 걸음을 나아가면 최종적으로 2ch^5가 누적되게 되는 것이죠! 같은 방법으로 2h로 한걸음을 가게 되면 c(2h)^5 = 32ch^5 가 쌓이게 되고, 우리는 그 두 값이 거의 같아지길 바라니까 그 차이인 30ch^5이 무시할 수 있을 만한 겂이 되는 h로 업데이트해나가면 되는 것이죠!

 

업데이트 된 걸음폭을 h' 라고하면, c(h')5 가 h에 대해 풀어보면 아래 그림에서 보이는 것처럼 쓸 수 있고, 따라서 h' 를 rho 값 적절하게 곱해줌으로써 업데이트를 할 수 있게 됩니다!

 

 

빨간 박스 안의 rho 값을 통해, 만약 rho가 1보다 크다는 것은 분모에 있는 두 값 사이의 차이가 충분히 작다는 것이고 따라서, 더 큰 걸음폭으로 걸어가도록 업데이트를 해주는 것이죠!

반대의 경우도 만약 rho가 1보다 작다면, 분모에 위치한 두 값 사이의 차이가 유의미하게 크다는 것이 되고, 앞선 업데이트 식을 통해 h'를 다시 계산해주어 프로세스를 반복하게 됩니다!

 

 

Adaptive RK4의 장점은 이러한 두번 계산된 값을 버리지 않고 이를 적절하게 이용해서 error의 order를 한 차수 더 줄일 수 있기 때문인데요! 아래와 같이 식을 업데이트하게 되면, 최종적으로 h^6 의 local truncation error 만을 가져가게 되는 것이죠!

 

 

제 부족한 실력이지만... 만약 지구의 공전운동을 adaptive RK4를 이용해 사용해보면 어떻게 될까요!!? (물론, 물리학과 수업에서는 편의상 dimensionless parameter들로 바꿔서 보다 편한 값을 쓰게 되지만! 여기서는 그런 과정 없이 실제 값을 써볼게요! 수업시간에 몇몇 학생들은 dimensionless 과정을 잘 따라오지 못하는 경우도 있더라구요! 이 설명은 또 다음기회에...!)

 

저는 보통 코딩을 하게되면, 크게 세가지 파트로 구성하는데요! 1. 패키지 불러오기, 2. 코딩하기, 3. 그림 그리기 로 하고 있어요! (저보다 초심자분들에겐 도움이 되길...)

def adaptive_rk4 부분을 살펴보면, while 구문을 통해 rho 값이 1보다 작아질 때까지 반복하여 원하는 걸음폭을 설정하고 그로부터 다음 스텝으로 나아가는 것을 확인해볼 수 있어요! (찬찬히 한번 따라가보면 좋을거 같네요!!)

 

 

지구의 공전운동을 그려보면 아래 오른쪽 그림과 같이 표현해볼 수 있어요! 오른쪽 끝이 태양으로부터의 원일점이며 여기를 초기값으로 준 상황입니다! RK4, adaptive RK4 모두 잘 풀어낸 걸 알 수 있고!

adaptive RK4에서의 h (걸음폭)을 시간에 대해 나타내면 (c) step size 그림과 같이 나타나게 됩니다! 해당 색상의 값이 급변하게 되는 이유로 h가 업데이트 되게 되어 전체적으로 h의 값이 시간에 따라 유동적을 바뀌면서 업데이트 되는 것을 확인할 수 있겠죠!!?

 

 

물론 python에서는 scipy 패기지 안의 scipy.integrate.ode 함수가 이미 이러한 adaptive step size를 지원하고 있고... 위의 내용을 잘 모르겠다!!? 싶으면 scipy.integrate.ode 함수를 쓰면 그냥 된다... 라고 생각하시면 될거 같아요!

 

 

참고에 적혀있는 것처럼 이러한 내용은 한양대 손승우교수님의 전산물리 수업자료를 바탕으로 하였는데요! python에선 ode함수가 있더라!!라고 알려준건 현재 조교중인 박상준선생님이 당시 학생일 때, 레포트 제출 시 저에게 알려준 내용입니다! (박상준선생님의 허락을 구해 아래 첨부합니다!)

 

 

 

네! 물론 c,c++에서도 아래와 같이 써볼 수 있구요! 보시면서 따라하면 충분히 이해되시리라 믿습니다~!!

 

 

오늘은 함수의 미분값이 특정부분에서 너무 stiff해서 step size를 adaptive하게 조절하며 numerical하게 풀어나가는 알고리듬인 Adaptive RK4에 대해 포스팅해보았습니다!!

틀린 부분이라던지... 좀더 많은 조언은 아래 이메일주소나 댓글로 알려주시면 저에게 많은 도움이 될 것 같아요!!

오늘 하루도 좋은 하루 되세요~!!! 저는 다음에 또 돌아오겠습니다~!