question for lyuka and the Ostrowski method
#1

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.

#2

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

#3

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.

#4

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.

#5

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

#6

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.

#7

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

#8

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



Lyuka

#9

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.

#10

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.

#11

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

#12

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

#13

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.

#14

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.

#15

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



Possibly Related Threads…
Thread Author Replies Views Last Post
  A fast Bernoulli Number method for the HP Prime Namir 16 5,361 11-22-2013, 04:46 PM
Last Post: Namir
  A brand new calculator benchmark: "middle square method seed test" Pier Aiello 25 6,761 09-13-2013, 01:58 PM
Last Post: Pier Aiello
  Downhill Simplex Method (Nelder Mead) for WP-34S Marcel Samek 0 866 07-07-2013, 03:05 AM
Last Post: Marcel Samek
  New Quadratic Lagrangian Interpolation Method Namir 2 1,316 07-20-2012, 04:32 PM
Last Post: Namir
  New variant for the Romberg Integration Method Namir 8 2,440 04-18-2012, 07:47 AM
Last Post: Nick_S
  question for lyuka and the Ostrowski method - bumped! hugh steers 7 2,286 08-27-2011, 11:47 AM
Last Post: Namir
  Lyuka and the Ostrowski method's for Root Seeking Namir 1 932 08-20-2011, 11:03 AM
Last Post: Lyuka
  An Attempt to Challenge the Bisection method for root findin Namir 6 2,029 08-14-2011, 05:57 PM
Last Post: Paul Dale
  Further to an earlier posting, Method to the madness! Geoff Quickfall 6 1,720 03-16-2010, 09:38 PM
Last Post: Geoff Quickfall
  New method of typing numbers? Saile (Brazil) 9 2,229 02-22-2009, 03:58 PM
Last Post: Maximilian Hohmann

Forum Jump: