Simpson算法及其推广形式
摘要：本文研究了辛普森公式的数值积分的计算方法问题，并且更进一步研究了变步长复化的辛普森公式和二重积分的辛普森公式的问题。首先是对
一维辛普森公式和变步长复化辛普森公式以及二维辛普森公式的推导及
其算法，进行误差分析，并且列举了实例。然后，对辛普森公式进行改
进，这里的改进最主要是对辛普森公式的代数精度进行提高，从而使辛
普森公式对积分的计算更加精确。另外，还研究了辛普森公式的推广形
式。最后，在结论的当中列举了一个例子。
关键词：辛普森公式算法改进推广形式二重积分的辛普森公式
Abstract：This paper first studies the calculation methods of the numerical integration in simpson formula, and then study of the long-simpson
formula and the double integral simpson formula problem. First, study the
algorithm and derived of one-dimensional simpson formula and
step-change in simpson formula, as well as two-dimensional simpson
formula, and then analysis the error. Finally , list the example. In this ,
improve the simpson formula. This improved the most important is to
incre ase the simpson formula’s accuracy of algebra. Besides, we study the
simpson formula’s promotion of forms. At the last, we list a example in
the conclusion.
Key word:The simpson formula, Algorithm, Improve, Promotion of forms, The simpson formula of the two-dimensional integral.
1 引言
辛普森公式主要的研究数值积分（numerical integration）的。何谓数值积分呢？其是求定积分的近似值的数值方法。即用被积函数的有限个抽样值的离散或加权平均近似值代替定积分的值。求某函数的定积分时，在多数情况下，被积函数的原函数很难用初等函数表达出来，因此能够借助微积分学的牛顿-莱布尼兹公式计算定积分的机会是不多的。另外，许多实际问题中的被积函数往往是列表函数或其他形式的非连续函数，对这类函数的定积分，也不能用不定积分方法求解。由于以上原因，数值积分的理论与方法一直是计算数学研究的基本课题。对微积分学做出杰出贡献的数学大师，如I.牛顿、L.欧拉、C.F.高斯等人也在数值积分这个领域做出了各自的贡献，并奠定了它的理论基础。
构造数值积分公式最通常的方法是用积分区间上的n 次插值多项式代替被积函数，由此导出的求积公式称为插值型求积公式。特别在节点分布等距的情形称为牛顿-柯茨公式，例如梯形公式与抛物线公式就是最基本的近似公式。但它们的精度较差。龙贝格算法是在区间逐次分半过程中，对梯形公式的近似值进行加权平均获得准确程度较高的积分近似值的一种方法，它具有公式简练、计算结果准确、使用方便、稳定性好等优点，因此在等距情形宜采用龙贝格求积公式。当用不等距节点进行计算时，常用高斯型求积公式计算，它在节点数目相同情况下，准确程度较高，稳定性好，而且还可以计算无穷积分。数值积分还是微分方程数值解法的重要依据。许多重要公式都可以用数值积分方程导出。
在数值积分的研究辛普森是一位很重要的数学家。辛普森（Thomas Simpson，公元1710年8月20日─公元1761年5月14日）是英国著名的数学家，他生于英格兰列斯特郡，并卒于同地。他的父亲是一位纺织工人，所以他主要靠自己力学成材，而他的第一份工作也是纺织。他对数学的兴趣最初是由一次日蚀所引发起的，他在一位占卜师的指导之下，他学会了算术和基本的代数。其后他放弃了纺织的工作，而当了一间学校的司阍。凭借他刻苦和持久的努力，他证明了他在数学方面的能力，以致于1735年他能够解决了数个有关微积分的问题，1737年他便开始撰写有关数学的文章。在1754年他成为了「淑女日记」一书的编辑，及后他到了伦敦的乌尔威治并出任数学教授一职，直至逝世。
辛普森最为人熟悉的贡献是他在插值法（Interpolation）及数值积分法（Numerical Method of Integration）方面，事实上他在概率方面也有一定的
工作，他在1740年推出了他的「机会的特性和法则」（The Nature and Laws of Chance ），而大部份他在这方面的结果也是建基于棣美弗早期的结果。另外，当时有一群讲师巡回在伦敦咖啡屋讲学，而辛普森是当中最突出的一位。他专研有关「误差理论」（Theory of Error ），并且意图证明数算平均数比单一观察较佳。Simpson 公式就是他的代表的定理。
在辛普森之后，很多后人都在其的基础上，不断完善simpson 公式，使其公式越来越完善。其中，华罗庚，王元的著作数《值积分及其应用》中，就是用数论的方法研究辛普森公式。在其研究中，假设()f x 是一个在[],αβ内定义了的函数，以后如果用到几次微商，便假设定()f x 有几次微商。我们用Euler 求和公式（详细参考该书第一章的内容）来推出普通数值积的Simpson 法。从而得出下面的内容：
命
12
33
s t r R R R =+ （其中t R ，r R 参照该书第二章的1
和2节） 则得
()1
1
01102246n n s n i i i i R f t dt y y y y n α
β
βα--+==?
?
-=-
+++ ???
∑∑?， 而且simpson 公式的余项[由（1）与（2）]为
()()3
22''03
''2201211233633611232n s n
b x b x R f x dx n n b x b x f x dx n n βαβααβαβαα????- ? ?--?????? ?=-++?+ ? ? ????
? ?
??
-??-??
????=+-+ ? ? ? ????????
??? 其中
()()()1
13
''2012112n l n l n
f t dt y y y n
b x f x dx n n α
β
βαβαβαα-=-??
-
++ ???
--??
?
???=-+ ? ? ????
???∑?? （1）
()3
1''2
122
11242n r R b x f x dx n n βαβαα--
-?-???
???
?=-++- ?
? ? ???
??????
?
（2） 如果()()IV f x a x b ≤≤存在，由部分积分可得
()()''220'''330122122n
n
b x b x f x dx n b x b x f x dx n
n βααβα
βαα??-????+-+= ? ? ??????
?-??-?
???-
+-+ ? ? ??????
???
（此处用了()33100,02b b ??
== ???
）
()()()()()()()'''4
42
440''''''2
44024412212219601222n
o
n IV n IV b x b x f x n
n b x b x f x dx
n n f f n b x b x f x dx
n
n b x b x n βαβααβαβααβααββαβααβα-??-?
???=-
+-+ ? ? ????
???-??-???
???++-+ ? ? ? ??????
???-=--??-??????++-+ ? ? ? ??????
???-??
=+ ?????0112960n
IV f x dx n βαα??-????--+ ? ? ???????
?
即
()5
440111232960n
IV s R b x b x f x dx n n βαβαα-??-??
????=+--?+ ?
? ? ???
??????
?、 定理1.1 如果()(),IV f x M x αβ≤≤≤，则 ()5
44
1802s
M
R n βα-≤
??
证. 由于 ()4
40432
1432
20432
432
1122960111111222(2412247201261236096011111122224122472012612360960n b x b x dx x x x x x x n dx x x x x x x ??+-- ???
???
?????+++ ? ? ? ??????? ?=-+-+-+-++ ? ???
???
?????+++ ? ? ? ??????? +-+-+-+-++ ??
??121664
)6
63180218021802
dx n n ?????=+= ???????
故得定理。 定理1.2 如果
()
''f x 是单调非负递减函数，则
()()3
''3
324s M
R f n βαα-≤
.
证. 由
()1122220
102212x x b x dx dx ??
=-+= ???
?? 及对1
02
η<<
，常有 ()2
20
2
20
2320
12211112212226312224108b x b x dx x x x x dx x x dx η
η
η
ηη????+- ? ????
???
????=-+++-++ ? ? ? ???????
??=
-=-≤ ?
??
?
??
由
()3
''22011232n
s R b x b x f x dx n n βαβαα-??-??
????=+-+ ?
? ? ???
?????
?? 以及第二中值公式可得
()()()
()
3
''2203''3
11232324s R f b x b x dx n M
f n
ξβααβαα-?????
?≤+- ? ? ??????
?-≤
?
即得定理1.2.
【】参考文献：华罗庚，王元， 数值积分及其应用 科学出版社 1963年 前人的经验都是很多的，这里只是列举了其中一个。现在就来介绍一下一些预备知识。首先，来介绍一下代数精度。
代数精度是用来衡量求积公式精确程度的一个概念，代数精度越高则求积公式就越准确。
定义1.1 如果一个一般的求积公式对于所以次数不大于m 的多项式()f x 都准确成立，即
()()0
n
b
k k a
k f x f x λ==∑?
而对于次数为m+1的某一个多项式()f x 并不准确成立，则称该求积公式具有m 次代数精度。
根据定义1.1，插值型求积公式()()0
n
b
k k a
k f x f x λ=≈∑?至少要具有n 次代数精度。
由于任一m 次多项式
()1110m m m m m p x a x a x a x a --=++???++ 是由k
x ()
0,1,2,,k m =???叠加而成的，因此可以利用函数
()k f x x =()0,1,2,,k m =???来判断求积公式的代数精度。
定理1.3 表明n+1个节点的求积公式是内插求积公式至少具有n 次代数精确度，也可能更高。
判断求积公式代数精度的步骤：
1.在内插求积公式两端分别代入最简单的代数多项式n x x Λ,,1进行计算，即代入f(x)=1, f(x)=x,…,f(x)=x n 计算。
2 .判断两端数值相等时的代数多项式的最高次数，则可以确定出代数精度。即对于任意不超过n 次的代数多项式都准确成立，而对某一个n ＋1次代数多项式不成立，则其代数精度为n 。
例1. 确定积分公式[])1()(3)(23
1)(211
1
f x f x f dx x f ++≈?-中的待定常数，并确
定其代数精度。
解 令f (x )＝1，
左边＝21)(1
11
1==??--dx dx x f 右边＝2)132(3
1
=++
令f (x )＝x
左边＝??--==11110)(xdx dx x f 右边＝[]13231
21++x x
令f(x)＝x 2， 左边＝321
12=
?-dx x 右边＝)132(3
12
2
21++x x 得：
[]???=-=???-==??????
?=
++=++1266
.06899.05266.02899.032)132(3
1013231
21
21222121x x x x x x x x 或得
令f (x )＝x 3，
左边＝ 右边＝0.6106或 0.3494 ∴左≠右 ∴该求积公式具有2次代数精度。 下面，再来研究一下拉格朗日插值公式。
设连续函数y = f （x ）在[a, b]上对给定n + 1个不同结点：
x 0, x 1, …, x n 分别取函数值
y 0, y 1, …, y n
其中 y i = f (x i ) i = 0, 1, 2,…, n 试构造一个次数不超过n 的插值多项式
n n n x a x a x a a x P ++++=Λ2210)(
使之满足条件
i i n y x P =)( i = 0, 1, 2,…, n
先求n 次多项式l k (x ) k = 0, 1,…, n
使
??
?≠==i
k i k x l i k ,
0,
1)( （3）
若做出这样的)(x l k ，则P n (x )的次数≤n ，另外，由（3）
i i n
k k k i n y x l y x P ==∑=)()(1 i = 0, 1, 2,…, n
即()n P x 满足插值条件n i y x P i i n ,,2,1,0,)(Λ==。于是问题归结为具体求
出基本插值多项式l k (x )。 根据（3）式，x k 以外所有的结点都是l k (x )的根，因此令
∏≠=+--=-----=n
k
j j j n k k k x x x x x x x x x x x x x l 01110)
()())(())(()(λλΛΛ
又由l k (x k ) = 1，得：
1
1
3=?
-dx x
)
())(())((1
1110n k k k k k k k x x x x x x x x x x -----=
+-ΛΛλ
所以有：
∏
≠=+-+---=----------=
n
k
j j j
k j n k k k k k k k n k k k x x x x x x x x x x x x x x x x x x x x x x x x x l 011101110)
())(())(()())(())(()(ΛΛΛΛ
代入（3）式即得()n P x 的表达式
k n k n
k
j j j
k j k n
k k n y x x x x y x l x P ∑∏
∑=≠==???
?
?
??--==000)()( （4）
（4）式称为拉格朗日插多项式。
我们把差f (x )－()n P x 称为用插值多项式()n P x 代替f (x )的余项，误差或
插值余项，记为：
)()()(x P x f x R n n -=
下面讨论余项的形式。
设x 是区间[a, b]中任意固定的数，若x 是插值节点x i ，则0)(=i n x R ，i =
0, 1, 2,…, n ，如果x 不是节点，考虑如下函数：
)()
()
()()()(11x R x t t P t f t n n n n ++-
-=ωω?
其中
∏=+-=n
k k n x x x 0
1)()(ω
由插值条件n i y x P i i n ,,2,1,0,)(Λ==知
0)(=i x ? i = 0, 1,…, n
而
0)()
()
()()()(11=-
-=++x R x x x P x f x n n n n ωω?
所以，?(t )在[a, b]上有n + 2个互异零点。应用罗尔定理，?’(t )在?(t )的每两个零点间有一个零点，即?’(t ) 在[a, b]上至少有n + 1个零点，进一步对?’(t )应用罗尔定理可知，?" (t )在
[a, b]上至少有n 个零点，继续上述讨论就可推得?(n +1)
(t )在[a, b]内至少有一
个零点，记之为ξ，即
0)()1(=+ξ?n
因为P n (t )为不高于n 次的多项式。所以
)!1()(0
)()1(1)1(+==+++n t t P n n n n ω
于是
)
()
()!1()()(01)
1()1(x x R n f
n n n n ++++-
==ωξξ?
移项、整理可得
)()!
1()
()(1)1(x n f x R n n n +++=ωξ
综上所述，得以下余项定理：
定理1.4： 设f (n )(x )在区间[a, b]上连续，f (n + 1)
(x )在[a, b]上存在，x 0, x 1,…,
x n 是[a, b]上互异的数，记插值问题n i y x P i i n ,,2,1,0,
)(Λ==的余项为
)()()(x P x f x R n n -=，那么，当x ∈ [a, b]时，有如下估计
)()!
1()
()(1)1(x n f x R n n n +++=ωξ
其中
∏=+-=n
j j n x x x 1
1)()(ω
2 辛普森公式及其算法
2.1 有限区间上的一维辛普森公式
2.1.1 利用插值多项式导出求积公式
设()f x 为被积函数，[],a b 为积分区间，01,,,n x x x ???为[],a b 内的n+1个互异的点，记()n L x 为相应的拉格朗日插值多项式，()n R x 为相应的插值余项，那么我们有
()()()n n f x I f R f =+ （2.1.1）
将上面的（2.1）式两边同事积分得
()()()b b b
n
n a
a
a
f x dx I f dx R f dx =+???
如果我们取
()()b
b
n a
a
f x dx I f dx ≈?
? （2.1.2）
那么截断误差为
()b
n a
R f dx ?
利用拉格朗日插值公式（4），我们有
()()()0n
n k k k I x f x LB x ==∑
其中()k LB x ，0,1,,k n =???，是基函数。从而上面的（2.1.2）式可以改成
()()()0
n
b
b k
k
a
a
k f x dx f x LB x dx =≈∑?? （2.1.3）
由于各()k f x 都是常数，()k LB x 是多项式函数，所以从理论上讲，上面的（2.1.3）的右边是可以计算的。记
()b
k k a w LB x dx =? （2.1.4）
由上面的（2.1.3）即得
()()0
n
b
k
k
a
k f x dx f x w
=≈∑? （2.1.5）
从而得到一种求数值积分的方法。上面的（2.1.4）和（2.1.5）被称为插值公式导出的求积公式。
回归第十章关于多项式插值的结论可知，由于任意次数不超过n 的多项式与它的任意n+1个基点的插值多项式恒等，再由求积公式的代数精度定义，我们立即得到：由n+1个基点的拉格朗日插值多项式所形成的求积公式的代数精度至少是n 。
为了强调上面的（2.1.5）是利用n+1个基点的插值多项式导出的求积公式，我们把公式中的各k w 改写成(),w n k ，从而有
()(),b
k a w n k LB x dx =? （2.1.6）
()()()0
,n
b
k
a
k f x dx f x w n k =≈∑? （2.1.7）
2.1.2 牛顿-科特斯求积公式的推导
虽然从理论上讲，可以利用上面的（2.1.6）和（2.1.7）完成()f x 在区间[],a b 上的数值积分的近似计算，但是利用（2.1.6）计算各(),w n k 仍然是一件比较麻烦的事。牛顿-科特斯求积公式的目的就是利用01,,,n x x x ???是区间[],a b 的一个等分分划的特点进一步简化（2.1.6）中的各(),w n k 的计算。
假如我们将积分区间[],a b 划分为n 等分，记
()
0,0,1,,k b a h n
x a x a kh k n x a th -?
=
??=??
=+=????
=+?
（2.1.8）
那么我们有
()(),0,1,,,,0,1,,0
j k j x x t j h j n
x x k j h k j n dx hdt x a t x b t n -=-=????
?
-=-=????
?
=??=?=?
=?=?
? （2.1.9）
再注意到对0,1,,k n =???，拉格朗日插值的基函数()k LB x 可表示为
()0n
j k j k j
j k
x x LB x x x =≠-=-∏
所以利用上面（2.1.8）和（2.1.9）给出的记号约定以及相应的积分变换，我们有
()000,n n b n j a j j k j j k
j k
x x t j
w n k dx dx x x k j
==≠≠--==--∏∏??
从而可以得到
()()
()000
,n
n
n
j j k
j j k
h
w n k t j dt k j =≠=≠=
--∏?∏ （2.1.10）
再注意到
()()()()()()
()011121!!
n
j j k
n k
k j k k n k k n k =≠--=-????--???-+????????
=--∏
所以由上面的（2.1.10）可以得到
()()()()00
1,!!n k
n n j j k
h w n k t j dt k n k -=≠-=--∏? 再利用
()
b a h n
-=
立即可以得到
()()()()()00
1,!!n k
n n j j k
w n k b a t j dt nk n k -=≠-=---∏? (2.1.11) 在上面的(2.1.11)中，记
()()()()00
1,!!n k
n n j j k
C n k t j dt nk n k -=≠-=--∏? （2.1.12） 则上面的（2.1.11）可表示为
()()(),,w n k b a C n k =-
从而前面的（2.1.7）可以表示为
()()()()0
,n
b
k
a
k f x dx b a C n k f x ==-∑? （2.1.13）
上面的（2.1.12）和（2.1.13）称为牛顿-科特斯求积公式，利用（2.1.12）得到的各(),C n k 称为牛顿-科特斯系数。虽然看上去利用（2.1.12）式计算为牛顿-科特斯系数也还是有些麻烦，但是与前面的（2.1.4）对比可知，各(),C n k 与具体问题无关，可以把它们“一次性地”计算出来，并制成表。这样，人们就可以通过查表得到这些值，从而可以免去这一部分的计算。这样一来，上面的（2.1.13）就真的称为机械求积公式了。
虽然在现代的计算机条件下牛顿-科特斯求积公式和牛顿-科特斯系数的作用已经很有限，但是牛顿-科特斯求积公式所体现出来的独具创造性的思维方式和大胆的探索精神为科学计算领域中的研究树立了光辉的榜样。
2.1.3 辛普森公式
在牛顿-科特斯求积公式中，如果取n=2时，那么k 可以取0，1，2，此时形成的的求积 公式就是辛普森公式，利用（2.1.10），我们可以得到
()()
()()()()
()()()()
()()2
2
1
1
00
1
0112,01220!2!6
142,10221!1!
6112,20122!0!
6
C t t dt C t t dt C t t dt -=--=
??-=--=
??-=
--=
????? 所以有
()()()()14
16626b
a
a b f x dx b a f a f f b ?+???=-+
+ ???????
?
（2.1.14） 【】参考文献：甄西丰编著 实用数值计算方法 清华大学出版社 2006年
2.1.4误差分析
应用Newton-Cotes 型求积分公式（2.1.13）计算定积分()()b
a I f f x =?dx
时，一方面由于它是由式（2.1.1）去掉余项()n E f 得到的，因而产生了离散误差()n E f ；另一方面，由于计算机的字长是有限的，函数值可能带有误差，并且计算()n I f 还会有舍入误差。
关于离散误差，我们有下面的定理。
定理2.1 设n 为偶数且()f x 在[],a b 上有n+2阶连续导数，则Newton-Cotes 型求积分公式（2.1.13）的离散误差为
()()()
()()()2320
12!
n n n
n h f E f t t t n dt n η++=
-???-+? (),a b η∈
若n 为奇数，且()f x 在[],a b 上有n+1阶连续导数，则
()()()
()()()1220
11!
n n n
n h f E f t t t n dt n ξ++=
-???-+?
(),a b ξ∈
特别，n=2时，Simpson 公式的离散误差为
()()()54290h E f f η=- 2
b a
h -=
由于Newton-Cotes 型求积分公式的误差公式证明比较繁杂，下面给出辛普森公式的误差的证明。
证明：首先考虑构造一个三次差值多项式()3p x 满足条件：
()()()()333,,22a b a b p a f a p b f b p f ++????=== ? ?????和''322a b a b p f ++????= ? ?????
。可以
证明满足上述条件的()3p x 与()f x 的误差为
()()()()()()2
434!2f a b f x p x x a x x b ξ+??
-=--- ???
（2.1.15）
其中a b ξ≤≤。
对（2.1.15）式两边分别进行积分得
()()()()()()2
434!2b b
b
a
a
a
f a b f x dx p x dx x a x x b dx ξ+??
-=--- ???
??
?
又因为()3p x 为三次代数多项式，而辛普森公式的代数精度为3.所以有
()()()()()()()3333141662614
16
626b
a
a b p x dx b a p a p p b a b b a f a f f b ?+?
??=-++ ?
??????
?+???=-+
+ ???????
?
从而有
()()()()()
()
()()2
414
1
6
6264!
2b
a
b
a
a b f x dx b a f a f f b f
a b x a x x b dx ξ?+???--+
+ ???????
+??
=--
- ???
?
?
据提设，()
()4f
ξ在[a,b]上连续，而()()2
02a b x a x x b +??---≤ ??
?，由积分中
值定理可知，在[a,b]中存在一个点η，使
(
)
()
()()()()()()()()()
2
42
45
44!
214!22880
b
a
b a f a b x a x x b dx
a b f x a x x b dx b a f ξηη+?
?--- ??
?+?
?=--- ??
?-=-
?
? 因此
()()()()()(
)
()
5
414166262880
b
a
a b f x dx b a f a f f b b a f η?+???--++ ???????
-=-
?
得证。
【】参考文献：韩国强编著，数值分析，华南理工大学出版社，2005年
2.2 有限区间的变步长复化辛普森公式
2.2.1 复化辛普森公式
假设()012212221,,,,,,,,k k n k x x x x x x x --??????把积分空间[],a b 划分为2n 等分，我们也可以认为是其中的()02221,,,,,n k x x x x -??????把区间划分为n 等份，并且21k x -就是第
k 个子空间()221,k k x x -????的中点。
记k I ，1,2,,k n =???为()f x 在第k 个子空间()221,k k x x -????应用辛普森公式所求
得的积分值，则有
()()
()()2122146k k k k b a I f x f x f x n ---?
?=
++?
?
记 21
k n
n k k S I ===∑
则有
()()()()1
20221211426n n n n k k k k b a S f x f x f x f x n --==-??
=+++????
∑∑ （2.2.1）
上面的（2.2.1）式称为复化辛普森公式，虽然我们也可以直接编写2n S 的计算机程序，但是没有必要那么急，而是我们要改进一下变步长复化辛普森公式的性能。
2.2.2 变步长复化辛普森公式
假设()012212221,,,,,,,,k k n k x x x x x x x --??????把积分空间[],a b 划分为2n 等分，那么我们可以得到下面的3个不同的积分值：
（1） 利用2k x ，0,1,2,,k n =???这n+1个点处的函数值和复化梯形公式计算出
n T ；
（2） 利用j x ，0,1,2,,2k n =???这2n+1个点处的函数值和复化梯形公式计算出
2n T ；
（3） 利用复化辛普森公式计算出2n S ； 显然，2n S ，n T ，2n T 之间应该存在一定的关系。 在这里，我们应该先知道复化梯形公式
()()()10111
22n n n j j b a T f x f x f x n -=??-=++?
???
∑ （2.2.2）
另外
()()()21202111
222n n n j j b a T f y f y f y n -=??-=++?
???
∑ （2.2.3） 这里的这两个式子的证明就不给出了。
按照（2.2.2）和（2.2.3）式方括号内的表示形式整理（2.2.1）式的方括号内的数学表达式，不难得到
()2126n b a S n
-=∑+∑ （2.2.4）
其中
()()()()1
1022121
1
2244n
n n k k k k f x f x f x f x --==∑=+++∑∑
()()()1
20221
2n n k k f x f x f x -=∑=++∑
上面的1∑，2∑还可以进行进一步调整为
()()()()()()2110211
2022111
42211222n n j j n n k k f x f x f x f x f x f x -=-=???∑=++?
?????????
∑=++???
???
∑∑ （2.2.5）
利用前面的（2.2.2）和（2.2.3）式以及上面的（2.2.5）,我们可以得到
12241
,6363
n n b a b a T T --∑=∑= （2.2.6） 再由（2.2.1）和（2.2.6）式即得到
2243n n n T T
S -= （2.2.7）
按照我们习惯记法，取2k n =，则有122k n +=，利用[]1S k +表示由（2.2.7）式所得到的2n S ，利用复化梯形公式得到的梯形序列[][]0,1,,T T ???我们有
[][][]()141/3S k T k T k +=+- （2.2.8）
所以由梯形序列[][]0,1,,T T ???我们可以得到了一个新的序列式[][]1,2,,S S ???称之为辛普森序列。
2.2.3 复化辛普森公式的算法实现 （1）复化辛普森公式的计算步骤
1）确定步长()/h b a N =-（N 为等份数）。
()122,0S f a h S =+=。
2）对1,2,,1k N =???-计算
()112,S S f a kh h =+++ ()22S S f a kh =++ 3） ()()12426S h f a S S f b =+++???? （2）算法流程图
复化辛普森公式的算法流程图见图1。
（3） 复化辛普森公式的matlab 程序
按照（2.2.1）编写复化辛普森求积函数（函数名：s_quad.m ） function I=S_quad(x,y) %复化辛普森求积公式，其中
%x 为向量，被积函数自变量的等距节点；
%y 为向量，被积函数在节点处的函数值。 n=length(x);m=length(y) if n~=m
error(‘the lengths of X and Y must be equal ’); return; end
if rem(n-1,2)~=0 I=T_quad(x,y); Return; End
N=(n-1)/2;h=(x(n)-x(1))/N;a=zeros(1,n); For k=1:N
a(2*k-1)= a(2*k-1)+1; a(2*k)= a(2*k)+4; a(2*k+1)= a(2*k+1)+1 end
I=h/6*sum(a*y');
下面给一个例子，用复化辛普森公式求积2
1
1x e dx --?，在积分区间中点与点之
间的间隔取为0.1。 解：输入
x=-1:0.1:1;y=exp(-x.^2); I=S_quad(x,y) 得到
I=1.4936
【】参考文献：薛毅编著，数值分析与实验，北京化工大学出版社，2005年
2.3 二重积分的辛普森公式
我们上面讨论的辛普森公式的求积分方法，稍经修改就可以用来计算重积分。
(),R
I f x y dxdy =?? （2.3.1）
的计算法，此处，假设积分区域为矩形域：
(){},|,R x y a x b c y d =≤≤≤≤。
给定正整数n 和m ，取步长
()(),h b a n k d c m =-=-。
把积分I 写成
(),b
d
a
c
I f x y dxdy =?
?
对积分 (),d c
f x y dy ?
视x 为常数，令,0,1,,2j y c jk j m =+=???。应用复合辛普森公式得
()()()()()()()
()102212114
4
4
,,2,4,,3,,,180
m m d
j j m c
j j k f x y dy f x y f x y f x y f x y d c k f x c d y μμ--==?
?=+++??
??
-?-
∈?∑∑?
于是
()()()()()()
10212121444
2,,334,,33,180
m b b j a a j m b b
j m a a j b
a
k k I f x y dx f x y dx
k k f x y dx f x y dx d c k f x dx y
μ-=-==+++-?-
?∑??∑???
再令,0,1,,2i x a ih i n =+=???，则对每一0,1,,2j m =???，有
()()()()()()()
()
102212114
4
4
,,2,4,,3,,,180
n n
b
j j i j i j n j a
i i j j j h f x y dy f x y f x y f x y f x y f y b a h
a b x
εε--==??
=+++??
??
?--
∈?∑∑?
从而
()()()()()()
()()
1002021011
111
2002221
11
1`1
2122211
1
[,2,4,9,2,4,8,2,n n i i i i m m n n j i j j j i m n m i j n j j i j hk
I f x y f x y f x y f x y f x y f x y f x y f x y --==---===---===≈+++++++∑∑∑∑∑∑∑∑