multitrapezia.cpp

常用数值算法的C++源码

#include "stdio.h"
#include "math.h"
#define  PI   3.141593
#define  SIZE 1025          //定义存储插值节点的数组大小

// 积分原函数
double f(double x,double y)
{
	double f = tan(pow(x,2)+pow(y,2));
	return f;
}

// 复化梯形公式计算二重积分
// x的积分区间为[a,b]
// y的积分区间为[c,d]
// 计算至相邻两次近似值之差的绝对值不超过err
double MultiTrapezia(double a,double b,double c,double d,double err)
{
	int    m    = 1;
	int    n    = 1;
    double h    = b-a;
	double k    = d-c;
	double sum1,sum2,sum3,sum4,sum5;
	double delta,result;
	double x1[SIZE],x2[SIZE],y1[SIZE],y2[SIZE],T1[SIZE],T2[SIZE],T3[SIZE],T4[SIZE];
	x1[0] = a;
	x1[1] = b;
	y1[0] = c;
	y1[1] = d;
	T1[0] = 0.25*h*k*(f(a,c)+f(a,d)+f(b,c)+f(b,d));

	for(int p=0;(2*m<SIZE)&&(2*n<SIZE);p++)
	{
		sum1 = 0;
	    sum2 = 0;
	    sum3 = 0;
	    sum4 = 0;
	    sum5 = 0;
		//计算插值节点
		for(int i=0;i<=2*m;i++) x2[i] = a+0.5*i*h;
		for(int j=0;j<=2*n;j++) y2[j] = c+0.5*j*k;
		//计算T_2m,2n
		for(i=0;i<m;i++) sum1 = sum1+f(x2[2*i+1],y1[0])+f(x2[2*i+1],y1[n]);
		for(j=0;j<n;j++) sum2 = sum2+f(x1[0],y2[2*j+1])+f(x1[m],y2[2*j+1]);
		for(i=0;i<m;i++)
			for(j=1;j<n;j++) sum3 = sum3+f(x2[2*i+1],y1[j]);
		for(i=1;i<m;i++)
			for(j=0;j<n;j++) sum4 = sum4+f(x1[i],y2[2*j+1]);
		for(i=0;i<m;i++)
			for(j=0;j<n;j++) sum5 = sum5+f(x2[2*i+1],y2[2*j+1]);
		T1[p+1] = 0.25*T1[p]+0.125*h*k*(sum1+sum2+2*(sum3+sum4+sum5));
        //判断是否达到精度要求
		delta = fabs(T1[p+1]-T1[p]);
		if(delta<=err) 
		{
			result = T1[p+1];
			break;
		}
		//计算T(1)_m,n
		T2[p] = 4*T1[p+1]/3-T1[p]/3;
		if(p-1>=0)
		{
			//判断是否达到精度要求
			delta = fabs(T2[p]-T2[p-1]);
			if(delta<=err) 
			{
				result = T2[p];
				break;
			}
            //计算T(2)_m,n
			T3[p-1] = 16*T2[p]/15-T2[p-1]/15;
			if(p-2>=0)
			{
				//判断是否达到精度要求
				delta = fabs(T3[p-1]-T3[p-2]);
				if(delta<=err)
				{
					result = T3[p-1];
					break;
				}
				//计算T(3)_m,n
				T4[p-2] = 64*T3[p-1]/63-T3[p-2]/63;
				if(p-3>=0)
				{
					//判断是否达到精度要求
					delta = fabs(T4[p-2]-T4[p-3]);
					if(delta<=err)
					{
						result = T4[p-2];
						break;
					}
				}
			}
		}
		//重新划分区间
		m = 2*m;
		n = 2*n;
		h = h/2;
		k = k/2;
		for(i=0;i<=m;i++) x1[i] = x2[i];
		for(j=0;j<=n;j++) y1[j] = y2[j];
	}
	if((2*m+1)>SIZE||(2*n+1)>SIZE)
	{
		printf("数据溢出\n");
		return NULL;
	}
	else return result;
}

void main()
{
	double result;
	result = MultiTrapezia(0,PI/3,0,PI/6,0.5*pow(10,-5));
	printf("积分结果为：%f\n",result);
}