北航数值分析-实习大作业2.doc

数值分析计算实习作业二1. 第二题 题目试用带双步位移的QR分解法求矩阵的全部特征值，并对其中的每一个实特征值求相应的特征向量。已知说明：1. 求矩阵特征值时，要求迭代的精度水平为2. 打印内容a) 算法的设计方案b) 全部源程序（要求注明主程序和每个子程序的功能）c) 矩阵经过拟上三角化后所得的矩阵d) 对矩阵进行QR分解方法结束后所得的矩阵Q、R以及乘积矩阵RQe) 矩阵的全部特征值，其中，如果是实数，则令f) 的相应于实特征值的特征向量3. 采用e型输出数据，并且至少显示12位有效数字2. 解答2.1算法的设计方案对矩阵A进行拟上三角化运算，得到一个拟上三角矩阵，再使用带双步位移的QR分解法就可以求出全部的特征值。对矩阵直接进行QR分解法，就可以得到矩阵Q、R及RQ的乘积。需要特别说明的是求矩阵实特征值对应的的特征向量：因为矩阵无重复特征值（观察前面的结果可知）。求某一特征值对应的特征向量时，的秩为n-1，对其进行“选主元的Gauss分解法”，当某一列选不出主元时（即k,k+1,n元素均为0），例如下图所示，记下列号，跳过此列，继续Gauss进行消元。消元完毕后，按从后向前的顺序求解解向量，当求解到上述记下的列号对应的解时，将其值赋值为1，直到计算完毕，再归一化就能得到实特征值对应的特征向量。2.2 源程序（见附录A）2.3 输出结果矩阵经过拟上三角化后所得的矩阵：-8.82751675883e-001 -9.93313649183e-002 -1.10334928599e+000-7.60044358564e-001 1.54910107991e-001 -1.94659186287e+000-8.78243638293e-002 -9.25588938718e-001 6.03259944053e-0011.51886095647e-001-2.34787836242e+000 2.37237010494e+000 1.81929082221e+0003.23780410155e-001 2.20579844032e-001 2.10269266255e+0001.81613808610e-001 1.27883908999e+000 -6.38057812440e-001-4-001-1.05568548241e-016 1.72827459997e+000 -1+000-1.24383926270e+000 -6.39975834174e-001 -2.00283307904e+0002.92494720612e-001 -6.41283006839e-001 9.78399762128e-0022.55776357416e-001-5.39338381277e-017 0.00000000000e+000 -1.29166953413e+000-1.11160351340e+000 1+000 -1.30735603002e+0001.80369917775e-001 -4.24638535837e-001 7.98895523930e-0021.60881992807e-0011.53346499662e-017 0.00000000000e+000 8.29733316381e-0171.56012629853e+000 8.12504939752e-001 4.42175683292e-001-3.58861612814e-002 4.69174231367e-001 -2.73659505009e-001-7.35933465775e-0021.30056273723e-016 0.00000000000e+000 -9.51733763689e-017-1.79893452016e-017 -7.70777375519e-001 -1.58305142574e+000-3.04284317680e-001 2.52871244603e-001 -6.70992540145e-0012.54461992908e-0011.61021672477e-016 0.00000000000e+000 1.57658260703e-0178.87400904240e-017 -6.45235846333e-017 -7.46345345694e-001-2.70836515702e-002 -9.48652189368e-001 1.19587108150e-0011.92926561795e-0021.36855018620e-016 0.00000000000e+000 1.12848923378e-016-1.05729133787e-017 2.46927074256e-017 0.00000000000e+000-7.70180137436e-001 -4.69762399062e-001 4.98825946801e-0011-001-2.78085130072e-017 0.00000000000e+000 -8.54810208228e-017-4.33818392824e-017 9.60328530894e-018 0.00000000000e+0000.00000000000e+000 7.01316709211e-001 1.58218068848e-0013.86259461423e-001-2.12460444005e-017 0.00000000000e+000 2.72921948474e-0172.33520544800e-016 -1.16926600276e-016 0.00000000000e+0000.00000000000e+000 0.00000000000e+000 4.84380760278e-0013.99277799518e-001对矩阵进行QR分解方法结束后所得的矩阵Q、R以及乘积矩阵RQ矩阵Q:-3.51926257953e-001 4.42759198224e-001 -6.95598251361e-0016.48620075365e-002 3.70971886190e-001 1.85584714361e-001-1.62894231963e-002 -1-001 -5.25537538372e-002-5.48658294357e-002-9.36027728736e-001 -1.66467918654e-001 2.61529954856e-001-2.43867172893e-002 -1.39477436089e-001 -6.97758539124e-0026.12447214296e-003 4.44050544314e-002 1.97590790973e-0022.06283697053e-002-4.20869709511e-017 -8.81052055469e-001 -3.98976279696e-0013.72030872848e-002 2.12779406409e-001 1.06446355722e-001-9.34317107976e-003 -6.77420046453e-002 -3.01434069867e-002-3-002-2-017 2.39575977129e-017 -5.37180680644e-001-1.23494585420e-001 -7.06315160872e-001 -3.53345636850e-0013.10143894826e-002 2.24867649160e-001 1.00060178353e-0011.04462274870e-0016.11345877564e-018 -6.81170463132e-018 3.85313956281e-0179.89223546862e-001 -1.23941173121e-001 -6.20035858983e-0025.44227283946e-003 3.94588163723e-002 1.75581335001e-0021.83305946291e-0025-017 -5.77714472780e-017 -5.44927201024e-0182-017 5.32361069026e-001 -6.73390034490e-0015.91058120587e-002 4.28542532387e-001 1.90690134319e-0011.99079449530e-0016.41944458360e-017 -7-017 4.88146744014e-0173.68833948323e-017 5.57266709004e-018 -6.05976150575e-001-9.16578303282e-002 -6.64558650897e-001 -2.95711087758e-001-3.08720746256e-0015.45599355978e-017 -6.07914733108e-017 8.28474669876e-017-6.79371484923e-017 -1.06190192156e-016 -2.47873538515e-0169.93339662512e-001 -9.69044031194e-002 -4.31199058447e-002-4.50169441118e-002-1.10864087707e-017 1.23526378012e-017 -4.28478929479e-0171.30872391235e-017 1.87794036527e-017 1.19276965418e-0161.97426148791e-017 5.41008800606e-001 -5.81754083823e-001-6.07348058054e-001-8.47015203309e-018 9.43756651499e-018 5.77456427978e-0181.37405591234e-016 1.45211729487e-016 9.91985492303e-017-4.51779279809e-017 -4.98694910963e-017 -7.22159133673e-0016.91726958888e-001矩阵R:2.50834274492e+000 -2+000 -1.31460907079e+000-3.55878749384e-002 -2.60985785039e-001 -1.28312184709e+000-1.39087861061e-001 -8.71289797216e-001 3.84936790297e-0013.35380289967e-0010.00000000000e+000 -1.96160327785e+000 2.40752372763e-0017.05471457282e-001 5.95720431828e-001 5.52697877468e-001-3.26820992441e-001 -5.76949866836e-002 2.87112933019e-001-8.89512875419e-0020.00000000000e+000 0.00000000000e+000 2.40453460199e+0001.70675809633e+000 -4.23956670409e-001 3.40533230581e+000-1.05001765585e-001 1.46225710273e+000 -6.68448746928e-001-4.02764620966e-0010.00000000000e+000 0.00000000000e+000 0.00000000000e+0001.57712208072e+000 6.39953513396e-001 3.46812787243e-001-5.70178664977e-002 4.01478805443e-001 -2.22247617631e-001-6.31705923644e-0020.00000000000e+000 0.00000000000e+000 0.00000000000e+000-1.87964619505e-016 -1.44784699777e+000 -1.41572400774e+000-2.80613904467e-001 -2.81791052189e-001 -4.61143488185e-0021.99662907996e-0010.00000000000e+000 0.00000000000e+000 0.00000000000e+000-9.40323553330e-017 -1.76632666819e-016 1.23164145154e+0001.61970100342e-001 1.96263827550e-001 5.35003562176e-001-1.50927342477e-0010.00000000000e+000 0.00000000000e+000 0.00000000000e+0008.25355059785e-018 1.55036704937e-017 1.27923107856e-017-7.75344191421e-001 -3.46451450882e-001 4.31222680350e-0011.23464369624e-0010.00000000000e+000 0.00000000000e+000 0.00000000000e+0005.98417879933e-017 1.12408272270e-016 9.27497494445e-0170.00000000000e+000 1.29631294061e+000 -4.28805331834e-0012.73733415816e-0010.00000000000e+000 0.00000000000e+000 0.00000000000e+0002.66280187565e-017 5.00187190719e-017 4.12711944393e-0170.00000000000e+000 -7.67971196511e-017 -6.70739644065e-001-4.84232012188e-0010.00000000000e+000 0.00000000000e+000 0.00000000000e+0002.77995048618e-017 5.22192671087e-017 4.30868995910e-0170.00000000000e+000 -8.01757697648e-017 -1.11022302463e-0167.16832392632e-002矩阵RQ乘积：1.16307441416e+000 2.63267093451e+000 -1.77279600327e+000-8.66889913852e-002 3.30050347105e-001 1.45516237121e+000-9.73065044859e-001 -4.87303117465e-001 -7.75641163049e-001-3.24920197911e-0011.83611506085e+000 1-001 -9.88038140313e-0015.58972569477e-001 4.69419006710e-002 -2.97847823701e-0011.61713057765e-002 6.93697770252e-001 1.36767057141e-0011.41909923152e-0021.19937255680e-016 -2.11852015353e+000 -1.87618974578e+000-5.40707194060e-001 1+000 -2.55032302022e+0001.69157793654e+000 1.22995161326e+000 1.38794777721e+0008.66750291724e-0019.22687822235e-018 -1.02807218346e-017 -8.47199512781e-0014.38291046832e-001 -1.00863219918e+000 -7.95937426150e-0014.76925886558e-001 4.07268308389e-001 4.09639049353e-0013.36337894086e-001-1.16824198777e-016 1.30167220401e-016 1.89831796776e-017-1.43224434245e+000 -5.74228490806e-001 1.21315147772e+000-3.45750862557e-001 -4.74985357312e-001 -3-001-4.29450701503e-0028.03128040945e-017 -8.94856937261e-017 4.41719790370e-017-1.57729087274e-016 6.55677959800e-001 -9.27525097446e-0012.52907984405e-001 6.90594921698e-001 -2.37443067582e-002-2.42978111978e-001-7.45016324353e-017 8.30108017939e-017 -8.87484776266e-0173-017 5.75545677225e-017 4.69840088488e-001-2.73077600953e-001 7.82129625980e-001 -9.58096493640e-0027.84623984132e-0027.31620980977e-017 -8-017 9.52044873069e-0174.77395921795e-017 -1.12781945424e-016 -4.35885624255e-0161.28767905894e+000 -3.57605890035e-001 -4.11672540881e-0033.91426821642e-0011-017 -1.28553757058e-017 1.16394943443e-017-2.91230374876e-017 -8.59482333643e-017 -1.68340747666e-016-6.41137456972e-017 -3.62876050354e-001 7.39898097535e-0017.24160830958e-002-6.07167934784e-019 6.76515338557e-019 -1.45194174703e-0175.80731181685e-017 7.23970610235e-018 -3.49640310258e-017-7.91872099409e-017 -2.90931754028e-017 -5-0024.95852290988e-002矩阵的全部特征值3.38303961744e+000 +0.00000000000e+000i-2.32349621021e+000 +8.93040517720e-001i-2.32349621021e+000 -8.93040517720e-001i1.57754855711e+000 +0.00000000000e+000i-1.48403982226e+000 +0.00000000000e+000i-9.80530956290e-001 +1-001i-9.80530956290e-001 -1-001i9.35588907819e-001 +0.00000000000e+000i6.36062787575e-001 +0.00000000000e+000i5.65048899350e-002 +0.00000000000e+000i的相应于实特征值的特征向量特征值3.38303961744e+000对应的特征向量为： -1.05221573436e-001 -2-001 -4.73017877776e-001 -2.60887478870e-001 -3.05805625140e-001 -2.58379783221e-001 8.73379363950e-002 4.05454033853e-001 5.09013113720e-001 2.40925044562e-001特征值1.57754855711e+000对应的特征向量为： 6.24962559917e-002 -1.12090045150e-002 -2.49666705179e-001 -1.31364486150e-001 -3.83541256120e-001 8-001 -1.24560352067e-001 -6.83561054721e-002 2.70587401777e-001 1.00519108146e-001特征值-1.48403982226e+000对应的特征向量为： -5.60118116800e-001 7.79341508770e-001 1.33780114857e-002 -2.77409287891e-001 3.00557588301e-003 -2.53483408960e-003 -2.06284849088e-002 -1.10134829107e-002 -1.22486166525e-002 3.23620930409e-002特征值9.35588907819e-001对应的特征向量为： 8.05651548269e-002 4.61113468803e-002 -1.50243791976e-002 -4.81208346184e-002 -3.53633892069e-001 2.08919077811e-001 -1.55745966395e-001 8.19670427271e-001 -3.50982512819e-001 2.88513852783e-002特征值6.36062787575e-001对应的特征向量为： 1.07037468397e-001 7.12345583256e-002 3.90237917138e-001 -4.46655513756e-002 -7.19034744275e-001 1.75815688440e-001 -2.26537961544e-001 3.76886150502e-001 2.95625408505e-001 2.25577972559e-002特征值5.65048899350e-002对应的特征向量为： -2.09017518582e-001 -2.00063796685e-001 3.89176061817e-001 -2.77939118707e-002 -3.93236547906e-001 -1.24703810952e-001 6.44810939558e-001 -3.02779853108e-001 -2.91095263314e-001 4.09436555812e-0023. 附录A （源程序）#include /*在编译环境gcc和VS2012下调试通过*/#include #include #include void QR(double *A,int n,double *Q,double *R);/*QR分解法*/int SBQR(double *A,int n,double *lambda);/*带双步位移的QR分解法*/void Hessenberg(double *A,int n,double *a);/*拟上三角化函数*/void vector(double *A,int n,double labmda,double *p);/*求非重根特征值lambda的特征向量*/void Print(double *A,int n);/*输出函数*/double *dmatrix(int m,int n);/*生成坐标从1开始的矩阵*/void free_dmatrix(double *a);/*释放由dmatrix申请的空间*/int main()int i,j,k,n=10;double *A,*Q,*R,*A_n_1,*RQ;double *lambda,*kesi;A=dmatrix(n,n);Q=dmatrix(n,n);R=dmatrix(n,n);RQ=dmatrix(n,n);A_n_1=dmatrix(n,n);lambda=dmatrix(n,2);kesi=dmatrix(n,n);/*初始化A*/for(i=1;i=n;i+)for(j=1;j=n;j+)if(i != j)Aij=sin(0.5*i+0.2*j);elseAij=1.5*cos(i+1.2*j);/*求A经过拟上三角化后得到的矩阵A(n-1)*/Hessenberg(A,n,A_n_1);printf(对A进行拟上三角化后得到的矩阵为：n);printf(A_n_1:n);Print(A_n_1,n);/*对A(n-1)矩阵进行QR后的Q、R矩阵以及矩阵乘积RQ*/QR(A_n_1,n,Q,R);printf(对上三角矩阵A(n-1)进行QR分解的结果：n);printf(矩阵Q:n);Print(Q,n);printf(矩阵R:n);Print(R,n);printf(矩阵乘积RQ：n);for(i=1;i=n;i+)/*计算RQ矩阵的乘积*/for(j=1;j=n;j+)RQij=0;for(k=1;k求出了全部特征值：n);printf(特征值：n);for(i=1;i未求出全部特征值：n);/*求实特征值对应的特征向量*/printf(全部实特征值对应的特征向量为：n);for(i=1;i=n;i+)if(fabs(lambdai2) 1e-12)/*确定特征值为实数*/vector(A,n,lambdai1,kesii);printf(特征值%.11e对应的特征向量为：n,lambdai1);for(k=1;k只能求非重实数根的特征向量！by lidezheng 日期：2013-11-1-*/void vector(double *A,int n,double lambda,double *x)int i,j,k,point,flag;double temp,m;double *a;a=dmatrix(n,n);/*初始化a,将A中的元素复制到a中*/for(i=1;i=n;i+)for(j=1;j=n;j+)aij=Aij;for(i=1;i=n;i+)aii-=lambda;/*对角线元素减去特征值*/*对a进行列主元素Gauss消元法*/point=n;for(k=1;kn的元素时候全为0*/for(i=k;i 1e-12)flag=0;break;if(flag = 1)point=k;continue;for(i=k,j=k;j=n;j+)/*选主元所在行号*/if(fabs(aik) fabs(ajk)i=j;if(i != k)/*不等则交换两行*/for(j=k;j=n;j+)temp=aij;aij=akj;akj=temp;for(i=k+1;i=n;i+)/*消元*/m=aik/akk;for(j=k+1;j=1;k-)if(k != point)temp=0;for(j=k+1;j=n;j+)temp-=akj*xj;xk=temp/akk;/*对特征向量归一化*/temp=0;for(i=1;i=n;i+)temp+=pow(xi,2);temp=sqrt(temp);for(i=1;i=n;i+)xi/=temp;free_dmatrix(a);/*-*函数名：SBQR(double *A,int n,double *labmda)*功 能：带双步位移的QR分解法，求全部特征值*参 数：*返回值：1-求出了全部特征值 0-求出了部分特征值*说 明：by lidezheng 日期：2013-10-31-*/int SBQR(double *A,int n,double *lambda)int i,j,l,r,L,k,m,flag,ret;double s,t,deta;double s12,s22;double dr,cr,hr,tr;double *ur,*vr,*pr,*qr,*wr;double *a,*B;a=dmatrix(n,n);/*对矩阵A进行拟上三角化运算*/Hessenberg(A,n,a);k=1;m=n;L=1000;/*双步位移QR循环*/while(1)if(m = 1 | fabs(amm-1) = 1)continue;s=am-1m-1+amm;t=am-1m-1*amm-amm-1*am-1m;deta=s*s-4*t;/*计算二次方程*/if(deta = 0)s10=(s+sqrt(deta)/2;s20=(s-sqrt(deta)/2;s11=s21=0;elses10=s20=s/2;s11=sqrt(-deta)/2;s21=-sqrt(-deta)/2;if(m =2 | fabs(am-1m-2) = 1e-12)lambdam1=s20;lambdam2=s21;lambdam-11=s10;lambdam-12=s11;m-=2;if(m = 0)ret=1;/*出口状态*/break;elsecontinue;elseif(k = L)ret=0;break;else/*计算s,t*/s=am-1m-1+amm;t=am-1m-1*amm-amm-1*am-1m;B=dmatrix(m,m);ur=(double *)malloc(n+1)*sizeof(double);vr=(double *)malloc(n+1)*sizeof(double);pr=(double *)malloc(n+1)*sizeof(double);qr=(double *)malloc(n+1)*sizeof(double);wr=(double *)malloc(n+1)*sizeof(double);/*计算m*m矩阵B*/for(i=1;i=m;i+)for(j=1;j=m;j+)Bij=0;for(l=1;l