辛普森求积公式
摘要
在工程实验及研究中，实际工作中，变量间未必都有线性关系，如服药后血药浓度与时间的关系；疾病疗效与疗程长短的关系；毒物剂量与致死率的关系等常呈曲线关系。曲线拟合是指选择适当的曲线类型来拟合观测数据，并用拟合的曲线方程分析两变量间的关系.可以说，曲线拟合模型与我们的生活生产密切相关.
本课题着重介绍曲线拟合模型及其应用，其中包括它的基本思想、模型的建立、以及具体应用.为了更好的了解曲线拟合模型，可以将它分为线性与非线性模型，在模型建立的基础上我们可以用最小二乘法来解决一些我们日常所应用的问题.
关键词曲线拟合；线性与非线性模型；最小二乘发
目录
引言 (1)
第一章曲线拟合 (2)
§1.1 基本思想及基本概念 (2)
§1.1.1 方法思想 (2)
§1.1.2几个基本概念 (2)
§1.2辛普森算法基本定义及其应用 (4)
§1.2.1辛普森求积公式的定义 (4)
§1.2.2辛普森求积公式的几何意义 (5)
§1.2.3辛普森求积公式的代数精度及其余项 (5)
§1.2.4辛普森公式的应用 (6)
第二章辛普森求积公式的拓展及其应用 (7)
§2.1 复化辛普森求积公式 (7)
§2.1.1问题的提出 (7)
§2.1.2复化辛普森公式及其分析 (7)
§2.1.3复化辛普森公式计算流程图 (8)
§2.1.4复化辛普森公式的应用 (9)
§2.2 变步长辛普森求积公式 (10)
§2.2.1变步长辛普森求积公式的导出过程 (10)
§2.2.2变步长辛普森求积公式的加速过程 (12)
§2.2.3变步长辛普森求积公式的算法流程图 (13)
§2.2.4变步长辛普森公式算法程序代码 (14)
§2.2.5变步长辛普森求积公式的应用 (14)
§2.2.6小结 (14)
§2.2.7数值求积公式在实际工程中的应用 (14)
参考文献 (16)
附录A (17)
附录B (18)
附录C (21)
引言
辛普森是英国数学家.1710年8月20日生于波士沃希；1761年5月14日卒于波士沃希.在定积分近似计算中，以他的姓来命名的“辛普森公式”，虽早在他之前牛顿的学生柯特斯（Cotes）和斯特林就已经得出了（包括一些更高阶的近似公式），但真正广泛地为人所知并加以应用，则是1743年辛普森重新发现之后的事了.辛普森的工作使牛顿的微积分学说得到了进一步完善.在我们的日常生活中计算积分与我们的生活生产密切相关.所以掌握数值积分方法是学生储备知识能量的武器.
数值积分的一个基本的计算策略，用易于积分的简单函数来逼近曲线)
y .
f
(x 简单曲线下面的面积近似等于)
f下面的面积.如果涉及初等函数的积分找不到其
(x
他由初等函数构成的解析表达式，或者只在一些离散的x点上知道函数的值，在多数情况下，被积函数的原函数很难用初等函数表达出来，因此能够借助微积分学的牛顿-莱布尼兹公式计算定积分的机会是不多的.那么就必须对定积分进行数值逼近.
数值积分实现是将整个闭区间]
f进行
(x
,
[b
a划分为N个小段，在每个小段上对)
低阶分段多项式逼近.对每个小段上的逼近多项式积分时，就得到基本公式.基本公式只涉及足够多的))
f
(x
x对来定义分段多项式的某一段，将此公式应用到N个小段并
,
(
把结果相加得到复合公式，或称为扩展公式. 在一个小段中节点的位置和数目决定了基本公式的很多重要特性.当节点均匀分布时，所有的积分公式称为牛顿—柯特斯公式.例如，梯形、辛普森、柯特斯求积公式等.
经典辛普森求积公式来源于Lagrange插值多项式的应用，它的代数精度高达3阶，其形变后的代数精度高达4阶，且二者都具有良好的稳定性与收敛性，从而提高了计算效率及准确度，是定积分近似计算常使用的方法，一直是理工科大学生必修的内容. 下面将给出具体辛普森求积公式的具体思想以及其算法程序设计并给出将其拓展后在实际工程问题中的应用.
第一章 辛普森求积公式的理论
实际问题当中常常需要计算积分，有些数值方法，如微分方程和积分方程的求解，也都和积分计算相联系.依据人们所熟知的微积分基本定理，对积分
dx x f I b
a ?=)(只要找到被积函数)(x f 的原函数)(X F ,便有下列牛顿-莱布尼茨公
式：?-=b a
a F
b F dx x f )()()(，但实际计算dx x f b
a
?)(往往遇到一些困难，如: 1))(x f 的
原函数不能用初等函数表示，故不能用牛顿-莱布尼茨公式计算.2) 虽然找到了
)(x f 的原函数, 但因表达式过于复杂而不便应用牛顿-莱布尼茨公式.3) )(x f 在
许多实际问题中是以列表函数的形式给出, 即仅仅知道其在一些节点处的函数值， 牛顿-莱布尼茨公式也不能直接运用，因此有必要研究积分的数值计算问题，数值积分是解决上述困难的一种有效方法.
§1.1基本思想及基本概念
§1.1.1 方法思想
由定积分中值定理：
b a a b f dx x f I b
a ≤≤-==?ξξ),
)(()(
可知: 积分可以通过被积函数在ξ处的值得到. 由于积分中值定理仅仅告诉我们ξ在一定条件下是存在的, 但并没有给出确定ξ的方法. 一个很自然的想法就是利用被积函数)(x f 在节点b x x x x a n =≤≤≤= 210处函数值的加权平均来替代(近似))(ξf , 按此思想有
)()(0
i n
i i b
a
x f A dx x f ∑?=≈ （1-1）
这就是数值求积的思想（有效地解决了本章开始提出的问题），权因子i A 和节点i x
n i ,,2,1,0 =的不同确定方法就对应不同的数值求积公式.
§1.1.2 几个基本概念
定义1.1 称形如(1-1)式的求积公式为机械求积公式,其中i A 仅节点的选择与
)(x f 无关，b x x x x a n =≤≤≤= 210称为求积节点，i A （n i ,,2,1,0 =）称为求积
系数.
定义1.2 如果某个求积公式对于次数不超过m 的多项式均能准确地成立,而对于
1+m 次多项式就不准确成立, 则称该求积公式具有m 次代数精度（或代数精确度）.
注1.1 a) m 越大近似程度越高，标志着使函数准确成立的“个数”越多，但代数精度不是唯一衡量标准.
b) 若机械求积公式的代数精度0≥m ，则有a b A n
i i -=∑=0.
c) 若机械求积公式的代数精度为m ，即当m x x x f ,,1)(=时，由（1.1）
式可得，对任意次数不超过m 的k 次多项式m k x P k ≤),(有)()(0
i k n
i i b
a
k x P A dx x P ∑?=≡.
d) 代精度的高低, 从侧面反映求积公式的精度高低.
定义1.3 称求积公式
∑==n
k k k n x f A I 0)(
为插值型求积公式，式中求积系数k A 通过插值基函数
.,1,0)
())(()()())(()()(110110n k x x x x x x x x x x x x x x x x x l n k k k k k k n k k k =--------=
+-+-
积分求得，即
.,,1,0,
)(n k dx x l A b
a k k ==? （1-2）
定理1.1 插值型求积公式的代数精度至少为n 次.
定义1.4 若节点将被积区间等分成n 等分, 即.,2,1,0,n i i n
a
b a x i =-+
=则相应的插值求积公式称为Newton-Cotes (牛顿-柯特斯)求积公式. 即等距节点情形下的插值求积公式称为牛顿-柯特斯公式, 相应的求积系数称为Cotes 系数.
常见的几个简单求积公式( Newton-Cotes 公式)，如表1-1所示：
表1-1 几种简单N-C 求积公式总结表
n
名称
形式
1=n
梯形求积
公式 )]()([2
)(b f a f a
b T dx x f b
a
+-=
≈?
2=n
辛普森求
积公式
)]()2
(4)([6)(b f b
a f a f a
b S dx x f b
a
+++-=
≈?
4=n
柯特斯求
积公式
)](7)(32)(12)(32)(7[90
)(321b f x f x f x f a f a
b C dx x f b
a
++++-=
≈?
其中.1,,1,,-=-=
+=n k n
a
b h kh a x k 注1.2 a ）8≥n 时，N-C 公式出现数值不稳定.
b ）n 为偶数时，N-C 公式的代数精度至少为1+n 次，n 为奇数时，N-C 公式的代数精度至少为n 次.
定义1.5 截断误差: 由
（1-3）
当1=n 时可得梯形求积公式的截断误差T R
]
,[)(,],[,)(12)("))((2
)("))((!
2)
("23b a C x f b a a b f dx b x a x f dx
b x a x f T I R b a b a T ∈∈--=--=--=-=??ηηηξ 类似的，可得当2=n ，4=n 时的截断误差
注1.3 从截断误差公式可知，当区间长度a b -较大时，求积公式误差较大.
§1.2辛普森算法基本定义及其应用
§1.2.1 辛普森求积公式的定义
设计积分区间],[b a 划分为n 等份，步长n
a
b -，选取等距节点kh a x k +=构造出的插值型求积公式
)()()
(k n k n x f c a b I -=
为牛顿—柯特斯（Newton-Cotes ）公式，式中)
(n k c 称为柯特斯系数.
根据插值型求积公式系数（1-2），引进变换th a x +=，则有
?∏?∏≠=-≠=---=---=n k
j j k n n n k
j j n k dt j t n k n nk dt j k j t a b h c 00
00)
()()!(!)1( 当2=n 时，由上式有
6
1
)2)(1(4120)
2(0
=--=?dt t t c dx x n f x f A dx x f f R f I f I b
a
n n i n
i i b
a
n n )()!
1()
()()(][][][1)1(0
?
∑?++=+=-==-ωξ
64)1(2120)2(1
=--=?dt t t c
6
1)1(4120)2(2
=-=?dt t t c 则相应的求积公式是辛普森求积公式：
)]()2
()([6)(b f b a f a f a b dx x f s b a ++-==? （1-4）
§1.2.2辛普森求积公式的几何意义
辛普森公式的几何意义就是用通过A,B,C 三点的抛物线)(x L y =代替)(x f y =所得曲边梯形面积，如图1.1所示.
§1.2.3辛普森求积公式的代数精度及其余项
由N-C 公式的特点知，当n 为偶数时N-C 公式的代数精度至少为1+n 次，由于Simpson 求积公式为2=n 时的N-C 公式，故它的代数精度至少为3次，
即3≥m
将4)(x x f =代入Simpson 公式（1-4）
左边5
5
54
a b dx x b
a
-==?
右边≠+++-=
))2
(4(644
4b b a a a b 左边 由此可知4)(x x f =使得Simpson 求积公式不准确成立，所以3=m 即Simpson 公式代数精度为3次
由N-C 公式的余项公式（1-3）知，当2=n 时可得辛普森求积公式的截断误差
y
x
O
)
(x L y =a
2
b
a +
b A
B
C
)
(x f y =图1.1 辛普森求积公式的几何意义图
],[)(],,[),()2
(1804)
4(4b a C x f b a f a b a b R s ∈∈---
=ηη （1-5） §1.2.4辛普森公式的应用
例1.1 用辛普森求积公式计算积分
dx x x
?+1
024.
由积分形式可知 2,1,0===n b a 用辛普森公式计算有下式
)]1()2
1(4)0([6141
02f f f dx x x s ++=+=?
其中2
4)(x x
x f +=
. 计算流程图
C 语言程序代码及其运算结果详见附录A 分析附录A 可知
111765.041
02=+?dx x x
开始
定义函数f （x ）
输入n ，a ，b 的值
计算h=（b-a ）/n
调用函数f （x ），计算s 的值
输出s 的值
结束
图1.2 例1.1流程图
第二章 辛普森求积公式的拓展及其应用
为了提高精度，通常在实际应用中往往采用将积分区间划分成若干个小区间，在各小区间上采用低次的求积公式，如：梯形公式或辛普森公式，然后再利用积分的可加性，把各区间上的积分加起来，便得到新的求积公式，这就是复化求积公式，本章重点介绍复化辛普森求积公式.
§2.1 复化辛普森求积公式
§2.1.1问题的提出
由截断误差可知，当区间长度a b -较大时，Newton-Cotes 求积公式的误差较大. 为构造更高精度的数值积分公式，可以采用分段低次多项式替代整体高次多项式，为此，利用积分关于区间具有可加性，将],[b a 区间上的积分，分成若干小区间上的积分，以此来减少积分区间长度引起的误差.这就引用了复化求积公式. 其基本思想是：先把积分区间分成一些长度较小的子区间，在每个子区间上使用低阶的牛顿-柯特斯公式，即利用
n
a
b h ih a x dx x f dx x f i n
i x x b
a
i
i -=
+==∑?
?
=-,,)()(1
1
并把小区间n i x x i i ,,2,1],,[1 =-上的积分dx x f i i x x ?
-1
)(用前面的方法近似求得，
由此即可得到相应的复化求积公式. 最常用的是复化梯形公式和复化辛普森公式，下面学习辛
普森求积公式.
§2.1.2复化辛普森公式及其分析
定义 2.1 将小区间n i x x i i ,,2,1],[1 =-上的积分分别用辛普森公式计算，即可得到复化辛普森公式
n
n
i i n i i i i i b
a
n
i S b f x f x f a f h
x f x f x f h
dx x f =+++=++≈∑∑?
∑=--=--=)]()(4)(2)([6)]
()(4)([6
)(1
1111212
1
其中2
2
1h x x i i =
=-. 另一种定义形式为：用分段二次插值函数代替，记1,2,1,0,2-==m k m n 在第k 段的两个小区间上，用三个结点))(,()),(,()),(,(2222121222++++k k k k k k x f x x f x x f x 作二次插值函数)(x s k ，然后积分，求m 段之和可得整个区间上的近似积分
m
a
b h x f x f x f x f h
s m k k m k k m n 2)
)(2)(4)()((31
1
2101220-=
+++=∑∑-=-=+ 称该求积公式为复化辛普森求积公式（抛物线公式）.
定理2.1 若],,[)(4b a C x f ∈则复化辛普森公式的截断误差为
b a f h a b S dx x f n b a
≤≤--=-?ξξ),()2
(180)()()4(4 且
0)],()([)2
1(1801)("'"'
)4(4
→-→
-?
h b f a f h S dx x f b
a
n
. 注 2.1 从误差公式可以看出当],[)(4b a C x f ∈时，n S 比n T 2的精度一般要高，但他们的计算量几乎一样.
注2.2 ○1n
S 属于机械型求积公式，但不属于插值型、也不属于N-C 求积公式. ○2n S 的代数精度为4次，具有稳定性和收敛性即][f I S n
→（∞→n 或∞→h ）.
§2.1.3复化辛普森公式计算流程图
为了减少计算工作量，优化程序设计，将复化辛普森公式
n
n
i i n i i i i i b a n i S b f x f x f a f h
x f x f x f h
dx x f =+++=++≈∑∑?∑=--=--=)]()(4)(2)([6)]()(4)([6)(1
1111
2121
改写为
])}2(])12([2{)]()([5.0[3}])12([2)2()]()([5.0{3}
])12([4)2(2)()({61
111111
∑∑∑∑∑==-=-=-++-++--=-+++++-=-+++++-=n i n i n i n i n
i n ih a f h i a f b f a f n a b h i a f ih a f b f a f n a b h i a f ih a f b f a f a b S 则于此相对应的辛普森流程图为：
§2.1.4 复化辛普森公式的应用
例2.1 用复化辛普森公式计算正弦积分的近似值
.
分析该积分可知,sin )(dx x
x
x f =
0=a ，1=b 则 125.08
1==-=n a b h 为步长
C 语言程序代码及其运算结果详见附录B 由此可知
94608.04=S
开始
输入A,B,N
H=（B-A ）/（2*N ）
S=0.5*（F （A ）- F （B ）），调用函数F
S=S+2*F[A+（2*I-1*H ）]+（F （A+2*I*H ）），调用函数F
S=（B-A ）/（3*N ）S
输出S
结束
I=1,N
定义函数F
++I
8sin 1
0==
?n dx
x x
S n 图2.1 复化辛普森算法流程图
例2.2 用复化辛普森公式计算定积分
84102=+?n dx x x
. 分析该积分可知2
4)(x x x f +=，0=a ，1=b 则125.081
==-=n a b h 为步长 C 语言程序代码及其运算结果详见附录B. 由此可知
11157.04=S
在利用插值型求积公式求积分时，为了提高精度有两种途径.一是提高积分区间上的插值多项式的阶数，从而也就提高了求积的阶数.但是，由于插值多项式的阶数越高，其逼近性质未必好（即精度未必能提高），因此，牛顿-柯特斯公式的阶数越高，其积分精度也未必越高，工程上一般只作到六阶牛顿-柯特斯公式（即龙贝格公式）为止.二是采用复化公式，尽量减小每个求积小区间的长度.在实际应用时，往往将两种方法混合使用，以便提高求积的精度.
§2.2 变步长辛普森求积公式
在数值积分中，精度是一个很重要的问题，如果误差太大，就没有实际意义.为了提高精度，通过需要在复化求积公式中尽量减少各细分小区间的长度，即减少步长
h .显然，如果步长h 取得太大，则精度就难以保证.但是，如果步长取得太小，则计算工作量就随之增大，并且，由于项数增加，其误差积累也就增大.因此，在采用复化公式求积时，关键的问题是合理地选择步长（即合理选择对整个积分区间的细分数），以便既能满足精度要求，又不至于引起过多的误差积累和过大的计算工作量.在实际计算过程中，通常采用变步长的求积法.
§2.2.1变步长辛普森求积公式的导出过程
变步长辛普森求积公式是建立在变步长梯形公式的基础上，同时它又是龙贝格算法导出的中间过程，我们知道, 若被积函数具有一定的光滑性, 则增加节点可以降低复化求积公式的截断误差.
这里需要解决的问题是增加节点后的复化求积方法能否充分利用已有的计算工作量. 譬如: 若将n T 作为?=b
a dx x f I )(的近似精度不够, 需减少步长(增加节点数)
计算相应的m T 来近似I , 当然我们想要充分利用已经求得的n T .为此, 设区间
n b a ],[等分后, 利用复化梯形公式已经求得n T 这一结果, 为了得到精度更高的数值
结果, 我们将原有的步长折半, 即把区间],[b a 分为n 2等分, 然后应用复化梯形公
式求得n T 2.下面将会看到这样既提高了精度, 又能充分利用已经求得的n T .事实上, 我们可以建立n T 与n T 2的下述递推关系. 设
n
a
b h x f x f h T n i i i n -=
+=∑-=+,
)]()([21
1 则
∑∑∑∑-=+-=+-=+++-=+=++=++=1
1
010111
02)
(221)(2)]()([4)]
()(2)([4
21212
1n i i n n k k n k k k i k k n h n x f h T x f h x f x f h x f x f x f h
T
其中n
a
b x x h k k -=-=-1 即，
∑-=+=1
2221n i n n
h T T 新增分点的函数值 注2.3 由上述公式可知在n T 的基础上计算n T 2只需调用n 次函数即可，最大限度地节省了n T 2的计算量.
加速公式的导出：由前面的误差分析，我们可以得到复化梯形公式n T 的截断误差为2)("12
h f a
b ξ--
，即 2)("12h f a
b T I n ξ--
=- 类似根据复化梯形公式n T 2的截断误差为2)2)(("12h
f a b η--，有 22)2
)(("12h
f a b T I n η-≈-
两式相比可得
4
1
2≈
-n T I ， 其中dx x f f I I b a ?==)()(
即
)(3
1
22n n n T T T I -≈- （2-1）
注2.4 ○
1公式（2-1）说明n T 2的误差可以近似地由n T 2与n T 表现， 这样就给出了
复化梯形公式估计误差的事后估计法.
○
2由公式（2-1）还可以得到校正公式（加速公式） n n n n n T T T T T I 3
1
34)(31222-=-+≈
数值实验结果表明，在一定条件下，上式计算出来的值比原来的n T 2好得多，上述公式称为梯形公式的加速公式.
梯形求积公式的实质：假设已知n T ，n T 2，则
n
k k k n k k k n k n k k k k k n n S x f x f x f h
x f x f h
x f x f h
x f x f h T T =++=+-
+++=-=+-=+-=-=+++∑∑∑)]()(4)([6
)]()([2
31)]}()([4
)]()([4{34313411
011010122
12121
即
n n n T T S 3
1342-=
上式表明n T 与n T 2通过上面公式处理后，可得精度更高的n S .即复化辛普森公式，这也是加速的实质.
§2.2.2变步长辛普森求积公式的加速过程
类似梯形加速公式的推导，由n S 的截断误差公式（1-5）可得
][151
2n n n S S S I -≈-
即
n n n n n S S S S S I 15
11516][15
1
222-=-+
≈
注2.5 ○
1上述两个公式分别称为复化辛普森公式估计误差的事后估计公式及复化辛普森公式的加速公式.
○
2类似地可以证明： n n n S S C 15
1
15162-=
○3在求得n C ，n
C 2的基础上，可以进一步加速得：龙贝格公式
n n n C C R 63
1
63642-=
§2.2.3变步长辛普森求积公式的算法流程图
开始 N=1，H=B-A
IP=F(A)+F(B) F
IC=0,X=A-H/2
K=1,N
X=X+H
IC=IC+F(x) F
I2=(4*IC+IP)*H/6
N=1
|I2-I| I1<=I2,IP=IP+2*IC N=N+N H=0.5*H Y N N I=2I 输出 结束 图2.2 变步长辛普森算法流程图 §2.2.4变步长辛普森公式算法程序代码 详见附录C §2.2.5变步长辛普森求积公式的应用 例2.3 用变步长辛普森求积公式计算定积分 dx x x ?+1024 取000001.0=ε. C 语言程序代码及其运算结果详见附录C. 分析结果可知 111572.041 02=+?dx x x §2.2.6小结 通过分析例1.1、2.2、2.3有下表2-1 表2-1 三种算法比较 算法名称 代数精度 积分形式 计算结果 余项 辛普森求积 3 dx x x ?+1 024 0.111765 111765.0)4ln 5(ln 2 1 -- 复化辛普森求积 4 dx x x ?+1 024 0.11157 11157.0)4ln 5(ln 2 1 -- 变步长辛普森求积 dx x x ?+1 024 0.111572 111572.0)4ln 5(ln 2 1 -- 由表2-1可以得出用变步长辛普森求积公式求得的结果偏离准确值的程度最小，即其计算结果最接近准确值，其次是复化辛普森求积方法，辛普森求积方法较前述两种方法误差较大.但三种算法均具有良好的稳定性与收敛性，从而提高了计算效率及准确度．在工程技术中有较为广泛的应用. §2.2.7 数值求积公式在实际工程中的应用 例 2.4人造地球卫星轨道可视为平面上的椭圆。我国第一颗人造地球卫星近地点距地球表面439Km ，远地点距地球表面2384Km ，地球半径为6371Km. 求：该问题的轨道长度. 模型 1）a ，b 分别是长半轴和短半轴； 2）焦距为c ； 3）地球半径为r ； 近地点和远地点与地球表面的距离分别是1S 和2S . r s s a 2221++= r s a c --=1 22c a b -= )(5.7782Km a = )(5.7721Km b = 椭圆的参数方程为 t a x cos =，)20(sin π≤≤=t t b y 弧长的公式?椭圆长度 dt t b t a L 21 2)cos sin (40 2222?+=π 利用数值积分公式计算卫星轨道长度问题，编写Mathematica 语句如下： }]2 ,0,{,][5.7721][5.7782*4[2222π t t Cos t Sin NIntegrate + 结果为:48707.4 图2.3 卫星轨道的示意图 参考文献 [1]李庆扬王能超易大义.数值分析基础[M].北京:清华大学出版,2008年. [2]徐士良.数值方法与计算机实现[M].北京:清华大学出版社，2006年. [3]封建湖聂玉峰王振海. 数值分析[M].西安:西北工业大学出版社，2003年. [4]朝伦巴根贾德彬. 数值计算方法[M].北京:中国水利水电出版社，2007年. [5]王立秋魏焕彩周学圣.工程数值分析题解[M].山东：山东大学出版社，2004年. [6]张韶华奚梅成陈效群 .数值计算方法与算法.北京：科学出版社，2006年. 附录A 本附录介绍例1.1 用辛普森求积公式计算过程中的变量说明、C语言程序代码、以及运行结果. 表变量说明表 变量名变量说明 a 积分上限 b 积分下限 n 区间等分后的个数 h 步长 s 结果 该积分用Simpson公式计算的C语言程序代码为： #include #include float f (float x) { float t; t=x/(4+pow(x,2)); return t;} main() {float n,a,b; float h,s; printf("putn="); scanf("%f",&n); printf("puta="); scanf("%f",&a); printf("putb="); scanf("%f",&b); printf("a=%f￡?b=%f,n=%f\n",a,b,n); h=(b-a)/n; printf("h=%f\n",h); s=(h/3)*(f(a)+4*f(a+h)+f(b)); printf("s=%f\n",s); } 其结果为： 图例1.1计算结果 实验四、复化梯形公式和复化Simpson公式的精度比较 （2学时） 一、实验目的与要求 1、熟悉复化Simpson公式和复化梯形公式的构造原理； 2、熟悉并掌握二者的余项表达式； 3、分别求出准确值，复化梯形的近似值，复化Simpson的近似值，并比较后两 者的精度； 4、从余项表达式，即误差曲线，来观察二者的精度，看哪个更接近于准确值。 二、实验内容： 对于函数 sin () x f x x =，试利用下表计算积分1 sin x I dx x =?。 表格如下： 注：分别利用复化梯形公式和复化Simpson公式计算，比较哪个精度更好。其中：积分的准确值0.9460831 I=。 三、实验步骤 1、熟悉理论知识，并编写相应的程序； 2、上机操作，从误差图形上观察误差，并与准确值相比较，看哪个精度更好； 3、得出结论，并整理实验报告。 四、实验注意事项 1、复化梯形公式，程序主体部分： for n=2:10 T(n)=0.5*T(n-1) for i=1:2^(n-2) T(n)=T(n)+(sin((2*i-1)/2^(n-1))/((2*i-1)/2^(n-1)))/2^(n-1)； end end 2、复化Simpson公式，程序主体部分： for i=1:10 n=2.^i x=0:1/n:1 f=sin(x)./x f(1)=1 s=0 for j=1:n/2 s=s+f(2*j) end t=0 for j=1:(n/2-1) t=t+f(2*j-1) end S(i)=1/3/n*(f(1)+4*s+2*t+f(n+1)) end 五．实验内容 复化梯形公式和复化辛普森公式的引入 复化梯形公式： 1 10[(()]2 n n k k k h T f x f x -+==+∑； 复化辛普森公式： 1 1102 [(4()()]6n n k k k k h S f x f x f x -++ ==++∑； 根据题意和复化梯形公式、复化辛普森公式的原理编辑程序求解代码如下： Matlab 代码 clc s=quad('sin(x)./x',0,1) p1=zeros(10,1); 数值分析实验报告—— 实验目的[1] 掌握复化梯形和辛普森数值积分法的基本原理和方法； [2] 编程MA TLAB程序实现复化梯形和辛普森数值积分 实验内容与步骤1.编程序实现复化梯形数值积分求积公式 function y=f(x) y=sqrt(x).*log(x); function T_n=F_H_T(a,b,n) h=(b-a)/n; for k=0:n x(k+1)=a+k*h; if x(k+1)==0 x(k+1)=10^(-10); end end T_1=h/2*(f(x(1))+f(x(n+1))); for i=2:n F(i)=h*f(x(i)); end T_2=sum(F); T_n=T_1+T_2; 实验内容与步骤运行结果： >> T_n=F_H_T(0,1,20) T_n = -0.4336 2.编程序实现复化辛普森数值积分求积公式 function y=f(x) y=sqrt(x).*log(x); function S_n=S_P_S(a,b,n) h=(b-a)/n; for k=0:n x(k+1)=a+k*h; x_k(k+1)=x(k+1)+1/2*h; if (x(k+1)==0)|(x_k(k+1)==0) x(k+1)=10^(-10); x_k(k+1)=10^(-10); end S_1=h/6*(f(x(1))+f(x(n+1))); for i=2:n F_1(i)=h/3*f(x(i)); end for j=1:n F_2(j)=2*h/3*f(x_k(j)); end S_2=sum(F_1)+sum(F_2); S_n=S_1+S_2; 运行结果： >> S_n=S_P_S(0,1,20) S_n = -0.4423 实验心得 通过此次实验的操作，我掌握了复合梯形公式和复合辛普森公式，对编程又有了新的突破！ 第四章共包含9个求积公式，1个余项公式。 1，机械求积公式 f x dx = A k f (x k )n k =0b a 2，插值求积公式 Ln x dx =b a [ l k (x )dx b a L (x k )n k =0] 3，梯形求积公式 f x dx = b ?a b a [f a +f b ] R n x =? b ?a 3f ′′ ξ 4，辛普森求积公式 f x dx = b ?a b a [f a +f (a +b )+f b ] R n x =? b ?a (b ?a )4f (4) ξ 5，复合梯形公式 f x dx =?b a [f a + f x k n?1k =1+f b ] h=(b-a)/n R n x =? b ?a h 2f ′′ ξ 6，复合辛普森公式 f x dx =?b a [f a +4 f x k +12 n?1k =0+2 f x k n?1k =1+f b ] h=(b-a)/n R n x =? b ?a (h )4f (4) ξ 7，高斯求积公式 ρ(x )f x dx = A k f (x k )n k =0b a 其中x k 为高斯点，n+1个节点对应2n+1级代数精度。 高斯点公式：ωn+1=（x-x 0）（x-x 1）…（x-x n ）= x n+1 + a 0x n + a 1x n-1+…+a n-1x+a n ，用 ρ(x )ωn +1 x φk (x )dx b a =0(k=0,…，n)求出待定系数a ，解方程ωn+1=0得高斯点。 重新代入 ρ(x )f x dx = A k f (x k )n k =0 b a 中求解方程组得到系数A 。 摘要 在工程实验及研究中，实际工作中，变量间未必都有线性关系，如服药后血药浓度与时间的关系；疾病疗效与疗程长短的关系；毒物剂量与致死率的关系等常呈曲线关系。曲线拟合是指选择适当的曲线类型来拟合观测数据，并用拟合的曲线方程分析两变量间的关系.可以说，曲线拟合模型与我们的生活生产密切相关. 本课题着重介绍曲线拟合模型及其应用，其中包括它的基本思想、模型的建立、以及具体应用.为了更好的了解曲线拟合模型，可以将它分为线性与非线性模型，在模型建立的基础上我们可以用最小二乘法来解决一些我们日常所应用的问题. 关键词曲线拟合；线性与非线性模型；最小二乘发 目录 引言 (1) 第一章曲线拟合 (2) §1.1 基本思想及基本概念 (2) §1.1.1 方法思想 (2) §1.1.2几个基本概念 (2) §1.2辛普森算法基本定义及其应用 (4) §1.2.1辛普森求积公式的定义 (4) §1.2.2辛普森求积公式的几何意义 (5) §1.2.3辛普森求积公式的代数精度及其余项 (5) §1.2.4辛普森公式的应用 (6) 第二章辛普森求积公式的拓展及其应用 (7) §2.1 复化辛普森求积公式 (7) §2.1.1问题的提出 (7) §2.1.2复化辛普森公式及其分析 (7) §2.1.3复化辛普森公式计算流程图 (8) §2.1.4复化辛普森公式的应用 (9) §2.2 变步长辛普森求积公式 (10) §2.2.1变步长辛普森求积公式的导出过程 (10) §2.2.2变步长辛普森求积公式的加速过程 (12) §2.2.3变步长辛普森求积公式的算法流程图 (13) §2.2.4变步长辛普森公式算法程序代码 (14) §2.2.5变步长辛普森求积公式的应用 (14) §2.2.6小结 (14) §2.2.7数值求积公式在实际工程中的应用 (14) 参考文献 (16) 附录A (17) 数值分析与实验课程设计 班级: 姓名: 学号: 08级应用数学《数值分析与实验（实践）》任务书 一、设计目的 通过《数值分析与实验（实践）》实践环节，掌握本门课程的众多数值解法和原理，并通过编写C语言或matlab程序，掌握各种基本算法在计算机中的具体表达方法，并逐一了解它们的优劣、稳定性以及收敛性。在熟练掌握C语言或matlab语言编程的基础上，编写算法和稳定性均佳、通用性强、可读性好，输入输出方便的程序，以解决实际中的一些科学计算问题。 二、设计教学内容 1、数值方法的稳定性； 2、禾U用牛顿法和割线法程序求出非线性方程的解，并比较它们之间的优 劣； 3、高斯消去法和列主元高斯消去法求解线性方程组； 雅克比法和高斯-赛德尔迭代法解方程组； 4、利用Lagrange插值多项式求未知点的近似值； 5、利用所给数据进行数据的多项式和可转化成多项式形式的函数拟合； 6、编写复化辛卜生公式和龙贝格算法，通过实际计算体会各种方法的精确 \ 度； 7、利用改进Euler方法和四阶Runge-Kutta方法求解初值问题的微分方程组； &利用幕法求矩阵按模最大的特征值及对应特征向量； \ （8个中选取1个） 二、设计时间 2011 —2012学年第1学期：第16周共计一周 教师签名: 2011年12月12日 、八 刖 数值计算方法是一种利用计算机解决数学 .言 问题的数值近似解方法， 特别是无法用人工过计算器计算的数学问题。数值计算方法常用于矩阵高次代数方程矩阵特征值与特征向量的数值解法，插值法，线性方程组迭代法，函数逼近，数值积分与微分，常微分方程初值问题数值解等。 作为数学与计算机之间的一条通道，数值计算的应用范围已十分广泛，作为用计算机解决实际问题的纽带，数值算法在求解线性方程组，曲线拟合、数值积分、数值微分，迭代方法、插值法、拟合法、最小二乘法等应用广泛。 数值计算方法是和计算机紧密相连的，现代计算机的出现为大规模的数值计 算创造了条件，集中而系统的研究适用于计算机的数值方法是十分必要的。数值计算方法是在数值计算实践和理论分析的基础上发展起来的。 通过数值计算方法与实验将有助于我们理解和掌握数值计算方法基本理论和相关软件的掌握，熟练求解一些数学模和运算，并提高我们的编程能力来解决实际问题。 关于复化辛普森(simpson)公式在线路坐标计算中的应用 天津西站项目部刘思传 摘要：本文里利用辛普森公式导证了线路坐标计算的公式，并在卡西欧FX-4800P计算器中编写了中边线坐标计算的源程序。 关键词：复化辛普森公式，线路坐标计算，曲率。 一.引言 随着我国道路建设等级和质量水平的飞速发展，公路、铁路建设的机械化和日产量日益提高，促使施工中在满足设计精度的前提下，尽可能快速、准确地进行测量放样和检查工作，本文线路曲率变化的特点，利用复化辛普森公式导证了线路坐标计算的通用公式，并利用卡西欧FX-4800P计算器编写了计算线路中边线坐标的源程序。 二.复化辛普森公式数学模型 把积分区间分成偶数等分，记，其中是节点总数，是积分子区间的总数。 记，，在每个区间上用辛普森数值积分公式计算，则得到复化辛普森公式，记为。 复化辛普森积分计算公式 而，称 (1) 式(1)即为辛普森复化公式。 三.线路坐标计算 2. 回旋曲线上点位坐标方位角的计算 如图1，设回旋曲线起点A 的曲率为A ρ，其里程为DK A ；回旋曲线终点B 的曲率为B ρ，其里程为DK B ，Ax ’'y 为以A 为坐标原点，以A 点切线为'x 轴的局部坐标系；Axy 为线路坐标系。 由此回旋曲线上各点曲率半径为R i 和该点离曲线起点的距离?i 成反比，故此任意点的曲率为 c l R i i i /1==ρ（=为常数）. （2） y ' Y B 图1 由式（2）可知，回旋曲线任意点的曲率按线性变化，由此回旋曲线上里程为DK i 点的曲率为 )(A i A B A B A i DK DK DK DK ---+=ρρρρ （3） 当曲线右偏时，取正；当曲线左偏时取负。在图1中有 ???????=== ?I A DK DK i i i dl dl dl R d ρβρβ1 （4） 将式（3）代入式（4）得 πρρβ180 *)(2A i A i i DK DK -+= （5） 若已知回旋曲线起点A 在线路坐标系下切线坐标方位角αA ，则里程为Dk i 点切线坐标方位角为 i A i βαα+= π180 （6） 将式(5)代入式（6）得 *)(2A i A i A i DK DK -++=ρραα π180 （7） 对于式(7) ，当，时，，则a i =a A ，式(7)变成计算直线段上任意点切线坐标方位角计算公式；当，时，， ，则式(7)代表圆曲线上任意点切线坐标方位角 计算公式。 可见，若已知曲线段起点和终点的曲率及起点的切线坐标方位角，式(7)便能计算任意线型点位切线坐标方位角。 3、回旋曲线点位坐标计算 由图1可得回旋曲线上点位在坐标系下坐标计算公式： SHANGHAI JIAO TONG UNIVERSITY 题目名称：复合梯形公式与复合辛普森公式对比学生: 学生学号: 班级: 学院(系): 目录 1.概述 (3) 2.问题提出 (3) 3.算法推导 (3) 4.算法框图 (3) 4.1复合梯形公式算法流程图 (3) 4.2 复合辛普森公式算法流程图 (3) 5.MATLAB源程序 (3) 6.结论与展望 (3) 图表目录 图4-1 复合梯形公式算法流程图 (1) 图4-2 复合辛普森公式算法流程图 (1) 图6-1 MATLAB计算结果 (3) 表2-1函数计算结果表 (3) 1.概述 梯形求积公式和辛普森求积公式分别是牛顿-科斯特公式中n=1和n=2时的情形。其中梯形求积公式可表示为 由于牛顿-科斯特公式在n≥8时不具有稳定性，故不可能通过提高阶的方法来提高求积精度。为了提高精度通常可把积分区间分成若干子区间（通常是等分），再在每个子区间上用低阶求积公式。这种方法称为复合求积法。 本文主要讨论复合梯形公式和复合辛普森公式在同一数学问题中的应用。首先给出了复合梯形公式和复合辛普森公式的推导过程以及其余项的表达形式，然 后用流程图的形式介绍算法思路，再运用MATLAB编写代码计算结果，最后对结果进行对比讨论。 希望通过两个算法在同一个算例中的应用对比，更好的理解和掌握复合梯形公式和复合辛普森公式的适用围和适用条件。并且能够熟悉MATLAB编程求解问题的流程，掌握编程化的思想方法。同时对两种方法的计算结果对比分析，讨论两种求积方法的计算精度。 2.问题提出 对于函数给出的函数表如下，试用复合梯形公式和复合辛普森公 式计算积分。 表2-1函数计算结果表 x f(x) 0 1 1/8 0.997397867081822 1/4 0.989615837018092 期中测试（下） 班级： 姓名： 学号： 分数： 一、填空题（20分） 1、计算积分1 ? 2、5个节点的牛顿- 3、求积公式 0()()b n k k k a f x dx A f x =≈∑? 4 、数值积分公式 ()1 1 ()29[(1)8(0)(1)]f x dx f f f -'≈-++? 5、求解一阶常微分方程初值问题00(,),()y f x y y x y '==的改进欧拉公式为 二、选择题（6分） 1、舍入误差是( A )产生的误差。 A. 只取有限位数 B ．模型准确值与用数值方法求得的准确值 C ． 观察与测量 D ．数学模型准确值与实际值 2、用 1+x 近似表示e x 所产生的误差是( C )误差。 A ． 模型 B ． 观测 C ． 截断 D ． 舍入 3、解线性方程组的主元素消去法中选择主元的目的是( A )。 A ．控制舍入误差 B ． 减小方法误差 C ．防止计算时溢出 D ． 简化计算 4、求解初值问题00(,),()y f x y y x y '==欧拉法的局部截断误差是( A )；中心欧拉法的局部截断误差是( B )； 改进欧拉法的局部截断误差是( B )；四阶龙格－库塔法的局部截断误差是( D ) A. 2()O h B. 3()O h C. 4()O h D. 5()O h 三、计算题（64分） 1、（10分）试分别推导复化梯形和复化辛普森求积公式。 证明：以积分 ()b f x dx 为例。将积分区间[,]a b 做n 等分，步长()/h b a n =-。 2、(10分) 求A 、B 使求积公式1 1 ()[(1)(1)][(0.5)(0.5)]f x dx A f f B f f -≈-++-+? 的代数精度尽量高, 并求其代数精度； 利用此公式求2 1 1 I dx x = ? (保留4位小数)。 解：2 ,,1)(x x x f =是精确成立，即 3、(12分) 取5个等距节点 ，分别用复化梯形公式和复化辛普森公式计算积分2 201 I dx = ?的近似值（保留4位小数）。 高等数值分析拉格朗日插值多项式切比雪夫高斯龙格现象复合梯形辛普森求积公式 解答： 1.拉格朗日插值函数： function y=lagrange (a,b,x) y=0; if length(a)==length(b) n=length(a); else disp('ERROR!length(a)!=length(b)') return; end for i=1:n k=1; for j=1:n if j~=i k=k.*(x-a(j))/(a(i)-a(j)); end end y=y+k*b(i); end 2.问题（a）： function Q_a m=100; n=10; x=-1:2/n:1; y=1./(1+9*x.^2); x0=-1:2/m:1; y0=lagrange(x,y,x0); y1=1./(1+9*x0.^2); plot(x0,y0,'--r'); hold on; plot(x0,y1,'-b'); end 3.问题（b）： function Q_b m=100; n=10; x=zeros(1,n+1); for i=1:n+1 x(i)=cos((2*i-1)*pi/(2*n+2)); end y=1./(1+9*x.^2); x0=-1:2/m:1; y0=lagrange(x,y,x0); y1=1./(1+9*x0.^2); plot(x0,y0,'--r'); hold on; plot(x0,y1,'-b'); end 4.问题（c）： main.m（m文件） figure(1) Q_a() figure(2) Q_b() syms x y=1/(1+9*x^2); I0=int(y,-1,1);%准确值 n=10; 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. 复化梯形公式与辛普森公式的matlab程序【程序代码】 cclc; disp('1.复化梯形公式求解'); disp('2.simpson公式求解'); disp('请进行选择：'); c=input(' '); if c==1 clc; disp('复化梯形公式'); disp('请输入积分下限'); a=input('a='); disp('请输入积分上限'); b=input('b='); disp('请输入等分的数目'); n=input('n='); h=(b-a)/n; s1=0; for i=1:n-1 s1=s1+fun1(i*h); end disp('复化梯形公式的结果：'); T=h/2*(fun1(a)+2*s1+fun1(b)) else if c==2 clc; disp('simpson公式'); disp('请输入积分下限'); a=input('a='); disp('请输入积分上限'); b=input('b='); disp('请输入等分的数目'); n=input('n='); h=(b-a)/n; s2=0; for i=0:n-1 s2=s2+fun1((i+0.5)*h); end disp('辛普森公式的结果：'); S=h/6*(fun1(a)+4*s2+2*s1+fun1(b)) end end disp('菜单选项'); disp('1.继续运算'); disp('2.退出程序！'); p=input(' '); if p==1 (fuhua); else if p==2 disp('正在退出，请稍候。。。'); 第三次实验 实验名称复合辛普森求积计算积分实验时间2012.05.06 姓名班级数应二班学号成绩 一、实验目的，内容 1.目的： 学习并理解复合辛普森求积计算积分的matlab实现。 2. 内容： 用matlab计算积分 1 4 ln 9 x xdx=- ? （精确值-0.4444），并求出达到。 二．代码 1. function Sn=ComSimpson(a,b,n) %复合辛普森求积 %f表示被积函数，本题中由f.m调用 %a,b分别表示积分上下限 %n表示区间分割次数 %sn表示该方法计算所返回的值 h=(b-a)/n; for k=0:n x(k+1)=a+k*h; x_k(k+1)=x(k+1)+1/2*h; if (x(k+1)==0)|(x_k(k+1)==0) x(k+1)=10^(-10); x_k(k+1)=10^(-10)；%误差精度 end end S1=h/6*(f(x(1))+f(x(n+1)));%复合辛普森求积公式第一部分for i=2:n F_1(i)=h/3*f(x(i)); end for j=1:n F_2(j)=2*h/3*f(x_k(j)); end S2=sum(F_1)+sum(F_2);%复合辛普森求积公式第二部分 Sn=S1+S2;%算得值 f的表达式由f.m（见附）文件调用。 附：function y=f(x) y=sqrt(x).*log(x); 三.数值结果: 在命令窗口输入如下指令时，出现如下结果 Sn=ComSimpson(0,1,300) Sn =-0.4438 Sn=ComSimpson(0,1,700) Sn =-0.4442 四．计算结果的分析: 当步长取值很小的时候，误差较大。当步长取得越细，即区间分割的越小时，实验计算值的精度越高，即越 趋近精确值。 五. 计算中出现的问题，解决方法及体会： 本实验过程中，了解了复合求积公式的几个构成，以及在matlab中的实现，深化了对于该问题的理解。一的个很简单的问题，哪怕就是一个小小符号的不注意，也不会得到结果。实验的严谨性，细腻性有待进一步提升。 对于实验，网上的程序有不少，关键是要看懂，弄懂，在实际的问题中能运用。借鉴吸收的同时，学会编程。多练习，才会有效果。 教师评语 指导教师：年月日 复合梯形公式与复合辛普森公式求积分 2010-12-26 09:37:23| 分类：数值分析| 标签：|字号大中小订阅 一实验目的 1. 掌握复合梯形公式与复合辛普森公式的基本思想。 2. 编程实现用复合梯形公式与复合辛普森公式求积分。 3. 熟悉matlab软件的使用。 二实验内容 1、用复合梯形公式计算积分I=4/(1+x2)dx ,求它0到1的积分。精确度为10-5.（0.00001） ，精确到 ●1计算公式 h=(b-a)/n h=h/2[(f(x0)+f(x1))+(f(x1)+f(x2))+(f(x2)+f(x3)+...+(f(xn-1)+f(xn)] 。 l1 算法分析 En=h2/12[f'(b)-f'(a)] 将区间[a,b]等分成n个小区间，在小区间上分别应用低次积分公式来构造公式，通过for循环来实现， 分的越细，越接近实际结果，精确度越高。 l2 源程序 function f1=fun4(x) %原函数 f1=4/(1+x^2); function ff=fun2(x) %函数对x求导 ff=-8*x/((1+x^2)^2); function f=tixing(a,b) %a，b是区间 a=0;b=1; disp('******复合梯形公式******') h=0.008; %h表示区间被等分成若干份后，每两个数的间距 m=(a:h:b); %形成一维矩阵，每两个数间的间隔是h n=length(m); %求上矩阵的长度，即元素个数 for i=1:n-1 D(i)=fun4(m(i))+fun4(m(i+1)); end R=h/2*sum(D); %积分结果 E=-(h^2)*(fun2(b)-fun2(a))/12; %余项，即精度 t=pi-R; [R;E;t] 实验结果讨论和分析 通过对h的值的改变，发现h值越小，即等分的区间越小，结果越精确，精确度越高。通过手算得到积分结果为π，实验结果为3.14158198692313，结果正确，可见复合梯形公式的精确度较高，运算次数 实验5 复合辛普森公式 李涛 201226100108 计自1201 一、实验目的 ● 用复合辛普森公式计算积分 dx x ? +48 2cos 1， 使误差不超过-4 10（注意所给积分特点，做出相应的处理后再计算） 二、实验步骤 1.算法原理 复合辛普森原理： 将区间],[b a 划分为n 等分，在每个子区间[]1,+k k x x 上采用辛普森公式，若记 ,2 1 21h x x k k + =+则得 ● ∑? -=== 1 )()(n k b a dx x f dx x f I ● ).()]()(4)([61 21f R x f x f x f h n n k k k k +++=∑-=+ 记 ● ∑-=+++=1 21)]()(4)([6n k k k k n x f x f x f h S ● ],)()(2)(4)([6101 1 21∑∑-=-=++++=n k n k k k b f x f x f a f h 称为复合辛普森求积公式，其余项为 ● .),(),()2(180)(10 1) 4(4∑-=+∈-=-=n k k k k k n n x x f h h S I f R ηη 于是当],[)(4b a C x f ∈时，与复合梯形公式相似有 ● ),(),()2 (180)() 4(4b a f h a b S I f R n n ∈-- =-=ηη 易知误差阶为4 h ，收敛性是显然的，实际上，只要],[)(b a C x f ∈则可得到收敛性，即 ● ? = ∞ →b a n n dx x f S )(lim 此外，由于n S 中求积公系数均为正数，故知辛普森公式计算稳定。 2.算法步骤 复合辛普森： 首先将区间],[b a 划分为n 等分，在每个子区间[]1,+k k x x 上采用辛普森公式，若记 ,2 1 21h x x k k + =+则得 ∑-=+++=1021)]()(4)([6n k k k k n x f x f x f h S ])()(2)(4)([6101 1 21∑∑-=-=++++=n k n k k k b f x f x f a f h 算法过程： 这里将辛普森公式写为Sn()函数，然后在Solve()函数里依次计算S1,S2,S4,S6.......当相邻的精度小于eps 时退出循环，则S2n 保存结果。 三.程序代码 #include 2. 编写用变步长复化辛普森公式计算积分()b a f x dx ? 的程序。 用上面编写的程序计算下列积分并分析计算结果 （1）0cos xdx π ? （2）220cos x x dx （3）?10dx x 程序： function S=bianfuhuasimpson(fx,a,b,eps,M) % 变步长复合simpson 求积公式 % 调用方式： S=fuhuasimpson(@fx,a,b,epsilon) % fx -- 求积函数（函数文件） % a, b -- 求积区间 % eps -- 计算精度 % M--最大允许输出划分数 n=1; h=(b-a)/n; T1=h*(feval(fx,a)-feval(fx,b))/2; Hn=h*feval(fx,(a+b)/2); S1=(T1+2*Hn)/3; n=2*n; % 最好与倒数第三行保持一致（变步长） while n<=M T2=(T1+Hn)/2; Hn=0; h=(b-a)/n; for j=1:n x(j)=a+(j-1/2)*h; y(j)=feval(fx,x(j)); Hn=Hn+y(j); end Hn=h*Hn; S2=(T2+2*Hn)/3; fprintf(' n=%2d S2=%-12.9f S2-S1=%-12.9f\n',n,S2,abs(S2-S1)); if abs(S2-S1) S=S2; % 达到下列条件之一，则运算终止: % (1).abs(S2-S1) 陕西科技大学 机械教改班 用C++的积分 其实积分的思想就是,微分—>求和—>取极限，如果是用纯手工法那就是先对一个函数微分，再求出它的面积，在取极限，因为我们的计算速度和计算量有限，现在有了计算机这个速度很快的机器，我们可以把微分后的每个小的面积加起来，为了满足精度，我们可以加大分区，即使实现不了微分出无限小的极限情况，我们也至少可以用有限次去接近他，下面我分析了四种不同的积分方法，和一个综合通用程序。 一．积分的基本思想 1、思路：微分—>求和—>取极限。 2、Newton —Leibniz 公式 ?-=b a a F b F dx x f ) ()()( 其中，)(x F 被积函数 )(x f 的原函数。 3、用计算机积分的思路 在积分区间内“微分—>求和—>控制精度”。因为计算机求和不可以取极限，也就是不可以无限次的加下去，所以要控制精度。 二．现有的理论 1、一阶求积公式---梯形公式 ?=+-=b a T b f a f a b dx x f )]()([2 )( 他只能精确计算被积函数为0、1次多项式时的积分。 2、二阶求积分公式——牛顿、科特斯公式 ?=+++-=b a S b f a b f a f a b dx x f )]()2(4)([6)( 他只能精确计算被积函数为0、1、2、3次多项式时的积分。 三．四种实现方法 1.复化矩形法 将积分区间[a,b]等分成n 个子区间: ],[],[],[],[],[112322110n n n n x x x x x x x x x x ---、、、 则h=(b-a)/n,区间端点值k x =a+kh 实习七 数值积分的辛普森方法 一、实习目的 1．掌握计算定积分近似值的辛普森方法； 2．理解复化辛普森求积公式。 二、相关知识 抛物线公式（辛普森公式） 将积分区间],[b a 作2n 等分：n i ih a x n a b h i 2,,2,1,0,,2 =+=-=，现在考察由分点22-k x 和k x 2形成的一个小区间],[222k k x x -，（12-k x 为中点），n k ,,2,1 =，在每一个 小区间],[222k k x x -上，作一条抛物线k k k x x y γβα++=2通过三点))(,(2222--k k x f x ， ))(,(1212--k k x f x 和))(,(22k k x f x ，这样就产生关于未知系数k α，k β和k γ的线性方程组 ?????=++=++=++------)() ()(222212122122222222k k k k k k k k k k k k k k k k k k x f x x x f x x x f x x γβαγβαγβα （7-1） 显然上述方程组有唯一解（由高等代数知识知）。 现在，以)(x f y =为顶的曲边梯形用以抛物线k k k x x y γβα++=2为顶的曲边梯形来 代替，其面积 dx x x dx x f k k n k x x k b a k k )()(1 2222γβα++≈∑??=-∑=--=n k k k x x 12226]4)(2)2([){(2222222222222222k k k k k k k k k k k k k x x k x x x x x x γβαγβα++++++++----- }222k k k k k x x γβα+++)]()(4)([6212221222k k k n k k k x f x f x f x x ++-=--=-∑ （7-2） 得抛物线公式，记为n S 2，化简后： {})()(4)(2)(4)(2)(4)(6212432102n n n x f x f x f x f x f x f x f n a b S +++++++-=- 在实际求解数值积分时，我们总是采用成倍加密节点的方法，就抛物线公式而言，若n S 2被认为精度不够，则接着计算n S 4，而精度是否达到要求，又以n n S S 24-是否足够小作为判 在公路中线坐标计算中，我们通常采用切线支距公式来计算曲线上各点的坐标。但当在不同的曲线上计算时就需用不同的计算公式，这为计算也带来不便。在设有缓和曲线的圆曲线半径较小或是卵形曲线上的坐标计算时，如公式选用不当就会出现较大计算误差，即便是能对切线支距公式进行多项展开，也会增加计算的难度。而用复化辛卜生公式不仅能解决不同曲线线型或直线上的坐标计算问题，而且用复化辛卜生公式计算完全是可逆的（即：可顺前进方向也可逆向计算），尤其在计算第二缓和曲线和卵形曲线时显得尤为方便。 用辛卜生公式计算坐标的精度可由人为或程序自行判断，其计算结果完全能保证坐标计算的精度要求。因此，可以说复化辛卜生公式是一个计算公路中线坐标的万能公式。下面本人就该公式在公路中线坐标计算中的具体应用进行实例解析。 一、复化辛卜生公式 式中： H=（Z i－Z A）/n (公式2) (公式3) Zi —待求点桩号 Z A—曲线元起点桩号 Z B—曲线元终点桩号 ρA—曲线元起点曲率 ρB—曲线元终点曲率 a i曲线上任意一点处切线方位角的计算方法有以下三种方法： 1．利用公式（3）求得曲率代入公式（2）计算 2．利用曲线元上已知起点和终点曲率用内插法求得曲率代入公式（2）计算 3．利用切线角公式计算 二、算例 例：已知雅（安）攀（枝花）高速公路西昌西宁立交A匝道一卵形曲线（卵形曲线相关参数见图一，其计算略。），相关设计数据见下表。现用辛卜生公式来计算卵形曲线中桩坐标。 图一 已知相关设计数据见下表： （一）由+271.881推算Zi=+223.715的坐标，n取2等分 用公式(3)、公式（2）计算+247.798处曲线及方位角: ρ+247.798=1÷75+(1÷50-1÷75)(247.798-271.881) ÷(223.715-271.881) =0.01666666666666667 a+247.798=71°24’18.5” +(0.016666667+1÷75)(247.798-271.881)×180÷π÷2 =50°42’26.37” 其它各点依次代入公式计算，结果见下表: 切线方位角图示1 将计算出的数据代入公式（1）求得+223.715中桩坐标如下： X=9880.438+(271.881-223.715)÷2÷6×(cos71°24’18.5”+4(cos61°37’52.22”+cos38°38’0.96”) 辛普森法則 - 數值積分(二) 在許多的實際問題中，常會遇到無法由積分方法求出的定積分。此類的定積分，我們可借由數值方法去求它的近似值，在此僅介紹二種較常用的方法，有興趣的讀者可閱讀數值分析(Numerical analysis)之書籍。 。 辛普森法則(Simpson ’s Rule) 若()f x 在[,]a b 上有定義，將區間[,]a b 分割為n 等分(取n 為偶數)，既 012n a x x x x b =<<<<= ，其中,0,1,,i x a i x i n =+??= ，b a x n -?=。 這裡我們想用過001122(,()),(,()),(,())x f x x f x x f x 三點的拋物線 2()g x x x αβγ=++來取代()f x 在02[,]x x 的定義，進而求出它的近似積分值1A (如圖6-3)，最後用連加的方式求得()f x 在[,]a b 上的近似積分。由假設我們有 20000220202 111122222()()()()22()()f x g x x x x x x x f x g x x x f x g x x x αβγαβγαβγαβγ ?==++??++????==++=++? ? ?????? ?==++? 圖6-3 令 220 0()()x x x x f x dx g x dx ≈?? 2 2x x x x dx αβγ=++? 2 323 2 x x x x x α β γ= + + 3322202020()()()3 2 x x x x x x α β γ= -+ -+- 2 2202 020022()4()322x x x x x x x x x αβγαβγαβγ????++?????=++++++++?? ? ? ? ??????????? 012[()4()()]3 x f x f x f x ?= ++ 所以 240 2 2 ()()n n b x x x a x x x f x dx f x dx -=+++? ??? 012234[()4()()][()4()()] 33 x x f x f x f x f x f x f x ??≈ ++++++ 21[()4()()]3 n n n x f x f x f x --?+++ 01234[()4()2()4()2()3 x f x f x f x f x f x ?= +++++ 212()4()()]n n n f x f x f x --+++ 若令01221[()4()2()2()4()()]3 n n n n x S f x f x f x f x f x f x --?= ++++++ ，且(4)()f x 在[,]a b 上連續，則我們可估計出辛普森法則的誤差值為 5 (4)4()()max |()|180b n n a a x b b a E f x dx S f x n ≤≤-=-≤? 例題1. 試用辛普森法則估計2 1 x e dx -?，取6n =。 解：令2 ()x f x e -=，1 6 x ?= ，則 所以 2 10 1 [140.972620.894840.778820.641240.49940.3679]18 x e dx -≈ +?+?+?+?+?+? 0.7468=复化梯形公式及复化辛普森公式的精度比较
编程MATLAB程序实现复化梯形和辛普森数值积分
9个求积公式
辛普森求积公式
数值分析与实验复化辛卜生公式龙贝格算法
关于辛普森(simpson)公式在线路坐标计算中的应用
复合梯形公式及复合辛普森公式对比
计算方法期中测试(二)答案
高等数值分析拉格朗日插值多项式切比雪夫高斯龙格现象复合梯形辛普森求积公式
辛普森公式
复化梯形公式,辛普森公式的matlab程序
复合辛普森求积
复合梯形公式与复合辛普森公式求积分
复合辛普森公式
变步长复化辛普森公式计算积分
复化梯形法复化矩形法变步长梯形变步长辛普森
数值积分的辛普森方法
复化辛普森公式应用
辛普森积分法则
- 变步长复化辛普森公式计算积分
- 复化辛普森求积公式的Matlab实现
- 复化辛普森公式和高斯求积公式方法计算积分matlab程序
- 复化求积公式
- 编程MATLAB程序实现复化梯形和辛普森数值积分精编版
- 复化求积公式
- 数值分析实习作业之变步长辛普森求积公式
- 数值分析与实验复化辛卜生公式龙贝格算法
- 辛普森求积公式分解
- 编程MATLAB程序实现复化梯形和辛普森数值积分
- 复化辛普森公式和高斯求积公式方法计算积分,matlab程序
- 复化辛普森求积公式
- 复化求积公式计算定积分
- 计算方法期中测试(二)答案
- 数值积分的辛普森方法
- 高斯积分公式第十一次课汇总.
- 复化求积公式求数值积分
- 利用MATLAB软件编写辛普森求积公式程序
- 辛普森公式
- (整理)复化梯形法复化矩形法变步长梯形变步长辛普森