2.cpp

基于龙贝格计算定积分的方法来计算无穷积分

#include<iostream> 
#include<math.h>  
using namespace std; 
double f(double x) 
{ 
	double y; 
	y=1/(x*x+1)/(x*x+4); 
	return y;  
} 
 
int main() 
{ 
	double F,h,T[2][100]={0};
	int m;//求第一列所用，第一列为梯形公式的逐次半分法求的  
	int i=0;  

	int a=0;
	int b=10000;
	double epxl=0.00001;
	int step=100;


	double qk=1;
	while(fabs(T[1][i]-qk)>=epxl)  
	{
		qk=T[1][i];
		m=1;			
		h=(b-a)/2; 
		T[0][0]=h*(f(a)+f(b)); 
		T[1][0]=T[0][0];
		do
		{
			for(int j=0;j<=i;j++) 
				T[0][j]=T[1][j]; 
			F=0;  
			for(int k=1;k<=pow(2,i);k++) 
				F+=f(a+(2*k-1)*h);  
			T[1][0]=T[0][0]/2+F*h; 
			for(i=1;i<=m;i++) 
			{ 
				T[1][i]=pow(4,i)*T[1][i-1]-T[0][i-1]; 
				T[1][i]/=(pow(4,i)-1); 
			} 
			i=i-1;  
			m+=1; 
			h/=2;  
		}while(fabs(T[1][i]-T[0][i-1])>=epxl);  
		b+=step;
	}
	cout<<"I="<<T[1][i]<<"\n"; 
	cout<<"此时的误差为："<<fabs(T[1][i]-qk)<<endl; 
	
	system("pause");  
	return 0;  
}