/[slib]/slib/root.scm
ViewVC logotype

Contents of /slib/root.scm

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.15 - (show annotations) (download)
Mon Apr 11 23:10:39 2011 UTC (12 years, 11 months ago) by jaffer
Branch: MAIN
CVS Tags: slib-3b6, slib-3b7, slib-3b4, slib-3b5, slib-3c1, HEAD
Changes since 1.14: +2 -4 lines
From Bill Schottstaedt:
(secant:find-root-1): Removed unused variable.

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))

savannah-hackers-public@gnu.org
ViewVC Help
Powered by ViewVC 1.1.26