2

以下は、 SICP 演習 1.29の私のコードです。この演習では、高次の手順を使用してシンプソンの規則を実装するよう求められますsum。integral元の手順よりも正確になるはずです。しかし、私のコードではそうではない理由がわかりません:

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

私のコードのいくつかの説明:

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

等しい

h/3 * (y_{0}
       + 4 * (y_{1} + y_{3} + ... + y_{n-1})
       + 2 * (y_{2} + y_{4} + ... + y_{n-2})
       + y_{n})

とsumを計算するために使用します。y_{1} + y_{3} + ... + y_{n-1}y_{2} + y_{4} + ... + y_{n-2}

ここに完全なコード：

#lang racket

(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 (integral f a b dx)
  (define (add-dx x) (+ x dx))
  (* (sum f (+ a (/ dx 2.0)) add-dx b)
     dx))

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

いくつかのテスト (正確な値は 0.25 である必要があります):

> (integral cube 0 1 0.01)
0.24998750000000042
> (integral cube 0 1 0.001)
0.249999875000001

> (simpson-integral cube 0 1.0 100)
0.23078806666666699
> (simpson-integral cube 0 1.0 1000)
0.24800798800666748
> (simpson-integral cube 0 1.0 10000)
0.2499999999999509
4

2 に答える 2

1

2項の作成方法に問題があります。偶数項 ( を掛ける) と奇数項 ( を掛ける) を交互に使用する方法4は正しくありません。sum現在の用語の偶数または奇数の性質を追跡するために追加のパラメーターを渡すことでこの問題を解決しました。他の方法もありますが、これは私にとってはうまくいき、精度が向上しました:

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

(define (simpson-integral f a b n)
  (let* ((h (/ (- b a) n))
         (term (lambda (x i)
                 (if (even? i)
                     (* 2.0 (f x))
                     (* 4.0 (f x)))))
         (next (lambda (x) (+ x h))))
    (* (+ (f a)
          (sum term a next b 1)
          (f b))
       (/ h 3.0))))

(simpson-integral cube 0 1 1000)
=> 0.2510004999999994
于 2013-10-31T14:36:20.103 に答える