1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146
| #include <assert.h>
#include <math.h>
#include <float.h>
/**
* Estimates the min location.
*
* @param a left border | of the range
* @param b right border | the min is seeked
* @param f function under investigation
* @param tol acceptable tolerance
*/
double fminbr(double a, double b, double (*f)(double x), double tol)
{
double x, v, w; /* Abscissae, descr. see above */
double fx; /* f(x) */
double fv; /* f(v) */
double fw; /* f(w) */
const double r = (3. - sqrt(5.0)) / 2; /* Gold section ratio */
assert( tol > 0 && b > a );
v = a + r * (b - a);
fv = f(v); /* First step - always gold section*/
x = v;
w = v;
fx = fv;
fw = fv;
for (;;) /* Main iteration loop */
{
double range = b - a; /* Range over which the minimum */
/* is seeked for */
double middle_range = (a + b) / 2;
double tol_act = FLT_EPSILON*fabs(x) + tol/3; /* Actual tolerance */
double new_step; /* Step at this iteration */
if ( fabs(x - middle_range) + range / 2 <= 2 * tol_act )
{
return x; /* Acceptable approx. is found */
}
/* Obtain the gold section step */
new_step = r * ( x < middle_range ? b - x : a - x );
/* Decide if the interpolation can be tried */
if ( fabs(x - w) >= tol_act ) /* If x and w are distinct */
{
/* interpolatiom may be tried */
double p; /* Interpolation step is calcula-*/
double q; /* ted as p/q; division operation*/
/* is delayed until last moment */
double t;
t = (x-w) * (fx-fv);
q = (x-v) * (fx-fw);
p = (x-v)*q - (x-w)*t;
q = 2*(q-t);
if ( q>(double)0 ) /* q was calculated with the op-*/
{
p = -p; /* posite sign; make q positive */
}
else /* and assign possible minus to */
{
q = -q; /* p */
}
if ( fabs(p) < fabs(new_step*q) && /* If x+p/q falls in [a,b]*/
p > q*(a-x+2*tol_act) && /* not too close to a and */
p < q*(b-x-2*tol_act) ) /* b, and isn't too large */
{
new_step = p/q; /* it is accepted */
}
/* If p/q is too large then the */
/* gold section procedure can */
/* reduce [a,b] range to more */
/* extent */
}
if ( fabs(new_step) < tol_act ) /* Adjust the step to be not less*/
{
if ( new_step > (double)0 ) /* than tolerance */
{
new_step = tol_act;
}
else
{
new_step = -tol_act;
}
}
/* Obtain the next approximation to min */
{
/* and reduce the enveloping range */
double t = x + new_step; /* Tentative point for the min */
double ft = (*f)(t);
if ( ft <= fx )
{
/* t is a better approximation */
if ( t < x ) /* Reduce the range so that */
{
b = x; /* t would fall within it */
}
else
{
a = x;
}
v = w;
w = x;
x = t; /* Assign the best approx to x */
fv=fw;
fw=fx;
fx=ft;
}
else /* x remains the better approx */
{
if ( t < x ) /* Reduce the range enclosing x */
{
a = t;
}
else
{
b = t;
}
if ( ft <= fw || w==x )
{
v = w;
w = t;
fv=fw;
fw=ft;
}
else if ( ft<=fv || v==x || v==w )
{
v = t;
fv=ft;
}
}
} /* ----- end-of-block ----- */
} /* ===== End of loop ===== */
} |
Partager