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
|
;;; "minimize.scm" finds minimum f(x) for x0 <= x <= x1.
;;; Author: Lars Arvestad
;;;
;;; This code is in the public domain.
;;@noindent
;;
;;The Golden Section Search
;;@footnote{David Kahaner, Cleve Moler, and Stephen Nash
;;@cite{Numerical Methods and Software}
;;Prentice-Hall, 1989, ISBN 0-13-627258-4}
;;algorithm finds minima of functions which
;;are expensive to compute or for which derivatives are not available.
;;Although optimum for the general case, convergence is slow,
;;requiring nearly 100 iterations for the example (x^3-2x-5).
;;
;;@noindent
;;
;;If the derivative is available, Newton-Raphson is probably a better
;;choice. If the function is inexpensive to compute, consider
;;approximating the derivative.
;;@body
;;
;;@var{x_0} are @var{x_1} real numbers. The (single argument)
;;procedure @var{f} is unimodal over the open interval (@var{x_0},
;;@var{x_1}). That is, there is exactly one point in the interval for
;;which the derivative of @var{f} is zero.
;;
;;@0 returns a pair (@var{x} . @var{f}(@var{x})) where @var{f}(@var{x})
;;is the minimum. The @var{prec} parameter is the stop criterion. If
;;@var{prec} is a positive number, then the iteration continues until
;;@var{x} is within @var{prec} from the true value. If @var{prec} is
;;a negative integer, then the procedure will iterate @var{-prec}
;;times or until convergence. If @var{prec} is a procedure of seven
;;arguments, @var{x0}, @var{x1}, @var{a}, @var{b}, @var{fa}, @var{fb},
;;and @var{count}, then the iterations will stop when the procedure
;;returns @code{#t}.
;;
;;Analytically, the minimum of x^3-2x-5 is 0.816497.
;;@example
;;(define func (lambda (x) (+ (* x (+ (* x x) -2)) -5)))
;;(golden-section-search func 0 1 (/ 10000))
;; ==> (816.4883855245578e-3 . -6.0886621077391165)
;;(golden-section-search func 0 1 -5)
;; ==> (819.6601125010515e-3 . -6.088637561916407)
;;(golden-section-search func 0 1
;; (lambda (a b c d e f g ) (= g 500)))
;; ==> (816.4965933140557e-3 . -6.088662107903635)
;;@end example
(define golden-section-search
(let ((gss 'golden-section-search:)
(r (/ (- (sqrt 5) 1) 2))) ; 1 / golden-section
(lambda (f x0 x1 prec)
(cond ((not (procedure? f)) (slib:error gss 'procedure? f))
((not (number? x0)) (slib:error gss 'number? x0))
((not (number? x1)) (slib:error gss 'number? x1))
((>= x0 x1) (slib:error gss x0 'not '< x1)))
(let ((stop?
(cond
((procedure? prec) prec)
((number? prec)
(if (>= prec 0)
(lambda (x0 x1 a b fa fb count) (<= (abs (- x1 x0)) prec))
(if (integer? prec)
(lambda (x0 x1 a b fa fb count) (>= count (- prec)))
(slib:error gss 'integer? prec))))
(else (slib:error gss 'procedure? prec))))
(a0 (+ x0 (* (- x1 x0) (- 1 r))))
(b0 (+ x0 (* (- x1 x0) r)))
(delta #f)
(fmax #f)
(fmin #f))
(let loop ((left x0)
(right x1)
(a a0)
(b b0)
(fa (f a0))
(fb (f b0))
(count 1))
(define finish
(lambda (x fx)
(if (> fx fmin) (slib:warn gss fx 'not 'min (list '> fmin)))
(if (and (> count 9) (or (eqv? x0 left) (eqv? x1 right)))
(slib:warn gss 'min 'not 'found))
(cons x fx)))
(case count
((1)
(set! fmax (max fa fb))
(set! fmin (min fa fb)))
((2)
(set! fmin (min fmin fa fb))
(if (= fmax fa fb) (slib:error gss 'flat? fmax)))
(else
(set! fmin (min fmin fa fb))))
(cond ((stop? left right a b fa fb count)
(if (< fa fb)
(finish a fa)
(finish b fb)))
((< fa fb)
(let ((a-next (+ left (* (- b left) (- 1 r)))))
(cond ((and delta (< delta (- b a)))
(finish a fa))
(else (set! delta (- b a))
(loop left b a-next a (f a-next) fa
(+ 1 count))))))
(else
(let ((b-next (+ a (* (- right a) r))))
(cond ((and delta (< delta (- b a)))
(finish b fb))
(else (set! delta (- b a))
(loop a right b b-next fb (f b-next)
(+ 1 count))))))))))))
|
