roots-one.lisp 14.2 KB
 liam committed Jan 16, 2008 1 2 ;; One-dimensional root solver. ;; Liam Healy  Liam Healy committed Mar 18, 2009 3 ;; Time-stamp: <2009-03-18 20:20:03EDT roots-one.lisp>  lhealy committed Jul 24, 2008 4 ;; $Id$  liam committed Jan 08, 2008 5 6 7  (in-package :gsl)  8 9 ;;; /usr/include/gsl/gsl_roots.h  liam committed Jan 08, 2008 10 11 12 13 ;;;;**************************************************************************** ;;;; Initialization ;;;;****************************************************************************  Liam Healy committed Dec 25, 2008 14 (defmobject one-dimensional-root-solver-f "gsl_root_fsolver"  Liam Healy committed Dec 23, 2008 15  ((type :pointer))  Liam Healy committed Jan 25, 2009 16 17  "one-dimensional root solver with function only" :initialize-suffix "set"  Liam Healy committed Feb 08, 2009 18  :initialize-args ((callback :pointer) (lower :double) (upper :double))  Liam Healy committed Mar 15, 2009 19  :callbacks (callback gsl-function (function))  Liam Healy committed Feb 10, 2009 20  :singular (function))  Liam Healy committed Feb 07, 2009 21   Liam Healy committed Dec 25, 2008 22 (defmobject one-dimensional-root-solver-fdf "gsl_root_fdfsolver"  Liam Healy committed Dec 23, 2008 23  ((type :pointer))  Liam Healy committed Jan 25, 2009 24 25  "one-dimensional root solver with function and derivative" :initialize-suffix "set"  Liam Healy committed Mar 15, 2009 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))))  Liam Healy committed Feb 08, 2009 37   Liam Healy committed Dec 26, 2008 38 (defmfun name ((solver one-dimensional-root-solver-f))  liam committed Jan 08, 2008 39  "gsl_root_fsolver_name"  Liam Healy committed Dec 26, 2008 40 41  (((mpointer solver) :pointer)) :definition :method  liam committed Jan 08, 2008 42  :c-return :string  liam committed Feb 18, 2008 43  :documentation ; FDL  liam committed Jan 08, 2008 44 45  "The name of the solver.")  Liam Healy committed Dec 26, 2008 46 (defmfun name ((solver one-dimensional-root-solver-fdf))  liam committed Jan 08, 2008 47  "gsl_root_fdfsolver_name"  Liam Healy committed Dec 26, 2008 48 49  (((mpointer solver) :pointer)) :definition :method  liam committed Jan 08, 2008 50  :c-return :string  liam committed Feb 18, 2008 51  :documentation ; FDL  liam committed Jan 08, 2008 52 53 54 55 56 57 58  "The name of the solver.") ;;;;**************************************************************************** ;;;; Iteration ;;;;**************************************************************************** ;; It appears that this is always returning :SUCCESS (0).  Liam Healy committed Jan 03, 2009 59 (defmfun iterate ((solver one-dimensional-root-solver-f))  liam committed Jan 08, 2008 60  "gsl_root_fsolver_iterate"  Liam Healy committed Dec 26, 2008 61  (((mpointer solver) :pointer))  Liam Healy committed Jan 03, 2009 62  :definition :method  Liam Healy committed Mar 15, 2009 63  :callback-object solver  liam committed Feb 18, 2008 64  :documentation ; FDL  Liam Healy committed Nov 16, 2008 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.")  liam committed Jan 08, 2008 71   Liam Healy committed Jan 03, 2009 72 (defmfun iterate ((solver one-dimensional-root-solver-fdf))  liam committed Jan 08, 2008 73  "gsl_root_fdfsolver_iterate"  Liam Healy committed Dec 26, 2008 74  (((mpointer solver) :pointer))  Liam Healy committed Jan 03, 2009 75  :definition :method  Liam Healy committed Mar 15, 2009 76  :callback-object solver  liam committed Feb 18, 2008 77  :documentation ; FDL  Liam Healy committed Nov 16, 2008 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.")  liam committed Jan 08, 2008 84   Liam Healy committed Jan 03, 2009 85 (defmfun solution ((solver one-dimensional-root-solver-f))  liam committed Jan 08, 2008 86  "gsl_root_fsolver_root"  Liam Healy committed Dec 26, 2008 87  (((mpointer solver) :pointer))  Liam Healy committed Jan 03, 2009 88  :definition :method  liam committed Jan 08, 2008 89  :c-return :double  liam committed Feb 18, 2008 90  :documentation ; FDL  liam committed Jan 08, 2008 91 92  "The current estimate of the root for the solver.")  Liam Healy committed Jan 03, 2009 93 (defmfun solution ((solver one-dimensional-root-solver-fdf))  liam committed Jan 08, 2008 94  "gsl_root_fdfsolver_root"  Liam Healy committed Dec 26, 2008 95  (((mpointer solver) :pointer))  Liam Healy committed Jan 03, 2009 96  :definition :method  liam committed Jan 08, 2008 97  :c-return :double  liam committed Feb 18, 2008 98  :documentation ; FDL  liam committed Jan 08, 2008 99 100  "The current estimate of the root for the solver.")  liam committed Feb 18, 2008 101 (defmfun fsolver-lower (solver)  liam committed Jan 08, 2008 102  "gsl_root_fsolver_x_lower"  Liam Healy committed Dec 26, 2008 103  (((mpointer solver) :pointer))  liam committed Jan 08, 2008 104  :c-return :double  liam committed Feb 18, 2008 105  :documentation ; FDL  liam committed Jan 08, 2008 106 107  "The lower end of the current bracketing interval for the solver.")  liam committed Feb 18, 2008 108 (defmfun fsolver-upper (solver)  liam committed Jan 08, 2008 109  "gsl_root_fsolver_x_upper"  Liam Healy committed Dec 26, 2008 110  (((mpointer solver) :pointer))  liam committed Jan 08, 2008 111  :c-return :double  liam committed Feb 18, 2008 112  :documentation ; FDL  liam committed Jan 08, 2008 113 114 115 116 117 118  "The upper end of the current bracketing interval for the solver.") ;;;;**************************************************************************** ;;;; Search stopping conditions ;;;;****************************************************************************  liam committed Feb 18, 2008 119 (defmfun root-test-interval (lower upper absolute-error relative-error)  liam committed Jan 08, 2008 120 121 122  "gsl_root_test_interval" ((lower :double) (upper :double) (absolute-error :double) (relative-error :double))  liam committed Jan 14, 2008 123  :c-return :success-continue ; GSL documentation not clear on this  liam committed Feb 18, 2008 124  :documentation ; FDL  liam committed Jan 08, 2008 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 committed Feb 18, 2008 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  liam committed Jan 08, 2008 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 committed Feb 18, 2008 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.")  liam committed Jan 08, 2008 140   liam committed Feb 18, 2008 141 (defmfun root-test-delta (x1 x0 absolute-error relative-error)  liam committed Jan 08, 2008 142 143 144  "gsl_root_test_delta" ((x1 :double) (x0 :double) (absolute-error :double) (relative-error :double))  liam committed Jan 14, 2008 145  :c-return :success-continue  liam committed Feb 18, 2008 146  :documentation ; FDL  liam committed Jan 08, 2008 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 committed Feb 18, 2008 154 (defmfun root-test-residual (f absolute-error)  liam committed Jan 08, 2008 155 156  "gsl_root_test_residual" ((f :double) (absolute-error :double))  liam committed Jan 14, 2008 157  :c-return :success-continue  liam committed Feb 18, 2008 158  :documentation ; FDL  liam committed Jan 08, 2008 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 committed Feb 18, 2008 164  for situations where the precise location of the root, x, is  liam committed Jan 08, 2008 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 committed Feb 16, 2009 172 (defmpar +bisection-fsolver+ "gsl_root_fsolver_bisection"  liam committed Feb 18, 2008 173  ;; FDL  liam committed Jan 08, 2008 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 committed Feb 16, 2009 188 (defmpar +false-position-fsolver+ "gsl_root_fsolver_falsepos"  liam committed Feb 18, 2008 189  ;; FDL  liam committed Jan 08, 2008 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 committed Feb 16, 2009 205 (defmpar +brent-fsolver+ "gsl_root_fsolver_brent"  liam committed Feb 18, 2008 206  ;; FDL  liam committed Jan 08, 2008 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 committed Feb 16, 2009 229 (defmpar +newton-fdfsolver+ "gsl_root_fdfsolver_newton"  liam committed Feb 18, 2008 230  ;; FDL  liam committed Jan 08, 2008 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 committed Feb 18, 2008 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  liam committed Jan 08, 2008 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 committed Feb 16, 2009 240 (defmpar +secant-fdfsolver+ "gsl_root_fdfsolver_secant"  liam committed Feb 18, 2008 241  ;; FDL  liam committed Jan 08, 2008 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 committed Feb 16, 2009 265 (defmpar +steffenson-fdfsolver+ "gsl_root_fdfsolver_steffenson"  liam committed Feb 18, 2008 266  ;; FDL  liam committed Jan 08, 2008 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)  liam committed Jan 13, 2008 290  (+ (* (+ (* a x) b) x) c))  liam committed Jan 08, 2008 291  (defun quadratic-derivative (x)  liam committed Jan 13, 2008 292  (+ (* 2 a x) b))  Liam Healy committed Jan 19, 2009 293 294 295  (defun quadratic-and-derivative (x) (values (+ (* (+ (* a x) b) x) c) (+ (* 2 a x) b))))  liam committed Jan 08, 2008 296   Liam Healy committed Jan 22, 2009 297 (defun roots-one-example-no-derivative  Liam Healy committed Feb 16, 2009 298  (&optional (method +brent-fsolver+) (print-steps t))  liam committed Jan 08, 2008 299  "Solving a quadratic, the example given in Sec. 32.10 of the GSL manual."  Liam Healy committed Dec 26, 2008 300  (let ((max-iter 50)  Liam Healy committed Dec 26, 2008 301  (solver  Liam Healy committed Feb 08, 2009 302  (make-one-dimensional-root-solver-f method 'quadratic 0.0d0 5.0d0)))  Liam Healy committed Jan 19, 2009 303 304  (when print-steps (format t "iter ~6t [lower ~24tupper] ~36troot ~44terr ~54terr(est)~&"))  liam committed Feb 03, 2008 305  (loop for iter from 0  Liam Healy committed Jan 03, 2009 306  for root = (solution solver)  Liam Healy committed Dec 26, 2008 307 308  for lower = (fsolver-lower solver) for upper = (fsolver-upper solver)  Liam Healy committed Jan 03, 2009 309  do (iterate solver)  Liam Healy committed Dec 26, 2008 310 311 312  while (and (< iter max-iter) (not (root-test-interval lower upper 0.0d0 0.001d0))) do  Liam Healy committed Jan 19, 2009 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))))  liam committed Jan 08, 2008 319   Liam Healy committed Jan 22, 2009 320 (defun roots-one-example-derivative  Liam Healy committed Feb 16, 2009 321  (&optional (method +newton-fdfsolver+) (print-steps t))  liam committed Jan 08, 2008 322  "Solving a quadratic, the example given in Sec. 32.10 of the GSL manual."  Liam Healy committed Dec 26, 2008 323 324 325  (let* ((max-iter 100) (initial 5.0d0) (solver (make-one-dimensional-root-solver-fdf  Liam Healy committed Jan 24, 2009 326  method  Liam Healy committed Mar 15, 2009 327  'quadratic 'quadratic-derivative 'quadratic-and-derivative  Liam Healy committed Jan 24, 2009 328  initial)))  Liam Healy committed Jan 19, 2009 329 330  (when print-steps (format t "iter ~6t ~8troot ~22terr ~34terr(est)~&"))  Liam Healy committed Dec 26, 2008 331 332  (loop for iter from 0 for oldroot = initial then root  Liam Healy committed Jan 03, 2009 333  for root = (progn (iterate solver) (solution solver))  Liam Healy committed Dec 26, 2008 334 335 336  while (and (< iter max-iter) (not (root-test-delta root oldroot 0.0d0 1.0d-5))) do  Liam Healy committed Jan 19, 2009 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))))  liam committed Jan 16, 2008 341   Liam Healy committed Jan 22, 2009 342 343 ;; To see step-by-step information as the solution progresses, make ;; the last argument T.  Liam Healy committed Jan 19, 2009 344 (save-test roots-one  Liam Healy committed Feb 16, 2009 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))`