数值分析北航大作业 QR分解求特征值

#include
#include

#define N 10
#define L 50
double a[N][N];

double Q[N][N]R[N][N]RQ[N][N];
double lamda[N]Re[N]Im[N];
double M[N][N]QK[N][N];
double biaozhi;
int mcnt;

void initial（void）;				//初始化矩阵A
void nishangsanjiaohua（void）;	//拟上三角化矩阵A
void QR_resolve（void）;			//QR分解拟上三角阵A（n-1）
void eig（void）;					//求矩阵的特征向量
void QR_resolveM（void）;			//对矩阵Mk进行QR分解
void solve（void）;				//解一元二次方程
void jiexiangliang（void）;		//求对应于实数特征值的特征向量



void main（void）
{
	double b[N][N];
	int ijl;

	initial（）;			//初始化矩阵A
	
//显示矩阵A
	printf（“初始矩阵A显示如下\n“）;
	for（i=0;i<=N-1;i++）
	{
		for（j=0;j<=N-1;j++）
		{
			printf（“a[%d][%d] = %.11e		“i+1j+1a[i][j]）;
		}
		printf（“\n“）;
	}

	nishangsanjiaohua（）;		//对矩阵A进行拟上三角化得到A（n-1）

//显示拟上三角化后的矩阵A（n-1）
	printf（“拟上三角化后矩阵A（n-1）显示如下\n“）;
	for（i=0;i<=N-1;i++）
	{
		for（j=0;j<=N-1;j++）
		{
			printf（“a[%d][%d] = %.11e		“i+1j+1a[i][j]）;
		}
		printf（“\n“）;
	}

	for（i=0;i<=N-1;i++）
		for（j=0;j<=N-1;j++）
			b[i][j]=a[i][j];

	QR_resolve（）;		//对A（n-1）进行QR分解，得矩阵Q，R，RQ


//显示矩阵Q
	printf（“对A（n-1）进行QR分解后矩阵Q显示如下\n“）;	
	for（i=0;i<=N-1;i++）
	{
		for（j=0;j<=N-1;j++）
			printf（“Q[%d][%d] = %.11e		“i+1j+1Q[i][j]）;
		printf（“\n“）;
	}

//显示矩阵R
	printf（“对A（n-1）进行QR分解后矩阵R显示如下\n“）;	
	for（i=0;i<=N-1;i++）
	{
		for（j=0;j<=N-1;j++）
		{
			R[i][j]=a[i][j];
			printf（“R[%d][%d] = %.11e		“i+1j+1R[i][j]）;
		}
		printf（“\n“）;
	}

//显示矩阵RQ
	printf（“对A（n-1）进行QR分解后矩阵RQ显示如下\n“）;
	for（i=0;i<=N-1;i++）
		for（j=0;j<=N-1;j++）
		{
			RQ[i][j]=0;
			for（l=0;l<=N-1;l++）
				RQ[i][j]+=R[i][l]*Q[l][j];
		}

	for（i=0;i<=N-1;i++）
		for（j=0;j<=N-1;j++）
		{
			printf（“RQ[%d][%d] = %.11e		“i+1j+1RQ[i][j]）;
		}

	for（i=0;i<=N-1;i++）
		for（j=0;j<=N-1;j++）
			a[i][j]=b[i][j];

//求出所有的特征值
	eig（）;
//显示特征值
	printf（“所有特征值显示如下:\n“）;
	for（i=0;i<=N-1;i++）
		printf（“λ%d（%.11e  %.11e）\n“i+1Re[i]Im[i]）;

	initial（）;		//初始化矩阵A，求A的特征向量

//求对应于实特征值的特征向量
	jiexiangliang（）;

}



//初始化矩阵A
void initial（void）
{
	int ij;
	for（i=0;i<=N-1;i++）
		for（j=0;j<=N-1;j++）
		{
			if（i==j）
				a[i][j]=1.5*cos（（i+1）+1.2*（j+1））;
			else
				a[i][j]=sin（0.5*（i+1）+0.2*（j+1））;
		}
}

//对矩阵A进行拟上三角化得到A（n-1）
void nishangsanjiaohua（void）
{
		
	int ijr;
	double sumcdht;
	double u[N]v[N]p[N]q[N]w[N];

	for（r=0;r<=N-3;r++）
	{
		for（i=r+2;i<=N-1;i++）
			if（fabs（a[i][r]）>=1e-12）break;
		if（i==N && fabs（a[N-1][r]）<=1e-12 ）
			continue;

		sum=0;
		for（i=r+1;i<=N-1;i++）
			sum+=a[i][r]*a[i][r];
		d=sqrt（sum）;
		if（a[r+1][r]>1e-12）
			c=-d;
		else
			c=d;
		h=c*c-c*a[r+1][r];

		for（i=0;i<=r;i++）
			u[i]=0;
		u[r+1]=a[r+1][r]-c;
		for（i=r+2;i<=N-1;i++）
			u[i]=a[i][r];

		for（i=0;i<=N-1;i++）
			v[i]=u[i]/h;
		for（i=0;i<=N-1;i++）
		{
			p[i]=0;
			q[i]=0;
			for（j=0;j<=N-1;j++）
			{
				p[i]+=a[j][i]*v[j];
				q[i]+=a[i][j]*v[j];
			}
		}

		t=0;
		for（i=0;i<=N-1;i++）
			t+=p[i]*v[i];
		for（i=0;i<=N-1;i++