原创龙贝格积分

 2008-4-13 22:52  2458 5 5 分类: 工业电子

/********************************************************************
* 用龙贝格法计算函数f(x)从a到b的积分值
* 输入: f--函数f(x)的指针
*       a--积分下限
*       b--积分上限
*       eps--计算精度
*       max_it--最大迭代次数
* 输出: 返回值为f(x)从a点到b点的积分值
*******************************************************************/
double romberg(double(*f)(double),double a,double b,
                          double eps,int max_it)
{
double *t,h;
int i,m,k;

if(!(t=(double *)malloc(max_it*sizeof(double)+1)))
return(ERROR_CODE);
h=b-a;
t[1]=h*((*f)(a)+(*f)(b))/2.0;
printf("%18.10e\n",t[1]);
for(k=2;k<max_it+1;k++)
{
double s,sm;
h*=0.5;
s=0.0;
for(i=0;i<pow(2,k-2);i++)
s+=(*f)(a+(2*i+1)*h);
sm=t[1];
t[1]=0.5*t[1]+h*s;
for(m=2;m<k+1;m++)
{
s=t[m];
t[m]=t[m-1]+(t[m-1]-sm)/(pow(4,m-1)-1);
if(m<k)
sm=s;
}
for(m=1;m<k+1;m++)
printf("%18.10e",t[m]);
printf("\n");
if(fabs(t[k]-sm)<eps)
{
sm=t[k];
free(t);
return(sm);
}
}
return(ERROR_CODE);
}

写原创有奖励！2024面包板原创奖励正在进行中