重积分_070460.cpp

利用计算机进行二重数值积分

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

long double f(long double x,long double y)
{
	return tan(x*x+y*y);
}


int ss(int r)
{
	int ssum=1;
	for(int i=1;i<=r;i++)
     ssum*=2;
	return ssum;
}

double T(int m,int n,long double a,long double b,long double c,long double d)
{
	long double Tsum=0;
	long double h=(b-a)/m;
	long double k=(d-c)/n;
	for(int i=0;i<=n-1;i++)
	{
		for(int j=0;j<=m-1;j++)
		{
		
			Tsum+=0.25*h*k*(f(a+i*h,c+j*k)+f(a+i*h,c+(j+1)*k)+f(a+(i+1)*h,c+j*k)+f(a+(i+1)*h,c+(j+1)*k));
		}
	}
	return Tsum;
	
}

double TT(long double T,int m,int n,long double a,long double b,long double c,long double d)
{
	long double TT1=0.0,TT2=0.0,TT3=0.0,TT4=0.0,TT5=0.0;
	for(int i=0;i<=m-1;i++)
		TT1+=f(a+(i+0.5)*(b-a)/m,c)+f(a+(i+0.5)*(b-a)/m,d);
	for(int j=0;j<=n-1;j++)
		TT2+=f(a,c+(j+0.5)*(d-c)/n)+f(b,c+(j+0.5)*(d-c)/n);
	for(int k1=0;k1<=m-1;k1++)
		for(int k2=1;k2<=n-1;k2++)
			TT3+=f(a+(k1+0.5)*(b-a)/m,c+k2*(d-c)/n);
	for(int l1=1;l1<=m-1;l1++)
		for(int l2=0;l2<=n-1;l2++)
			TT4+=f(a+l1*(b-a)/m,c+(l2+0.5)*(d-c)/n);
	for(int m1=0;m1<=m-1;m1++)
		for(int m2=0;m2<=n-1;m2++)
			TT5+=f(a+(m1+0.5)*(b-a)/m,c+(m2+0.5)*(d-c)/n);
	return 0.25*T+0.125*((b-a)/m)*((d-c)/n)*(TT1+TT2+2*(TT3+TT4+TT5));
}
long double waitui(long double t[20],long double t1[20],long double t2[20],long double t3[20],
				   long double a,long double b,long double c,long double d,int z,int x,int m,int n,double e)
{
	
	long double eps,eps1;
	int r=4;
    while(r<50)
	{
       t[r]=TT(t[r-1],ss(r-1)*m,ss(r-1)*n,a,b,c,d);
	   t1[r-1]=t[r]*4/3-t[r-1]/3;
	   t2[r-2]=t1[r-1]*16.0/15-t1[r-2]/15;
	   t3[r-3]=t2[r-2]*64.0/63-t2[r-3]/63;
	   eps=(t3[r-3]-t3[r-2])/255;
	   eps1=(t2[r-2]-t2[r-1])/63;
	   z=ss(r+1)*m;x=ss(r+1)*n;
	   r++;
	   if(fabs(eps)>e&&fabs(eps1)>e)
	   {
		  cout<<"在区间["<<a<<","<<b<<"],["<<c<<","<<d<<"]分别各等分"<<z<<","<<x<<endl;
		  return t3[r-4];}
	}
}

void main()
{	
    long double a=0.0,b=0.523598775,c=0,d=1.047197551,e,h,k;
	long double t[20],t1[20],t2[20],t3[20];
	int m=1,n=1,z=0,x=0;

    cout<<"输入区间划分初始间隔h,k的值:"<<endl;
	cin>>h>>k;
	cout<<"输入允许误差:"<<endl;
	cin>>e;
	m=(int)((b-a)/h);                                                                                                                                          
	n=(int)((d-c)/k);

	t[0]=T(m,n,a,b,c,d);
	t[1]=TT(t[0],m,n,a,b,c,d);
	t1[0]=t[1]*4.0/3-t[0]/3;
	if((fabs(t[1]-t[0])/3)<e)
	{cout<<"积分结果:"<<setprecision(10)<<t1[0]<<endl;z=ss(2)*m;x=ss(2)*n;
	cout<<"在区间["<<a<<","<<b<<"],["<<c<<","<<d<<"]分别各等分"<<z<<","<<x<<endl;}
	else{
		t[2]=TT(t[1],ss(1)*m,ss(1)*n,a,b,c,d);
		t1[1]=t[2]*4.0/3-t[1]/3;
		t2[0]=t1[1]*16.0/15-t1[0]/15;
		if((fabs(t1[1]-t1[0])/15)<e)
		{cout<<"积分结果:"<<setprecision(10)<<t2[0]<<endl;z=ss(3)*m;x=ss(3)*n;
	    cout<<"在区间["<<a<<","<<b<<"],["<<c<<","<<d<<"]分别各等分"<<z<<","<<x<<endl;}
		else{
		    t[3]=TT(t[2],ss(2)*m,ss(2)*n,a,b,c,d);
			t1[2]=t[3]*4.0/3-t[2]/3;
		    t2[1]=t1[2]*16.0/15-t1[1]/15;
			t3[0]=t2[1]*64.0/63-t2[0]/63;
			if((fabs(t2[1]-t2[0])/63)<e)
			{cout<<"积分结果:"<<setprecision(10)<<t3[0]<<endl;z=ss(4)*m;x=ss(4)*n;
			cout<<"在区间["<<a<<","<<b<<"],["<<c<<","<<d<<"]分别各等分"<<z<<","<<x<<endl;}
			else 
				cout<<"积分结果:"<<setprecision(10)<<waitui(t,t1,t2,t3,a,b,c,d,z,x,m,n,e)<<endl;
		}
	}
	
}