123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246 |
- (declare (usual-integrations))
- (define* (zbrent f x1 x2 #:optional tol itmax)
- (if (default-object? tol) (set! tol 0))
- (if (default-object? itmax) (set! itmax 100))
- (let ((fa (f x1)) (fb (f x2)))
- (if (or (and (> fa 0) (> fb 0)) (and (< fa 0) (< fb 0)))
- (error "Root must be bracketed in zbrent" x1 x2))
- (let lp ((iter 0) (a x1) (fa fa)
- (b x2) (fb fb)
- (c x1) (fc fa)
- (d (- x2 x1)) (e (- x2 x1)))
- (cond ((fix:= iter itmax)
- (list #f b iter))
- ((or (and (> fb 0) (> fc 0))
- (and (< fb 0) (< fc 0)))
- (let ((u (- b a)))
- (lp iter a fa b fb a fa u u)))
- ((< (abs fc) (abs fb))
- (lp iter b fb c fc a fa d e))
- (else
- (let ((tol1
- (+ (* *machine-epsilon* (abs b))
- (/ tol 2.0)))
- (xm (/ (- c b) 2.0)))
- (define (next1 p q)
- (let ((p (abs p))
- (q (if (> p 0.0) (- q) q)))
- (let ((min1 (- (* 3.0 xm q) (abs (* tol1 q))))
- (min2 (abs (* e q))))
- (if (< (* 2.0 p) (min min1 min2))
- (next2 d (/ p q))
- (next2 d xm)))))
- (define (next2 e d)
- (let ((nb
- (if (> (abs d) tol1)
- (+ b d)
- (+ b
- (if (> xm 0.0)
- (abs tol1)
- (- (abs tol1)))))))
- (lp (fix:+ iter 1) b fb nb (f nb) c fc d e)))
- (cond ((or (<= (abs xm) tol1) (= fb 0.0))
- (list #t b iter))
- ((and (>= (abs e) tol1) (> (abs fa) (abs fb)))
- (let ((s (/ fb fa)))
- (if (= a c)
- (next1 (* 2.0 xm s) (- 1.0 s))
- (let ((q (/ fa fc)) (r (/ fb fc)))
- (next1 (* s
- (- (* 2.0 xm q (- q r))
- (* (- b a) (- r 1.0))))
- (* (- q 1.0)
- (- r 1.0)
- (- s 1.0)))))))
- (else
- (next2 d xm)))))))))
|