question for lyuka and the Ostrowski method - bumped!



#2

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.


#3

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.


#4

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.

#5

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.


#6

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.


#7

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


#8

Your suggestion should work too.

Question about the last equation? Should be SIGN of X sub-A, without the ABS function?

Namir

#9

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 1,645 11-22-2013, 04:46 PM
Last Post: Namir
  A brand new calculator benchmark: "middle square method seed test" Pier Aiello 25 2,185 09-13-2013, 01:58 PM
Last Post: Pier Aiello
  Downhill Simplex Method (Nelder Mead) for WP-34S Marcel Samek 0 269 07-07-2013, 03:05 AM
Last Post: Marcel Samek
  New Quadratic Lagrangian Interpolation Method Namir 2 423 07-20-2012, 04:32 PM
Last Post: Namir
  New variant for the Romberg Integration Method Namir 8 852 04-18-2012, 07:47 AM
Last Post: Nick_S
  question for lyuka and the Ostrowski method hugh steers 14 1,059 08-22-2011, 04:40 PM
Last Post: Lyuka
  Lyuka and the Ostrowski method's for Root Seeking Namir 1 306 08-20-2011, 11:03 AM
Last Post: Lyuka
  An Attempt to Challenge the Bisection method for root findin Namir 6 684 08-14-2011, 05:57 PM
Last Post: Paul Dale
  Further to an earlier posting, Method to the madness! Geoff Quickfall 6 574 03-16-2010, 09:38 PM
Last Post: Geoff Quickfall
  New method of typing numbers? Saile (Brazil) 9 784 02-22-2009, 03:58 PM
Last Post: Maximilian Hohmann

Forum Jump: