From: Matthias Felleisen (matthias at ccs.neu.edu) Date: Sat Jan 21 17:24:01 EST 2012 Previous message: [racket] adaptive simpson integration Next message: [racket] FFI: Casting pointers Messages sorted by: [date] [thread] [subject] [author]

```Thanks for the code. Sadly I misread the article on Simpson's error estimate, so the contract is all wrong :-(

#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))))

(provide/contract
;; deliver the indeterminate adaptive Simpson integral for f
(->i ((f (-> real? real?)) (L real?) (R  (L) (and/c real? (>/c L))))
(#:epsilon (ε real?))
(r real?))]
[step (-> real? real?)])

;; -----------------------------------------------------------------------------
;; implementation

(define (adaptive-simpson f L R #:epsilon [ε .000000001])
(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* 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))))

(define (mid L R) (/ (+ L R) 2.))

;; simplistic prototype: many calls to f per step
;; (asr f L R ε whole)
(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* 15))]
[else (define epsilon1 (/ ε 2))
(+ (asr f L M epsilon1 left*) (asr f M R epsilon1 right*))]))

;; -----------------------------------------------------------------------------
;; examples
(define const=5 (lambda (x) 5))
(check-equal? (adaptive-simpson const=5 0 1) 5.0)
(check-equal? (adaptive-simpson const=5 0 10) 50.0)

(check-equal? (adaptive-simpson values 0 1) .5)

(define step (lambda (x) (if (< x 1) (sin x) 1)))