# 附录 C 数值计算技术

## C.1 辛普森积分法

$$\frac{h}{3}[(y0+y_n)+4(y_1+y_3+\dots +y{n-1})+2(y2+y_4+\dots +y{n-2})]$$

(define integrate-simpson
(lambda (f a b n)
;...


    ;...
(unless (even? n) (set! n (+ n 1)))
;...


    ;...
(let* ((h (/ (- b a) n))
(h*2 (* h 2))
(n/2 (/ n 2))
;...


           ;...
(sum-every-other-ordinate-starting-from
(lambda (x0 num-ordinates)
(let loop ((x x0) (i 0) (r 0))
(if (>= i num-ordinates) r
(loop (+ x h*2)
(+ i 1)
(+ r (f x)))))))
;...


           ;...
(y0+yn (+ (f a) (f b)))
(y1+y3+...+y.n-1
(sum-every-other-ordinate-starting-from
(+ a h) n/2))
(y2+y4+...+y.n-2
(sum-every-other-ordinate-starting-from
(+ a h*2) (- n/2 1))))
(* 1/3 h
(+ y0+yn
(* 4.0 y1+y3+...+y.n-1)
(* 2.0 y2+y4+...+y.n-2))))))


$$\phi(x) = \frac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{2}}$$

(define *pi* (* 4 (atan 1)))

(define phi
(lambda (x)
(* (/ 1 (sqrt (* 2 *pi*)))
(exp (- (* 1/2 (* x x)))))))


(integrate-simpson phi 0 1 10)
(integrate-simpson phi 0 2 10)
(integrate-simpson phi 0 3 10)


## C.2 自适应区间长度

(define integrate-adaptive-simpson-first-try
(lambda (f a b e)
(let loop ((n 4)
(iprev (integrate-simpson f a b 2)))
(let ((icurr (integrate-simpson f a b n)))
(if (<= (abs (- icurr iprev)) e)
icurr
(loop (+ n 2)))))))


(define integrate-adaptive-simpson-second-try
(lambda (f a b e)
(let integrate-segment ((a a) (b b) (e e))
(let ((i2 (integrate-simpson f a b 2))
(i4 (integrate-simpson f a b 4)))
(if (<= (abs (- i2 i4)) e)
i4
(let ((c (/ (+ a b) 2))
(e (/ e 2)))
(+ (integrate-segment a c e)
(integrate-segment c b e))))))))


(define integrate-adaptive-simpson
(lambda (f a b e)
(let* ((h (/ (- b a) 4))
(mid.a.b (+ a (* 2 h))))
(let integrate-segment ((x0 a)
(x2 mid.a.b)
(x4 b)
(y0 (f a))
(y2 (f mid.a.b))
(y4 (f b))
(h h)
(e e))
(let* ((x1 (+ x0 h))
(x3 (+ x2 h))
(y1 (f x1))
(y3 (f x3))
(i2 (* 2/3 h (+ y0 y4 (* 4.0 y2))))
(i4 (* 1/3 h (+ y0 y4 (* 4.0 (+ y1 y3))
(* 2.0 y2)))))
(if (<= (abs (- i2 i4)) e)
i4
(let ((h (/ h 2)) (e (/ e 2)))
(+ (integrate-segment
x0 x1 x2 y0 y1 y2 h e)
(integrate-segment
x2 x3 x4 y2 y3 y4 h e)))))))))


integrate-segment现在显式地设置了四个h大小的区间，五个纵坐标y0,y1,y2,y3,y4。积分i4用到了所有的坐标值，i2的区间大小是两倍的h，故只用到了y0,y2,y4。很容易看出i2i4用到的和符合辛普森公式中的和。

(integrate-simpson          exp 0 20 10)
(integrate-simpson          exp 0 20 20)
(integrate-simpson          exp 0 20 40)
(- (exp 20) 1)


## C.3 广义积分（反常积分）

$$\Gamma(n)=\int_{0}^{\infty} x^{n-1}e^{-x}dx$$

a. $\Gamma(1)=1$ b. 对n>0， $\Gamma(n+1)=n\Gamma(n)$

(define gamma-1-to-2
(lambda (n e)
(unless (< 1 n 2)
(error 'gamma-1-to-2 "argument outside (1, 2)"))
;...


    ;...
(let ((gamma-integrand
(let ((n-1 (- n 1)))
(lambda (x)
(* (expt x n-1)
(exp (- x))))))
;...


$$\int0^{\infty}g(x)dx \approx \int_0^{x_c}g(x)dx + \int{x_c}^{\infty}t(x)dx$$

          ;...
(e/100 (/ e 100)))
(let loop ((xc 4) (yc (gamma-integrand 4)))
(let* ((tail-integrand
(lambda (x)
(* yc (exp (- (- x xc))))))
(x1 (* 2 xc))
(y1 (gamma-integrand x1))
(y1-estimated (tail-integrand x1)))
(if (<= (abs (- y1 y1-estimated)) e/100)
gamma-integrand
0 xc e/100)
yc)
(loop x1 y1)))))))


(define gamma
(lambda (n e)
(cond ((< n 1) (/ (gamma (+ n 1) e) n))
((= n 1) 1)
((< 1 n 2) (gamma-1-to-2 n e))
(else (let ((n-1 (- n 1)))
(* n-1 (gamma n-1 e)))))))


(gamma 3/2 .001)
(* 1/2 (sqrt *pi*))


[1]. 想了解为什么这种近似是合理的，请参考任意一种基础的微积分教程。 [2]. $\phi$ 是服从正态或高斯分布的随机变量的概率密度函数。其均值为0而方差为1。定积分 $\int_0^z\phi(x)dx$ 是该随机变量在0到z之间取值的概率。然而你并不需要了解这么多也可以理解这个例子！ [3]. 如果Scheme没有atan过程，我们可以用我们的数值积分过程来得到积分 $\int_0^1(1+x^2)^{-1}dx$ 的值，即 $\pi/4$。 [4]. 把常量——如phi中的 (/ 1 (sqrt (* 2 *pi*)))——提取到被积分函数外面，可以加速integrate-simpson中纵坐标的计算。 [5]. 对大于0的实数n来说$\Gamma(n)$是“减小后阶乘”函数（把正整数n映射到(n-1)!）的一个扩展。