北航数值分析大作业第二题精解

目标：使用带双步位移的QR 分解法求矩阵10*10[]ij A a =的全部特征值，并对其中的每一个实特征值求相应的特征向量。已知：sin(0.50.2)()

1.5cos( 1.2)(){i j i j ij i j i j a +≠+== (i,j=1,2, (10)

算法：

以上是程序运作的逻辑，其中具体的函数的算法，大部分都是数值分析课本上的逻辑，在这里特别写出矩阵A 的实特征值对应的一个特征向量的求法：

()[]()()

[]()[]()111111I 00000

i n n n B A I gause i n Q A I u Bu u λλ-⨯-⨯-=-⨯-⎛⎫

⎪-=−−−−→=−−−−−−→= ⎪⎝

⎭

选主元的消元

检查知无重特征值

由于=0i A I λ-，因此在经过选主元的高斯消元以后，i A I λ-即B 的最后一行必然为零，左上方变为n-1阶单位矩阵[]()()11I n n -⨯-，右上方变为n-1阶向量[]()11n Q ⨯-，然后令n u 1=-，则

()1,2,,1j j u Q j n ==⋅⋅⋅-。

这样即求出所有A所有实特征值对应的一个特征向量。

#include

#include

#include

#define N 10

#define E 1.0e-12

#define MAX 10000

//以下是符号函数

double sgn(double a)

{

double z;

if(a>E) z=1;

else z=-1;

return z;

}

//以下是矩阵的拟三角分解

void nishangsanjiaodiv(double A[N][N])

{

int i,j,k;

int m=0;

double d,c,h,t;

double u[N],p[N],q[N],w[N];

for(i=0;i

{

for(j=i+2;j

if(m==(N-2-i)) continue;

for(j=i+1,d=0;j

d=sqrt(d);

c=-1*sgn(A[i+1][i])*d;

h=c*c-c*A[i+1][i];

for(j=i+2;j

for(j=0;j

u[i+1]=A[i+1][i]-c;

for(j=0;j

{for(k=i+1,p[j]=0;k

for(j=0;j

{for(k=i+1,q[j]=0;k

for(j=0,t=0;j

for(j=0;j

for(j=0;j

{

for(k=0;k

}

}

}

//以下是矩阵的QR分解

void qrdiv(double A[N][N],double Q[N][N],double R[N][N])

{

int i,j,k;

//int m=0;

double d,c,h;

double u[N],w[N],p[N];

for(i=0;i

{for(j=0;j

for(i=0;i

{for(j=0;j

for(i=0;i

{

//for(j=i+1;j

//if(m==(N-1-i)) continue;

for(j=i,d=0;j

d=sqrt(d);

c=-1*sgn(R[i][i])*d;

h=c*c-c*R[i][i];

for(j=i+1;j

for(j=0;j

u[i]=R[i][i]-c;

for(j=0;j

for(j=0;j

for(j=0;j

{for(k=i,p[j]=0;k

for(j=0;j

{

for(k=0;k

}

}

}//矩阵的QR分解