roots-one.lisp 14.2 KB
Newer Older
1 2
;; One-dimensional root solver.
;; Liam Healy 
3
;; Time-stamp: <2009-03-18 20:20:03EDT roots-one.lisp>
lhealy's avatar
lhealy committed
4
;; $Id$
5 6 7

(in-package :gsl)

8 9
;;; /usr/include/gsl/gsl_roots.h

10 11 12 13
;;;;****************************************************************************
;;;; Initialization
;;;;****************************************************************************

14
(defmobject one-dimensional-root-solver-f "gsl_root_fsolver"
15
  ((type :pointer))
Liam Healy's avatar
Liam Healy committed
16 17
  "one-dimensional root solver with function only"
  :initialize-suffix "set"
18
  :initialize-args ((callback :pointer) (lower :double) (upper :double))
19
  :callbacks (callback gsl-function (function))
20
  :singular (function))
21

22
(defmobject one-dimensional-root-solver-fdf "gsl_root_fdfsolver"
23
  ((type :pointer))
Liam Healy's avatar
Liam Healy committed
24 25
  "one-dimensional root solver with function and derivative"
  :initialize-suffix "set"
26 27 28 29 30 31 32 33 34 35 36
  :initialize-args ((callback :pointer) (root-guess :double))
  :callbacks
  (callback gsl-function-fdf
	    (function)
	    (df)
	    (fdf :void :double (:double :cvector 1) (:double :cvector 1)))
  :arglists-function
  (lambda (set)
    `((type &optional (function nil ,set) df fdf root-guess)
      (:type type)
      (:functions (list function df fdf) :root-guess root-guess))))
37

38
(defmfun name ((solver one-dimensional-root-solver-f))
39
  "gsl_root_fsolver_name"
40 41
  (((mpointer solver) :pointer))
  :definition :method
42
  :c-return :string
liam's avatar
liam committed
43
  :documentation			; FDL
44 45
  "The name of the solver.")

46
(defmfun name ((solver one-dimensional-root-solver-fdf))
47
  "gsl_root_fdfsolver_name"
48 49
  (((mpointer solver) :pointer))
  :definition :method
50
  :c-return :string
liam's avatar
liam committed
51
  :documentation			; FDL
52 53 54 55 56 57 58
  "The name of the solver.")

;;;;****************************************************************************
;;;; Iteration
;;;;****************************************************************************

;; It appears that this is always returning :SUCCESS (0).
59
(defmfun iterate ((solver one-dimensional-root-solver-f))
60
  "gsl_root_fsolver_iterate"
61
  (((mpointer solver) :pointer))
62
  :definition :method
63
  :callback-object solver
liam's avatar
liam committed
64
  :documentation			; FDL
65 66 67 68 69 70
  "Perform a single iteration of the solver.  The following errors may
   be signalled: 'bad-function-supplied, the iteration encountered a
   singular point where the function or its derivative evaluated to
   infinity or NaN, or 'gsl-division-by-zero, the derivative of the
   function vanished at the iteration point, preventing the algorithm
   from continuing without a division by zero.")
71

72
(defmfun iterate ((solver one-dimensional-root-solver-fdf))
73
  "gsl_root_fdfsolver_iterate"
74
  (((mpointer solver) :pointer))
75
  :definition :method
76
  :callback-object solver
liam's avatar
liam committed
77
  :documentation			; FDL
78 79 80 81 82 83
  "Perform a single iteration of the solver.  The following errors may
   be signalled: 'bad-function-supplied, the iteration encountered a
   singular point where the function or its derivative evaluated to
   infinity or NaN, or 'gsl-division-by-zero, the derivative of the
   function vanished at the iteration point, preventing the algorithm
   from continuing without a division by zero.")
84

85
(defmfun solution ((solver one-dimensional-root-solver-f))
86
  "gsl_root_fsolver_root"
87
  (((mpointer solver) :pointer))
88
  :definition :method
89
  :c-return :double
liam's avatar
liam committed
90
  :documentation			; FDL
91 92
  "The current estimate of the root for the solver.")

93
(defmfun solution ((solver one-dimensional-root-solver-fdf))
94
  "gsl_root_fdfsolver_root"
95
  (((mpointer solver) :pointer))
96
  :definition :method
97
  :c-return :double
liam's avatar
liam committed
98
  :documentation			; FDL
99 100
  "The current estimate of the root for the solver.")

liam's avatar
liam committed
101
(defmfun fsolver-lower (solver)
102
  "gsl_root_fsolver_x_lower"
103
  (((mpointer solver) :pointer))
104
  :c-return :double
liam's avatar
liam committed
105
  :documentation			; FDL
106 107
  "The lower end of the current bracketing interval for the solver.")

liam's avatar
liam committed
108
(defmfun fsolver-upper (solver)
109
  "gsl_root_fsolver_x_upper"
110
  (((mpointer solver) :pointer))
111
  :c-return :double
liam's avatar
liam committed
112
  :documentation 			; FDL
113 114 115 116 117 118
  "The upper end of the current bracketing interval for the solver.")

;;;;****************************************************************************
;;;; Search stopping conditions
;;;;****************************************************************************

liam's avatar
liam committed
119
(defmfun root-test-interval (lower upper absolute-error relative-error)
120 121 122
  "gsl_root_test_interval"
  ((lower :double) (upper :double)
   (absolute-error :double) (relative-error :double))
123
  :c-return :success-continue	 ; GSL documentation not clear on this
liam's avatar
liam committed
124
  :documentation			; FDL
125 126 127 128 129
  "Test for the convergence of the interval [lower,upper]
   with absolute error absolute-error and relative error
   relative-error.  This returns T
   if the following condition is achieved,
   |a - b| < epsabs + epsrel min(|a|,|b|) 
liam's avatar
liam committed
130 131 132
   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
133 134 135 136
   ensures that the relative error is accurately estimated for roots close
   to the origin.

   This condition on the interval also implies that any estimate of the
liam's avatar
liam committed
137 138 139
   root r in the interval satisfies the same condition with respect
   to the true root r^*, |r - r^*| < epsabs + epsrel r^*
   assuming that the true root r^* is contained within the interval.")
140

liam's avatar
liam committed
141
(defmfun root-test-delta (x1 x0 absolute-error relative-error)
142 143 144
  "gsl_root_test_delta"
  ((x1 :double) (x0 :double)
   (absolute-error :double) (relative-error :double))
145
  :c-return :success-continue
liam's avatar
liam committed
146
  :documentation			; FDL
147 148 149 150 151 152 153
  "Test for the convergence of the sequence ... x0, x1
   with absolute error absolute-error and relative error
   relative-error.  The test returns T if the following
   condition is achieved,
   |x_1 - x_0| < epsabs + epsrel |x_1|
   and returns NIL otherwise.")

liam's avatar
liam committed
154
(defmfun root-test-residual (f absolute-error)
155 156
  "gsl_root_test_residual"
  ((f :double) (absolute-error :double))
157
  :c-return :success-continue
liam's avatar
liam committed
158
  :documentation			; FDL
159 160 161 162 163
  "Tests the residual value f against the absolute
   error bound absolute-error.  The test returns T if the
   following condition is achieved,
   |f| < epsabs
   and returns NIL otherwise.  This criterion is suitable
liam's avatar
liam committed
164
   for situations where the precise location of the root, x, is
165 166 167 168 169 170 171
   unimportant provided a value can be found where the residual,
   |f(x)|, is small enough.")

;;;;****************************************************************************
;;;; Root bracketing algorithms
;;;;****************************************************************************

Liam Healy's avatar
Liam Healy committed
172
(defmpar +bisection-fsolver+ "gsl_root_fsolver_bisection"
liam's avatar
liam committed
173
  ;; FDL
174 175 176 177 178 179 180 181 182 183 184 185 186 187
  "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.

   On each iteration, the interval is bisected and the value of the
   function at the midpoint is calculated.  The sign of this value is used
   to determine which half of the interval does not contain a root.  That
   half is discarded to give a new, smaller interval containing the
   root.  This procedure can be continued indefinitely until the interval is
   sufficiently small.

   At any time the current estimate of the root is taken as the midpoint of
   the interval.")

Liam Healy's avatar
Liam Healy committed
188
(defmpar +false-position-fsolver+ "gsl_root_fsolver_falsepos"
liam's avatar
liam committed
189
  ;; FDL
190 191 192 193 194 195 196 197 198 199 200 201 202 203 204
  "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.

   On each iteration a line is drawn between the endpoints (a,f(a))
   and (b,f(b)) and the point where this line crosses the
   x-axis taken as a ``midpoint''.  The value of the function at
   this point is calculated and its sign is used to determine which side of
   the interval does not contain a root.  That side is discarded to give a
   new, smaller interval containing the root.  This procedure can be
   continued indefinitely until the interval is sufficiently small.

   The best estimate of the root is taken from the linear interpolation of
   the interval on the current iteration.")

Liam Healy's avatar
Liam Healy committed
205
(defmpar +brent-fsolver+ "gsl_root_fsolver_brent"
liam's avatar
liam committed
206
  ;; FDL
207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228
  "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.

   On each iteration Brent's method approximates the function using an
   interpolating curve.  On the first iteration this is a linear
   interpolation of the two endpoints.  For subsequent iterations the
   algorithm uses an inverse quadratic fit to the last three points, for
   higher accuracy.  The intercept of the interpolating curve with the
   x-axis is taken as a guess for the root.  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 bisection
   step.

   The best estimate of the root is taken from the most recent
   interpolation or bisection.")

;;;;****************************************************************************
;;;; Root finding algorithms using derivatives
;;;;****************************************************************************

Liam Healy's avatar
Liam Healy committed
229
(defmpar +newton-fdfsolver+ "gsl_root_fdfsolver_newton"
liam's avatar
liam committed
230
  ;; FDL
231 232
  "Newton's Method is the standard root-polishing algorithm.  The algorithm
   begins with an initial guess for the location of the root.  On each
liam's avatar
liam committed
233 234
   iteration, a line tangent to the function f is drawn at that
   position.  The point where this line crosses the x-axis becomes
235 236 237 238 239
   the new guess.  The iteration is defined by the following sequence,
   x_{i+1} = x_i - f(x_i) / f'(x_i)
   Newton's method converges quadratically for single roots, and linearly
   for multiple roots.")

Liam Healy's avatar
Liam Healy committed
240
(defmpar +secant-fdfsolver+ "gsl_root_fdfsolver_secant"
liam's avatar
liam committed
241
  ;; FDL
242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264
  "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
   the derivative to compute a first step,
   x_1 = x_0 - f(x_0)/f'(x_0)
   Subsequent iterations avoid the evaluation of the derivative by
   replacing it with a numerical estimate, the slope of the line through
   the previous two points,
   x_{i+1} = x_i - f(x_i) / f'_{est}
    where
   f'_{est} =  f(x_{i}) - f(x_{i-1}) / x_i - x_{i-1}
   When the derivative does not change significantly in the vicinity of the
   root the secant method gives a useful saving.  Asymptotically the secant
   method is faster than Newton's method whenever the cost of evaluating
   the derivative is more than 0.44 times the cost of evaluating the
   function itself.  As with all methods of computing a numerical
   derivative the estimate can suffer from cancellation errors if the
   separation of the points becomes too small.

   On single roots, the method has a convergence of order (1 + \sqrt
   5)/2 (approximately 1.62).  It converges linearly for multiple
   roots.")

Liam Healy's avatar
Liam Healy committed
265
(defmpar +steffenson-fdfsolver+ "gsl_root_fdfsolver_steffenson"
liam's avatar
liam committed
266
  ;; FDL
267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289
  "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
   then the acceleration procedure generates a new sequence R_i,
   R_i = x_i - (x_{i+1} - x_i)^2 / (x_{i+2} - 2 x_{i+1} + x_i)
   which converges faster than the original sequence under reasonable
   conditions.  The new sequence requires three terms before it can produce
   its first value so the method returns accelerated values on the second
   and subsequent iterations.  On the first iteration it returns the
   ordinary Newton estimate.  The Newton iterate is also returned if the
   denominator of the acceleration term ever becomes zero.

   As with all acceleration procedures this method can become unstable if
   the function is not well-behaved.")

;;;;****************************************************************************
;;;; Examples
;;;;****************************************************************************

;;; This is the example given in Sec. 32.10.

(let ((a 1.0d0) (b 0.0d0) (c -5.0d0))
  (defun quadratic (x)
290
    (+ (* (+ (* a x) b) x) c))
291
  (defun quadratic-derivative (x)
292
    (+ (* 2 a x) b))
293 294 295
  (defun quadratic-and-derivative (x)
    (values (+ (* (+ (* a x) b) x) c)
	    (+ (* 2 a x) b))))
296

297
(defun roots-one-example-no-derivative
Liam Healy's avatar
Liam Healy committed
298
    (&optional (method +brent-fsolver+) (print-steps t))
299
  "Solving a quadratic, the example given in Sec. 32.10 of the GSL manual."
300
  (let ((max-iter 50)
301
	(solver
302
	 (make-one-dimensional-root-solver-f method 'quadratic 0.0d0 5.0d0)))
303 304
    (when print-steps
      (format t "iter ~6t   [lower ~24tupper] ~36troot ~44terr ~54terr(est)~&"))
305
    (loop for iter from 0
306
       for root = (solution solver)
307 308
       for lower = (fsolver-lower solver)
       for upper = (fsolver-upper solver)
309
       do (iterate solver)
310 311 312
       while  (and (< iter max-iter)
		   (not (root-test-interval lower upper 0.0d0 0.001d0)))
       do
313 314 315 316 317 318
       (when print-steps
	 (format t "~d~6t~10,6f~18t~10,6f~28t~12,9f ~44t~10,4g ~10,4g~&"
		 iter lower upper
		 root (- root (sqrt 5.0d0))
		 (- upper lower)))
       finally (return root))))
319

320
(defun roots-one-example-derivative
Liam Healy's avatar
Liam Healy committed
321
    (&optional (method +newton-fdfsolver+) (print-steps t))
322
  "Solving a quadratic, the example given in Sec. 32.10 of the GSL manual."
323 324 325
  (let* ((max-iter 100)
	 (initial 5.0d0)
	 (solver (make-one-dimensional-root-solver-fdf
326
		  method
327
		  'quadratic 'quadratic-derivative 'quadratic-and-derivative
328
		  initial)))
329 330
    (when print-steps
      (format t "iter ~6t ~8troot ~22terr ~34terr(est)~&"))
331 332
    (loop for iter from 0
       for oldroot = initial then root
333
       for root = (progn (iterate solver) (solution solver))
334 335 336
       while (and (< iter max-iter)
		  (not (root-test-delta root oldroot 0.0d0 1.0d-5)))
       do
337 338 339 340
       (when print-steps
	 (format t "~d~6t~10,8g ~18t~10,6g~34t~10,6g~&"
		 iter root (- root (sqrt 5.0d0)) (- root oldroot)))
       finally (return root))))
341

342 343
;; To see step-by-step information as the solution progresses, make
;; the last argument T.
344
(save-test roots-one
Liam Healy's avatar
Liam Healy committed
345 346 347 348 349 350
 (roots-one-example-no-derivative +bisection-fsolver+ nil)
 (roots-one-example-no-derivative +false-position-fsolver+ nil)
 (roots-one-example-no-derivative +brent-fsolver+ nil)
 (roots-one-example-derivative +newton-fdfsolver+ nil)
 (roots-one-example-derivative +secant-fdfsolver+ nil)
 (roots-one-example-derivative +steffenson-fdfsolver+ nil))