<< Previous exercise (1.28) | sicp-solutions | Next exercise (1.30) >>
(define (round-to-next-even x)
(+ x (remainder x 2)))
(define (simpson f a b n)
(define fixed-n (round-to-next-even n))
(define h (/ (- b a) fixed-n))
(define (simpson-term k)
(define y (f (+ a (* k h))))
(if (or (= k 0) (= k fixed-n))
(* 1 y)
(if (even? k)
(* 2 y)
(* 4 y))))
(* (/ h 3) (sum simpson-term 0 inc fixed-n)))
This is a similar solution, complete with testing code. There is no rounding of n to the next even number, because the exercise assumes n to be even.
;; ex 1.29 (define (cube x) (* x x x)) (define (inc n) (+ n 1)) (define (sum term a next b) (if (> a b) 0 (+ (term a) (sum term (next a) next b)))) (define (simpson-integral f a b n) (define h (/ (- b a) n)) (define (yk k) (f (+ a (* h k)))) (define (simpson-term k) (* (cond ((or (= k 0) (= k n)) 1) ((odd? k) 4) (else 2)) (yk k))) (* (/ h 3) (sum simpson-term 0 inc n))) ;; Testing (simpson-integral cube 0 1 100) (simpson-integral cube 0 1 1000)
There is a third way which approaches the solution by breaking the problem into four parts: (f y0), (f yn) and two sums,one over even k and another over odd k.
(define (another-simpson-integral f a b n)
(define h (/ (- b a) n))
(define (add-2h x) (+ x (* 2 h)))
(* (/ h 3.0) (+ (f a)
(* 4.0 (sum f (+ a h) add-2h b))
(* 2.0 (sum f (add-2h a) add-2h (- b h)))
(f b))))
Here's a version that sums over pairs of terms (2 y_k + 4 y_k+1). No conditionals or special cases are needed anywhere, but there's an extra term [f(b) - f(a)] to be added to the final count.
(define (simpson f a b n)
(define h (/ (- b a) n))
(define (y k)
(f (+ a (* k h))))
(define (ypair k)
(+ (* 2 (y k))
(* 4 (y (+ k 1)))))
(define (add-2 k)
(+ k 2))
(* (/ h 3) (+ (sum ypair 0 add-2 (- n 1))
(- (f b) (f a)))))
(define (sim-integral f a b n)
(define h (/ (- a b) n))
(define (y k) (f (+ a (* k h))))
(define (coeff k) (if (is-even? k) 2 4))
(define (part-term k) (* (coeff k) (y k)))
(define part-value (sum part-term 1 inc (- n 1)))
(* (/ h 3) (+ (y 0) (y n) part-value)))
Another solution in GNU Guile:
;;;; Under Creative Commons Attribution-ShareAlike 4.0 ;;;; International. See ;;;; <https://creativecommons.org/licenses/by-sa/4.0/>. (define-module (net ricketyspace sicp one twentynine) #:export (simpson)) (define (sum term a next b) (if (> a b) 0 (+ (term a) (sum term (next a) next b)))) (define (simpson f a b n) (let* ((h (/ (- b a) (* 1.0 n))) (y (lambda (k) (+ a (* k h)))) (ce (lambda (k) ;coefficient (cond ((or (= k 0) (= k n)) 1) ((even? k) 2) (else 4)))) (term (lambda (k) (* (ce k) (f (y k))))) (next (lambda (k) (1+ k)))) (* (/ h 3.0) (sum term 0 next n)))) ;;; Guile REPL ;;; scheme@(guile-user)> ,use (net ricketyspace sicp one twentynine) ;;; scheme@(guile-user)> (define (cube x) (* x x x)) ;;; scheme@(guile-user)> (simpson cube 0 1 100) ;;; $23 = 0.24999999999999992 ;;; scheme@(guile-user)> (simpson cube 0 1 1000) ;;; $24 = 0.2500000000000003
While testing the integral example presented before the ex 1.29 I had some problems for small dx [;Aborting!: maximum recursion depth exceeded]. So I tried a iterative version :
(define (sum-iter ans foo a next b)
(if (> a b)
ans
(sum-iter (+ ans (foo a)) foo (next a) next b)))
Using sum-iter in Simpson Integral definition:
(define (simpson-integral foo a b n)
(define h (/ (- b a) n))
(define (y k) (foo (+ a (* k h))))
(define (intersimp k) (* (cond ((or (= k 0) (= k n)) 1)
((even? k) 2)
(else 4))
(y k)))
(define (inc a) (+ a 1))
(* (sum-iter 0 intersimp 0 inc n ) (/ h 3.0))
)
QUESTION
Simpson’s Rule is a more accurate method of numerical integration than the method illustrated above. Using Simpson’s Rule, the integral of a function f between a and b is approximated as
h/3 (y_0 + 4 y_1 + 2 y_2 + 4 y_3 + 2 y_4 + … + 2 y_{n − 2} + 4 y_{n − 1} + y_n)
where h = (b − a)/n, for some even integer n, and y_k = f(a + kh). (Increasing n increases the accuracy of the approximation.) Define a procedure that takes as arguments f, a, b, and n and returns the value of the integral, computed using Simpson’s Rule. Use your procedure to integrate cube between 0 and 1 (with n = 100 and n = 1000), and compare the results to those of the integral procedure shown above.
ANSWER
The integral of cube between 0 and 1 computed using the composite Simpson’s rule is exact because the error of that method is proportional to the fourth derivative of the function to integrate, so for cube the error is zero. The integral of cube between 0 and 1 computed using the composite midpoint rule is not exact because the error of that method is proportional to the second derivate of the function to integrate, so for cube the error is nonzero.