Runge-Kutta 알고리듬, 예제 by 바죠

출처: 위키피디아

Runge-Kutta 방법은, 미분방정식들을 동시에 적분하는 수치 알고리듬이다. 아주 널리 알려진 유용한 계산 방법이다.
"
first-order differential equation set,
self-starting,
독일 사람들이 만듬.
"

In numerical analysis, the Runge–Kutta methods are an important family of implicit and explicit iterative methods for the approximation of solutions of ordinary differential equations. These techniques were developed around 1900 by the German mathematicians C. Runge and M.W. Kutta.
http://en.wikipedia.org/wiki/Runge-Kutta


http://gams.nist.gov/serve.cgi/Module/ODE/RKSUITE.F90/12230/
dderkf.txt

      SUBROUTINE DDERKF (DF, NEQ, T, Y, TOUT, INFO, RTOL, ATOL, IDID,
     +   RWORK, LRW, IWORK, LIW, RPAR, IPAR)
C***BEGIN PROLOGUE  DDERKF
C***PURPOSE  Solve an initial value problem in ordinary differential
C            equations using a Runge-Kutta-Fehlberg scheme.
C***LIBRARY   SLATEC (DEPAC)
C***CATEGORY  I1A1A
C***TYPE      DOUBLE PRECISION (DERKF-S, DDERKF-D)
C***KEYWORDS  DEPAC, INITIAL VALUE PROBLEMS, ODE,
C             ORDINARY DIFFERENTIAL EQUATIONS, RKF,
C             RUNGE-KUTTA-FEHLBERG METHODS
C***AUTHOR  Watts, H. A., (SNLA)
C           Shampine, L. F., (SNLA)
C***DESCRIPTION



 DF -- Provide a subroutine of the form
C                               DF(X,U,UPRIME,RPAR,IPAR)
C             to define the system of first order differential equations
C             which is to be solved.  For the given values of X and the
C             vector  U(*)=(U(1),U(2),...,U(NEQ)) , the subroutine must
C             evaluate the NEQ components of the system of differential
C             equations  DU/DX=DF(X,U)  and store the derivatives in the
C             array UPRIME(*), that is,  UPRIME(I) = * DU(I)/DX *  for
C             equations I=1,...,NEQ.  X와 U는 각각 변화되면 안됩니다.

cf. 실제 응용의 한 예: 분자구조 이완
http://incredible.egloos.com/3141336

!234567890
       module rk_df
       implicit none
       public
       save
       integer neq,info(15),idid,lrw,liw
       real*8 t, tout
       integer, allocatable ::  iwork(:),ipar(:)
       real*8, allocatable :: y(:), rwork(:), rtol(:), atol(:), rpar(:)

       contains

       subroutine initial
       implicit none

      t=0.0d0
       tout=0.1d0

       neq=3
       allocate(y(neq))

       y(1)=0.d0
       y(2)=1.d0
       y(3)=1.d0

      info(1)=0
       lrw= 33+7*neq ; liw= 34
       allocate(iwork(liw)) ; allocate(rwork(lrw))
       allocate(atol(neq),rtol(neq))
       rtol=1.d-10 ; atol=1.d-10
       allocate(rpar(neq)) ; allocate(ipar(neq))
       end subroutine initial

       subroutine final
       implicit none
       deallocate(y) ; deallocate(rwork) ; deallocate(iwork)
       deallocate(atol,rtol) ; deallocate(rpar) ; deallocate(ipar)
       end subroutine final

!      DF(X,U,UPRIME,RPAR,IPAR)  x, u가 이 서브루틴에서 바뀌면 안된다. uprime()만 계산해 주세요.
       subroutine df(x,u,uprime,rpar_z,ipar_z)
       implicit none
       real*8 x,u(1)
       real*8 rpar_z(1)
       integer ipar_z(1)
       real*8 uprime(1)

       uprime(1)=u(2)*u(3)
       uprime(2)=-u(1)*u(3)
       uprime(3)=-0.522d0*u(1)*u(2)

       end subroutine df

       end module rk_df
!
       program rk_test
       USE rk_df, ONLY : neq, t, y, tout, info, rtol, atol, idid,  rwork, lrw, iwork, liw, rpar, ipar
       USE rk_df, ONLY : df,initial,final
       implicit none

       call initial



       call dderkf(df, neq, t, y, tout, info, rtol, atol, idid,  rwork, lrw, iwork, liw, rpar, ipar)
!      write(6,*) t,' t'
       write(6,*) tout,' tout'
       write(6,*) y,' y'
       write(6,'(i5,1x,a4)') idid,'idid'
       write(6,'(15i5,1x,a4)') info,'info'

       tout=tout+0.1d0
       call dderkf(df, neq, t, y, tout, info, rtol, atol, idid,  rwork, lrw, iwork, liw, rpar, ipar)
!      write(6,*) t,' t'
       write(6,*) tout,' tout'
       write(6,*) y,' y'
       write(6,'(i5,1x,a4)') idid,'idid'
       write(6,'(15i5,1x,a4)') info,'info'


       tout=tout+0.1d0
       call dderkf(df, neq, t, y, tout, info, rtol, atol, idid,  rwork, lrw, iwork, liw, rpar, ipar)
!      write(6,*) t,' t'
       write(6,*) tout,' tout'
       write(6,*) y,' y'
       write(6,'(i5,1x,a4)') idid,'idid'
       write(6,'(15i5,1x,a4)') info,'info'


       call final

       end program rk_test


함수이름을 서브루틴의 한 인수로 넘겨서 계산한 적은 여러번 있었다. 서브루틴 이름을 직접 넘겨서 계산하기는 처음이다. 물론, 문법적으로 하자가 없다.
pgf90 rk_test.f90 dderk.f
a.out

   0.1000000000000000       tout
   9.9747046229319682E-002   0.9950128274411856        0.9973998069856306     
  y
    2 idid
    1    0    0    0    0    0    0    0    0    0    0    0    0    0    0 info
   0.2000000000000000       tout
   0.1979932758824953        0.9802033782404827        0.9897155661785535     
  y
    2 idid
    1    0    0    0    0    0    0    0    0    0    0    0    0    0    0 info
   0.3000000000000000       tout
   0.2933201731290916        0.9560142656098026        0.9772864626591783     
  y
    2 idid
    1    0    0    0    0    0    0    0    0    0    0    0    0    0    0 info


----------------------------------------------------------------

http://www.math.tulane.edu/~amedovik/



----------------------------------------------------------------
dv/dt= a
dx/dt= v
neq=2 (입자의 수 1개, 1차원의 경우)
로 두고 Runge-Kutta 방법으로 뉴튼의 운동방정식을 적분할 수도 있다. 그렇게 하는 사람 거의 없다.
MD 방법에서는, Predictor-Corrector 방법이나, velocity Verlet 방법을 사용한다.
그 이유는, Runge-Kutta 방법을 사용하면, MD 한 단계, 한 스텦을 진행하는 데, 상대적으로 많은 힘 계산을 요구한다. 물론 정확하다는 장점은 있을 수 있다. 통상, 힘 계산 부분은 계산량이 많다. 분자동역학에서는 이 부분이 전체 계산에 차지하는 비중이 너무나도 크다. 이 방법을 고집할 수 없다.

분자동역학에 이용될 경우, 에너지 드리프트가 발생한다. 국소적으로 표가 나는 것이 아니라 상대적으로 긴 시간이 흐르면 에너지 드리프트가 발생한다. 이점이 Verlet 알고리듬과 대비되는 점이다. Sympletic integrator로서 MD에서 꼭 사용해야 하는 이유가 있는것이다.

http://en.wikipedia.org/wiki/Verlet_integration

cf. LBFGS 를 이용한 함수 최적화 방법:
http://incredible.egloos.com/3042973
Numerical Analysis - Numerical Methods Project
http://math.fullerton.edu/mathews/numerical.html


------------------------------------------------------------
http://numericalmethods.eng.usf.edu/ebooks/runge4th_08ode_ebook.htm
runge4th_08ode_ebook.pdf

핑백

덧글

댓글 입력 영역

최근 포토로그



MathJax