感冒,鼻子流黄绿色鼻涕:龙贝格程序例题

// 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行没有直接关系，可以用两个数组存储计算结果，前面算过的数不符合要求就
//直接删掉了。这样可以节约内存。这两个数组的长度是变化的，随着数表行数逐渐加长。