mymath.h

自编数学函数调用库

//————————————————————————————————————————————
//本文件为自己所编的数学函数调用库，附加有屏蔽的主程序以说明调用方法；
#ifndef mymath_h
#define mymath_h
#include<math.h>
#include"Point.h"
//template<class double>
double f(double (*F)(double x),double x,double a,double b){
	return F(((b-a)*x+a+b)/2);
}
//template<class double>
double GL(double (*F)(double x),double a,double b,int n){//高斯勒让德积分算法
	double I;
	switch (n){
	case 0:
		I=f(F,0,a,b)*2.0;break;
	case 1:
		I=f(F,0.5773503,a,b)*1.0+f(F,-0.5773503,a,b)*1.0;
			break;	
	case 2:
		I=f(F,0.7745967,a,b)*0.5555556+f(F,-0.7745967,a,b)*0.5555556+f(F,0,a,b)*0.8888889;
			break;
	case 3:
		I=f(F,0.8611363,a,b)*0.3478548+f(F,-0.8611363,a,b)*0.3478548+
			f(F,0.3398810,a,b)*0.6521452+f(F,-0.3398810,a,b)*0.6521452;
			break;
	case 4:
		I=f(F,0.9061793,a,b)*0.2369269+f(F,-0.9061793,a,b)*0.2369269+
		f(F,0,a,b)*0.5688889+f(F,0.5384693,a,b)*0.4786287+f(F,-0.5384693,a,b)*0.4786287;
		break;
	}
	return (b-a)*I/2;
}

/*void main(){
	int i; 
	double a=2.5,b=8.4;
    //double a=-1,b=1;
	for(i=0;i<=4;i++)
		cout<<GL(Fun,a,b,i)<<endl;
	cout<<sqrt(15)/5;
}*/
double LBG(double (*f)(double ),double a,double b){//龙贝格积分算法
	int i,k,t,n,m,p;
	double h,temp=0;
	double A[15],B,C;
	double e=1.0e-7;
	k=0;
	h=b-a;
	A[0]=(f(a)+f(b))*h/2;
	C=1;
	while(C>e){
		B=A[0];
		k=k+1;
		h=h/2;
		p=(int)pow(2,k-1);
		for(i=1;i<=p;i++)
			temp+=f(a+(2*i-1)*h);
		A[k]=A[k-1]/2+h*temp;
		temp=0;
		t=k-1;
		m=4;
		n=1;
		do{
			A[t]=(m*A[t+1]-A[t])/(m-1);
			t=t-1;m=4*m;
			n++;
		}
		while(t+1);
		C=fabs(B-A[0]);
	}
	return A[0];
}
double LBG2(double (*f)(double,int ),int j,double a,double b){//多加一个参数的龙贝格积分算法
	int i,k,t,n,m,p;
	double h,temp=0;
	double A[15],B,C;
	double e=1.0e-7;
	k=0;
	h=b-a;
	A[0]=(f(a,j)+f(b,j))*h/2;
	C=1;
	while(C>e){
		B=A[0];
		k=k+1;
		h=h/2;
		p=(int)pow(2,k-1);
		for(i=1;i<=p;i++)
			temp+=f(a+(2*i-1)*h,j);
		A[k]=A[k-1]/2+h*temp;
		temp=0;
		t=k-1;
		m=4;
		n=1;
		do{
			A[t]=(m*A[t+1]-A[t])/(m-1);
			t=t-1;m=4*m;
			n++;
		}
		while(t+1);
		C=fabs(B-A[0]);
	}
	return A[0];
}
template <class Type>//LU分解法解方程组模板函数，计数从0开始
Type* LU(Type (*a)[50],Type* b,int n){
	int i,j,r,k;
	Type temp=0;
	Type* x=new Type[n];
	for(i=1;i<n;i++)//LU分解:
		a[i][0]=a[i][0]/a[0][0];
	for(r=1;r<n;r++){
		for(i=r;i<n;i++)
		{for(k=0;k<r;k++)
		temp+=a[r][k]*a[k][i];
		a[r][i]=a[r][i]-temp;
		temp=0;
		}	
		if(r!=n-1)
		for(i=r+1;i<n;i++)
		{for(k=0;k<r;k++)
		temp+=a[i][k]*a[k][r];
		a[i][r]=(a[i][r]-temp)/a[r][r];
		temp=0;
		}		
	}
	x[0]=b[0];//回代:
	for(i=1;i<n;i++){
		for(k=0;k<i;k++)temp+=a[i][k]*x[k];
		x[i]=b[i]-temp;
		temp=0;
	}
	x[n-1]=x[n-1]/a[n-1][n-1];
	for(i=n-2;i>=0;i--){
		for(k=i+1;k<n;k++)temp+=a[i][k]*x[k];
		x[i]=(x[i]-temp)/a[i][i];
		temp=0;
	}
	cout<<"矩阵L:"<<endl;//打印LU:
	for(i=0;i<n;i++)
	{for(j=0;j<i;j++)cout<<a[i][j]<<" ";cout<<1<<endl;}
	cout<<"矩阵U:"<<endl;
	for(i=0;i<n;i++)
	{for(k=0;k<i;k++)cout<<" ";
	 for(j=i;j<n;j++)cout<<a[i][j]<<" ";
	 cout<<endl;}/**/
	return x;
}
/*void main(){
	int i;
	double A[5][5]={3.24,-2.18,5.09,-2.37,1.21,
	                0.73,3.85,-6.23,4.80,-5.93,
					2.88,5.73,-7.02,-9.17,3.58,
					2.10,3.02,-0.78,3.85,-6.00,
					1.20,-4.13,6.48,0.00,-3.24};
	double b[5]={28.36,-36.00,24.48,-16.23,4.34};
	double *x;
	x=LU(A,b,5);
	for(i=0;i<5;i++)
		cout<<"x["<<i+1<<"]="<<x[i]<<endl;
}*/
template <class Type>//追赶法解方程组，计数从1开始
void Chase(Type x[],Type a[],Type b[],Type c[],Type* f,int n){
	int i;
	Type* d;d=new Type[n+1];
	Type* y;y=new Type[n+1];
	d[1]=c[1]/b[1];
	for(i=2;i<n;i++)
		d[i]=c[i]/(b[i]-a[i]*d[i-1]);
	y[1]=f[1]/b[1];
	for(i=2;i<=n;i++)
		y[i]=(f[i]-a[i]*y[i-1])/(b[i]-a[i]*d[i-1]);
	x[n]=y[n];
	for(i=n-1;i>=1;i--)
		x[i]=y[i]-d[i]*x[i+1];
	delete []d;
	delete []y;
}
/*void main(){
	int i;
	double x[4]={0,0,0,0};
	double a[4]={0,0,2,2};
	double b[4]={0,1,1,1};
	double c[4]={2,2,0,0};
	double f[4]={0,3,5,3};
    Chase( x, a, b, c,f,3);
	for(i=1;i<=3;i++)
		cout<<"x["<<i<<"]="<<x[i]<<endl;
}*/

template <class Type>//高斯塞德尔迭代法求解方程组，计数从0开始
void GS(Type (*a)[50],Type* b,Type x[],int n,int m){
	int i,j,k;
	Type temp1=0,temp2=0;
	for(k=0;k<m;k++){
		for(i=0;i<n;i++){
		for(j=0;j<i;j++)temp1+=a[i][j]*x[j];
		for(j=i+1;j<n;j++)temp2+=a[i][j]*x[j];
		x[i]=(b[i]-temp1-temp2)/a[i][i];
		temp1=temp2=0;
		}
	}
}
/*void main(){
	int i;
	/*double A[3][3]={8,-3,2,4,11,-1,6,3,12};
	double b[3]={20,33,36};
	double x[3]={0,0,0};
	double A[3][3]={6,-1,-1,-1,6,-1,-1,-1,6};
	double b[3]={11.33,32,42};
	double x[3]={0,0,0};
	GS(A,b,x,3,15);
	for(i=0;i<3;i++)
		cout<<"x["<<i+1<<"]="<<x[i]<<endl;
}*/

double f(Point p){
	return sqrt(4*p.x*p.x-p.y*p.y);
}
double Guass(double (*f)(Point x),Point i,Point j,Point k){//3点公式
	double s=((i.x-j.x)*(i.y-k.y)-(i.x-k.x)*(i.y-j.y))/2;
	Point p1,p2,p3;//下面这一条语句它是先执行后面的加号
	p1=0.666666666666667*i+j*0.166666666666667+k*0.166666666666667;
	p2=0.666666666666667*j+i*0.166666666666667+k*0.166666666666667;
	p3=0.666666666666667*k+i*0.166666666666667+j*0.166666666666667;
	return s*(f(p1)+f(p2)+f(p3))/3;
}
/*void main(){
    cout.precision(4);
	Point p1(0,0),p2(1,0),p3(1,1);
	double result;
	result=Guass(f,p1,p2,p3);
	cout<<result<<endl;
	cout<<(3.1415926535897932/3+sqrt(3)/2)/3;	
}*/
double f(Point p,int r,int l){
	return sqrt(4*p.x*p.x-p.y*p.y);
}
double Guass3(double (*f)(Point x,int ,int ),int r,int l,Point i,Point j,Point k){//3点公式
	double s=((i.x-j.x)*(i.y-k.y)-(i.x-k.x)*(i.y-j.y))/2;
	Point p1,p2,p3;//下面这一条语句它是先执行后面的加号
	p1=0.666666666666667*i+j*0.166666666666667+k*0.166666666666667;
	p2=0.666666666666667*j+i*0.166666666666667+k*0.166666666666667;
	p3=0.666666666666667*k+i*0.166666666666667+j*0.166666666666667;
	return s*(f(p1,r,l)+f(p2,r,l)+f(p3,r,l))/3;
}
double Guass4(double (*f)(Point x,int ,int ),int r,int l,Point i,Point j,Point k,Point m){
	return Guass3(f,r,l,i,j,k)+Guass3(f,r,l,k,m,i);
}
void main(){
    cout.precision(4);
	int r,l;
	Point p1(0,0),p2(1,0),p3(1,1),p4(0,1);
	double result;
	result=Guass4(f,r,l,p1,p2,p3,p4);
	cout<<result<<endl;
	//cout<<(3.1415926535897932/3+sqrt(3)/2)/3;
	cout<<double(1)/12;
}/**/
#endif