0

From textbooks I found that the tolerance for Golden Section Search method should be set to $\sqrt{\epsilon}$, where $\epsilon$ - is the machine epsilon. This can be derived from Taylor series. So, in other words, setting tolerance less then $\sqrt{\epsilon}$ is useless, since the computer wouldn't differ the values, whose differences are within that tolerance.

In my program, I used double format ( $\epsilon = 2.2*10^{-16}$), meaning for tolerance $\approx 10^{-8}$. However, when I tried to use tolerance levels like $10^{-11}$ and $10^{-10}$, I got two corresponding to those tolerances abscissas $x1$ and $x2$ for minimum of the same function. But |x1 - x2| turned out to be $10^{-12}$, which is not even close to machine epsilon. Could you please explain why is that? According to textbooks the difference should be on the order of machine epsilon. Thanks in advance. Here is my full program.

#include<iostream>

#include<iomanip>

#include<cmath> using namespace std;

double function(double x) { const double bohrRadius = 5.29e-11; const double C = 1 / (81 * sqrt(3 * M_PI)) * pow(1/bohrRadius , 2/3); return C * (27 - 18 * x + 2 * x * x) * exp(-x / 3); }

void bracket(double (*func)(double x), double x0, double step, double &xleft, double &xright) {

double f0, f1, f2; double x1, x2; double factor = 1.5; bool FLAG = false; f0 = func(x0); x1 = x0 + step; f1 = func(x1); if (f1 > f0) { step = - step; x1 = x0 + step; f1 = func(x1); if (f1 > f0) { xleft = x1; xright = x0 - step; return; } } do { step = factor * step; x2 = x1 + step; f2 = func(x2); if (f2 > f1) { xleft = x2; xright = x0; FLAG = true; return; } x0 = x1; x1 = x2; f0 = f1; f1 = f2; }while (!FLAG); }

double findExtremum(double (*func)(double x), double a, double b, double tol) {

double R = 0.618033989; double f1, f2; double x1, x2; x1 = b - R * (b - a); x2 = a + R * (b - a); f1 = func(x1); f2 = func(x2);

do {

if (f1 > f2) {
  a = x1;
  x1 = x2;
  f1 = f2;
  x2 = a + R * (b - a);
  f2 = func(x2);
}
else {

  b = x2;
  x2 = x1;
  f2 = f1;
  x1 = b - R * (b - a);
  f1 = func(x1);
}

}while ( abs(x1 - x2) > tol);

return (x1 + x2)/2; }

int main() { double Xmin1, Xmin2; double xleft = 0; double xright = 0; double step = 1; double initialGuess = 2; bracket(function, initialGuess, step, xleft, xright); if (xleft < xright) { Xmin1 = findExtremum(function, xleft, xright, 1e-10); Xmin2 = findExtremum(function, xleft, xright, 1e-11); double diff = Xmin1 - Xmin2; cout << diff << endl; } else { Xmin1 = findExtremum(function, xright, xleft, 1e-10); Xmin2 = findExtremum(function, xright, xleft, 1e-11); double diff = Xmin1 - Xmin2; cout << diff << endl; } }

atah1991
  • 3
  • 1
    Are you finding the minimum of a function or a root? It would be good to give the function and your results. It should be $10^{-16}$, not $\epsilon^{-16}$ and some other places. – Ross Millikan Mar 15 '15 at 23:06
  • Thanks for your comments. I am trying to find the minimum of the function – atah1991 Mar 15 '15 at 23:21

1 Answers1

0

If you are trying to find the minimum, the point isn't that the $x$ values will be too close, but that the $y$ values will be. Let $f(x)=1+(x-\pi)^2$ If we try to minimize this, the minimum will be at $(\pi,1)$, but for $x$ values between $\pi - \sqrt \epsilon$ and $\pi + \sqrt \epsilon$ the function value will be indistinguishable from $1$ In your problem, try taking $f(x1)$ and $f(x2)$ and see how close they are. This also gives a clue why we say "of order $\sqrt \epsilon$. If we had $f(x)=1+10^{-10}(x-\pi)^2$ the uncertainty in $x$ would be much larger. If we had $f(x)=(x-\pi)^2$ we could find the minimum to great accuracy.

Ross Millikan
  • 374,822
  • If this is true, then why the temination is on |x1 - x2| < $\sqrt{\epsilon}$, but not on the values of functions theirselve? – atah1991 Mar 15 '15 at 23:49
  • I tried what you said and the values are really very close ( within machine epsilon). Thank you, but the previous question still holds. I guess, we can terminate in loop using the values of functions and compare them with $\epsilon$ but not its square root. The only reason I see, that it is not cheap in terms of operations. I mean calling the function takes some time. Correct me if I am wrong or if I am missing something. – atah1991 Mar 15 '15 at 23:55
  • Because the heuristic is that when the $x$ values are closer than $\sqrt \epsilon,$ the $y$ values will be closer than $\epsilon$. The examples I gave show that this is not always true. You are right that we usually expect function calls to be expensive, but you have presumably already called it at $x1$ and $x2$ so know the value. The question is whether it is worth calling at the interior point or not. – Ross Millikan Mar 16 '15 at 00:23