mymath.h

数值计算方法中常用方法的程序实现

	assert(a.size()==b.size());
	assert((N+1)==c.size());
	
	std::size_t n = a.size();
	//构造系数矩阵
	std::vector<std::vector<double> > matrix;
	std::vector<double> temp(N+2,0);
	for(std::size_t i=0; i!=N+1; ++i)
		matrix.push_back(temp);
	//确定系数矩阵的值
	for(std::size_t i=0; i<=N; ++i)
		for(std::size_t j=0; j<=N; ++j)
		{
			double sum=0;
			for(std::size_t k=0; k!=n; ++k)
				sum+=b[k]*pow(a[k].getX(),(i+j)*1.0);
			matrix[i][j]=sum;
		}
	for(std::size_t i=0; i<=N; ++i)
	{
		double s=0;
		for(std::size_t k=0; k!=n; ++k)
			s+=b[k]*pow(a[k].getX(),i*1.0)*a[k].getY();
		matrix[i][N+1]=s;
	}
	//改进平方根法计算
	SqrtMethod(matrix, c);
}
//--------------------------------
//以下积分用到的函数，题目改变，只需要进入f改动函数原型就可以了
double f(double x)
{
	return pow(x,3.0)-2*x-5;
}
//以下程序中，lower表示积分下限，upper表示积分上限，n表示积分域被等分的段数，
//(*f)为一个指向函数f的指针，该指针含有一个double形参，并返回一个double值
//复化梯形求积公式求积分
double TrapeziumIntegral(double lower, double upper, std::size_t n, double (*f)(double))
{
	double h=(upper-lower)/n;
	double sum=((*f)(lower)+(*f)(upper))/2;  
	double s=0;
	for(std::size_t i=1; i<n; ++i)
		s+=(*f)(lower+i*(upper-lower)/n);
	sum+=s;
	return sum*h;	
}
//-----------------------------------
//复化辛普森求积公式求积分
double SimpsonIntergral(double lower, double upper, std::size_t n, double (*f)(double))
{
	double h=(upper-lower)/2/n;
	double sum=(*f)(lower)+(*f)(upper);   
	double s1=0;
	double s2=0;
	for(std::size_t i=1; i<=n; ++i)
		s1+=(*f)(lower+(2*i-1)*(upper-lower)/2/n);
	for(std::size_t i=1; i<n; ++i)
		s2+=(*f)(lower+i*(upper-lower)/n);
	sum+=4*s1+2*s2;
	return sum*h/3;
}
//-----------------------------------
//龙贝格积分
double RombergIntergral(double lower, double upper)
{
	std::vector<double> T;
	std::vector<double> S;
	std::vector<double> C;
	std::vector<double> R;
	std::size_t k=0;
	while(1)
	{
		T.push_back(TrapeziumIntegral(lower,upper,(unsigned)(pow(2.0,k*1.0)),f));
		if(k>=1) S.push_back(T[k]*4/3-T[k-1]/3);
		if(k>=2) C.push_back(S[k-1]*16/15-S[k-2]/15);
		if(k>=3) R.push_back(C[k-2]*64/63-C[k-3]/63);
		if(k>3&&fabs(TrapeziumIntegral(lower,upper,(unsigned)(pow(2.0,k*1.0))+1,f)-T[k])<EPS) break;
		k++;
	}
	return R[R.size()-1];
}
//------------------------------------
//不等距点上给定函数的数值积分
void PointIntgrl(const std::vector<double>& X, const std::vector<double>& Y, std::vector<double>& SUM)
{
	std::size_t n=X.size();
	assert(X.size()==Y.size());
	
	std::vector<double> S(n,0);
	std::vector<double> A(S);
	std::vector<double> B(S);
	std::vector<double> C(S);
	std::vector<double> F(S);
	std::vector<double> W(S);
	std::vector<double> SB(S);
	std::vector<double> G(S);
	std::vector<double> EM(S);
	
	for(std::size_t i=1; i!=n; ++i)
		S[i]=X[i]-X[i-1];
	
	for(std::size_t i=1; i!=n-1; ++i)
	{
		A[i]=S[i]/6;
		B[i]=(S[i]+S[i+1])/3;
		C[i]=S[i+1]/6;
		F[i]=(Y[i+1]-Y[i])/S[i+1]-(Y[i]-Y[i-1])/S[i];
	}
	
	A[n-1]=-0.5;
	B[0]=1.0;
	B[n-1]=1.0;
	C[0]=-0.5;
	F[0]=0;
	F[n-1]=0;
	W[0]=B[0];
	SB[0]=C[0]/W[0];
	G[0]=0;
	
	for(std::size_t i=1; i!=n; ++i)
	{
		W[i]=B[i]-A[i]*SB[i-1];
		SB[i]=C[i]/W[i];
		G[i]=(F[i]-A[i]*G[i-1])/W[i];
	}
	
	EM[n-1]=G[n-1];
	for(std::size_t i=1; i!=n; ++i)
		EM[n-i-1]=G[n-i-1]-SB[n-i-1]*EM[n-i];
	
	SUM[0]=0;
	for(std::size_t i=1; i!=n; ++i)
		SUM[i]=SUM[i-1]+S[i]*(Y[i]+Y[i-1])/2-pow(S[i],3.0)*(EM[i]+EM[i-1])/24;
}
//------------------------------------------
void SPLINT(const std::vector<double>& X, const std::vector<double>& Y, const std::vector<double>& Z, std::vector<double>& YINT)
{
	std::size_t n=X.size();
	assert(X.size()==Y.size());
	std::size_t m=Z.size();
	
	std::vector<double> S(n,0);
	std::vector<double> A(S);
	std::vector<double> B(S);
	std::vector<double> C(S);
	std::vector<double> F(S);
	std::vector<double> W(S);
	std::vector<double> SB(S);
	std::vector<double> G(S);
	std::vector<double> EM(S);
	for(std::size_t i=1; i!=n; ++i)
		S[i]=X[i]-X[i-1]; 
	
	for(std::size_t i=1; i!=n-1; ++i)
	{
		A[i]=S[i]/6;
		B[i]=(S[i]+S[i+1])/3;
		C[i]=S[i+1]/6;
		F[i]=(Y[i+1]-Y[i])/S[i+1]-(Y[i]-Y[i-1])/S[i];
	}
	
	A[n-1]=-0.5;
	B[0]=1.0;
	B[n-1]=1.0;
	C[0]=-0.5;
	F[0]=0;
	F[n-1]=0;
	W[0]=B[0];
	SB[0]=C[0]/W[0];
	G[0]=0;
	
	for(std::size_t i=1; i!=n; ++i)
	{
		W[i]=B[i]-A[i]*SB[i-1];
		SB[i]=C[i]/W[i];
		G[i]=(F[i]-A[i]*G[i-1])/W[i];
	}
	
	EM[n-1]=G[n-1];
	for(std::size_t i=1; i!=n; ++i)
		EM[n-i-1]=G[n-i-1]-SB[n-i-1]*EM[n-i];
		
	for(std::size_t i=0;i!=m;++i)
	{
		std::size_t k=1;
		if(Z[i]-X[0]<0)
		{
			if(Z[i]<(1.1*X[0]-0.1*X[1]))
				YINT[i]=EM[k-1]*pow((X[k]-Z[i]),3.0)/6.0/S[k]+EM[k]*pow(Z[i]-X[k-1],3.0)/6.0/S[k]+(Y[k]/S[k]-EM[i]*S[i]/6.0)*(Z[i]-X[k-1])
				+(Y[k-1]/S[k]-EM[k-1]*S[k]/6.0)*(X[k]-Z[i]);
		}
		
		else if(Z[i]==X[0])	YINT[i]=Y[0];
		
		else if(Z[i]-X[0]>0)
		{
		//	bool sign=false;
			while(1)
			{
				if(Z[i]-X[k]<0)
				{
					YINT[i]=EM[k-1]*pow((X[k]-Z[i]),3.0)/6.0/S[k]+EM[k]
							*pow(Z[i]-X[k-1],3.0)/6.0/S[k]+(Y[k]/S[k]-EM[i]
							*S[i]/6.0)*(Z[i]-X[k-1])
							+(Y[k-1]/S[k]-EM[k-1]*S[k]/6.0)*(X[k]-Z[i]);
					break;
				}
				else if(Z[i]-X[k]==0) { YINT[i]=Y[k]; break; }
				else if(Z[i]-X[k]>0)
				{
					k++;
					if(k>n&&Z[i]>(1.1*X[n-1]-0.1*X[n-2]))
					{
						k=n;
						YINT[i]=EM[k-1]*pow((X[k]-Z[i]),3.0)/6.0/S[k]+EM[k]
								*pow(Z[i]-X[k-1],3.0)/6.0/S[k]+(Y[k]/S[k]-EM[i]
								*S[i]/6.0)*(Z[i]-X[k-1])
								+(Y[k-1]/S[k]-EM[k-1]*S[k]/6.0)*(X[k]-Z[i]);
					}
				}
			}
		}
	}
}
//----------------------------------------------
//方程求根
//简化牛顿下山法
//参数为一个函数指针来定义所求函数，x为函数的自变量，
//c为系数，通常为x0点的导数，m为松弛因子，如不收敛，可以换为一半
double SNTon(double (*f)(double), double x, double c, double m=1)
{
	while(1)
	{
		double X=x-m*(*f)(x)/c;
		if(fabs(X-x)<EPS) break;
		else x=X;
	}
	return x;
}
//-----------------------------------
//弦割法(2点)
double XGMethod(double (*f)(double), double x0, double x1)
{
	while(1)
	{
		double X=x1-(*f)(x1)*(x1-x0)/((*f)(x1)-(*f)(x0));
		if(fabs(X-x1)<EPS) break;
		else { x0=x1; x1=X; }
	}
	return x1;
}
//-----------------------------------
//线性方程组的迭代解法，适用于稀疏矩阵
//雅可比迭代法
// a为增广矩阵，x为所求
void Jacobi(const std::vector<std::vector<double> >& a, std::vector<double>& x)
{
	assert(a.size()==x.size());
	assert(a[0].size()==a.size()+1);
	
	std::size_t n=a.size();
	std::vector<double> s(n,0);
	std::vector<double> X(n,0);
	
	while(1)
	{
		double D=0;
		for(std::size_t i=0; i!=n; ++i)
			for(std::size_t j=0; j!=n; ++j)
				if(i!=j) s[i]+=a[i][j]*x[j];
		for(std::size_t i=0; i!=n; ++i)
		{
			X[i]=(a[i][n]-s[i])/a[i][i];
			D+=fabs(X[i]-x[i]);
		}
		if(D<EPS) break;
		else { x.swap(X); s.assign(n,0); }
	}
}
//-------------------------------------
//SOR迭代法解线性方程组
//m控制收敛速度的松弛因子，当m为1时，为塞德尔迭代
void SOR(const std::vector<std::vector<double> >& a, std::vector<double>& x, double m=1)
{
	assert(a.size()==x.size());
	assert(a[0].size()==a.size()+1);
	
	std::size_t n=a.size();
	std::vector<double> s(n,0);
	
	while(1)
	{
		double D=0;
		std::vector<double> X(x);
		for(std::size_t i=0; i!=n; ++i)
		{
			for(std::size_t j=0; j!=n; ++j)
				if(i!=j) s[i]+=a[i][j]*x[j];
			x[i]=m*(a[i][n]-s[i])/a[i][i]+(1-m)*x[i];
			D+=fabs(X[i]-x[i]);
		}
		if(D<EPS) break;
		else s.assign(n,0);
	}
}
//------------------------------------
//矩阵相乘
//c只能为所有元素全为零的向量
void MultiplyMatrix(const std::vector<std::vector<double> >& a,
					const std::vector<std::vector<double> >& b,
					std::vector<std::vector<double> >& c)
{
	std::size_t m1=a.size();
	std::size_t n1=a[0].size();
	std::size_t m2=b.size();
	std::size_t n2=b[0].size();
	//参数检查
	assert(n1==m2);
	assert(m1==c.size());
	assert(n2==c[0].size());
	for(std::size_t i=0; i!=m1; ++i)
		for(std::size_t j=0; j!=n2; ++j)
			assert(c[i][j]==0);
	//计算主体
	for(std::size_t i=0; i!=m1; ++i)
		for(std::size_t j=0; j!=n2; ++j)
			for(std::size_t k=0; k!=n1; k++)
			c[i][j]+=a[i][k]*b[k][j];
}
//-----------------------------------
//向量中模最大的最大元素
double MaxElem(const std::vector<double>& a)
{
	std::size_t n=a.size();
	double temp=a[0];
	for(std::size_t j=0; j!=n; ++j)
		if(fabs(temp)<fabs(a[j])) temp=a[j];
	return temp;
}
//矩阵中模最大的元素
double MaxElem(const std::vector<std::vector<double> >& a)
{
	std::size_t m=a.size();
	std::size_t n=a[0].size();
	double temp=a[0][0];
	for(std::size_t i=0; i!=m; ++i)
		for(std::size_t j=0; j!=n; ++j)
			if(fabs(temp)<fabs(a[i][j])) temp=a[i][j];
	return temp;
}
//-----------------------------------
//乘幂法
//求最大的特征值及对应的特征向量
void PWMethod(const std::vector<std::vector<double> >& a)
{
	std::size_t m=a.size();
	std::size_t n=a[0].size();
	//初始化
	std::vector<double> temp1(1,1);
	std::vector<double> temp2(1,0);
	std::vector<std::vector<double> > z;
	std::vector<std::vector<double> > y;
	for(std::size_t i=0; i!=n; ++i)
		z.push_back(temp1);
	for(std::size_t i=0; i!=m; ++i)
		y.push_back(temp2);
	double s; //最大的特征值
	//循环开始
	while(1)
	{
		double D=0;
		std::vector<double> Z;
		for(std::size_t i=0; i!=m; ++i)
			Z.push_back(z[i][0]);
		MultiplyMatrix(a,z,y);
		s=MaxElem(y);
		for(std::size_t i=0; i!=m; ++i)
		{
			z[i][0]=y[i][0]/s;
			D+=fabs(Z[i]-z[i][0]);
		}
		if(D<EPS) break;
		else //别忘了y向量中的元素要置零
			for(std::size_t i=0; i!=m; ++i)
				y[i][0]=0;
	}
	std::cout<<std::setprecision(10);
	std::cout<<"模最大的特征值: "<<s<<std::endl;
	std::cout<<"对应的特征向量: "<<std::endl;
	for(std::size_t i=0; i!=m; ++i)
		std::cout<<z[i][0]<<std::endl;
}
//------------------------------------
//反幂法
//主要用于知道某特征值的近似值，求该值的精确值和对应的特征向量
void InversePWMethod(const std::vector<std::vector<double> >& a, double p)
{
	assert(a.size()==a[0].size());
	std::size_t n=a.size();
	std::vector<double> temp(n,0);
	std::vector<std::vector<double> > r;
	std::vector<std::vector<double> > l;
	std::vector<std::vector<double> > A;
	for(std::size_t i=0; i!=n; ++i)
	{
		r.push_back(temp);
		l.push_back(temp);
		A.push_back(temp);
	}
	//计算A-pI
	for(std::size_t i=0; i!=n; ++i)
		for(std::size_t j=0; j!=n; ++j)
		{
			if(i==j) A[i][j]=a[i][j]-p;
			else A[i][j]=a[i][j];
		}
	//LU分解
	for(std::size_t i=1; i<=n; ++i)
	{
		//计算R
		for(std::size_t j=i; j<=n; ++j)
		{
			double m1=0;
			for(std::size_t k=1; k<=i-1; ++k)
				m1+=l[i-1][k-1]*r[k-1][j-1];
			r[i-1][j-1]=A[i-1][j-1]-m1;
		}
		//计算l
		for(std::size_t j=i+1; j<=n; ++j)
		{
			double m3=0;
			for(std::size_t k=1; k<=i-1; ++k)
				m3+=l[j-1][k-1]*r[k-1][i-1];
			l[j-1][i-1]=(A[j-1][i-1]-m3)/r[i-1][i-1];
		}
		l[i-1][i-1]=1;  //L对角线元素都变为1
	}
	//循环主体
	std::vector<double> u(n,1);
	std::vector<double> y(n,0);