Special Topics II: Algorithms and Numerical Recipes of Integration
Numerical integration has been an active field in mathematics even before the introduction of computers. This comes from the simple fact that not all integrations can be carried out analytically, and numerical methods, however tedious without the help of a computer, become the only way to solve the problem. Furthermore, the results of many common integrals, such as those for the error function erf(x) and gamma function ()x Γ, do not exist as analytical functions and their values must be found numerically. This is quite different from differentiation. For most of the functions encountered in physics and engineering, differentiations can usually be carried out by analytical means. However, there are exceptions. On many occasions, we encounter functions that are given numerically. In this case, the functional form is not available and the derivatives can only be obtained using numerical methods. Furthermore, numerical differentiation is a good way to be introduced to the wide world of finite difference methods for solving a large variety of problems, including differential equations. 2-1 Numerical Integration
One of the more common forms of integration may be represented by
[,]()b
a b a
I f x dx =? （2-1）
This is a standard one-dimensional definite integral with both upper and lower limits of the integration specified. For simplicity, we can assume the integrated f(x) is greater than or equal to zero everwhere in the interval x=[a, b]. Under such conditions, the integral [,]a b I may be interpreted as the area bound above by f(x), below by the x axis, on the left by x=a, and on the right by x=b, as illustrated schematically in Fig. 2-1 by the dotted area. This form of integral is known in mathematics as the (definite) Riemann integral and is the only form with which we shall be dealing here. Certain types of indefinite integrals can also be evaluated using computers. The subject belongs to symbolic manipulation, or computer algebra, and we shall not be concerned with it here.
数值积分的使用可以用以速度v 沿x 方向运动的粒子所经过的距离为例进行说明。若v 为常数，则在时刻t=a 与时刻t=b 之间经过的距离可以简单表为
()d v b a =?-
另一方面，如果速度作为时间的函数而变化，则上式由下面的积分代替
()b
a
d v t dt =? （2-2）
若v(t)可以表为某一解析函数，则积分通常不需借助于数值方法即可解决。然而在许多情况下，v(t)只能以数据表格或其他不能以解析方法积分的形式给出。
作为另一个例子，考虑一辆在公路上行驶的汽车。如果我们只能用一分钟抽样一次的方法得到其速度，那么它在比如说10分钟的时间里经过了多少路程呢？假设抽样的结果就是表2-1的第二列所列出的数值。在这里，经过的路程就是速度对时间的积分这一概念仍然是正确的，只不过因为v(t)不是已知函数，（2-2）式中的积分不能以解析方式写出。然而，通过使用下面的方法，我们可以得到路程的一个非常好的“估计”。由于我们不知道任意时刻的准确速度，我们就以某一分钟的抽样结果作为那一分钟的平均速度t v 。于是，在那一分钟内经过的路程成为60t v ?。表中的第三列给出了这一数据。通过将列中的十个结果相加，我们就得到了十分钟内经过的路程。这基本上就是数值积分的精神。我们无需了解积分变量取到所有值时的被积函数，取而代之，我们将积分区间分解为若干个子区间。在每个子区间内，我们用常数或某种极其简单的函数如线性函数来近似代替被积函数。然后，将每一个部分的贡献相加，就得到了整个区间上的值，正如我们在汽车行驶路程的例子中所做的那样。
Table 2-1: Distance traveled in ten minutes.
用数值方法对积分求值常被称为数值求积法，其理由可在2-4节明显看出。虽然绝大多数方法遵循上一段所述的基本精神，对精确度与效率的要求导致了各种改进。
内容提要：数值积分的精神：将积分区间分解为若干个子区间，在每个子区间内，用常数或某种极其简单的函数如线性函数来近似代替被积函数。然后，将每一个部分的贡献相加，以得到整个区间上的值。 2-2 矩形与梯形法则
作为对（2-1）式数值积分的基本思想的一个简单应用，我们将x=a 与x=b 之间的区间分割为一些更小的子区间，其中第i 个子区间始于1i x x -=，终于i x x =。为简单起见，我们由所有子区间具有相同大小开始
1i i h x x -=-
总共有N 个这样的子区间，即
Nh b a =-
我们希望计算的曲线下面积现在被分割为了N 个窄带，每一个的宽度都为h 。若h 足够小，我们可将每一窄条简化为某种简单的形状。在矩形法则中，我们采用的是矩形，矩形的高则取为f(x)在子区间内的某个合理平均值。在梯形法则中，我们采用在1i x x -=处高为
1()i f x -，在i x x =处高为()i f x 的梯形。当取h 趋近于0的极限时，这两种法则是等效的，
都给出积分的准确值。
矩形法则 为了应用矩形法则，我们需要每个子区间内f(x)的平均值。该值由下式给出
11()i
i x i x f f x dx h
-=
? 若对所有N 个子区间都得到i f 的准确值，我们就得到结果
[,]1
()N
b
a b i a
i I f x dx h f ===∑?
然而，这在应用上是不现实的。矩形法则的目的在于通过对i f 的近似，用最小的计算量高效地求出积分。
对于一个缓慢变化的函数，在小区间内可以用函数在区间中点上的值很好的近似代替其平均值
1/2()i i f f x -≈
其中1/211
()2
i i i x x x --≡
+。为简化符号，在不致引起混淆的情况下令 ()i i f f x ≡ （2-3）
即i f 为f(x)在i x x =处的值。从而f(x)在区间[a,b]上的积分为
[,]1/21
N
a b i i I h f -==∑ （2-4）
注意由于我们已经使0x a =，f(x)在第一个子区间[a,a+h]上的近似应表示为
1/21/2()f f x ≡。
（2-4）式给出的数值积分方法就是矩形法则。容易看出，使用这一方法的误差随着子区间数目N 的增大而减小。考虑子区间[1i x -,i x ]的贡献
11
[,]1/2()i
i i i x x x i x I f x dx hf ---=≈?
（2-5）
我们可以通过在中点1/2i x -的邻域内对函数f(x)作泰勒级数展开来验证这一近似的精确性。
2
1/21/21/21/21/2
(3)
31/21/211()()()1!2!
1()3!
i i i i i i i f x f f x x f x x f x x -------'''=+
-+-+
-+ （2-6）
其中我们用1/2i f -'，1/2i f -''和(3)
1/2i f -分别表示f(x)在1/2i x x -=处的一阶，
二阶和三阶导数。 利用（2-6）式，积分在子区间内的值可以表示为
1
1
1111/21/21/22
1/21/2(3)31/21/21
()()1!
1()2!
1()3!
i
i
i i i i i i i i x x x i i i x x x x i i x x i i x f x dx f dx f x x dx f x x dx f x x dx ------------'=+
-''+
-+-+?
????
（2-7）
其中第一项是（2-5）式中用到的近似，第二项则由于其中的积分等于零而消除。从而，
矩形法则在宽度为h 的单个子区间内的最高阶误差由第三项给出
11
1[,]1/2
321/21/21/2
()1()2!24
i
i i i i i x x x i x x i i i x I f x dx hf h f x x dx f -------?≡-''''≈-=?
?
在整个[a,b]区间上的总误差则通过将所有N 个子区间的贡献相加得到
32[,][,]
2
1()()()()2424b
a b a b a
b a b a I f x dx I h f f N ξξ--''''?=-≈=? 其中我们利用了Nh=(b-a)，并且取()f ξ''为f(x)在[a,b]上的二阶导数的均值。 梯形法则 作为矩形法则的一种替代方法，我们用一个梯形来近似代替每个子区间的面积，梯形的四个顶点分别位于1(,0)i x -，11(,)i i x f --，(,)i i x f 与(,0)i x 。梯形的面积为
1[,]1()2
i i x x i i h
I f f --=
+ （2-8） 图2-2示出了两种法则的差别。
对于整个[a,b]区间，积分值由所有窄带的和给出
[,]
101211111
()()222
N a b i i N N i I h f f h f f f f f --==+=+++++∑ （2-9） 这是我们从一开始就可以预见到的非常合理的结果。与其他值相比，0f 与N f 的贡献只具有一半的重要性，这是由于它们位于区间的两个端点。由Euler-Maclaurin 积分求和公式，此结果还跟着有（例如参见Abramowitz and Stegun[1]）
01212
2(21)(21)220011(){}22
()()2!(2)!
b
N N a
k k k k N N f x dx h f f f f f B B h f f h f f k ---=+++++''------?
（2-10）
系数2k B 为伯努利数，其中
11
2
B =-
而其他所有奇数阶项为零。前几个偶数阶项的值如下
01B = 216B =
4130B =- 6142
B = 更一般的，它们由产生函数
01!
k
k t k t t B e k ∞
==-∑ （2-11） 给出。从k B 直到60B 的明确值可以在Abramowitz and Stegun[1]的第810页找到。同时应注意，（2-11）式说明梯形法则的误差由h 的偶数次方给出。
梯形法则的优点在于f(x)的值只在格点（网点）上，即01,,,N x x x x = 处求出。由问题2-1可以看出，此方法中由有限步长引起的误差的预期值是矩形法则中的2倍。而另一方面，梯形法则的简洁性则有助于进行进一步的改进。在3-5节中我们将看到，此方法可以与外插法一起使用而构成一种非常高效的数值积分算法。框图2-1列出了梯形法则算法的主要步骤。
该算法还包含了一种方法，即由一个小的子区间数开始，通过循环迭代逐步增加数目，而无需在格点上重复计算被积函数的值。稍后我们将看到，这在利用此方法作为更高级技巧——例如外插法——的核心，以近似计算积分在h 趋近于零时的值时十分有用。
内容提要：将曲线下面积分割为N 个等宽窄条，将每一窄条简化为某种简单的形状。在矩形法则中采用的是矩形，在梯形法则中采用的则是梯形。矩形法则在小区间内用函数在区间中点上的值近似代替其平均值。梯形法则中梯形的四个顶点分别位于1(,0)i x -，
11(,)i i x f --，(,)i i x f 与(,0)i x 。矩形法则的误差比梯形法则小，但梯形法则的简洁性有助于
进行进一步的改进。
2-3 Simpson 法则
在上一节中我们看到，矩形法则与梯形法则所伴随的误差都正比于1/N 2。为了提高精确度，同时避免大量增加计算量，我们可以应用Simpson 法则。在（2-7）式所给出的积分的泰勒展开式中，可以看到，若积分的上下限相对于展开的中心点对称，则含有f(x)的奇数阶导数的项都将等于零。利用这一性质，我们可以在相邻的两个子区间内对面积作泰勒级数展开。若展开是在两子区间合并后的中点i x 的邻域内，则可得
1
111
1
1111111
[,]2
(3)335(4)()111'()''()()...
1!2!3!12
20''0()
2!3
i i i i i i i i i i i i x x x x x x x x i i i i i i i x x x x i i i I f x dx
f dx f x x dx f x x dx f x x dx f h f h O h f +-+-++++----==+-+-+-+=?++?++?
???? （2-12）
可见，如果将含有f(x)二阶导数的项包含进数值积分中，可使积分的精确度提高两个级次，达到4
1/N 。
此目的可通过以下方法达到。利用2-7节中将介绍的中心差分，可将f(x)的二阶导数近似表为
211221
''()|(2)i i x x i i i d f f x f f f dx h
=--≡≈-+ （2-13）
将此结果代入（2-12）式，则子区间[xi-1,xi+1]内的积分近似为
115(4)[,]115(4)111
2(2)()
3
141
()()
333
i i x x i i i i i i i i i I hf h f f f O h f h f f f O h f -+----=+-++=+++ （2-14）
对于整个[a,b]区间，结果为
[,]0123211
(42424)3
a b N N N I h f f f f f f f --=+++++++ （2-15）
我们已经假设区间的数目N 为偶数。如图2-3所示，奇数格点的贡献是偶数格点贡献的2倍。这种权重上的差异来自我们为修正基本方法所得的一级结果所引入的f(x)二阶导数的贡献。两端点0f 和N f 的权重则仅为偶数点的一半。
除提高精确度以外，（2-15）式给出的Simpson 法则的另一优点在于其自然引出了一种算法，即通过迭代使积分达到所需要的精确度。原则上，通过估计式（2-12）中省略的高阶项的贡献大小，我们就可以得出数值计算中的误差。但是，由于这需要知道一些高阶导数（如（2-12）式中的(4)
()f
ξ）的平均值，这样做并不容易。比较相继的两次迭代结果是较为简
便的一种方法，尽管较为不确定。如果计算值I[a,b]之差小于预设的允许值ε，则很可能已经达到了所需的精确度。
对于Simpson 法则，该算法可通过如下方法实施。假设我们在x=a 和x=b 之间对函数f(x)进行积分，且要求达到精确度ε。由（2-15）式可见，对积分的贡献可分为三个部分，即端点区，奇数点与偶数点，
01353124642
()()()d N
o N N e N N S N f f S N f f f f f S N f f f f f ----≡+≡+++++≡+++++ 其中()k f f a kh ≡+，()/h b a N =-。式（2-15）可用d S ，o S 及e S 以如下形式写出
[,]1
(){()2()4()}3
a b d e o I N h S N S N S N =++ （2-16）
由此可见，对于一个新的N ，我们可利用上一次循环中已有的d S ，o S 及e S 的值。
为便于进行迭代计算，应从N=1N 开始，其中1N 应为某个较小而合理的区间数目，例如取1N =6。这是迭代计算的第一步。在下一轮迭代中，将区间的数目翻倍，而令步长h 减半，即令N=2N =21N 。以这种方式改变子区间数目的好处在于，对于一个新的N ，来自两个端点的贡献（存储为d S ）是不变的，而偶数点与上一次迭代计算过的偶数点完全相同。也就是说，在新的迭代中的e S 的值就是上一次循环的o S 与e S 之和。所需计算的仅仅是函数在（新的）奇数点上的值。从而我们得到结果
211()()()e e o S N S N S N =+ 2/221
()((21))N o i S N f a i h ==
+-∑
其中2()/h b a N =-。现在我们就可以由（2-16）式得到I[a,b]在N=N 2时的新值了。如
果该值与前一循环所得结果之差小于ε，则可认为计算结果已收敛；否则即再次令子区间数目加倍，使N=N 3=2N 2。然后依此类推，直到计算结果收敛或者达到最大迭代次数。框图2-2给出了算法的主要步骤。
举一个简单的例子就可以说明Simpson 法则的威力。考虑积分
/2
sin I xdx π=?
众所周知，其结果为cos(0)=1。为了用数值方法求出其值，我们由N=6开始迭代。除了分别位于0和π/2的两个端点，我们还分别在i x =/12π，/4π及5/12π处有三个奇数点，在i x =/6π，/3π处有两个偶数点。在第一次循环中，求出了被积函数在这七个点上的值。表2-2给出了积分间隔h 的大小及所得的积分值。
Table 2-2: Example of the accuracy of Simpson’s rule for different N.
在这个简单的例子中，我们实际上可以用（2-12）式来估计对任一给定的N 的预期误差。在任一给定的子区间中，在Simpson 法则已包含的二阶导数项之后的第一个非零项是
(4)512()()4!5
i i N f x h ?=
? 对于整个积分，我们必须将此结果乘以N 并求出f(x)在[a,b]区间内的四阶导数的平均值
(4)()f ξ
(4)412
()()()4!5
N f h b a ξ?=
?- 由于正弦函数的偶数阶导数仍然是正弦函数，我们得到|(4)
()1f ξ≤。在所讨论的例子
中，(b-a)=/2π。从而在误差估计中我们可以认为(4)
()()b a f ξ-与1具有相同量级。由此
得到
8210-?
其对应于不同的h 的具体数值在表2-2的最后一列中给出。
在第二列中给出了数值积分的结果。将这些结果与积分的准确值进行比较即可看出，我们的误差估计具有正确的数量级。如果将子区间的数目翻倍并使N=12，我们期望得到一百万分之一的精确度。表中给出的数据同时显示，在单精度计算中，在N=24之后继续进行迭代已经没有意义，因为此时算法所带来的误差已经可比拟于计算机的截断（舍入）误差了。
Simpson 法则还有若干变种。例如，可改变不同格点的贡献的相对权重，以达到比（2-15）式更高阶的精确度。同时，利用外插方法还可提高算法的效率，对此本书将在3-5节讨论Romberg 积分方法时予以详细介绍。
内容提要：应用Simpson 法则可以提高数值积分精确度，同时避免大量增加计算量。将f(x)的二阶导数近似表出，代入子区间内积分的泰勒展开式，对于整个[a,b]区间求和。引入f(x)二阶导数带来了格点权重上的差异，并且引出了一种迭代算法，可达到预设的精确度。 2-4 高斯求积法
在梯形法则与Simpson 法则中所用的近似可以看作是把每个子区间中的被积函数分别替换为二次与三次多项式。用多项式近似代替被积函数的想法同样可以应用于整个[a,b]区间。这就是高斯求积法——数值积分的高斯方法——的基本思想。
这一方法有若干种不同的形式，其中典型的一种是高斯-勒让德积分。一般地，在给定区间[a,b]上的函数f(x)可以以一组完备的多项式()l P x 为基础展开
()()n
l l l f x P x α==∑ （2-17）
其中展开系数l α的定义与所用的多项式()l P x 有关，稍后将给出例子。对于每一个多项式，下标l 表示其次数，即变量x 的最高幂次，在（2-17）式中由求和上限n 给出。
一般来说，采用一组正交的多项式来展开函数更为方便。原则上可使用任意正交的多项式，在这里我们采用将在4-2节中讨论的勒让德多项式，其正交性表示为
1
,1
2
()()21
k l k l P x P x dx l δ+-=
-?
（2-18） 由于勒让德多项式仅在区间[-l ,+l ]上有定义，我们需要将任意一个积分的上下限由[a,b]（其中a,b 有限）变为[-l ,+l ]，即作如下代换
22
b a b a
x x -+→
+
（2-19）
例如，（2-1）式给出的积分应改写为
1[,]1()()222
b
a b a
b a b a b a
I f x dx f x dx +---+=→
+?? 在以下的讨论中，我们将仅限于讨论
1
[1,1]1
()I f x dx +-+-=? （2-20）
并以此积分作为原型。
考虑被积函数可以近似表为一个（2n-1）次多项式函数的情况
21()()n f x p x -≈ （2-21）
换句话说，如果我们以（2-17）式给出的形式将f(x)按勒让德多项式展开，我们只需用到最高为l =(2n-1)次的()l P x 。如果近似足够好的话，将f(x)的积分用2n 个参数表达出来是可能的，而高斯求积法为我们提供了一种达到最精确值的方法。也就是说，我们希望将积分写为这样的形式
[1,1]1
()n
i i i I f x ω-+==∑ （2-22）
其中待定的横坐标i x 与权重因子i ω应能使我们得到所用计算量下可能达到的最高精确度。作为例子，图2-4示出了使用四阶高斯-勒让德求积法计算的积分
1
1
x e dx +-?
。
横坐标与权重因子 为了确定（2-22）式中的横坐标与权重因子，我们暂时假设被积函数是(2n+1)次的多项式，也就是说（2-21）式是一个准确的表达式而不是一个近似。利用多项式分解可将21()n p x -分解为两项之和
2111()()()()n n n n p x p x P x q x ---=+ （2-23）
其中P n (x)是n 阶勒让德多项式。等式右边的另外两个量1()n p x -和1()n q x -是(n-1)次多项式。
在从-1到+1进行积分时，（2-23）式右边第一项消失了。这是多项式的正交性以如下方式造成的。由于1()n p x -被构造为(n-1)次多项式，因而可以用直到(n-1)阶的勒让德多项式将其展开为
1
10
()()n n k k k p x a P x --==∑ （2-24）
由于此表达式不包含()n P x ，我们得到
1
1
1
11
1
()()()()0n n n k k n k p x P x dx a P x P x dx -++---===∑?
?
其中我们交换了积分和求和的顺序，并且利用了（2-18）式给出的()l P x 的正交性。由此得到等式
1
1
2111
1
()()n n p x dx q x dx ++----=?
?
此外，作为一个n 次多项式，函数()n P x 在区间[-1,+1]内有n 个零点。设这些零点分别位于x=1x ,2x ,…,n x 。在这些点上有关系
211()()n i n i p x q x --= （2-25）
这可以由（2-23）式看出。
类似于（2-24）式，我们同样可以将多项式1()n q x -按勒让德多项式展开
1
10
()()n n k k k q x P x ω--==∑
其中k ω是多项式()k P x 的权重。结合（2-25）式，我们在()n P x 的根上得到关系
1
210
()()n n i k k i k p x P x ω--==∑
或用矩阵形式表为
1
n i ik k k p P ω-==∑ （2-26）
其中使用了符号
21()i n i p p x -≡ 以及 ()ik k i p p x ≡
注意()k i P x 代表了在()n P x 零点处的直到k=(n-1)阶的勒让德多项式的值。i x 的值独立于多项式21()n p x -。
内容提要：在整个[a,b]区间用多项式近似代替被积函数，此即高斯求积法的主要思想。采用正交的勒让德多项式来展开函数，希望将积分写为如[1,1]1
()n
i
i
i I f x ω-+==
∑的形式。通过
求待定的横坐标i x 与权重因子i ω得到所用计算量下可能达到的最高精确度。 表2-3：高斯-勒让德积分的横坐标与权重因素
式（2-26）的关系可以转化为以i p 为变量求k ω值的形式：
n
1k i ik i 1
p {P }ω-==∑ (2-27)
此时P -1 是矩阵{P ik }的逆矩阵， 即：
1
ik kj ij k
{P
}P δ-=∑ 。这个展开式有助于我们接下来
求积分的值：如果认为式（2-21）是精确值，即有：
11
n 1
1
[1,1]2n q k k -1
k 0
1
1
I f (x)dx p (x)dx P (x)dx ω++-+-+-=--===∑??? （2-28）
取k=0， 得P 0（x ）=1，并且根据勒让德多项式的正交性，我们得到
1
k k0k01
2
P (x)dx 22k 1
δδ+-=
=+?
n 项求和中唯一的非零积分。于是，式（2-28）的最终结果为：
n
1[1,1]0k ko k 1
I 22p {P }ω--+===∑ （2-29）
最后的等式从式（2-27）得来。
为了数值积分计算，式（2-29）可以改写为我们更加熟悉的形式。已知
i 2n 1i i p p (x )f (x )-≡≈
表 2-4：1
x 1
e dx +-?
的高斯勒让德积分
f(x)的积分以式（2-22）的形式表达即为：
n
[1,1]i i i 1
I f (x )ω-+==∑
其中i ω是i f (x )在Pn(x)在零点处的权重，并给出
1i i02{P }ω=
我们并不关心如何逆变换多项式矩阵k i P {P (x )}=。5-2节将讨论其通用方法。现在我们把得到的结果写成：
i 2`2
i n i 2
(1x )[P (x )]
ω=
- （2-30） 其中`n i n i d
P (x )P (x)|x x dx
=
= P(x i )的值可以用4-2节所述方式计算得出。我们还可以在《Abramowitz and Stegun 》中查阅到几张扩展表格，内含在不同阶的高斯积分方法时的i x 和i ω值。表2-3为，阶数是n=4，5 和12时相对应的值。
高斯-勒让德积分的例子。我们来计算一下这个积分作为例子：
+1
x -1
I=e dx ? 式（2-30）
表2-4为多个计算量的值(其中，2栏和3栏为Xi 点上被积函数的值)；表2-3为横坐标和权重系数。与1x 111
e dx e e 2.350402387
+--=-=?
的准确值比较时，我们发现，即使当n= 4 逼
近时，用高斯－勒让德积分法所取得的值仍然精确到7个有效数字。
理论上，指数函数是自变量x 的无穷级数23
x
x x x e e ......
1!2!3!=++++
当我们用n=4做高斯－勒让德积分法时，我们把函数近似为一个最高级是(2n-1)=7级的函数。在这个范围内，我们忽略任何幂高于X8的项。也就是说，如果阶n=5时，结果会更精确。
在这种情况下，最高级数为X9的项也包括在内。 从表2-4的第二部分中我们看到，这时的计算结果为2.3504023。本质上，这个值是单精度计算可获的最佳值。 要想在不是1±的范围内得到相同的值，我们可以对变量做出如下调整：
a
1
x ay a
1
e dx a e dy ++--=?
? （2-32）
现在这个积分与式（2-31）的形式相同，故可以进行相同的计算。方框2-3为具体算法。如我们所预想的，当a 超过1时，计算的精度会降低。
现在，我们来做一个有关高斯积分法精确度的有趣练习——计算积分
1
81
x dx +-?
。当n=4时，
所得结果为0.21061227。比之于精确值2/9，它有一个5%的误差。造成这种较大误差的原因很明显，这个方法把被积函数近似为了（2n-1）=7 阶的多项式。要是我们代入n=5的话，我们便得到与单精度计算的确值相同的结果----0.22222221。这就是高斯积分法的重要局限：为了获得某个精度所需多项式的阶数很大程度上取决于被积函数的性质。总体上来说，我们不易确定某个函数和精度所需的多项式逼近的阶数。另外，这种方法本身也没有给出一个计算误差的方法。所以，如果不做独立的验算，很难保证所获值的可靠性。
高斯-勒让德积分：以级数求和的形式来求解数值积分，并举例说明了当n 取不同值时对积分精确度的影响。高斯积分法的重要局限：为了获得某个精度所需要的多项式的阶数很大程度上取决于被积函数的性质。 2-5 蒙特卡罗积分
用讨论过的数值积分法，被积函数f(x)只能在横轴的规律间隔点上求值。对于矩形，梯形，和森逊法则，f （x ）的求值是在N+1个等间距的点X 0, X 1, …..Xn 上进行的。在做高斯－勒让德积分法时，点集Xi 取自于某个阶数的勒让德多项式在 [－1，＋1]内的解。本节我们采用不同的思路，利用数据采样的方法来实现这样做得好处是，它简化处理了多维积分和某些非正常积分（在书的下一部分我们会回过来讨论这一问题）。 对于一维常态被积函数（在这部分我们将继续作为例子使用），这种方法却不比我们见过的其他方法有效了。
为了讨论的方便，我们选取积分区间为[-1，+1]。这样做是因为多数均匀分布随机数生成器都设定在这一范围。对于那些下限不为0或上限不为1的积分来说，我们很容易做出下面的转换：
1
x (x a)b a
→
--
这样一来，实际的数值积分运算将在[0,1]区间进行。这就与式（2-19）的情况就差不多了（除了这里的转换只限于积分在[-1，+1]）。
蒙特卡罗积分法不要求提前选取网格点，而是选取随机样本。如果样本确实是随机的话，点X 1,X 2……, n x 会以相同概率分布在[-1,+1]区间内。如果有N 个这样的点，那么相邻的两个点之间的平均距离应为:
arg h 1/(N 1)1/N l e
N
=-→ （2-33）
图2-5：直方图显示了在[0，1]区间内用7-2栏内的减数法生成的5000个标准随机数的分布 此时积分的值为
N
1
[0,1]
i 0
N i 1
1I f (x)dx f (x )N →∞==∑→? （2-34） 在x 轴上选取点的时候，我们可以用一个平均分布随机数生成器（如：7-1节讨论的生成器）。 如果样本确实是随机的，那么在区间1
1[x x,x x]22
-+ 内选出的点的数量与x 的值是无关的。换句话说，若是把[0,1]划分为若干个的等大的子区间，每个子区间内点的数量在统计涨落中是一样的。 统计学的基本概念告诉我们，这种情况只出现在样本够大的时候（如表2-5所示的一个有5000个随机点的样本）。
一个存在于蒙特卡罗法中迭代法的问题是怎么知道N 在什么时候大到使得式（2-34）有效。为了判定有效性，我们可以用一种实际的方法，并采取迭代策略。这与我们在处理辛普森法则时做的一样。我们先取N1个点的样本并且代入式(2-34)求积分的值。我们再用相同的步
骤取另外N1个不同的点并且也代入式(2-34)求积分的值。如果两次运算的结果在一定的可以接受的误差范围E内相同，那么两个样本就有可能足够大了。若还要提高精度，我们可以结合两个样本的结果，把它们的平均作为最终结果。而当误差太大了的时候，我们可以把两组合并为一组样本,其点数为N2 =2N1，同时取N2个不同的点，再把两个新的样本相互比较。反复这个步骤直到理想的精度，或者在反复的次数超过了计算允许的最大迭代数的时候，终止这一步骤。
在辛普森法则的例子中，我们先取一个有少量点的样本(N= 2~6)。我们知道，我们的方法需要的随机样本应该有均匀分布的点，但是如此小的样本中的点是不可能这样分布的。
严谨的人可能想从一个更大的样本出发，如：2
1N ~10 。我们可以以标准概率函数为例来检验随着抽样数据点的增加蒙特卡罗积分的收敛性
2
x
t x
e dt +--?
（2-35）
这个我们以后再§6.1
中讨论。为了简化起见，我们应该首先考虑当x=1的情况。此外，由于被积函数关于x=0对称，我们可以把积分形式改为
2
1
t 0
e dt -
这使我们能够使用无需做任何变换的[0，1]区间内标准分布的随机数生成器。这个值和标准正态分布曲线下方所包围的面积相等，并且在标准的平均值标准差之内。根据表中的数据我们得到数值
2x
t x
1
e dt-)=0.68268952
+--
表2-5：