matrix.cpp

本文的重点在于通过一个实例来解决特殊的条形矩阵特征值得解算.

#include "Matrix.h"
#include "Max&Min.h"
//转化后的矩阵的初始化。
void InitMatrix(double A[M][N],double b,double c)  
{
	int i = 0;
	int j = 0;
	//输入矩阵第一行。
	A[i][j] = 0; j++;
	A[i][j] = 0; j++;
	for (j=2; j<N; j++)
	{
		A[i][j] = c;
	}
	//第一行输入完毕，下面输入第二行。
	i++; j = 0; 
	A[i][j] = 0; j++;
	for (j=1;j<N;j++)
	{
		A[i][j] = b;
	}
	//第二行输入完毕，下面输入第三行。
	i++; j = 0;
	for (j=0;j<N;j++)
	{
		A[i][j] = (1.64 - 0.024 * (j+1)) * sin(0.2 * (j+1)) - 0.64 * exp(0.1 / (j+1));
	}
	//第三行输入完毕，下面输入第四行。
    i++; j = 0;
	for (j=0;j<N-1;j++)
	{
		A[i][j] = b;
	}
	A[i][j] = 0;
	//第四行输入完毕，下面输入第五行。
	i++; j= 0;
	for (j=0;j<N-2;j++)
	{
		A[i][j] = c;
	}
	A[i][j] = 0; j++;
	A[i][j] = 0;
	//矩阵输入完毕

}

//求按模为最大的特征值。
void AbsMaxEigenval(double A[M][N],double &MaxEigenval)
{   
	
	int i,j,k;
	int r = 2;
	int s = 2; //定义上下带宽。
	int time = 10000; //设定迭代次数。
	double preci = 1e-12;  //设定迭代精度。
	double p;
	double temp;
	double Y[N] = {0};
	double U[N] = {0};
	
	//设定向量U的初值。
	for (i=0;i<N;i++)
	{
		U[i] = 1.0;
	}                       
	for (k=0;k<time;k++)
	{
		p = 0;
		for (i=0;i<N;i++)
		{
			p = U[i]*U[i] + p;    
		}
		for (i=0;i<N;i++)
		{
			Y[i] = U[i]/sqrt(p);
		}
		for (i=0;i<N;i++)
		{
			U[i] = 0;
			for (j = Max(i-r,0) ;j<=Min(i+s,N-1); j++)
			{
				U[i] = U[i] + A[i-j+s][j] * Y[j];
			}
		}
		temp = MaxEigenval;
		MaxEigenval = 0;
		for (i=0;i<N;i++)
		{
			MaxEigenval = MaxEigenval + Y[i]*U[i];
		}
		//若迭代精度和迭代次数达到要求，则迭代成功。
		if(fabs((MaxEigenval - temp) / MaxEigenval) < preci && k > 20)
		{
			return ;
			
		}
	
	}
	cout<<"failed!"<<endl;
	return;

}   

//求按模为最小的特征值。
void AbsMinEigenval(double A[M][N],double &MinEigenval)
{
	int i,k;
	int r = 2;
	int s = 2; //定义上下带宽。
	int time = 10000; //设定迭代次数。
	double preci = 1e-12;  //设定迭代精度。
	double p;
	double temp;
	double Y[N] = {0};
	double U[N] = {0};
	
	//设定向量U的初值。
	for (i=0;i<N;i++)
	{
		U[i] = 1.0;
	} 
	for (k=0;k<time;k++)
	{
		p = 0;
		for (i=0;i<N;i++)
		{
			p = U[i]*U[i] + p;    
		}
		for (i=0;i<N;i++)
		{
			Y[i] = U[i]/sqrt(p);
		}
		Doolittle(A,U,Y);//用Doolittle方法求U。
		temp = MinEigenval;
		MinEigenval = 0;
		for (i=0;i<N;i++)
		{
			MinEigenval = MinEigenval + Y[i]*U[i];
		}
		//若迭代精度和迭代次数达到要求，则迭代成功。
		if(fabs((MinEigenval - temp) / MinEigenval) < preci && k > 20)
		{
			MinEigenval = 1 / MinEigenval;
			return ;	
		}
	}
	cout<<"failed!"<<endl;
	return;
}



//用Doolittle方法解带状方程组。
void Doolittle(double A[M][N],double x[N], double b[N])
{
	int i,j,k,t;
	int s = 2;
	int r = 2;
	double D[M][N];
	double p[N];
	for (i=0;i<M;i++)
	{
		for (j=0;j<N;j++)
		{
			D[i][j] = A[i][j];
		}
	}
	for (j=0;j<N;j++)
	{
		p[j] = b[j];
	}
	//完全按照书上的算法进行计算。
	for (k=0;k<=N-1;k++)
	{
		for (j=k;j<=Min(k+s,N-1);j++)
		{
			for (t=Max(0,Max(k-r,j-s));t<=k-1;t++)
			{
				D[k-j+s][j] = D[k-j+s][j] -  D[k-t+s][t] * D[t-j+s][j];

			}
		}
	
		for (i=k+1;i<=Min(k+r,N-1);i++)
		{
			for (t=Max(0,Max(i-r,k-s));t<=k-1;t++)
			{
				D[i-k+s][k] = D[i-k+s][k] - D[i-t+s][t] * D[t-k+s][k];
				
			}
			D[i-k+s][k] /= D[s][k];
		}
	}
	for (i=1;i<N;i++)
	{
		for (t=Max(0,i-r);t<=i-1;t++)
		{
			p[i] = p[i] - D[i-t+s][t] * p[t];
		}
	}

	x[N-1] = p[N-1] / D[s][N-1];
	for (i=N-2;i>=0;i--)
	{
		x[i] = p[i];
		for (t=i+1;t<=Min(i+s,N-1);t++)
		{
			x[i] = x[i] - D[i-t+s][t] * x[t];
		}
		x[i] = x[i] / D[s][i];
	}
	return;
}

//求带状矩阵的行列式。
double Det(double A[M][N])
{   
	int i,j,k,t;
	int s = 2;
	int r = 2;
	double D[M][N];
	double det;
	for (i=0;i<M;i++)
	{
		for (j=0;j<N;j++)
		{
			D[i][j] = A[i][j];
		}
	}

	//先作Doolittle分解。
	for (k=0;k<=N-1;k++)
	{
		for (j=k;j<=Min(k+s,N-1);j++)
		{
			for (t=Max(0,Max(k-r,j-s));t<=k-1;t++)
			{
				D[k-j+s][j] = D[k-j+s][j] -  D[k-t+s][t] * D[t-j+s][j];

			}
		}
	
		for (i=k+1;i<=Min(k+r,N-1);i++)
		{
			for (t=Max(0,Max(i-r,k-s));t<=k-1;t++)
			{
				D[i-k+s][k] = D[i-k+s][k] - D[i-t+s][t] * D[t-k+s][k];
				
			}
			D[i-k+s][k] /= D[s][k];
		}
	}
	
	
    //分解后矩阵的第S+1行的连乘积就是该矩阵的行列式值。


	 det = 1.0;
	for (j=0;j<=N-1;j++)
	{
		det = det * D[s][j];
	}

	return det;
	
}

//求最大和最小的特征值。
void MaMiEigenval(double A[M][N],double &MaxEigenval,double &MinEigenval)
{
	int i;
	int s = 2;
	int r = 2;
	double k = 1.0;
	AbsMaxEigenval(A,k);//先求绝对值最大的特征根。
	//下面分情况讨论。
	//如果求得的根为正，把它赋给最大的特征值，以k为偏移量，
	//重新得到的绝对值最大的特征值与k之和即为矩阵的最小特征值
	if (k>0)
	{
		MaxEigenval = k;
		for (i=0;i<N;i++)
		{
			A[s][i] = A[s][i] - MaxEigenval;
		}
		AbsMaxEigenval(A,MinEigenval);
		MinEigenval = MinEigenval + MaxEigenval;
		//再恢复矩阵。
		for (i=0;i<N;i++)
		{
			A[s][i] = A[s][i] + MaxEigenval;
		}
	
	}
	//否则以k为偏移量，把它赋给最小的特征值，以k为偏移量，
	//重新得到的绝对值最大的特征值与k之和即为矩阵的最大特征值
	else
	{
		MinEigenval = k;
		for (i=0;i<N;i++)
		{
			A[s][i] = A[s][i] - MinEigenval;
		}
		AbsMaxEigenval(A,MaxEigenval);
		MaxEigenval = MaxEigenval +MinEigenval;
		//再恢复矩阵。
		for (i=0;i<N;i++)
		{
			A[s][i] = A[s][i] + MinEigenval;
		}
		

	}
	
}