3 September 2021 11:22:11.610 AM bisection_rc_test(): FORTRAN90 version. Test bisection_rc(). test01() Demonstrate bisection_rc() on a simple example. The function is evaluated in a separate routine. I X FX DX 1 0.00000 1.00000 1.00000 2 1.00000 -0.459698 1.00000 3 0.500000 0.377583 0.500000 4 0.750000 -0.183111E-01 0.250000 5 0.625000 0.185963 0.125000 6 0.687500 0.853349E-01 0.625000E-01 7 0.718750 0.338794E-01 0.312500E-01 8 0.734375 0.787473E-02 0.156250E-01 9 0.742188 -0.519571E-02 0.781250E-02 10 0.738281 0.134515E-02 0.390625E-02 11 0.740234 -0.192387E-02 0.195312E-02 12 0.739258 -0.289009E-03 0.976562E-03 13 0.738770 0.528158E-03 0.488281E-03 14 0.739014 0.119597E-03 0.244141E-03 15 0.739136 -0.847007E-04 0.122070E-03 16 0.739075 0.174493E-04 0.610352E-04 17 0.739105 -0.336253E-04 0.305176E-04 18 0.739090 -0.808791E-05 0.152588E-04 19 0.739082 0.468074E-05 0.762939E-05 20 0.739086 -0.170358E-05 0.381470E-05 21 0.739084 0.148858E-05 0.190735E-05 22 0.739085 -0.107502E-06 0.953674E-06 Interval is tiny. A = 0.739084 F(A) = 0.148858E-05 X = 0.739085 F(X) = -0.107502E-06 B = 0.739086 F(B) = -0.170358E-05 test02(): Demonstrate bisection_rc() on a simple example. The function is evaluated within this routine. I X FX DX 1 0.00000 5.00000 1.00000 2 1.00000 -3.13768 1.00000 3 0.500000 -3.03503 0.500000 4 0.250000 4.98958 0.250000 5 0.375000 -2.71136 0.125000 6 0.312500 3.47923 0.625000E-01 7 0.343750 -2.34926 0.312500E-01 8 0.328125 0.872882 0.156250E-01 9 0.335938 -0.922331 0.781250E-02 10 0.332031 -0.384975E-01 0.390625E-02 11 0.330078 0.418288 0.195312E-02 12 0.331055 0.189594 0.976562E-03 13 0.331543 0.754015E-01 0.488281E-03 14 0.331787 0.184063E-01 0.244141E-03 15 0.331909 -0.100581E-01 0.122070E-03 16 0.331848 0.417113E-02 0.610352E-04 17 0.331879 -0.294424E-02 0.305176E-04 18 0.331863 0.613256E-03 0.152588E-04 19 0.331871 -0.116554E-02 0.762939E-05 20 0.331867 -0.276155E-03 0.381470E-05 21 0.331865 0.168548E-03 0.190735E-05 22 0.331866 -0.538042E-04 0.953674E-06 23 0.331866 0.573715E-04 0.476837E-06 24 0.331866 0.178359E-05 0.238419E-06 25 0.331866 -0.260103E-04 0.119209E-06 26 0.331866 -0.121134E-04 0.596046E-07 27 0.331866 -0.516490E-05 0.298023E-07 28 0.331866 -0.169065E-05 0.149012E-07 29 0.331866 0.464671E-07 0.745058E-08 30 0.331866 -0.822093E-06 0.372529E-08 Reached iteration limit. A = 0.331866 F(A) = 0.464671E-07 X = 0.331866 F(X) = -0.822093E-06 B = 0.331866 F(B) = -0.169065E-05 test03() Demonstrate bisection_rc() on a probability example. The cardioid probability density function has a cumulative density function of the form: CDF(X) = ( pi + x - alpha + 2 beta * sin ( x - alpha ) ) / ( 2 * pi ) where alpha and beta are parameters, and x is a value in the range -pi <= x <= +pi. CDF(X) is the probability that a random sample will have a value less than or equal to X. As X moves from -pi to +pi, the CDF rises from 0 (no probability) to 1 (certain probability). Assuming that: * ALPHA = 0.00000 * BETA = 0.250000 determine the value X where the Cardioid CDF is exactly 0.75. I X FX DX 1 -3.14159 -0.750000 6.28319 2 3.14159 0.250000 6.28319 3 0.00000 -0.250000 3.14159 4 1.57080 0.795775E-01 1.57080 5 0.785398 -0.687302E-01 0.785398 6 1.17810 0.110200E-01 0.392699 7 0.981748 -0.275838E-01 0.196350 8 1.07992 -0.794394E-02 0.981748E-01 9 1.12901 0.162468E-02 0.490874E-01 10 1.10447 -0.313822E-02 0.245437E-01 11 1.11674 -0.751383E-03 0.122718E-01 12 1.12287 0.438000E-03 0.613592E-02 13 1.11981 -0.156355E-03 0.306796E-02 14 1.12134 0.140907E-03 0.153398E-02 15 1.12057 -0.770286E-05 0.766990E-03 16 1.12096 0.666073E-04 0.383495E-03 17 1.12076 0.294535E-04 0.191748E-03 18 1.12067 0.108757E-04 0.958738E-04 19 1.12062 0.158648E-05 0.479369E-04 20 1.12060 -0.305817E-05 0.239684E-04 21 1.12061 -0.735842E-06 0.119842E-04 Function is small. A = 1.12060 F(A) = -0.305817E-05 X = 1.12061 F(X) = -0.735842E-06 B = 1.12062 F(B) = 0.158648E-05 Look at the actual cardioid CDF value now: Cardioid( 1.12061 ) = 0.749999 test04(): The freezing pipe problem. At the beginning of a cold spell, the soil is at a uniform temperature of Ti. The cold spell applies a uniform air temperature of Tc, which begins to cool the soil. As a function of depth x and time t, the soil temperature will now cool down as: ( T(x,t) - Tc ) / ( Ti - Tc ) = erf ( 0.5 * x / sqrt ( alpha * t ) ). where: Ti = 20 degrees centigrade, Tc = -15 degrees centigrade, alpha = 0.000000138 meter^2 / second, thermal conductivity; and erf() is the error function. Water freezes at 0 degrees centigrade. What depth x in meters must a water pipe be buried so that it will not freeze even if this cold snap lasts for 60 days? I X FX DX 1 0.00000 -15.0000 1000.00 2 1000.00 20.0000 1000.00 3 500.000 20.0000 500.000 4 250.000 20.0000 250.000 5 125.000 20.0000 125.000 6 62.5000 20.0000 62.5000 7 31.2500 20.0000 31.2500 8 15.6250 20.0000 15.6250 9 7.81250 20.0000 7.81250 10 3.90625 19.9618 3.90625 11 1.95312 16.4124 1.95312 12 0.976562 5.50088 0.976562 13 0.488281 -3.90920 0.488281 14 0.732422 1.08845 0.244141 15 0.610352 -1.34538 0.122070 16 0.671387 -0.111044 0.610352E-01 17 0.701904 0.493193 0.305176E-01 18 0.686646 0.192180 0.152588E-01 19 0.679016 0.408426E-01 0.762939E-02 20 0.675201 -0.350325E-01 0.381470E-02 21 0.677109 0.292217E-02 0.190735E-02 22 0.676155 -0.160509E-01 0.953674E-03 23 0.676632 -0.656328E-02 0.476837E-03 24 0.676870 -0.182029E-02 0.238419E-03 25 0.676990 0.551011E-03 0.119209E-03 26 0.676930 -0.634621E-03 0.596046E-04 27 0.676960 -0.418011E-04 0.298023E-04 28 0.676975 0.254606E-03 0.149012E-04 29 0.676967 0.106403E-03 0.745058E-05 30 0.676963 0.323008E-04 0.372529E-05 Reached iteration limit. A = 0.676960 , F(A) = -0.418011E-04 X = 0.676963 , F(X) = 0.323008E-04 B = 0.676967 , F(B) = 0.106403E-03 test05() Use bisection_rc() on the Kepler equation. 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. Given eccentricity E = 0.100000 Given angle M = 24.8511 (degrees) = 0.433733 (radians) Given E and M, find corresponding X. I X FX DX 1 0.00000 -0.433733 3.14159 2 3.14159 2.70786 3.14159 3 1.57080 1.03706 1.57080 4 0.785398 0.280954 0.785398 5 0.392699 -0.793026E-01 0.392699 6 0.589049 0.997583E-01 0.196350 7 0.490874 0.100008E-01 0.981748E-01 8 0.441786 -0.347024E-01 0.490874E-01 9 0.466330 -0.123643E-01 0.245437E-01 10 0.478602 -0.118521E-02 0.122718E-01 11 0.484738 0.440694E-02 0.613592E-02 12 0.481670 0.161065E-02 0.306796E-02 13 0.480136 0.212664E-03 0.153398E-02 14 0.479369 -0.486286E-03 0.766990E-03 15 0.479752 -0.136814E-03 0.383495E-03 16 0.479944 0.379243E-04 0.191748E-03 17 0.479848 -0.494451E-04 0.958738E-04 18 0.479896 -0.576046E-05 0.479369E-04 19 0.479920 0.160819E-04 0.239684E-04 20 0.479908 0.516072E-05 0.119842E-04 21 0.479902 -0.299868E-06 0.599211E-05 22 0.479905 0.243043E-05 0.299606E-05 23 0.479904 0.106528E-05 0.149803E-05 24 0.479903 0.382706E-06 0.749014E-06 25 0.479903 0.414190E-07 0.374507E-06 26 0.479902 -0.129224E-06 0.187254E-06 27 0.479903 -0.439027E-07 0.936268E-07 28 0.479903 -0.124187E-08 0.468134E-07 29 0.479903 0.200886E-07 0.234067E-07 30 0.479903 0.942335E-08 0.117033E-07 Reached iteration limit. In Radians: A = 0.479903 , F(A) = -0.124187E-08 X = 0.479903 , F(X) = 0.942335E-08 B = 0.479903 , F(B) = 0.200886E-07 In Degrees: A = 27.4964 , F(A) = -0.124187E-08 X = 27.4964 , F(X) = 0.942335E-08 B = 27.4964 , F(B) = 0.200886E-07 bisection_rc_test(): Normal end of execution. 3 September 2021 11:22:11.610 AM