北航研究生数值分析作业第二题.docx

北航研究生数值分析作业第二题.docx

《北航研究生数值分析作业第二题.docx》由会员分享，可在线阅读，更多相关《北航研究生数值分析作业第二题.docx（18页珍藏版）》请在冰豆网上搜索。

北航研究生数值分析作业第二题

北航研究生数值分析作业第二题：

一、算法设计方案

1.按照题目给出的矩阵定义对矩阵A赋初值：

对应的函数为a_init（）；

2.对矩阵A进行householder变换，使其拟上三角化：

对应的函数为householder（）；

3.输出拟上三角化后的A：

对应的函数为aout（int）；

4.对拟上三角化后的矩阵A使用带双步位移的QR分解法逐次迭代（最大迭代次数L=500），逐个求出其特征值，对应的函数为eigen_a（）；中间包含两个子程序：

calc_mk（）和qr_analyze（），分别用来计算矩阵Mk和对Mk进行QR分解并得到Ak+1;

5.输出QR分解过程完毕后的A及求得的特征向量：

对应的函数为aout（）和eigenvalout（）；

6.对于在第三步中求得的每个实特征值，使用带原点平移的反幂法求出其对应的特征向量，对应的函数为eigenvec（）；其中包含一个解方程（A-μI）=yk-1的程序段。

这部分也用迭代完成，仍然将最大迭代次数L设置为500；

7.输出矩阵A的特征向量，结束计算：

对应的函数为eigenvecout（）。

算法编译环境：

vlsualc++6.0

二、源程序如下：

#include

#include

#defineN10//矩阵阶数；

#defineEPSL1.0e-12//迭代的精度水平；

#defineL500//迭代最大次数；

#defineOUTPUTMODE1//输出格式：

0--输出至屏幕，1--输出至文件

doublea[N][N],a2[N][N],eigen[N][N];//声明矩阵A；

doublesa_re[N]={0},sa_im[N]={0};//声明矩阵的特征值数组；

doubleu_init[N]={2,1,2,1,2,1,2,1,2,1};//定义反幂法中使用的初始向量u；

//主程序开始；

intmain（）

{

FILE*p;

voida_init（）;

voidhouseholder（）;

voidequal_zero（doublematrix[N][N],int）;

voideigenvec（）;

inteigen_a（）;

voidaout（int）;

voideigenvalout（int）;

voideigenvecout（int）;

if（OUTPUTMODE）

{

p=fopen（"Result.txt","w+"）;

fprintf（p,"计算结果：

\n"）;

fclose（p）;

}

a_init（）;//对矩阵A进行初始化；

householder（）;//对矩阵A进行拟上三角化；

equal_zero（a,N）;//对矩阵A的元素进行归零处理，消除误差；

aout（OUTPUTMODE）;//输出A；

if（eigen_a（））printf（"迭代超过最大次数，特征值求解结果可能不正确。

\n"）;//求矩阵A的特征值；

equal_zero（a,N）;//对矩阵A的元素进行归零处理，消除误差；

aout（OUTPUTMODE）;//输出A；

eigenvalout（OUTPUTMODE）;//输出矩阵的特征值；

eigenvec（）;//求矩阵A的特征向量；

eigenvecout（OUTPUTMODE）;//输出矩阵A的特征向量；

getchar（）;

}

voida_init（）

{

inti,j;

for（i=1;i<=N;i++）

{

for（j=1;j<=N;j++）

a2[i-1][j-1]=a[i-1][j-1]=sin（0.5*i+0.2*j）;

}

for（i=1;i<=N;i++）

a2[i-1][i-1]=a[i-1][i-1]=cos（i+1.2*i）*1.5;//这里使用if语句反而更慢，所以赋值赋了两次。

}

voidhouseholder（）//对矩阵进行拟上三角化；

{

intr,i,j;

doubletmp_ir,tmp_dr,tmp_pr,tmp_qr,tmp_tr,dr,cr,hr,tr;

doubleur[N]={0},pr[N]={0},qr[N]={0},wr[N]={0};

intsgn2（double）;

for（r=0;r

{

//第一步；

tmp_ir=0;

for（i=r+2;i

if（tmp_ir==N-2-r）continue;

else

{

//第二步；

tmp_dr=0;

for（i=r+1;i

dr=sqrt（tmp_dr）;

cr=-1.0*sgn2（a[r+1][r]）*dr;

hr=cr*cr-cr*a[r+1][r];

//第三步；

for（i=0;i

for（i=r+2;i

ur[r+1]=a[r+1][r]-cr;

//第四步；

for（i=0;i

{

tmp_pr=0;

tmp_qr=0;

for（j=0;j

{

tmp_pr=tmp_pr+a[j][i]*ur[j];

tmp_qr=tmp_qr+a[i][j]*ur[j];

}

pr[i]=tmp_pr/hr;

qr[i]=tmp_qr/hr;

}

tmp_tr=0;

for（i=0;i

tr=tmp_tr/hr;

for（i=0;i

for（i=0;i

{

for（j=0;j

{

a[i][j]=a[i][j]-wr[i]*ur[j]-ur[i]*pr[j];

}

}

}

//第五步：

（继续）

}

}

intsgn2（doublea）//求cr时用到的sgn子程序

{

if（a>=0）return1;

elsereturn-1;

}

voidequal_zero（doublematrix[N][N],intrank）//对矩阵进行归零处理；

{

inti,j;

for（i=0;i

{

for（j=0;j

}

}

inteigen_a（）//计算A的特征值；

{

intsnum=N-1,m=N-1,flag=0,flag_7to4=0,step4_cont=0,k=1;

doublemk[N][N],det_dk,stmpb,s1_re,s1_im,s2_re,s2_im,ms,mt,b24ac;

voidcalc_mk（doublemk[N][N],int,double,double）;

voidqr_analyze（doublemk[N][N],int）;

while

（1）

{

//第三步；

if（（-EPSL

flag_7to4）

{

sa_re[snum]=a[m][m],snum--,m--,flag=1;

}

else{}

//第四步；

flag_7to4=0;

if（flag）

{

flag=0;

if（m==0）

{

sa_re[snum]=a[0][0];

return0;

}

elseif（m==-1）return0;

elsestep4_cont=1;

}

//第五步；

else

{

det_dk=a[m-1][m-1]*a[m][m]-a[m-1][m]*a[m][m-1];

stmpb=-1.0*（a[m-1][m-1]+a[m][m]）;

b24ac=stmpb*stmpb-4.0*det_dk;

if（b24ac<0）

{

b24ac=-b24ac;

s1_re=-0.5*stmpb;

s1_im=0.5*sqrt（b24ac）;

s2_re=-0.5*stmpb;

s2_im=-0.5*sqrt（b24ac）;

}

else

{

s1_re=（-1.0*stmpb+sqrt（b24ac））/2;

s1_im=0;

s2_re=（-1.0*stmpb-sqrt（b24ac））/2;

s2_im=0;

}

}

if（step4_cont）

{

step4_cont=0;

continue;

}

//第六步；

if（m==1）

{

sa_re[1]=s2_re;

sa_im[1]=s2_im;

sa_re[0]=s1_re;

sa_im[0]=s1_im;

return0;

}

//第七步；

else

{

if（-EPSL

{

sa_re[snum]=s2_re,sa_im[snum]=s2_im,snum--,m--;

sa_re[snum]=s1_re,sa_im[snum]=s1_im,snum--,m--;

flag=flag_7to4=1;

}

else

{

//第八步；

if（k==L）return1;

//第九步；

ms=a[m-1][m-1]+a[m][m];

mt=a[m-1][m-1]*a[m][m]-a[m-1][m]*a[m][m-1];

calc_mk（mk,m,ms,mt）;//计算矩阵Mk；

qr_analyze（mk,m）;//对Mk进行QR分解并更新Ak；

//第十步；

k++;

}

}

}

}

voidcalc_mk（doublemk[N][N],intm,doublems,doublemt）//函数eigen_a（）中计算Mk的子程序；

{

inti,j,k;

doubletmp_sum;

for（i=0;i<=m;i++）

{

for（j=0;j<=m;j++）

{

tmp_sum=0;

for（k=0;k<=m;k++）tmp_sum+=a[i][k]*a[k][j];

mk[i][j]=tmp_sum-ms*a[i][j];

}

}

for（i=0;i<=m;i++）mk[i][i]+=mt;

}

voidqr_analyze（doublemk[N][N],intm）//函数eigen_a（）中对Mk进行QR分解并更新Ak的子程序

{

inttmp_ir,r,i,j;

doubleur[N],vr[N],pr[N],qr[N],wr[N],tmp_dr,tmp_vr,tmp_pr,tmp_qr,tmp_tr,dr,cr,hr,tr;

for（r=0;r

{

//第一步；

tmp_ir=0;

for（i=r+1;i<=m;i++）

tmp_ir=tmp_ir+（mk[i][r]==0）;

if（tmp_ir==m-r）

continue;

else

{

//第二步；

tmp_dr=0;

for（i=r;i<=m;i++）tmp_dr=tmp_dr+mk[i][r]*mk[i][r];

dr=sqrt（tmp_dr）;

cr=-1.0*sgn2（mk[r][r]）*dr;

hr=cr*cr-cr*mk[r][r];

//第三步；

for（i=0;i

for（i=r+1;i<=m;i++）ur[i]=mk[i][r];

ur[r]=mk[r][r]-cr;

//第四步；

for（i=0;i<=m;i++）

{

tmp_vr=0;

tmp_pr=0;

tmp_qr=0;

for（j=0;j<=m;j++）

{

tmp_vr=tmp_vr+mk[j][i]*ur[j];

tmp_pr=tmp_pr+a[j][i]*ur[j];

tmp_qr=tmp_qr+a[i][j]*ur[j];

}

vr[i]=tmp_vr/hr;

pr[i]=tmp_pr/hr;

qr[i]=tmp_qr/hr;

}

tmp_tr=0;

for（i=0;i<=m;i++）tmp_tr=tmp_tr+pr[i]*ur[i];

tr=tmp_tr/hr;

for（i=0;i<=m;i++）wr[i]=qr[i]-tr*ur[i];

for（i=0;i<=m;i++）

{

for（j=0;j<=m;j++）

{

mk[i][j]=mk[i][j]-ur[i]*vr[j];

a[i][j]=a[i][j]-wr[i]*ur[j]-ur[i]*pr[j];

}

}

}

//第五步：

（继续）

}

}

voideigenvec（）//计算矩阵A的特征向量；

{

inti,j,k,s=0,ik=0,n=0;

doubleu[N],y[N],mik,beta=0,beta_tmp,yita,tmp_yita,tmp;

doubleabs（double）;

voidaout（）;

for（i=0;i

for（s=0;s

{

if（!

sa_im[s]）

{

n=0;

do

{

for（i=0;i

{

for（j=0;j

}

tmp_yita=0;

for（i=0;i

yita=sqrt（tmp_yita）;

for（i=0;i

for（i=0;i

if（!

n）beta=beta_tmp;

else

{

if（abs（（beta_tmp-beta）/beta）

{

for（i=0;i

break;

}

elsebeta=beta_tmp;

}

//解方程组（A-uI）uk=yk；

for（k=0;k

{

ik=k;

for（i=k;i

for（j=k;j

{

tmp=a[k][j];a[k][j]=a[ik][j];a[ik][j]=tmp;

}

tmp=y[k];y[k]=y[ik];y[ik]=tmp;

for（i=k+1;i

{

mik=a[i][k]/a[k][k];

for（j=k;j

y[i]-=mik*y[k];

}

}

u[N-1]=y[N-1]/a[N-1][N-1];

for（k=N-2;k>=0;k--）

{

tmp=0;

for（j=k+1;j

u[k]=（y[k]-tmp）/a[k][k];

}

beta_tmp=0;

for（i=0;i

}

while（n++<=L）;

}

else{}

}

}

doubleabs（doublex）//求绝对值；

{

if（x<0）return-x;

returnx;

}

voidaout（intopmode）//打印A；opmode=1:

到文件，0:

到屏幕；

{

inti,j;

FILE*p;

if（opmode）

{

p=fopen（"Result.txt","at+"）;

fprintf（p,"Matrix（A）=\n{\n"）;

for（i=0;i

{

fprintf（p,"{"）;

for（j=0;j

fprintf（p,"%.12le",a[i][N-1]）;

fprintf（p,"},\n"）;

}

fprintf（p,"}\n\n"）;

fclose（p）;

printf（"矩阵拟上三角化完成。

\n"）;

}

else

{

printf（"Matrix（A）=\n{\n"）;

for（i=0;i

{

printf（"{"）;

for（j=0;j

printf（"%.12le\n",a[i][N-1]）;

printf（"},\n"）;

}

printf（"}\n\n"）;

}

}

voideigenvalout（intopmode）//输出矩阵的特征值；opmode=1:

到文件，0:

到屏幕；

{

inti,j;

FILE*p;

if（opmode）

{

p=fopen（"Result.txt","at+"）;

for（i=0;i

fclose（p）;

printf（"矩阵特征值计算完毕。

\n"）;

}

else

{

for（i=0;i

printf（"特征值%d=（%.12le,%.12le）\n",i+1,sa_re[i],sa_im[i]）;

}

}

voideigenvecout（intopmode）//输出矩阵的特征值；opmode=1:

到文件，0:

到屏幕；

{

inti,j;

FILE*p;

if（opmode）

{

p=fopen（"Result.txt","at+"）;

for（i=0;i

{

if（!

sa_im[i]）

{

fprintf（p,"\n"）;

fprintf（p,"特征向量%d:

",i）;

for（j=0;j

}

}

fclose（p）;

printf（"矩阵特征向量计算完毕。

\n"）;

}

else

{

for（i=0;i

{

if（!

sa_im[i]）

{

printf（"\n"）;

printf（"特征向量%d:

",i+1）;

for（j=0;j

}

}

}

}

三、程序运行结果

计算结果：

Matrix（A）=

{

{-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-00