数值分析计算实习作业二

数值分析

计算实习题二

学号： 姓名： 院系：

2015年11月28日

数值分析计算实习作业二

一、算法

采用拟上三角化的双步位移的QR分解法计算特征值，并利用列主元素Guass消去法求实数特征值对应的特征向量，定义全局变量n=10。 1 主程序

在定义一维和二维数组时将数组的大小定为n+1（大小为11），这样在使用单个元素时的角标可以与矩阵或向量元素的角标相对应。

先输入矩阵A的各个元素；

保留数组A，定义数组B的元素与A相同，即B=A，后面的矩阵变换用数组B来进行；

对B进行拟上三角化； 输出拟上三角化后的矩阵B；

对拟上三角化后的B进行QR分解，输出矩阵Q、R、RQ； 对拟上三角化的数组矩阵B进行QR分解； 输出A的全部特征值；

判断所有特征值，如果特征值为实数（即储存特征值虚部的数组中元素为0），利用列主元素高斯消去法求其特征向量；

输出实特征值对应的特征向量。 2 拟上三角化程序

对于r=1,2，……，n-2循环执行

（1）求第r列的第r+2个元素到第n个元素的平方和，如果和为0，则跳出该次循环直接进行下一次r的循环。如果和不为0，则继续执行。 （2）计算第r列从第r+1到第n个元素的平方和的平方根

drir1(B[i][r])n2 当B[r+1][r]小于等于0时，取cr=dr，否则取cr=-dr hr=cr*cr-cr*B[r+1][r]

（3）令ur(0,...,0,B[r1][r]cr,B[r2][r],...,B[n][r])T （4）计算

1

数值分析计算实习作业二

prBT(ur)/(hr)qrB(ur)/(hr)tr(pr)(ur)/(hr)wrqr(tr)(ur)BB(wr)(ur)T(ur)(pr)T（5）执行下一次循环。

3 拟上三角化后的矩阵的QR分解

初始令Q=I。

对于r=1,2，……，n-1循环执行

（1）求第r列的第r+1个元素到第n个元素的平方和，如果和为0，则跳出该次循环直接进行下一次r的循环。如果和不为0，则继续执行。 （2）计算第r列从第r到第n个元素的平方和的平方根

dr

ir1(B[i][r])n2 当B[r+1][r]小于等于0时，取cr=dr，否则取cr=-dr hr=cr*cr-cr*B[r+1][r]

（3）令ur(0,...,0,B[r1][r]cr,B[r2][r],...,B[n][r])T （4）计算

wrQ(ur)QQ(wr)*(ur)T/(hr)prB(ur)/(hr)BB(ur)(pr)TT

（5）执行下一次循环。

循环结束后的Q和B就是分解后的Q和R，之后再求RQ。 4 QR分解过程

给定精度水平sigma=e-12，给定迭代的次数限制L，令k=1，m=n。利用goto和if判断语句在下面几个步骤中跳转。

loop3：如果fabs(B[m][m-1])loop4：如果m=1，则得到特征值B[1][1]放入特征值数组中，跳转到loop9，如果m=2，则求二阶子阵特征值放入特征值数组中，跳转到loop9，否则跳转到

2

数值分析计算实习作业二

loop3。

loop5：如果fabs(B[m-1][m-2])<=sigma，则对右下角相邻的一个二阶子阵执行求特征值函数得到两个特征值放入数组中，并且令m=m-2，跳转到loop4；否则跳转到loop6。

loop6：如果k=L，则跳出整个循环，输出”未得到全部特征值；否则跳转到loop7。

loop7：计算

sB[m1][m1]B[m][m]tB[m1][m1]*B[m][m]B[m][m1]*B[m1][m] MB*Bs*Bt*I之后对M执行QR分解，同时得到新的经过QR分解后的矩阵B，跳转到loop8。

loop8：令k=k+1，跳转到loop3。 loop9：结束整个循环。 最后输出特征值。 5对M的QR分解

对于r=1,2，……，m循环执行

（1）将M矩阵的第r列的第r+1个元素到第m个元素的绝对值相加，如果和为0，则跳出该次循环直接进行下一次r的循环。如果和不为0，则继续执行。 （2）计算矩阵M第r列从第r到第m个元素的平方和的平方根

dr(M[i][r])irm2 当M[r][r]小于等于0时，取cr=dr，否则取cr=-dr hr=cr*cr-cr*M[r][r]

（3）令ur(0,...,0,M[r][r]cr,M[r1][r],...,M[m][r])T （4）计算

3

数值分析计算实习作业二

vrMT(ur)/(hr)MM(ur)*(vr)TprBT(ur)/(hr)qrB(ur)/(hr)tr(pr)T(ur)/(hr)wrqr(tr)(ur)BB(wr)(ur)T(ur)(pr)T（5）执行下一次循环。

最终得到的矩阵B便是对M进行QR分解之后随之变动的矩阵B。 6二阶子阵的特征值

此函数的输入量为m，计算

bB[m1][m1]B[m][m]

cB[m1][m1]*B[m][m]B[m][m1]*B[m1][m]

之后判断一元二次函数的是否小于零，从而得到两个特征值的实部和虚部分别放入特征值数组的相应位置。 7列主元素高斯消去法

参照课本的步骤，需要注意的是因为AI为奇异，所以在求特征向量的元素时令最后一个元素为1，再利用课本上的步骤倒推其他元素的值。

二、程序

#include #include #include

#define L 2500//定义迭代次数 #define n 10

double A[n+1][n+1]={0},B[n+1][n+1]={0},M[n+1][n+1]={0},X[n+1][n+1]={0}; double lamr[n+1]={0},lamc[n+1]={0};//lamr储存特征值的实部，lamc储存特征值的虚部

double sigma=1e-12;//定义精度 void hess();//拟上三角化函数 void AQR();//A的QR分解函数

void QRmeth();//带双步位移的QR分解函数 void MQRfenjie(int m);//M的QR分解函数

void gauss(int y);//列主元素gauss消去法求特征向量函数

/*****************以下为主程序*****************/

4

数值分析计算实习作业二

void main() {

int i,j,k;

for(i=1;i<=n;i++)//输入矩阵A {

for(j=1;j<=n;j++) {

if(i!=j) A[i][j]=sin(0.5*i+0.2*j);

else if(i==j) A[i][j]=1.52*cos(i+1.2*j); } }

for(i=1;i<=n;i++)//令矩阵B=A，对B做拟上三角化 {

for(j=1;j<=n;j++) {

B[i][j]=A[i][j]; } }

hess();//拟上三角化矩阵B

AQR();//将拟上三角化后的矩阵B进行QR分解 QRmeth();//带双步位移的QR分解

for(k=1;k<=n;k++)//判断特征值是否为实数，并求实数的特征向量 {

if(lamc[k]!=0) continue; else {

for(i=1;i<=n;i++)//求特征向量方程的矩阵系数 {

for(j=1;j<=n;j++) {

if(i==j) A[i][j]=A[i][j]-lamr[k]; else A[i][j]=A[i][j]; } }

gauss(k);//gauss消去法求特征向量 } } }

/*****************以下为拟上三角化程序*****************/ void hess()//拟上三角化函数 {

int r=0,i=0,j=0; double temp1,sum=0;

double dr=0,cr=0,hr=0,tr=0;

5

数值分析计算实习作业二

double ur[n+1]={0},pr[n+1]={0},qr[n+1]={0},wr[n+1]={0}; for(r=1;r<=n-2;r++) {

for(temp1=0.0,i=r+2;i<=n;i++) temp1+=B[i][r]*B[i][r]; if(temp1==0) continue; else {

temp1+=B[r+1][r]*B[r+1][r]; dr=sqrt(temp1);

if(B[r+1][r]>0) cr=-dr; else cr=dr;

hr=cr*cr-cr*B[r+1][r]; for(i=1;i<=r;i++) ur[i]=0.0; ur[r+1]=B[r+1][r]-cr; for(i=r+2;i<=n;i++) ur[i]=A[i][r];

for(j=1;j<=n;j++) {

for(pr[j]=0.0,i=1;i<=n;i++) pr[j]+=B[i][j]*ur[i]/hr; }

for(tr=0.0,i=1;i<=n;i++) {

tr+=pr[i]*ur[i]/hr; }

for(i=1;i<=n;i++) {

for(qr[i]=0.0,j=1;j<=n;j++) qr[i]+=B[i][j]*ur[j]/hr; }

for(i=1;i<=n;i++) {

wr[i]=qr[i]-tr*ur[i]; }

for(i=1;i<=n;i++) {

for(j=1;j<=n;j++) B[i][j]-=(wr[i]*ur[j]+ur[i]*pr[j]); } } }

for(i=1;i<=n;i++)

6

数值分析计算实习作业二

{

for(j=1;j<=n;j++)

if(fabs(B[i][j])printf(\"拟上三角化后A(n-1):\\n\");

for(i=1;i<=n;i++)//输出拟上三角化函数 {

for(j=1;j<=n;j++) {

printf(\"%1.12e,\ if(j%3==0) printf(\"\\n\"); }

printf(\"\\n\"); } }

/*****************以下为QR分解程序*****************/ void AQR()//QR分解函数 {

double

ur[n+1],wr[n+1],pr[n+1],cr,dr,hr,R[n+1][n+1],Q[n+1][n+1]={0},rq[n+1][n+1]={0}; int i,j,r,k;

for(i=1;i<=n;i++)//令矩阵R等于矩阵B，对R进行QR分解 for(j=1;j<=n;j++) R[i][j]=B[i][j];

for(i=1;i<=n;i++) Q[i][i]=1;//Q初始为单位阵

for(r=1;r<=n-1;r++) {

for(dr=0.0,i=r;i<=n;i++) dr+=R[i][r]*R[i][r]; dr=sqrt(dr);

if (dr==0) continue; else {

if(R[r][r]>0) cr=-dr; else cr=dr;

hr=cr*cr-cr*R[r][r]; for(i=1;iur[r]=R[r][r]-cr; for(i=r+1;i<=n;i++) ur[i]=R[i][r]; for(i=1;i<=n;i++)

7

数值分析计算实习作业二

{

for (pr[i]=0.0,j=1;j<=n;j++) pr[i]+=R[j][i]*ur[j]/hr; }

for(i=1;i<=n;i++) {

for(j=1;j<=n;j++)

R[i][j]=R[i][j]-ur[i]*pr[j]; }

for(i=1;ifor(wr[i]=0.0,j=1;j<=n;j++) wr[i]+=Q[i][j]*ur[j]; }

for(i=1;i<=n;i++) {

for(j=1;j<=n;j++)

Q[i][j]-=wr[i]*ur[j]/hr; } } }

for(i=1;i<=n;i++)//求RQ的矩阵值 {

for(j=1;j<=n;j++) {

for(rq[i][j]=0.0,k=1;k<=n;k++) rq[i][j]+=R[i][k]*Q[k][j]; } }

printf(\"生成的Q阵：\\n\"); for(i=1;i<=n;i++) {

for(j=1;j<=n;j++)

if(fabs(Q[i][j])for(i=1;i<=n;i++)//输出矩阵Q {

for(j=1;j<=n;j++) {

printf(\"%1.12e,\ if(j%3==0) printf(\"\\n\"); }

printf(\"\\n\"); }

printf(\"生成的R阵：\\n\");

8

数值分析计算实习作业二

for(i=1;i<=n;i++) {

for(j=1;j<=n;j++)

if(fabs(R[i][j])for(i=1;i<=n;i++)//输出矩阵R {

for(j=1;j<=n;j++) {

printf(\"%1.12e,\ if(j%3==0) printf(\"\\n\"); }

printf(\"\\n\"); }

printf(\"生成的RQ阵：\\n\"); for(i=1;i<=n;i++) {

for(j=1;j<=n;j++)

if(fabs(rq[i][j])for(i=1;i<=n;i++)//输出矩阵RQ {

for(j=1;j<=n;j++) {

printf(\"%1.12e,\ if(j%3==0) printf(\"\\n\"); }

printf(\"\\n\"); } }

/*****************以下为带双步位移的QR分解的程序*****************/ void QRmeth()//带双步位移的QR分解 {

int i,j,k=1,m=n,s; int nR=0,nC=0;

double s1,t,x,sum,I[n+1][n+1]={0}; for(i=1;i<=n;i++) I[i][i]=1;

loop3:if(fabs(B[m][m-1])<=sigma) {

lamr[m]=B[m][m]; m--;

goto loop4; }

else goto loop5;

9

数值分析计算实习作业二

loop4:if(m==1) {

lamr[m]=B[1][1]; goto loop9; }

else if(m==2)//求二阶子阵的特征值 {

s1=B[m-1][m-1]+B[m][m];

t=B[m-1][m-1]*B[m][m]-B[m][m-1]*B[m-1][m]; x=s1*s1-4*t; if(x>=0) {

lamr[m-1]=(s1+sqrt(x))/2; lamr[m]=(s1-sqrt(x))/2; } else {

lamr[m-1]=s1/2;

lamc[m-1]=sqrt(-x)/2; lamr[m]=s1/2;

lamc[m]=-sqrt(-x)/2; }

goto loop9; }

else goto loop3; loop5:{

if(fabs(B[m-1][m-2])<=sigma) {

s1=B[m-1][m-1]+B[m][m];

t=B[m-1][m-1]*B[m][m]-B[m][m-1]*B[m-1][m]; x=s1*s1-4*t; if(x>=0) {

lamr[m-1]=(s1+sqrt(x))/2; lamr[m]=(s1-sqrt(x))/2; } else {

lamr[m-1]=s1/2;

lamc[m-1]=sqrt(-x)/2; lamr[m]=s1/2;

lamc[m]=-sqrt(-x)/2; } m--;

10

数值分析计算实习作业二

m--; goto loop4; }

else goto loop6; } loop6:{

if(k==L) printf(\"计算终止，未能得到全部特征值\\n\"); else goto loop7; } loop7:{

s1=B[m-1][m-1]+B[m][m];

t=B[m-1][m-1]*B[m][m]-B[m][m-1]*B[m-1][m];

for(i=1;i<=m;i++)//生成M for(j=1;j<=m;j++) {

for(sum=0.0,s=1;s<=m;s++) sum+=B[i][s]*B[s][j];

M[i][j]=sum-s1*B[i][j]+t*I[i][j]; }

MQRfenjie(m);//M进行QR分解，求出新的矩阵B goto loop8; }

loop8:{k++;goto loop3;} loop9:{ ; }

printf(\"矩阵的全部特征值为：\\n\");

for(i=1;i<=n;i++)//对于实数特征值和复数特征值两种不同的输出方法 {

if(lamc[i]==0) printf(\"lam%d=%13.11le\\n\

else printf(\"lam%d=(%13.11le,%13.11le)\\n\} }

/*****************以下为对M的QR分解的程序*****************/ void MQRfenjie(int m)//QR分解函数 {

int r=0,i=0,j=0; double temp1,sum=0;

double dr=0,cr=0,hr=0,tr=0;

double ur[n+1]={0},vr[n+1]={0},pr[n+1]={0},qr[n+1]={0},wr[n+1]={0}; for(r=1;r<=m-1;r++) {

11

数值分析计算实习作业二

for(temp1=0.0,i=r+1;i<=m;i++) temp1+=M[i][r]*M[i][r]; if(temp1==0) continue; else {

temp1+=M[r][r]*M[r][r]; dr=sqrt(temp1);

if(M[r][r]>0) cr=-dr; else cr=dr;

hr=cr*cr-cr*M[r][r]; for(i=1;iur[r]=M[r][r]-cr; for(i=r+1;i<=m;i++) ur[i]=M[i][r]; for(j=1;j<=m;j++)

for(vr[j]=0.0,i=1;i<=m;i++) vr[j]+=M[i][j]*ur[i]/hr; for(i=1;i<=m;i++) for(j=1;j<=m;j++) M[i][j]-=ur[i]*vr[j]; for(j=1;j<=m;j++)

for(pr[j]=0.0,i=1;i<=m;i++) pr[j]+=B[i][j]*ur[i]/hr; for(i=1;i<=m;i++)

for(qr[i]=0.0,j=1;j<=m;j++) qr[i]+=B[i][j]*ur[j]/hr; for(tr=0.0,i=1;i<=m;i++) tr+=pr[i]*ur[i]/hr; for(i=1;i<=m;i++) wr[i]=qr[i]-tr*ur[i]; for(i=1;i<=m;i++) for(j=1;j<=m;j++)

B[i][j]-=(wr[i]*ur[j]+ur[i]*pr[j]); } } }

/**************以下为列主元素高消去法求特征向量程序**************/ void gauss(int y)//列主元的高斯消元法求解特征向量 {

int k,i,j,p;

double temp1,temp2,b[n+1]={0},max,m; for(k=1;k<=n-1;k++)//换行过程 {

12

数值分析计算实习作业二

max=A[k][k];p=k; for(i=k;i<=n;i++) {

if(fabs(A[i][k])>max) p=i; }

for(j=k;j<=n;j++) {

temp1=A[k][j]; A[k][j]=A[p][j]; A[p][j]=temp1; }

for(i=k+1;i<=n;i++)//消元过程 {

m=A[i][k]/A[k][k]; for(j=k+1;j<=n;j++) {

A[i][j]=A[i][j]-m*A[k][j]; }

b[i]=b[i]-m*b[k]; } }

X[y][n]=1;//因系数矩阵奇异，故无关特征向量只有n-1个，所以取最后一个元素为1

for(k=n-1;k>=1;k--) {

temp2=0;

for(j=k+1;j<=n;j++) {

temp2+=A[k][j]*X[y][j]; }

X[y][k]=(b[k]-temp2)/A[k][k]; }

printf(\"对应特征值%1.12e的特征向量为:\\n\ for(i=1;i<=n;i++) printf(\"%1.12e\\n\}

三、结果

拟上三角化的结果

13

数值分析计算实习作业二

14

数值分析计算实习作业二

QR分解后的Q矩阵

15

数值分析计算实习作业二

QR分解后的R矩阵

16

数值分析计算实习作业二

QR分解后的RQ矩阵

矩阵的全部特征值（带括号的表示复数，逗号前表示实部，逗号后表示虚部）

17

数值分析计算实习作业二

实数特征值对应的特征向量

18

数值分析计算实习作业二

四、讨论

1、对于复数特征值的储存方法，采用两个数组，一个储存实部，一个储存虚部，因为复数特征值只在二阶子阵的特征值求解中产生，很容易判断特征值是否为复数，也很容易分别求出复数特征值的实部和虚部，因此可以将复数特征值储存在两个数组中。

2、关于特征向量的求解，需要求解AIX0这个方程组，对于之前讲的迭代法，高斯消去法等线性方程组的解法前提都是系数矩阵非奇异，而这个方程组的系数矩阵显然是奇异的（否则X只有零解了），所以在采用高斯消去法时在求解过程中先规定最后一个元素的值为1，然后再倒推出其余各元素。

19