北航数值分析第二次大作业

数值分析第二次大作业1、 算法的设计方案首先构造矩阵 A，利用Householder矩阵对矩阵A作相似变换，把A化为拟上三角矩阵A(n-1)，其算法如书上第61页所示。使用QR分解法对矩阵A(n-1)进行QR分解，其算法如书上第59页所示，进而可以求出所得矩阵的Q、R、RQ阵。然后对A(n-1)进行带双步位移的QR分解求解矩阵的全部特征值。这里采用这么几步进行：第一步：判断是否，若不是，则进入第四步。若是，则为特征值，m=m-1，若m=1，进入第二步，若m=2进入第三步，否则进入第四步。第二步：m=1，则为特征值，转向结束步。第三步：m=2，则可以求出A的两个特征值和，转向结束步。第四步：判断是否，若不是，进入第五步。若是，则得到A的两个特征值和，令m=m-2，若m=1，进入第二步，若m=2进入第三步，否则进入第一步。第五步：判断是否达到循环上限，若达到上限，则结束，否则进入第六步。第六步：对A进行双步移位QR分解，这里的算法如书上第64页所示，分解之后转向计数步。计数步：对循环次数进行计数，并转向第一步。结束步：显示所求得的特征值。最后对实特征值利用列主元高斯消元法求解其对应的特征向量，其算法见书第17页。2、 对程序的几点说明1、 对A求拟上三角化矩阵A(n-1)，这里采用一个单独的函数NSSjiao( )来求出。2、 对矩阵A(n-1)进行QR分解采用一个单独的函数QRfenjie( )求出。3、 对矩阵A(n-1)进行带双步位移的QR分解采用一个单独的函数SBQR( )求出。4、 对实特征值利用列主元高斯消元法求解其对应的特征向量在这里采用一个单独的函数vector( )来实现。5、 所有的结果都用C语言文件指令使之显示在一个txt文件中，方便结果的复制与粘贴。A(n-1)输出到NSSjiao.txt中，A(n-1) QR分解所得矩阵的Q、R、RQ阵输出到q.txt中，A的特征值输出到tezhenzhi.txt中，A实特征值对应的特征向量输出到vector.txt中。三、源程序如下：#include stdio.h#include math.h#define N 10+1#define n 10#define L 1000double aNN;double ReN=0;double ImN=0;int sgn(double a); /符号函数void NSSjiao(); /拟上三角化函数void QRfenjie(); /QR分解函数void SBQR(); /双步QR分解函数void vector(); /求特征向量的函数void main() int i,j;for(i=1;i=n;i+)for(j=1;j=0)return(1); else return(-1);void NSSjiao() int r,i,j; int cnt=0; double dr=0; double cr=0; double hr=0; double tr=0; double urN=0; double prN=0; double qrN=0; double wrN=0;for(r=1;r=n-2;r+)for(i=r+2;i=n;i+)if(air=0)cnt+; if(cnt!=n-r-1)for(i=r+1;i=n;i+) dr+=air*air; dr=sqrt(dr);cr=-sgn(ar+1r)*dr;hr=cr*cr-cr*ar+1r; for(i=1;i=r;i+) uri=0; for(i=r+2;i=n;i+) uri=air; urr+1=ar+1r-cr; for(i=1;i=n;i+) for(j=1;j=n;j+) pri+=aji*urj/hr; for(i=1;i=n;i+) for(j=1;j=n;j+) qri+=aij*urj/hr; for(i=1;i=n;i+) tr+=pri*uri/hr; for(i=1;i=n;i+) wri=qri-tr*uri; for(i=1;i=n;i+) for(j=1;j=n;j+) aij-=(wri*urj+uri*prj); dr=0; cr=0; hr=0; tr=0; for(i=1;i=n;i+) pri=0; qri=0; FILE *fp;fp=fopen(NSSjiao.txt,a);fprintf(fp,A矩阵拟上三角化后为：n);for(i=1;i=n;i+)for(j=1;j=n;j+)fprintf(fp,%13.11et,aij);fprintf(fp,n);fclose(fp);void QRfenjie() int r,i,j,t; int cnt=0; double dr=0; double cr=0; double hr=0; double urN=0; double wrN=0; double prN=0; double QNN=0; double RQNN=0; for(i=1;i=n;i+)for(j=1;j=n;j+)if(i=j)Qij=1; /Q阵初始化 for(r=1;r=n-1;r+) for(i=r+1;i=n;i+) if(air=0)cnt+;if(cnt!=n-r) for(i=r;i=n;i+) dr+=air*air; dr=sqrt(dr);cr=-sgn(arr)*dr;hr=cr*cr-cr*arr; for(i=1;ir;i+) uri=0; for(i=r+1;i=n;i+) uri=air; urr=arr-cr; for(i=1;i=n;i+) for(j=1;j=n;j+) wri+=Qij*urj; for(i=1;i=n;i+) for(j=1;j=n;j+) Qij-=wri*urj/hr; for(i=1;i=n;i+) for(j=1;j=n;j+) pri+=aji*urj/hr; for(i=1;i=n;i+) for(j=1;j=n;j+) aij-=uri*prj; dr=0; cr=0; hr=0; for(i=1;i=n;i+) wri=0; pri=0; FILE *fp;fp=fopen(q.txt,a);fprintf(fp,A(n-1)矩阵QR分解后Q阵为:n);for(i=1;i=n;i+)for(j=1;j=n;j+)fprintf(fp,%13.11et,Qij);fprintf(fp,n);fclose(fp);fp=fopen(Q.txt,a);fprintf(fp,A(n-1)矩阵QR分解后R阵为：n);for(i=1;i=n;i+)for(j=1;j=n;j+)fprintf(fp,%13.11et,aij);fprintf(fp,n);fclose(fp);fp=fopen(Q.txt,a);fprintf(fp,A(n-1)矩阵QR分解后RQ阵为：n);for(i=1;i=n;i+)for(j=1;j=n;j+)for(t=1;t=n;t+)RQij+=ait*Qtj;fprintf(fp,%13.11et,RQij);fprintf(fp,n);fclose(fp);void SBQR() int i,j;int m,k;int r;int cnt;int count;double s=0,t=0;double delta=0;double dr=0;double cr=0;double hr=0;double tr=0;double urN=0;double wrN=0;double prN=0;double vrN=0;double qrN=0;double INN=0;double aaNN=0;double cNN=0;double bNN=0;double MNN=0;for(i=1;i=n;i+)for(j=1;j=n;j+) if(i!=j) aij=sin(0.5*i+0.2*j); else aij=1.5*cos(i+1.2*j); /赋值，矩阵初始化NSSjiao();m=n;count=1;goto loop1;loop1:if(fabs(amm-1)=1e-12) /降一阶的判断Rem=amm;Imm=0;m-; if(m=1) goto loop2; if(m=2) goto loop3; else goto loop4;else goto loop4;loop2: /m=1Rem=a11;Imm=0;goto loopend;loop3: /m=2 s=am-1m-1+amm; t=am-1m-1*amm-amm-1*am-1m; delta=s*s-4*t; if(delta=0) Rem=s/2; Imm=sqrt(fabs(delta)/2; m-; Rem=s/2; Imm=sqrt(fabs(delta)/2; m-; else Rem=s/2+sqrt(fabs(delta)/2; Imm=0; m-; Rem=-s/2-sqrt(fabs(delta)/2; Imm=0; m-; goto loopend;loop4:if(fabs(am-1m-2)=1e-12) /降二阶的判断 s=am-1m-1+amm; t=am-1m-1*amm-amm-1*am-1m; delta=s*s-4*t;if(delta=0) Rem=s/2; Imm=sqrt(fabs(delta)/2; m-; Rem=s/2; Imm=-sqrt(fabs(delta)/2; m-; else Rem=s/2+sqrt(fabs(delta)/2; Imm=0; m-; Rem=s/2-sqrt(fabs(delta)/2; Imm=0; m-; if(m=1)goto loop2;else if(m=2)goto loop3;else goto loop1;else goto loop5;loop5: if(count=L)goto loopend;else goto loop6;loop6: s=am-1m-1+amm;t=am-1m-1*amm-amm-1*am-1m;for(i=1;i=m;i+)for(j=1;j=m;j+)if(i=j)Iij=1;for(i=1;i=m;i+)for(j=1;j=m;j+)for(k=1;k=m;k+)aaij+=aik*akj;for(i=1;i=m;i+)for(j=1;j=m;j+)Mij=aaij-s*aij+t*Iij;bij=Mij;cij=aij;for(i=1;i=m;i+)for(j=1;j=m;j+)aaij=0;for(r=1;r=m-1;r+) for(i=r+1;i=m;i+) if(bir=0) cnt+; if(cnt!=m-r) for(i=r;i=m;i+) dr+=bir*bir;dr=sqrt(dr);cr=-sgn(brr)*dr;hr=cr*cr-cr*brr; for(i=1;ir;i+)uri=0;for(i=r+1;i=m;i+)uri=bir;urr=brr-cr;for(i=1;i=m;i+)for(j=1;j=m;j+)vri+=bji*urj/hr;for(i=1;i=m;i+)for(j=1;j=m;j+)bij-=uri*vrj;for(i=1;i=m;i+)for(j=1;j=m;j+)pri+=cji*urj/hr;for(i=1;i=m;i+)for(j=1;j=m;j+)qri+=cij*urj/hr;for(i=1;i=m;i+)tr+=pri*uri/hr;for(i=1;i=m;i+)wri=qri-tr*uri;for(i=1;i=m;i+)for(j=1;j=m;j+)cij-=(wri*urj+uri*prj);for(i=1;i=m;i+)for(j=1;j=m;j+)aij=cij; dr=0; cr=0; hr=0; tr=0;for(i=1;i=m;i+) pri=0; qri=0; vri=0; goto loopcount;loopcount: count+; goto loop1;loopend: FILE *fp;fp=fopen(tezhenzhi.txt,a);fprintf(fp,矩阵A的特征值为:n);for(i=1;i=0)fprintf(fp,%13.11e+%13.11e*in,Rei,Imi);elsefprintf(fp,%13.11e%13.11e*in,Rei,Imi);fclose(fp);void vector() /列主元高斯消元法求解特征向量int count;int k,i,j,ik;double m=0;double xN=0;double temp=0;double sum=0;for(count=1;count=n;count+)if (Imcount=0)for (i=1;i=n;i+)for (j=1;j=n;j+)if(i!=j)aij=sin(0.5*i+0.2*j);else aij=1.5*cos(i+1.2*j)-Recount;for (k=1;k=n;k+)ik=k;for(i=k;ifabs(aikk)ik=i;for (j=k;j=n;j+)temp=akj;akj=aikj;aikj=temp;for (i=k+1;i=n;i+)m=aik/akk;for (j=k+1;j=1;k-)temp=0;for (j=k+1;j=n;j+)temp+=akj*xj;xk=-temp/akk;FILE *fp;fp=fopen(vector.txt,a);fprintf(fp,属于特征值%13.11e的特征向量为n,Recount);for(i=1;i=n;i+)fprintf(fp,%13.11et,xi);fprintf(fp,n,0);fclose(fp);3、 运行结果如下对A进行拟上三角化得到的矩阵A(n-1):-8.82751675883e-001-9.93313649183e-002-1.10334928599e+000-7.60044358564e-0011.54910107991e-001-1.94659186287e+000-8.78243638293e-002-9.25588938718e-0016.03259944053e-0011.51886095647e-001-2.34787836242e+0002.37237010494e+0001.81929082221e+0003.23780410155e-0012.20579844032e-0012.10269266255e+0001.81613808610e-0011.27883908999e+000-6.38057812440e-001-4-0010.00000000000e+0001.72827459997e+000-1+000-1.24383926270e+000-6.39975834174e-001-2.00283307904e+0002.92494720612e-001-6.41283006839e-0019.78399762128e-0022.55776357416e-0010.00000000000e+0000.00000000000e+000-1.29166953413e+000-1.11160351340e+0001+000-1.30735603002e+0001.80369917775e-001-4.24638535837e-0017.98895523930e-0021.60881992807e-0010.00000000000e+0000.00000000000e+000-8.29733316381e-0171.56012629853e+0008.12504939751e-0014.42175683292e-001-3.58861612814e-0024.69174231367e-001-2.73659505009e-001-7.35933465775e-0020.00000000000e+0000.00000000000e+0009.51733763689e-0170.00000000000e+000-7.70777375519e-001-1.58305142574e+000-3.04284317680e-0012.52871244603e-001-6.70992540145e-0012.54461992908e-0010.00000000000e+0000.00000000000e+000-1.57658260703e-0170.00000000000e+000-6.45235846333e-017-7.46345345694e-001-2.70836515702e-002-9.48652189368e-0011.19587108150e-0011.92926561795e-0020.00000000000e+0000.00000000000e+000-1.12848923378e-0160.00000000000e+0002.46927074256e-0171.62316355137e-016-7.70180137436e-001-4.69762399062e-0014.98825946801e-0011-0010.00000000000e+0000.00000000000e+0008.54810208228e-0170.00000000000e+0009.60328530894e-0186.30104169532e-0170.00000000000e+0007.01316709211e-0011.58218068848e-0013.86259461423e-0010.00000000000e+0000.00000000000e+000-2.72921948474e-0170.00000000000e+000-1.16926600276e-016-3.56137515328e-0170.00000000000e+0000.00000000000e+0004.84380760278e-0013.99277799518e-001A(n-1)矩阵QR后Q阵为：-3.51926257953e-0014.42759198224e-001-6.95598251361e-0016.48620075365e-0023.70971886190e-0011.85584714361e-001-1.62894231963e-002-1-001-5.25537538372e-002-5.48658294357e-002-9.36027728736e-001-1.66467918654e-0012.61529954856e-001-2.43867172893e-002-1.39477436089e-001-6.97758539124e-0026.12447214296e-0034.44050544314e-0021.97590790973e-0022.06283697053e-0020.00000000000e+000-8.81052055469e-001-3.98976279696e-0013.72030872848e-0022.12779406409e-0011.06446355722e-001-9.34317107976e-003-6.77420046453e-002-3.01434069867e-002-3-0020.00000000000e+0000.00000000000e+000-5.37180680644e-001-1.23494585420e-001-7.06315160872e-001-3.53345636850e-0013.10143894826e-0022.24867649160e-0011.00060178353e-0011.04462274870e-0010.00000000000e+0000.00000000000e+000-3.45070233422e-0179.89223546862e-001-1.23941173121e-001-6.20035858983e-0025.44227283946e-0033.94588163723e-0021.75581335001e-0021.83305946291e-0020.00000000000e+0000.00000000000e+0003.95807888521e-017-4.28342438788e-0175.32361069026e-001-6.73390034490e-0015.91058120587e-0024.28542532387e-0011.90690134319e-0011.99079449530e-0010.00000000000e+0000.00000000000e+000-6.55670584123e-0187.09565284549e-0184.96214253920e-017-6.05976150575e-001-9.16578303282e-002-6.64558650897e-001-2.95711087758e-001-3.08720746256e-0010.00000000000e+0000.00000000000e+000-4.69317111444e-0175.07893960458e-0171.91367990217e-0172.69244300483e-0169.93339662512e-001-9.69044031194e-002-4.31199058447e-002-4.50169441118e-0020.00000000000e+0000.00000000000e+0003.55499233623e-017-3.84720499853e-017-3.40472688637e-017-7.54341472248e-017-5.42101159583e-0185.41008800606e-001-5.81754083823e-001-6.07348058054e-0010.00000000000e+0000.00000000000e+000-1-0171.22832726408e-0178.95117787782e-0171.01897981674e-016-1.04758476565e-0171.02297400980e-017-7.22159133673e-0016.91726958888e-001A(n-1)矩阵QR分解后R阵为：2.50834274492e+000-2+000-1.31460907079e+000-3.55878749384e-002-2.60985785039e-001-1.28312184709e+000-1.39087861061e-001-8.71289797216e-0013.84936790297e-0013.35380289967e-0010.00000000000e+000-1.96160327785e+0002.40752372763e-0017.05471457282e-0015.95720431828e-0015.52697877468e-001-3.26820992441e-001-5.76949866836e-0022.87112933019e-001-8.89512875419e-0020.00000000000e+0000.00000000000e+0002.40453460199e+0001.70675809633e+000-4.23956670409e-0013.40533230581e+000-1.05001765585e-0011.46225710273e+000-6.68448746928e-001-4.02764620966e-0010.00000000000e+0000.00000000000e+0000.00000000000e+0001.57712208072e+0006.39953513396e-0013.46812787243e-001-5.70178664977e-0024.01478805443e-001-2.22247617631e-001-6.31705923644e-0020.00000000000e+0000.00000000000e+0000.00000000000e+0000.00000000000e+000-1.44784699777e+000-1.41572400774e+000-2.80613904467e-001-2.81791052189e-001-4.61143488185e-0021.99662907996e-0010.00000000000e+0000.00000000000e+0000.00000000000e+0000.00000000000e+0000.00000000000e+0001.23164145154e+0001.61970100342e-0011.96263827550e-0015.35003562176e-001-1.50927342477e-0010.00000000000e+0000.00000000000e+0000.00000000000e+0000.00000000000e+0000.00000000000e+0000.00000000000e+000-7.75344191421e-001-3.46451450882e-0014.31222680350e-0011.23464369624e-0010.00000000000e+0000.00000000000e+0000.00000000000e+0000.00000000000e+0000.00000000000e+0000.00000000000e+0000.00000000000e+0001.29631294061e+000-4.28805331834e-0012.73733415816e-0010.00000000000e+0000.00000000000e+0000.00000000000e+0000.00000000000e+0000.00000000000e+0000.00000000000e+0000.00000000000e+0000.00000000000e+000-6.70739644065e-001-4.84232012188e-0010.00000000000e+0000.00000000000e+0000.00000000000e+0000.00000000000e+0000.00000000000e+0000.00000000000e+0000.00000000000e+0000.00000000000e+0001.11022302463e-0167.16832392632e-002A(n-1)矩阵QR分解后RQ阵为：1.16307441416e+0002.63267093451e+000-1.77279600327e+000-8.66889913852e-0023.30050347105e-0011.45516237121e+000-9.73065044859e-001-4.87303117465e-001-7.75641163049e-001-3.24920197911e-0011.83611506085e+0001-001-9.88038140313e-0015.58972569477e-0014.69419006710e-002-2.97847823701e-0011.61713057765e-0026.93697770252e-0011.36767057141e-0011.41909923152e-0020.00000000000e+000-2.11852015353e+000-1.87618974578e+000-5.40707194060e-0011+000-2.55032302022e+0001.69157793654e+0001.22995161326e+0001.38794777721e+0008.66750291724e-0010.00000000000e+0000.00000000000e+000-8.47199512781e-0014.38291046832e-001-1.00863219918e+000-7.95937426150e-0014.76925886558e-0014.07268308389e-0014.09639049353e-0013.36337894086e-0010.00000000000e+0000.00000000000e+0005.08466026679e-018-1.43224434245e+000-5.74228490806e-0011.21315147772e+000-3.45750862557e-001-4.74985357312e-001-3-001-4.29450701503e-0020.00000000000e+0000.00000000000e+0005.92087592796e-017-6.40755909193e-0176.55677959800e-001-9.27525097446e-0012.52907984405e-0016.90594921698e-001-2.37443067582e-002-2.42978111978e-0010.00000000000e+0000.00000000000e+0003.52718385158e-017-3.81711071673e-017-4.87340949313e-0174.69840088488e-001-2.73077600953e-0017.82129625980e-001-9.58096493640e-0027.84623984132e-0020.00000000000e+0000.00000000000e+000-7.91891382103e-0178.56983136776e-0176.39092955976e-0174.09264318023e-0161.28767905894e+000-3.57605890035e-001-4.11672540881e-0033.91426821642e-0010.00000000000e+0000.00000000000e+000-1.83485631758e-0171.98567752865e-017-2.05076157533e-0171.25440835582e-0188.70882817833e-018-3.62876050354e-0017.39898097535e-0017.24160830958e-0020.00000000000e+0000.00000000000e+000-8-0198.80504771643e-0196.41649425503e-0187.30437740079e-018-7.50942694044e-0196.07973436028e-017-5-0024.95852290988e-002矩阵A的特征值为:3.38303961744e+000+0.00000000000e+000*i-2.32349621021e+000-8.93040517720e-001*i-2.32349621021e+000+8.93040517720e-001*i1.57754855711e+000+0.00000000000e+000*i-1.48403982226e+000+0.00000000000e