I have coded up the adaptive Simpson integration algorithm below, including a
contract. The algorithm has been mentioned recently on 'users'. Any comments
from the numerically inclined? -- Matthias
#lang racket
(require rackunit)
;; -----------------------------------------------------------------------------
;; interface, with contract
;; Simpson's rule for approximating an integral
(define (simpson f L R)
(* (/ (- R L) 6) (+ (f L) (* 4 (f (mid L R))) (f R))))
(define (mid L R) (/ (+ L R) 2.))
;; adaptive Simpson limits approximate error to 15 ε
(define (simpson-error-estimate f L R ε)
[define M (mid L R)]
(< (abs (+ (simpson f L M) (+ simpson f M R) (- (simpson f L R)))) (* 15 ε)))
(provide/contract
[adaptive-simpson
(->i ((f (-> real? real?))) (#:epsilon (ε real?))
(r (f ε)
(->i ((L real?) (R (L) (and/c real? (>/c L)))) () (r real?)
;; adds six function calls to f, which seems tolerable
#:post (L R) (simpson-error-estimate f L R ε))))])
;; -----------------------------------------------------------------------------
;; implementation
(define (adaptive-simpson f #:epsilon [ε .000000001])
(lambda (L R)
(define f@L (f L))
(define f@R (f R))
(define-values (M f@M whole) (simpson-1call-to-f f L f@L R f@R))
(asr f L f@L R f@R ε whole M f@M)))
;; computationally efficient: 2 function calls per step
(define (asr f L f@L R f@R ε whole M f@M)
(define-values (leftM f@leftM left*) (simpson-1call-to-f f L f@L M f@M))
(define-values (rightM f@rightM right*) (simpson-1call-to-f f M f@M R f@R))
(define delta* (- (+ left* right*) whole))
(cond
[(<= (abs delta*) (* 15 ε)) (+ left* right* (/ delta* whole 15))]
[else (define epsilon1 (/ ε 2))
(+ (asr f L f@L M f@M epsilon1 left* leftM f@leftM)
(asr f M f@M R f@R epsilon1 right* rightM f@rightM))]))
(define (simpson-1call-to-f f L f@L R f@R)
(define M (mid L R))
(define f@M (f M))
(values M f@M (* (/ (abs (- R L)) 6) (+ f@L (* 4 f@M) f@R))))
;; simplistic prototype: many calls to f per step
(define (asr.v0 f L R ε whole)
(define M (mid L R))
(define left* (simpson f L M))
(define right* (simpson f M R))
(define delta* (- (+ left* right*) whole))
(cond
[(<= (abs delta*) (* 15 ε)) (+ left* right* (/ delta* whole 15))]
[else (define epsilon1 (/ ε 2))
(+ (asr.v0 f L M epsilon1 left*) (asr.v0 f M R epsilon1 right*))]))
;; -----------------------------------------------------------------------------
;; examples
(define int-const=5 (adaptive-simpson (lambda (x) 5)))
(check-equal? (int-const=5 0 1) 5.0)
(check-equal? (int-const=5 0 10) 50.0)
(define int-identity (adaptive-simpson values))
(check-equal? (int-identity 0 1) .5)
(define step (lambda (x) (if (< x 1) (sin x) 1)))
(define int-step (adaptive-simpson step))
(int-step 0 2)
____________________
Racket Users list:
http://lists.racket-lang.org/users
