Here an implementation of the 4th-order Runge-Kutta method based on the subroutine for F90, call rk4, supplied in the book Numerical Recipes in Fortran 90 by Press WH, Teukolsky SA, Vetterling WT and Flannery BP is presented.
This post is devoted to the driver program to integrate a set of two first order ODE as an initial value problem. For reference, the reader is advised to find rk4 in chapter B16 of that book.
The reader is also advised to compare with my raw implementation of the 4th-order Runge-Kutta method based on an algorithm from the book of Numerical methods for Engineers by Chapra SC and Canale RP in the following link,
Contents
______________________________
1. Some comments on rk4 from Numerical Recipes
In comparison with the one I presented in a previous post this routine presents several advantages that worth to mention,
- it uses double precision for calculations. It means that it has a better handling of numerical errors (if such),
- it uses interface blocks. This accounts for usage of variables between the main program and the external functions (that can be located in other files as well),
- it uses a function called assert_eq3 for size checking between the number of variables, ODE and output solutions of the routine. If this check is not tru then ane error message is printed and the program is terminated.
2. The code
The dirver program is not too complicated. In fact it was adapted from that presented by the same authors of the source book Numerical Recipes for in Fortran 90 called: Numerical Recipes. Example Book [Fortran] 2nd ed.. You may look at their driver program, xrk4, in chapter 16.
The proposed driver program in this post is quite different since it was thought for a test with Bessel functions of the first kind. So, formatting of output text and results was changed. Also, the output was conviniently modified so that a table of results could be easily used for comparison with those in the example extracted from the book of Chapra and Canale.
In my case I used four different files:
- main.f90
- MySubroutines.f90. In this file rk4, load and derivs are placed using a module.
- MyTools.f90. In this file the function assert_eq3 is placed using a module. You may find the code for this function in the appendix of the Numerical Recipes in Fortran 90.
- PrecisionConstantsStuff.f90. In this file, the different parameters used for precision handling are declared. This file is similar to nr.f90 presented in the appendix of Numerical Recipes in Fortran 90.
The reader may use the following code for his/her own PrecisionConstantsStuff.f90 file.
module PrecisionConstantsStuff
!
INTEGER, PARAMETER :: LPI = SELECTED_INT_KIND(10)
!
INTEGER, PARAMETER:: SP = KIND(1.0)
!
INTEGER, PARAMETER:: DP = KIND(1.0D0)
!
end module PrecisionConstantsStuff
2.1 Main program
The driver program was also named xrk4. Here it is,
program xrk4
USE MySubroutines
USE PrecisionConstantsStuff
implicit none
! Driver for routine rk4
integer, parameter:: N = 2 ! Number of ODEs
integer, parameter:: NS = 4 ! Number of steps for integration
integer(I4B) :: i,j ! Loop counters
real(DP) :: h, x
REAL(DP) :: dydx(n),yout(n),y(n)
real(dp) :: xi, xf ! Initial and final values of integration range
!
xi = 0.0_dp
xf = 2.0_dp
!Step size calculation
h = (xf-xi)/NS
!
x = xi
!
call Load(y)
!
! Printing results
PRINT*, '4th order Runge-Kutta pogram'
PRINT*, 'Numerical recipes routine'
!
PRINT*, '================================RESULTS=============================='
PRINT*, ' x y1 y2'
PRINT*, '======================================================================'
PRINT*, x, y(1), y(2)
!
call derivs(x,y,dydx)
!
do i=1,NS
x = x + h
call rk4(y,dydx,x,h,yout,derivs)
y(1) = yout(1)
y(2) = yout(2)
call derivs(x,y,dydx)
PRINT*, x, yout(1),yout(2)
end do
!
PRINT*, '======================================================================'
!
endprogram xrk4
2.2 Subroutine load
This subroutine is very similar to those used in previous programs. in the book of Numerical Recipes in Fortran 90 a similar routine is used. This is the routine,
subroutine Load(y)
use PrecisionConstantsStuff
implicit none
real(dp), intent(out):: y(1:2)
!
! In this routine you may store all known boundary conditions at xi
y(1) = 4.0_dp
y(2) = 6.0_dp
!
return
end subroutine Load
2.3 Subroutine derivs
subroutine derivs(x,y,dydx)
USE PrecisionConstantsStuff
implicit none
REAL(DP), INTENT(IN) :: x
REAL(DP), DIMENSION(:), INTENT(IN) :: y
REAL(DP), DIMENSION(:), INTENT(OUT) :: dydx
dydx(1)=-0.5_dp*y(1)
dydx(2)=4.0_dp-0.3_dp*y(2)-0.1_dp*y(1)
return
end subroutine
3. Comparison of results
Here is shown the output of the program. As you can see, the program works well.
 |
| Fig. 01 Comparison of results for the solution of two ODE. Notice that xrk4 is very precise. |
This is the end of the post. The Numerical Recipes in Fortran 90 has several other subroutines to couple with rk4. I shall treat those subroutines in the near future.
Any question? Write in the comments and I shall try to help.
Perhaps the following post may be of interest for you
- Up and running with Numerical Recipes in Fortran 90
- Bracketing roots (zbrac) - Numerical Recipes in F90
- An improvement for bracketing roots (zbrak) - Numerical recipes in F90
- The bisection method (rtbis) - Numerical Recipes in F90
- The false position method (rtfls) - Numerical recipes in F90
- The 4th-order Runge-Kutta method for systems of ode in F90
- Dealing with higher order ODE when using the 4th-order Runge-Kutta mehtod
- Boundary value problem by a shooting method - Case #1 Linear interpolation
=========
Ildebrando
No comments:
Post a Comment