The Museum of HP Calculators

HP Forum Archive 20

[ Return to Index | Top of Index ]

question for lyuka and the Ostrowski method - bumped!
Message #1 Posted by hugh steers on 25 Aug 2011, 7:54 p.m.

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!
Message #2 Posted by Namir on 25 Aug 2011, 8:43 p.m.,
in response to message #1 by hugh steers

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!
Message #3 Posted by Lyuka on 28 Aug 2011, 12:52 a.m.,
in response to message #2 by Namir

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.

      
Re: question for lyuka and the Ostrowski method - bumped!
Message #4 Posted by Lyuka on 25 Aug 2011, 10:01 p.m.,
in response to message #1 by hugh steers

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!
Message #5 Posted by Namir on 26 Aug 2011, 9:27 a.m.,
in response to message #4 by Lyuka

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!
Message #6 Posted by C.Ret on 26 Aug 2011, 11:49 a.m.,
in response to message #5 by Namir

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!
Message #7 Posted by Namir on 26 Aug 2011, 2:44 p.m.,
in response to message #6 by C.Ret

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!
Message #8 Posted by Namir on 27 Aug 2011, 11:47 a.m.,
in response to message #1 by hugh steers

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.

[ Return to Index | Top of Index ]

Go back to the main exhibit hall