chap5.h

romberger 积分

#include<iostream.h>
#include<math.h>
#include<fstream.h>
#include<iomanip.h>


int const max_step=100; //定义一个常数，代表最大循环步长
                        //在本次程序能达到的显示精度之下，100作为最大循环步长已经完全足够了！

class romberg 
{
public:
	romberg(long double eps);	   //类的构造函数，传递参数为由用户输入椭圆的参数
		
	long double f(long double x);  //求被积函数f(x)的函数值
	long double computeT(long double aa, long double bb, long int n); 
	                               //计算梯形积分，其中aa和bb为区间的两个端点，n为等分数
		
	void output();	               //用romberg积分法进行计算，并将结果输出在屏幕上
	
private:
	int n;			               //区间等分数（或者说节点数-1）
	long double a,b;               //被积区间的两个端点
	long double epsilon;           //积分能达到的精度
	
	long double R[max_step][2];    //数组R，用于romberg积分中的存储
};


romberg::romberg(long double eps)  //构造函数
{
	n=1;                           //初始步长从1开始
	
	//本程序先求出1/4的椭圆周长，然后再乘以4得出总的周长
	//因此积分区间为0到pi/2
	a=0,b=3.1415926/2;             
    
	//将精度按照用户输入的值进行设置
	//由于先求出的是1/4周长，因此精度也相应的变为原来精度的1/4
	epsilon=eps/4;                 

	for(int i=0;i<2;i++)            //将数组清零
		for(int j=0;j<max_step;j++)
			R[j][i]=0;
}


long double romberg::f(long double x)         //求被积函数f(x)的函数值
{
	long double y=sqrt((pow(sin(x),2.0)*3+1));
	return y;
}


long double romberg::computeT(long double aa, long double bb, long int n)
//计算梯形积分（用复化梯形积分公式求得）
//传递的参数为：aa和bb代表积分区间的两个端点值，n为区间等分数
//返回n等分的梯形积分值
{
	long double h=(bb-aa)/n;          //h为将区间n等分后的步长
	long double sum=0 ;               //sum用于累加积分值
	for(int i=1;i<n;i++)
		sum+=f(aa+i*h);
	
	sum+=(f(aa)+f(bb))/2;             //不要忘记加上两个端点函数值的一半
	return(h*sum);                    //返回求出的梯形积分值，要乘以步长，即(bb-aa)/n
}



void romberg::output()                //用romberg积分法进行计算，并将结果输出在屏幕上
{
	cout<<setiosflags(ios::fixed);
	cout<<setprecision(14);           //结果显示到小数点后8位
	
	R[0][1]=computeT(a,b,n);          //先计算1等分的时候的梯形积分值
	n*=2;                             

	for(int m=2;m<max_step;m++)
	{
		int i;

		//从1到m-1开始循环
		for(i=0;i<m-1;i++)
		{
			R[i][0]=R[i][1];
			
			//这里的输出是在调试程序的时候进行检测的,也可以用来检测收敛的速度
			//cout<<R[i][0]<<endl;
		}
        
	
       	R[0][1]=computeT(a,b,n);
		n*=2;
		
		for(i=1;i<m;i++)
		//外推公式
		{
			R[i][1]=R[i-1][1]+(R[i-1][1]-R[i-1][0])/(pow(4,i)-1);
		}

		//这里的输出是在调试程序的时候进行检测的,也可以用来检测收敛的速度
		/*for(i=0;i<m;i++)
		{
			cout<<R[i][1]<<endl;
		}
		*/

				
        //当最后一次的积分值与前次之差的绝对值小于要求的精度时，返回最后一次积分值 
		if( ( (R[m-1][1]-R[m-2][0]) < epsilon ) && ( (R[m-1][1]-R[m-2][0]) > (-1)*epsilon )) 
		{
			cout<<"在设置的精度下求得椭圆周长为： "
				<<4*R[m-1][1]<<endl;   //注意，由于求出来得是1/4周长值，所以返回的时候要再乘以4
			cout<<endl;
			return;
		}
		//cout<<endl;
	}

    //若循环达到最大步长数仍旧没有求出所要的解时，报告无解的信息
	cout<<"无解！请重新输入有效数位数"<<endl;
}