diff options
Diffstat (limited to 'root.scm')
| -rw-r--r-- | root.scm | 149 |
1 files changed, 149 insertions, 0 deletions
diff --git a/root.scm b/root.scm new file mode 100644 index 0000000..5ba78c1 --- /dev/null +++ b/root.scm @@ -0,0 +1,149 @@ +;;;"root.scm" Newton's and Laguerre's methods for finding roots. +;Copyright (C) 1996 Aubrey Jaffer +; +;Permission to copy this software, to redistribute it, and to use it +;for any purpose is granted, subject to the following restrictions and +;understandings. +; +;1. Any copy made of this software must include this copyright notice +;in full. +; +;2. I have made no warrantee or representation that the operation of +;this software will be error-free, and I am under no obligation to +;provide any services, by way of maintenance, update, or otherwise. +; +;3. In conjunction with products arising from the use of this +;material, there shall be no use of my name in any advertising, +;promotional, or sales literature without prior written consent in +;each case. + +;;;; Newton's Method explained in: +;;; D. E. Knuth, "The Art of Computer Programming", Vol 2 / +;;; Seminumerical Algorithms, Reading Massachusetts, Addison-Wesley +;;; Publishing Company, 2nd Edition, p. 510 + +(define (newton:find-integer-root f df/dx x_0) + (let loop ((x x_0) (fx (f x_0))) + (cond + ((zero? fx) x) + (else + (let ((df (df/dx x))) + (cond + ((zero? df) #f) ; stuck at local min/max + (else + (let* ((delta (quotient (+ fx (quotient df 2)) df)) + (next-x (cond ((not (zero? delta)) (- x delta)) + ((positive? fx) (- x 1)) + (else (- x -1)))) + (next-fx (f next-x))) + (cond ((>= (abs next-fx) (abs fx)) x) + (else (loop next-x next-fx))))))))))) + +(define (integer-sqrt y) + (newton:find-integer-root (lambda (x) (- (* x x) y)) + (lambda (x) (* 2 x)) + (ash 1 (quotient (integer-length y) 2)))) + +(define (newton:find-root f df/dx x_0 prec) + (if (and (negative? prec) (integer? prec)) + (let loop ((x x_0) (fx (f x_0)) (count prec)) + (cond ((zero? count) x) + (else (let ((df (df/dx x))) + (cond ((zero? df) #f) ; stuck at local min/max + (else (let* ((next-x (- x (/ fx df))) + (next-fx (f next-x))) + (cond ((= next-x x) x) + ((> (abs next-fx) (abs fx)) #f) + (else (loop next-x next-fx + (+ 1 count))))))))))) + (let loop ((x x_0) (fx (f x_0))) + (cond ((< (abs fx) prec) x) + (else (let ((df (df/dx x))) + (cond ((zero? df) #f) ; stuck at local min/max + (else (let* ((next-x (- x (/ fx df))) + (next-fx (f next-x))) + (cond ((= next-x x) x) + ((> (abs next-fx) (abs fx)) #f) + (else (loop next-x next-fx)))))))))))) + +;;; H. J. Orchard, "The Laguerre Method for Finding the Zeros of +;;; Polynomials", IEEE Transactions on Circuits and Systems, Vol. 36, +;;; No. 11, November 1989, pp 1377-1381. + +(define (laguerre:find-root f df/dz ddf/dz^2 z_0 prec) + (if (and (negative? prec) (integer? prec)) + (let loop ((z z_0) (fz (f z_0)) (count prec)) + (cond ((zero? count) z) + (else + (let* ((df (df/dz z)) + (ddf (ddf/dz^2 z)) + (disc (sqrt (- (* df df) (* fz ddf))))) + (if (zero? disc) + #f + (let* ((next-z + (- z (/ fz (if (negative? (+ (* (real-part df) + (real-part disc)) + (* (imag-part df) + (imag-part disc)))) + (- disc) disc)))) + (next-fz (f next-z))) + (cond ((>= (magnitude next-fz) (magnitude fz)) z) + (else (loop next-z next-fz (+ 1 count)))))))))) + (let loop ((z z_0) (fz (f z_0)) (delta-z #f)) + (cond ((< (magnitude fz) prec) z) + (else + (let* ((df (df/dz z)) + (ddf (ddf/dz^2 z)) + (disc (sqrt (- (* df df) (* fz ddf))))) + (print 'disc disc) + (if (zero? disc) + #f + (let* ((next-z + (- z (/ fz (if (negative? (+ (* (real-part df) + (real-part disc)) + (* (imag-part df) + (imag-part disc)))) + (- disc) disc)))) + (next-delta-z (magnitude (- next-z z)))) + (print 'next-z next-z ) + (print '(f next-z) (f next-z)) + (print 'delta-z delta-z 'next-delta-z next-delta-z) + (cond ((zero? next-delta-z) z) + ((and delta-z (>= next-delta-z delta-z)) z) + (else + (loop next-z (f next-z) next-delta-z))))))))))) + +(define (laguerre:find-polynomial-root deg f df/dz ddf/dz^2 z_0 prec) + (if (and (negative? prec) (integer? prec)) + (let loop ((z z_0) (fz (f z_0)) (count prec)) + (cond ((zero? count) z) + (else + (let* ((df (df/dz z)) + (ddf (ddf/dz^2 z)) + (tmp (* (+ deg -1) df)) + (sqrt-H (sqrt (- (* tmp tmp) (* deg (+ deg -1) fz ddf)))) + (df+sqrt-H (+ df sqrt-H)) + (df-sqrt-H (- df sqrt-H)) + (next-z + (- z (/ (* deg fz) + (if (>= (magnitude df+sqrt-H) + (magnitude df-sqrt-H)) + df+sqrt-H + df-sqrt-H))))) + (loop next-z (f next-z) (+ 1 count)))))) + (let loop ((z z_0) (fz (f z_0))) + (cond ((< (magnitude fz) prec) z) + (else + (let* ((df (df/dz z)) + (ddf (ddf/dz^2 z)) + (tmp (* (+ deg -1) df)) + (sqrt-H (sqrt (- (* tmp tmp) (* deg (+ deg -1) fz ddf)))) + (df+sqrt-H (+ df sqrt-H)) + (df-sqrt-H (- df sqrt-H)) + (next-z + (- z (/ (* deg fz) + (if (>= (magnitude df+sqrt-H) + (magnitude df-sqrt-H)) + df+sqrt-H + df-sqrt-H))))) + (loop next-z (f next-z)))))))) |