龙贝格公式求积分的c语言程序,小生提问c语言编的龙贝格算法求积分?大家看看哪儿错了,运行结果不对...

该楼层疑似违规已被系统折叠 隐藏此楼查看此楼

应该是在算T[k][0]的地方，但是找不出来. 求高人指点迷津

#include

#include

#include

#define PI 3.1415926536

#define R 6371

#define H 2384

#define h1 439

#define a (2*R+H+h1)/2

#define c (H-h1)/2//有关常数值

#define FUN(x) pow(1 - pow(c/a,2) * pow(sin (x),2),0.5)//积分函数

#define N 20 //最大加速次数

#define FAC 4*a//积分函数前面的系数

#define ZERO 1.564640219//第一次用梯形公式时积分函数值

void main()

{

double A = 0,B = PI/2;//积分区间

double T[N][N];

int k,m;

int i,j;

T[0][0] = ZERO;

for(k = 1;k

{

double t = 0,h = (B - A) * pow(0.5,k);

double x;

int n = (int)pow(2,k-1);

#ifdef DEBUG_M

printf("%d\t %lf\n",n,h);

#endif

for(m = 0;m < n;m++)

{

x = h + 2 * h * m;

t += h * FUN(x);

}

T[k][0] = 0.5 * T[k-1][0] + t;

}//求出梯形值T((B-A)*pow(2,k))，即T[k][0]

for(k = 0;k

{

for(m = 1;m<=k;m++)

T[k][m] = pow(4,m)/(pow(4,m)-1) * T[k][m-1] - 1/(pow(4,m)-1) * T[k-1][m-1];

if(fabs(T[k][k] - T[k-1][k-1]) <=0.0001)

break;

}//逐个求出各行其他元素T[i][j]

for(i = 0;i<=k;i++)

for(j = 0;j<=i;j++)

T[i][j] *= FAC;

for(i = 0;i<=k;i++)

{

for(j = 0;j<=i;j++)

printf("%.2lf\t\t",T[i][j]);

printf("\n");

}

}