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 ===== */
} |
Fonction minimum qui ne marche pas je ne sais pourquoi
)
Répondre avec citation
.
+ 
Partager