RK4

MODULE Global_data  !Symbolic names for kind types of single- and double-precision reals:   INTEGER, PARAMETER :: SP = KIND(1.0)   INTEGER, PARAMETER :: DP = KIND(1.0D0)   !Frequently used mathematical constants (with precision to spare):   REAL(DP), PARAMETER :: Pi=3.141592653589793238462643383279502884197_dp   ! order of the problem  INTEGER, PARAMETER :: no_of_equations=2  ! parameters for the equation  REAL(DP) :: y_initial(no_of_equations)  REAL(DP) :: x_start,x_finish,h,kappaEND MODULE Global_dataMODULE RungeKutta_method! here we just need dp, kappa, h and no_of_equations! if we set wp=>sp we can change to single precision numbers  USE global_data, ONLY : wp=>dp,h,kappa,no_of_equations  IMPLICIT NONECONTAINS  FUNCTION func(x,y)!    we assign y to have 'no_of_equations's elements throughout to keep the program consistent    REAL(wp) , INTENT(IN) :: x,y(no_of_equations)    REAL(wp) :: func(no_of_equations)!   The ODE here is y'' + kappa y' + x y = 0!   so that the set of first ODEs is!    y_1' = y_2!    y_2' = - kappa y_2 - x y_1! where we write y_1 = y and y_2 = y'    func(1) = y(2)    func(2) = -kappa*y(2) - x*y(1)  END FUNCTION func  FUNCTION rungekutta4(x,y)    REAL(wp) , INTENT(IN) :: x,y(no_of_equations)    REAL(wp),DIMENSION(no_of_equations) :: rungekutta4,k1,k2,k3,k4    ! Note here that k1,k2,... etc are two element arrays.    k1 = h*func(x,y)    k2 = h*func(x+0.5_wp*h,y+0.5_wp*k1)    k3 = h*func(x+0.5_wp*h,y+0.5_wp*k2)    k4 = h*func(x+h,y+k3)    rungekutta4 = y + (k1+2._wp*(k2+k3) + k4)/6._wp  END FUNCTION rungekutta4END MODULE RungeKutta_methodPROGRAM euler  ! we need all the variables from global data  USE global_data, ONLY : wp=>dp,kappa,h,x_start,x_finish,y_initial,no_of_equations  ! RungeKutta_method gives us our rungekutta method function  USE RungeKutta_method  IMPLICIT NONE    INTEGER :: i,no_of_steps  REAL(wp) :: x,y(no_of_equations)    ! Setup initial conditions  ! set x_0 = 0  x_start = 0._wp  ! set x_n = 1  x_finish = 1._wp  ! Set y(x=0) = 0  y_initial(1)=0._wp    ! Set y'(x=0) = 1  y_initial(2)=1._wp  ! Set the number of steps  no_of_steps = 10  ! Set the parameter kappa from the ODE  kappa = 5._wp  ! Set the step size  h = (x_finish-x_start)/dble(no_of_steps)  ! Open a file to write to and output some info to screen  WRITE(6,*)' Run with ::',no_of_steps,' steps.'  OPEN(unit=10,file="runge_kutta.dat")  ! Set the start values of x and y  x = x_start  y = y_initial  WRITE(6,'(3(E15.8,1X))')x,y  WRITE(10,'(3(E15.8,1X))')x,y      ! Loop over the required no. of steps to reach x_finish  DO i=1,no_of_steps! update the values of x and y     y = rungekutta4(x,y)     x = x + h     ! pick out only 10 value to print at.     if(mod(i,no_of_steps/10)==0)THEN     ! print to screen and file         WRITE(6,'(3(E15.8,1X))')x,y         WRITE(10,'(3(E15.8,1X))')x,y     END if          END DO  CLOSE(unit=10)  END PROGRAM euler