练习 1.29

我们先将题目给出的条件翻译成表达式,将公式主体的分析放到后面。

变量 h 的定义由条件 \(h = (b - a)/n\) 给出,将它翻译成表达式:

(define h (/ (- b a)
             n))

变量 ab 是自由的,因为等会 h 就要放到辛普森函数主体内,所以,不用担心这些自由变量。

继续将下一个条件 \(y_k = f(a + kh)\) 翻译成函数:

(define (y k)
    (f (+ a (* k h))))

函数 \(f\) ,变量 \(a\)\(h\) 都是自由变量。

现在将注意力放到近似值公式上:

\(\frac{h}{3}[ y_0 + 4y_1 + 2y_2 + 4y_3 + 2y_4 + \dotsi + 2y_{n-2} + 4y_{n-1} + y_n ]\)

在公式的最外边,大括号之外,乘上了 \(\frac{h}{3}\) ,我们同样可以在自己的函数上也用 (* (/ h 3) ...) 来达到同样的效果。

在公式的内部,大括号包围的部分,计算多个 \(y_k\) 之和,根据下标 \(k\) 的不同, \(y_k\) 有三个不同的乘法因子:

  • \(k = 0\)\(k = n\) 时,乘法因子为 \(1\)
  • \(k\) 为奇数时,乘法因子为 \(4\)
  • \(k\) 为偶数时,乘法因子为 \(2\)

根据这三条规则,可以给出函数 factor ,它接受一个参数 \(k\) ,返回相应的 \(y_k\) 的乘法因子:

(define (factor k)
    (cond ((or (= k 0) (= k n))
            1)
          ((odd? k)
            4)
          (else
            2)))

和之前的几个定义一样, factor 内部也有一个自由变量 \(n\)

有了 factor 函数,我们就可以使用表达式 (* (factor k) (y k)) 计算出大括号内的各个加法项。

综合以上这些条件,现在可以写出完整的辛普森函数了:

;;; 29-simpson.scm

(load "p38-sum.scm")

(define (simpson f a b n)
    
    (define h (/ (- b a) n))

    (define (y k)
        (f (+ a (* k h))))

    (define (factor k)
        (cond ((or (= k 0) (= k n))
                1)
              ((odd? k)
                4)
              (else
                2)))
    
    (define (term k)
        (* (factor k)
           (y k)))

    (define (next k)
        (+ k 1))

    (if (not (even? n))
        (error "n can't be odd")
        (* (/ h 3)
           (sum term (exact->inexact 0) next n))))

sum 的定义照抄书本 38 页的代码:

;;; p38-sum.scm

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

辛普森函数检查了输入参数 n 的情况,确保它是一个偶数,否则引发一个错误。

测试

分别将 n 设为 1001000 ,用 simpson 函数求出 cube01 的积分:

1 ]=> (load "29-simpson.scm")

;Loading "29-simpson.scm"...
;  Loading "p38-sum.scm"... done
;... done
;Value: simpson

1 ]=> (define (cube x)
          (* x x x))

;Value: cube

1 ]=> (simpson cube 0 1 100)

;Value: .24999999999999992

1 ]=> (simpson cube 0 1 1000)

;Value: .2500000000000003

可以看到,和书本 39 页 integral 函数计算出的积分相比, simpson 的计算结果精度更高。

讨论

blog comments powered by Disqus