question for lyuka and the Ostrowski method - bumped! - Printable Version +- HP Forums (https://archived.hpcalc.org/museumforum) +-- Forum: HP Museum Forums (https://archived.hpcalc.org/museumforum/forum-1.html) +--- Forum: Old HP Forum Archives (https://archived.hpcalc.org/museumforum/forum-2.html) +--- Thread: question for lyuka and the Ostrowski method - bumped! (/thread-191841.html) question for lyuka and the Ostrowski method - bumped! - hugh steers - 08-25-2011 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. Re: question for lyuka and the Ostrowski method - bumped! - Namir - 08-25-2011 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. Re: question for lyuka and the Ostrowski method - bumped! - Lyuka - 08-25-2011 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. Re: question for lyuka and the Ostrowski method - bumped! - Namir - 08-26-2011 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. Re: question for lyuka and the Ostrowski method - bumped! - C.Ret - 08-26-2011 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 ``` Re: question for lyuka and the Ostrowski method - bumped! - Namir - 08-26-2011 Your suggestion should work too. Question about the last equation? Should be SIGN of X sub-A, without the ABS function? Namir Re: question for lyuka and the Ostrowski method - bumped! - Namir - 08-27-2011 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. Re: question for lyuka and the Ostrowski method - bumped! - Lyuka - 08-28-2011 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.