重庆分公司,新征程启航
为企业提供网站建设、域名注册、服务器等服务
我给一楼加的注释以及修改:
创新互联公司专注于企业成都全网营销推广、网站重做改版、虎丘网站定制设计、自适应品牌网站建设、H5高端网站建设、商城开发、集团公司官网建设、成都外贸网站制作、高端网站制作、响应式网页设计等建站业务,价格优惠性价比高,为虎丘等各大城市提供网站开发制作服务。
#includestdio.h
#includemath.h
#define ARRAYBOUND 10001
void main()
{
int i = 0; //辅助变量,最常见那种
int n = 0; //将所求定积分函数曲线在x轴方向,平均分成n等分;n越大,结果越精确;不过限于此算法限制nARRAYBOUND,否则溢出.
float x[ARRAYBOUND];//ARRAYBOUND维浮点数组,存放离散的x坐标值
float y[ARRAYBOUND];//ARRAYBOUND维浮点数组,存放每个x坐标对应的函数值;x[i],y[i]满足y[i]=f(x[i]),f是你要求定积分的函数
float x0 = 0.0; //定积分下限
float xn = 0.0; //定积分上限
float h = 0.0; //面积微元宽度
float J = 0.0; //辅助变量
/*f=x^3*/ //这里说明要求定积分的是函数f(x)=x*x*x;(y等于x的立方,x^3是vb的写法)
// printf("input x0,xn,n:");
printf("请分别输入下限(x0),上限(xn),精度(n):");
scanf("%f",x0);
scanf("%f",xn);
scanf("%d",n);
h=(xn-x0)/n;//将函数图形在x方向平分成n份,h是每个面积微元的宽度
x[0]=x0; //将积分下限赋值给x[0]
for(i=0;i=n nARRAYBOUND;i++)
{
x[i]=x[0]+i*h; //计算n个离散的横坐标值,存入x[]数组
y[i]=(float)pow(x[i],3);//计算n个横坐标对应的函数值,存入y[]数组。在此可以改变要求积分的函数
}
// J=0.0;
for(i=0;in;i++)
{
//J=J+y[i]+y[i+1];
J+=y[i];//将所有纵坐标值代数相加,存入J
}
//J=J*h/2.0;
J=J*h;//所有微元面积一次求解,因为∑h*y[i]=h*∑y[i];
printf("\nn=%d \n所求定积分值是: %f\n",n,J);
}
我将//J=J+y[i]+y[i+1]改为J+=y[i];将//J=J*h/2.0;改为J=J*h只是帮助lz理解
其实,这两种表达在理论上是等价的,不过我发现修改后,在n同样大小的情况下,结果的精度有一点点下降,还真不知为什么???
这样的话lz应该能理解了吧,其实一楼的算法还有不少值得改进的地方,希望lz能有所突破!!
看不出有什么编译或执行的问题。但您的integral函数中的循环控制有点小问题,应该是
for ( i=a;ib;i+=trace )
您现在的写法,多加了一次计算,不知道这个对运算结果会不会造成影响。
积分分为两种,数值积分,公式积分。
公式积分:部分函数可以直接用公式求得其不定积分函数。C语言中可以直接用积分公式写出其积分函数。
数值积分:按照积分的定义,设置积分范围的步长,用梯形面积累加求得其积分。
以【f(x)=x*sin(x) 从1到2的积分】为例:
#include math.h
#include stdio.h
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;
}
这是辛普森积分法。
给你写了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();
}
给一个辛普森求积分的例子。
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
基本是这样的,用梯形发求定积分,对应于一个积分式就要有一段程序,不过你可以改变程序的一小部分来改变你所要求的积分式。
以c为例:求f(x)=xsinx从1到2的积分
#include math.h
float integral(float(*fun)(float x),float a,float b,int,n)
{float 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;/*返回积分值*/
}
float f(float x)
{return(x*sinx) /*修改此处可以改变被积函数*/
}
main()
{float y;
y=integral(f,1.0,2.0,150);/*修改此处可以改变积分上下限和步长*/
printf("y=%f\n",y);
}