北航 数值分析第二次大作业(带双步位移的QR方法)

一、算法设计方案：
按题目要求，本程序运用带双步位移的QR方法求解给定矩阵的特征值，并对每一实特征值，求解其相应的特征向量。

总体思路：
1）初始化矩阵
首先需要将需要求解的矩阵输入程序。

为了防止矩阵在后面的计算中被破坏保存A[][]。

2）对给定的矩阵进行拟上三角化
为了尽量减少计算量，提高程序的运行效率，在对矩阵进行QR分解之前，先进行拟上三角化。

由于矩阵的QR 分解不改变矩阵的结构，所以具有拟上三角形状的矩阵的QR分解可以减少大量的计算量。

这里用函数
void QuasiTriangularization()来实现，函数形参为double型N维方阵double a[][N]。

3）对拟上三角化后的矩阵进行QR分解
对拟上三角化的矩阵进行QR分解会大大减小计算量。

用子程序void QR_decomposition()来实现，将Q、R设为形参，然后将计算出来的结果传入Q和R，然后求出RQ乘积。

4）对拟上三角化后的矩阵进行带双步位移的QR分解
为了加速收敛，对QR分解引入双步位移，适当选取位移量，可以避免进行复数运算。

为了进一步减少计算量，在每次进行QR分解之前，先判断是否可以直接得到矩阵的一个特征值或者通过简单的运算得到矩阵的一对特征值。

若可以，则得到特征值，同时对矩阵进行降阶处理；若不可以，则进行QR分解。

这里用函数intTwoStepDisplacement_QR()来实现。

这是用来存储计算得到的特征值的二维数组。

考虑到特征值可能为复数，因此将所有特征值均当成复数处理。

此函数中，QR分解部分用子函数void QR_decompositionMk()实现。

这里形参有三个，分别用来传递引入双步位移后的Mk阵，A矩阵，以及当前目标矩阵的维数m。

5）计算特征向量
得到特征值后，计算实特征值相应的特征向量。

这里判断所得特征值的虚数部分是否为零。

求实特征值的特征向量采用求解相应的方程组((A-λI)x=0)的方法。

因此先初始化矩阵Array，计算(A-λI)，再求解方程组。

方程组的求解采用列主元的高斯消去法，由函数intGauss_column(double a[][N],double b[],double X[])实现。

由于此给定矩阵的特殊性，其没有重根，所有对应于每一特征值只有一个特征向量，因此可以用高斯消去法求解此奇异的线性方程组。

首先用高斯消去将矩阵(A-λI)化为上三角阵，其最后一行必定全为零。

因此在反代的过程中，只要令x[]的最后一个元素为“1”，即可得到方程组的一个基础解系，从而也就是一个特征向量。

6）输出有关结果
根据题目要求，需要输出拟上三角化后的矩阵、QR分解结束后的矩阵、矩阵全部特征值及对应实特征值的特征向量。

由于输出结果要求都要保留12位有效数字，所以将结果输出到文件result.txt中更有利于数据的打印。

程序中构造了两个输出函数专门来解决不同输出结果的问题，void print_lambda(complex lambda[]);void print_matrix(double mat[][N])。

二、程序源代码：
#include "stdafx.h"
#include "stdlib.h"
#include "math.h"
#define L 100 //定义最大迭代次数
#define N 10 //定义矩阵大小
#define err 1e-12 //定义误差限
//定义一个结构体来表示复数
typedefstruct complex{
double rpart;
double ipart;
};
FILE * pReadFile;
//主函数
int _tmain(intargc, _TCHAR* argv[])
{
inti,j,t;
double y[N],X[N],a[N][N],A[N][N],B[N][N],Q[N][N],R[N][N],RQ[N][N];
struct complex lambda[N];
//声明要调用的函数
void QuasiTriangularization(double a[][N]);
void QR_decomposition(double A[][N],double Q[][N],double R[][N]);
void QR_decompositionMk(double mk[][N],int m, double a[][N]);
void print_lambda(complex lambda[]);
void print_matrix(double mat[][N]);
void multi_matrix(double a[][N],double b[][N],double c[][N]);
intTwoStepDisplacement_QR(double a[][N],complex lambda[]);
intGauss_column(double a[][N],double b[],double X[]);
//矩阵初始化
for (i = 0; i < N; i++){
for (j = 0; j < N; j++){
if (i == j){
a[i][j] = 1.5 * cos((i+1) + 1.2 * (j+1));
A[i][j] = a[i][j];
}
else{
a[i][j] = sin(0.5 * (i+1) + 0.2 * (j+1));
A[i][j] = a[i][j];
}
}
}
for (i = 0; i < N; i++){
y[i] = 0;
}
//对矩阵a[][]拟上三角化
QuasiTriangularization(a);
//打开文件result.txt
pReadFile = fopen( "result.txt", "w" );
//打印结果到文件result.txt中
fprintf(pReadFile,"拟上三角化之后的矩阵A[%d][%d]:\n",N,N);
//printf("拟上三角化之后的矩阵A[%d][%d]:\n",N,N);
print_matrix(a);
//对拟上三角化后的矩阵a[][]，进行QR分解
QR_decomposition(a,Q,R);
fprintf(pReadFile,"Q[%d][%d]:\n",N,N);
//printf("Q[%d][%d]:\n",N,N);
print_matrix(Q);
fprintf(pReadFile,"R[%d][%d]:\n",N,N);
//printf("R[%d][%d]:\n",N,N);
print_matrix(R);
multi_matrix(R,Q,RQ);
fprintf(pReadFile,"RQ[%d][%d]:\n",N,N);
//printf("RQ[%d][%d]:\n",N,N);
print_matrix(RQ);
//双步位移QR分解求全部特征值
TwoStepDisplacement_QR(a,lambda);
//利用校正的列主元素高斯消元法求每个实特征值对应的特征向量
for (t = 0; t < N; t++){
if (lambda[t].ipart == 0){
for (i = 0; i < N; i++){
for (j = 0; j < N; j++){
if (i == j)
B[i][j] = A[i][j] - lambda[1].rpart;
else
B[i][j] = A[i][j];
}
}
Gauss_column(B,y,X);
fprintf(pReadFile,"\n矩阵与特征值λ=%.12e对应的特征向量为X[%d]：\n",lambda[t].rpart,N);
//printf("\n矩阵与特征值λ=%.12e对应的特征向量为X[%d]：\n",lambda[t].rpart,N);
for (i = 0; i < N; i++){
fprintf(pReadFile,"%.12e ",i,X[i]);
//printf("X[%d]: %.12e ",i,X[i]);
}
}
}
fclose( pReadFile );
scanf("%d",&i);
return 0;
}//主函数
/************************************* 矩阵的拟上三角化
输入n阶方阵a[][]，将a[][]拟上三角化
无返回值
**************************************/ void QuasiTriangularization(double a[][N]){
intr,i,j;
double tr,hr,cr,dr,sum;
double ur[N],pr[N],qr[N],wr[N];
for (r = 0; r < N-2; r++){
//判断a[i][r](i=r+2,r+3,...,n-2)是否全为零
sum = 0;//变量sum使用前清零
for (i = r+2; i < N; i++){
sum = sum || a[i][r];
}
//如果不是全部都是零，则计算
if (sum){
//计算dr
sum = 0;
for (i = r+1; i < N; i++){
sum += a[i][r] * a[i][r];
}
dr = sqrt(sum);
//计算cr
if (a[r+1][r] > 0)
cr = -dr;
else
cr = dr;
//计算hr
hr = cr * cr - cr * a[r+1][r];
//取向量ur[]
for (i = 0; i < N; i++){
if (i < r+1)
ur[i] = 0;
else
if (i == r+1)
ur[i] = a[i][r] - cr;
else
ur[i] = a[i][r];
}
//计算向量qr[]
for (i = 0; i < N; i++){
sum = 0;
for (j = r+1; j < N; j++)
sum += a[i][j] * ur[j];
qr[i] = sum / hr;
}
//计算向量pr[]
for (i = 0; i < N; i++){
sum = 0;
for (j = r+1; j < N; j++)
sum += a[j][i] * ur[j];
pr[i] = sum / hr;
}
//计算tr
sum = 0;
for (i = r+1; i < N; i++)
sum += pr[i] * ur[i];
tr = sum / hr;
//计算wr[]
for (i = 0; i < N; i++){
if (i < r+1)
wr[i] = qr[i];
else
wr[i] = qr[i] - tr * ur[i];
}
//计算新产生的矩阵a[][]
for (i = 0; i < N; i++)
for (j = 0; j < N; j++)
a[i][j] = a[i][j] - wr[i] * ur[j] - ur[i] * pr[j];
}
}
}
/*************************************
矩阵的QR分解
将A[][]分解为Q[][]*R[][]
无返回值
**************************************/
void QR_decomposition(double A[][N],double Q[][N],double R[][N]){ intr,i,j;
double tr,hr,cr,dr,sum;
double ur[N],pr[N],wr[N];
//取矩阵R1[][]为A[][]
for (i = 0; i < N; i++)
for (j = 0; j < N; j++)
R[i][j] = A[i][j];
//取矩阵Q1[][]为单位矩阵
for (i = 0; i < N; i++)
for (j = 0; j < N; j++){
if (i == j)
Q[i][j] = 1;
else
Q[i][j] = 0;
}
//
for (r = 0; r < N-1; r++){
//判断R[i][r](i=r+1,r+2,...,N-1)是否全为零
sum = 0;//变量sum使用前清零
for (i = r+1; i < N; i++){
sum = sum || R[i][r];
}
if (sum){
//计算dr
sum = 0;
for (i = r; i < N; i++){
sum += R[i][r] * R[i][r];
}
dr = sqrt(sum);
//计算cr
if (R[r][r] > 0)
cr = -dr;
else
cr = dr;
//计算hr
hr = cr * cr - cr * R[r][r];
//取向量ur[]
for (i = 0; i < N; i++){
if (i < r)
ur[i] = 0;
else
if (i == r)
ur[i] = R[i][r] - cr;
else
ur[i] = R[i][r];
}
//计算向量wr[]
for (i = 0; i < N; i++){
sum = 0;
for (j = r; j < N; j++)
sum += Q[i][j] * ur[j];
wr[i] = sum;
}
//计算新的矩阵Qr+1[][],存储在Q[][]里面
for (i = 0; i < N; i++)
for (j = 0; j < N; j++)
Q[i][j] = Q[i][j] - wr[i] * ur[j] / hr;
//计算向量pr[]
for (i = 0; i < N; i++){
sum = 0;
for (j = r; j < N; j++)
sum += R[j][i] * ur[j];
pr[i] = sum / hr;
}
//计算新产生的R[][]
for (i = 0; i < N; i++)
for (j = 0; j < N; j++)
R[i][j] = R[i][j] - ur[i] * pr[j];
}
}
}
/*************************************
Mk[][]矩阵的QR分解及求Ak+1[][]的计算
形参m为每次降阶之后的值
无返回值
**************************************/
void QR_decompositionMk(double mk[][N],int m, double a[][N]){ intr,i,j;
double tr,hr,cr,dr,sum;
double ur[N],pr[N],qr[N],wr[N],vr[N],br[N][N],Cr[N][N];
//取矩阵br[][]
for (i = 0; i <= m; i++)
for (j = 0; j <= m; j++)
br[i][j] = mk[i][j];
//取矩阵Cr[][]
for (i = 0; i <= m; i++)
for (j = 0; j <= m; j++)
Cr[i][j] = a[i][j];
//
for (r = 0; r <= m-1; r++){
//判断br[i][r](i=r+1,r+2,...,m)是否全为零
sum = 0;//变量sum使用前清零
for (i = r+1; i <= m; i++){
sum = sum || br[i][r];
}
if (sum){
//计算dr
sum = 0;
for (i = r; i <= m; i++){
sum += br[i][r] * br[i][r];
}
dr = sqrt(sum);
//计算cr
if (br[r][r] > 0)
cr = -dr;
else
cr = dr;
//计算hr
hr = cr * cr - cr * br[r][r];
//取向量ur[]
for (i = 0; i <= m; i++){
if (i < r)
ur[i] = 0;
else
if (i == r)
ur[i] = br[i][r] - cr;
else
ur[i] = br[i][r];
}
//计算向量vr[]
for (i = 0; i <= m; i++){
sum = 0;
for (j = r; j <= m; j++)
sum += br[j][i] * ur[j];
vr[i] = sum / hr;
}
//计算新的矩阵Br+1[][],存储在Br[][]里面
for (i = 0; i <= m; i++)
for (j = 0; j <= m; j++)
br[i][j] = br[i][j] - ur[i] * vr[j];
//计算向量pr[]
for (i = 0; i <= m; i++){
sum = 0;
for (j = r; j <= m; j++)
sum += Cr[j][i] * ur[j];
pr[i] = sum / hr;
}
//计算向量qr[]
for (i = 0; i <= m; i++){
sum = 0;
for (j = r; j <= m; j++)
sum += Cr[i][j] * ur[j];
qr[i] = sum / hr;
}
//计算tr
sum = 0;
for (i = r; i <= m; i++)
sum += pr[i] * ur[i];
tr = sum / hr;
//计算向量wr[]
for (i = 0; i <= m; i++){
if (i < r)
wr[i] = qr[i];
else
wr[i] = qr[i] - tr * ur[i];
}
//计算新产生的矩阵Cr[][]
for (i = 0; i <= m; i++)
for (j = 0; j <= m; j++)
Cr[i][j] = Cr[i][j] - wr[i] * ur[j] - ur[i] * pr[j]; }
}
//将计算出的Cr[][]矩阵中的值赋给矩阵A[][]，得到新的矩阵。

for (i = 0; i <= m; i++)
for (j = 0; j <= m; j++)
a[i][j] = Cr[i][j];
}
/*******************************************
高斯列主元素消元法求矩阵实特征值的特征向量
等到的特征向量放在X[]里面
返回值为整型
********************************************/ intGauss_column(double a[][N],double b[],double X[])
{
inti,j,k,row_no;
double temp,a_ik,m_ik,sum = 0;
//消元过程
for (k = 0; k < N - 1; k++){
//选行号，并记录
row_no = k;
a_ik = abs(a[k][k]);
for(i = k + 1; i < N; i++){
if (abs(a[i][k]) >a_ik){
a_ik = abs(a[i][k]);
row_no = i;
}
}
//交换刚刚选择的最大行与第K行的所有元素
if (row_no != k){//如果不是第K行的元素最大，则交换行for (j = k; j < N; j++){
temp = a[k][j];
a[k][j] = a[row_no][j];
a[row_no][j] = temp;
}
temp = b[k];
b[k] = b[row_no];
b[row_no] = temp;
}
if(a[k][k] == 0){
printf("error");
return(0);
}
for(i = k+1; i < N; i++){
m_ik = a[i][k] / a[k][k];
a[i][k] = 0;
for(j = k+1; j < N; j++){
a[i][j] = a[i][j] - m_ik * a[k][j];
}
b[i] = b[i] - m_ik * b[k];
}
}
//回代过程，由于a的行列式值为零，所以有无穷多个解。

所以不妨设x[N-1] := 1 X[N-1] = 1;//b[N-1] / a[N-1][N-1];
for(k = N-2; k >= 0; k--){
sum = 0;
for(j = k+1; j < N; j++)
sum += X[j] * a[k][j];
X[k] = (b[k] - sum) / a[k][k];
}
return 0;
}
/*******************************************
打印特征向量的值到result文件中
无返回值
********************************************/
void print_lambda(complex lambda[]){
int i;
fprintf(pReadFile,"已经求出所有的%d个特征值λ：\n",N);
//printf("已经求出所有的%d个特征值λ：\n",N);
for (i = 0; i < N; i++){
if (lambda[i].ipart != 0){
if (lambda[i].ipart> 0)
fprintf(pReadFile,"%.12e + %.12e i\n",lambda[i].rpart,lambda[i].ipart);
//printf("%.12e + %.12e i\n",lambda[i].rpart,lambda[i].ipart);
else
fprintf(pReadFile,"%.12e + %.12e i\n",lambda[i].rpart,-lambda[i].ipart);
//printf("%.12e - %.12e i\n",lambda[i].rpart,-lambda[i].ipart);
}
else
fprintf(pReadFile,"%.12e\n",lambda[i].rpart);
//printf("%.12e\n",lambda[i].rpart);
}
fprintf(pReadFile,"\n");
//printf("\n");
}
/*******************************************
打印参量传入矩阵到result文件中
无返回值
********************************************/
void print_matrix(double mat[][N]){
inti,j;
for (i = 0; i < N; i++){
for (j = 0; j < N; j++)
fprintf(pReadFile,"%.12e ",mat[i][j]);
//printf("%.12e ",mat[i][j]);
//printf("\n");
fprintf(pReadFile,"\n");
}
}
/*******************************************
带双步位移的QR分解
将所有求出的特征值传入参量lambda[]中
返回值为整型
********************************************/ intTwoStepDisplacement_QR(double a[][N],complex lambda[]){
inti,j,m,l,flag = 1,count = 0,k = 0;
double mk[N][N],detD,b,s,t,sum,delta;
struct complex s1,s2;
m = N - 1;
//(3)到(9)的循环
while (1)
{
k++;
//(3)到(7)的循环
while (1)
{
//(3)到(4)的循环
while (1)
{
if (flag){
//(3)
if (abs(a[m][m-1]) <= err){
lambda[count].rpart = a[m][m];
lambda[count].ipart = 0;
count++;
m--;
}
else
break;
}
flag = 1;
//(4)
if (m == 0){
lambda[count].rpart = a[m][m];
lambda[count].ipart = 0;
print_lambda(lambda);
return 0;
}
else
if (m == -1){
print_lambda(lambda);
return 0;
}
}
//(5)
detD = a[m-1][m-1] * a[m][m] - a[m-1][m] * a[m][m-1];
b = -(a[m-1][m-1] + a[m][m]);
delta = b * b - 4 * detD;
if (delta >= 0){
s1.rpart = -(b + sqrt(delta)) / 2;
s1.ipart = 0;
s2.rpart = -(b - sqrt(delta)) / 2;
s2.ipart = 0;
}
else{
s1.rpart = -b / 2;
s1.ipart = -sqrt(-delta) / 2;
s2.rpart = -b / 2;
s2.ipart = sqrt(-delta) / 2;
}
//(6)
if (m == 1){
lambda[count] = s1;
count++;
lambda[count] = s2;
print_lambda(lambda);
return 0;
}
//(7)
if (abs(a[m-1][m-2]) <= err){
m = m - 2;
lambda[count] = s1;
count++;
lambda[count] = s2;
count++;
}
else
break;
flag = 0;
}
//(8)
if (k == L){
fprintf(pReadFile,"在规定的迭代次数L=%d 内未能得到全部特征值",L);
//printf("在规定的迭代次数L=%d 内未能得到全部特征值",L);
//scanf("%d",&i);
return 0;
}
//(9)
s = a[m-1][m-1] + a[m][m];
t = a[m-1][m-1] * a[m][m] - a[m][m-1] * a[m-1][m];
//计算矩阵Mk
for (i = 0; i <= m; i ++){
for (j = 0; j <= m; j++){
sum = 0;
for (l = 0; l <= m; l++)
sum += a[i][l] * a[l][j];
if (i == j)
mk[i][j] = sum - s * a[i][j] + t;
else
mk[i][j] = sum - s * a[i][j];
}
}
QR_decompositionMk(mk,m,a);
}
}
/*******************************************
矩阵乘法
a[][]与b[][]的乘积传入c[][]中
无返回值
********************************************/ void multi_matrix(double a[][N],double b[][N],double c[][N]){ inti,j,t;
double sum;
for (i = 0; i < N; i++){
for (j = 0; j < N; j++){
sum = 0;
for (t = 0; t < N; t++)
sum += a[i][t] * b[t][j];
c[i][j] = sum;
}
}
}
三、运算结果：
1.拟上三角化之后的矩阵A[10][10]:
-8.827516758830e-001 -9.933136491826e-002 -1.103349285994e+000 -7.600443585637e-001 1.549101079914e-001 -1.946591862872e+000 -8.782436382927e-002 -9.255889387184e-001 6.032599440534e-001 1.518860956469e-001 -2.347878362416e+000 2.372370104937e+000 1.819290822208e+000 3.237804101546e-001 2.205798440320e-001
2.102692662546e+000 1.816138086098e-001 1.278839089990e+000 -6.380578124405e-001 -4.154075603804e-001
0.000000000000e+000 1.728274599968e+000 -1.171467642785e+000 -1.243839262699e+000 -6.399758341743e-001
-2.002833079037e+000 2.924947206124e-001 -6.412830068395e-001 9.783997621285e-002 2.557763574160e-001
0.000000000000e+000 0.000000000000e+000 -1.291669534130e+000 -1.111603513396e+000 1.171346824096e+000
-1.307356030021e+000 1.803699177750e-001 -4.246385358369e-001 7.988955239304e-002 1.608819928069e-001
0.000000000000e+0000.000000000000e+0000.000000000000e+000 1.560126298527e+000 8.125049397515e-001
4.421756832923e-001 -3.588616128137e-002 4.691742313671e-001 -2.736595050092e-001 -7.359334657750e-002
0.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+000 -7.707773755194e-001
-1.583051425742e+000 -3.042843176799e-001 2.528712446035e-001 -6.709925401449e-001 2.544619929082e-001
0.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+000
-7.463453456938e-001 -2.708365157019e-002 -9.486521893682e-001 1.195871081495e-001 1.929265617952e-002
0.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000
e+000 -7.701801374364e-001 -4.697623990618e-001 4.988259468008e-001 1.137691603776e-001
0.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000
e+0000.000000000000e+0007.013167092107e-001 1.582180688475e-001 3.862594614233e-001
0.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000
e+0000.000000000000e+0000.000000000000e+0004.843807602783e-001 3.992777995177e-001
2.QR分解之后拟上三角化之后的矩阵A[10][10]的结果Q[10][10]:
-3.519262579534e-001 4.427591982245e-001 -6.955982513607e-001 6.486200753651e-002 3.709718861896e-001
1.855847143605e-001 -1.628942319628e-002 -1.181053169648e-001 -5.255375383724e-002 -5.486582943568e-002
-9.360277287361e-001 -1.664679186545e-001 2.615299548560e-001 -2.438671728934e-002 -1.394774360893e-001 -6.977585391241e-002 6.124472142963e-003 4.440505443139e-002 1.975907909728e-002 2.062836970533e-002
0.000000000000e+000 -8.810520554692e-001 -3.989762796959e-001 3.720308728479e-002 2.127794064090e-001
1.064463557221e-001 -9.343171079758e-003 -6.774200464527e-002 -3.014340698675e-002 -3.146955080444e-002
0.000000000000e+0000.000000000000e+000 -5.371806806439e-001 -1.234945854205e-001 -7.063151608719e-001
-3.533456368496e-001 3.101438948264e-002 2.248676491598e-001 1.000601783527e-001 1.044622748702e-001
0.000000000000e+0000.000000000000e+0000.000000000000e+000 9.892235468621e-001 -1.239411731211e-001
-6.200358589825e-002 5.442272839461e-003 3.945881637235e-002 1.755813350011e-002 1.833059462907e-002
0.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+000 5.323610690264e-001
-6.733900344896e-001 5.910581205868e-002 4.285425323867e-001 1.906901343193e-001 1.990794495295e-001
0.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+000
-6.0597********e-001 -9.165783032819e-002 -6.645586508974e-001 -2.957110877580e-001 -3.0872********e-001
0.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000
e+000 9.933396625117e-001 -9.690440311939e-002 -4.311990584470e-002 -4.501694411183e-002
0.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000
e+0000.000000000000e+000 5.410088006061e-001 -5.817540838226e-001 -6.0734********e-001
0.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000
e+0000.000000000000e+0000.000000000000e+000 -7.221591336735e-001 6.917269588876e-001
3.QR分解之后拟上三角化之后的矩阵A[10][10]的结果R[10][10]:
2.508342744917e+000 -2.185646885493e+000 -1.314609070786e+000 -
3.558787493835e-002 -2.609857850388e-001 -1.283121847090e+000 -1.390878610606e-001 -8.712897972161e-001 3.849367902971e-001 3.353802899665e-001
0.000000000000e+000 -1.961603277854e+000 2.407523727633e-001 7.054714572823e-001 5.957204318279e-001 5.526978774676e-001 -3.268209924413e-001 -5.769498668364e-002 2.871129330189e-001 -8.895128754189e-002
0.000000000000e+000 0.000000000000e+000 2.404534601993e+000 1.706758096328e+000 -4.239566704091e-001 3.405332305815e+000 -1.050017655852e-001 1.462257102734e+000 -6.684487469283e-001 -4.027*********e-001
0.000000000000e+000 0.000000000000e+0000.000000000000e+000 1.577122080722e+000 6.399535133956e-001 3.468127872427e-001 -5.701786649768e-002 4.014788054433e-001 -2.222476176311e-001 -6.317059236442e-002
0.000000000000e+000 0.000000000000e+0000.000000000000e+0000.000000000000e+000 -1.447846997770e+000 -1.415724007744e+000 -2.806139044665e-001 -2.817910521892e-001 -4.611434881851e-002 1.996629079956e-001
0.000000000000e+000 0.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+000
1.231641451542e+000 1.619701003419e-001 1.962638275504e-001 5.350035621760e-001 -1.509273424767e-001
0.000000000000e+000
0.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+000-7.753441914209e-00
1 -3.464514508821e-001 4.312226803504e-001 1.234643696237e-001
0.000000000000e+000
0.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+00
01.296312940612e+000 -4.288053318338e-001 2.737334158165e-001
0.000000000000e+000
0.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+00 00.000000000000e+000 -6.707396440648e-001 -4.842320121884e-001
0.000000000000e+000
0.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+00 00.000000000000e+0000.000000000000e+000 7.168323926323e-002
4.QR分解之后拟上三角化之后的矩阵A[10][10]的结果Q、R乘积矩阵为:
1.163074414164e+000
2.632670934508e+000 -1.772796003272e+000 -8.668899138521e-002
3.300503471047e-001 1.455162371214e+000 -9.730650448593e-001 -
4.873031174655e-001 -7.756411630489e-001 -3.249201979113e-001
1.836115060851e+000 1.144286420080e-001 -9.880381403133e-001 5.589725694767e-001 4.694190067101e-002 -
2.978478237007e-001 1.617130577649e-002 6.936977702522e-001 1.367670571405e-001 1.419099231519e-002
0.000000000000e+000 -2.118520153533e+000 -1.876189745783e+000 -5.407071940597e-001 1.171538359721e+000 -2.550323020223e+000 1.691577936540e+000 1.229951613262e+000 1.387947777212e+000 8.667502917242e-001
0.000000000000e+0000.000000000000e+000 -8.471995127808e-001 4.382910468318e-001 -1.008632199185e+000 -7.959374261495e-001 4.769258865577e-001 4.072683083890e-001 4.096390493527e-001 3.363378940862e-001
0.000000000000e+0000.000000000000e+0000.000000000000e+000 -1.432244342447e+000 -5.742284908055e-001
1.213151477723e+000 -3.457508625575e-001 -4.749853573124e-001 -3.176158274191e-001 -4.294507015032e-002
0.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+000 6.556779598004e-001 -9.275250974463e-001 2.529079844053e-001 6.905949216976e-001 -2.374430675823e-002 -2.429781119781e-001
0.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+000
4.698400884876e-001 -2.730776009527e-001 7.821296259798e-001 -9.580964936399e-002 7.846239841323e-002
0.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000 e+000 1.287679058937e+000 -3.576058900348e-001 -4.116725408806e-003 3.914268216423e-001
0.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000 e+0000.000000000000e+000-3.628760503545e-001 7.398980975354e-0017.241608309576e-002
0.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000
e+0000.000000000000e+0000.000000000000e+000 -5.176670596524e-002 4.958522909877e-002
5.矩阵的所有特征值λ为:
λ1 =5.650488993501e-002
λ2 =6.360627875745e-001
λ3 =9.355889078188e-001
λ4 =-9.805309562902e-001 + 1.139489127430e-001*i
λ5 =-9.805309562902e-001 + 1.139489127430e-001*i
λ6 =-1.484039822259e+000
λ7 =1.577548557113e+000
λ8 =-2.323496210212e+000 + 8.930405177200e-001*i
λ9 =-2.323496210212e+000 + 8.930405177200e-001*i
λ10=3.383039617436e+000
6.矩阵的所有实特征值对应的特征向量为:
矩阵与特征值λ1=5.650488993501e-002对应的特征向量为X1[10]：
-5.105003830620e+000 -4.886319842355e+000 9.505161576142e+000 -6.788331788207e-001 -9.604334110489e+000 -3.0457********e+000 1.574873885600e+001 -7.395037126270e+000 -7.109654943654e+000 1.000000000000e+000
矩阵与特征值λ2=6.360627875745e-001对应的特征向量为X2[10]：
4.745031936534e+000 3.157868541747e+000 1.729946912419e+001 -1.980049331453e+000
-3.187521973521e+001 7.794009603193e+000 -1.004255685844e+001 1.670757770479e+001
1.310524273052e+001 1.000000000000e+000
矩阵与特征值λ3=9.355889078188e-001对应的特征向量为X3[10]：
2.792418944526e+000 1.598236841510e+000 -5.207507040906e-001 -1.667886451700e+000
-1.225708535858e+001 7.241214790770e+000 -5.398214501428e+000 2.841008912972e+001 -1.216518754415e+0011.000000000000e+000
矩阵与特征值λ6=-1.484039822259e+000对应的特征向量为X6[10]：
-1.730784582110e+001 2.408192534965e+001 4.133852365114e-001 -8.572044074529e+000
9.287334657923e-002 -7.832726042765e-002 -6.374274025708e-001 -3.403204760828e-001
-3.784865409365e-001 1.000000000000e+000
矩阵与特征值λ7=1.577548557113e+000对应的特征向量为X7[10]：
6.217350824575e-001 -1.115111815225e-001 -2.483773580801e+000 -1.306860840419e+000
-3.815605442529e+000 8.117305509387e+000 -1.239170883674e+000 -6.800309586190e-001
2.691900144837e+000 1.000000000000e+000
矩阵与特征值λ10=3.383039617436e+000对应的特征向量为X10[10]：
-4.367398732970e-001 -9.063574014237e-001 -1.963340418326e+000 -1.082857447819e+000 -1.269297783865e+000 -1.0724********e+000 3.625108238697e-001 1.682905297743e+000
2.112744711308e+000 1.000000000000e+000。