这个是我网上找的关于龙贝格的代码,给你参考参考,因为我没听过这个名词··帮不到你什么忙···
你的代码中
程序代码:
for(i=1;i<n;i++)
{for(j=1;j<=n-i+1;j++)
T[i+1][j-1]=(pow(4,i)*T[i][j]-T[i][j-1])/(pow(4,i)-1);
}
printf("%f",T[i][0]);
这里有问题,最后一次大循环时红色部分已经越界了,而且到printf那里 i的值已经等于n了,
printf("%f",T[i][0]);肯定会越界,因为不了解这个算法,所以也不知道怎么改~~~LZ看看吧···
下面是我网上找的有关代码 ,LZ试试···
程序代码:
#include<iostream>
#include<cmath>
using namespace std;
#define f(x) (4.0/(1+x*x)) //举例函数
#define epsilon 0.0001 //精度
#define MAXREPT 10 //迭代次数,到最后仍达不到精度要求,则输出T(m=10).
double Romberg(double aa, double bb)
{ //aa,bb 积分上下限
int m, n;//m控制迭代次数, 而n控制复化梯形积分的分点数. n=2^m
double h, x;
double s, q;
double ep; //精度要求
double *y = new double[MAXREPT];//为节省空间,只需一维数组
//每次循环依次存储Romberg计算表的每行元素,以供计算下一行,算完后更新
double p;//p总是指示待计算元素的前一个元素(同一行)
//迭代初值
h = bb - aa;
y[0] = h*(f(aa) + f(bb))/2.0;
m = 1;
n = 1;
ep = epsilon + 1.0;
//迭代计算
while ((ep >= epsilon) && (m < MAXREPT))
{
//复化积分公式求T2n(Romberg计算表中的第一列),n初始为1,以后倍增
p = 0.0;
for (int i=0; i<n; i++)//求Hn
{
x = aa + (i+0.5)*h;
p = p + f(x);
}
p = (y[0] + h*p)/2.0;//求T2n = 1/2(Tn+Hn),用p指示
//求第m行元素,根据Romberg计算表本行的前一个元素(p指示),
//和上一行左上角元素(y[k-1]指示)求得.
s = 1.0;
for (int k=1; k<=m; k++)
{
s = 4.0*s;
q = (s*p - y[k-1])/(s - 1.0);
y[k-1] = p;
p = q;
}
p = fabs(q - y[m-1]);
m = m + 1;
y[m-1] = q;
n = n + n; h = h/2.0;
}
return (q);
}
int main()
{
double a,b;
cout<<"Romberg积分,请输入积分范围a,b:"<<endl;
cin>>a>>b;
cout<<"积分结果:"<<Romberg(a, b)<<endl;
system("pause");
return 0;
}