数据结构算法:龙贝格积分法
发布时间：2010/11/10 15:25:17 来源：www.xue.net 编辑：城市总裁吧
估计很多人对龙贝格积分法不是很连接吧，下面我一起来看看龙贝格积分法代码：
#include <math.h>
#include <stdio.h>
#define eps 0.0000001
#define max 20
double f(double x)
{
if(x==0)
return 1;
else
return (sin(x)/x);
}
void romberg(double a,double b)
{
double t[max][4]={0},h=1.0,e=1.0+eps;
double fnew;
int i,j,k=1,m;
t[0][0]=h*(f(a)+f(b))/2.0;
while((k<max)&&(e>eps))
{
fnew=0;
for(i=1;i<=(int)(pow(2,k-1));i++)
fnew=fnew+f(a+(i-0.5)*h);
t[k][0]=(t[k-1][0]+h*fnew)/2.0;
for(m=1;m<=k;m++)
{
if(m>3)
break;
t[k][m]=(pow(4,m)*t[k][m-1]-t[k-1][m-1])/(pow(4,m)-1);
}
if(k>=4)
e=fabs(t[k][3]-t[k-1][3]);
k++;
h=h/2.0;
}
if(k>max)
printf(\"method failed.\\n\");
else
{
printf(\"\\tt\\t\\ts\\t\\tc\\t\\tr\\n\");
for(i=0;i<k;i++)
{
printf(\"k=%d\\t\",i);
for(j=0;j<4;j++)
if(i>=j)
printf(\"%0.9lf\\t\",t[i][j]);
printf(\"\\n\");
}
}
}
int main()
{
double a=0.0,b=1.0;
romberg(a,b);
return 0;
}
