附錄 C 數值計算技術
================
遞歸（包括循環）與Scheme的算術基本過程結合可以實現各種數值計算技術。作為一個例子，我們來實現辛普森法則，這是一個用來計算定積分的數值解的過程。
## C.1 辛普森積分法
函數f(x)在區間[a,b]上的定積分可以看作是f(x)曲線下方從x=a到x=b的區域的面積。也就是說，我們把f(x)的曲線繪制在xy平面上，然后找到由該曲線，x軸，x=a和x=b所圍成區域的面積即是積分值。
根據辛普森法則，我們把積分區間[a,b]劃分為n個相等的區間，n是一個偶數（n越大，近似的效果就越好）。區間邊界在x軸上形成了n+1個點，即：$x_0, x_1, \ldots, x_i, x_{i+1}, \ldots, x_n$ ，其中 $x_0=a, x_n=b$ 。每個小區間的長度是h=(b-a)/n，這樣 $x_i=a+i*h$ ，我們然后計算f(x)在區間端點的縱坐標值，即，其中 $y_i=f(x_i)=f(x+i*h)$。辛普森法則用下列算式模擬f(x)在a到b之間的定積分：
$$\frac{h}{3}[(y_0+y_n)+4(y_1+y_3+\dots +y_{n-1})+2(y_2+y_4+\dots +y_{n-2})]$$
我們定義一個過程`integrate-simpson`，該過程接受四個參數：積分函數`f`，積分限`a`和`b`，以及劃分區間的數目`n`。
```scheme
(define integrate-simpson
(lambda (f a b n)
;...
```
首先我們需要在`integrate-simpson`函數里做的是保證`n`為偶數，如果不是的話，我們直接把`n`加上1.
```scheme
;...
(unless (even? n) (set! n (+ n 1)))
;...
```
接下來我們把區間長度保存在變量`h`中。我們引入兩個變量`h*2`和`n/2`來保存`h`的兩倍和`n`的一半，因為我們會在后面的計算過程中經常用到這兩個變量。
```scheme
;...
(let* ((h (/ (- b a) n))
(h*2 (* h 2))
(n/2 (/ n 2))
;...
```
我們注意到 $y_1 + y_3 + \cdots + y_{n?1}$ 與 $y_2 + y_4 + \cdots + y_{n?2}$ 都要把每個縱坐標進行相加。所以我們定義一個本地過程`sum-every-other-ordinate-starting-from`來捕獲公共的循環過程。通過把這個循環抽象為一個過程，我們可以避免重復寫這個循環。這樣不僅減少了勞動量，而且減少了錯誤發生的可能，因為我們只需要調試一個地方即可。
過程`sum-every-other-ordinate-starting-from`接受兩個參數：起始縱坐標和被相加的縱坐標的數量。
```scheme
;...
(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)))))))
;...
```
我們現在可以計算著三個縱坐標的和，然后把它們拼起來得到最后的結果。注意 $y_1+y_3+ \cdots + y_{n?1}$中有n/2項，在 $y_2 + y_4 + \cdots + y_{n?2}$ 中有(n/2)-1項。
```scheme
;...
(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))))))
```
現在我們來用`integrate-simpson`來求下面函數的定積分：
$$\phi(x) = \frac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{2}}$$
我們首先需要用Scheme的S表達式來定義 $\phi$ 。
```scheme
(define *pi* (* 4 (atan 1)))
(define phi
(lambda (x)
(* (/ 1 (sqrt (* 2 *pi*)))
(exp (- (* 1/2 (* x x)))))))
```
注意我們用 $tan^{-1}1 = \frac{\pi}{4}$ 來定義`*pi*`。
下面的調用分別計算了從0到1,2,3的積分值。都使用了10個區間。
```scheme
(integrate-simpson phi 0 1 10)
(integrate-simpson phi 0 2 10)
(integrate-simpson phi 0 3 10)
```
如果精確到小數點后四位，上面的值應該分別是0.3413 0.4772 0.4987。可以看出我們實現的辛普森積分法確實獲得了相當精確的值！
## C.2 自適應區間長度
每次都指定區間數目感覺不是很方便。對某個積分來說足夠好的`n`可能對另一個積分來說差太多。這種情況下，最好指定一個可以接受的誤差`e`，然后讓程序計算到底需要分多少個區間。完成該任務的典型方法是讓程序通過增加n來得到更好的結果，直到連續兩次結果之間的誤差小于`e`。因此：
```scheme
(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)))))))
```
這里我們連續兩次計算辛普森積分（用我們最初定義的過程`integrate-simpson`），`n`從2,4，。。。注意`n`必須是偶數。當當前`n`的積分值`icurr`與前一次`n`的積分值`iprev`的差小于`e`時，我們返回`icurr`。
這種方法的問題是我們沒有考慮對于某個函數來說可能只有某一段或多段能從增長的區間中獲益。而對于函數的其他段而言，區間增長只會增加計算量，而不會讓整體的結果更好。對于一個增長的適應過程而言，我們可以把積分拆成相鄰的幾段，讓每段的精度獨立的增長。
```scheme
(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))))))))
```
初始段是從`a`到`b`，為了找到一段的積分，我們用兩個最小的區間數目2和4來計算辛普森積分`i2`和`i4`。如果誤差在`e`以內，就返回`i4`。如果沒有的話我們就把區間分成兩份，遞歸計算每段的積分并相加。通常，同一層的不同段以它們自己的節奏來匯聚。注意當我們積分半段時，允許的誤差也要減半，這樣才不會產生精度丟失。
這個過程中仍然存在一些低效之處：積分`i4`重新計算了計算`i2`時用到的三個縱坐標，而且每個半段的積分都重新計算了`i2`和`i4`用到的三個縱坐標。我們可以通過顯式保存`i2`和`i4`用到的和并在命名let中傳更多參數來解決這種低效率。這樣有利于在`integrate-segment`內部和連續調用`integrate-segment`時共享一些信息：
```scheme
(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`。很容易看出`i2`和`i4`用到的和符合辛普森公式中的和。
比較下面對積分 $\int_{0}^{20}e^{x}dx$ 的近似：
```scheme
(integrate-simpson exp 0 20 10)
(integrate-simpson exp 0 20 20)
(integrate-simpson exp 0 20 40)
(integrate-adaptive-simpson exp 0 20 .001)
(- (exp 20) 1)
```
可以分析出最后一個是正確的答案。看看你能不能找到一個最小的`n`（如果設得太小會算得超級慢。。。）讓`(integrate-simpson exp 0 20 n)`返回一個和`integrate-adaptive-simpson`算出的差不多的答案？
## C.3 廣義積分（反常積分）
辛普森積分法不能直接用來計算廣義積分（這種積分的被積函數在某個點的值無窮大或者積分區間的端點無窮大）。然而這個積分法仍然可以用于部分積分，而剩下的部分用其他辦法來獲得近似值。比如，考慮伽瑪函數 $\Gamma$ 。對n>0， $\Gamma(n)$ 被定義為下面的積分（積分上限為無窮）：
$$\Gamma(n)=\int_{0}^{\infty} x^{n-1}e^{-x}dx$$
從上式可以看出兩個結論：
a. $\Gamma(1)=1$
b. 對n>0， $\Gamma(n+1)=n\Gamma(n)$
這就意味著如果我們知道 $\Gamma$ 在區間(1,2)上的值，我們就可以知道任何n>0的 $\Gamma(n)$ 。實際上，如果我們放寬條件n>0，我們可以用結論b來把 $\Gamma(n)$ 擴展到n ≤ 0，而函數在n ≤ 0時會發散。
我們首先實現一個Scheme過程`gamma-1-to-2`，其參數n在區間(1,2)內。`gamma-1-to-2`的第二個參數`e`是精確度。
```scheme
(define gamma-1-to-2
(lambda (n e)
(unless (< 1 n 2)
(error 'gamma-1-to-2 "argument outside (1, 2)"))
;...
```
我們引入一個局部變量`gamma-integrand`來保存 $\Gamma$ 中的被積函數 $g(x)=x^{n-1}e^x$ ：
```scheme
;...
(let ((gamma-integrand
(let ((n-1 (- n 1)))
(lambda (x)
(* (expt x n-1)
(exp (- x))))))
;...
```
我們現在需要讓g(x)從0積分到 $\infty$ ，首先我們沒法定義一個“無窮”的區間表示；因此我們用辛普森公式只積分一個部分，如 $[0,x_c]$ （c的意思是截取(cut)），對于剩下的部位，“尾部”，區間 $[x_c,\infty]$ ，我們用一個“尾部”被積函數t(x)來近似g(x)，這樣更加容易分析和處理。事實上，可以很容易看出對于足夠大的 $x_c$ ，我們可以把g(x)替換為一個遞減的指數函數 $t(x)=y_ce^{-(x-x_c)}$ ，其中 $y_c=g(x_c)$ ，因此：
$$\int_0^{\infty}g(x)dx \approx \int_0^{x_c}g(x)dx + \int_{x_c}^{\infty}t(x)dx$$
前一個積分可以用辛普森公式來解出，后一個就是 $y_c$ 。為了找 $x_c$ ，我們從一個比較小的值（如`4`）開始，然后通過每次擴大一倍來改進積分結果，直到在 $2x_c$ 處的縱坐標（即 $g(2x_c)$ ）在一個特定的由“尾部”的被積函數縱坐標誤差之內。對于辛普森積分和尾部的積分計算，我們都需要誤差為`e/100`，就是比`e`再小兩個數量級，這樣對總體的誤差不會有太大影響：
```scheme
;...
(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)
(+ (integrate-adaptive-simpson
gamma-integrand
0 xc e/100)
yc)
(loop x1 y1)))))))
```
我們現在可以寫一個更通用的過程`gamma`來返回任意`n`對應的 $\Gamma(n)$ ：
```scheme
(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(\frac{3}{2})$ ：
```scheme
(gamma 3/2 .001)
(* 1/2 (sqrt *pi*))
```
第二個值是理論上的正確答案。（這是由于 $\Gamma(\frac{3}{2})=\frac{1}{2}\Gamma(\frac{1}{2})$ ，而可以證明 $\Gamma(\frac{1}{2})$ 的值是 $\pi^{\frac{1}{2}}$ ）你可以通過更改過程`gamma`的第二個參數（誤差）讓結果達到任何你想要的近似程度。
--------------------------------------------
[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)!）的一個擴展。
