C语言用辛普森公式求sinx在0到π上的定积分的源程序

用梯形法估算，再用辛普森法。

fsimpf 积分函数

a,b 积分下上限，旅卖eps 精度。

#include<stdlib.h>

#include <math.h>

double fsimpf(double x)

{

return sin(x)

}

double fsimp(double a,double b,double eps)

{

int n,k

double h,t1,t2,s1,s2,ep,p,x

n=1h=b-a

t1=h*(fsimpf(a)+fsimpf(b))/2.0

s1=t1

ep=eps+1.0

while (ep>=eps)

{

p=0.0

for (k=0k<=n-1k++)

{

x=a+(k+0.5)*h

p=p+fsimpf(x)

}

t2=(t1+h*p)/2.0

s2=(4.0*t2-t1)/3.0

ep=fabs(s2-s1)

t1=t2s1=s2n=n+nh=h/2.0

}

return(s2)

}

void main()

{

double a,b,eps,t

a=0.0b=3.141592653589793238eps=0.0000001

a definite integral by Simpson Method.

t=fsimp(a,b,eps)

printf("拆姿逗%e\n",t)

printf("册昌\n Press any key to quit...")

getch()

}

辛普森（Simpson）公式是牛顿-科特斯公式当n=2时的情形，也称为三点公式。利用区间二等分的三个点来进行积分插值。其科特斯系数分别为1/6，4/6，1/6。

中文名

辛普森公式

别称

三点公式

应用学科

立体几何

科特斯系数

1/6，含戚桐4/6，1/6

应用

立体几何中用来求拟柱体体积的公式。

公式内容

设拟柱体的高（两底仔春面α，β间的距离）为H，如果用平行于底面的平面γ去截该图形，所得到的截面面积是平面γ与平面α之间距离h的不超过3次的函数，那么该拟柱体的体积V为

V = H (S_1 + 4S_0 + S_2) /6.

式中，S_1和S_2是两底面的面积，S_0是中截面的面积（即平面γ与平面α之间距离h=H/2时得到的截面的面积）。

事实上，不光是拟柱体，其他符合条件（所有顶点都在两个平行平面上、用平行于底面的平面去截该图形时所得到的截面面积是该平面与一底之间距离的不超过3次的函数）的立体图形也可以利用该公式求体积。

计算实例

例1：计算底面积为S、高为h的谈坦柱体的体积。

解：此题中S_1 = S_0 = S_2 = S，H = h，所以V = H (S_1 + 4S_0 + S_2) /6 = h (S + 4S + S) /6 = S h。

例2：计算底面积为S、高为h的锥体的体积。

解：此题中S_1 = S，S_0 = S /4，S_2 = 0，H = h，所以V = H (S_1 + 4S_0 + S_2) /6 = h (S + 4S /4 + 0) /6 = S h /3。

例3：计算半径为r的球的体积。

解：此题中S_1 = S_2 = 0，S_0 = πr^2，H = 2r，所以V = H (S_1 + 4S_0 + S_2) /6 = 2r (0 + 4πr^2 + 0) /6 = 4πr^3 /3。

公式证明

只需要证明根据公式算出来的体积和用积分算出来的体积相等即可。

设截面面积是截面高h的不超过3次的函数：f(h)= ah^3 + bh^2 + ch + d。

那么，

利用积分计算体积，可以得到（积分限为0~h）：

V = ∫ f(x) dx

= ah^4 /4 + bh^3 /3 + ch^2 /2 +dh；

利用公式计算体积，可以得到：

V = H (S_1 + 4S_0 + S_2) /6

= h ( f(0) + 4f(h/2) + f(h) ) /6

= h [ d + 4 (ah^3 /8 + bh^2 /4 + ch /2 + d) + (ah^3 + bh^2 + ch + d) ]/6

= ah^4 /4 + bh^3 /3 + ch^2 /2 +dh。

因此两式相等，公式得证。

Remark：当函数f(h)次数高于或等于4次时，公式一般不成立。这只需验证f(h)=h^4时公式不成立即可。

---