Commit 058e342b authored by liam's avatar liam

Add minimization-one for univariate minimization.


git-svn-id: svn+ssh://pop/opt/space/mathematics/gsl/trunk@3260 a3d8a0fb-c1db-0310-ace7-a616afeb9e30
parent dede295f
......@@ -3,7 +3,7 @@
; description: Definition of GSLL system
; date:
; author: Liam Healy
; modified: Mon Jan 7 2008 - 00:13
; modified: Tue Jan 8 2008 - 22:52
;********************************************************
;;; $Id: $
......@@ -168,12 +168,10 @@
(:file "types")
(:file "lookup")
(:file "evaluation")))
;; Once common callback definitions are moved out of
;; numerical-integration to something in init, the
;; numerical-integration dependency won't be necessary.
(:file "numerical-differentiation" :depends-on (init general))
(:file "chebyshev" :depends-on (init general))
(:file "series-acceleration" :depends-on (init))
(:file "wavelet" :depends-on (init data))
(:file "hankel" :depends-on (init data))
(:file "roots-one" :depends-on (init))))
(:file "roots-one" :depends-on (init general))
(:file "minimization-one" :depends-on (init general))))
......@@ -307,7 +307,7 @@ svn checkout svn://common-lisp.net/project/gsll/subversion/trunk</pre>
<tr>
<td><a
href="One-dimensional-Minimization.html">One-dimensional Minimization</a></td>
<td>Not started</td>
<td>Done</td>
</tr>
<tr>
<td><a href="Multidimensional-Root_002dFinding.html">
......
;********************************************************
; file: minimization-one.lisp
; description: Univariate minimization
; date: Tue Jan 8 2008 - 21:02
; author: Liam Healy
; modified: Tue Jan 8 2008 - 22:50
;********************************************************
;;; $Id: $
(in-package :gsl)
;;;;****************************************************************************
;;;; Initialization
;;;;****************************************************************************
(defun-gsl allocate-fminimizer (type)
"gsl_min_fminimizer_alloc"
((type :pointer))
:c-return :pointer
:documentation
"Allocate an instance of a minimizer of the given type.")
(defun-gsl set-fminimizer (minimizer function minimum lower upper)
"gsl_min_fminimizer_set"
((minimizer :pointer) (function :pointer)
(minimum :double) (lower :double) (upper :double))
:documentation
"Set, or reset, an existing minimizer to use the
function and the initial search interval [lower,
upper], with a guess for the location of the minimum.")
(defun-gsl set-fminimizer-with-values
(minimizer function x-minimum x-lower x-upper
f-minimum f-lower f-upper)
"gsl_min_fminimizer_set_with_values"
((minimizer :pointer) (function :pointer)
(x-minimum :double) (x-lower :double) (x-upper :double)
(f-minimum :double) (f-lower :double) (f-upper :double))
:documentation
"Set, or reset, an existing minimizer to use the
function and the initial search interval [lower,
upper], with a guess for the location of the minimum, using
supplied rather than computed values of the function.")
(defun-gsl free-fminimizer (minimizer)
"gsl_min_fminimizer_free"
((minimizer :pointer))
:c-return :void
:documentation
"Free all the memory associated with the minimizer.")
(defun-gsl fminimizer-name (minimizer)
"gsl_min_fminimizer_name"
((minimizer :pointer))
:c-return :string
:documentation
"The name of the minimizer.")
(export '(with-fminimizer))
(defmacro with-fminimizer
((minimizer minimizer-type function minimum lower upper) &body body)
"Create and initialize an fminimizer for one-dimensional problems,
and clean up afterwards."
`(let ((,minimizer (allocate-fminimizer ,minimizer-type)))
(unwind-protect
(progn
(set-fminimizer ,minimizer ,function ,minimum ,lower ,upper)
,@body)
(free-fminimizer ,minimizer))))
;;;;****************************************************************************
;;;; Iteration
;;;;****************************************************************************
(defun-gsl iterate-fminimizer (minimizer)
"gsl_min_fminimizer_iterate"
((minimizer :pointer))
:c-return :success-continue
:documentation
"Perform a single iteration of the minimizer. The following
errors may be signalled: :EBADFUNC,
the iteration encountered a singular point where the function or its
derivative evaluated to infinity or NaN, or
:FAILURE, the algorithm could not improve the current best approximation or
bounding interval.")
(defun-gsl fminimizer-x-minimum (minimizer)
"gsl_min_fminimizer_x_lower"
((minimizer :pointer))
:c-return :double
:documentation
"The current estimate of the position of the minimum for the minimizer.")
(defun-gsl fminimizer-x-lower (minimizer)
"gsl_min_fminimizer_x_lower"
((minimizer :pointer))
:c-return :double
:documentation
"The current lower bound of the interval for the minimizer.")
(defun-gsl fminimizer-x-upper (minimizer)
"gsl_min_fminimizer_x_upper"
((minimizer :pointer))
:c-return :double
:documentation
"The current upper bound of the interval for the minimizer.")
(defun-gsl fminimizer-f-minimum (minimizer)
"gsl_min_fminimizer_f_minimum"
((minimizer :pointer))
:c-return :double
:documentation
"The value of the function at the current estimate of the minimum for the
minimizer.")
(defun-gsl fminimizer-f-lower (minimizer)
"gsl_min_fminimizer_f_lower"
((minimizer :pointer))
:c-return :double
:documentation
"The value of the function at the current estimate of the lower bound
for the minimizer.")
(defun-gsl fminimizer-f-upper (minimizer)
"gsl_min_fminimizer_f_upper"
((minimizer :pointer))
:c-return :double
:documentation
"The value of the function at the current estimate of the upper bound
for the minimizer.")
;;;;****************************************************************************
;;;; Stopping parameters
;;;;****************************************************************************
(defun-gsl min-test-interval (lower upper absolute-error relative-error)
"gsl_min_test_interval"
((lower :double) (upper :double)
(absolute-error :double) (relative-error :double))
:c-return :success-failure
:documentation
"Test for the convergence of the interval [lower,upper]
with absolute error and relative error specified.
The test returns T if the following condition is achieved:
|a - b| < epsabs + epsrel min(|a|,|b|)
when the interval x = [a,b] does not include the origin. If the
interval includes the origin then min(|a|,|b|) is replaced by
zero (which is the minimum value of |x| over the interval). This
ensures that the relative error is accurately estimated for minima close
to the origin.
This condition on the interval also implies that any estimate of the
minimum x_m in the interval satisfies the same condition with respect
to the true minimum x_m^*,
|x_m - x_m^*| < epsabs + epsrel x_m^*
assuming that the true minimum x_m^* is contained within the interval.")
;;;;****************************************************************************
;;;; Minimization algorithms
;;;;****************************************************************************
(defvariable *golden-section-fminimizer* "gsl_min_fminimizer_goldensection"
"The golden section algorithm is the simplest method of bracketing
the minimum of a function. It is the slowest algorithm provided by the
library, with linear convergence.
On each iteration, the algorithm first compares the subintervals from
the endpoints to the current minimum. The larger subinterval is divided
in a golden section (using the famous ratio (3-sqrt 5)/2 =
0.3189660...) and the value of the function at this new point is
calculated. The new value is used with the constraint f(a') >
f(x') < f(b') to a select new interval containing the minimum, by
discarding the least useful point. This procedure can be continued
indefinitely until the interval is sufficiently small. Choosing the
golden section as the bisection ratio can be shown to provide the
fastest convergence for this type of algorithm.")
(defvariable *brent-fminimizer* "gsl_min_fminimizer_brent"
"The Brent minimization algorithm combines a parabolic
interpolation with the golden section algorithm. This produces a fast
algorithm which is still robust.
The outline of the algorithm can be summarized as follows: on each
iteration Brent's method approximates the function using an
interpolating parabola through three existing points. The minimum of the
parabola is taken as a guess for the minimum. If it lies within the
bounds of the current interval then the interpolating point is accepted,
and used to generate a smaller interval. If the interpolating point is
not accepted then the algorithm falls back to an ordinary golden section
step. The full details of Brent's method include some additional checks
to improve convergence.")
;;;;****************************************************************************
;;;; Example
;;;;****************************************************************************
;;; This is the example given in Sec. 33.8. The results are different
;;; than given there.
(defun-scalar minimization-one-fn (x)
(1+ (cos x)))
(defun minimization-one-example ()
"Solving a minimum, the example given in Sec. 33.8 of the GSL manual."
(let ((max-iter 100))
(with-fminimizer
(minimizer *brent-fminimizer* minimization-one-fn 2.0d0 0.0d0 6.0d0)
(format t "~&iter ~6t [lower ~24tupper] ~36tmin ~44tmin err ~54tupper-lower")
(loop for iter from 0
for min = (fminimizer-x-minimum minimizer)
for lower = (fminimizer-x-lower minimizer)
for upper = (fminimizer-x-upper minimizer)
do (iterate-fminimizer minimizer)
while (and (< iter max-iter)
;; abs and rel error swapped in example?
(not (min-test-interval lower upper 0.001d0 0.0d0)))
do
(format t "~&~d~6t~10,6f~18t~10,6f~28t~12,9f ~44t~10,4g ~10,4g"
iter lower upper
min (- min pi)
(- upper lower))))))
......@@ -3,7 +3,7 @@
; description: One-dimensional root solver.
; date: Sun Dec 9 2007 - 15:47
; author: Liam Healy
; modified: Mon Jan 7 2008 - 22:25
; modified: Tue Jan 8 2008 - 22:53
;********************************************************
;;; $Id: $
......@@ -204,7 +204,7 @@
;;;; Search stopping conditions
;;;;****************************************************************************
(defun-gsl test-interval (lower upper absolute-error relative-error)
(defun-gsl root-test-interval (lower upper absolute-error relative-error)
"gsl_root_test_interval"
((lower :double) (upper :double)
(absolute-error :double) (relative-error :double))
......@@ -227,7 +227,7 @@
|r - r^*| < epsabs + epsrel r^*
assuming that the true root @math{r^*} is contained within the interval.")
(defun-gsl test-delta (x1 x0 absolute-error relative-error)
(defun-gsl root-test-delta (x1 x0 absolute-error relative-error)
"gsl_root_test_delta"
((x1 :double) (x0 :double)
(absolute-error :double) (relative-error :double))
......@@ -240,7 +240,7 @@
|x_1 - x_0| < epsabs + epsrel |x_1|
and returns NIL otherwise.")
(defun-gsl test-residual (f absolute-error)
(defun-gsl root-test-residual (f absolute-error)
"gsl_root_test_residual"
((f :double) (absolute-error :double))
:c-return :success-failure
......@@ -258,7 +258,7 @@
;;;; Root bracketing algorithms
;;;;****************************************************************************
(defvariable bisection-fsolver "gsl_root_fsolver_bisection"
(defvariable *bisection-fsolver* "gsl_root_fsolver_bisection"
"The bisection algorithm is the simplest method of bracketing the
roots of a function. It is the slowest algorithm provided by
the library, with linear convergence.
......@@ -273,7 +273,7 @@
At any time the current estimate of the root is taken as the midpoint of
the interval.")
(defvariable false-position-fsolver "gsl_root_fsolver_falsepos"
(defvariable *false-position-fsolver* "gsl_root_fsolver_falsepos"
"The false position algorithm is a method of finding roots based on
linear interpolation. Its convergence is linear, but it is usually
faster than bisection.
......@@ -289,7 +289,7 @@
The best estimate of the root is taken from the linear interpolation of
the interval on the current iteration.")
(defvariable brent-fsolver "gsl_root_fsolver_brent"
(defvariable *brent-fsolver* "gsl_root_fsolver_brent"
"The Brent-Dekker method (referred to here as Brent's method)
combines an interpolation strategy with the bisection algorithm. This
produces a fast algorithm which is still robust.
......@@ -312,7 +312,7 @@
;;;; Root finding algorithms using derivatives
;;;;****************************************************************************
(defvariable newton-fdfsolver "gsl_root_fdfsolver_newton"
(defvariable *newton-fdfsolver* "gsl_root_fdfsolver_newton"
"Newton's Method is the standard root-polishing algorithm. The algorithm
begins with an initial guess for the location of the root. On each
iteration, a line tangent to the function @math{f} is drawn at that
......@@ -322,7 +322,7 @@
Newton's method converges quadratically for single roots, and linearly
for multiple roots.")
(defvariable secant-fdfsolver "gsl_root_fdfsolver_secant"
(defvariable *secant-fdfsolver* "gsl_root_fdfsolver_secant"
"The secant method is a simplified version of Newton's method which does
not require the computation of the derivative on every step.
On its first iteration the algorithm begins with Newton's method, using
......@@ -346,7 +346,7 @@
5)/2 (approximately 1.62). It converges linearly for multiple
roots.")
(defvariable steffenson-fdfsolver "gsl_root_fdfsolver_steffenson"
(defvariable *steffenson-fdfsolver* "gsl_root_fdfsolver_steffenson"
"The Steffenson method provides the fastest convergence of all the
routines. It combines the basic Newton algorithm with an Aitken
``delta-squared'' acceleration. If the Newton iterates are x_i
......@@ -383,7 +383,7 @@
(defun roots-one-example ()
"Solving a quadratic, the example given in Sec. 32.10 of the GSL manual."
(let ((max-iter 50))
(with-fsolver (solver brent-fsolver quadratic 0.0d0 5.0d0)
(with-fsolver (solver *brent-fsolver* quadratic 0.0d0 5.0d0)
(format t "~&iter ~6t [lower ~24tupper] ~36troot ~44terr ~54terr(est)")
(loop for iter from 0
for root = (fsolver-root solver)
......@@ -391,7 +391,7 @@
for upper = (fsolver-upper solver)
do (iterate-fsolver solver)
while (and (< iter max-iter)
(not (test-interval lower upper 0.0d0 0.001d0)))
(not (root-test-interval lower upper 0.0d0 0.001d0)))
do
(format t "~&~d~6t~10,6f~18t~10,6f~28t~12,9f ~44t~10,4g ~10,4g"
iter lower upper
......@@ -416,14 +416,14 @@
(let ((max-iter 100)
(initial 5.0d0))
(with-fdfsolver
(solver newton-fdfsolver *quadratic* initial)
(solver *newton-fdfsolver* *quadratic* initial)
(format t "~&iter ~6t ~36troot ~44terr ~54terr(est)")
(loop for iter from 0
for oldroot = initial then root
for root = (fdfsolver-root solver)
do (iterate-fdfsolver solver)
while (and (< iter max-iter)
(not (test-delta root oldroot 0.0d0 1.0d-3)))
(not (root-test-delta root oldroot 0.0d0 1.0d-3)))
do
(format t "~&~d~6t~g ~18t~10,4g ~10,4g"
iter root (- root (sqrt 5.0d0)) (- root oldroot))))))
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment