The Museum of HP Calculators

HP Forum Archive 20

[ Return to Index | Top of Index ]

question for lyuka and the Ostrowski method
Message #1 Posted by hugh steers on 16 Aug 2011, 7:13 p.m.

hi,

wanted to try out this method. took the code here,

http://www.finetune.jp/~lyuka/technote/fract/ost.html

and tried one of my tests. this is smooth test from the temperature industry.

double f(double t)
{
    /* example from pressure */
    double 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 - 1e-5; // value = 1e-5
}

not very pretty i know, but it should have a real world solution to a value of 1e-5 (hardcoded) between -275.15 and 100, say.

i do this: 

int main()
{
    double          a, b, root;
    int             r;

    printf("vwat\n");
    r= ost(a=-271.15, b=100.0, f, EPSILON, ITERATION, &root);
    printf("main: %g, %g, f(%.18g) = %g, %d\n", a, b, root, f(root), r);
    return 0;
}

and it gives me -256.

here's the results from my other methods,

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
2 iterations
solve vwat=1e-5 in [-273.150000, 100.000000] with ostRoot = -256.086474

secant blows, but brent and ridder are fine. ost is the -256 as mentioned.

is there a way to fix for this, otherwise the method has failed.

thanks for any info, -- hugh.

      
Re: question for lyuka and the Ostrowski method
Message #2 Posted by Namir on 16 Aug 2011, 10:12 p.m.,
in response to message #1 by hugh steers

Hugh,

Thanks for the nice link. Interestingly the author uses my favorite test function f(x) = e^x - 3*x^2. Small world! :-)

I am going to look into the code and write a version for Excel VBA.

Namir

            
Re: question for lyuka and the Ostrowski method
Message #3 Posted by Lyuka on 16 Aug 2011, 11:44 p.m.,
in response to message #2 by Namir

Hi, Namir

It's not a small world, I've chosen the function because you often use the function as an example :-)

Best regards,
Lyuka

                  
Re: question for lyuka and the Ostrowski method
Message #4 Posted by Namir on 17 Aug 2011, 1:11 a.m.,
in response to message #3 by Lyuka

Hi Lyuka,

Thank you for a wonderful web page and links. The C and 42S source codes (plus the .raw file for the 42s) are very very nice.

Namir

Edited: 17 Aug 2011, 1:11 a.m.

            
Re: question for lyuka and the Ostrowski method
Message #5 Posted by Namir on 17 Aug 2011, 7:58 a.m.,
in response to message #2 by Namir

Lyuka's article about the Ostrowski algorithm comes in an interesting time. For the last week, I have been playing with basic root-seeking algorithms trying to use random numbers as a tool to tweak and improve these algorithms. Here is the pseudo-code for a method that is interesting. It is based on Bisection and uses random numbers:

Given: 

1. Function Fx(x)=0
2. Function Rand() that generates uniform random numbers in interval (0,1).
3. Interval {A,B] containing the root of function Fx(x)=0, such that A < B.
4. The Factor value used to tweak root guesses. Optimum value is around 0.75.  
5. Toler is smallest size of interval [A.B] used to obtain the root.

FA = Fx(A)
FB = Fx(B)
bFlag = False ' compensation mode flag
Do
  H = Abs(B - A) / 2
  ' do we need to compensate for the last iteration's failure
  ' to shrink the size of interval [A,B] by more than half?
  If bFlag Then
    ' compare the absolute function values at A and B to determine
    ' which end of the interval we will shift the new guess 
    ' towards, in hope of getting closer to the root
    If Abs(FA) < Abs(FB) Then
      ' closer to A
      X = A + H * Rand() * Factor
    Else
      ' closer to B
      X = B - H * Rand() * Factor
    End If
  Else
    ' get a random value inside interval {A, B] replacing the
    ' typical Bisection step of calculating X as (A+B)/2
    '
    ' The new value of X reflects success using random numbers
    ' if it shrinks the size of the interval [A,B] by more than
    ' half.
    X = A + (B - A) * Rand()
  End If
  FX = Fx(X)
  If FA * FX > 0 Then
    A = X
    FA = FX
  Else
    B = X
    FB = FX
  End If
  ' did the last iteration fail to shrink the size of
  ' the interval [A,B] by more than half?
  ' If yes, then random selection has failed and we need to
  ' compensate in the next iteration 
  bFlag = Abs(B - A) > H
Loop Until Abs(A - B) <= Toler 
Refined guess for root = X

I have been running solutions in batches of 10000 iterations. The above algorithm does better than the Bisection method between 20% to 50% of the time and saving between 1 to 8 iterations (in significant frequencies). The results depend on the interval [A,B] and where the root is located within the interval. If the root is near the initial boundaries of interval [A,B] then the above algorithm shows more success. Using the "Factor" is the latest aspect of tweaking the algorithm. I am still studying the effect of using Factor and its optimum values.

The above algorithm is not the greatest improvement on Bisection. Like I said before I am looking at using random numbers in tweaking basic algorithms I have used the approach with Newton and Halley's method. So far no success.

Namir

Edited: 17 Aug 2011, 8:15 a.m.

      
Re: question for lyuka and the Ostrowski method
Message #6 Posted by Lyuka on 16 Aug 2011, 11:40 p.m.,
in response to message #1 by hugh steers

Hi hugh,

Thank you for your interest on the page " Ostrowski's method " and "SOLVE alternative for the HP 42S " (necessary not to be OT ;-)

Take a look at the return code, it would be '3' for your initial guesses of {-271.15, 100.0}, which means

Quote:
3 : The search concentrated in a local flat region of the function. f(x_{n+2}) == f(x_{n+1})

It would occur mainly because of too wide set of initial guesses.
If you try that of {-271.15, 0} or {-100, 100},
the solve will converge to the correct answer of -106.310044029008282.

Quote:
main: -271.15, 100, f(-256.651643757074339) = -1e-05, 3
main: -271.15, 0, f(-106.310044029008282) = 1.52466e-20, 1
main: -100, 100, f(-106.310044029008282) = 1.52466e-20, 1

Best regards,
Lyuka

Edited: 17 Aug 2011, 1:38 a.m.

            
Re: question for lyuka and the Ostrowski method
Message #7 Posted by Namir on 17 Aug 2011, 1:48 a.m.,
in response to message #6 by Lyuka

I like your C code because it examines the function values for the initial guesses and is able to exit if any of these guesses are at or close to the root. This is what I call the "deluxe version" of an implementation where all scenarios are handled.

Namir

                  
Re: question for lyuka and the Ostrowski method
Message #8 Posted by Lyuka on 17 Aug 2011, 2:05 a.m.,
in response to message #7 by Namir

It's my pleasure, thank you for saying that :-)

Lyuka

                        
Re: question for lyuka and the Ostrowski method
Message #9 Posted by Namir on 17 Aug 2011, 8:56 a.m.,
in response to message #8 by Lyuka

Hi,

I was looking at the following part of your C code:

   if (0.0 == (f - e) * (f - d)) {
        c = (a + 3 * b) * 0.25;
        if (0.0 == (f = (*func)(c))) { *x = c; return 0; }
        if (0.0 == (f - e) * (f - d)) { *x = c; return 7; }
    }

How about something like:

   while (0.0 == (f - e) * (f - d)) {
        c = a + (b-a)*rand()/RAND_MAX;
        if (0.0 == (f = (*func)(c))) { *x = c; return 0; }
    }

The suggested code calculates c as a random value in the interval [a,b].

Namir

                              
Re: question for lyuka and the Ostrowski method
Message #10 Posted by Lyuka on 17 Aug 2011, 9:39 a.m.,
in response to message #9 by Namir

Hi,

It seems good an idea, however the iteration limit must be used.
If the condition kept false several times, it must be concerned that it would be occurring because of too narrow initial guesses, or not.
In such case, the third guess 'c' would better be put outside of the 'a' to 'b' region.
This kind of 'pivot' strategy is difficult to decide.

Best regards,
Lyuka

                                    
Re: question for lyuka and the Ostrowski method
Message #11 Posted by Namir on 17 Aug 2011, 10:37 a.m.,
in response to message #10 by Lyuka

You are right. After I posted the message I thought teh while loop perhaps needed a counter? Also the range for calculating c may need to be wider. You are not obligated to have a < c < b.

                                          
Re: question for lyuka and the Ostrowski method
Message #12 Posted by Lyuka on 18 Aug 2011, 4:53 p.m.,
in response to message #11 by Namir

Hi Namir,

Inspired by your hints of using random() function,
I modified the _ost() function to deal with the case f(x_n+2) == f(x_n+1)
by pivoting the value to the random position between two previous values.

Quote:
if (0 == (f - e)) {
    if (p++ > MAXPVT) { *x = c; return 3; }
    c = a + (b - a) * (1.0 + (double)random()) / M2E31;
    f = (*func)(c);
    continue;
}

It seems to work well, also for the case of hugh steer's temperature related function.

Best regards,
Lyuka

Edited: 18 Aug 2011, 5:10 p.m.

      
Re: question for lyuka and the Ostrowski method
Message #13 Posted by Namir on 16 Aug 2011, 11:43 p.m.,
in response to message #1 by hugh steers

Hugh,

I got to play with Ostrowski's method using f(x)=e^x-3*x^2 and comparing it with Newton's method and the faster Halley's method. Since Ostrowski's method uses two initial guesses (and calculates a third guess as the average of the first two), I used the average of these guesses for the Newton and Halley methods, to even up the score a bit. In this case, Ostrowski's method did better than Newton but not better than Halley.

If I use the first guess for Newton and Halley methods, then depending on that guess and the second guess (used only by Ostrowski's method) I get a few cases where Ostrowski's method outperforms the other two!

Thanks for the fantastic article (and thanks for Luyka for writing that web page). It is a true gem in numerical analysis and for a root-seeking junkie like me!!

Namir

Edited: 17 Aug 2011, 1:23 a.m.

            
Re: question for lyuka and the Ostrowski method
Message #14 Posted by hugh steers on 17 Aug 2011, 7:12 a.m.,
in response to message #13 by Namir

all credit goes to lyuka for the article. i followed his post earlier (seems to have disappeared) to investigate this method. yes, i am interested in root methods too. mostly i use ridders method and find it very stable. the code is also small. i wanted to try this method out and compare it.

to lyuka, thanks. it looks like the method finds my function flat. it is true that the curve has a flattened end. this has fooled other methods too. the annoying thing is that this particular example comes from a real world case.

now im thinking if there's a hybrid method involving ostrowski and other ideas.

-- hugh.

      
Re: question for lyuka and the Ostrowski method
Message #15 Posted by Lyuka on 22 Aug 2011, 4:40 p.m.,
in response to message #1 by hugh steers

Hi hugh steers,

I've been tried to make my implementation of Ostrowski's method to be able to solve the water saturation pressure related equations, as you shown above, for various initial guesses.

Now, the latest version of the program (Rev.2.8) seems to work well for such kind of equations with almost any initial guesses within the range of [-273.15 : 100], except when the two guesses are at flat region near -273.15.

I've tested the program using the water saturation pressures equations by Sonntag and Wexler as well, to solve it for the pressure to be 1e-5.

Quote:
double
s(x)
    double          x;
{
    double          t;
    double          p = 1e-5;

    t = x + 273.15;
    if (t < 0) return 0;
    return exp(((16.635794 + 2.433502 * log(t) + 1.673952e-5 * t * t) - 2.711193e-2 * t)
               - 6096.9385 / t) * 100.0 - p;
}

double
w(x)
    double          x;
{
    double          t;
    double          p = 1e-5;

    t = x + 273.15;
    if (t < 0) return 0;
    return exp(((18.87643845 + 2.858487 * log(t))
        + t * (-2.8354721e-2 + t * (1.7838301e-5 + t * (-8.4150417e-10 + t * 4.4412543e-13))))
        + (-6017.0128 -2991.2729 / t) / t) - p;
}

Here is the results of test run of the program.

Quote:
main: s(t) is 'P_water_sat_Sonntag'
main:   273.15,      100, s(-125.417920181) = 2.37609682364e-17, r=1
main: -125.986,   40.366, s(-125.417920181) = -4.23516473627e-20, r=1
main:  24.7879,  19.0635, s(-125.417920181) = 2.87991202066e-20, r=1
main: -199.434,  67.0312, s(-125.417920181) = -1.22277676266e-17, r=1
main: -169.498,  13.5149, s(-125.417920181) = 2.87991202066e-20, r=1
main: -38.4868, -95.0093, s(-125.417920181) = 2.87991202066e-20, r=1
main: -81.5745, -137.031, s(-125.417920181) = -8.49743452686e-18, r=1
main:  68.7282,  82.1745, s(-125.417920181) = 5.92584249899e-18, r=1
main: -5.49065, -35.9342, s(-125.417920181) = 2.87991202066e-20, r=1
main: -46.6596, -220.311, s(-125.417920181) = -1.64154985178e-18, r=1
main: -221.942, -182.517, s(-125.417920181) = -1.13502414932e-19, r=1
main: -214.685,  26.9286, s(-125.417920181) = -1.13502414932e-19, r=1
main:  99.5987, -232.548, s(-125.417920181) = 2.87991202066e-20, r=1
main: -162.686, -44.5434, s(-125.417920181) = -1.13502414932e-19, r=1
main: -88.9695, -77.5122, s(-125.417920181) = 1.30781887056e-18, r=1
main: -163.997,   89.841, s(-125.417920181) = -6.43745039913e-18, r=1
main:  14.1433, -76.5951, s(-125.417920181) = 2.87991202066e-20, r=1
main:  59.5242, -123.805, s(-125.417920181) = -4.23516473627e-20, r=1
main:  -141.63, -167.431, s(-125.417920181) = -4.23516473627e-20, r=1
main:  69.7847,  28.2524, s(-125.417920181) = -4.23516473627e-20, r=1
main:  81.0914, -247.121, s(-125.417920181) = -2.31578807779e-18, r=1

main: w(t) is 'P_water_sat_Wexler'
main:  -273.15,      100, w(-125.418952599) = -2.69247289386e-19, r=1
main: -142.961,  59.0403, w(-125.418952599) = -4.11348646816e-19, r=1
main: -102.359, -265.678, w(-125.418952599) = -7.61196812136e-17, r=1
main:  89.0421, -184.236, w(-125.418952599) = -2.05504751383e-20, r=1
main: -71.7384, -173.644, w(-125.418952599) = -1.37059029852e-18, r=1
main:  10.5368, -133.142, w(-125.418952599) = -1.62664240278e-19, r=1
main: -23.9889, -81.8974, w(-125.418952599) = 5.0500203577e-20, r=1
main: -258.493, -74.7811, w(-125.418952599) = -2.05504751383e-20, r=1
main:  74.1817,  74.5643, w(-125.418952599) = 2.07554174068e-18, r=1
main: -167.066, -4.12663, w(-125.418952599) = -1.62664240278e-19, r=1
main: -141.037, -34.3419, w(-125.418952599) = -7.31096139779e-19, r=1
main: -108.925, -211.217, w(-125.418952599) = 5.0500203577e-20, r=1
main: -149.885,  36.2664, w(-125.418952599) = 5.0500203577e-20, r=1
main: -142.413,  60.2119, w(-125.418952599) = -1.62664240278e-19, r=1
main: -53.4989,  83.7561, w(-125.418952599) = 2.09049253395e-17, r=1
main:  47.2651,  -27.877, w(-125.418952599) = -2.05504751383e-20, r=1
main:  71.6293, -109.128, w(-125.418952599) = -1.62664240278e-19, r=1
main:  30.8803, -124.473, w(-125.418952599) = 9.97751739957e-17, r=1
main:  66.7792, -17.8339, w(-125.418952599) = -1.62664240278e-19, r=1
main: -192.615, -93.1086, w(-125.418952599) = -1.9819578349e-19, r=1
main:  55.6183, -218.051, w(-125.418952599) = -4.60355789388e-18, r=1

Best regards,
Lyuka

[ Return to Index | Top of Index ]

Go back to the main exhibit hall