romberg.cpp

计算方法经典算法(弦截法二分法牛顿迭代法龙贝格算法)

//用龙贝格算法求a到b间函数2*x*x*x的积分


#include<iostream.h>
#include<math.h>
#include<stdlib.h>

double f(double x)
{
	return 2*x*x*x;
}


void rombeg(double a,double b,double l)//参数分别为积分的上限下限和精度
{
	double t[10][10];                  //t[10][10]用来保存T[]，S[]，C[]，R[]
	double h,m;      
	int k=1;
	h=b-a;
	t[0][0]=(f(a)+f(b))*h/2;
	while(1){
		double index=0;
		for(int i=1;i<=pow(2,k-1);i++)
		{
			m=a+(i-0.5)*h;
			index+=f(m);
		}
		t[0][k]=(t[0][k-1]+h*index)/2;
		for(int j=1;j<=k;j++)
			t[j][k-j]=(pow(4,j)*t[j-1][k-j+1]-t[j-1][k-j])/(pow(4,j)-1);
		if(fabs(t[k][0]-t[k-1][0])<l)
		{
			cout<<"I="<<t[k][0]<<endl;    //输出中间结果
			for(int s1=0;s1<=k;s1++)      //
				for(int s2=s1;s2<k;s2++)  //
					cout<<t[s1][s2]<<endl;//输出中间结果
			exit(1);
		}
		else
		{   k=k+1;h=h/2;}              //精度不够，继续
	}
}


void main()
{
	cout<<"please input a,b,l:"<<endl;
	double a,b,l;
	cin>>a>>b>>l;
	rombeg(a,b,l);
}