[racket] adaptive simpson integration

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 at L (f L))
    (define f at R (f R))
    (define-values (M f at M whole) (simpson-1call-to-f f L f at L R f at R))
    (asr f L f at L R f at R ε whole M f at M)))

;; computationally efficient: 2 function calls per step 
(define (asr f L f at L R f at R ε whole M f at M)
  (define-values (leftM  f at leftM  left*)  (simpson-1call-to-f f L f at L M f at M))
  (define-values (rightM f at rightM right*) (simpson-1call-to-f f M f at M R f at R))
  (define delta* (- (+ left* right*) whole))
  (cond
    [(<= (abs delta*) (* 15 ε)) (+ left* right* (/ delta* whole 15))]
    [else (define epsilon1 (/ ε 2))
          (+ (asr f L f at L M f at M epsilon1 left*  leftM  f at leftM) 
             (asr f M f at M R f at R epsilon1 right* rightM f at rightM))]))

(define (simpson-1call-to-f f L f at L R f at R)
  (define M (mid L R))
  (define f at M (f M))
  (values M f at M (* (/ (abs (- R L)) 6) (+ f at L (* 4 f at M) f at 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)