codepad
[
create a new paste
]
login
|
about
Language:
C
C++
D
Haskell
Lua
OCaml
PHP
Perl
Plain Text
Python
Ruby
Scheme
Tcl
; numerical integration (define-syntax fold-of (syntax-rules (range in is) ((_ "z" f b e) (set! b (f b e))) ((_ "z" f b e (v range fst pst stp) c ...) (let* ((x fst) (p pst) (s stp) (le? (if (positive? s) <= >=))) (do ((v x (+ v s))) ((le? p v) b) (fold-of "z" f b e c ...)))) ((_ "z" f b e (v range fst pst) c ...) (let* ((x fst) (p pst) (s (if (< x p) 1 -1))) (fold-of "z" f b e (v range x p s) c ...))) ((_ "z" f b e (v range pst) c ...) (fold-of "z" f b e (v range 0 pst) c ...)) ((_ "z" f b e (x in xs) c ...) (do ((t xs (cdr t))) ((null? t) b) (let ((x (car t))) (fold-of "z" f b e c ...)))) ((_ "z" f b e (x is y) c ...) (let ((x y)) (fold-of "z" f b e c ...))) ((_ "z" f b e p? c ...) (if p? (fold-of "z" f b e c ...))) ((_ f i e c ...) (let ((b i)) (fold-of "z" f b e c ...))))) (define-syntax sum-of (syntax-rules () ((_ arg ...) (fold-of + 0 arg ...)))) (define (int-rect f a b n) (let* ((w (/ (- b a) n)) (mid (lambda (i) (+ a (/ w -2) (* i w)))) (g (lambda (i) (* w (f (mid i)))))) (sum-of (g i) (i range 1 (+ n 1))))) (define (int-trap f a b n) (let* ((w (/ (- b a) n)) (lo (lambda (i) (+ a (* (- i 1) w)))) (hi (lambda (i) (+ a (* i w)))) (g (lambda (i) (* w 1/2 (+ (f (lo i)) (f (hi i))))))) (sum-of (g i) (i range 1 (+ n 1))))) (define (int-simp f a b n) (let* ((w (/ (- b a) n)) (lo (lambda (i) (+ a (* (- i 1) w)))) (mid (lambda (i) (+ a (/ w -2) (* i w)))) (hi (lambda (i) (+ a (* i w)))) (g (lambda (i) (* w 1/6 (+ (f (lo i)) (* 4 (f (mid i))) (f (hi i))))))) (sum-of (g i) (i range 1 (+ n 1))))) (define (integrate method f a b . epsilon) (let ((g5 (method f a b 5)) (g10 (method f a b 10)) (mid (/ (+ a b) 2)) (eps (if (null? epsilon) #e1e-7 (car epsilon)))) (if (< (abs (- g10 g5)) eps) g10 (+ (integrate method f a mid eps) (integrate method f mid b eps))))) (define (cube x) (* x x x)) (display (int-rect cube 0 1 10000.)) (newline) (display (int-trap cube 0 1 10000.)) (newline) (display (int-simp cube 0 1 10000.)) (newline) (define (approx-pi n) (inexact->exact (round (integrate int-simp (lambda (x) (/ (log x))) 2 n 1e-7)))) (display (approx-pi 1e21)) (newline)
Private
[
?
]
Run code
Submit