复化梯形公式+二分+外推-计算二重积分

//Calculate.cpp:定义控制台应用程序的入口点。//如需要重新使用本程序，只需要修改程序开始处的积分区间和积分函数即可。#includestdafx.hintT_Rank=-1;//记录当前计算的阶数doubleT_Matrix[16][16];#definePI3.14159265358979323846#definex00.0#definexmPI/3#definey00.0#defineynPI/6#defineeps00.5e-5doublef(doublex,doubley){doubleSum;Sum=tan(x*x+y*y);returnSum;}doubleT(doublea,doubleb,intm,doublec,doubled,intn)//复化梯形公式计算二重积分{doubleSum=0;//用于对小区间的计算结果进行求和doubleh;//x轴上[a,b]区间经过m等分后，每个区间的大小doublek;//y轴上[c,d]区间经过n等分后，每个区间的大小inti;//x轴的每个区间的序号intj;//y轴的每个区间的序号h=(b-a)/m;k=(d-c)/n;for(i=0;im;i++)//复化梯形计算公式{for(j=0;jn;j++){Sum=Sum+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);}}Sum=Sum*h*k/4;T_Rank++;T_Matrix[0][0]=Sum;//保存最初计算结果returnSum;}doubleT_Dich(doublea,doubleb,intm0,doublec,doubled,intn0,ints)//二分法{doubleSum=0;//二分后的计算结果doubletemp0=0;//二分前的计算结果doubletemp1=0;doubletemp2=0;doubletemp3=0;doubletemp4=0;doubletemp5=0;intm;//[a,b]进行m等分intn;//]c,d]进行n等分doubleh;//x轴上[a,b]区间经过m等分后，每个区间的大小doublek;//y轴上[c,d]区间经过n等分后，每个区间的大小inti;//x轴的每个区间的序号intj;//y轴的每个区间的序号if(s1||T_Ranks-1)//如果上一步的二分结果不存在，直接返回0return0;else//如果上一步的二分结果存在，将其赋值给temp0temp0=T_Matrix[s-1][0];m=m0*pow(2.0,s-1);//根据二分的阶数和初始m0，计算当前的mn=n0*pow(2.0,s-1);//根据二分的阶数和初始n0，计算当前的nh=(b-a)/m;k=(d-c)/n;for(i=0;im;i++)temp1+=f(a+i*h+0.5*h,c)+f(a+i*h+0.5*h,d);for(j=0;jn;j++)temp2+=f(a,c+j*k+0.5*k)+f(b,c+j*k+0.5*k);for(i=0;im;i++){for(j=1;jn;j++)temp3+=f(a+i*h+0.5*h,c+j*k);}for(i=1;im;i++){for(j=0;jn;j++)temp4+=f(a+i*h,c+j*k+0.5*k);}for(i=0;im;i++){for(j=0;jn;j++)temp5+=f(a+i*h+0.5*h,c+j*k+0.5*k);}Sum=temp0*0.25+h*k*0.125*(temp1+temp2+2*(temp3+temp4+temp5));T_Rank++;T_Matrix[s][0]=Sum;//保存计算结果returnSum;}voidT_Extra(inti,intj)//外推法{T_Matrix[i][j]=T_Matrix[i+1][j-1]*pow(4.0,j)/(pow(4.0,j)-1.0)-T_Matrix[i][j-1]/(pow(4.0,j)-1.0);}doublecalc(doublea,doubleb,doublec,doubled,doubleeps)//基于二分法和外推法，根据误差，进行计算{//二重积分函数为f(x,y)，积分区间为a=x=b,c=y=d,误差为einti;intj;T_Rank=-1;for(i=0;i16;i++){for(j=0;j16;j++)T_Matrix[i][j]=0;}//对储存中间结果的矩阵进行清零intm=1;//将x轴的区间1等分intn=1;//将y轴的区间1等分doubleSum=0;T(a,b,m,c,d,n);//先计算最初的Tmn(f)T_Dich(a,b,m,c,d,n,T_Rank+1);//二分法中间值while(fabs(T_Matrix[0][T_Rank-1]-T_Matrix[1][T_Rank-1])eps){T_Dich(a,b,m,c,d,n,T_Rank+1);//二分法中间值T_Extra(0,T_Rank-1);//外推for(i=1;i=T_Rank-1;i++)T_Extra(T_Rank-i,i);//外推}returnT_Matrix[1][T_Rank-1];}int_tmain(intargc,_TCHAR*argv[]){doubles;//二重积分的结果doublea,b,c,d,eps;//x轴区间[a,b]，y轴区间[c,d]，误差epsclock_tstart,finish;//程序执行的开始时间和结束时间doublet;//程序执行时间chardisp[16];//printf(%f\n,pow(4.0,1));a=x0;b=xm;c=y0;d=yn;eps=eps0;start=clock();s=calc(a,b,c,d,eps);finish=clock();t=1.0*(finish-start)/CLOCKS_PER_SEC;printf(Thevaluesatisfiedis%f\n,s);printf(caculatetimeis%fseconds.\n,t);std::ofstreamoutfile(process.txt);//将程序执行结果输出outfile复化梯形+二分法+外推法计算二重积分。第0列为二分的结果，其他列均为外推的结果：\n;inti,j;for(i=0;i=T_Rank;i++)//向输出文件中，打印矩阵{for(j=0;jT_Rank;j++){if(jT_Rank-i)break;outfileT_Matrix[i][j]\t;}outfilestd::endl;}outfileThevaluesatisfiedissstd::endl;outfile计算时间为：tstd::endl;outfile.close();//关闭输出文件printf(Pressentertocontinue...);getchar();return0;}