c程序romberg求积分

一剑行天下 posted @ 2008年11月07日 16:13 in 数值分析 , 1879 阅读

#include<stdio.h>
#include<math.h>
double T[100][100];
double f(double x)
{
double f1=0;
if(x==0)
f1=0;
else
f1=pow(x,1.5);
return f1;
}
main()
{
double a=0,b=0,e=0;
double ff=0;
int i=0,j=0;
int m=0,k=1;
double h,T0;
printf("请输入初始点a和末尾点b的值:\n");
scanf("%lf,%lf",&a,&b);
printf("请输入误差e\n");
scanf("%lf",&e);
h=(b-a)/2;
T0=T[0][0]=h*(f(a)+f(b));
printf("i=%d 时T=%.7lf\n",m-1,T[0][0]);

do
{
ff=0;
m++;
for(k=1;k<=(pow(2,m-1));k++)
{
   ff=ff+f(a+(2*k-1)*h);
}
T[0][m]=T0/2 +h*ff;
printf("t[0][%d]=%lf\n",m,T[0][m]);
for(j=1;j<=m;j++)
{
   T[j][m-j]=(pow(4,j) *T[j-1][m-j+1] - T[j-1][m-j] ) / (pow(4,j)-1);
   printf("m=%d 时T[j][m-j]=T[%d][%d]=%.6lf\n",m,j,m-j,T[j][m-j]);
}
T0=T[0][m];
h=h/2;
}while (fabs((T[m][0]-T[m-1][0]))>=e);
printf("对于要满足的误差e＝%.6lf进行了%d次加速运算\n",e,m);
printf("T[%d][0]=%.7lf\n",m,T[j][0]);

}
/* 运行结果：
请输入初始点a和末尾点b的值:
0,1
请输入误差e
0.000007
i=-1 时T=0.5000000
t[0][1]=0.426777
m=1 时T[j][m-j]=T[1][0]=0.402369
t[0][2]=0.407018
m=2 时T[j][m-j]=T[1][1]=0.400432
m=2 时T[j][m-j]=T[2][0]=0.400303
t[0][3]=0.401812
m=3 时T[j][m-j]=T[1][2]=0.400077
m=3 时T[j][m-j]=T[2][1]=0.400054
m=3 时T[j][m-j]=T[3][0]=0.400050
t[0][4]=0.400463
m=4 时T[j][m-j]=T[1][3]=0.400014
m=4 时T[j][m-j]=T[2][2]=0.400009
m=4 时T[j][m-j]=T[3][1]=0.400009
m=4 时T[j][m-j]=T[4][0]=0.400009
t[0][5]=0.400118
m=5 时T[j][m-j]=T[1][4]=0.400002
m=5 时T[j][m-j]=T[2][3]=0.400002
m=5 时T[j][m-j]=T[3][2]=0.400002
m=5 时T[j][m-j]=T[4][1]=0.400002
m=5 时T[j][m-j]=T[5][0]=0.400002
t[0][6]=0.400030
m=6 时T[j][m-j]=T[1][5]=0.400000
m=6 时T[j][m-j]=T[2][4]=0.400000
m=6 时T[j][m-j]=T[3][3]=0.400000
m=6 时T[j][m-j]=T[4][2]=0.400000
m=6 时T[j][m-j]=T[5][1]=0.400000
m=6 时T[j][m-j]=T[6][0]=0.400000
对于要满足的误差e＝0.000007进行了6次加速运算
T[6][0]=0.0000000

一剑行天下's Blog

Happy coding

一剑行天下

分类

最新评论

最新留言

链接

功能

c程序romberg求积分