Exercise1-28 <---> Exercise1-30

Exercise 1.29

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 $$ \frac{h}{3}[y_0 + 4y_1 + 2y_2 + 4y_3 + 2y_4+ \ldots + 2y_{n-2}+4y_{n-1}+y_n]$$ where h = (b - a)/n, for some even integer n, and yk=f(a + hk). (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.

Ru: Правило Симпсона - более точный метод численного интегрирования, чем представленный выше. С помощью правила Симпсона интеграл функции f между a и b приближенно вычисляется в виде $$ \frac{h}{3}[y_0+4y_1+2y_2+4y_3+2y_4+ \ldots +2y_{n-2}+4y_{n-1}+y_n]$$ где h = (b - a)/n, для какого-то четного целого числа n, а yk=f(ak + h). (Увеличение n повышает точность приближенного вычисления.) Определите процедуру, которая принимает в качестве аргументов f, a, b и n, и возвращает значение интеграла, вычисленное по правилу Симпсона. С помощью этой процедуры проинтегрируйте cube между 0 и 1 (с n = 100 и n = 1000) и сравните результаты с процедурой integral, приведенной выше.

Note: in the solutions below, "sum" refers to the function defined in exercise 1-30.

Scheme solution: 

(define (inc n) (+ n 1))
(define (even? n) (= (remainder n 2) 0))
(define (simpson f a b n)
  (define h (/ (- b a) n))
  (define (multiplier k) ;; to determine coeff near y[k]
    (cond ((= 0 k) 1)
          ((= n k) 1)
          ((even? k) 2)
          (else 4)))
  (define (term k) ;; sum item
    (* (multiplier k)
       (f (+ a (* k h)))))
  (* (/ h 3.0)
     (sum term 0 inc n)))

Another variation: 

(define (inc1 n) (+ n 1))
(define (inc2 n) (+ n 2))

(define (simpson f a b n)
  (let* ((h (/ (- b a) n))
         (term (lambda(k) (f (+ a (* k h))))))
    (* (/ h 3.0)
       (+ (sum term 0 inc1 n)
          (sum term 1 inc1 (- n 1))
          (* 2.0 (sum term 1 inc2 (- n 1)))))))

Haskell solution: (Note: there is also a Prelude function called "sum"; so make sure to define it from exercise 1.30) 

simpson f a b n = (h/3) * sum term 0 (1+) n
  where h = (b-a) / fromIntegral n
        term k = multiplier k * f (a + fromIntegral k * h) -- sum item
        multiplier 0 = 1             -- to determine coeff near y[k]
        multiplier k | k == n    = 1
                     | even k    = 2
                     | otherwise = 4

Another variation: 

simpson f a b n = (h/3) * (sum term 0 (1+) n
                           + sum term 1  (1+) (n-1)
                           + 2 * sum term 1 (2+) (n-1))
  where
    h = (b-a) / fromIntegral n
    term k = f (a + fromIntegral k * h)

OCaml solution: 

let simpson f a b n =
  let h = (b -. a) /. float n in
  let multiplier = function  (* to determine coeff near y[k] *)
     0                  -> 1.
   | k when k = n       -> 1.
   | k when k mod 2 = 0 -> 2.
   | _                  -> 4.
  in
  let term k = multiplier k *. f (a +. float k *. h) in  (* sum item *)
    (h /. 3.) *. sum term 0 succ n

Another variation: 

let simpson f a b n =
  let h = (b -. a) /. float n in
  let term k = f (a +. float k *. h) in
    (h /. 3.) *. (sum term 0 succ n
                  +. sum term 1 succ (n-1)
                  +. 2. *. sum term 1 ((+) 2) (n-1))

Standard ML solution: 

fun simpson (f, a, b, n) = let
  val h = (b - a) / real n
  fun multiplier 0 = 1.0            (* to determine coeff near y[k] *)
    | multiplier k =
        if      k = n       then 1.0
        else if k mod 2 = 0 then 2.0
        else                     4.0
  fun term k = multiplier k * f (a + real k * h)  (* sum item *)
in
  (h / 3.0) * sum (term, 0, fn x => x + 1, n)
end

Another variation: 

fun simpson (f, a, b, n) = let
  val h = (b - a) / real n
  fun term k = f (a + real k * h)
in
  (h / 3.0) * (sum (term, 0, fn x => x + 1, n)
                + sum (term, 1, fn x => x + 1, n-1)
                + 2.0 * sum (term, 1, fn x => x + 2, n-1))
end

Exercise1-28 <---> Exercise1-30

Comments

(define (simpson-integral f a b n)
  (define h (/ (- b a) n))
  (define (addtwo x) (+ x 2.0))
  (define (fs i) (f (+ a (* h i))))
  (* (/ h 3)
     (+ (fs 0)
        (* (sum fs 1 addtwo (- n 1)) 4)
        (* (sum fs 2 addtwo (- n 2)) 2)
        (fs n))))
Posted by nx at 2009-01-05 06:56:29 X

You have "yk=f(ak + h)", that should be "yk = f(a + kh)"

I'm pretty sure that all of the solutions are right, however.

Posted by 24-196-171-251 at 2009-03-10 10:09:27 X

:) :)) :( ;) :\ |) X-( B) Markup

Exercise1-29 (last edited 2011-07-24 20:17:30 by NickDaly)