mymath.h

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

	std::vector<double> z(n,0);
	double s;
	while(1)
	{
		//计算ry=u
		y[n-1]=u[n-1]/r[n-1][n-1];
		for(int i=n-1; i!=0; i--)
		{	
			double m=0;
			for(std::size_t j=i+1;j<=n;j++)
				m+=r[i-1][j-1]*y[j-1];
			y[i-1]=(u[i-1]-m)/r[i-1][i-1];
		}
		std::vector<double> Z(z);
		//选最大元素
		s=MaxElem(y);
		for(std::size_t i=0; i!=n; ++i)
			z[i]=y[i]/s;
		//判断循环是否中止
		double D=0;
		for(std::size_t i=0; i!=n; ++i)
			D+=fabs(Z[i]-z[i]);
		if(D<EPS) break;
		else 
		{	//解lu=z
			u[0]=z[0];
			for(std::size_t i=1; i!=n; ++i)
			{	
				double m=0;
				for(std::size_t j=0; j!=i; ++j)
					m+=l[i][j]*u[j];
				u[i]=z[i]-m;
			}
		}
	}
	std::cout<<std::setprecision(10);
	std::cout<<"特征值: "<<p+1/s<<std::endl;
	std::cout<<"对应的特征向量: "<<std::endl;
	for(std::size_t i=0; i!=n; ++i)
		std::cout<<z[i]<<std::endl;
}
//-------------------------------------
//雅可比法
//用于求解实对称矩阵的特征值和特征向量
void JacobiMatrix(const std::vector<std::vector<double> >& a)
{
	//初始化
	std::vector<std::vector<double> > A(a); //多维向量也支持一维向量一样的拷贝
	std::size_t n=A.size();
	std::vector<double> temp(n,0);
	std::vector<std::vector<double> > I;   //临时矩阵，代表单位矩阵
	for(std::size_t i=0; i!=n; ++i)
		I.push_back(temp);
	for(std::size_t i=0; i!=n; ++i)
		I[i][i]=1;
	std::vector<std::vector<double> > P(I);
	//循环开始
	while(1)
	{
		//计算非对角线绝对值最大的元素及其矩阵中的坐标
		std::size_t ii=0;
		std::size_t jj=1;
		double MAX=0;
		for(std::size_t i=0; i!=n-1; ++i)
		{
			for(std::size_t j=i+1; j!=n; ++j)
				if(fabs(MAX)<fabs(A[i][j])) 
				{
					MAX=A[i][j];
					ii=i;
					jj=j;
				}
		}
		//循环中止条件
		if(fabs(MAX)<EPS) break;
		//计算cos和sin
		double c;  //cos
		double s;  //sin
		if(A[ii][ii]==A[jj][jj])
		{
			if(MAX>0)
			{
				c=pow(2.0,0.5)/2;
				s=c;
			}
			if(MAX<0)
			{
				s=pow(2.0,0.5)/2;
				c=-s;
			}
		}
		else 
		{
			c=cos(atan(2*MAX/(A[ii][ii]-A[jj][jj]))/2);
			s=sin(atan(2*MAX/(A[ii][ii]-A[jj][jj]))/2);
		}
		//计算A=P1AP2和P
		std::vector<std::vector<double> > A1;  //临时矩阵
		std::vector<std::vector<double> > A2;  //临时矩阵
		std::vector<std::vector<double> > Ptemp;  //临时矩阵
		for(std::size_t i=0; i!=n; ++i)
		{
			A1.push_back(temp);
			A2.push_back(temp);
			Ptemp.push_back(temp);
		}
		std::vector<std::vector<double> > P1(I);   
		std::vector<std::vector<double> > P2(I);   
		P1[ii][ii]=c;
		P1[jj][jj]=c;
		P1[ii][jj]=s;
		P1[jj][ii]=-s;
		P2[ii][ii]=c;
		P2[jj][jj]=c;
		P2[ii][jj]=-s;
		P2[jj][ii]=s;
		MultiplyMatrix(P1,A,A1);
		MultiplyMatrix(A1,P2,A2);
		MultiplyMatrix(P,P2,Ptemp);
		for(std::size_t i=0; i!=n; ++i)
			for(std::size_t j=0; j!=n; ++j)
			{
				A[i][j]=A2[i][j];
				P[i][j]=Ptemp[i][j];
			}
	}
	//打印结果
	std::cout<<std::fixed;
	std::cout<<"特征值: "<<std::endl;
	for(std::size_t i=0; i!=n; ++i)
		std::cout<<A[i][i]<<std::endl;
	std::cout<<"特征向量: "<<std::endl;
	for(std::size_t i=0; i!=n; ++i)
	{
		for(std::size_t j=0; j!=n; ++j)
			std::cout<<P[i][j]<<'\t';
		std::cout<<std::endl;
	}
}
//---------------------------------------
//符号判断
int sign(double x)
{
    if(x>0) return 1;
    else if(x<0) return -1;
    else return 0;
}
//---------------------------------------
//镜射矩阵的转化(方阵)
//a为原始矩阵，A为转化后的矩阵
void HouseHolder(const std::vector<std::vector<double> >& a, std::vector<std::vector<double> >& b)
{
	assert(a.size()==a[0].size());
	assert(b.size()==b[0].size());
	assert(a.size()==a.size());
	
	std::size_t n=a.size();
	std::vector<std::vector<double> > A(a);
	//循环
	for(std::size_t j=0; j!=n-2; ++j)
	{
		//计算c和rou
		double c=0;
		double rou=0;
		for(std::size_t i=j+1; i!=n; ++i)
			c+=A[i][j]*A[i][j];
		c=-sign(A[j+1][j])*pow(c,0.5);
		rou=pow(2*c*(c-A[j+1][j]),0.5);
		//计算u
		std::vector<std::vector<double> >u;
		std::vector<double> temp1(1,0);
		for(std::size_t k=0; k!=n; ++k)
			u.push_back(temp1);
		u[j+1][0]=(A[j+1][j]-c)/rou;
		for(std::size_t k=j+2; k!=n; ++k)
			u[k][0]=A[k][j]/rou;
		//计算uT
		std::vector<std::vector<double> > uT;
		std::vector<double> temp2;
		for(std::size_t k=0; k!=n; ++k)
			temp2.push_back(u[k][0]);
		uT.push_back(temp2);
		//计算H及HT    //H=I-2*u*uT
		std::vector<std::vector<double> > H;
		std::vector<std::vector<double> > HT;
		std::vector<std::vector<double> > uuT; //u*uT
		std::vector<double> temp3(n,0);
		for(std::size_t k=0; k!=n; ++k)
		{
			uuT.push_back(temp3);
			H.push_back(temp3);
			HT.push_back(temp3);
		}
		std::vector<std::vector<double> > I;   //单位矩阵
		for(std::size_t i=0; i!=n; ++i)
			I.push_back(temp3);
		for(std::size_t i=0; i!=n; ++i)
			I[i][i]=1;
		MultiplyMatrix(u,uT,uuT);
		for(std::size_t k=0; k!=n; ++k)
			for(std::size_t i=0; i!=n; ++i)
			{
				H[k][i]=I[k][i]-uuT[k][i]*2;
				HT[i][k]=H[k][i];
			}
		//计算矩阵b
		std::vector<std::vector<double> > A1;  //临时矩阵
		for(std::size_t i=0; i!=n; ++i)
			A1.push_back(temp3);
		MultiplyMatrix(H,A,A1);
		MultiplyMatrix(A1,HT,b);
		//如果循环没有完毕，则要把b中元素拷贝到A中，然后b中元素要清零
		if(j==n-3) break;
		else
			for(std::size_t i=0; i!=n; ++i)
				for(std::size_t k=0; k!=n; ++k)
				{
					A[i][k]=b[i][k];
					b[i][k]=0;
				}
	}
}
//------------------------------------------
//拟上三角阵(海森伯格矩阵)的QR算法
void QR(const std::vector<std::vector<double> >& a)
{
	assert(a.size()==a[0].size());
	std::size_t n=a.size();
	std::vector<std::vector<double> > A(a);
	while(1)
	{
		double u=A[n-1][n-1];  //位移量u
		for(std::size_t i=0; i!=n; ++i)
			A[i][i]=A[i][i]-u;;         //
		
		std::vector<double> Cta;
		std::vector<double> s;
		std::vector<double> c;
		for(std::size_t i=0; i!=n-1; ++i)
		{
			Cta.push_back(atan(A[i+1][i]/A[i][i]));
			s.push_back(sin(Cta[i]));
			c.push_back(cos(Cta[i]));
		//解法一：
			//左乘旋转矩阵P
			for(std::size_t j=i; j!=n; ++j)
			{
				double z=A[i][j]*c[i]+A[i+1][j]*s[i];
				A[i+1][j]=-A[i][j]*s[i]+A[i+1][j]*c[i];
				A[i][j]=z;
			}
		}
		//右乘旋转矩阵P
		for(std::size_t i=0; i!=n-1; ++i)
			for(std::size_t j=0; j<=i+1; ++j)
			{
				double z=A[j][i]*c[i]+A[j][i+1]*s[i];
				A[j][i+1]=-A[j][i]*s[i]+A[j][i+1]*c[i];
				A[j][i]=z;
			}
		
		//解法二：
			//左乘旋转矩阵P
		/*	std::vector<double> temp(n,0);
			std::vector<std::vector<double> > A1;   //临时矩阵
			std::vector<std::vector<double> > P1;   //临时矩阵
			for(std::size_t j=0; j!=n; ++j)
			{
				P1.push_back(temp);
				A1.push_back(temp);
			}
			for(std::size_t j=0; j!=n; ++j)
				P1[j][j]=1;
			P1[i][i]=c[i];
			P1[i+1][i+1]=c[i];
			P1[i][i+1]=s[i];
			P1[i+1][i]=-s[i];
			MultiplyMatrix(P1,A,A1);
			for(std::size_t j=0; j!=n; ++j)
				for(std::size_t k=0; k!=n; ++k)
					A[j][k]=A1[j][k];
		}
		//右乘旋转矩阵P
		for(std::size_t i=0; i!=n-1; ++i)
		{
			std::vector<double> temp(n,0);
			std::vector<std::vector<double> > A1;   //临时矩阵
			std::vector<std::vector<double> > P1;   //临时矩阵
			for(std::size_t j=0; j!=n; ++j)
			{
				P1.push_back(temp);
				A1.push_back(temp);
			}
			for(std::size_t j=0; j!=n; ++j)
				P1[j][j]=1;
			P1[i][i]=c[i];
			P1[i+1][i+1]=c[i];
			P1[i][i+1]=-s[i];
			P1[i+1][i]=s[i];
			MultiplyMatrix(A,P1,A1);
			for(std::size_t j=0; j!=n; ++j)
				for(std::size_t k=0; k!=n; ++k)
					A[j][k]=A1[j][k];
		}*/
		//------------------
		for(std::size_t i=0; i!=n; ++i)
			A[i][i]+=u;    //位移量u
		//判断收敛条件
		double D=0;
		for(std::size_t i=0; i!=n; ++i)
			for(std::size_t j=0; j!=n; ++j)
				if(i!=j) D+=fabs(A[i][j]);
		if(D<1) break;
		
		for(size_t i=0; i!=3; ++i)
		{	
			for(size_t j=0; j!=3; ++j)
				std::cout<<A[i][j]<<' ';
			std::cout<<std::endl;
		}
		system("pause");
	}
	//应用反幂法求特征值对应的特征向量
	for(std::size_t i=0; i!=n; ++i)
		InversePWMethod(a,A[i][i]);
}
//-------------------------------------------
//常微分问题的数值解法
//函数原型y'=f(x,y)
double Func(double x, double y)
{
	assert(x!=-0.5);
	return -0.9*y/(1+2*x);
}
//-------------------------------------------
//简单迭代法
//low下限，high上限，h步长,y0初值，Func函数指针
//当步长为上下限之差时，即为求一个点处的迭代值
void SimpleIterative(double low, double high, double h, double y0, double (*Func)(double,double))
{
	assert(high>low);
	std::size_t n=static_cast<unsigned>((high-low)/h);
	//欧拉法计算初值
	std::vector<double> y;
	y.push_back(y0);
	for(std::size_t i=1; i<=n; ++i) 
		y.push_back(y[i-1]+h*(*Func)(low+i*h,y[i-1]));
	//迭代
	while(1)
	{
		std::vector<double> temp(n+1,0);
		temp[0]=y0;
		for(std::size_t i=0; i<=n; ++i)
			temp[i+1]=temp[i]+h*((*Func)(low+i*h,y[i])+(*Func)(low+(i+1)*h,y[i+1]))/2;
		double s=0;
		for(std::size_t i=0; i<=n; ++i)
			s+=fabs(temp[i]-y[i]);
		if(s<EPS) break;
		else 
			for(std::size_t i=0; i<=n; ++i)
				y[i]=temp[i];
	}
	std::cout<<std::setprecision(10);
	for(std::size_t i=0; i<=n; ++i)
		std::cout<<y[i]<<std::endl;
}
//-----------------------------------------------
//龙格-库塔法(四级四阶)
void R_K(double low, double high, double h, double y0, double (*Func)(double,double))
{
	assert(high>low);
	std::size_t n=static_cast<unsigned>((high-low)/h);
	std::vector<double> y;
	y.push_back(y0);
	for(std::size_t i=0; i!=n; ++i)
	{
		double K1=h*(*Func)(low+i*h,y[i]);
		double K2=h*((*Func)(low+i*h+h/2,y[i]+K1/2));
		double K3=h*((*Func)(low+i*h+h/2,y[i]+K2/2));
		double K4=h*((*Func)(low+i*h+h,y[i]+K3));
		y.push_back(y[i]+(K1+2*K2+2*K3+K4)/6);
	}
	for(std::size_t i=0; i<=n; ++i)
		std::cout<<std::setprecision(10)<<y[i]<<std::endl;
}
//------------------------------------------------
//修正哈明预测-校正公式
void EmmendHamming(double low, double high, double h, 
					double y0, double y1, double y2, double y3, 
					double (*Func)(double,double))
{
	assert(high>low);
	std::size_t n=static_cast<unsigned>((high-low)/h);
	std::vector<double> y(n+1,0);
	y[0]=y0;
	y[1]=y1;
	y[2]=y2;;
	y[3]=y3;
	std::vector<double> f(n+1,0);
	f[0]=(*Func)(low,y[0]);
	f[1]=(*Func)(low+h,y[1]);
	f[2]=(*Func)(low+2*h,y[2]);
	std::vector<double> p(n+1,0);
	std::vector<double> m(p);
	std::vector<double> c(p);
	for(std::size_t i=3; i!=n; ++i)
	{
		f[i]=(*Func)(low+i*h,y[i]);
		p[i+1]=y[i-3]+4*h*(2*f[i]-f[i-1]+2*f[i-2])/3;
		m[i+1]=p[i+1]+(c[i]-p[i])*112/121;
		c[i+1]=(9*y[i]-y[i-2])/8+h*((*Func)(low+(i+1)*h,m[i+1])+2*f[i]-f[i-1])*3/8;
		y[i+1]=c[i+1]-(c[i+1]-p[i+1])*9/121;
	}
	for(std::size_t i=0; i<=n; ++i)
		std::cout<<std::setprecision(10)<<y[i]<<std::endl;
}
//------------------------------------------
//文件读取 
int DateInput(std::string ss, std::vector<std::vector<double> >& v)
{
	std::ifstream in(ss.c_str());
	if(!in)
	{
		std::cerr<<"error: unable to open input file: "<<in<<std::endl;
		return -1;
	}
	for(std::string s; getline(in,s); )
	{
		double x;
		std::vector<double> a;
		for(std::istringstream sin(s); sin>>x; a.push_back(x));
		v.push_back(a);
	}
	return 1;
}
//--------------------------------------------

#endif