All pastes #2055902 Raw Edit

Mine

public text v1 · immutable
#2055902 ·published 2011-05-09 20:53 UTC
rendered paste body
(define abscissa car)
(define weight cdr)

(define coefficients (list
  (cons (- (sqrt (/ (+ 15.0 (* 2.0 (sqrt 30.0))) 35.0))) (/ (- 90.0 (* 5.0 (sqrt 30.0))) 180.0))
  (cons (- (sqrt (/ (- 15.0 (* 2.0 (sqrt 30.0))) 35.0))) (/ (+ 90.0 (* 5.0 (sqrt 30.0))) 180.0))
  (cons (sqrt (/ (- 15.0 (* 2.0 (sqrt 30.0))) 35.0)) (/ (+ 90.0 (* 5.0 (sqrt 30.0))) 180.0))
  (cons (sqrt (/ (+ 15.0 (* 2.0 (sqrt 30.0))) 35.0)) (/ (- 90.0 (* 5.0 (sqrt 30.0))) 180.0))
))

(define min-iteration-count 10)
(define max-iteration-count 50)
(define absolute-accuracy 0.05)
(define relative-accuracy 0.1)

(define (integrate func imin imax)

    (define (stage n)
      (let* (
        (step (/ (- imax imin) n))
        (halfStep (/ step 2.0)))
        (define (i-iter i midPoint sum)
          (define (j-iter coeffs sum)
            (if (null? coeffs)
                sum
                (let ((coeff-pair (car coeffs)))
                     (j-iter (cdr coeffs)
                             (+ sum
                                (* (weight coeff-pair)
                                   (func (+ midPoint
                                            (* halfStep
                                               (abscissa coeff-pair))))))))))
          (if (zero? i)
              sum
              (i-iter (1- i)
                      (+ midPoint step)
                      (j-iter coefficients sum))))
      
        (* halfStep (i-iter n (+ imin halfStep) 0.0))))

    (define (integrate-iter i n oldt)
      (if (zero? i)
          (exit) ; TODO: error continuation or something
          (let* ((t (stage n))
                 (delta (abs (- t oldt)))
                 (limit (max absolute-accuracy
                             (* relative-accuracy
                                (+ (abs oldt) (abs t))
                                0.5))))
                    (if (and (>= (1+ i) min-iteration-count) (<= delta limit))
                        t
                        (let* ((ratio (min 4
                                           (expt (/ delta limit)
                                                 (/ 0.5 (length coefficients)))))
                               (newn (max (truncate (* ratio n))
                                          (1+ n))))
                              (integrate-iter (1- i) newn t))))))
    
    (integrate-iter max-iteration-count 2 (stage 1)))

(display "f(x) = x*2 from x=0 to 3: ")
(display (integrate (lambda (x) (* x 2)) 0.0 3.0))
(display "\nf(x) = sin(x) from x=0 to 4: ")
(display (integrate sin 0.0 4.0))
(display "\n")