1 |
;;;"root.scm" Newton's and Laguerre's methods for finding roots. |
2 |
;Copyright (C) 1996, 1997 Aubrey Jaffer |
3 |
; |
4 |
;Permission to copy this software, to modify it, to redistribute it, |
5 |
;to distribute modified versions, and to use it for any purpose is |
6 |
;granted, subject to the following restrictions and understandings. |
7 |
; |
8 |
;1. Any copy made of this software must include this copyright notice |
9 |
;in full. |
10 |
; |
11 |
;2. I have made no warranty or representation that the operation of |
12 |
;this software will be error-free, and I am under no obligation to |
13 |
;provide any services, by way of maintenance, update, or otherwise. |
14 |
; |
15 |
;3. In conjunction with products arising from the use of this |
16 |
;material, there shall be no use of my name in any advertising, |
17 |
;promotional, or sales literature without prior written consent in |
18 |
;each case. |
19 |
|
20 |
(require 'logical) |
21 |
|
22 |
;;;; Newton's Method explained in: |
23 |
;;; D. E. Knuth, "The Art of Computer Programming", Vol 2 / |
24 |
;;; Seminumerical Algorithms, Reading Massachusetts, Addison-Wesley |
25 |
;;; Publishing Company, 2nd Edition, p. 510 |
26 |
;@ |
27 |
(define (newton:find-integer-root f df/dx x_0) |
28 |
(let loop ((x x_0) (fx (f x_0))) |
29 |
(cond |
30 |
((zero? fx) x) |
31 |
(else |
32 |
(let ((df (df/dx x))) |
33 |
(cond |
34 |
((zero? df) #f) ; stuck at local min/max |
35 |
(else |
36 |
(let* ((delta (quotient (+ fx (quotient df 2)) df)) |
37 |
(next-x (cond ((not (zero? delta)) (- x delta)) |
38 |
((positive? fx) (- x 1)) |
39 |
(else (- x -1)))) |
40 |
(next-fx (f next-x))) |
41 |
(cond ((>= (abs next-fx) (abs fx)) x) |
42 |
(else (loop next-x next-fx))))))))))) |
43 |
|
44 |
;;(define (integer-sqrt y) |
45 |
;; (newton:find-integer-root (lambda (x) (- (* x x) y)) |
46 |
;; (lambda (x) (* 2 x)) |
47 |
;; (ash 1 (quotient (integer-length y) 2)))) |
48 |
|
49 |
;@ |
50 |
(define (newton:find-root f df/dx x_0 prec) |
51 |
(if (and (negative? prec) (integer? prec)) |
52 |
(let loop ((x x_0) (fx (f x_0)) (count prec)) |
53 |
(cond ((zero? count) x) |
54 |
(else (let ((df (df/dx x))) |
55 |
(cond ((zero? df) #f) ; stuck at local min/max |
56 |
(else (let* ((next-x (- x (/ fx df))) |
57 |
(next-fx (f next-x))) |
58 |
(cond ((= next-x x) x) |
59 |
((> (abs next-fx) (abs fx)) #f) |
60 |
(else (loop next-x next-fx |
61 |
(+ 1 count))))))))))) |
62 |
(let loop ((x x_0) (fx (f x_0))) |
63 |
(cond ((< (abs fx) prec) x) |
64 |
(else (let ((df (df/dx x))) |
65 |
(cond ((zero? df) #f) ; stuck at local min/max |
66 |
(else (let* ((next-x (- x (/ fx df))) |
67 |
(next-fx (f next-x))) |
68 |
(cond ((= next-x x) x) |
69 |
((> (abs next-fx) (abs fx)) #f) |
70 |
(else (loop next-x next-fx)))))))))))) |
71 |
|
72 |
;;; H. J. Orchard, "The Laguerre Method for Finding the Zeros of |
73 |
;;; Polynomials", IEEE Transactions on Circuits and Systems, Vol. 36, |
74 |
;;; No. 11, November 1989, pp 1377-1381. |
75 |
;@ |
76 |
(define (laguerre:find-root f df/dz ddf/dz^2 z_0 prec) |
77 |
(if (and (negative? prec) (integer? prec)) |
78 |
(let loop ((z z_0) (fz (f z_0)) (count prec)) |
79 |
(cond ((zero? count) z) |
80 |
(else |
81 |
(let* ((df (df/dz z)) |
82 |
(ddf (ddf/dz^2 z)) |
83 |
(disc (sqrt (- (* df df) (* fz ddf))))) |
84 |
(if (zero? disc) |
85 |
#f |
86 |
(let* ((next-z |
87 |
(- z (/ fz (if (negative? (+ (* (real-part df) |
88 |
(real-part disc)) |
89 |
(* (imag-part df) |
90 |
(imag-part disc)))) |
91 |
(- disc) disc)))) |
92 |
(next-fz (f next-z))) |
93 |
(cond ((>= (magnitude next-fz) (magnitude fz)) z) |
94 |
(else (loop next-z next-fz (+ 1 count)))))))))) |
95 |
(let loop ((z z_0) (fz (f z_0)) (delta-z #f)) |
96 |
(cond ((< (magnitude fz) prec) z) |
97 |
(else |
98 |
(let* ((df (df/dz z)) |
99 |
(ddf (ddf/dz^2 z)) |
100 |
(disc (sqrt (- (* df df) (* fz ddf))))) |
101 |
;;(print 'disc disc) |
102 |
(if (zero? disc) |
103 |
#f |
104 |
(let* ((next-z |
105 |
(- z (/ fz (if (negative? (+ (* (real-part df) |
106 |
(real-part disc)) |
107 |
(* (imag-part df) |
108 |
(imag-part disc)))) |
109 |
(- disc) disc)))) |
110 |
(next-delta-z (magnitude (- next-z z)))) |
111 |
;;(print 'next-z next-z ) |
112 |
;;(print '(f next-z) (f next-z)) |
113 |
;;(print 'delta-z delta-z 'next-delta-z next-delta-z) |
114 |
(cond ((zero? next-delta-z) z) |
115 |
((and delta-z (>= next-delta-z delta-z)) z) |
116 |
(else |
117 |
(loop next-z (f next-z) next-delta-z))))))))))) |
118 |
;@ |
119 |
(define (laguerre:find-polynomial-root deg f df/dz ddf/dz^2 z_0 prec) |
120 |
(if (and (negative? prec) (integer? prec)) |
121 |
(let loop ((z z_0) (fz (f z_0)) (count prec)) |
122 |
(cond ((zero? count) z) |
123 |
(else |
124 |
(let* ((df (df/dz z)) |
125 |
(ddf (ddf/dz^2 z)) |
126 |
(tmp (* (+ deg -1) df)) |
127 |
(sqrt-H (sqrt (- (* tmp tmp) (* deg (+ deg -1) fz ddf)))) |
128 |
(df+sqrt-H (+ df sqrt-H)) |
129 |
(df-sqrt-H (- df sqrt-H)) |
130 |
(next-z |
131 |
(- z (/ (* deg fz) |
132 |
(if (>= (magnitude df+sqrt-H) |
133 |
(magnitude df-sqrt-H)) |
134 |
df+sqrt-H |
135 |
df-sqrt-H))))) |
136 |
(loop next-z (f next-z) (+ 1 count)))))) |
137 |
(let loop ((z z_0) (fz (f z_0))) |
138 |
(cond ((< (magnitude fz) prec) z) |
139 |
(else |
140 |
(let* ((df (df/dz z)) |
141 |
(ddf (ddf/dz^2 z)) |
142 |
(tmp (* (+ deg -1) df)) |
143 |
(sqrt-H (sqrt (- (* tmp tmp) (* deg (+ deg -1) fz ddf)))) |
144 |
(df+sqrt-H (+ df sqrt-H)) |
145 |
(df-sqrt-H (- df sqrt-H)) |
146 |
(next-z |
147 |
(- z (/ (* deg fz) |
148 |
(if (>= (magnitude df+sqrt-H) |
149 |
(magnitude df-sqrt-H)) |
150 |
df+sqrt-H |
151 |
df-sqrt-H))))) |
152 |
(loop next-z (f next-z)))))))) |
153 |
|
154 |
(define (secant:find-root-1 f x0 x1 prec must-bracket?) |
155 |
(letrec ((stop? |
156 |
(cond ((procedure? prec) prec) |
157 |
((and (integer? prec) (negative? prec)) |
158 |
(lambda (x0 f0 x1 f1 count) |
159 |
(>= count (- prec)))) |
160 |
(else |
161 |
(lambda (x0 f0 x1 f1 count) |
162 |
(and (< (abs f0) prec) |
163 |
(< (abs f1) prec)))))) |
164 |
(bracket-iter |
165 |
(lambda (xlo flo glo xhi fhi ghi count) |
166 |
(define (step xnew fnew) |
167 |
(cond ((or (= xnew xlo) |
168 |
(= xnew xhi)) |
169 |
(let ((xmid (+ xlo (* 1/2 (- xhi xlo))))) |
170 |
(if (= xnew xmid) |
171 |
xmid |
172 |
(step xmid (f xmid))))) |
173 |
((positive? fnew) |
174 |
(bracket-iter xlo flo (if glo (* 0.5 glo) 1) |
175 |
xnew fnew #f |
176 |
(+ count 1))) |
177 |
(else |
178 |
(bracket-iter xnew fnew #f |
179 |
xhi fhi (if ghi (* 0.5 ghi) 1) |
180 |
(+ count 1))))) |
181 |
(if (stop? xlo flo xhi fhi count) |
182 |
(if (> (abs flo) (abs fhi)) xhi xlo) |
183 |
(let* ((fflo (if glo (* glo flo) flo)) |
184 |
(ffhi (if ghi (* ghi fhi) fhi)) |
185 |
(del (- (/ fflo (- ffhi fflo)))) |
186 |
(xnew (+ xlo (* del (- xhi xlo)))) |
187 |
(fnew (f xnew))) |
188 |
(step xnew fnew)))))) |
189 |
(let ((f0 (f x0)) |
190 |
(f1 (f x1))) |
191 |
(cond ((<= f0 0 f1) |
192 |
(bracket-iter x0 f0 #f x1 f1 #f 0)) |
193 |
((<= f1 0 f0) |
194 |
(bracket-iter x1 f1 #f x0 f0 #f 0)) |
195 |
(must-bracket? #f) |
196 |
(else |
197 |
(let secant-iter ((x0 x0) |
198 |
(f0 f0) |
199 |
(x1 x1) |
200 |
(f1 f1) |
201 |
(count 0)) |
202 |
(cond ((stop? x0 f0 x1 f1 count) |
203 |
(if (> (abs f0) (abs f1)) x1 x0)) |
204 |
((<= f0 0 f1) |
205 |
(bracket-iter x0 f0 #f x1 f1 #f count)) |
206 |
((>= f0 0 f1) |
207 |
(bracket-iter x1 f1 #f x0 f0 #f count)) |
208 |
((= f0 f1) #f) |
209 |
(else |
210 |
(let ((xnew (+ x0 (* (- (/ f0 (- f1 f0))) (- x1 x0))))) |
211 |
(secant-iter x1 f1 xnew (f xnew) (+ count 1))))))))))) |
212 |
;@ |
213 |
(define (secant:find-root f x0 x1 prec) |
214 |
(secant:find-root-1 f x0 x1 prec #f)) |
215 |
(define (secant:find-bracketed-root f x0 x1 prec) |
216 |
(secant:find-root-1 f x0 x1 prec #t)) |