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

 ▼ hugh steers Senior Member Posts: 536 Threads: 56 Joined: Jul 2005 08-25-2011, 07:54 PM Hello Lyuka, Namir and any interested parties. Sorry for the delay (been painting). I manually bumped this thread as the original is way down the list. Lyuka, thanks for the much improved Ostrowski method code. I have been trying it out. Firstly, it gets good results in almost all cases, beating many other methods. Also, you have fixed the case for the water equation. i'll just cover this again for detail, ```here is the water vapor pressure from temperature: Value Vwat(Value t) { Value 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; } ``` which you were solving. i mention this again just to add some background and also to point out the way it is normally solved. The inverse is called Twat :-) and proceeds usually with a baked in method of newton using both sides: ```Value Twat(Value e) { Value td = 0.0001; Value ta = 1e-7; Value T = 100; int c = 0; if (e >= 0) for (;;) { Value dt = 2 * td * (Vwat(T) - e) / (Vwat(T + td) - Vwat(T - td)); T -= dt; ++c; if (fabs(dt) < ta) break; } printf("twat count = %d\n", c); return T; } ``` which gives good results for this special case. now to give some results of your latest code against others in my test: ```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 ost iterations = 19 solve vwat=1e-5 in [-273.150000, 100.000000] with ostRoot = -106.310045 root not bracketed solve exp(6x-x^4-1)=1 in [0.000000, 2.000000] with brentRoot = 0.000000 100 iterations solve exp(6x-x^4-1)=1 in [0.000000, 2.000000] with secantRoot = -1.#IND00 4 iterations solve exp(6x-x^4-1)=1 in [0.000000, 2.000000] with ridderRoot = 1.757772 ost iterations = 9 solve exp(6x-x^4-1)=1 in [0.000000, 2.000000] with ostRoot = 1.757772 root not bracketed solve sin(x)-x=1 in [0.000000, 6.000000] with brentRoot = 0.000000 6 iterations solve sin(x)-x=1 in [0.000000, 6.000000] with secantRoot = -1.934563 9 iterations solve sin(x)-x=1 in [0.000000, 6.000000] with ridderRoot = -1.934563 ost iterations = 9 solve sin(x)-x=1 in [0.000000, 6.000000] with ostRoot = -1.934563 1 iterations solve (x-1)^2=1 in [-2.000000, 2.000000] with brentRoot = 2.000000 0 iterations solve (x-1)^2=1 in [-2.000000, 2.000000] with secantRoot = 2.000000 0 iterations solve (x-1)^2=1 in [-2.000000, 2.000000] with ridderRoot = -1.#QNAN0 solve (x-1)^2=1 in [-2.000000, 2.000000] with ostRoot = 2.000000 root not bracketed solve tan(x)-1/x=0 in [1.000000, 4.000000] with brentRoot = 1.000000 13 iterations solve tan(x)-1/x=0 in [1.000000, 4.000000] with secantRoot = -3.425618 5 iterations solve tan(x)-1/x=0 in [1.000000, 4.000000] with ridderRoot = 3.425618 ost iterations = 6 solve tan(x)-1/x=0 in [1.000000, 4.000000] with ostRoot = 3.425618 29 iterations solve (x-1)^5=0 in [0.000000, 1.500000] with brentRoot = 1.000000 100 iterations solve (x-1)^5=0 in [0.000000, 1.500000] with secantRoot = 1.000000 23 iterations solve (x-1)^5=0 in [0.000000, 1.500000] with ridderRoot = 1.000000 ost iterations = 8 solve (x-1)^5=0 in [0.000000, 1.500000] with ostRoot = 1.020561 29 iterations solve x^2-10000x+1=0 in [-1.000000, 1.000000] with brentRoot = 0.000100 3 iterations solve x^2-10000x+1=0 in [-1.000000, 1.000000] with secantRoot = 0.000100 2 iterations solve x^2-10000x+1=0 in [-1.000000, 1.000000] with ridderRoot = 0.000100 ost iterations = 2 solve x^2-10000x+1=0 in [-1.000000, 1.000000] with ostRoot = 0.000100 ``` you can see the improved Ostrowski method gets good results, but often with a few more iterations than ridders. ridders is fooled by my (x-1)^2-1 quadratic when i gave bad bounds. bounds are generally a problem with black box solvers since how is one to know good bounds to start? anyhow, Ostrowski struggles with a very flat curve at (x-1)^5. the other methods get 1.00000 but Ostrowski stops at 1.020561, which is not so good. i havent tried Namir's 4th order Ostrowski method. would that be better at this? thanks for the input, -- hugh. ▼ Namir Posting Freak Posts: 2,247 Threads: 200 Joined: Jun 2005 08-25-2011, 08:43 PM I got somewhat different results for (x-1)^5. When the initial guess is 1.5, tolerance for the root is 1e-7, and tolerance for the function value is 1e-15, I get: ```1. Newton 5635 iterations and 11271 function calls to get root = 1.00190837280905. 2. Halley 4763 iterations and 14290 function calls to get root = 1.00176867211809. 3. Simple Ostrowski 18 iterations and 22 function calls to get root = 0.999148573. 4. Ostrowski 4th Order (Grau and Diaz-Barrero) 71 iterations and 76 function calls to get root = 1.00048060997239. 5. Ostrowski (Luyka) 28 iterations and 31 function calls to get root = 1.00096671. ``` I modified Luyka's implementation to take one initial guess and calculate the second and third one based on that one initial guess--this way it is on par with the other methods. All the above methods are close to the solution of 1 by about 0.0001. The equation (x-1)^5 seems to pose a problem for methods that use first and second derivatives. Namir Edited: 25 Aug 2011, 8:44 p.m. ▼ Lyuka Member Posts: 169 Threads: 12 Joined: Aug 2007 08-28-2011, 12:52 AM Hi, I modified the code (Rev.3.01) to cope with the very flat region near the multiple root. When the solver encountered to that region, it will temporary change the strategy to use extrapolation. The behavior of the program for (x-1)^5 == 0 is as follows. (epsilon = 1e-50) Quote: ```main: (x - 1)^5 == 0 ost:c=0.75, f=-9.765625000000e-04, i=0, p=0, xp=0 ost:c=0.761904761905, f=-7.651622719419e-04, i=1, p=0, xp=0 ost:c=0.803189200351, f=-2.952872052788e-04, i=2, p=0, xp=0 ost:c=0.850693877302, f=-7.419729755354e-05, i=3, p=0, xp=0 ost:c=0.878293985111, f=-2.670300955364e-05, i=4, p=0, xp=0 ost:c=0.902610176346, f=-8.761286801440e-06, i=5, p=0, xp=0 ost:c=0.92234070912, f=-2.824666003999e-06, i=6, p=0, xp=0 ost:c=0.937768603186, f=-9.333569060579e-07, i=7, p=0, xp=0 ost:c=0.950216234015, f=-3.058008812555e-07, i=8, p=0, xp=0 ost:c=0.960176299828, f=-1.001631669682e-07, i=9, p=0, xp=0 ost:c=0.968133427318, f=-3.286069794729e-08, i=10, p=0, xp=0 ost:c=0.974504259003, f=-1.077303803554e-08, i=11, p=0, xp=0 ost:c=0.979601206998, f=-3.532013499117e-09, i=12, p=0, xp=0 ost:c=0.983678883914, f=-1.158108552378e-09, i=13, p=0, xp=0 ost:c=0.986941586119, f=-3.797100974770e-10, i=14, p=0, xp=0 ost:c=0.989552033388, f=-1.244969983086e-10, i=15, p=0, xp=0 ost:c=0.999942690022, f=-6.182316366438e-22, i=16, p=0, xp=3 ost:c=0.999942690022, f=-6.182316273200e-22, i=17, p=0, xp=2 ost:c=0.999954164675, f=-2.023026512788e-22, i=18, p=0, xp=1 ost:c=0.99996505284, f=-5.212660555528e-23, i=19, p=0, xp=0 ost:c=0.999971418708, f=-1.907257352170e-23, i=20, p=0, xp=0 ost:c=0.999977183259, f=-6.183980300985e-24, i=21, p=0, xp=0 ost:c=0.999981796202, f=-1.998987508836e-24, i=22, p=0, xp=0 ost:c=0.999985410203, f=-6.610681203429e-25, i=23, p=0, xp=0 ost:c=1.00000024196, f=8.292946741195e-34, i=24, p=0, xp=3 ost:c=1.00000024196, f=8.292935820281e-34, i=25, p=0, xp=2 ost:c=1.00000019372, f=2.728495698200e-34, i=26, p=0, xp=1 ost:c=1.0000001476, f=7.005533461945e-35, i=27, p=0, xp=0 ost:c=1.00000012073, f=2.564963532939e-35, i=28, p=0, xp=0 ost:c=1.00000009639, f=8.318656491505e-36, i=29, p=0, xp=0 ost:c=1.00000007689, f=2.688370094878e-36, i=30, p=0, xp=0 ost:c=1.00000006163, f=8.891095035563e-37, i=31, p=0, xp=0 ost:c=0.999999998879, f=-1.771758485393e-45, i=32, p=0, xp=3 ost:c=0.999999998879, f=-1.771754099344e-45, i=33, p=0, xp=2 ost:c=0.999999999102, f=-5.831575405038e-46, i=34, p=0, xp=1 ost:c=0.999999999301, f=-1.674142801941e-46, i=35, p=0, xp=0 ost:c=0.999999999434, f=-5.832025164875e-47, i=36, p=0, xp=0 ost:c=0.999999999547, f=-1.901053635115e-47, i=37, p=0, xp=0 ost:c=0.999999999638, f=-6.183968655566e-48, i=38, p=0, xp=0 ost:c=0.999999999995, f=-3.345012594248e-57, i=39, p=0, xp=3 ost:c=0.999999999995, f=-3.345012594248e-57, i=40, p=0, xp=2 ost:c=0.999999999238, f=-2.560955415462e-46, i=41, p=1, xp=2 ost:c=0.999999999661, f=-4.464993380546e-48, i=42, p=1, xp=1 ost:c=0.999999999679, f=-3.407016434088e-48, i=43, p=1, xp=0 ost:c=0.999999999747, f=-1.038342380366e-48, i=44, p=1, xp=0 ost:c=0.999999999803, f=-2.969976269079e-49, i=45, p=1, xp=0 ost:c=0.99999999984, f=-1.049595464097e-49, i=46, p=1, xp=0 ost:c=0.999999999872, f=-3.400293066647e-50, i=47, p=1, xp=0 ost:c=0.999999999898, f=-1.106478409920e-50, i=48, p=1, xp=0 ost:c=0.999999999918, f=-3.648288759905e-51, i=49, p=1, xp=0 ``` Regards, Lyuka /* Rev.3.01 : typo in the source code corrected. */ Edited: 29 Aug 2011, 4:14 p.m. Lyuka Member Posts: 169 Threads: 12 Joined: Aug 2007 08-25-2011, 10:01 PM Hi hugh steers, Namir and all who interested. I've modified the code (Rev.2.94) to struggle to escape from the flat region or extremes, this morning. Now, it can solve the water saturation pressure equation started at the very flat region e.g. [-273.15, -273.0]. Please try the latest code (Rev.2.94). Quote: ```/*s(x) : Water saturation pressures (Sonntag)*/ double s(x) double x; { double t; double p = 1e-5; t = x + 273.15; if (t < 0) return -p; return exp(((16.635794 + 2.433502 * log(t) + 0.00001673952 * t * t) - 0.02711193 * t) - 6096.9385 / t) - p; } main: s(t) is 'P_water_sat_Sonntag' ost:a=-273.15, b=-273.1, c=-273.125, f=-1.000000000000e-05, i=0, p=0 ost:a=-273.15, b=-273.225, c=-273.206078387, f=-1.000000000000e-05, i=1, p=1 ost:a=-273.15, b=-273.0375, c=-273.064372309, f=-1.000000000000e-05, i=2, p=2 ost:a=-273.15, b=-273.31875, c=-273.295157305, f=-1.000000000000e-05, i=3, p=3 ost:a=-273.15, b=-272.896875, c=-273.090425146, f=-1.000000000000e-05, i=4, p=4 ost:a=-273.15, b=-273.5296875, c=-273.285100373, f=-1.000000000000e-05, i=5, p=5 ost:a=-273.15, b=-272.58046875, c=-272.731564881, f=-1.000000000000e-05, i=6, p=6 ost:a=-273.15, b=-274.004296875, c=-273.411032864, f=-1.000000000000e-05, i=7, p=7 ost:a=-273.15, b=-271.868554688, c=-272.448762741, f=-1.000000000000e-05, i=8, p=8 ost:a=-273.15, b=-275.072167969, c=-274.073068154, f=-1.000000000000e-05, i=9, p=9 ost:a=-273.15, b=-270.266748047, c=-271.383252597, f=-1.000000000000e-05, i=10, p=10 ost:a=-273.15, b=-277.47487793, c=-274.800747146, f=-1.000000000000e-05, i=11, p=11 ost:a=-273.15, b=-266.662683105, c=-269.830272597, f=-1.000000000000e-05, i=12, p=12 ost:a=-273.15, b=-282.880975342, c=-281.866044434, f=-1.000000000000e-05, i=13, p=13 ost:a=-273.15, b=-258.553536987, c=-260.536164572, f=-1.000000000000e-05, i=14, p=14 ost:a=-273.15, b=-295.044694519, c=-286.697293227, f=-1.000000000000e-05, i=15, p=15 ost:a=-273.15, b=-240.307958221, c=-250.48456364, f=-1.000000000000e-05, i=16, p=16 ost:a=-273.15, b=-322.413062668, c=-282.332745028, f=-1.000000000000e-05, i=17, p=17 ost:a=-273.15, b=-199.255405998, c=-229.286334046, f=-1.000000000000e-05, i=18, p=18 ost:a=-273.15, b=-383.991891003, c=-281.658556097, f=-1.000000000000e-05, i=19, p=19 ost:a=-273.15, b=-106.887163496, c=-227.42340973, f=-1.000000000000e-05, i=20, p=20 ost:a=-273.15, b=-106.887163496, c=-242.794125415, f=-1.000000000000e-05, i=21, p=21 ost:a=-273.15, b=-106.887163496, c=-145.766953091, f=-9.999848970363e-06, i=22, p=22 ost:a=-106.887163496, b=-145.766953091, c=-106.88706331, f=-1.155857857677e-06, i=23, p=22 ost:a=-145.766953091, b=-106.88706331, c=-106.338424292, f=-6.002266450796e-08, i=24, p=22 ost:a=-106.88706331, b=-106.338424292, c=-106.311453476, f=-2.988926674806e-09, i=25, p=22 ost:a=-106.338424292, b=-106.311453476, c=-106.310044116, f=-1.842629348946e-13, i=26, p=22 ost:a=-106.311453476, b=-106.310044116, c=-106.310044029, f=1.496866036476e-20, i=27, p=22 main: -273.15, -273.1, s(-106.310044029) = 1.49686603648e-20, r=1 ``` And also, it can find a root of (x-1)^5 in 27 iterations. Quote: ```/* f6(x) = (x-1)^5 */ ost:a=0, b=1.5, c=0.75, f=-9.765625000000e-04, i=0, p=0 ost:a=1.5, b=0.75, c=0.761904761905, f=-7.651622719419e-04, i=1, p=0 ost:a=0.75, b=0.761904761905, c=0.803189200351, f=-2.952872052788e-04, i=2, p=0 ost:a=0.761904761905, b=0.803189200351, c=0.850693877302, f=-7.419729755354e-05, i=3, p=0 ost:a=0.803189200351, b=0.850693877302, c=0.878293985111, f=-2.670300955364e-05, i=4, p=0 ost:a=0.850693877302, b=0.878293985111, c=0.902610176346, f=-8.761286801440e-06, i=5, p=0 ost:a=0.878293985111, b=0.902610176346, c=0.92234070912, f=-2.824666003999e-06, i=6, p=0 ost:a=0.902610176346, b=0.92234070912, c=0.937768603186, f=-9.333569060579e-07, i=7, p=0 ost:a=0.92234070912, b=0.937768603186, c=0.950216234015, f=-3.058008812555e-07, i=8, p=0 ost:a=0.937768603186, b=0.950216234015, c=0.960176299828, f=-1.001631669682e-07, i=9, p=0 ost:a=0.950216234015, b=0.960176299828, c=0.968133427318, f=-3.286069794729e-08, i=10, p=0 ost:a=0.960176299828, b=0.968133427318, c=0.974504259003, f=-1.077303803554e-08, i=11, p=0 ost:a=0.968133427318, b=0.974504259003, c=0.979601206998, f=-3.532013499117e-09, i=12, p=0 ost:a=0.974504259003, b=0.979601206998, c=0.983678883914, f=-1.158108552378e-09, i=13, p=0 ost:a=0.979601206998, b=0.983678883914, c=0.986941586119, f=-3.797100974770e-10, i=14, p=0 ost:a=0.983678883914, b=0.986941586119, c=0.989552033388, f=-1.244969983086e-10, i=15, p=0 ost:a=0.986941586119, b=0.989552033388, c=0.991640629072, f=-4.081952835848e-11, i=16, p=0 ost:a=0.989552033388, b=0.991640629072, c=0.993311710543, f=-1.338367252555e-11, i=17, p=0 ost:a=0.991640629072, b=0.993311710543, c=0.994648733138, f=-4.388165984146e-12, i=18, p=0 ost:a=0.993311710543, b=0.994648733138, c=0.995718477996, f=-1.438768435661e-12, i=19, p=0 ost:a=0.994648733138, b=0.995718477996, c=0.996574375736, f=-4.717355769657e-13, i=20, p=0 ost:a=0.995718477996, b=0.996574375736, c=0.99725917517, f=-1.546701119566e-13, i=21, p=0 ost:a=0.996574375736, b=0.99725917517, c=0.997807079768, f=-5.071240091293e-14, i=22, p=0 ost:a=0.99725917517, b=0.997807079768, c=0.998245455502, f=-1.662730778483e-14, i=23, p=0 ost:a=0.997807079768, b=0.998245455502, c=0.998596197729, f=-5.451671807322e-15, i=24, p=0 ost:a=0.998245455502, b=0.998596197729, c=0.99887682483, f=-1.787464685796e-15, i=25, p=0 ost:a=0.998596197729, b=0.99887682483, c=0.999101353168, f=-5.860642591876e-16, i=26, p=0 ost:a=0.99887682483, b=0.999101353168, c=0.999280997168, f=-1.921555815013e-16, i=27, p=0 main:f6: 0, 1.5, s(0.999280997168) = -1.92155581501e-16, r=1 ``` Regards, Lyuka Edited: 25 Aug 2011, 11:59 p.m. ▼ Namir Posting Freak Posts: 2,247 Threads: 200 Joined: Jun 2005 08-26-2011, 09:27 AM I changed your implementation so that it uses the following pseudo-code to use one guess, A, and then generate the other two guesses, B and C: ```Given the guess A for root: Shift = 0.1 ' or any other value near 0.1 ' obtain guess B from the value of A If A > 0 Then B = (1 + Shift) * A ElseIf A < 0 Then B = A + Shift * Abs(A) Else B = Shift End If FA = Fx(A) ' evaluate function at A FB = Fx(B) ' evaluate function at B C = (A + B) / 2 ' calculate C as mid-value between A and B FC = Fx(C) ' evaluate function at C ``` The above modification seems to work well for all the cases and initial guesses I used. Namir Edited: 26 Aug 2011, 10:38 a.m. ▼ C.Ret Senior Member Posts: 260 Threads: 0 Joined: Oct 2008 08-26-2011, 11:49 AM Hi, May I suggest a small modification to your last implementation which appears to me quite clever since it is using the Shift parameter to automatically create the second "user guess";It is a good idea. I we may discuss the relation between this shift amplitude and the system accuracy eps used as zero approximation. It may be an implicit relation between these two parameters and the resulting number of iterations. But, you can go one step forward, your shift parameter may be used for the computation of the third point as well. For clarity, let H = Shift / 2 . Now we have the three initial values xA, xC and xB equidistant : xA = user input. xC = xA+H xB = xC + H = xA + 2.H Please account for my following modification of your small modification from the original implementation: The algorithm may now be : ```INPUT “System accuracy (zero approximation)”; eps INPUT “User guess (first root approximation)”; XA INPUT “Initial step (or initial scale in % )”; shift LET H = shift * ( ABS(XA) + 1 – SIGN(ABS(XA)) ) / 2 ‘ where SIGN(x)returns -1, 0 or +1 respectively for x negative, null or positive. ‘ where ABS returns absolute value of its argument x. LET XC = XA + 1*H LET XB = XA + 2*H LET FA = FN Fx(XA) ' evaluate function at A LET FB = FN Fx(XB) ' evaluate function at B LET FC = FN Fx(XC) ' evaluate function at C ... ``` Note, that in the code H may be negative, depending on how you are considering intervals in negative parts. The following computation of H also works, but negative user's guess have to be greater (on the right) than expected negative root, when positive guess have to be less (on the left) than expected positive root. ```LET H = shift * ( XA + 1 – SIGN(ABS(XA)) ) / 2 ``` ▼ Namir Posting Freak Posts: 2,247 Threads: 200 Joined: Jun 2005 08-26-2011, 02:44 PM Your suggestion should work too. Question about the last equation? Should be SIGN of X sub-A, without the ABS function? Namir Namir Posting Freak Posts: 2,247 Threads: 200 Joined: Jun 2005 08-27-2011, 11:47 AM I just wanted to comment about the methods that enhance Ostrowski's method (itself a refinement to Newton's method). These methods calculate 2, 3, and even 4 intermediate refinements in each step. You can think of these algorithms as having additional sub-steps. That is why counting function calls seems more relevant than counting iterations (of course it's good to know both counts). It seems that these methods are vulnerable to equations with non-smooth curves (such as (x-1)^5). Namir Edited: 27 Aug 2011, 11:48 a.m.

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

Forum Jump: 