3.txt

romberg算法的C++实现 很好用

// Romberg.cpp : Defines the entry point for the console application. 
// 

#include "stdafx.h" 

const double e = 1.0E-8; 
const double end = 1.0E-6; 

//积分上下限 
const double a = 0.0; 
const double b = 1.0; 

//被积函数 
double function(long double x) 
{ 
if(abs(x)<e) return 1.0; 
else return sin(x)/x; 
} 

int min(int x, int y) 
{ 
return x>y?y:x; 
} 

int _tmain(int argc, _TCHAR* argv[]) 
{ 
double * temp; 
double * t = new double[1]; 
int m = 1; 
int j = 0; 
double h; 
t[0] = (b-a)/2*(function(a)+function(b));//t_0[0] 
while(true) 
{ 
temp = new double[m];//创建数组temp[m] 
int k = sizeof(double)*min(4,m); 
::memcpy(temp,t,k);//将t数组的数拷到temp数组 
delete []t; 
t = new double[m+1];//当前行比上一行多一个数 
h = (b-a)/pow(2.0,m);//步长 
t[0] = temp[0]/2; 
for(int i = 1;i<=pow(2,m-1);i++) 
{ 
t[0] += h*function(a+(2*i-1)*h);//相当于t_k[0]; 
} 
for(int i = 1;i<=min(3,m);i++) 
{ 
t[i] = (pow(4.0,i)*t[i-1]-temp[i-1])/(pow(4.0,i)-1);//t_m[i] 
} 

//产生4行数后判断精度是否满足要求 
if(m>3) 
{ 
if(abs(t[3]-temp[3])/t[3]<end)//判断第4列（八阶收敛） 
{ 
//满足精度要求，输出结果，跳出循环 
printf("%d\n",m);//循环次数 
printf("%2.15f\n",t[3]);//积分结果 
// std::cout<<t[3]<<std::endl;//这是C++的输出函数。 
break; 
} 
} 
delete []temp; 
m++; 
} 
delete []t; 
return 0; 
} 
//由于第k+1行的数可以由第k行的数求出，与1，2,...,k-1行没有直接关系，可以用两个数组存储计算结果，前面算过的数不符合要求就 
//直接删掉了。这样可以节约内存。这两个数组的长度是变化的，随着数表行数逐渐加长。