integralellipticek.cpp

用C++或C平台实现的计算源码

		xt=&GaussPoints12[0];
		ceo=&GaussPtCeof12[0];
		break;
	case 13:
		xt=&GaussPoints13[0];
		ceo=&GaussPtCeof13[0];
		break;
	case 14:
		xt=&GaussPoints14[0];
		ceo=&GaussPtCeof14[0];
		break;
	case 15:
		xt=&GaussPoints15[0];
		ceo=&GaussPtCeof15[0];
		break;
	default:
		std::cout<<" number of Gauss points must be between 2 and 15!! "<<std::endl;
		return 0;
		break;
	}
	preVal=0;
	m=1;
	while (true) {
		val=0.0; 
		dx=(b-a)/m;
		for(j=0;j<m;j++){
			a0=a+j*dx;
			b0=a+(j+1)*dx;
			valtmp=0;
			for(i=0;i<nGauss;i++) 
			{ 
				normxt=(xt[i]*(b0-a0)+a0+b0)/2.0; 
				valtmp+=ceo[i]*fun(x0,normxt); 
			}
			val+=valtmp;
		}
		val=val*dx/2;
		if (fabs(val-preVal)/fabs(val)<eps) break;
		preVal=val;
		++m;
		if(m>mIter) {
			std::cout<<" !!! number of iter have already more than"<<mIter<<std::endl;
			break;
		}
	}
	return val;
}

double IntegralGaussLegendre2D(double(*fun)(double,double),double xa,double xb,double ya,double yb,double eps=1e-7,int nGauss=10) 
{ 
	double preVal,val,valtmp,normxt;
	double *xt,*ceo;
	int i,j,m;
	double a0,b0,dx;
	const int mIter=10;

	switch(nGauss) {
	case 2:
		xt=&GaussPoints2[0];
		ceo=&GaussPtCeof2[0];
		break;
	case 3:
		xt=&GaussPoints3[0];
		ceo=&GaussPtCeof3[0];
		break;
	case 4:
		xt=&GaussPoints4[0];
		ceo=&GaussPtCeof4[0];
		break;
	case 5:
		xt=&GaussPoints5[0];
		ceo=&GaussPtCeof5[0];
		break;
	case 6:
		xt=&GaussPoints6[0];
		ceo=&GaussPtCeof6[0];
		break;
	case 7:
		xt=&GaussPoints7[0];
		ceo=&GaussPtCeof7[0];
		break;
	case 8:
		xt=&GaussPoints8[0];
		ceo=&GaussPtCeof8[0];
		break;
	case 9:
		xt=&GaussPoints9[0];
		ceo=&GaussPtCeof9[0];
		break;
	case 10:
		xt=&GaussPoints10[0];
		ceo=&GaussPtCeof10[0];
		break;
	case 11:
		xt=&GaussPoints11[0];
		ceo=&GaussPtCeof11[0];
		break;
	case 12:
		xt=&GaussPoints12[0];
		ceo=&GaussPtCeof12[0];
		break;
	case 13:
		xt=&GaussPoints13[0];
		ceo=&GaussPtCeof13[0];
		break;
	case 14:
		xt=&GaussPoints14[0];
		ceo=&GaussPtCeof14[0];
		break;
	case 15:
		xt=&GaussPoints15[0];
		ceo=&GaussPtCeof15[0];
		break;
	default:
		std::cout<<" number of Gauss points must be between 2 and 15!! "<<std::endl;
		return 0;
		break;
	}
	preVal=0;
	m=1;
	while (true) {
		val=0.0; 
		dx=(xb-xa)/m;
		for(j=0;j<m;j++){
			a0=xa+j*dx;
			b0=xa+(j+1)*dx;
			valtmp=0;
			for(i=0;i<nGauss;i++) 
			{ 
				normxt=(xt[i]*(b0-a0)+a0+b0)/2.0; 
				valtmp+=ceo[i]*IntegralGaussLegendre2D_1(fun,normxt,ya,yb); 
			}
			val+=valtmp;
		}
		val=val*dx/2;
		if (fabs(val-preVal)/fabs(val)<eps) break;
		preVal=val;
		++m;
		if(m>mIter) {
			std::cout<<" !!! number of iter have already more than"<<mIter<<std::endl;
			break;
		}
	}
	return val;
}

double IntegralSimpsonStandard(double(*fun)(double),double a,double b,double eps=1e-5) 
{ 
	double preval,val,valtmp;
	int i,k,n;
	double dx,xa,xh,xb;
	const int mIter=5;

	k=0;
	preval=0;
	while (true) {
		n=pow(2,k);
		dx=(b-a)/n;
		valtmp=0;
		for(i=0;i<n;i++){
			xa=a+i*dx;
			xh=xa+0.5*dx;
			xb=xa+dx;
			valtmp+=fun(xa)+4*fun(xh)+fun(xb);
		}
		val=valtmp*dx/6;
		if (fabs(val-preval)/fabs(val)<eps) break;
		++k;
		preval=val;
		if (k>mIter) {
			std::cout<<" !!! number of iter have already more than"<<mIter<<std::endl;
			break;
		}		
	}
	return val;
}

double IntegralEllipticFirstK(double k,double eps=1e-5,int algo=3) 
{ 
	if(fabs(k)>=1){
		std::cout<<k<<" must be less than 1"<<std::endl;
		return 0;
	}
	int i,j,n;
	double sn,valtmp,tmpold,tmpeven,preval;

	if (algo==1) {
		n=1;
		preval=0.5*PI;
		while (true) {
			sn=0;
			for(i=1;i<=n;i++){
				tmpold=1;	tmpeven=1;
				for(j=1;j<=i;j++){
					tmpold*=(2*j-1);
					tmpeven*=2*j;
				}
				valtmp=tmpold/tmpeven;
				sn+=powf(k,2*i)*powf(valtmp,2);
			}
			sn=PI*0.5*(1+sn);
			if (fabs(preval-sn)/fabs(sn)<eps) break;
			preval=sn;
			++n;
		}
		return sn;
	}
	double a0,b0,c0,ai,bi;
	if (algo==2) {
		preval=PI/2;
		a0=1;
		b0=sqrt(1-k*k);
		c0=k;
		for(i=1;;i++){
			ai=(a0+b0)/2;
			bi=sqrt(a0*b0);
			sn=PI*0.5/ai;
			if (fabs(sn-preval)/fabs(sn)<eps) break;
			preval=sn;
			a0=ai;
			b0=bi;
		}
		return sn;
	}
	double kk;
	if (algo==3) {
		if ((k+1e-5)>=1) {
			return log(4./sqrt(1-k*k));
		}
		preval=PI/2;
		valtmp=1;
		kk=k;
		for(i=1;;i++){
			kk=sqrt(1-kk*kk);
			kk=(1-kk)/(1+kk);
			valtmp*=(1+kk);
			if (kk<1e-5) {
				valtmp*=PI/2.;
				break;
			}
		}
		return valtmp;
	}
	return 0;
}

double IntegralEllipticSecondE(double k,double eps=1e-5,int algo=3) 
{ 
	if(fabs(k)>=1){
		std::cout<<k<<" must be less than 1"<<std::endl;
		return 0;
	}
	int i,j,n;
	double sn,valtmp,tmpold,tmpeven,preval;

	if (algo==1) {
		n=1;
		preval=0.5*PI;
		while (true) {
			sn=0;
			for(i=1;i<=n;i++){
				tmpold=1;	tmpeven=1;
				for(j=1;j<=i;j++){
					tmpold*=(2*j-1);
					tmpeven*=2*j;
				}
				valtmp=tmpold/tmpeven;
				sn+=powf(k,2*i)*powf(valtmp,2)/(2*i-1);
			}
			sn=PI*0.5*(1-sn);
			if (fabs(preval-sn)/fabs(sn)<eps) break;
			preval=sn;
			++n;
		}
		return sn;
	}
	double a0,b0,c0,ai,bi,ci,citmp,K,E;
	if (algo==2) {
		preval=PI/2;
		a0=1;
		b0=sqrt(1-k*k);
		c0=k;
		citmp=0;
		for(i=1;;i++){
			ai=(a0+b0)/2;
			bi=sqrt(a0*b0);
			ci=(a0-b0)/2;
			K=PI*0.5/ai;
			sn=K*(1-0.5*(k*k+citmp));
			if (fabs(sn-preval)/fabs(sn)<eps) break;
			preval=sn;
			for(j=1;j<=i;j++){
				citmp+=powf(2,j)*powf(ci,2);
			}
			a0=ai;
			b0=bi;
			c0=ci;
		}
		return sn;
	}
	double kk,kk1;
	if (algo==3) {
		if ((k+1e-5)>=1) {
			return log(4./sqrt(1-k*k));
		}
		preval=PI/2;
		valtmp=1;
		kk=k;
		for(i=1;;i++){
			kk1=sqrt(1-kk*kk);
			kk=(1-kk1)/(1+kk1);
			if (kk<1e-5) {
				K=E=0.5*PI;
				valtmp=(E+kk1*K)/(1+kk1);
				break;
			}
		}
		return valtmp;
	}
	return 0;
}


double IntegralSimpsonQuick(double(*fun)(double),
				double a,double b,double eps=1e-5) 
{ 
	double valtmp;
	int i,k,n;
	double dx,t1,t2,s1,s2,x;
	const int mIter=5;

	k=0;
	n=1;
	dx=(b-a);
	t1=dx*(fun(a)+fun(b))/2;
	s1=t1;
	while (true) {
		valtmp=0;
		for(i=0;i<n;i++){
			x=a+(i+0.5)*dx;
			valtmp+=fun(x);
		}
		t2=(t1+dx*valtmp)/2.;
		s2=(4*t2-t1)/3.;
		if (fabs(s2-s1)<eps) break;
		t1=t2;
		n=n+n;
		dx=dx/2.;
		++k;
		if (k>mIter) {
			std::cout<<" !!! number of iter have already more than"<<mIter<<std::endl;
			break;
		}		
	}	
	return s2;
}

double funtest(double x)
{
	//return x*x;
	return(log(9*x+0.5)/(cos(x)+1/x));
}

double funtest2D(double x,double y)
{
	//return x*x;
	return(x*x+y*y);
}

void main()
{

	FILE *pfile;
	pfile=fopen("data1.txt","w+");
	double a,b,val,dk,pho,x1,y1,z1,k;
	int i;
	a=0;
	b=1;
	IntegralGaussLegendre2D(funtest2D,a,b,a,b);


	x1=10,y1=10,z1=10;
	a=1;
	pho=sqrt(x1*x1+y1*y1+z1*z1);
	k=sqrt(4*a*pho/(z1*z1+(a+pho)*(a+pho)));

	dk=1./100;
	for(i=0;i<100;i++)
	{
		k=i*dk;
		val=IntegralEllipticFirstK(k,1e-5,1);
		fprintf(pfile,"%10.5e\t%10.5e\n",k,val);
	}
	fclose(pfile);
	

/*
	a=0.3;
	b=1.2;
	val=IntegralGaussLegendre(funtest,a,b,1e-7,5);
	val=IntegralSimpsonQuick(funtest,a,b,1e-7);
	std::cout<<"result:"<<val<<endl;
*/
}