Van Wijngaarden-Dekker-Brent Method


sham_IIT
01-05-2009, 11:59 PM
Can I use Brent's method for one-dimensional function whose roots lies between [0 1] ?? Considering a=0, b=1, fa > 0.0 and fb > 0.0.
When I tried to use Brent's algorithm for such function, root is obtained quite faster than bisection method but while execution it doesnot accept interpolation means not a single time it interpolates and enter into the check
if (2.0p < (min1 < min2 ? min1 : min2))
{
e=d;
d=p/q;
}
whereas it always goes to the else part
else
{
d=xm;
e=d;
}

can anybody clear my doubts:
1. Is is correct to use Brent's algorithm for such functions ??
2. Why its taking less computational time than the time taken to find the root using only bisection method ??

davekw7x
01-07-2009, 11:17 AM
Can I...

What function and interval are you using?

The bisection method and Brent's method work for a certain class of functions if you give the algorithm a bracketed interval (the sign of the function at the end points must be different). I note, in particular, that the NR functions rtbis() and zbrent() exit with an error if the function value at both endpoints is greater than zero (or both less than zero).


Assuming that there was a misprint when you said "fa > 0.0 and fb > 0.0," can you give us some other information?

How did you measure the timing?

How did you determine what branch was being taken in zbrent()?

Did you modify the distribution functions? Did you use a debugger to set a breakpoint and/or single-step through the loop? Or what?

Also, you might tell us what compiler and operating system you are using. Sometimes this makes a difference to people who would like to help.

Maybe you can post your test program so that we might discover some specifics rather than dealing in vague generalities.


Regards,

Dave

sham_IIT
01-13-2009, 03:50 AM
Function I have is one-dimensional cost function.
optimum value of the variable should range between 0 to 1.
I TRIED TO OBTAIN THE OPTIMUM VALUE OF THE VARIABLE USING 'zbrent'.................

Initial values I took x1 = 0 and x2 = 1
I used xacc = 1.0e-8,
zbrent function i used as :
/////////////////////////////////////////////
float zbrent(const float x1,const float x2,const float tol)
{
const int ITMAX= 500;
int iter;
float a=x1,b=x2,c=x2,d,e,min1,min2;
float fa,fb,fc,p,q,r,s,tol1,xm;
func(a,fa);
func(b,fb);

//if ((fa>0.0 && fb>0.0) || (fa<0.0 && fb<0.0))
//cout<<"root must be bracketed";

fc=fb;
for (iter=0;iter<ITMAX;iter++)
{
if ((fb>0.0 && fc>0.0) || (fb<0.0 && fc<0.0))
{
c=a;
fc=fa;
e=d=b-a;
}
if (fabs(fc) < fabs(fb))
{
a=b;
b=c;
c=a;
fa=fb;
fb=fc;
fc=fa;
}
tol1=2.0*EPS*fabs(b)+0.5*tol;
xm = 0.5*(c-b);
if (fabs(xm) <= tol1 || fb == 0.0)
return b;
if (fabs(e) >= tol1 && fabs(fa) > fabs(fb))
{
s=fb/fa;
if (a==c)
{
p=2.0*xm*s;
q=1.0-s;
}
else
{
q=fa/fc;
r=fb/fc;
p=s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0));
q=(q-1.0)*(r-1.0)*(s-1.0);
}
if (p>0.0)
q=-q;
p=fabs(p);
min1 = 3.0*xm*q-fabs(tol1*q);
min2=fabs(e*q);
if (2.0*p < (min1<min2?min1:min2))
{
e=d;
d=p/q;
}
else
{
d=xm;
e=d;
}
}
else
{
d=xm;
e=d;
}
a=b;
fa=fb;
if (fabs(d) > tol1)
b+=d;
else
b+=SIGN(tol1,xm);
func(b,fb);
}
cout<<"maximum no of iterations exceeded in zbrent";
return 0.0;
}
////////////////////////////////////////////////
***NOTE: the error for (fa>0.0 && fb>0.0) or (fa<0.0 && fb<0.0)
I have removed.....

Now if I call,
temp = zbrent(x1,x4,xacc);

while execution not a single time it interpolates and enter into the check
if (2.0p < (min1 < min2 ? min1 : min2))
{
e=d;
d=p/q;
}

whereas it always goes to the else part

else
{
d=xm;
e=d;
}
i.e. it bisects...
The answer it gives ranges between 0 to 1 and also better than given by goldensection.

As per my idea zbent works when there is a change in sign of the function between the interval [x1,x2] right?? Then by removing the error part (if fa>0 and fb>0) how come I'm getting the answer. Even if I change the set of data of my cost function still I'm getting the answer ranging between [0,1].......

Can you please give me some explanations regarding the above problem.

davekw7x
01-13-2009, 09:56 AM
one-dimensional cost functionWait a minute! Are you trying to find a zero of the function or the minimum value of the function? The zbrent() and bisection (rtbis()) functions are for finding real zeros of functions. You indicated that you are looking for a zero of a function.

...zbent works when there is a change in sign of the function between the interval [x1,x2] right?? Then by removing the error part (if fa>0 and fb>0) how come...Let me get this straight:
You are defeating the test and then feeding the method a function that would fail the test and wondering why it worked after all? Maybe you just got lucky. (Or did it work? I'm not sure what you are getting at here. Maybe I am missing the whole point.)

Anyhow...

If you have a continuous function that has the same sign at both ends of an interval, then the possibilities are:

No real roots on the interval.
An even number of real roots on the interval (taking multiplicity into account).


If a real root exists on such an interval (that is, an interval for which the sign of the function is the same at both endpoints), then bisection and/or Brent may or may not find one of the real roots.



Assuming that you are really looking for a zero of the function:
Why not use zbrak() to try to discover actual sub-intervals with endpoints for which the function takes on opposite signs and then try rtbis() and/or zbrent() find a root on such an interval?



On the other hand, if you are looking for a minimum value (not a zero):You could try mnbrak() then golden() and/or and brent().




Regards,

Dave

sham_IIT
01-14-2009, 11:57 PM
I m trying to minimize a function which is dependent on a single variable say 'g'. Its not zero of a function but rather I want to get the value of 'g' (whose value ranges from 1 to 0) for which the function is minimum. Ultimately its minimization of a function....

davekw7x
01-15-2009, 08:44 AM
...trying to minimize a function... ...
In your first post you indicated that you obtained a root of the function (which I took to mean a zero, not a minimum value). I think that the lines of code you showed in that first post were from the Numerical Recipes Version 2 zbrent() function, which is for finding a zero of a function, not a minimum value.

I would say that it is not reasonable to depend on zbrent() (rtbis() either) to find a zero of a function unless you can feed it an interval for which the function takes on different signs at the endpoints. I mean, it is not impossible for it to work (after you disable the test for the condition that makes it most likely to work), but I wouldn't depend on it.

There are other ways.

Chapter 10 of the Numerical Recipes text contains descriptions of several approaches to finding minimum values of functions.

Regards,

Dave

sham_IIT
01-15-2009, 10:49 PM
Thanks for your kind suggestions....
I can understand that 'zbrent' should be used if the function values at the end point are of opposite sign. But still was not understanding how come its giving me answer when I remove the check...Anyways I'll try to use other algorithms suggested by you.

Regards
Sham