question for lyuka and the Ostrowski method « Next Oldest | Next Newest »

 ▼ hugh steers Senior Member Posts: 536 Threads: 56 Joined: Jul 2005 08-16-2011, 07:13 PM hi, wanted to try out this method. took the code here, http://www.finetune.jp/~lyuka/technote/fract/ost.html and tried one of my tests. this is smooth test from the temperature industry. ```double f(double t) { /* example from pressure */ double x; t = t + 273.15; x = 0; if (t > 0) { x = -6096.9385 / t + 16.635794 - 0.02711193 * t + 0.00001673952 * t * t + 2.433502 * log(t); x = exp(x); } return x - 1e-5; // value = 1e-5 } ``` not very pretty i know, but it should have a real world solution to a value of 1e-5 (hardcoded) between -275.15 and 100, say. i do this: ```int main() { double a, b, root; int r; printf("vwat\n"); r= ost(a=-271.15, b=100.0, f, EPSILON, ITERATION, &root); printf("main: %g, %g, f(%.18g) = %g, %d\n", a, b, root, f(root), r); return 0; } ``` and it gives me -256. here's the results from my other methods, ```37 iterations solve vwat=1e-5 in [-273.150000, 100.000000] with brentRoot = -106.310044 100 iterations solve vwat=1e-5 in [-273.150000, 100.000000] with secantRoot = -1.#IND00 8 iterations solve vwat=1e-5 in [-273.150000, 100.000000] with ridderRoot = -106.310044 2 iterations solve vwat=1e-5 in [-273.150000, 100.000000] with ostRoot = -256.086474 ``` secant blows, but brent and ridder are fine. ost is the -256 as mentioned. is there a way to fix for this, otherwise the method has failed. thanks for any info, -- hugh. ▼ Namir Posting Freak Posts: 2,247 Threads: 200 Joined: Jun 2005 08-16-2011, 10:12 PM Hugh, Thanks for the nice link. Interestingly the author uses my favorite test function f(x) = e^x - 3*x^2. Small world! :-) I am going to look into the code and write a version for Excel VBA. Namir ▼ Lyuka Member Posts: 169 Threads: 12 Joined: Aug 2007 08-16-2011, 11:44 PM Hi, Namir It's not a small world, I've chosen the function because you often use the function as an example :-) Best regards, Lyuka ▼ Namir Posting Freak Posts: 2,247 Threads: 200 Joined: Jun 2005 08-17-2011, 01:11 AM Hi Lyuka, Thank you for a wonderful web page and links. The C and 42S source codes (plus the .raw file for the 42s) are very very nice. Namir Edited: 17 Aug 2011, 1:11 a.m. Namir Posting Freak Posts: 2,247 Threads: 200 Joined: Jun 2005 08-17-2011, 07:58 AM Lyuka's article about the Ostrowski algorithm comes in an interesting time. For the last week, I have been playing with basic root-seeking algorithms trying to use random numbers as a tool to tweak and improve these algorithms. Here is the pseudo-code for a method that is interesting. It is based on Bisection and uses random numbers: ```Given: 1. Function Fx(x)=0 2. Function Rand() that generates uniform random numbers in interval (0,1). 3. Interval {A,B] containing the root of function Fx(x)=0, such that A < B. 4. The Factor value used to tweak root guesses. Optimum value is around 0.75. 5. Toler is smallest size of interval [A.B] used to obtain the root. FA = Fx(A) FB = Fx(B) bFlag = False ' compensation mode flag Do H = Abs(B - A) / 2 ' do we need to compensate for the last iteration's failure ' to shrink the size of interval [A,B] by more than half? If bFlag Then ' compare the absolute function values at A and B to determine ' which end of the interval we will shift the new guess ' towards, in hope of getting closer to the root If Abs(FA) < Abs(FB) Then ' closer to A X = A + H * Rand() * Factor Else ' closer to B X = B - H * Rand() * Factor End If Else ' get a random value inside interval {A, B] replacing the ' typical Bisection step of calculating X as (A+B)/2 ' ' The new value of X reflects success using random numbers ' if it shrinks the size of the interval [A,B] by more than ' half. X = A + (B - A) * Rand() End If FX = Fx(X) If FA * FX > 0 Then A = X FA = FX Else B = X FB = FX End If ' did the last iteration fail to shrink the size of ' the interval [A,B] by more than half? ' If yes, then random selection has failed and we need to ' compensate in the next iteration bFlag = Abs(B - A) > H Loop Until Abs(A - B) <= Toler Refined guess for root = X ``` I have been running solutions in batches of 10000 iterations. The above algorithm does better than the Bisection method between 20% to 50% of the time and saving between 1 to 8 iterations (in significant frequencies). The results depend on the interval [A,B] and where the root is located within the interval. If the root is near the initial boundaries of interval [A,B] then the above algorithm shows more success. Using the "Factor" is the latest aspect of tweaking the algorithm. I am still studying the effect of using Factor and its optimum values. The above algorithm is not the greatest improvement on Bisection. Like I said before I am looking at using random numbers in tweaking basic algorithms I have used the approach with Newton and Halley's method. So far no success. Namir Edited: 17 Aug 2011, 8:15 a.m. Lyuka Member Posts: 169 Threads: 12 Joined: Aug 2007 08-16-2011, 11:40 PM Hi hugh, Thank you for your interest on the page " Ostrowski's method " and "SOLVE alternative for the HP 42S " (necessary not to be OT ;-) Take a look at the return code, it would be '3' for your initial guesses of {-271.15, 100.0}, which means Quote: 3 : The search concentrated in a local flat region of the function. f(x_{n+2}) == f(x_{n+1}) It would occur mainly because of too wide set of initial guesses. If you try that of {-271.15, 0} or {-100, 100}, the solve will converge to the correct answer of -106.310044029008282. Quote:```main: -271.15, 100, f(-256.651643757074339) = -1e-05, 3 main: -271.15, 0, f(-106.310044029008282) = 1.52466e-20, 1 main: -100, 100, f(-106.310044029008282) = 1.52466e-20, 1 ``` Best regards, Lyuka Edited: 17 Aug 2011, 1:38 a.m. ▼ Namir Posting Freak Posts: 2,247 Threads: 200 Joined: Jun 2005 08-17-2011, 01:48 AM I like your C code because it examines the function values for the initial guesses and is able to exit if any of these guesses are at or close to the root. This is what I call the "deluxe version" of an implementation where all scenarios are handled. Namir ▼ Lyuka Member Posts: 169 Threads: 12 Joined: Aug 2007 08-17-2011, 02:05 AM It's my pleasure, thank you for saying that :-) Lyuka ▼ Namir Posting Freak Posts: 2,247 Threads: 200 Joined: Jun 2005 08-17-2011, 08:56 AM Hi, I was looking at the following part of your C code: ``` if (0.0 == (f - e) * (f - d)) { c = (a + 3 * b) * 0.25; if (0.0 == (f = (*func)(c))) { *x = c; return 0; } if (0.0 == (f - e) * (f - d)) { *x = c; return 7; } } ``` How about something like: ``` while (0.0 == (f - e) * (f - d)) { c = a + (b-a)*rand()/RAND_MAX; if (0.0 == (f = (*func)(c))) { *x = c; return 0; } } ``` The suggested code calculates c as a random value in the interval [a,b]. Namir ▼ Lyuka Member Posts: 169 Threads: 12 Joined: Aug 2007 08-17-2011, 09:39 AM Hi, It seems good an idea, however the iteration limit must be used. If the condition kept false several times, it must be concerned that it would be occurring because of too narrow initial guesses, or not. In such case, the third guess 'c' would better be put outside of the 'a' to 'b' region. This kind of 'pivot' strategy is difficult to decide. Best regards, Lyuka ▼ Namir Posting Freak Posts: 2,247 Threads: 200 Joined: Jun 2005 08-17-2011, 10:37 AM You are right. After I posted the message I thought teh while loop perhaps needed a counter? Also the range for calculating c may need to be wider. You are not obligated to have a < c < b. ▼ Lyuka Member Posts: 169 Threads: 12 Joined: Aug 2007 08-18-2011, 04:53 PM Hi Namir, Inspired by your hints of using random() function, I modified the _ost() function to deal with the case f(x_n+2) == f(x_n+1) by pivoting the value to the random position between two previous values. Quote:```if (0 == (f - e)) { if (p++ > MAXPVT) { *x = c; return 3; } c = a + (b - a) * (1.0 + (double)random()) / M2E31; f = (*func)(c); continue; } ``` It seems to work well, also for the case of hugh steer's temperature related function. Best regards, Lyuka Edited: 18 Aug 2011, 5:10 p.m. Namir Posting Freak Posts: 2,247 Threads: 200 Joined: Jun 2005 08-16-2011, 11:43 PM Hugh, I got to play with Ostrowski's method using f(x)=e^x-3*x^2 and comparing it with Newton's method and the faster Halley's method. Since Ostrowski's method uses two initial guesses (and calculates a third guess as the average of the first two), I used the average of these guesses for the Newton and Halley methods, to even up the score a bit. In this case, Ostrowski's method did better than Newton but not better than Halley. If I use the first guess for Newton and Halley methods, then depending on that guess and the second guess (used only by Ostrowski's method) I get a few cases where Ostrowski's method outperforms the other two! Thanks for the fantastic article (and thanks for Luyka for writing that web page). It is a true gem in numerical analysis and for a root-seeking junkie like me!! Namir Edited: 17 Aug 2011, 1:23 a.m. ▼ hugh steers Senior Member Posts: 536 Threads: 56 Joined: Jul 2005 08-17-2011, 07:12 AM all credit goes to lyuka for the article. i followed his post earlier (seems to have disappeared) to investigate this method. yes, i am interested in root methods too. mostly i use ridders method and find it very stable. the code is also small. i wanted to try this method out and compare it. to lyuka, thanks. it looks like the method finds my function flat. it is true that the curve has a flattened end. this has fooled other methods too. the annoying thing is that this particular example comes from a real world case. now im thinking if there's a hybrid method involving ostrowski and other ideas. -- hugh. Lyuka Member Posts: 169 Threads: 12 Joined: Aug 2007 08-22-2011, 04:40 PM Hi hugh steers, I've been tried to make my implementation of Ostrowski's method to be able to solve the water saturation pressure related equations, as you shown above, for various initial guesses. Now, the latest version of the program (Rev.2.8) seems to work well for such kind of equations with almost any initial guesses within the range of [-273.15 : 100], except when the two guesses are at flat region near -273.15. I've tested the program using the water saturation pressures equations by Sonntag and Wexler as well, to solve it for the pressure to be 1e-5. Quote: ```double s(x) double x; { double t; double p = 1e-5; t = x + 273.15; if (t < 0) return 0; return exp(((16.635794 + 2.433502 * log(t) + 1.673952e-5 * t * t) - 2.711193e-2 * t) - 6096.9385 / t) * 100.0 - p; } double w(x) double x; { double t; double p = 1e-5; t = x + 273.15; if (t < 0) return 0; return exp(((18.87643845 + 2.858487 * log(t)) + t * (-2.8354721e-2 + t * (1.7838301e-5 + t * (-8.4150417e-10 + t * 4.4412543e-13)))) + (-6017.0128 -2991.2729 / t) / t) - p; } ``` Here is the results of test run of the program. Quote: ```main: s(t) is 'P_water_sat_Sonntag' main: 273.15, 100, s(-125.417920181) = 2.37609682364e-17, r=1 main: -125.986, 40.366, s(-125.417920181) = -4.23516473627e-20, r=1 main: 24.7879, 19.0635, s(-125.417920181) = 2.87991202066e-20, r=1 main: -199.434, 67.0312, s(-125.417920181) = -1.22277676266e-17, r=1 main: -169.498, 13.5149, s(-125.417920181) = 2.87991202066e-20, r=1 main: -38.4868, -95.0093, s(-125.417920181) = 2.87991202066e-20, r=1 main: -81.5745, -137.031, s(-125.417920181) = -8.49743452686e-18, r=1 main: 68.7282, 82.1745, s(-125.417920181) = 5.92584249899e-18, r=1 main: -5.49065, -35.9342, s(-125.417920181) = 2.87991202066e-20, r=1 main: -46.6596, -220.311, s(-125.417920181) = -1.64154985178e-18, r=1 main: -221.942, -182.517, s(-125.417920181) = -1.13502414932e-19, r=1 main: -214.685, 26.9286, s(-125.417920181) = -1.13502414932e-19, r=1 main: 99.5987, -232.548, s(-125.417920181) = 2.87991202066e-20, r=1 main: -162.686, -44.5434, s(-125.417920181) = -1.13502414932e-19, r=1 main: -88.9695, -77.5122, s(-125.417920181) = 1.30781887056e-18, r=1 main: -163.997, 89.841, s(-125.417920181) = -6.43745039913e-18, r=1 main: 14.1433, -76.5951, s(-125.417920181) = 2.87991202066e-20, r=1 main: 59.5242, -123.805, s(-125.417920181) = -4.23516473627e-20, r=1 main: -141.63, -167.431, s(-125.417920181) = -4.23516473627e-20, r=1 main: 69.7847, 28.2524, s(-125.417920181) = -4.23516473627e-20, r=1 main: 81.0914, -247.121, s(-125.417920181) = -2.31578807779e-18, r=1 main: w(t) is 'P_water_sat_Wexler' main: -273.15, 100, w(-125.418952599) = -2.69247289386e-19, r=1 main: -142.961, 59.0403, w(-125.418952599) = -4.11348646816e-19, r=1 main: -102.359, -265.678, w(-125.418952599) = -7.61196812136e-17, r=1 main: 89.0421, -184.236, w(-125.418952599) = -2.05504751383e-20, r=1 main: -71.7384, -173.644, w(-125.418952599) = -1.37059029852e-18, r=1 main: 10.5368, -133.142, w(-125.418952599) = -1.62664240278e-19, r=1 main: -23.9889, -81.8974, w(-125.418952599) = 5.0500203577e-20, r=1 main: -258.493, -74.7811, w(-125.418952599) = -2.05504751383e-20, r=1 main: 74.1817, 74.5643, w(-125.418952599) = 2.07554174068e-18, r=1 main: -167.066, -4.12663, w(-125.418952599) = -1.62664240278e-19, r=1 main: -141.037, -34.3419, w(-125.418952599) = -7.31096139779e-19, r=1 main: -108.925, -211.217, w(-125.418952599) = 5.0500203577e-20, r=1 main: -149.885, 36.2664, w(-125.418952599) = 5.0500203577e-20, r=1 main: -142.413, 60.2119, w(-125.418952599) = -1.62664240278e-19, r=1 main: -53.4989, 83.7561, w(-125.418952599) = 2.09049253395e-17, r=1 main: 47.2651, -27.877, w(-125.418952599) = -2.05504751383e-20, r=1 main: 71.6293, -109.128, w(-125.418952599) = -1.62664240278e-19, r=1 main: 30.8803, -124.473, w(-125.418952599) = 9.97751739957e-17, r=1 main: 66.7792, -17.8339, w(-125.418952599) = -1.62664240278e-19, r=1 main: -192.615, -93.1086, w(-125.418952599) = -1.9819578349e-19, r=1 main: 55.6183, -218.051, w(-125.418952599) = -4.60355789388e-18, r=1 ``` Best regards, Lyuka

 Possibly Related Threads... Thread Author Replies Views Last Post A fast Bernoulli Number method for the HP Prime Namir 16 2,211 11-22-2013, 04:46 PM Last Post: Namir A brand new calculator benchmark: "middle square method seed test" Pier Aiello 25 3,051 09-13-2013, 01:58 PM Last Post: Pier Aiello Downhill Simplex Method (Nelder Mead) for WP-34S Marcel Samek 0 397 07-07-2013, 03:05 AM Last Post: Marcel Samek New Quadratic Lagrangian Interpolation Method Namir 2 637 07-20-2012, 04:32 PM Last Post: Namir New variant for the Romberg Integration Method Namir 8 1,188 04-18-2012, 07:47 AM Last Post: Nick_S question for lyuka and the Ostrowski method - bumped! hugh steers 7 1,067 08-27-2011, 11:47 AM Last Post: Namir Lyuka and the Ostrowski method's for Root Seeking Namir 1 444 08-20-2011, 11:03 AM Last Post: Lyuka An Attempt to Challenge the Bisection method for root findin Namir 6 930 08-14-2011, 05:57 PM Last Post: Paul Dale Further to an earlier posting, Method to the madness! Geoff Quickfall 6 784 03-16-2010, 09:38 PM Last Post: Geoff Quickfall New method of typing numbers? Saile (Brazil) 9 1,057 02-22-2009, 03:58 PM Last Post: Maximilian Hohmann

Forum Jump: 