1 //
2 // minsearch.c
3 //
4 // search to find minimum of a function using parabolic
5 // interpolation
6 //
7
8 #include <stdio.h>
9 #include <stdlib.h>
10 #include <math.h>
11 #include <assert.h>
12
function(double _x)13 double function(double _x) {
14 #if 0
15 // -sin(x) minimum at x = pi/2
16 double y = 1.0f-sin(_x);
17 #else
18 // exp(-(x-pi/2)^2)
19 double t = _x - M_PI/2;
20 double y = -exp(-t*t);
21 #endif
22 return y;
23 }
24
main()25 int main() {
26 // options
27 unsigned int num_iterations = 10;
28 double x0 = 0.3f; // lower value
29 double x1; // middle value (ignored)
30 double x2 = 3.1f; // upper value
31 double tol = 1e-6f; // error tolerance
32
33 //
34 double y0 = function(x0);
35 double y1;
36 double y2 = function(x2);
37
38 double x_hat;
39 double del = 0.0f;
40 unsigned int i;
41 for (i=0; i<num_iterations; i++) {
42 // choose center point between [x0,x2]
43 x1 = 0.5f*(x0 + x2);
44 y1 = function(0.5f*(x0+x2));
45
46 #if 0
47 printf("%4u : (%8.3f,%8.3f,%8.3f) (%8.3f,%8.3f,%8.3f)\n",
48 i,
49 x0, x1, x2,
50 y0, y1, y2);
51 #endif
52
53 // ensure values are reasonable
54 assert(x2 > x1);
55 assert(x1 > x0);
56
57 // compute minimum
58 double t0 = y0*(x1*x1 - x2*x2) +
59 y1*(x2*x2 - x0*x0) +
60 y2*(x0*x0 - x1*x1);
61
62 double t1 = y0*(x1 - x2) +
63 y1*(x2 - x0) +
64 y2*(x0 - x1);
65
66 x_hat = 0.5f * t0 / t1;
67 del = x_hat - x1;
68
69 // print
70 printf("%4u : %12.8f, e=%12.4e\n", i, x_hat, x_hat-M_PI/2);
71
72 if (x_hat > x1) {
73 // new minimum
74 x0 = x1;
75 y0 = y1;
76 } else {
77 // new maximum
78 x2 = x1;
79 y2 = y1;
80 }
81
82 if (fabs(del) < tol)
83 break;
84 }
85
86 printf("estimated minimum : %f (%f)\n", x_hat, M_PI/2);
87
88 return 0;
89 }
90