123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221 |
- (declare (usual-integrations))
- (define (bisect-2 f x0 x1 eps)
- (let loop ((x0 x0) (fx0 (f x0)) (x1 x1) (fx1 (f x1)))
- (if (= fx0 0.0)
- x0
- (if (= fx1 0.0)
- x1
- (if (> (* fx1 fx0) 0.0)
- (error "root not bounded")
- (let ((xm (/ (+ x0 x1) 2.0)))
- (if (close-enuf? x0 x1 eps)
- xm
- (let ((fxm (f xm)))
- (if (< (* fx1 fxm) 0.0)
- (loop xm fxm x1 fx1)
- (loop x0 fx0 xm fxm))))))))))
- (define (bisect-fp f x0 x1 eps)
- (let loop ((x0 x0) (fx0 (f x0)) (x1 x1) (fx1 (f x1)))
- (if (= fx0 0.0)
- x0
- (if (= fx1 0.0)
- x1
- (if (> (* fx1 fx0) 0.0)
- (error "root not bounded")
- (let ((xm (/ (- (* fx1 x0) (* fx0 x1)) (- fx1 fx0))))
- (if (close-enuf? x0 x1 eps)
- xm
- (let ((fxm (f xm)))
- (if (< (* fx1 fxm) 0.0)
- (loop xm fxm x1 fx1)
- (loop x0 fx0 xm fxm))))))))))
- (define *bisect-break* 60)
- (define *bisect-wallp* #f)
- (define* (bisect f x0 x1 eps #:optional n-break)
- (let ((n-break (if (default-object? n-break) *bisect-break* n-break)))
- (let loop ((x0 x0) (fx0 (f x0)) (x1 x1) (fx1 (f x1)) (iter 0))
- (if *bisect-wallp* (write-line (list x0 x1)))
- (if (= fx0 0.0)
- x0
- (if (= fx1 0.0)
- x1
- (if (> (* fx1 fx0) 0.0)
- (error "root not bounded")
- (let ((xm (if (< iter n-break)
- (/ (+ x0 x1) 2.)
- (/ (- (* fx1 x0) (* fx0 x1)) (- fx1 fx0)))))
- (if (close-enuf? x0 x1 eps)
- xm
- (let ((fxm (f xm)))
- (if (< (* fx1 fxm) 0.0)
- (loop xm fxm x1 fx1 (fix:+ iter 1))
- (loop x0 fx0 xm fxm (fix:+ iter 1))))))))))))
- (define (find-a-root f x0 x1 dx eps continue failure)
- (define (find x0 x1)
- (if (> (abs (- x0 x1)) dx)
- (let ((f1 (f x1)) (f0 (f x0)))
- (if (< (* f0 f1) 0)
- (continue (bisect f x0 x1 eps))
- (let ((xm (/ (+ x0 x1) 2)))
- (find x0 xm)
- (find xm x1))))
- failure))
- (find x0 x1))
- (define (search-for-roots f x0 x1 eps small)
- (define (find-roots x0 x1)
- (let ((f1 (f x1)) (f0 (f x0)))
- (if (< (abs (- x1 x0)) small)
- (if (< (* f0 f1) 0)
- (list (bisect f x0 x1 eps))
- '())
- (let ((xm (/ (+ x0 x1) 2)))
- (append (find-roots x0 xm)
- (find-roots xm x1))))))
- (find-roots x0 x1))
|