北航 数值分析报告第二次大作业带双步位移地QR方法Word格式文档下载.docx
- 文档编号:20707515
- 上传时间:2023-01-25
- 格式:DOCX
- 页数:24
- 大小:22.21KB
北航 数值分析报告第二次大作业带双步位移地QR方法Word格式文档下载.docx
《北航 数值分析报告第二次大作业带双步位移地QR方法Word格式文档下载.docx》由会员分享,可在线阅读,更多相关《北航 数值分析报告第二次大作业带双步位移地QR方法Word格式文档下载.docx(24页珍藏版)》请在冰豆网上搜索。
由于此给定矩阵的特殊性,其没有重根,所有对应于每一特征值只有一个特征向量,因此可以用高斯消去法求解此奇异的线性方程组。
首先用高斯消去将矩阵(A-λI)化为上三角阵,其最后一行必定全为零。
因此在反代的过程中,只要令x[]的最后一个元素为“1”,即可得到方程组的一个基础解系,从而也就是一个特征向量。
6)输出有关结果
根据题目要求,需要输出拟上三角化后的矩阵、QR分解结束后的矩阵、矩阵全部特征值及对应实特征值的特征向量。
由于输出结果要求都要保留12位有效数字,所以将结果输出到文件result.txt中更有利于数据的打印。
程序中构造了两个输出函数专门来解决不同输出结果的问题,voidprint_lambda(complexlambda[]);
voidprint_matrix(doublemat[][N])。
二、程序源代码:
#include"
stdafx.h"
stdlib.h"
math.h"
#defineL100//定义最大迭代次数
#defineN10//定义矩阵大小
#defineerr1e-12//定义误差限
//定义一个结构体来表示复数
typedefstructcomplex{
doublerpart;
doubleipart;
};
FILE*pReadFile;
//主函数
int_tmain(intargc,_TCHAR*argv[])
{
inti,j,t;
doubley[N],X[N],a[N][N],A[N][N],B[N][N],Q[N][N],R[N][N],RQ[N][N];
structcomplexlambda[N];
//声明要调用的函数
voidQuasiTriangularization(doublea[][N]);
voidQR_decomposition(doubleA[][N],doubleQ[][N],doubleR[][N]);
voidQR_decompositionMk(doublemk[][N],intm,doublea[][N]);
voidprint_lambda(complexlambda[]);
voidprint_matrix(doublemat[][N]);
voidmulti_matrix(doublea[][N],doubleb[][N],doublec[][N]);
intTwoStepDisplacement_QR(doublea[][N],complexlambda[]);
intGauss_column(doublea[][N],doubleb[],doubleX[]);
//矩阵初始化
for(i=0;
i<
N;
i++){
for(j=0;
j<
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));
}
}
y[i]=0;
//对矩阵a[][]拟上三角化
QuasiTriangularization(a);
//打开文件result.txt
pReadFile=fopen("
result.txt"
"
w"
);
//打印结果到文件result.txt中
fprintf(pReadFile,"
拟上三角化之后的矩阵A[%d][%d]:
\n"
N,N);
//printf("
print_matrix(a);
//对拟上三角化后的矩阵a[][],进行QR分解
QR_decomposition(a,Q,R);
Q[%d][%d]:
print_matrix(Q);
R[%d][%d]:
print_matrix(R);
multi_matrix(R,Q,RQ);
RQ[%d][%d]:
print_matrix(RQ);
//双步位移QR分解求全部特征值
TwoStepDisplacement_QR(a,lambda);
//利用校正的列主元素高斯消元法求每个实特征值对应的特征向量
for(t=0;
t<
t++){
if(lambda[t].ipart==0){
for(i=0;
for(j=0;
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]:
lambda[t].rpart,N);
//printf("
fprintf(pReadFile,"
%.12e"
i,X[i]);
//printf("
X[%d]:
%.12e"
fclose(pReadFile);
scanf("
%d"
&
i);
return0;
}//主函数
/*************************************
矩阵的拟上三角化
输入n阶方阵a[][],将a[][]拟上三角化
无返回值
**************************************/
voidQuasiTriangularization(doublea[][N]){
intr,i,j;
doubletr,hr,cr,dr,sum;
doubleur[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;
sum=sum||a[i][r];
//如果不是全部都是零,则计算
if(sum){
//计算dr
sum=0;
for(i=r+1;
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[]
if(i<
r+1)
ur[i]=0;
else
if(i==r+1)
ur[i]=a[i][r]-cr;
ur[i]=a[i][r];
//计算向量qr[]
sum=0;
for(j=r+1;
j++)
sum+=a[i][j]*ur[j];
qr[i]=sum/hr;
//计算向量pr[]
sum+=a[j][i]*ur[j];
pr[i]=sum/hr;
//计算tr
i++)
sum+=pr[i]*ur[i];
tr=sum/hr;
//计算wr[]
wr[i]=qr[i];
else
wr[i]=qr[i]-tr*ur[i];
//计算新产生的矩阵a[][]
i++)
j++)
a[i][j]=a[i][j]-wr[i]*ur[j]-ur[i]*pr[j];
}
矩阵的QR分解
将A[][]分解为Q[][]*R[][]
voidQR_decomposition(doubleA[][N],doubleQ[][N],doubleR[][N]){
doubleur[N],pr[N],wr[N];
//取矩阵R1[][]为A[][]
R[i][j]=A[i][j];
//取矩阵Q1[][]为单位矩阵
if(i==j)
Q[i][j]=1;
else
Q[i][j]=0;
//
N-1;
//判断R[i][r](i=r+1,r+2,...,N-1)是否全为零
for(i=r+1;
sum=sum||R[i][r];
for(i=r;
sum+=R[i][r]*R[i][r];
if(R[r][r]>
hr=cr*cr-cr*R[r][r];
r)
if(i==r)
ur[i]=R[i][r]-cr;
ur[i]=R[i][r];
//计算向量wr[]
for(j=r;
sum+=Q[i][j]*ur[j];
wr[i]=sum;
//计算新的矩阵Qr+1[][],存储在Q[][]里面
Q[i][j]=Q[i][j]-wr[i]*ur[j]/hr;
sum+=R[j][i]*ur[j];
//计算新产生的R[][]
R[i][j]=R[i][j]-ur[i]*pr[j];
}
Mk[][]矩阵的QR分解及求Ak+1[][]的计算
形参m为每次降阶之后的值
voidQR_decompositionMk(doublemk[][N],intm,doublea[][N]){
doubleur[N],pr[N],qr[N],wr[N],vr[N],br[N][N],Cr[N][N];
//取矩阵br[][]
=m;
br[i][j]=mk[i][j];
//取矩阵Cr[][]
Cr[i][j]=a[i][j];
=m-1;
//判断br[i][r](i=r+1,r+2,...,m)是否全为零
sum=sum||br[i][r];
sum+=br[i][r]*br[i][r];
if(br[r][r]>
hr=cr*cr-cr*br[r][r];
ur[i]=br[i][r]-cr;
ur[i]=br[i][r];
//计算向量vr[]
sum+=br[j][i]*ur[j];
vr[i]=sum/hr;
//计算新的矩阵Br+1[][],存储在Br[][]里面
br[i][j]=br[i][j]-ur[i]*vr[j];
sum+=Cr[j][i]*ur[j];
sum+=Cr[i][j]*ur[j];
//计算新产生的矩阵Cr[][]
Cr[i][j]=Cr[i][j]-wr[i]*ur[j]-ur[i]*pr[j];
//将计算出的Cr[][]矩阵中的值赋给矩阵A[][],得到新的矩阵。
a[i][j]=Cr[i][j];
/*******************************************
高斯列主元素消元法求矩阵实特征值的特征向量
等到的特征向量放在X[]里面
返回值为整型
********************************************/
intGauss_column(doublea[][N],doubleb[],doubleX[])
inti,j,k,row_no;
doubletemp,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;
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;
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;
m_ik=a[i][k]/a[k][k];
a[i][k]=0;
for(j=k+1;
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 北航 数值分析报告第二次大作业带双步位移地QR方法 数值 分析 报告 第二次 作业 带双步 位移 QR 方法