17 January 2015 08:35:27 PM BISECTION_RC_PRB: C++ version. Test the BISECTION_RC library. TEST01 Demonstrate BISECTION_RC on a simple example. The function is evaluated in a separate routine. I X FX DX 1 0 1 1 2 1 -0.459698 1 3 0.5 0.377583 0.5 4 0.75 -0.0183111 0.25 5 0.625 0.185963 0.125 6 0.6875 0.0853349 0.0625 7 0.71875 0.0338794 0.03125 8 0.734375 0.00787473 0.015625 9 0.742188 -0.00519571 0.0078125 10 0.738281 0.00134515 0.00390625 11 0.740234 -0.00192387 0.00195312 12 0.739258 -0.000289009 0.000976562 13 0.73877 0.000528158 0.000488281 14 0.739014 0.000119597 0.000244141 15 0.739136 -8.47007e-05 0.00012207 16 0.739075 1.74493e-05 6.10352e-05 17 0.739105 -3.36253e-05 3.05176e-05 18 0.73909 -8.08791e-06 1.52588e-05 19 0.739082 4.68074e-06 7.62939e-06 20 0.739086 -1.70358e-06 3.8147e-06 21 0.739084 1.48858e-06 1.90735e-06 22 0.739085 -1.07502e-07 9.53674e-07 Interval is tiny. A = 0.739084 F(A) = 1.48858e-06 X = 0.739085 F(X) = -1.07502e-07 B = 0.739086 F(B) = -1.70358e-06 TEST02 Demonstrate BISECTION_RC on a simple example. The function is evaluated within this routine. I X FX DX 1 0 5 1 2 1 -3.13768 1 3 0.5 -3.03503 0.5 4 0.25 4.98958 0.25 5 0.375 -2.71136 0.125 6 0.3125 3.47923 0.0625 7 0.34375 -2.34926 0.03125 8 0.328125 0.872882 0.015625 9 0.335938 -0.922331 0.0078125 10 0.332031 -0.0384975 0.00390625 11 0.330078 0.418288 0.00195312 12 0.331055 0.189594 0.000976562 13 0.331543 0.0754015 0.000488281 14 0.331787 0.0184063 0.000244141 15 0.331909 -0.0100581 0.00012207 16 0.331848 0.00417113 6.10352e-05 17 0.331879 -0.00294424 3.05176e-05 18 0.331863 0.000613256 1.52588e-05 19 0.331871 -0.00116554 7.62939e-06 20 0.331867 -0.000276155 3.8147e-06 21 0.331865 0.000168548 1.90735e-06 22 0.331866 -5.38042e-05 9.53674e-07 23 0.331866 5.73715e-05 4.76837e-07 24 0.331866 1.78359e-06 2.38419e-07 25 0.331866 -2.60103e-05 1.19209e-07 26 0.331866 -1.21134e-05 5.96046e-08 27 0.331866 -5.1649e-06 2.98023e-08 28 0.331866 -1.69065e-06 1.49012e-08 29 0.331866 4.64671e-08 7.45058e-09 30 0.331866 -8.22093e-07 3.72529e-09 Reached iteration limit. A = 0.331866, F(A) = 4.64671e-08 X = 0.331866, F(X) = -8.22093e-07 B = 0.331866, F(B) = -1.69065e-06 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 * BETA = 0.25 determine the value X where the Cardioid CDF is exactly 0.75. I X FX DX 1 -3.14159 -0.75 6.28319 2 3.14159 0.25 6.28319 3 0 -0.25 3.14159 4 1.5708 0.0795775 1.5708 5 0.785398 -0.0687302 0.785398 6 1.1781 0.01102 0.392699 7 0.981748 -0.0275838 0.19635 8 1.07992 -0.00794394 0.0981748 9 1.12901 0.00162468 0.0490874 10 1.10447 -0.00313822 0.0245437 11 1.11674 -0.000751383 0.0122718 12 1.12287 0.000438 0.00613592 13 1.11981 -0.000156355 0.00306796 14 1.12134 0.000140907 0.00153398 15 1.12057 -7.70286e-06 0.00076699 Function is small. A = 1.11981, F(A) = -0.000156355 X = 1.12057, F(X) = -7.70286e-06 B = 1.12134, F(B) = 0.000140907 Look at the actual cardioid CDF value now: Cardioid(1.12057) = 0.749992 BISECTION_RC_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 -15 1000 2 1000 20 1000 3 500 20 500 4 250 20 250 5 125 20 125 6 62.5 20 62.5 7 31.25 20 31.25 8 15.625 20 15.625 9 7.8125 20 7.8125 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.9092 0.488281 14 0.732422 1.08845 0.244141 15 0.610352 -1.34538 0.12207 16 0.671387 -0.111044 0.0610352 17 0.701904 0.493193 0.0305176 18 0.686646 0.19218 0.0152588 19 0.679016 0.0408426 0.00762939 20 0.675201 -0.0350325 0.0038147 21 0.677109 0.00292217 0.00190735 22 0.676155 -0.0160509 0.000953674 23 0.676632 -0.00656328 0.000476837 24 0.67687 -0.00182029 0.000238419 25 0.67699 0.000551011 0.000119209 26 0.67693 -0.000634621 5.96046e-05 27 0.67696 -4.18011e-05 2.98023e-05 28 0.676975 0.000254606 1.49012e-05 29 0.676967 0.000106403 7.45058e-06 30 0.676963 3.23008e-05 3.72529e-06 Reached iteration limit. A = 0.67696, F(A) = -4.18011e-05 X = 0.676963, F(X) = 3.23008e-05 B = 0.676967, F(B) = 0.000106403 TEST05 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.1 Given angle M = 24.8511 (degrees) = 0.433733 (radians) Given E and M, find corresponding X. I X FX DX 1 0 -0.433733 3.14159 2 3.14159 2.70786 3.14159 3 1.5708 1.03706 1.5708 4 0.785398 0.280954 0.785398 5 0.392699 -0.0793026 0.392699 6 0.589049 0.0997583 0.19635 7 0.490874 0.0100008 0.0981748 8 0.441786 -0.0347024 0.0490874 9 0.46633 -0.0123643 0.0245437 10 0.478602 -0.00118521 0.0122718 11 0.484738 0.00440694 0.00613592 12 0.48167 0.00161065 0.00306796 13 0.480136 0.000212664 0.00153398 14 0.479369 -0.000486286 0.00076699 15 0.479752 -0.000136814 0.000383495 16 0.479944 3.79243e-05 0.000191748 17 0.479848 -4.94451e-05 9.58738e-05 18 0.479896 -5.76046e-06 4.79369e-05 19 0.47992 1.60819e-05 2.39684e-05 20 0.479908 5.16072e-06 1.19842e-05 21 0.479902 -2.99868e-07 5.99211e-06 22 0.479905 2.43043e-06 2.99606e-06 23 0.479904 1.06528e-06 1.49803e-06 24 0.479903 3.82706e-07 7.49014e-07 25 0.479903 4.1419e-08 3.74507e-07 26 0.479902 -1.29224e-07 1.87254e-07 27 0.479903 -4.39027e-08 9.36268e-08 28 0.479903 -1.24187e-09 4.68134e-08 29 0.479903 2.00886e-08 2.34067e-08 30 0.479903 9.42335e-09 1.17033e-08 Reached iteration limit. In Radians: A = 0.479903, F(A) = -1.24187e-09 X = 0.479903, F(X) = 9.42335e-09 B = 0.479903, F(B) = 2.00886e-08 In Degrees: A = 27.4964, F(A) = -1.24187e-09 X = 27.4964, F(X) = 9.42335e-09 B = 27.4964, F(B) = 2.00886e-08 BISECTION_RC_PRB: Normal end of execution. 17 January 2015 08:35:27 PM