!****************************************************************************** !* Integrate a System of Ordinary Differential Equations By the * !* Runge-Kutta-Fehlberg method (simple or double precision) * !* -------------------------------------------------------------------------- * !* REFERENCE: H A Watts and L F Shampine, * !* Sandia Laboratories, * !* Albuquerque, New Mexico. * !* -------------------------------------------------------------------------- * !* SAMPLE RUN: * !* * !* March 11 2005 9:03:41.625 AM * !* * !* PROGRAM TRKF45 * !* Demonstrate the RKF45 ODE integrator. * !* * !* TEST01 * !* Solve a scalar equation: * !* * !* Y' = 0.25 * Y * ( 1 - Y / 20 ) * !* * !* T Y Y_Exact Error * !* * !* 0.00000 1.00000 1.00000 0.00000 * !* 4.00000 2.50322 2.50322 -0.452995E-05 * !* 8.00000 5.60008 5.60009 -0.104904E-04 * !* 12.0000 10.2777 10.2777 -0.114441E-04 * !* 16.0000 14.8368 14.8368 -0.858307E-05 * !* 20.0000 17.7302 17.7302 0.00000 * !* * !* TEST02 * !* Solve a vector equation: * !* * !* Y'(1) = Y(2) * !* Y'(2) = - Y(1) * !* * !* T Y1 Y2 * !* * !* 0.00000 1.00000 0.00000 * !* 0.523599 0.866026 -0.500000 * !* 1.04720 0.500002 -0.866027 * !* 1.57080 0.599724E-06 -1.00000 * !* 2.09440 -0.500001 -0.866028 * !* 2.61799 -0.866028 -0.500002 * !* 3.14159 -1.00001 -0.801131E-06 * !* 3.66519 -0.866031 0.500002 * !* 4.18879 -0.500004 0.866031 * !* 4.71239 -0.110995E-05 1.00001 * !* 5.23599 0.500003 0.866033 * !* 5.75959 0.866032 0.500005 * !* 6.28319 1.00001 0.100422E-05 * !* * !* TEST03 * !* Solve a scalar equation in double precision: * !* * !* Y' = 0.25 * Y * ( 1 - Y / 20 ) * !* * !* T Y Y_Exact Error * !* * !* 0.00000 1.00000 1.00000 0.00000 * !* 4.00000 2.50322 2.50322 -0.135049E-08 * !* 8.00000 5.60009 5.60009 -0.372025E-08 * !* 12.0000 10.2777 10.2777 0.248324E-09 * !* 16.0000 14.8368 14.8368 -0.957398E-08 * !* 20.0000 17.7302 17.7302 -0.139216E-08 * !* * !* TEST04 * !* Solve a vector equation in double precision: * !* * !* Y'(1) = Y(2) * !* Y'(2) = - Y(1) * !* * !* T Y1 Y2 * !* * !* 0.00000 1.00000 0.00000 * !* 0.523599 0.866025 -0.500000 * !* 1.04720 0.500000 -0.866025 * !* 1.57080 0.183435E-08 -1.00000 * !* 2.09440 -0.500000 -0.866025 * !* 2.61799 -0.866025 -0.500000 * !* 3.14159 -1.00000 -0.366759E-08 * !* 3.66519 -0.866025 0.500000 * !* 4.18879 -0.500000 0.866025 * !* 4.71239 -0.550445E-08 1.00000 * !* 5.23599 0.500000 0.866025 * !* 5.75959 0.866025 0.500000 * !* 6.28319 1.00000 0.734112E-08 * !* * !* TEST05 * !* Solve a vector equation in simple precision: * !* * !* Y'(1) = Y(2) * !* Y'(2) = Y(3) * !* Y'(3) = Y(4) * !* Y'(4) = Y(5) * !* Y'(5) = 45 * y((3) * Y(4) * Y(5) - 40 * Y(4)**3 / (9 * Y(3)**2) * !* * !* T Y1 Y2 Y3 Y4 Y5 * !* 0.00000 1.00000 1.00000 1.00000 1.00000 1.00000 * !* 0.136364 1.14610 1.14609 1.14587 1.14068 1.05604 * !* 0.272727 1.31353 1.31340 1.31128 1.28532 1.05248 * !* 0.409091 1.50538 1.50460 1.49612 1.42333 0.952094 * !* 0.545455 1.72508 1.72223 1.69844 1.53859 0.711111 * !* 0.681818 1.97638 1.96840 1.91370 1.60897 0.288088 * !* 0.818182 2.26328 2.24438 2.13400 1.60781 -0.339182 * !* 0.954545 2.58984 2.55011 2.34770 1.50801 -1.15027 * !* 1.09091 2.96003 2.88369 2.53985 1.28946 -2.06094 * !* 1.22727 3.37739 3.24105 2.69376 0.948201 -2.92046 * !* 1.36364 3.84475 3.61589 2.79372 0.503786 -3.54302 * !* 1.50000 4.36396 4.00000 2.82843 -0.673361E-06 -3.77124 * !* * !* -------------------------------------------------------------------------- * !* * !* NOTE: To link with file rkf45.f90. * !****************************************************************************** program trkf45 implicit none call timestamp ( ) !optional print *,' ' print *,' PROGRAM TRKF45' print *,' Demonstrate the RKF45 ODE integrator.' call test01 print *,' ' pause call test02 print *,' ' pause call test03 print *,' ' pause call test04 print *,' ' pause call test05 print *,' ' pause stop end !of main program subroutine test01 !-------------------------------------------------------------------- ! ! TEST01 solves a scalar ODE. ! !-------------------------------------------------------------------- implicit none ! integer, parameter :: neqn = 1 ! real abserr external f_01 integer i_step integer iflag integer iwork(5) integer n_step real relerr real t real t_out real t_start real t_stop real work(6*neqn+3) real y(neqn) real y_exact_01 ! write ( *, '(a)') ' ' write ( *, '(a)' ) 'TEST01' write ( *, '(a)' ) ' Solve a scalar equation:' write ( *, '(a)') ' ' write ( *, '(a)' ) ' Y'' = 0.25 * Y * ( 1 - Y / 20 )' write ( *, '(a)') ' ' abserr = 0.00001E+00 relerr = 0.00001E+00 iflag = 1 t_start = 0.0E+00 t_stop = 20.0E+00 n_step = 5 t_out = 0.0E+00 y(1) = 1.0E+00 write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' T Y Y_Exact Error' write ( *, '(a)' ) ' ' write ( *, '(4g14.6)' ) t_out, y(1), y_exact_01(t_out), y(1) - y_exact_01(t_out) do i_step = 1, n_step t = ( ( n_step - i_step + 1 ) * t_start & + ( i_step - 1 ) * t_stop ) / real ( n_step ) t_out = ( ( n_step - i_step ) * t_start & + ( i_step ) * t_stop ) / real ( n_step ) call rkf45 ( f_01, neqn, y, t, t_out, relerr, abserr, iflag, work, iwork ) write ( *, '(4g14.6)' ) t_out, y(1), y_exact_01(t_out), & y(1) - y_exact_01(t_out) end do return end subroutine f_01 ( t, y, yp ) !-------------------------------------------------------------------- ! ! F_01 evaluates the derivative for the ODE. ! !-------------------------------------------------------------------- implicit none ! real t real y(1) real yp(1) ! yp(1) = 0.25E+00 * y(1) * ( 1.0E+00 - y(1) / 20.0E+00 ) return end function y_exact_01 ( t ) !-------------------------------------------------------------------- ! ! Y_EXACT_01 evaluates the exact solution of the ODE. ! !-------------------------------------------------------------------- implicit none ! real t real y_exact_01 ! y_exact_01 = 20.0E+00 / ( 1.0E+00 + 19.0E+00 * exp ( - 0.25E+00 * t ) ) return end subroutine test02 !-------------------------------------------------------------------- ! ! TEST02 solves a vector ODE. ! !-------------------------------------------------------------------- implicit none ! integer, parameter :: neqn = 2 ! real abserr external f_02 integer i_step integer iflag integer iwork(5) integer n_step real relerr real t real t_out real t_start real t_stop real work(6*neqn+3) real y(neqn) ! write ( *, '(a)') ' ' write ( *, '(a)' ) 'TEST02' write ( *, '(a)' ) ' Solve a vector equation:' write ( *, '(a)') ' ' write ( *, '(a)' ) ' Y''(1) = Y(2)' write ( *, '(a)' ) ' Y''(2) = - Y(1)' write ( *, '(a)') ' ' abserr = 0.00001E+00 relerr = 0.00001E+00 iflag = 1 t_start = 0.0E+00 t_stop = 2.0E+00 * 3.14159265E+00 n_step = 12 t_out = 0.0E+00 y(1) = 1.0E+00 y(2) = 0.0E+00 write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' T Y1 Y2' write ( *, '(a)' ) ' ' write ( *, '(4g14.6)' ) t_out, y(1), y(2) do i_step = 1, n_step t = ( ( n_step - i_step + 1 ) * t_start & + ( i_step - 1 ) * t_stop ) / real ( n_step ) t_out = ( ( n_step - i_step ) * t_start & + ( i_step ) * t_stop ) / real ( n_step ) call rkf45 ( f_02, neqn, y, t, t_out, relerr, abserr, iflag, work, iwork ) write ( *, '(4g14.6)' ) t_out, y(1), y(2) end do return end subroutine f_02 ( t, y, yp ) !-------------------------------------------------------------------- ! ! F_02 evaluates the derivative for the ODE. ! !-------------------------------------------------------------------- implicit none ! real t real y(2) real yp(2) ! yp(1) = y(2) yp(2) = - y(1) return end subroutine test03 !-------------------------------------------------------------------- ! ! TEST03 solves a scalar ODE in double precision. ! !-------------------------------------------------------------------- implicit none ! integer, parameter :: neqn = 1 ! double precision abserr external f_03 integer i_step integer iflag integer iwork(5) integer n_step double precision relerr double precision t double precision t_out double precision t_start double precision t_stop double precision work(6*neqn+3) double precision y(neqn) double precision y_exact_03 ! write ( *, '(a)') ' ' write ( *, '(a)' ) 'TEST03' write ( *, '(a)' ) ' Solve a scalar equation in double precision:' write ( *, '(a)') ' ' write ( *, '(a)' ) ' Y'' = 0.25 * Y * ( 1 - Y / 20 )' write ( *, '(a)') ' ' abserr = 0.000000001D+00 relerr = 0.000000001D+00 iflag = 1 t_start = 0.0D+00 t_stop = 20.0D+00 n_step = 5 t_out = 0.0D+00 y(1) = 1.0D+00 write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' T Y Y_Exact Error' write ( *, '(a)' ) ' ' write ( *, '(4g14.6)' ) t_out, y(1), y_exact_03(t_out), & y(1) - y_exact_03(t_out) do i_step = 1, n_step t = ( ( n_step - i_step + 1 ) * t_start & + ( i_step - 1 ) * t_stop ) / dble ( n_step ) t_out = ( ( n_step - i_step ) * t_start & + ( i_step ) * t_stop ) / dble ( n_step ) call rkf45_d ( f_03, neqn, y, t, t_out, relerr, abserr, iflag, work, iwork ) write ( *, '(4g14.6)' ) t_out, y(1), y_exact_03(t_out), & y(1) - y_exact_03(t_out) end do return end subroutine f_03 ( t, y, yp ) !-------------------------------------------------------------------- ! ! F_03 evaluates the derivative for the ODE. ! !-------------------------------------------------------------------- implicit none ! double precision t double precision y(1) double precision yp(1) ! yp(1) = 0.25D+00 * y(1) * ( 1.0D+00 - y(1) / 20.0D+00 ) return end function y_exact_03 ( t ) !------------------------------------------------------------------------ ! ! Y_EXACT_03 evaluates the exact solution of the ODE. ! !------------------------------------------------------------------------ implicit none ! double precision t double precision y_exact_03 ! y_exact_03 = 20.0D+00 / ( 1.0D+00 + 19.0D+00 * exp ( - 0.25D+00 * t ) ) return end subroutine test04 !-------------------------------------------------------------------- ! ! TEST04 solves a vector ODE in double precision. ! !-------------------------------------------------------------------- implicit none ! integer, parameter :: neqn = 2 ! double precision abserr external f_04 integer i_step integer iflag integer iwork(5) integer n_step double precision relerr double precision t double precision t_out double precision t_start double precision t_stop double precision work(6*neqn+3) double precision y(neqn) ! write ( *, '(a)') ' ' write ( *, '(a)' ) 'TEST04' write ( *, '(a)' ) ' Solve a vector equation in double precision:' write ( *, '(a)') ' ' write ( *, '(a)' ) ' Y''(1) = Y(2)' write ( *, '(a)' ) ' Y''(2) = - Y(1)' write ( *, '(a)') ' ' abserr = 0.000000001D+00 relerr = 0.000000001D+00 iflag = 1 t_start = 0.0D+00 t_stop = 2.0D+00 * 3.14159265D+00 n_step = 12 t_out = 0.0D+00 y(1) = 1.0D+00 y(2) = 0.0D+00 write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' T Y1 Y2' write ( *, '(a)' ) ' ' write ( *, '(4g14.6)' ) t_out, y(1), y(2) do i_step = 1, n_step t = ( ( n_step - i_step + 1 ) * t_start & + ( i_step - 1 ) * t_stop ) / dble ( n_step ) t_out = ( ( n_step - i_step ) * t_start & + ( i_step ) * t_stop ) / dble ( n_step ) call rkf45_d ( f_04, neqn, y, t, t_out, relerr, abserr, iflag, work, iwork ) write ( *, '(4g14.6)' ) t_out, y(1), y(2) end do return end subroutine f_04 ( t, y, yp ) !-------------------------------------------------------------------- ! ! F_04 evaluates the derivative for the ODE. ! !-------------------------------------------------------------------- implicit none ! double precision t double precision y(2) double precision yp(2) ! yp(1) = y(2) yp(2) = - y(1) return end subroutine test05 !-------------------------------------------------------------------- ! ! TEST05 solves a vector ODE in simple precision. ! !-------------------------------------------------------------------- implicit none ! integer, parameter :: neqn = 5 ! real abserr external f_05 integer i_step integer iflag integer iwork(5) integer n_step real relerr real t real t_out real t_start real t_stop real work(6*neqn+3) real y(neqn) ! write ( *, '(a)') ' ' write ( *, '(a)' ) 'TEST05' write ( *, '(a)' ) ' Solve a vector equation in simple precision:' write ( *, '(a)') ' ' write ( *, '(a)' ) ' Y''(1) = Y(2)' write ( *, '(a)' ) ' Y''(2) = Y(3)' write ( *, '(a)' ) ' Y''(3) = Y(4)' write ( *, '(a)' ) ' Y''(4) = Y(5)' write ( *, '(a)' ) ' Y''(5) = 45 * y((3) * Y(4) * Y(5) - 40 * Y(4)**3 / (9 * Y(3)**2)' write ( *, '(a)') ' ' abserr = 0.000001E+00 relerr = 0.000001E+00 iflag = 1 t_start = 0.0E+00 t_stop = 1.5E+00 n_step = 11 t_out = 0.0E+00 y(1) = 1.0E+00 y(2) = 1.0E+00 y(3) = 1.0E+00 y(4) = 1.0E+00 y(5) = 1.0E+00 write ( *, '(a)' ) ' T Y1 Y2 Y3 Y4 Y5' write ( *, '(6g13.6)' ) t_out, y(1), y(2), y(3), y(4), y(5) do i_step = 1, n_step t = ( ( n_step - i_step + 1 ) * t_start & + ( i_step - 1 ) * t_stop ) / dble ( n_step ) t_out = ( ( n_step - i_step ) * t_start & + ( i_step ) * t_stop ) / dble ( n_step ) call rkf45 ( f_05, neqn, y, t, t_out, relerr, abserr, iflag, work, iwork ) write ( *, '(6g13.6)' ) t_out, y(1), y(2), y(3), y(4), y(5) end do return end subroutine f_05 ( t, y, yp ) !-------------------------------------------------------------------- ! ! F_05 evaluates the derivative for the ODE. ! !-------------------------------------------------------------------- implicit none ! real t real y(5) real yp(5) ! yp(1) = y(2) yp(2) = y(3) yp(3) = y(4) yp(4) = y(5) yp(5) = (45 * y(3) * y(4) * y(5) - 40 * y(4)**3) / (9 * Y(3)**2) return end !end of file trkf45.f90 (J-P Moreau Release 1.1, 03/15/05) !Release 1.1: added TEST #5.