Решение упражнения 1.29 из SICP

16 October, 2007 (20:14) | Решения упражнений

Итак, мы переходим к упражнениям раздела 1.3, в котором речь идет о процедурах высших порядков. Поддержка процедур высших порядков и вообще работа с процедурами как с данными является очень мощным средством построения абстракций. В упражнениях мы как раз этими абстракциями и займемся.

Реализация формулы Симпсона для вычисления определенного интеграла выглядит следующим образом:

(define (simpson f a b n) 
  (define h (/ (- b a) n)) 
  (define (g k) 
    (define (c k) 
      (cond ((= k 0) 1) 
            ((= k n) 1) 
            ((even? k) 2) 
            (else 4))) 
    (* (c k) (f (+ a (* k h))))) 
  (define (inc k) (+ k 1)) 
  (/ (* (sum g 0 inc n) h) 3))

Заметим, что мы переиспользовали процедуру sum, которая задает схему вычислений, в которую можно уложить и метод Симпсона.

Интегрируя с помощью этой процедуры функцию f(x)=x3 между 0 и 1 с n=100 и n=1000, получим:

> (simpson cube 0 1 100) 
1/4 
> (simpson cube 0 1 1000) 
1/4

То есть видно, что формула Симпсона дает более точный результат, чем формула центральных прямоугольников, использовавшаяся в процедуре integral из основного текста раздела.

Comments

Comment from thror
Date: February 14, 2008, 11:59 pm

Все, у кого не встречал используют формулу Симпсона так… А я использую всегда так. Обратите внимание—корректность определения сохраняется, а чередование коэффициентов получается само собой

(define (simpson f a b n)
  (let ((h (/ (- b a) n)))
    (/ (* h (sum (lambda (x) (+ (f (- x h)) (* 4 (f x)) (f (+ x h))))
                 (+ a h)
                 (lambda (x) (+ x h h))
                 b))
       3)))

–thror

Comment from Sergey Khenkin
Date: February 15, 2008, 9:55 am

А я когда-то встречал и такой вариант тоже.
Я думаю, дело в том, что авторы подводят все формулы под один вид i=1..n ci f(xi) и записывают определение ci, i=1..n отдельно.

Comment from Sergey Khenkin
Date: February 15, 2008, 10:05 am

Возможно, еще одним основанием для использования формулы именно в классическом виде является то, что при этом интегрируемая функция в кажой точке вычисляется ровно один раз (что может быть важно в плане эффективности для сложных функций). В предложенном же варианте функция в среднем вычисляется 1.5 раза в каждой точке.

Мне, кстати, очень нравится такая короткая запись. Симпатично выглядит.

Comment from gorilych
Date: July 15, 2008, 5:19 am

Как вариант:

(define (simpson f a b n)
(define h (/ (- b a) n))
(define dx (* 2 h))
(define (next x) (+ x dx))
(* (/ h 3)
(+ (f a)
(f b)
(* 2 (sum f (+ a dx) next (- b dx)))
(* 4 (sum f (+ a h) next (- b h))))))

Кроме того, формула Симпсона использует кубическое приближение, поэтому в данном случае интегрирования кубической функции не просто “более точный результат” а именно _точный_ результат. В чём можно убедиться, используя n = 2:

1 ]=> (simpson cube 0 1 2)

;Value: 1/4

Comment from psv
Date: March 19, 2009, 12:11 am

полученный алгоритм имеет какие то специэффекты


(define (integral-simps f a b n)
  (define h (/ (- b a) n))
  (define (step2 x) (+ x h h))
  (/ (* h
        (+ (f a)
           (f b)
           (* 2. (sum f (+ a h h) step2 (- b h h)))
           (* 4. (sum f (+ a h)   step2 (- b h)))))
     3.))

guile> (integral-simps cube 0. 1. 60)
0.218833847736626
guile> (integral-simps cube 0. 1. 70)
0.249999999999999
guile> (integral-simps cube 0. 1. 72)
0.249999999999999
guile> (integral-simps cube 0. 1. 74)
0.249999999999999
guile> (integral-simps cube 0. 1. 76)
0.22504211555825





Comment from Irv
Date: January 7, 2016, 10:22 pm

(define (cube x) (* x x x))

(define (sum term a next b)
  (if (> a b)
      0
      (+ (term a) (sum term (next a) next b))))

(define (rectangles-integral f a b dx)
  (define (next x)
    (+ x dx))
  (* dx (sum f (+ a (/ dx 2)) next b)))

(define (simphson-integral f a b n)
  (define h (/ (- b a) n))
  (define (y i)
    (f (+ a (* i h))))
  (define (iter i)
    (if (> n i)
        (+
         (* (if (even? i) 2 4) (y i))
         (iter (+ 1 i)))           
        (y n)))
  (* (/ h 3.0) (+ (y 0) (iter 1))))

(define (good-simphson-integral f a b n)
  (define h (/ (- b a) n))  
  (define (next x) (+ x h))
  (define (term x)
    (define k (/ (- x a) h))
    (* (cond [(= k 0) 1]
             [(= k n) 1]
             [else (if (even? k) 2 4)])
       (f x)))
  (* (/ h 3.0) (sum term a next b)))

(rectangles-integral cube 0 1 .001)
(simphson-integral cube 0 1 1000)
(good-simphson-integral cube 0 1 1000)

Write a comment