189 8069 5689
给一个辛普森求积分的例子。
创新互联公司云计算的互联网服务提供商,拥有超过13年的服务器租用、眉山服务器托管、云服务器、虚拟主机、网站系统开发经验，已先后获得国家工业和信息化部颁发的互联网数据中心业务许可证。专业提供云主机、虚拟主机、主机域名、VPS主机、云服务器、香港云服务器、免备案服务器等。
1. 关于辛普森法（Simpson‘s Rule）[1]
辛普森法是数值方法求积分的一种方法，它与梯形法和矩形法不同，辛普森法用二次曲线去逼近实际的积分，如下图所示：
（图片来自维基百科，感谢作者的无私奉献）
2. 代码示例
根据[1]中所给出的伪代码，不难写出C语言版本的辛普森法求积分代码。
#include stdio.h
#include math.h
typedef double (*funcptr)(double a);
double integrate_simpson(funcptr f, float lower, float upper, int n)
{
int i;
float a, b;
float h;
double result;
printf("%f %f %d\n", lower, upper, n);
h = (upper - lower)/n;
printf("%f\n", h);
result = 0;
a = lower;
b = a+h;
for (i=0;in;i++)
{
double tmp;
tmp = (b-a)/6 * ( f(a) + 4*f( (a+b)/2 ) + f(b) );
printf("%f %f %f\n", a, b, tmp);
result = result + tmp;
a = b;
b = a+h;
}
return result;
}
double f(double a)
{
return exp(a);
}
void main(void)
{
double re;
re = integrate_simpson(f, 0, 1, 100);
printf("%f\n",re);
}
3. 输出样式
moose@debian:~$ ./a.out
0.000000 1.000000 100
0.010000
0.000000 0.010000 0.010050
0.010000 0.020000 0.010151
0.020000 0.030000 0.010253
0.030000 0.040000 0.010356
0.040000 0.050000 0.010460
0.050000 0.060000 0.010565
0.060000 0.070000 0.010672
0.070000 0.080000 0.010779
0.080000 0.090000 0.010887
0.090000 0.100000 0.010997
0.100000 0.110000 0.011107
0.110000 0.120000 0.011219
0.120000 0.130000 0.011332
0.130000 0.140000 0.011445
0.140000 0.150000 0.011560
0.150000 0.160000 0.011677
0.160000 0.170000 0.011794
0.170000 0.180000 0.011913
0.180000 0.190000 0.012032
0.190000 0.200000 0.012153
0.200000 0.210000 0.012275
0.210000 0.220000 0.012399
0.220000 0.230000 0.012523
0.230000 0.240000 0.012649
0.240000 0.250000 0.012776
0.250000 0.260000 0.012905
0.260000 0.270000 0.013034
0.270000 0.280000 0.013165
0.280000 0.290000 0.013298
0.290000 0.300000 0.013431
0.300000 0.310000 0.013566
0.310000 0.320000 0.013703
0.320000 0.330000 0.013840
0.330000 0.340000 0.013979
0.340000 0.350000 0.014120
0.350000 0.360000 0.014262
0.360000 0.370000 0.014405
0.370000 0.380000 0.014550
0.380000 0.390000 0.014696
0.390000 0.400000 0.014844
0.400000 0.410000 0.014993
0.410000 0.420000 0.015144
0.420000 0.430000 0.015296
0.430000 0.440000 0.015450
0.440000 0.450000 0.015605
0.450000 0.460000 0.015762
0.460000 0.470000 0.015920
0.470000 0.480000 0.016080
0.480000 0.490000 0.016242
0.490000 0.500000 0.016405
0.500000 0.510000 0.016570
0.510000 0.520000 0.016736
0.520000 0.530000 0.016905
0.530000 0.540000 0.017075
0.540000 0.550000 0.017246
0.550000 0.560000 0.017419
0.560000 0.570000 0.017595
0.570000 0.580000 0.017771
0.580000 0.590000 0.017950
0.590000 0.600000 0.018130
0.600000 0.610000 0.018313
0.610000 0.620000 0.018497
0.620000 0.630000 0.018683
0.630000 0.640000 0.018870
0.640000 0.650000 0.019060
0.650000 0.660000 0.019251
0.660000 0.670000 0.019445
0.670000 0.680000 0.019640
0.680000 0.690000 0.019838
0.690000 0.700000 0.020037
0.700000 0.710000 0.020239
0.710000 0.720000 0.020442
0.720000 0.730000 0.020647
0.730000 0.740000 0.020855
0.740000 0.750000 0.021064
0.750000 0.760000 0.021276
0.760000 0.770000 0.021490
0.770000 0.780000 0.021706
0.780000 0.790000 0.021924
0.790000 0.800000 0.022144
0.800000 0.810000 0.022367
0.810000 0.820000 0.022592
0.820000 0.830000 0.022819
0.830000 0.839999 0.023048
0.839999 0.849999 0.023280
0.849999 0.859999 0.023514
0.859999 0.869999 0.023750
0.869999 0.879999 0.023989
0.879999 0.889999 0.024230
0.889999 0.899999 0.024473
0.899999 0.909999 0.024719
0.909999 0.919999 0.024968
0.919999 0.929999 0.025219
0.929999 0.939999 0.025472
0.939999 0.949999 0.025728
0.949999 0.959999 0.025987
0.959999 0.969999 0.026248
0.969999 0.979999 0.026512
0.979999 0.989999 0.026778
0.989999 0.999999 0.027047
1.718280
另：
matlab中用内置的int求积分的结果：
1.718281828459046
matlab中用前述simpson法求积分的结果：
1.718281828465257
4. 注意
本例中直接调用了标准math库里的exp()，实际上可以通过改写函数f()，将exp也做数值求解。
参考资料：
1. en.wikipedia.org/wiki/Simpson's_rule
这是辛普森积分法。
给你写了fun_1( ），fun_2(），请自己添加另外几个被积函数。
调用方法 t=fsimp(a,b,eps,fun_i);
a,b --上下限，eps -- 迭代精度要求。
#includestdio.h
#includestdlib.h
#include math.h
double fun_1(double x)
{
return 1.0 + x ;
}
double fun_2(double x)
{
return 2.0 * x + 3.0 ;
}
double fsimp(double a,double b,double eps, double (*P)(double))
{
int n,k;
double h,t1,t2,s1,s2,ep,p,x;
n=1; h=b-a;
t1=h*(P(a)+P(b))/2.0;
s1=t1;
ep=eps+1.0;
while (ep=eps)
{
p=0.0;
for (k=0;k=n-1;k++)
{
x=a+(k+0.5)*h;
p=p+P(x);
}
t2=(t1+h*p)/2.0;
s2=(4.0*t2-t1)/3.0;
ep=fabs(s2-s1);
t1=t2; s1=s2; n=n+n; h=h/2.0;
}
return(s2);
}
void main()
{
double a,b,eps,t;
a=0.0; b=3.141592653589793238; eps=0.0000001;
// a definite integral by Simpson Method.
t=fsimp(a,b,eps,fun_1);
printf("%g\n",t);
t=fsimp(a,b,eps,fun_2);
printf("%g\n",t);
// ...
printf("\n Press any key to quit...");
getch();
}
对于一重定积分来说其求解可以使用梯形法进行求解，计算公式如下所示：
其中，f(x)为被积函数，为横坐标的两点间的间隔，越小，则计算出的结果越精确。
对于求解此类问题可以使用C语言中的回调函数编写通用的计算函数，代码如下：
#include stdio.h
#include stdlib.h
#includemath.h
//功能：返回f(x)在积分区间[a,b]的值
//参数：FunCallBack 指向用于计算f(x)的函数
// a 积分区间的起始值
// b 积分区间的结束值
// dx 横坐标的间隔数，越小计算结果越准确
double Calculate(double (*FunCallBack)(double x),
double a,double b,double dx)
{
double doui;
double total = 0; //保存最后的计算结果
for (doui = a; doui = b; doui += dx)
{
total += FunCallBack(doui)*dx;
}
return total;
}
double f2(double x)
{
return x*x;
}
double f(double x)
{
return x;
}
double f3(double x)
{
return x*x*x ;
}
int main()
{
double total;
total = (Calculate(f, 2, 3, 0.000001));
printf("total = %lf\n", total);
total = (Calculate(f2, 2, 3, 0.000001));
printf("total = %lf\n", total);
total = (Calculate(f3, 2, 3, 0.000001));
printf("total = %lf\n", total);
return 0 ;
}
其中，函数f，f2，f3为自行编写的关于x的被积函数。
运行结果：
total = 2.500000
total = 6.333331
total = 16.249991
#include stdio.h
#define RES (1e-6)
double integ(double a,double b,double f(double))
{
double sum;
for(sum=0;ab;a+=RES)
{
sum+=f(a)*RES;
}
return sum;
}
double f(double x)
{
return x*x;
}
int main()
{
printf("%lf\n",integ(0,0.1,f));
return 0;
}
#include
#include
double integral(double(*fun)(double x),double a,double b,int,n){
double s,h,y;
int i;
s=(fun(a)+fun(b))/2;
h=(b-a)/n; /*积分步长*/
for(i=1;in;i++)
s=s+fun(a+i*h);
y=s*h;
return y;/*返回积分值*/
}
double f(double x){
return(x*sinx) /*修改此处可以改变被积函数*/
}
int main(){
double y;
y=integral(f,1.0,2.0,150);/*修改此处可以改变积分上下限和步数,步长=(上限-下限)/步数*/
printf("y=%f\n",y);
return 0;
}
int main()
问题就是出在数据类型上的选用上，precision=0.0000001时已经超过了float的数据范围，所以导致数据截断后precision=0.000000，从而程序在计算积分时可能陷入死循环，应该采用double型数据类型。其实不推荐楼主用如此多的define语句，程序的可读性和风格应该重于编程员的劳动度。。。
还有楼主对自然对数e的define也已经超过了计算机的可识别范围。。您那样精确的定义e并不会在结果上获得更加精确地结果，其实反倒会起到相反的作用，要知道与其用一个这样可能导致内存出错以及必定会导致数据截断的变量来实现精度的提高远远不如采用一个更精确的积分算法，而且c语言提供了自然数e为底的指数函数~而且貌似您的积分算法是不准确的，梯形积分的定义并非如此，其再两端的函数值应该只取1/2.希望您多加细心~
如果不介意的话，就是你的precision应该改为step~这样会能更加准备的表达了这个变量的作用，在你的程序中precision变量其实是积分步长~在数值计算方法中积分精度的控制往往不是通过细化步长来表达，而是通过后一个积分值-前一个积分值precision 这样来实现精度控制~呵呵