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.

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.

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.

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.

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 x_A, x_C and x_B equidistant :

x_A = user input.

x_C = x_A+H

x_B = x_C + H = x_A + 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

Your suggestion should work too.

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

Namir

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.

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.