program main !*****************************************************************************80 ! !! sde_test() tests sde(). ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 21 September 2012 ! ! Author: ! ! John Burkardt ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) call timestamp ( ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'sde_test():' write ( *, '(a)' ) ' Fortran90 version' write ( *, '(a)' ) ' Test sde().' call test01 ( ) call test02 ( ) call test03 ( ) call test04 ( ) call test05 ( ) call test06 ( ) call test07 ( ) call test08 ( ) call test09 ( ) call test10 ( ) call test11 ( ) ! ! Terminate. ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SDE_TEST' write ( *, '(a)' ) ' Normal end of execution.' write ( *, '(a)' ) ' ' call timestamp ( ) stop 0 end subroutine test01 ( ) !*****************************************************************************80 ! !! TEST01 tests BPATH. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 14 September 2012 ! ! Author: ! ! John Burkardt ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer, parameter :: n = 500 real ( kind = rk ) w(0:n) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST01:' write ( *, '(a)' ) ' BPATH generates a sample Brownian motion path' call bpath ( n, w ) call bpath_gnuplot ( n, w ) return end subroutine test02 ( ) !*****************************************************************************80 ! !! TEST02 tests BPATH_AVERAGE. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 14 September 2012 ! ! Author: ! ! John Burkardt ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer, parameter :: m = 1000 integer, parameter :: n = 500 real ( kind = rk ) error real ( kind = rk ), allocatable :: u(:,:) real ( kind = rk ) umean(0:n) allocate ( u(1:m,0:n) ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST02:' write ( *, '(a)' ) ' BPATH_AVERAGE generates many Brownian paths' write ( *, '(a)' ) ' and averages them.' call bpath_average ( m, n, u, umean, error ) call bpath_average_gnuplot ( m, n, u, umean ) deallocate ( u ) return end subroutine test03 ( ) !*****************************************************************************80 ! !! TEST03 tests CHAIN. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 15 September 2012 ! ! Author: ! ! John Burkardt ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer, parameter :: n = 200 real ( kind = rk ) diff real ( kind = rk ) vem(0:n) real ( kind = rk ) xem(0:n) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST03:' write ( *, '(a)' ) ' CHAIN solves a stochastic differential equation for' write ( *, '(a)' ) ' a function of a stochastic variable X.' write ( *, '(a)' ) ' We can solve for X(t), and then evaluate V(X(t)).' write ( *, '(a)' ) ' Or, we can apply the stochastic chain rule to derive an' write ( *, '(a)' ) ' an SDE for V, and solve that.' call chain ( n, xem, vem, diff ) write ( *, '(a)' ) ' ' write ( *, '(a,g14.6)' ) ' Maximum | Sqrt(X) - V | = ', diff call chain_gnuplot ( n, xem, vem ) return end subroutine test04 ( ) !*****************************************************************************80 ! !! TEST04 tests EM. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 15 September 2012 ! ! Author: ! ! John Burkardt ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer, parameter :: n = 256 real ( kind = rk ) diff real ( kind = rk ) t(0:n) real ( kind = rk ) t2(0:n/4) real ( kind = rk ) xtrue(0:n) real ( kind = rk ) xem(0:n/4) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST04:' write ( *, '(a)' ) ' EM solves a stochastic differential equation' write ( *, '(a)' ) ' using the Euler-Maruyama method.' call em ( n, t, xtrue, t2, xem, diff ) write ( *, '(a)' ) ' ' write ( *, '(a,g14.6)' ) ' | Exact X(T) - EM X(T) | = ', diff call em_gnuplot ( n, t, xtrue, t2, xem ) return end subroutine test05 ( ) !*****************************************************************************80 ! !! TEST05 tests EMSTRONG. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 16 September 2012 ! ! Author: ! ! John Burkardt ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer, parameter :: m = 100 integer, parameter :: n = 512 integer, parameter :: p_max = 6 real ( kind = rk ) dtvals(p_max) real ( kind = rk ) xerr(p_max) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST05:' write ( *, '(a)' ) ' EMSTRONG investigates the strong convergence' write ( *, '(a)' ) ' of the Euler-Maruyama method.' call emstrong ( m, n, p_max, dtvals, xerr ) call emstrong_gnuplot ( p_max, dtvals, xerr ) return end subroutine test06 ( ) !*****************************************************************************80 ! !! TEST06 tests EMWEAK. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 17 September 2012 ! ! Author: ! ! John Burkardt ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer, parameter :: m = 50000 integer, parameter :: p_max = 5 real ( kind = rk ) dtvals(p_max) integer method real ( kind = rk ) xerr(p_max) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST06:' write ( *, '(a)' ) ' EMWEAK investigates the weak convergence' write ( *, '(a)' ) ' of the Euler-Maruyama method.' method = 0 call emweak ( method, m, p_max, dtvals, xerr ) call emweak_gnuplot ( p_max, dtvals, xerr, method ) method = 1 call emweak ( method, m, p_max, dtvals, xerr ) call emweak_gnuplot ( p_max, dtvals, xerr, method ) return end subroutine test07 ( ) !*****************************************************************************80 ! !! TEST07 tests MILSTRONG. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 18 September 2012 ! ! Author: ! ! John Burkardt ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer, parameter :: p_max = 4 real ( kind = rk ) dtvals(p_max) real ( kind = rk ) xerr(p_max) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST07:' write ( *, '(a)' ) ' MILSTRONG investigates the strong convergence' write ( *, '(a)' ) ' of the Milstein method.' call milstrong ( p_max, dtvals, xerr ) call milstrong_gnuplot ( p_max, dtvals, xerr ) return end subroutine test08 ( ) !*****************************************************************************80 ! !! TEST08 tests STAB_ASYMPTOTIC. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 18 September 2012 ! ! Author: ! ! John Burkardt ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer, parameter :: n = 1000 integer, parameter :: p_max = 3 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST08:' write ( *, '(a)' ) ' STAB_ASYMPTOTIC investigates the asymptotic' write ( *, '(a)' ) ' stability of the Euler-Maruyama method.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' For technical reasons, the plotting is done' write ( *, '(a)' ) ' in the same routine as the computations.' call stab_asymptotic ( n, p_max ) return end subroutine test09 ( ) !*****************************************************************************80 ! !! TEST09 tests STAB_MEANSQUARE. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 26 September 2012 ! ! Author: ! ! John Burkardt ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST09:' write ( *, '(a)' ) ' STAB_MEANSQUARE investigates the mean square' write ( *, '(a)' ) ' stability of the Euler-Maruyama method.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' For technical reasons, the plotting is done' write ( *, '(a)' ) ' in the same routine as the computations.' call stab_meansquare ( ) return end subroutine test10 ( ) !*****************************************************************************80 ! !! TEST10 tests STOCHASTIC_INTEGRAL_ITO. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 21 September 2012 ! ! Author: ! ! John Burkardt ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ) error real ( kind = rk ) estimate real ( kind = rk ) exact integer i integer n write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST10:' write ( *, '(a)' ) ' Estimate the Ito integral of W(t) dW over [0,1].' write ( *, '(a)' ) ' ' write ( *, '(a)' ) & ' Abs Rel' write ( *, '(a)' ) & ' N Exact Estimate Error Error' write ( *, '(a)' ) ' ' n = 100 do i = 1, 7 call stochastic_integral_ito ( n, estimate, exact, error ) write ( *, '(2x,i8,2x,g16.8,2x,g16.8,2x,g10.2,2x,g10.2)' ) & n, exact, estimate, error, error / exact n = n * 4 end do return end subroutine test11 ( ) !*****************************************************************************80 ! !! TEST11 tests STOCHASTIC_INTEGRAL_STRAT. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 21 September 2012 ! ! Author: ! ! John Burkardt ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ) error real ( kind = rk ) estimate real ( kind = rk ) exact integer i integer n write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST11:' write ( *, '(a)' ) & ' Estimate the Stratonovich integral of W(t) dW over [0,1].' write ( *, '(a)' ) ' ' write ( *, '(a)' ) & ' Abs Rel' write ( *, '(a)' ) & ' N Exact Estimate Error Error' write ( *, '(a)' ) ' ' n = 100 do i = 1, 7 call stochastic_integral_strat ( n, estimate, exact, error ) write ( *, '(2x,i8,2x,g16.8,2x,g16.8,2x,g10.2,2x,g10.2)' ) & n, exact, estimate, error, error / exact n = n * 4 end do return end subroutine timestamp ( ) !*****************************************************************************80 ! !! timestamp() prints the current YMDHMS date as a time stamp. ! ! Example: ! ! 31 May 2001 9:45:54.872 AM ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 18 May 2013 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! None ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) character ( len = 8 ) ampm integer d integer h integer m integer mm character ( len = 9 ), parameter, dimension(12) :: month = (/ & 'January ', 'February ', 'March ', 'April ', & 'May ', 'June ', 'July ', 'August ', & 'September', 'October ', 'November ', 'December ' /) integer n integer s integer values(8) integer y call date_and_time ( values = values ) y = values(1) m = values(2) d = values(3) h = values(5) n = values(6) s = values(7) mm = values(8) if ( h < 12 ) then ampm = 'AM' else if ( h == 12 ) then if ( n == 0 .and. s == 0 ) then ampm = 'Noon' else ampm = 'PM' end if else h = h - 12 if ( h < 12 ) then ampm = 'PM' else if ( h == 12 ) then if ( n == 0 .and. s == 0 ) then ampm = 'Midnight' else ampm = 'AM' end if end if end if write ( *, '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) & d, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, trim ( ampm ) return end