数值分析第二次作业.docx
《数值分析第二次作业.docx》由会员分享，可在线阅读，更多相关《数值分析第二次作业.docx（25页珍藏版）》请在冰点文库上搜索。
数值分析第二次作业
数值分析大作业
（二）
院系：
学号：
姓名：
一、算法设计方案：
1、矩阵A的存储和特征值的数据结构设计
A是一个10×10的方阵，且由给定的元素计算公式知A是非对称的矩阵，因此采取二维数组来存储数组A。
考虑到特征值为复数的情况，采用结构体来存储特征值的实部和虚部，当特征值为实数时，虚部赋值为零。
2、求解过程的分析
首先利用HouseHolder矩阵将A拟上三角化，得到A（n-1）。
然后利用带双步位移的QR分解法求解矩阵A的特征值。
利用选列主元素的高斯消去法求解A的实特征值对应的特征向量，将等式两边移项后，相当于求解齐次线性方程组。
在计算前首先要求齐次方程组的系数矩阵的秩，从而确定自由量,赋予适当的初值，从而可以求得特征值对应的一个特征值。
利用QR分解的基本方法可以得到A（n-1）进行QR分解后的Q、R和RQ。
二、源程序
#include
#include
#defineN10
#defineepsilon1E-12
#defineL10000
typedefstructComp
{
doublereal;
doubleimage;
}Comp;
voidgetRealMatrix（intn,doubleA[][n],doubleQ[][n],doubleregulate）;
intisContainThisValue（intm,intmainIndex[],intindex）;
intmaxElementIndex（intm,intn,intcolNum,doublearray[][n],intmainIndex[]）;
voidmainElementElimination（intn,doubleequationArray[][n],doublesolveArray[],double）;
voidqrMethod（int,doublea[][]）;
intsgn（doublerealNum）;
voidoutputMatrix（intm,intn,doubleA[][n]）;
voidcalcMatrix（intm,intn,doubleA[][n],doubleM[][n],doubles,doublet）;
voidproductMatrix（intm,intn,doubleA[][n],doubleM[][n]）;
voidhessenbergMatrix（intn,doublematrix_a[][]）;
voidsolvePowerEquation（doubleb,doublec,ComptwoDArr[2]）;
voidqrIterateFormula（intm,intn,doubleM[][n],doubleA[][n]）;
voidinitialMatrix（intn,doubleA[][]）;
doublecrossProduct（double*p_1,double*p_2,intstartIndex,intendIndex,intincrement）;
voidtranspositionMatirx（intm,intn,doublematrix[][]）;
voiditerateFormula（intm,intn,doublematrix[][],doublew[],doubleu[],doublep[]）;
voidqrWithTwoStepMove（intn,doublematrix[][],doublematrix_Q[][],ComplamdaValue[],doublee,doubleiterCount）;
intmain（）
{
doubleMatrix_Q[N][N],Matrix_A[N][N];
Complamdas[N];
doublevector[N]={0};
inti,j;
//初始化矩阵
initialMatrix（N,Matrix_A）;
//拟上三角化化
hessenbergMatrix（N,Matrix_A）;
//输出拟上三角化后的矩阵
printf（"A拟上三角化的矩阵：
\n"）;
outputMatrix（N,N,Matrix_A）;
//双步位移的QR
qrWithTwoStepMove（N,Matrix_A,Matrix_Q,lamdas,epsilon,L）;
//求解Q、R、RQ
initialMatrix（N,Matrix_A）;
hessenbergMatrix（N,Matrix_A）;
qrMethod（N,Matrix_A）;
for（i=0;i<=N-1;i++）
{
printf（"%.12E,%.12E\n",lamdas[i].real,lamdas[i].image）;
}
//求解实特征值对应的特征向量
initialMatrix（N,Matrix_A）;
for（i=0;i<=N-1;i++）
{
if（lamdas[i].image!
=0）
continue;
else
{
getRealMatrix（N,Matrix_A,Matrix_Q,lamdas[i].real）;
mainElementElimination（N,Matrix_Q,vector,1）;
printf（"%.12E的特征向量为\n",lamdas[i].real）;
for（j=0;j<=N-1;j++）
printf（"%.12E",vector[j]）;
printf（"\n"）;
}
}
return0;
}
//输出矩阵
voidoutputMatrix（intm,intn,doubleA[][n]）
{
inti,j;
for（i=0;i<=m-1;i++）
{
for（j=0;j<=m-1;j++）
{
printf（"%0.12E",A[i][j]）;
}
printf（"\n"）;
}
printf（"\n"）;
}
//初始化矩阵的元素
voidinitialMatrix（intn,doubleA[][n]）
{
inti,j;
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））;
}
}
//对矩阵进行拟上三角化
voidhessenbergMatrix（intn,doublematrix_a[][n]）
{
intr,i;
doubleu[n],p[n],w[n],c,d,h;
for（r=0;r<=n-3;r++）
{
if（crossProduct（&matrix_a[r+2][r],&matrix_a[r+2][r],r+2,n-1,n）==0）
continue;
d=sqrt（crossProduct（&matrix_a[r+1][r],&matrix_a[r+1][r],r+1,n-1,n））;
c=matrix_a[r+1][r]==0?
d:
（-sgn（matrix_a[r+1][r]）*d）;
h=c*c-c*matrix_a[r+1][r];
for（i=0;i<=n-1;i++）
{
u[i]=（i<（r+1））?
0:
（i==（r+1））?
（matrix_a[r+1][r]-c）:
matrix_a[i][r];
}
transpositionMatirx（n,n,matrix_a）;
for（i=0;i<=n-1;i++）
{
p[i]=crossProduct（&matrix_a[0][0]+n*i,u,0,n-1,1）/h;
}
transpositionMatirx（n,n,matrix_a）;
for（i=0;i<=n-1;i++）
{
w[i]=（crossProduct（&matrix_a[0][0]+n*i,u,0,n-1,1）
-crossProduct（p,u,0,n-1,1）*u[i]）/h;
}
iterateFormula（n,n,matrix_a,w,u,p）;
}
}
//双步位移的QR算法实现
voidqrWithTwoStepMove（intn,doublematrix[][n],doublematrix_Q[][n],ComplamdaValue[],doublee,doubleiterCount）
{
intk,m,i=0,index=0;
doubles,t;
Compsolver[2];
k=1;
m=n-1;
while
（1）
{
if（fabs（matrix[m][m-1]）<=e）
{
lamdaValue[index].real=matrix[m][m];
lamdaValue[index].image=0;
index++;
m=m-1;
if（m<2）
{
if（m==0）
{
lamdaValue[index].real=matrix[m][m];
lamdaValue[index].image=0;
index++;
}
elseif（m==1）
{
s=matrix[m-1][m-1]+matrix[m][m];
t=matrix[m-1][m-1]*matrix[m][m]-matrix[m][m-1]*matrix[m-1][m];
solvePowerEquation（s,t,solver）;
for（i=0;i<=1;i++,index++）
lamdaValue[index]=solver[i];
}
return;
}
}
else
{
if（fabs（matrix[m-1][m-2]）<=e）
{
s=matrix[m-1][m-1]+matrix[m][m];
t=matrix[m-1][m-1]*matrix[m][m]-matrix[m][m-1]*matrix[m-1][m];
solvePowerEquation（s,t,solver）;
for（i=0;i<=1;i++,index++）
lamdaValue[index]=solver[i];
m=m-2;
if（m<2）
{
if（m==0）
{
lamdaValue[index].real=matrix[m][m];
lamdaValue[index].image=0;
index++;
}
elseif（m==1）
{
s=matrix[m-1][m-1]+matrix[m][m];
t=matrix[m-1][m-1]*matrix[m][m]-matrix[m][m-1]*matrix[m-1][m];
solvePowerEquation（s,t,solver）;
for（i=0;i<=1;i++,index++）
lamdaValue[index]=solver[i];
}
return;
}
}
else
{
if（k==iterCount）
return;
else
{
s=matrix[m-1][m-1]+matrix[m][m];
t=matrix[m-1][m-1]*matrix[m][m]-matrix[m][m-1]*matrix[m-1][m];
calcMatrix（m,n,matrix,matrix_Q,s,t）;
qrIterateFormula（m+1,n,matrix_Q,matrix）;
k=k+1;
}
}
}
}
}
//矩阵拟上三角化的迭代式子
voiditerateFormula（intm,intn,doublematrix[][n],doublew[],doubleu[],doublep[]）
{
inti,j;
for（i=0;i<=m-1;i++）
for（j=0;j<=m-1;j++）
{
matrix[i][j]=matrix[i][j]-w[i]*u[j]-u[i]*p[j];
}
}
//计算向量的内积，startIndex是起始下标，endIndex是结束下标，increment是元素增量
doublecrossProduct（double*p_1,double*p_2,intstartIndex,intendIndex,intincrement）
{
doublesum=0;
inti;
if（p_1==p_2）
{
for（i=startIndex;i<=endIndex;i++）
{
sum+=（*p_1）*（*p_1）;
p_1=p_1+increment;
}
}
else
{
for（i=startIndex;i<=endIndex;i++）
{
sum+=（*p_1）*（*p_2）;
p_1=p_1+increment;
p_2=p_2+increment;
}
}
returnsum;
}
//矩阵转置
voidtranspositionMatirx（intm,intn,doublematrix[][n]）
{
inti,j;
doublechangeVar;
for（i=0;i<=m-1;i++）
{
for（j=i;j<=m-1;j++）
{
if（i!
=j）
{
changeVar=matrix[i][j];
matrix[i][j]=matrix[j][i];
matrix[j][i]=changeVar;
}
}
}
}
//符号函数
intsgn（doublerealNum）
{
returnrealNum==0?
0:
（realNum>0?
1:
-1）;
}
//代入的矩阵
voidcalcMatrix（intm,intn,doubleA[][n],doubleM[][n],doubles,doublet）
{
inti,j;
productMatrix（m,n,A,M）;
for（i=0;i<=m;i++）
for（j=0;j<=m;j++）
{
M[i][j]=M[i][j]-s*A[i][j]+（（i==j）?
t:
0）;
}
}
//计算A平方
voidproductMatrix（intm,intn,doubleA[][n],doubleM[][n]）
{
inti,j,r;
doublesum;
for（i=0;i<=m;i++）
for（j=0;j<=m;j++）
{
sum=0;
for（r=0;r<=m;r++）
sum+=A[i][r]*A[r][j];
M[i][j]=sum;
}
}
//计算lamda
voidsolvePowerEquation（doubleb,doublec,ComptwoDArr[2]）
{
doubledelta=b*b-4*c;
if（delta<0）
{
twoDArr[0].real=b/2;
twoDArr[0].image=sqrt（fabs（delta））/2;
twoDArr[1].real=b/2;
twoDArr[1].image=-sqrt（fabs（delta））/2;
}
else
{
twoDArr[0].real=（b+sqrt（delta））/2;
twoDArr[0].image=0;
twoDArr[1].real=（b-sqrt（delta））/2;
twoDArr[1].image=0;
}
}
//双步移位QR中的迭代公式
voidqrIterateFormula（intm,intn,doubleM[][n],doubleA[][n]）
{
intr,i;
doublec,d,h,u[m],v[m],p[m],w[m];
for（r=0;r<=m-2;r++）
{
if（crossProduct（&（M[r+1][r]）,&（M[r+1][r]）,r+1,m-1,n）==0）
{
continue;
}
else
{
d=sqrt（crossProduct（&M[r][r],&M[r][r],r,m-1,n））;
c=（M[r][r]==0）?
d:
（-sgn（M[r][r]）*d）;
h=c*c-c*M[r][r];
for（i=0;i<=m-1;i++）
{
u[i]=（i0:
（i==r）?
（M[r][r]-c）:
M[i][r];
}
transpositionMatirx（m,n,M）;
for（i=0;i<=m-1;i++）
{
v[i]=crossProduct（&M[0][0]+n*i,u,0,m-1,1）/h;
}
transpositionMatirx（m,n,M）;
for（i=0;i<=m-1;i++）
p[i]=0;
iterateFormula（m,n,M,u,v,p）;
//第二部分计算
transpositionMatirx（m,n,A）;
for（i=0;i<=m-1;i++）
{
p[i]=crossProduct（&A[0][0]+n*i,u,0,m-1,1）/h;
}
transpositionMatirx（m,n,A）;
for（i=0;i<=m-1;i++）
{
w[i]=（crossProduct（&A[0][0]+n*i,u,0,m-1,1）
-crossProduct（p,u,0,m-1,1）*u[i]）/h;
}
iterateFormula（m,n,A,w,u,p）;
}
}
}
//基本QR方法
voidqrMethod（intn,doubleA[][n]）
{
intr,i,j;
doublec,d,h,sum;
doubleu[n],w[n],p[n];
doubleQ[n][n];
doubleR[n][n];
//初始化为单位阵
for（i=0;i<=n-1;i++）
for（j=0;j<=n-1;j++）
{
if（i==j）
{
Q[i][j]=1;
R[i][j]=0;
}
else
{
Q[i][j]=0;
R[i][j]=0;
}
}
for（r=0;r<=n-1;r++）
{
//判断第r列元素从r+1行开始是否全部为0
if（crossProduct（&A[r+1][r],&A[r+1][r],r+1,n-1,n）==0）
continue;
else
{
d=sqrt（crossProduct（&A[r][r],&A[r][r],r,n-1,n））;
c=（A[r][r]==0）?
d:
（-sgn（A[r][r]）*d）;
h=c*c-c*A[r][r];
//初始化向量u
for（i=0;i<=n-1;i++）
{
u[i]=（i0:
（i==r）?
（A[r][r]-c）:
A[i][r];
}
for（i=0;i<=n-1;i++）
{
w[i]=crossProduct（&Q[0][0]+n*i,u,0,n-1,1）/h;
}
//辅助作用
for（i=0;i<=n-1;i++）
p[i]=0;
//迭代出Q
iterateFormula（n,n,Q,w,u,p）;
transpositionMatirx（n,n,A）;
for（i=0;i<=n-1;i++）
{
p[i]=crossProduct（&A[0][0]+n*i,u,0,n-1,1）/h;
}
//辅助作用
for（i=0;i<=n-1;i++）
w[i]=0;
transpositionMatirx（n,n,A）;
iterateFormula（n,n,A,w,u,p）;
}
}
printf（"正交阵Q:
\n"）;
outputMatrix（n,n,Q）;
printf（"上三角阵R:
\n"）;
outputMatrix（n,n,A）;
for（i=0;i<=n-1;i++）
{
for（j=0;j<=n-1;j++）
{
sum=0;
for（r=i;r<=n-1;r++）
sum+=A[i][r]*Q[r][j];
R[i][j]=sum;
}
}
printf（"输出RQ：
\n"）;
outputMatrix（N,N,R）;
}
//得到齐次线性方程组的系数矩阵
voidgetRealMatrix（intn,doubleA[][n],doubleQ[][n],doubleregulate）
{
inti,j;
for（i=0;i<=n-1;i++）
for（j=0;j<=n-1;j++）
{
升级会员 