/[sapaclisp]/sapaclisp/basic-math.lisp
ViewVC logotype

Contents of /sapaclisp/basic-math.lisp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1.1.1 - (show annotations) (vendor branch)
Mon May 9 20:24:31 2005 UTC (8 years, 11 months ago) by mmommer
Branch: MAIN, T
CVS Tags: initial_checkin, HEAD
Changes since 1.1: +0 -0 lines
Initial Checkin
1 ;;;-*- Mode: LISP; Package: :SAPA; Syntax: COMMON-LISP -*-
2 ;-----------------------------------------------------------------------------
3 ; (c) 1993, Donald B. Percival <dbp@apl.washington.edu>
4 ;
5 ; This code is licensed under the terms of the modified BSD license
6 ; ("sans advertising clause"). See the file COPYING for details.
7 ;
8 ; Comments about this software can be addressed to dbp@apl.washington.edu
9 ;-----------------------------------------------------------------------------
10 ;
11 ; basic-math.lisp
12 ;
13 ; a collection of Lisp functions for certain basic mathematical operations ...
14 ; Note: before compiling and loading basic-math.lisp,
15 ; you should compile and load
16 ; sapa-package.lisp
17 ;
18 ;-----------------------------------------------------------------------------
19 (if (not (find-package :SAPA))
20 (defpackage :SAPA (:USE :COMMON-LISP)))
21
22 (in-package :SAPA)
23
24 (export '(;;; functions related to the gamma function and its derivatives ...
25 log-of-gamma
26 factorial
27 digamma
28 trigamma
29
30 ;;; functions for dealing with polynomials ...
31 evaluate-polynomial
32 multiply-2-polynomials
33 multiply-polynomials
34 zeros-of-polynomial
35
36 ;;; functions for finding the root of a function ...
37 Newton-Raphson
38 bisection-with-Newton-Raphson
39 secant-method
40
41 ;;; functions for numerical integation ...
42 simple-numerical-integration
43 Gauss-Legendre-quadrature
44 ))
45
46 ;-----------------------------------------------------------------------------
47 ;;; used by the next set of functions ...
48 (defparameter +coefficents-for-log-of-gamma+
49 (list 76.18009172947146d0 -86.50532032941677d0 24.01409824083091d0
50 -1.231739572450155d0 0.1208650973866179d-2 -0.5395239384953d-5))
51
52 ;-----------------------------------------------------------------------------
53 ;;; used by the next set of functions ...
54 (defconstant +log-pi+ (log pi))
55 (defconstant +euler-constant+ 0.577215664901532860606512d0)
56
57 ;-----------------------------------------------------------------------------
58 ;-----------------------------------------------------------------------------
59 ;;; The functions log-of-gamma
60 ;;; factorial
61 ;;; digamma
62 ;;; trigamma
63 ;;; compute various quantities that are related to the gamma function
64 ;;; (log-of-gamma, factorial), its first derivative (digamma) and
65 ;;; its second derivative (trigamma).
66 ;-----------------------------------------------------------------------------
67 ;-----------------------------------------------------------------------------
68 (defun log-of-gamma (xx)
69 "given xx (a real or complex valued number
70 whose real part is greater than 0),
71 returns the log (base e) of gamma(xx)
72 ---
73 Note: based upon the discussion in Section 6.1,
74 Numerical Recipes, Second Edition"
75 (assert (plusp (realpart xx)))
76 (if (< (realpart xx) 1)
77 ;;; use reflection formula (6.1.4)
78 (let ((1-xx (- 1.0 xx)))
79 (- (+ +log-pi+ (log 1-xx))
80 (+ (log-of-gamma (1+ 1-xx)) (log (sin (* pi 1-xx))))))
81 ;;; since Re(xx) > 1, use approximation due to Lanczos
82 (let* ((x (1- xx))
83 (tmp-1 (+ x 5.5))
84 (tmp-2 (- (* (+ x 0.5) (log tmp-1)) tmp-1))
85 (ser 1.000000000190015d0))
86 (dolist (a-coefficient +coefficents-for-log-of-gamma+)
87 (incf x)
88 (incf ser (/ a-coefficient x)))
89 (+ tmp-2 (log (* 2.5066282746310005 ser))))))
90
91 ;;; (log-of-gamma 6.0) ;==> 4.787491742782046
92 ;;; above should agree with
93 ;;; (log (* 2 3 4 5)) ;==> 4.787491742782046
94 ;;; and it does!
95 ;;; (log-of-gamma 0.5) ;==> 0.5723649429246563
96 ;;; above should agree with
97 ;;; (* 0.5 (log pi)) ;==> 0.5723649429247001
98 ;;; and it does!
99 ;;; (exp (log-of-gamma 1.755)) ;==> 0.9202092223790562
100 ;;; above should agree with entry in row 2, column 2,
101 ;;; page 270 of Abramowitz and Stegun, which has
102 ;;; 0.9202092224
103 ;;; and it does!
104
105 ;-----------------------------------------------------------------------------
106 (defun factorial (k)
107 "given an integer k, returns k!"
108 (assert (and (integerp k) (not (minusp k))))
109 (cond
110 ((> k 15) ; arbitrary choice ...
111 (exp (log-of-gamma (1+ k))))
112 ((or (= k 1) (= k 0))
113 1)
114 (t
115 (* k (factorial (1- k))))))
116
117 #|
118 (factorial 0) ;==> 1
119 (factorial 1) ;==> 1
120 (factorial 6) ;==> 720
121 (factorial 25) ;==> 1.5511210043610457E+25
122 |#
123
124 ;-----------------------------------------------------------------------------
125 (defun digamma
126 (x
127 &key
128 (x-recursion 8.5)
129 (x-small 1.0e-5))
130 "given
131 [1] x (required)
132 ==> a positive number
133 [2] x-recursion (keyword; 8.5)
134 ==> for noninteger x,
135 if x-small < x < x-recursion,
136 recursive formula is used
137 [3] x-small (keyword; 1.0e-5)
138 ==> if x <= x-small,
139 small x approximation used (default )
140 returns
141 [1] value of digamma (psi, derivative of log(gamma)) function at x
142 ---
143 Note: see Abramowitz and Stegun's equation 6.3.2;
144 expansion 6.3.18 plus recurrence 6.3.5;
145 the small x formula is Equation (5) from
146 Algorithm AS 103 -- Psi (Digamma) Function -- by J. M. Bernardo
147 in Applied Statistics, Vol. 25, No. 3, 1976, pp. 315--317;
148 note that better accuracy can be obtained by increasing
149 x-recursion --- 10.0 to 100.0 are probably useful values"
150 (assert (plusp x))
151 (cond
152 ((integerp x)
153 (let ((sum (- +euler-constant+)))
154 (dotimes (i (1- x) sum)
155 (incf sum (/ (1+ i))))))
156 ((<= x x-small)
157 (- (+ +euler-constant+ (/ x))))
158 ((< x x-recursion)
159 (- (digamma (1+ x) :x-recursion x-recursion) (/ x)))
160 (t
161 (let* ((x2 (* x x))
162 (x4 (* x2 x2)))
163 (- (log x)
164 (/ (* 2.0 x))
165 (/ (* 12.0 x2))
166 (/ (* -120.0 x4))
167 (/ (* 252.0 x2 x4)))))))
168
169 #|
170 (let ((x*1000 1500)
171 x)
172 (dotimes (i 15)
173 (setf x (float (/ x*1000 1000)))
174 (format t "~&~5,3F ~12,10F" x (digamma x))
175 (incf x*1000 5)))
176 ;==>
177 1.500 0.0364899738
178 1.505 0.0411536541
179 1.510 0.0457967894
180 1.515 0.0504195526
181 1.520 0.0550221144
182 1.525 0.0596046437
183 1.530 0.0641673072
184 1.535 0.0687102696
185 1.540 0.0732336935
186 1.545 0.0777377398
187 1.550 0.0822225674
188 1.555 0.0866883332
189 1.560 0.0911351924
190 1.565 0.0955632983
191 1.570 0.0999728023
192 ;;; good agreement with 4th column of top of page 269, Abramowitz and Stegun
193 |#
194
195 ;-----------------------------------------------------------------------------
196 (defun trigamma
197 (x
198 &key
199 (x-recursion 2.0))
200 "given
201 [1] x (required)
202 ==> a positive number
203 [2] x-recursion (keyword; 2.0)
204 ==> if x < x-recursion,
205 recursive formula is used
206 returns
207 [1] value of trigamma function at x
208 ---
209 Note: expansion 6.4.12 plus recurrence 6.4.6
210 of Abramowitz and Stegun;
211 better accuracy can be obtained by increasing
212 x-recursion ---10.0 to 30.0 is probably a useful upper limit"
213 (assert (plusp x))
214 (if (< x x-recursion)
215 (+ (trigamma (1+ x) :x-recursion x-recursion) (/ (* x x)))
216 (let* ((x2 (* x x))
217 (x3 (* x x2))
218 (x5 (* x2 x3))
219 (x7 (* x2 x5)))
220 (+ (/ x)
221 (/ (* 2.0 x2))
222 (/ (* 6.0 x3))
223 (/ (* -30.0 x5))
224 (/ (* 42.0 x7))
225 (/ (* -30.0 x2 x7))))))
226 #|
227 (let ((x*1000 1500)
228 x)
229 (dotimes (i 15)
230 (setf x (float (/ x*1000 1000)))
231 (format t "~&~5,3F ~12,10F" x (trigamma x))
232 (incf x*1000 5)))
233 ;==>
234 1.500 0.9348000492
235 1.505 0.9306736517
236 1.510 0.9265820503
237 1.515 0.9225248209
238 1.520 0.9185015462
239 1.525 0.9145118155
240 1.530 0.9105552242
241 1.535 0.9066313746
242 1.540 0.9027398746
243 1.545 0.8988803384
244 1.550 0.8950523864
245 1.555 0.8912556443
246 1.560 0.8874897441
247 1.565 0.8837543230
248 1.570 0.8800490239
249 ;;; 5 place agreement with 5th column of top of page 269, Abramowitz and Stegun
250 |#
251
252 ;-----------------------------------------------------------------------------
253 ;-----------------------------------------------------------------------------
254 ;;; The functions evaluate-polynomial
255 ;;; multiply-2-polynomials
256 ;;; multiply-polynomials
257 ;;; zeros-of-polynomial
258 ;;; evaluate, multiply and find the roots of polynomials whose coefficients
259 ;;; are represented by sequences of numbers.
260 ;-----------------------------------------------------------------------------
261 ;-----------------------------------------------------------------------------
262 (defun evaluate-polynomial
263 (polynomial
264 z
265 &key
266 (degree-of-polynomial (1- (length polynomial)))
267 (number-of-derivatives 0)
268 (complex-values (make-array (1+ number-of-derivatives)))
269 (bounds-to-be-computed nil)
270 (bounds (if bounds-to-be-computed
271 (make-array (1+ number-of-derivatives)))))
272 "given
273 [1] polynomial (required)
274 ==> sequence of length n+1 with real or
275 complex-valued numbers giving the n+1
276 coefficients of a polynomial of nth order;
277 (elt polynomial i) = coefficient of z^(n-i)
278 [2] z (required)
279 ==> complex point at which the polynomial
280 to be evaluated
281 [3] degree-of-polynomial (keyword; (1- (length polynomial)))
282 ==> degree of polynomial
283 [4] number-of-derivatives (keyword; 0)
284 ==> number of derivatives to be evaluated
285 [5] complex-values (keyword; array of length (1+ number-of-derivatives))
286 <== sequence with values of polynomial
287 and derivatives
288 [6] bounds-to-be-computed (keyword; nil)
289 ==> if t, error bounds are computed
290 for values in complex-values
291 [7] bounds (keyword; array of length (1+ number-of-derivatives))
292 <== sequence with error bounds
293 returns
294 [1] the sequence complex-values containing the value
295 of the polynomial and its derivatives
296 [2] bounds, a sequence containing optional error bounds
297 ---
298 Note: this is a Lisp version of cmlib routine cpevl;
299 d1 was specified in cpevl as (expt 2 (- 1 (i1mach 11))),
300 where (i1mach 11) ==> THE NUMBER OF BASE-2 DIGITS (SINGLE PRECISION).
301 I have taken this to be equivalent to machine epsilon.
302 If this is in fact NOT the case, then the optional error bounds
303 might not be computed correctly."
304 (let ((d1 single-float-epsilon)
305 ci cim1)
306 (dotimes (j (1+ degree-of-polynomial))
307 (dotimes (i (min (1+ number-of-derivatives)
308 (1+ (- degree-of-polynomial j))))
309 (setf ci (if (plusp j) (elt complex-values i) 0.0))
310 (setf cim1 (if (plusp i)
311 (elt complex-values (1- i)) (elt polynomial j)))
312 (setf (elt complex-values i) (+ cim1 (* z ci)))
313 (if bounds-to-be-computed
314 (let* ((bi (if (plusp j) (elt bounds i) 0.0))
315 (bim1 (if (plusp i) (elt bounds (1- i)) 0.0))
316 (tf (+ bi (* (complex-of-absolutes ci)
317 (+ (* 3.0 d1) (* 4.0 d1 d1)))))
318 (r (realpart (* (complex-of-absolutes z)
319 (complex (realpart tf) (- (imagpart tf))))))
320 (s (imagpart (* (complex-of-absolutes z) tf))))
321 (setf (elt bounds i)
322 (if (plusp j)
323 (* (1+ (* 8 d1))
324 (+ bim1
325 (* d1 (complex-of-absolutes cim1))
326 (complex r s)))
327 0.0)))))))
328 (values complex-values bounds))
329
330 #|
331 (evaluate-polynomial #(1 2 3) 1)
332 ;==> #(6.0)
333 ; nil
334 (evaluate-polynomial #(1 2 3) 2)
335 ;==> #(11.0)
336 ; nil
337 |#
338
339 ;-----------------------------------------------------------------------------
340 (defun multiply-2-polynomials
341 (polynomial-1
342 polynomial-2
343 &key
344 (degree-of-1
345 (1- (length polynomial-1)))
346 (degree-of-2
347 (1- (length polynomial-2)))
348 (product-polynomial
349 (make-array (+ degree-of-1 degree-of-2 1))))
350 "given
351 [1] polynomial-1 (required)
352 ==> sequence of length n+1 with real or
353 complex-valued numbers giving the n+1
354 coefficients of a polynomial of nth order;
355 (elt polynomial-1 i) = coefficient of z^(n-i)
356 [2] polynomial-2 (required)
357 ==> another sequence representing another polynomial
358 [3] degree-of-1 (keyword; (1- (length polynomial-1)))
359 ==> degree of first polynomial
360 [4] degree-of-2 (keyword; (1- (length polynomial-2)))
361 ==> degree of second polynomial
362 [5] product-polynomial (keyword; array of length (+ degree-of-1 degree-of-2 1))
363 <== a sequence with product of two polynomials
364 returns
365 [1] product-polynomial, an sequence with coefficients
366 of polynomial given by the product of
367 polynomial-1 and polynomial-2
368 ---
369 Note: this routine works for a polynomial p represented either by
370 p(0) + p(1)*z + ... + p(degp)*z^degp
371 or by
372 p(0)*z**degp + p(1)*z**(degp-1) + ... + p(degp)"
373 (let (k)
374 (dotimes (i (+ degree-of-1 degree-of-2 1) product-polynomial)
375 (setf (elt product-polynomial i) 0.0)
376 (dotimes (j (1+ degree-of-1))
377 (setf k (- i j))
378 (if (and (>= k 0) (<= k degree-of-2))
379 (incf (elt product-polynomial i)
380 (* (elt polynomial-1 j) (elt polynomial-2 k))))))))
381
382 #|
383 (multiply-2-polynomials #(1 2 3) #(4 3 2 1))
384 ;==>#(4.0 11.0 20.0 14.0 8.0 3.0)
385 (multiply-2-polynomials #(4 3 2 1) #(1 2 3))
386 ;==>#(4.0 11.0 20.0 14.0 8.0 3.0)
387 (multiply-2-polynomials #(4 3 2 1) #(2))
388 ;==>#(8.0 6.0 4.0 2.0)
389 |#
390
391 ;-----------------------------------------------------------------------------
392 (defun multiply-polynomials (&rest polys)
393 "given
394 [1] sequences representing any number of polynomials (required)
395 returns
396 [1] a sequence representing their product
397 ---
398 Note: this function was written by Andrew G. Bruce"
399 (let ((n (length polys)))
400 (cond ((= n 0) polys)
401 ((= n 1) (first polys))
402 ((> n 1) (reduce #'multiply-2-polynomials polys)))))
403
404 #|
405 (multiply-2-polynomials #(1 2) #(4 3))
406 ;==> #(4.0 11.0 6.0)
407 (multiply-2-polynomials #(4.0 11.0 6.0) #(1 2 3))
408 ;==> #(4.0 19.0 40.0 45.0 18.0)
409 (multiply-polynomials #(1 2) #(4 3) #(1 2 3))
410 ;==> #(4.0 19.0 40.0 45.0 18.0)
411 (multiply-polynomials #(1) #(4) #(2))
412 ;==> #(8.0)
413 |#
414
415 ;-----------------------------------------------------------------------------
416 (defun zeros-of-polynomial
417 (polynomial
418 &key
419 (degree-of-polynomial (1- (length polynomial)))
420 (the-roots (make-array degree-of-polynomial))
421 (maximum-number-of-iterations 25))
422 "given
423 [1] polynomial (required)
424 ==> sequence with coefficients of a polynomial;
425 (elt polynomial 0) must be nonzero
426 [2] degree-of-polynomial (keyword; (1- (length polynomial)))
427 ==> degree of polynomial
428 [3] the-roots (keyword; array of length degree-of-polynomial)
429 <== number of derivatives to be evaluated
430 [4] maximum-number-of-iterations (keyword; 25)
431 ==> maximum number of iterations
432 returns
433 [1] t or nil, where t indicates that all went well,
434 whereas nil indicates that convergence did not occur
435 at end of specificed number of iterations
436 [2] the-roots, a vector with the required roots
437 of the polynomial"
438 (cond
439 ;;; if this is a first degree, we're done
440 ((= degree-of-polynomial 1)
441 (setf (elt the-roots 0) (- (/ (elt polynomial 1) (elt polynomial 0))))
442 (values t the-roots)) ;t = all went well
443 ;;; If constant term is zero, we know one root and can get
444 ;;; the others by dealing with a polynomial of one less degree.
445 ((zerop (abs (elt polynomial degree-of-polynomial)))
446 (setf (elt the-roots (1- degree-of-polynomial)) 0.0)
447 (zeros-of-polynomial polynomial
448 :degree-of-polynomial
449 (1- degree-of-polynomial)
450 :the-roots
451 the-roots
452 :maximum-number-of-iterations
453 maximum-number-of-iterations))
454 ;;; Here we have to do some honest work: we first generate
455 ;;; some initial estimates for the roots and then the final estimates.
456 (t
457 (let ((initial-estimates (make-array degree-of-polynomial))
458 (scratch (make-array (* 2 (1+ degree-of-polynomial)))))
459 (initial-estimates-for-zeros-of-polynomial
460 polynomial
461 :degree-of-polynomial degree-of-polynomial
462 :initial-estimates initial-estimates
463 :scratch scratch)
464 ;;; If this function return nil, convergence was NOT found
465 (if (zeros-of-polynomial-with-initial-estimates
466 polynomial initial-estimates
467 :degree-of-polynomial degree-of-polynomial
468 :scratch scratch
469 :final-results the-roots
470 :maximum-number-of-iterations maximum-number-of-iterations)
471 (values t the-roots) ;convergence occurred
472 (values nil the-roots) ;did not converge
473 )))))
474
475 #|
476 (zeros-of-polynomial #(1 -2))
477 ;==> t
478 ; #(2)
479 (zeros-of-polynomial #(1 -2 0))
480 ;==> t
481 ; #(2 0.0)
482 (multiply-polynomials #(1 -2) #(1 -3) #(1 5))
483 ;==> #(1.0 0.0 -19.0 30.0)
484 (zeros-of-polynomial #(1.0 0.0 -19.0 30.0))
485 ;==>t
486 ; #(#c(2.9999999999999996 -2.2190775908547207E-22)
487 #c(-5.0 1.9700224365534028E-19)
488 #c(1.9999999999999998 2.219061434983382E-22))
489 |#
490
491 ;-----------------------------------------------------------------------------
492 ;-----------------------------------------------------------------------------
493 ;;; The functions Newton-Raphson
494 ;;; bisection-with-Newton-Raphson
495 ;;; secant-method
496 ;;; can be used to find the zero (root) of a function in a specfied interval.
497 ;-----------------------------------------------------------------------------
498 ;-----------------------------------------------------------------------------
499 (defun Newton-Raphson
500 (f
501 f-prime
502 x-left
503 x-right
504 &key
505 (accuracy (* 10.0 single-float-epsilon))
506 (maximum-number-of-iterations 20)
507 (prevent-bracket-jumping-p t))
508 "given
509 [1] f (required)
510 ==> a function with a single argument
511 [2] f-prime (required)
512 ==> another function with a single argument,
513 this one being the first derivative of f
514 [3] x-left (required)
515 ==> left-hand bracket for the desired root;
516 i.e., left-hand bracket <= desired root
517 [4] x-right (required)
518 ==> right-hand bracket for the desired root;
519 i.e., desired root <= right-hand bracket
520 [5] accuracy (keyword; (* 10.0 single-float-epsilon))
521 ==> desired relative accuracy for computed rood
522 [6] maximum-number-of-iterations (keyword; 20)
523 ==> maximum number of iterations
524 [7] prevent-bracket-jumping-p (keyword; t)
525 ==> if t, allows Newton-Raphson to continue
526 if it jumps out of the interval
527 [x-left, x-right];
528 if nil, jumping out of the interval
529 causes an error to be signaled
530 returns
531 [1] a root of f in [x-left, x-right];
532 i.e., a value x in that interval
533 such that f(x) = 0
534 [2] the number of iterations required
535 ---
536 Note: this function is based loosely on rtnewt,
537 Section 9.4, Numerical Recipes, Second Edition"
538 (assert (< x-left x-right))
539 (let ((x (* 0.5 (+ x-left x-right)))
540 delta-x denom-for-accuracy-test)
541 (dotimes (j maximum-number-of-iterations
542 (if (not (cerror "returns solution so far"
543 "exceeding maximum number of iterations"))
544 (values x maximum-number-of-iterations)))
545 (setf delta-x (/ (funcall f x) (funcall f-prime x)))
546 (setf denom-for-accuracy-test (+ (abs x)
547 (abs (decf x delta-x))))
548 (cond
549 (prevent-bracket-jumping-p
550 (if (< x x-left) (setf x x-left))
551 (if (> x x-right) (setf x x-right))
552 (if (< (/ (abs delta-x) denom-for-accuracy-test) accuracy)
553 (return (values x (1+ j)))))
554 ((<= x-left x x-right)
555 (if (< (/ (abs delta-x) denom-for-accuracy-test) accuracy)
556 (return (values x (1+ j)))))
557 (t
558 (error "jumped out of brackets")
559 )))))
560
561 #|
562 (defun exp-2 (x) (- (exp x) 2.0))
563
564 (Newton-Raphson
565 #'exp-2
566 #'exp ;first derivative of exp-2
567 0.0
568 20.0)
569
570 ;==> 0.6931471805599453 ;the root
571 ; 15 ;number of iterations needed
572 (exp-2 0.6931471805599453)
573 ;==> 0.0
574 |#
575
576 ;-----------------------------------------------------------------------------
577 (defun bisection-with-Newton-Raphson
578 (f
579 f-prime
580 x-left
581 x-right
582 &key
583 (accuracy (* 10.0 single-float-epsilon))
584 (maximum-number-of-iterations 100))
585 "given
586 [1] f (required)
587 ==> a function with a single argument
588 [2] f-prime (required)
589 ==> another function with a single argument,
590 this one being the first derivative of f
591 [3] x-left (required)
592 ==> left-hand bracket for the desired root;
593 i.e., left-hand bracket <= desired root
594 [4] x-right (required)
595 ==> right-hand bracket for the desired root;
596 i.e., desired root <= right-hand bracket
597 [5] accuracy (keyword; (* 10.0 single-float-epsilon))
598 ==> desired relative accuracy for computed rood
599 [6] maximum-number-of-iterations (keyword; 100)
600 ==> maximum number of iterations
601 returns
602 [1] a root of f in [x-left, x-right];
603 i.e., a value x in that interval
604 such that f(x) = 0
605 [2] the number of iterations required
606 ---
607 Note: this function is based loosely on rtsafe,
608 Section 9.4, Numerical Recipes, Second Edition"
609 (let ((f-low (funcall f x-left))
610 (f-high (funcall f x-right))
611 df f-new
612 x-low x-high
613 rtsafe dxold dx temp)
614 (when (>= (* f-low f-high) 0.0)
615 (cond ((zerop f-low)
616 (values x-left 0))
617 ((zerop f-high)
618 (values x-right 0))
619 (t
620 (error "root not bracketed"))))
621 (cond ((< f-low 0.0)
622 (setf x-low x-left
623 x-high x-right))
624 (t
625 (setf x-high x-left
626 x-low x-right)))
627 (setf rtsafe (* 0.5 (+ x-left x-right))
628 dxold (abs (- x-right x-left))
629 dx dxold
630 f-new (funcall f rtsafe)
631 df (funcall f-prime rtsafe))
632 (dotimes (j maximum-number-of-iterations
633 (if (not (cerror "returns solution so far"
634 "exceeding maximum number of iterations"))
635 (values rtsafe maximum-number-of-iterations)))
636 (cond ((or (>= (* (- (* (- rtsafe x-high) df) f-new)
637 (- (* (- rtsafe x-low) df) f-new))
638 0.0)
639 (> (abs (* 2.0 f-new)) (abs (* dxold df))))
640 (setf dxold dx
641 dx (* 0.5 (- x-high x-low))
642 rtsafe (+ x-low dx))
643 (when (= x-low rtsafe) (return (values rtsafe (1+ j)))))
644 (t
645 (setf dxold dx
646 dx (/ f-new df)
647 temp rtsafe
648 rtsafe (- rtsafe dx))
649 (when (= temp rtsafe) (return (values rtsafe (1+ j))))))
650 (when (< (abs dx) accuracy) (return (values rtsafe (1+ j))))
651 (setf f-new (funcall f rtsafe)
652 df (funcall f-prime rtsafe))
653 (if (< f-new 0.0)
654 (setf x-low rtsafe)
655 (setf x-high rtsafe)))))
656
657 #|
658 (defun exp-2 (x)
659 (- (exp x) 2.0))
660
661 (bisection-with-Newton-Raphson
662 #'exp-2 #'exp 0.0 20.0
663 :maximum-number-of-iterations 100)
664 ;==> 0.6931471805599447
665 ; 60
666 ;;; good agreement with Newton-Raphson, but takes a lot longer!
667 |#
668
669 ;-----------------------------------------------------------------------------
670 (defun secant-method
671 (f
672 x-left
673 x-right
674 &key
675 (accuracy (* 10.0 single-float-epsilon))
676 (maximum-number-of-iterations 50))
677 "given
678 [1] f (required)
679 ==> a function with a single argument
680 [2] x-left (required)
681 ==> left-hand bracket for the desired root;
682 i.e., left-hand bracket <= desired root
683 [3] x-right (required)
684 ==> right-hand bracket for the desired root;
685 i.e., desired root <= right-hand bracket
686 [4] accuracy (keyword; (* 10.0 single-float-epsilon))
687 ==> desired relative accuracy for computed rood
688 [5] maximum-number-of-iterations (keyword; 50)
689 ==> maximum number of iterations
690 returns
691 [1] a root of f in [x-left, x-right];
692 i.e., a value x in that interval
693 such that f(x) = 0
694 [2] the number of iterations required
695 ---
696 Note: this function is based loosely on rtsec,
697 Section 9.2, Numerical Recipes, Second Edition"
698 (let ((f-left (funcall f x-left))
699 (f-right (funcall f x-right))
700 x-mid f-mid approx-f-mid approx-f-prime-mid delta-x x-new f-new
701 denom-for-accuracy-test)
702 (dotimes (j maximum-number-of-iterations
703 (if (not (cerror "returns solution so far"
704 "exceeding maximum number of iterations"))
705 (values x-new maximum-number-of-iterations)))
706 (setf x-mid (* 0.5 (+ x-left x-right))
707 f-mid (funcall f x-mid)
708 approx-f-mid (* 0.5 (+ f-left f-right))
709 approx-f-prime-mid (/ (- f-right f-left)
710 (- x-right x-left))
711 delta-x (/ approx-f-mid approx-f-prime-mid)
712 x-new (- x-mid delta-x)
713 f-new (funcall f x-new))
714 (setf denom-for-accuracy-test (+ (abs x-mid) (abs x-new)))
715 (if (or (zerop f-new)
716 (< (/ (abs delta-x) denom-for-accuracy-test) accuracy))
717 (return (values x-new (1+ j))))
718 (if (>= (* f-mid f-left) 0)
719 (setf x-left x-mid
720 f-left f-mid))
721 (if (>= (* f-mid f-right) 0)
722 (setf x-right x-mid
723 f-right f-mid))
724 (if (>= (* f-new f-left) 0)
725 (setf x-left x-new
726 f-left f-new)
727 (setf x-right x-new
728 f-right f-new))
729 )))
730
731 #|
732 (defun exp-2 (x)
733 (- (exp x) 2.0))
734
735 (secant-method #'exp-2 0.0 20.0)
736 ;==> 0.6931471805599453
737 ; 13
738 |#
739
740 ;-----------------------------------------------------------------------------
741 ;-----------------------------------------------------------------------------
742 ;;; The functions simple-numerical-integration
743 ;;; Gauss-Legendre-quadrature
744 ;;; are used for numerical integration.
745 ;-----------------------------------------------------------------------------
746 ;-----------------------------------------------------------------------------
747 (defun simple-numerical-integration
748 (f
749 a
750 b
751 &key
752 (accuracy 1.0e-6)
753 (maximum-number-of-iterations 20))
754 "given
755 [1] f (required)
756 ==> a function with a single argument
757 [2] a (required)
758 ==> left-hand limit for numerical integration
759 [3] b (required)
760 ==> right-hand limit for numerical integration
761 i.e., a < b
762 [4] accuracy (keyword; 1.0e-6)
763 ==> desired relative accuracy for computed rood
764 [5] maximum-number-of-iterations (keyword; 20)
765 ==> maximum number of iterations
766 returns
767 [1] the integral of f over the interval [a, b]
768 [2] the number of iterations required
769 ---
770 Note: this function is based on qtrap,
771 Section 4.2, Numerical Recipes, Second Edition"
772 (let ((s-old nil)
773 s-new)
774 (dotimes (i maximum-number-of-iterations
775 (if (not (cerror "returns solution so far"
776 "exceeding maximum number of iterations"))
777 (values s-old maximum-number-of-iterations)))
778 (setf s-new (trapezoidal-rule f a b s-old i))
779 (if (and s-old (< (abs (- s-new s-old))
780 (* accuracy (abs s-old))))
781 (return (values s-new (1+ i))))
782 (setf s-old s-new))))
783
784 #|
785 ;;; first test of trapezoidal-rule
786 (- (exp 1) 1)
787 ;==> 1.718281828459045
788
789 (let (s)
790 (dotimes (i 20)
791 (setf s (print (trapezoidal-rule #'exp 0.0 1.0 s i)))))
792 ;==> 1.718281828459572 (great agreement, but takes a long time!)
793
794 ;;; second test
795 (- (exp pi) (exp 1.0))
796 ;==> 20.42241080432022
797
798 (let (s)
799 (dotimes (i 10)
800 (setf s (print (trapezoidal-rule #'exp 1.0 pi s i)))))
801 ;==> 20.42244057984665
802
803 ;;; now test simple-numerical-integration on these two cases
804 (simple-numerical-integration #'exp 0.0 1.0)
805 ;==> 1.7182823746860927 (good agreement)
806 ; 10
807
808 (simple-numerical-integration #'exp 1.0 pi)
809 ;==> 20.422412665291343
810 ; 12
811 |#
812
813 ;-----------------------------------------------------------------------------
814 (defun Gauss-Legendre-quadrature
815 (lower-limit
816 upper-limit
817 N)
818 "given
819 [1] lower-limit (required)
820 ==> lower limit of integration
821 [2] upper-limit (required)
822 ==> upper limit of integration
823 [3] N (required)
824 ==> number of points to be computed
825 in Gauss-Legendre quadrature
826 returns
827 [1] an N-dimensional vector of abscissas points
828 [2] an N-dimensional vector of weights
829 ---
830 Note: this function is based on gauleg,
831 Section 4.5, Numerical Recipes, Second Edition"
832 (let ((abscissas (make-array N))
833 (weights (make-array N))
834 (eps 3.0d-14)
835 (m (truncate (+ N 1) 2))
836 (N+half (+ N 0.5))
837 (xm (* 0.5 (+ upper-limit lower-limit)))
838 (xl (* 0.5 (- upper-limit lower-limit)))
839 z pp)
840 ;;; Loop over desired roots.
841 (dotimes (i m (values abscissas weights))
842 (setf z (cos (* pi (/ (+ i 0.75) N+half))))
843 :label-1
844 (do ((p1 1.0 1.0)
845 (p2 0.0 0.0)
846 (z1 nil)
847 p3)
848 ;;; since z1 is initially nil,
849 ;;; we cannot exit before the first iteration
850 ((and z1 (<= (abs (- z z1)) eps)))
851 (dotimes (jm1 n)
852 (setf p3 p2
853 p2 p1
854 p1 (/ (- (* (1+ (* 2.0 jm1)) z p2) (* jm1 p3))
855 (1+ jm1))))
856 (setf pp (* n (/ (- (* z p1) p2) (- (* z z) 1.0)))
857 z1 z
858 z (- z1 (/ p1 pp))))
859 (setf (aref abscissas i) (- xm (* xl z))
860 (aref abscissas (- n (1+ i))) (+ xm (* xl z))
861 (aref weights i) (* 2.0 (/ xl (* (- 1.0 (* z z)) pp pp)))
862 (aref weights (- n (1+ i))) (aref weights i)))))
863
864 #|
865 (multiple-value-bind (abscissas weights)
866 (gauss-legendre-quadrature -1.0 1.0 32)
867 (dotimes (i (length abscissas))
868 (format t "~& ~2D ~13,10F ~13,10F"
869 (1+ i)
870 (aref abscissas i)
871 (aref weights i))))
872
873 1 -0.9972638618 0.0070186100
874 2 -0.9856115115 0.0162743947
875 3 -0.9647622556 0.0253920653
876 4 -0.9349060759 0.0342738629
877 5 -0.8963211558 0.0428358980
878 6 -0.8493676137 0.0509980593
879 7 -0.7944837960 0.0586840935
880 8 -0.7321821187 0.0658222228
881 9 -0.6630442669 0.0723457941
882 10 -0.5877157572 0.0781938958
883 11 -0.5068999089 0.0833119242
884 12 -0.4213512761 0.0876520930
885 13 -0.3318686023 0.0911738787
886 14 -0.2392873623 0.0938443991
887 15 -0.1444719616 0.0956387201
888 16 -0.0483076657 0.0965400885
889 17 0.0483076657 0.0965400885
890 18 0.1444719616 0.0956387201
891 19 0.2392873623 0.0938443991
892 20 0.3318686023 0.0911738787
893 21 0.4213512761 0.0876520930
894 22 0.5068999089 0.0833119242
895 23 0.5877157572 0.0781938958
896 24 0.6630442669 0.0723457941
897 25 0.7321821187 0.0658222228
898 26 0.7944837960 0.0586840935
899 27 0.8493676137 0.0509980593
900 28 0.8963211558 0.0428358980
901 29 0.9349060759 0.0342738629
902 30 0.9647622556 0.0253920653
903 31 0.9856115115 0.0162743947
904 32 0.9972638618 0.0070186100
905
906 ;;; excellent agreement with top of page 917 of Abramowitz and Stegun
907 |#
908
909 ;-----------------------------------------------------------------------------
910 ;-----------------------------------------------------------------------------
911 ;;; Everything below here consists of internal symbols in the SAPA package
912 ;;; and should be regarded as "dirty laundry" ...
913 ;-----------------------------------------------------------------------------
914 ;-----------------------------------------------------------------------------
915 (defun zeros-of-polynomial-with-initial-estimates
916 (polynomial initial-estimates
917 &key
918 (degree-of-polynomial
919 (1- (length polynomial)))
920 (scratch (make-array (1+ degree-of-polynomial)))
921 (final-results initial-estimates)
922 (maximum-number-of-iterations 25))
923 (if (< (length initial-estimates) degree-of-polynomial)
924 (error (format nil "too few initial estimates (~D<~D)"
925 (length initial-estimates) degree-of-polynomial)))
926 (cond
927 ((not (eq final-results initial-estimates))
928 (if (< (length final-results) degree-of-polynomial)
929 (error (format nil "not enough space for results (~D<~D)"
930 (length final-results) degree-of-polynomial)))
931 (dotimes (i degree-of-polynomial)
932 (setf (elt final-results i) (elt initial-estimates i)))))
933 (let ((value-of-polynomial (make-array 1))
934 (error-bound-for-value (make-array 1))
935 (number-of-roots-found 0)
936 temp)
937 (dotimes (number-of-iterations (* maximum-number-of-iterations
938 degree-of-polynomial))
939 (dotimes (i degree-of-polynomial)
940 (cond
941 ((or (zerop number-of-iterations) (plusp (abs (elt scratch i))))
942 ;;; evaluate polynomial at point (elt final-results i)
943 ;;; and stuff result into (elt value-of-polynomial 0)
944 (evaluate-polynomial polynomial (elt final-results i)
945 :degree-of-polynomial degree-of-polynomial
946 :complex-values value-of-polynomial
947 :bounds-to-be-computed t
948 :bounds error-bound-for-value)
949 (cond
950 ((> (+ (abs (realpart (elt value-of-polynomial 0)))
951 (abs (imagpart (elt value-of-polynomial 0))))
952 (+ (abs (realpart (elt error-bound-for-value 0)))
953 (abs (imagpart (elt error-bound-for-value 0)))))
954 (setf temp (elt polynomial 0))
955 (dotimes (j degree-of-polynomial)
956 (if (not (= j i))
957 (setf temp (* temp (- (elt final-results i)
958 (elt final-results j))))))
959 (setf (elt scratch i) (/ (elt value-of-polynomial 0) temp)))
960 (t
961 (setf (elt scratch i) 0.0)
962 (incf number-of-roots-found))))))
963 (dotimes (i degree-of-polynomial)
964 (decf (elt final-results i) (elt scratch i)))
965 ;;; if exit is done here, the routine will return t;
966 ;;; otherwise, the routine will quit when the maximum number of
967 ;;; interations have been made and return nil.
968 (if (= number-of-roots-found degree-of-polynomial) (return t)))))
969
970 ;-----------------------------------------------------------------------------
971 (defun initial-estimates-for-zeros-of-polynomial
972 (polynomial
973 &key
974 (degree-of-polynomial (1- (length polynomial)))
975 (initial-estimates (make-array degree-of-polynomial))
976 (scratch (make-array (* 2 (1+ degree-of-polynomial)))))
977 (let ((imax (1+ degree-of-polynomial))
978 (pn (make-array 1))
979 (scratch-offset
980 (make-array (+ degree-of-polynomial 2)
981 :displaced-to scratch
982 :displaced-index-offset degree-of-polynomial))
983 (temp (- (/ (elt polynomial 1)
984 (* (elt polynomial 0) degree-of-polynomial))))
985 x u v)
986 (evaluate-polynomial polynomial temp
987 :degree-of-polynomial degree-of-polynomial
988 :number-of-derivatives degree-of-polynomial
989 :complex-values scratch)
990 (setf (elt scratch degree-of-polynomial)
991 (abs (elt scratch degree-of-polynomial)))
992 (dotimes (i degree-of-polynomial)
993 (setf (elt scratch (+ degree-of-polynomial i 1))
994 (- (abs (elt scratch (- degree-of-polynomial i 1)))))
995 (if (< (realpart (elt scratch (+ degree-of-polynomial i 1)))
996 (realpart (elt scratch imax)))
997 (setf imax (+ degree-of-polynomial i 1))))
998 (setf x (expt (- (/ (realpart (elt scratch imax))
999 (realpart (elt scratch degree-of-polynomial))))
1000 (/ (float (- imax degree-of-polynomial)))))
1001 (loop
1002 (setf x (* 2.0 x))
1003 (evaluate-polynomial scratch-offset x
1004 :degree-of-polynomial degree-of-polynomial
1005 :complex-values pn)
1006 (if (>= (realpart (elt pn 0)) 0.0) (return)))
1007 (setf u (* 0.5 x))
1008 (setf v x)
1009 (loop
1010 (setf x (* 0.5 (+ u v)))
1011 (evaluate-polynomial scratch-offset x
1012 :degree-of-polynomial degree-of-polynomial
1013 :complex-values pn)
1014 (if (plusp (realpart (elt pn 0)))
1015 (setf v x)
1016 (setf u x))
1017 (if (<= (- v u) (* 0.001 (1+ v))) (return)))
1018 (dotimes (i degree-of-polynomial)
1019 (setf u (* (/ pi degree-of-polynomial) (+ 0.5 (* 2.0 i))))
1020 (setf (elt initial-estimates i)
1021 (+ temp (* (max x (* 0.001 (abs temp)))
1022 (complex (cos u) (sin u))))))))
1023
1024 ;-----------------------------------------------------------------------------
1025 (defun complex-of-absolutes (z)
1026 (complex (abs (realpart z)) (abs (imagpart z))))
1027
1028 ;-----------------------------------------------------------------------------
1029 ;;; trapezoidal-rule is based upon trapzd, Section 4.2, Numerical Recipes,
1030 ;;; Second Edition. It integrates f from a to b;
1031 ;;; on the first pass, should set n = 0, in which case s is ignored;
1032 ;;; for n = 1, 2, ..., calculates refined estimate of integral using s
1033 ;;; returned at previous stage; i.e.,
1034 ;;; at each stage n=m, function returns updated value of s,
1035 ;;; which should be used to call the function at stage n=m+1
1036 (defun trapezoidal-rule (f a b s n)
1037 (cond
1038 ((zerop n)
1039 (* 0.5 (- b a) (+ (funcall f a)
1040 (funcall f b))))
1041 (t
1042 (let* ((it (expt 2 (- n 1)))
1043 (del (/ (- b a) it))
1044 (x (+ a (* 0.5 del)))
1045 (sum 0.0))
1046 (dotimes (j it (* 0.5 (+ s (/ (* sum (- b a)) it))))
1047 (incf sum (funcall f x))
1048 (incf x del))))))
1049

  ViewVC Help
Powered by ViewVC 1.1.5