Posts: 536
Threads: 56
Joined: Jul 2005
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.
Posts: 2,247
Threads: 200
Joined: Jun 2005
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
Posts: 169
Threads: 12
Joined: Aug 2007
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.
Posts: 2,247
Threads: 200
Joined: Jun 2005
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.
Posts: 169
Threads: 12
Joined: Aug 2007
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
Posts: 2,247
Threads: 200
Joined: Jun 2005
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.
Posts: 2,247
Threads: 200
Joined: Jun 2005
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
Posts: 169
Threads: 12
Joined: Aug 2007
It's my pleasure, thank you for saying that :-)
Lyuka
Posts: 536
Threads: 56
Joined: Jul 2005
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.
Posts: 2,247
Threads: 200
Joined: Jun 2005
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.
Posts: 2,247
Threads: 200
Joined: Jun 2005
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
Posts: 169
Threads: 12
Joined: Aug 2007
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
Posts: 2,247
Threads: 200
Joined: Jun 2005
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.
Posts: 169
Threads: 12
Joined: Aug 2007
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.
Posts: 169
Threads: 12
Joined: Aug 2007
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
|