Exercise 1.29 -- 辛普森积分法

来看看这道题：

 

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

 

 

写道

where h = (b - a)/n, for some even integer n, and yk = 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.

 

大概意思就是： 辛普森积分法是一个很牛逼的东西，函数f在a和b之间的值，可以用下面的近似公式来计算：

 

其中，n是一个偶数，h=(b-a)/n, y_k=f(a+kh);n的值越大，积分的结果越精确。需要你写一函数，计算立方运算从0到1的积分，然后与书中的积分函数integral做比较。

 

integral的实现

先看一看wiki上对积分的解释：

就是将一条曲线（当然，可积的）的覆盖面积，变成一个个小矩形的面积的和，而矩形的高，就是矩形中央的x对应的f(x) -- 这个其实无所谓选哪个，只要矩形足够小就可以。 所以，有这么一个公式：

 

 为了能够将“一个个小矩形面积的和累加”，我们需要先写一个sigma来做累加操作：

(defun sum (term a next b)
  (if (> a b)
      0
    (+ (funcall term a)
       (sum term (funcall next a) next b))))

 函数sum的参数a、b相当于上面公式中的a和b，term则相当于上面的f，next则是一个取“步长”的方法。

 

我们要计算立方的积分，那么term是最容易想到的：

(defun cube (x) (* x x x))

 

下面是next，next可以看作是从a开始，每次加上一个很小的数dx：

(defun add-dx (x) (+ x dx))

 所以有了曲线面积：

 

(defun integral (f a b dx)
  (defun add-dx (x) (+ x dx))
  (* (sum f
          (+ a (/ dx 2.0))
          #'add-dx
          b)
     dx))

 我们以0.01为步长（将曲线分为100个小矩形），来算一下它的面积：

(integral #'cube 0 1 0.01)

 

得到的结果是这样的：

写道

0.24998708

 

Simpson's rule

练习题1.29给出了另外一种计算积分的方法，而不是使用前面的矩形面积叠加。简单的说，是这样的：

 

 

在练习题里面，变换成了这样：

 

如果你把“y_0+4y_1+2y_2..."变换一下：

 

写道

y_0 + 4y_1 + y_2
y_2 + 4y_3 + y_4
y_4 + 4y_5 + y_6
...
y_n-2 + 4y_n-1 + y_n

 就和辛普森法则是一样的了。

 

由于题目中规定n一定是偶数，参照前面介绍的方法，我们只需要定义好“a", "b", “term”和next，计算积分了。

 

不论如何，a自然就是我们的起始值了。

term按照上面的定义，

 

term(a) = f(a) + 4 * f(a+h) + f(a+2h)

 

变成方法，就是：

    (defun simp-term (a)
      (+ (funcall f a)
         (* 4 (funcall f (+ a h)))
         (funcall f (+ a (* 2 h)))))

 

 

再看一下上面列出的公式，y的下标从0开始，每次增加2，然后到n-2结束，所以，对于h=(b-a)/n,

b的值应该是 b-2h

next的计算方法应该是 x + 2h

 

把前面的步骤放到一起执行以下：

 

(defun simp-integral (f a b n)
  (let ((h (/ (- b a) n)))
    (defun add-dx (x) (+ x (* 2 h)))
    (defun simp-term (a)
      (+ (funcall f a)
         (* 4 (funcall f (+ a h)))
         (funcall f (+ a (* 2 h)))))
    (* (/ h 3.0)
       (sum #'simp-term a #'add-dx (- b (* 2 h))))
    ))

(simp-integral #'cube 0 1 100)

 

看看结果：

0.24999999

 

与前面的积分计算相比，精度提高了一些。