program main !*****************************************************************************80 ! !! bisection_rc_test() tests bisection_rc(). ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 03 September 2021 ! ! Author: ! ! John Burkardt ! implicit none call timestamp ( ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'bisection_rc_test():' write ( *, '(a)' ) ' FORTRAN90 version.' write ( *, '(a)' ) ' Test bisection_rc().' call test01 ( ) call test02 ( ) call test03 ( ) call test04 ( ) call test05 ( ) ! ! Terminate. ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'bisection_rc_test():' write ( *, '(a)' ) ' Normal end of execution.' write ( *, '(a)' ) ' ' call timestamp ( ) stop 0 end subroutine test01 ( ) !*****************************************************************************80 ! !! test01() tests bisection_rc(), evaluating the function in a separate routine. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 03 September 2021 ! ! Author: ! ! John Burkardt ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ) a real ( kind = rk ) b real ( kind = rk ) dx real ( kind = rk ) dx_tol real ( kind = rk ), external :: f01 real ( kind = rk ) fx real ( kind = rk ) fx_tol integer it integer it_max integer job real ( kind = rk ) x write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'test01()' write ( *, '(a)' ) ' Demonstrate bisection_rc() on a simple example.' write ( *, '(a)' ) ' The function is evaluated in a separate routine.' fx_tol = 1.0D-08 dx_tol = 1.0D-06 it = 0 it_max = 30 job = 0 a = 0.0D+00 b = 1.0D+00 write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' I X FX DX' write ( *, '(a)' ) ' ' do call bisection_rc ( a, b, x, fx, job ) if ( job < 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Error return.' exit end if it = it + 1 fx = f01 ( x ) if ( it <= 2 ) then dx = abs ( b - a ) else dx = 0.5D+00 * abs ( b - a ) end if write ( *, '(2x,i4,2x,g14.6,2x,g14.6,2x,g14.6)' ) it, x, fx, dx if ( abs ( fx ) <= fx_tol ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Function is small.' exit end if if ( dx <= dx_tol ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Interval is tiny.' exit end if if ( it_max <= it ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Reached iteration limit.' exit end if end do write ( *, '(a)' ) ' ' write ( *, '(a,g14.6,a,g14.6)' ) ' A = ', a, ' F(A) = ', f01 ( a ) write ( *, '(a,g14.6,a,g14.6)' ) ' X = ', x, ' F(X) = ', f01 ( x ) write ( *, '(a,g14.6,a,g14.6)' ) ' B = ', b, ' F(B) = ', f01 ( b ) return end function f01 ( x ) !*****************************************************************************80 ! !! f01() evaluates the function f(x) = cos ( x ) - x which is zero around 0.74 ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 03 September 2021 ! ! Author: ! ! John Burkardt ! ! Input: ! ! real ( kind = rk ) X, the evaluation point. ! ! Output: ! ! real ( kind = rk ) F01, the function value. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ) f01 real ( kind = rk ) x f01 = cos ( x ) - x return end subroutine test02 ( ) !*****************************************************************************80 ! !! test02() tests bisection_rc(), evaluating the function within the routine. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 03 September 2021 ! ! Author: ! ! John Burkardt ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ) a real ( kind = rk ) b real ( kind = rk ) dx real ( kind = rk ) dx_tol real ( kind = rk ) fx real ( kind = rk ) fx_tol integer it integer it_max integer job real ( kind = rk ) x write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'test02():' write ( *, '(a)' ) ' Demonstrate bisection_rc() on a simple example.' write ( *, '(a)' ) ' The function is evaluated within this routine.' fx_tol = 1.0D-09 dx_tol = 1.0D-09 it = 0 it_max = 30 job = 0 a = 0.0D+00 b = 1.0D+00 write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' I X FX DX' write ( *, '(a)' ) ' ' do call bisection_rc ( a, b, x, fx, job ) if ( job < 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Error return.' exit end if it = it + 1 fx = cos ( 100.0D+00 * x ) - 4.0D+00 * erf ( 30.0D+00 * x - 10.0D+00 ) if ( it <= 2 ) then dx = abs ( b - a ) else dx = 0.5 * abs ( b - a ) end if write ( *, '(2x,i4,2x,g14.6,2x,g14.6,2x,g14.6)' ) it, x, fx, dx if ( abs ( fx ) <= fx_tol ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Function is small.' exit end if if ( dx <= dx_tol ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Interval is tiny.' exit end if if ( it_max <= it ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Reached iteration limit.' exit end if end do write ( *, '(a)' ) ' ' fx = cos ( 100.0D+00 * a ) - 4.0D+00 * erf ( 30.0D+00 * a - 10.0D+00 ) write ( *, '(a,g14.6,a,g14.6)' ) ' A = ', a, ' F(A) = ', fx fx = cos ( 100.0D+00 * x ) - 4.0D+00 * erf ( 30.0D+00 * x - 10.0D+00 ) write ( *, '(a,g14.6,a,g14.6)' ) ' X = ', x, ' F(X) = ', fx fx = cos ( 100.0D+00 * b ) - 4.0D+00 * erf ( 30.0D+00 * b - 10.0D+00 ) write ( *, '(a,g14.6,a,g14.6)' ) ' B = ', b, ' F(B) = ', fx return end subroutine test03 ( ) !*****************************************************************************80 ! !! test03() tests bisection_rc(), to invert the cardioid CDF. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 03 September 2021 ! ! Author: ! ! John Burkardt ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ) a real ( kind = rk ), parameter :: alpha = 0.0D+00 real ( kind = rk ) b real ( kind = rk ), parameter :: beta = 0.25D+00 real ( kind = rk ) cdf real ( kind = rk ) dx real ( kind = rk ) dx_tol real ( kind = rk ) fx real ( kind = rk ) fx_tol integer it integer it_max integer job real ( kind = rk ), parameter :: r8_pi = 3.141592653589793D+00 real ( kind = rk ) x write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'test03()' write ( *, '(a)' ) ' Demonstrate bisection_rc() on a probability example.' write ( *, '(a)' ) '' write ( *, '(a)' ) ' The cardioid probability density function has a' write ( *, '(a)' ) ' cumulative density function of the form:' write ( *, '(a)' ) ' CDF(X) = ( pi + x - alpha + 2 beta * sin ( x - alpha ) ) / ( 2 * pi )' write ( *, '(a)' ) ' where alpha and beta are parameters, and x is a value' write ( *, '(a)' ) ' in the range -pi <= x <= +pi.' write ( *, '(a)' ) '' write ( *, '(a)' ) ' CDF(X) is the probability that a random sample will have' write ( *, '(a)' ) ' a value less than or equal to X.' write ( *, '(a)' ) '' write ( *, '(a)' ) ' As X moves from -pi to +pi,' write ( *, '(a)' ) ' the CDF rises from 0 (no probability)' write ( *, '(a)' ) ' to 1 (certain probability).' write ( *, '(a)' ) '' write ( *, '(a)' ) ' Assuming that:' write ( *, '(a,g14.6)' ) ' * ALPHA = ', alpha write ( *, '(a,g14.6)' ) ' * BETA = ', beta write ( *, '(a)' ) ' determine the value X where the Cardioid CDF is exactly 0.75.' fx_tol = 1.0D-06 dx_tol = 1.0D-08 it = 0 it_max = 30 job = 0 a = - r8_pi b = + r8_pi write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' I X FX DX' write ( *, '(a)' ) ' ' do call bisection_rc ( a, b, x, fx, job ) if ( job < 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Error return.' exit end if it = it + 1 cdf = ( r8_pi + x - alpha + 2.0D+00 * beta * sin ( x - alpha ) ) / ( 2.0D+00 * r8_pi ) fx = cdf - 0.75D+00 if ( it <= 2 ) then dx = abs ( b - a ) else dx = 0.5 * abs ( b - a ) end if write ( *, '(2x,i4,2x,g14.6,2x,g14.6,2x,g14.6)' ) it, x, fx, dx if ( abs ( fx ) <= fx_tol ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Function is small.' exit end if if ( dx <= dx_tol ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Interval is tiny.' exit end if if ( it_max <= it ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Reached iteration limit.' exit end if end do write ( *, '(a)' ) ' ' cdf = ( r8_pi + a - alpha + 2.0D+00 * beta * sin ( a - alpha ) ) / ( 2.0D+00 * r8_pi ) fx = cdf - 0.75D+00 write ( *, '(a,g14.6,a,g14.6)' ) ' A = ', a, ' F(A) = ', fx cdf = ( r8_pi + x - alpha + 2.0D+00 * beta * sin ( x - alpha ) ) / ( 2.0D+00 * r8_pi ) fx = cdf - 0.75D+00 write ( *, '(a,g14.6,a,g14.6)' ) ' X = ', x, ' F(X) = ', fx cdf = ( r8_pi + b - alpha + 2.0D+00 * beta * sin ( b - alpha ) ) / ( 2.0D+00 * r8_pi ) fx = cdf - 0.75D+00 write ( *, '(a,g14.6,a,g14.6)' ) ' B = ', b, ' F(B) = ', fx write ( *, '(a)' ) '' write ( *, '(a)' ) ' Look at the actual cardioid CDF value now:' write ( *, '(a)' ) '' cdf = ( r8_pi + x - alpha + 2.0D+00 * beta * sin ( x - alpha ) ) / ( 2.0D+00 * r8_pi ) write ( *, '(a,g12.6,a,g14.6)' ) ' Cardioid(', x, ') = ', cdf return end subroutine test04 ( ) !*****************************************************************************80 ! !! test04() tests bisection_rc() for the pipe freezing problem. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 03 September 2021 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Cleve Moler, ! Numerical Computing with MATLAB, ! SIAM, 2004, ! ISBN13: 978-0-898716-60-3, ! LC: QA297.M625, ! ebook: https://www.mathworks.com/moler/chapters.html ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ) a real ( kind = rk ) alpha real ( kind = rk ) b real ( kind = rk ) dx real ( kind = rk ) dx_tol real ( kind = rk ) fx real ( kind = rk ) fx_tol integer it integer it_max integer job real ( kind = rk ) t real ( kind = rk ) tc real ( kind = rk ) ti real ( kind = rk ) x write ( *, '(a)' ) '' write ( *, '(a)' ) 'test04():' write ( *, '(a)' ) ' The freezing pipe problem.' write ( *, '(a)' ) '' write ( *, '(a)' ) ' At the beginning of a cold spell, the soil is at a uniform' write ( *, '(a)' ) ' temperature of Ti. The cold spell applies a uniform air' write ( *, '(a)' ) ' temperature of Tc, which begins to cool the soil.' write ( *, '(a)' ) ' As a function of depth x and time t, the soil temperature' write ( *, '(a)' ) ' will now cool down as:' write ( *, '(a)' ) ' ( T(x,t) - Tc ) / ( Ti - Tc ) = erf ( 0.5 * x / sqrt ( alpha * t ) ).' write ( *, '(a)' ) ' where:' write ( *, '(a)' ) ' Ti = 20 degrees centigrade,' write ( *, '(a)' ) ' Tc = -15 degrees centigrade,' write ( *, '(a)' ) ' alpha = 0.000000138 meter^2 / second, thermal conductivity;' write ( *, '(a)' ) ' and erf() is the error function.' write ( *, '(a)' ) ' Water freezes at 0 degrees centigrade.' write ( *, '(a)' ) '' write ( *, '(a)' ) ' What depth x in meters must a water pipe be buried so that it will' write ( *, '(a)' ) ' not freeze even if this cold snap lasts for 60 days?' ! ! Problem parameters. ! ti = 20.0D+00 tc = -15.0D+00 t = 60.0D+00 * 24.0D+00 * 60.0D+00 * 60.0D+00 alpha = 0.000000138D+00 ! ! Iteration parameters. ! fx_tol = 1.0D-09 dx_tol = 1.0D-09 it = 0 it_max = 30 job = 0 fx = 0.0D+00 ! ! Initial guess for interval. ! a = 0.0D+00 b = 1000.0D+00 write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' I X FX DX' write ( *, '(a)' ) ' ' do call bisection_rc ( a, b, x, fx, job ) if ( job < 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Error return.' exit end if it = it + 1 fx = tc + ( ti - tc ) * erf ( 0.5D+00 * x / sqrt ( alpha * t ) ) if ( it <= 2 ) then dx = abs ( b - a ) else dx = 0.5 * abs ( b - a ) end if write ( *, '(2x,i4,2x,g14.6,2x,g14.6,2x,g14.6)' ) it, x, fx, dx if ( abs ( fx ) <= fx_tol ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Function is small.' exit end if if ( dx <= dx_tol ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Interval is tiny.' exit end if if ( it_max <= it ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Reached iteration limit.' exit end if end do write ( *, '(a)' ) '' fx = tc + ( ti - tc ) * erf ( 0.5D+00 * a / sqrt ( alpha * t ) ) write ( *, '(a,g14.6,a,g14.6)' ) ' A = ', a, ', F(A) = ', fx fx = tc + ( ti - tc ) * erf ( 0.5D+00 * x / sqrt ( alpha * t ) ) write ( *, '(a,g14.6,a,g14.6)' ) ' X = ', x, ', F(X) = ', fx fx = tc + ( ti - tc ) * erf ( 0.5D+00 * b / sqrt ( alpha * t ) ) write ( *, '(a,g14.6,a,g14.6)' ) ' B = ', b, ', F(B) = ', fx return end subroutine test05 ( ) !*****************************************************************************80 ! !! test05() tests bisection_rc() for Kepler's problem. ! ! Discussion: ! ! Kepler's equation has the form: ! ! X = M + E * sin ( X ) ! ! X represents the eccentric anomaly of a planet, the angle between the ! perihelion (the point on the orbit nearest to the sun) through the sun ! to the center of the ellipse, and the line from the center of the ellipse ! to the planet. ! ! There are two parameters, E and M: ! ! * E is the eccentricity of the orbit, which should be between 0 and 1.0; ! ! * M is the angle from the perihelion made by a fictitious planet traveling ! on a circular orbit centered at the sun, and traveling at a constant ! angular velocity equal to the average angular velocity of the true ! planet. M is usually between 0 and 180 degrees, but can have any value. ! ! For convenience, X and M are measured in degrees. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 03 September 2021 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Cleve Moler, ! Numerical Computing with MATLAB, ! SIAM, 2004, ! ISBN13: 978-0-898716-60-3, ! LC: QA297.M625, ! ebook: https://www.mathworks.com/moler/chapters.html ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ) ad real ( kind = rk ) ar real ( kind = rk ) bd real ( kind = rk ) br real ( kind = rk ) dx real ( kind = rk ) dx_tol real ( kind = rk ) e real ( kind = rk ) fx real ( kind = rk ) fx_tol integer it integer it_max integer job real ( kind = rk ) md real ( kind = rk ) mr real ( kind = rk ), parameter :: r8_pi = 3.141592653589793D+00 real ( kind = rk ) xd real ( kind = rk ) xr write ( *, '(a)' ) '' write ( *, '(a)' ) 'test05()' write ( *, '(a)' ) ' Use bisection_rc() on the Kepler equation.' write ( *, '(a)' ) '' write ( *, '(a)' ) ' Kepler''s equation has the form'; write ( *, '(a)' ) '' write ( *, '(a)' ) ' X = M + E * sin ( X )' write ( *, '(a)' ) '' write ( *, '(a)' ) ' X represents the eccentric anomaly of a planet, the angle between the' write ( *, '(a)' ) ' perihelion (the point on the orbit nearest to the sun) through the sun' write ( *, '(a)' ) ' to the center of the ellipse, and the line from the center of the ellipse' write ( *, '(a)' ) ' to the planet.' write ( *, '(a)' ) '' write ( *, '(a)' ) ' There are two parameters, E and M:' write ( *, '(a)' ) '' write ( *, '(a)' ) ' * E is the eccentricity of the orbit, which should be between 0 and 1.0;' write ( *, '(a)' ) '' write ( *, '(a)' ) ' * M is the angle from the perihelion made by a fictitious planet traveling' write ( *, '(a)' ) ' on a circular orbit centered at the sun, and traveling at a constant' write ( *, '(a)' ) ' angular velocity equal to the average angular velocity of the true' write ( *, '(a)' ) ' planet. M is usually between 0 and 180 degrees, but can have any value.' write ( *, '(a)' ) '' write ( *, '(a)' ) ' For convenience, X and M are measured in degrees.' ! ! Problem parameters. ! md = 24.851090D+00 mr = md * r8_pi / 180.0D+00 e = 0.1D+00 write ( *, '(a)' ) '' write ( *, '(a,g14.6)' ) ' Given eccentricity E = ', e write ( *, '(a,g14.6,a)' ) ' Given angle M = ', md, ' (degrees)' write ( *, '(a,g14.6,a)' ) ' = ', mr, ' (radians)' write ( *, '(a)' ) '' write ( *, '(a)' ) ' Given E and M, find corresponding X.' ! ! Iteration parameters. ! fx_tol = 1.0D-09 dx_tol = 1.0D-09 it = 0 it_max = 30 job = 0 fx = 0.0D+00 ! ! Initial guess for interval. ! ad = 0.0D+00 bd = 180.0D+00 ar = ad * r8_pi / 180.0D+00 br = bd * r8_pi / 180.0D+00 write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' I X FX DX' write ( *, '(a)' ) ' ' do call bisection_rc ( ar, br, xr, fx, job ) if ( job < 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Error return.' exit end if it = it + 1 fx = xr - mr - e * sin ( xr ) if ( it <= 2 ) then dx = abs ( br - ar ) else dx = 0.5 * abs ( br - ar ) end if write ( *, '(2x,i4,2x,g14.6,2x,g14.6,2x,g14.6)' ) it, xr, fx, dx if ( abs ( fx ) <= fx_tol ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Function is small.' exit end if if ( dx <= dx_tol ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Interval is tiny.' exit end if if ( it_max <= it ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Reached iteration limit.' exit end if end do write ( *, '(a)' ) '' write ( *, '(a)' ) ' In Radians:' write ( *, '(a)' ) '' fx = ar - mr - e * sin ( ar ); write ( *, '(a,g14.6,a,g14.6)' ) ' A = ', ar, ', F(A) = ', fx fx = xr - mr - e * sin ( xr ); write ( *, '(a,g14.6,a,g14.6)' ) ' X = ', xr, ', F(X) = ', fx fx = br - mr - e * sin ( br ); write ( *, '(a,g14.6,a,g14.6)' ) ' B = ', br, ', F(B) = ', fx ad = ar * 180.0D+00 / r8_pi xd = xr * 180.0D+00 / r8_pi bd = br * 180.0D+00 / r8_pi write ( *, '(a)' ) '' write ( *, '(a)' ) ' In Degrees:' write ( *, '(a)' ) '' fx = ( ad - md ) * r8_pi / 180.0D+00 - e * sin ( ad * r8_pi / 180.0D+00 ); write ( *, '(a,g14.6,a,g14.6)' ) ' A = ', ad, ', F(A) = ', fx fx = ( xd - md ) * r8_pi / 180.0D+00 - e * sin ( xd * r8_pi / 180.0D+00 ); write ( *, '(a,g14.6,a,g14.6)' ) ' X = ', xd, ', F(X) = ', fx fx = ( bd - md ) * r8_pi / 180.0D+00 - e * sin ( bd * r8_pi / 180.0D+00); write ( *, '(a,g14.6,a,g14.6)' ) ' B = ', bd, ', F(B) = ', fx return end