《数值分析大作业1.docx》由会员分享，可在线阅读，更多相关《数值分析大作业1.docx（25页珍藏版）》请在三一文库上搜索。
1、数值分析 A计算实习题目：幂法和反幂法求矩阵特征值院（系）名称机械工程及自动化学院专 业名称航空宇航制造工程班级SY12073 班学号SY1207415姓名陈俊2012 年 10 月目录第一章概述11.1数值分析 A计算实习概述11.2 题目概述1第二章算法的设计方案32.1 题目分析32.2 算法设计3第三章全部源程序6第四章特征值、条件数和行列式值的输出154.1 结果输出15第五章关于迭代初始向量的选取讨论165.1 迭代初始向量的选取对结果的影响165.2 讨论16北京航空航天大学数值分析 A计算实习报告第一章 概述1.1数值分析 A计算实习概述1）目的：训练运用计算机进行科学与工程计
2、算的能力。2）要求：（1）独立进行算法设计、程序设计和上机运算，并得出正确的结果；（2）编制程序时全部采用双精度，要求按题目的要求设计输出，并执行打印。1.2 题目概述已知：设有 501501 的矩阵 a1bc babc2 cba3bcA501501cbabc499cba500bcba501其中：0.1a = (1.64 - 0.024i )sin(0.2i ) - 0.64e, b = 0.16, c = -0.064 ，矩阵 A 的特征值iili (i = 1, 2, 501)满足 l1 l2 l501，|ls|=min|li|；试求：（1） l1 ， l501 和 ls 的值；（2）A
3、的与数 m= l+ kl501 - l1最接近的特征值 l (k =1,239) ；k140ik（3）A 的（谱范数）条件数 cond(A)2 和行列式 detA；说明：（1）在所有的算法中，凡是要给出精度水平的，都取 e =10-12 ；1北京航空航天大学数值分析 A计算实习报告（2）选择算法时，应使 A 的所有零元素都不存储；（3）打印以下内容：算法的设计方案；全部源程序（要求注明主程序和每个子程序的功能）特征值和（k=1,2,39）；讨论迭代初始向量的选取对计算结果的影响，并说明原因；（4）采用 e 型输出所有计算结果，并且至少显示 12 为有效数字。2北京航空航天大学数值分析 A计算实
4、习报告第二章 算法的设计方案2.1 题目分析由题目所求问题出发，结合题目已知条件方阵 A 为实对称矩阵，并且l l l=l112501 ，根据已有知识：对非奇异的实对称矩阵 A 有 cond ( A)2ln（ l1 和 ln 分别是矩阵 A 的模为最大和模为最小的特征值）；由定义可知，lS 相当于为矩阵 A 的按模最小的特征值；数 m= l + kl- l1051直接以 l和 l为系数因k1401501子；可以发现，l1 和 l501 的求解是本题的求解关键，可以选择幂法求解 l1 和 l501 的数值，用反幂法求解方阵 A 的按模最小的特征值 lS 。由于反幂法中要求解线性方程 Auk =
5、yk -1 并且题目要求求解行列式的值 detA，考虑到精度要求问题，选择三角分解法对矩阵 A 进行 LU 分解，则行列式值 det A = 1501 U 主对角线元素；矩阵 A 中与 mk 最接近的特征值 lik 可以用反幂法求解矩阵（ A - mk I ）按模最小的特征值 l ，则 lik = l1 + uk 。2.2 算法设计1）方阵 A 的存储方阵 A 为 501 阶实对称方阵，且上下半带宽均为 2，属于大型带状矩阵，为了节省存储量，方阵 A 的带外元素不给存储，仅存 A 的带内元素。为此，设置一个二维数组 Store_C（m，n），用于存放方阵 A 的带内元素，其中m=2+2+1=5
6、，n=501。数组 C 的第 j 列存放方阵 A 的第 j 列带内元素，并使方阵 A 的主对角线元素存放在数组 C 的第 3 行中。数组 C 存放完方阵 A 的带内元素后，多出的单元取零，这些零元共有（1+2）+（1+2）=6 个。在数组 C 中检索方阵 A 的带内元素的方法是：A 的带内元素 aii =C 中的元素 ci - j +2, j ；3北京航空航天大学数值分析 A计算实习报告2）幂法的算法说明及程序流程图（1）幂法的算法说明A. 取初始向量uiB. 初始向量的2范数C. 将初始向量单位化，得到初始迭代向量D. 算法迭代E. 求出特征值（2）幂法程序流程图入口任选初始迭代向量U0h=
7、k -1u Tu( k -1) k-1yk -1 = uk -1 / hk-1uk = Ayk -1bk= y Tuk -1 kK=1,2，b k - b k -1/b k eNY出口图 1幂法程序流程图3）反幂法的算法说明和程序流程图（1）反幂法算法说明反幂法的流程与幂法的迭代流程基本相同，但是在反幂法的迭代中， mk的求解由幂法中的 uk = Ayk -1 变为通过解线性方程 Amk = yk -1 得到，并且迭代4北京航空航天大学数值分析 A计算实习报告b- bk -1 e 变为b k-1 - bk-11 e ，最后得到的特征值取倒精度判断由幂法的kbkb -1k数即为所要求的按模最小的
8、特征值的数值。（2）反幂法的程序流程图图 2反幂法程序流程图4） l1 、 l501 和 lS 的求解因为方阵 A 为实对称矩阵，并且 l1 l2 l501 ，所以矩阵 A 的按模最大的特征值必为 l1 或 l501 ，采用如下步骤求解 l1 和 l501 ： 采用幂法求解矩阵 A 中按模最大的特征值 lM 1 ； 矩阵 B = A - lM I ，采用幂法求解矩阵 B 中按模最大的特征值 lM 2 ； lM 3 = lM 1 + lM 2 ；比较 lM 1 和 lM 3 二者的数值大小，则数值大的即为l501 ，数值小的即为 l1 ；5北京航空航天大学数值分析 A计算实习报告 lS 的求解：
9、lS 为矩阵 A 的按模最小的特征值，采用反幂法求解。求解结果：本题利用上述求解步骤调用幂法函数和反幂法函数求解得出l1 ， l501 和的数值分别为l1 = -10.700113615 ， l501 = 9.7246341， lS = -0.0055579 ；5）矩阵 A 中与数 m= l + kl501 - l1, k = 1, 2, 39 最接近的特征值 l的k140ik求解根据已经求解得出 l和 l的数值，代入式 m= l + kl501 - l1中，利用1501k140VC+ 中 的for循 环 语 句 求 解 出 mk , k = 1, 2,39 的 数 值 ， 矩 阵S t o
10、r e=_AD- mk I ，即有 Store_D2i = Store_D2i-mk ，对其采用反幂法求解按模最小的特征值 l ，则与 mk 最接近的对应的特征值 lik = l1 + mk ；6）矩阵 A 的条件数 cond(A)2 和行列式 detA 的求解因为 A 是非奇异的实对称矩阵，所以有 cond ( A)2 = l1 （ l 和 l 分别是ln1n矩阵 A 的模为最大和模为最小的特征值）；由已经求解出的 l1 ， l501 和 lS 的数值（ l1 = -10.7001，l501 = 9.7246 ，lS = -0.0055579 ；），可以判断得出按模最大和按模最小的特征值分别
11、为 l1 = -10.7001 ， lS = -0.0055579 将其带入l1-10.7001cond ( A)=1925.204 ，即为所求矩阵 A 的条件数。行列式2ln0.0055579detA 的求解，对矩阵 A 先进行 LU 分解，然后利用 VC+的 for 循环语句，求取 U 的对角线元素乘积，所得数值即为即为矩阵 A 的行列式值。求解得出的 detA=2.772786141752E+118；第三章 全部源程序/ shuzhi1.cpp : 定义控制台应用程序的入口点。6北京航空航天大学数值分析 A计算实习报告/#include stdafx.h#include math.h#d
12、efinee1.0e-12#define N 501/设定精度水平const double b=0.16;const double c=-0.064;/定义常量b、c/*/* 函数：输入稀疏矩阵A的主对角线元素，将A存储于二维数组C中(为了节约存储空间)*/ /*/void Structure(double Diagonal_A , double Store_CN )int i;Store_C00=Store_C01=Store_C10=Store_C4N-1=Store_C4N-2=St ore_C3N-1=0;for (i=2;iN;i+)Store_C0i=c;for (i=1;iN;i
13、+)Store_C1i=b;for (i=0;iN;i+)Store_C2i=Diagonal_Ai;for (i=0;iN-1;i+)Store_C3i=b;for (i=0;iN-2;i+)Store_C4i=c;7北京航空航天大学数值分析 A计算实习报告/*/* 函数：初始化迭代初始向量U0 */*/void Init_U(double U)int i;for (i=0;iN;i+)Ui=1;/*/ /* 函数：幂法求出A的按模最大的特征值*/ /*/ double Mifa(double Store_CN,double U) double Max;double Last_Eigenva
14、lue;double Pre_Eigenvalue=0.0;double Error_level;double Temp_UN;int i,r;double sign;doMax=0.0;for (i=0;iN;i+)if (Maxfabs(Ui)Max=fabs(Ui);sign=Ui/Max;r=i;for (i=0;iN;i+)Ui=Ui/Max;Temp_U0=Store_C20*U0+Store_C11*U1+Store_C02*U2;8北京航空航天大学数值分析 A计算实习报告Temp_U1=Store_C30*U0+Store_C21*U1+Store_C12*U2+Store_C
15、 03*U3;Temp_UN-2=Store_C4N-4*UN-4+Store_C3N-3*UN-3+Store_C2N-2* UN-2+Store_C1N-1*UN-1;Temp_UN-1=Store_C4N-3*UN-3+Store_C3N-2*UN-2+Store_C2N-1* UN-1;for (i=2;iN-2;i+)Temp_Ui=Store_C4i-2*Ui-2+Store_C3i-1*Ui-1+Store_C2i*Ui+Store_C1i+1*Ui+1+Store_C0i+2*Ui+2;Last_Eigenvalue=Pre_Eigenvalue;Pre_Eigenvalue=
16、sign*Temp_Ur;for (i=0;ie); return Pre_Eigenvalue;/*/* 函数：用来返回个整型值中最大的*/*/int Max3(int a,int b,int c)int temp;if(a=b)temp=a;else temp=b;if(tempb)temp=b;return temp;/*/ /* 函数：将矩阵进行LU分解，并存储在原二维数组中*/ /*/ void Decomposition_LU(double Store_CN) int i, j, k, t;for (k=0; kN; k+)/对矩阵A(存储在二维数组C中)进行LU分解for (j=
17、k; j=Min2(k+2,N-1); j+)/求U矩阵,存储在数组C对应位置Store_Ck-j+2j = Store_Ck-j+2j;for (t=Max3(0,k-2,j-2); t=k-1; t+)Store_Ck-j+2j-=Store_Ck-t+2t * Store_Ct-j+2j;for (i=k+1; i=Min2(k+2,N-1); i+)/求L矩阵,存储在数组C对应位置Store_Ci-k+2k = Store_Ci-k+2k;for (t=Max3(0,i-2,k-2); te)Last_Eigenvalue=Pre_Eigenvalue;Modele_Value=0.0
18、;for (i=0;iN;i+)Modele_Value+=Ui*Ui;Modele_Value=sqrt(Modele_Value);for (i=0;iN;i+)Ui=Ui/Modele_Value;for (i=0;iN;i+)Ti=Ui;for (t=Max3(0,i-2,0);t=0;i-)Temp_Ui=Ti;for (t=i+1;t=Min2(i+2,N-1);t+)Temp_Ui-=Store_Ci-t+2t*Temp_Ut;Temp_Ui/=Store_C2i;Pre_Eigenvalue=0.0;for (i=0;iN;i+)Pre_Eigenvalue+=Ui*Temp_
19、Ui;11北京航空航天大学数值分析 A计算实习报告for (i=0;iN;i+)Ui=Temp_Ui;Error_level=fabs(1/Pre_Eigenvalue-1/Last_Eigenvalue)/fabs(1/Pre_Eigenvalue);return Pre_Eigenvalue;/*/ /* 函数：给定偏移量U_pan，返回矩阵（A - U_pan*E) */ /*/ void Change(double SN,double U_pan) int i;for (i=0;iN;i+)S2i-=U_pan;/*/*主函数*/*/int _tmain(int argc, _TCHA
20、R* argv)double Diagonal_AN;/矩阵A的主对角线元素double Store_C5N,Store_D5N;/二维数组存储矩阵Adouble UN;/定义迭代初始向量double Eigenvalue,Eigenvalue_Temp,Eigenvalue_S,Eigenvalue_1,Eigenvalue_N;/特征值double Cond_Num;/条件数double Deter_Value=1.0;/行列式值double U_inter39,Eigenvalue_inter39; /线性数列UK、及相应的特征值 int i,k;for (i=1;iEigenvalue
21、)/输出最大、最小的特征值printf(最大的特征值是：);printf(%.12en,Eigenvalue_Temp);Eigenvalue_N=Eigenvalue_Temp;printf(最小的特征值是：);printf(%.12en,Eigenvalue);Eigenvalue_1=Eigenvalue;elseprintf(最大的特征值是：);printf(%.12en,Eigenvalue);Eigenvalue_N=Eigenvalue;printf(最小的特征值是：);printf(%.12en,Eigenvalue_Temp);Eigenvalue_1=Eigenvalue_
22、Temp;Structure(Diagonal_A, Store_C); Decomposition_LU(Store_C);/将矩阵A进行三角分解（LU分解）Init_U(U);Eigenvalue_S=1/Fan_Mifa(Store_C,U);/反幂法求出A的按模最小的特征值printf(按模最小的特征值是：);printf(%.12en,Eigenvalue_S);Cond_Num=fabs(Eigenvalue/Eigenvalue_S); /计算矩阵A的条件数cond(A)2printf(矩阵A的条件数cond(A)2是：);printf(%.12en,Cond_Num);13北京
23、航空航天大学数值分析 A计算实习报告for (i=0;iN;i+)/计算矩阵A的行列式的值Deter_Value=Deter_Value*Store_C2i;printf(矩阵A的行列式值det(A)是：);printf(%.12en,Deter_Value);for (k=1;k40;k+)/计算A的与线性数列UK最接近的特征值U_interk-1=Eigenvalue_1+k*(Eigenvalue_N-Eigenvalue_1)/40;Structure(Diagonal_A, Store_D);Change(Store_D,U_interk-1);Init_U(U);Eigenvalu
24、e_interk-1=1/Fan_Mifa(Store_D,U)+U_interk-1; printf(数U%d = %.12e , 与其最接近的特征值是：%.12en,k,U_interk-1,Eigenvalue_interk-1);return 0;14北京航空航天大学数值分析 A计算实习报告第四章 特征值、条件数和行列式值的输出4.1 结果输出1) 特征值 l1 , l501, ls ；2) cond(A)2 和 detA ；3) 与 uk 最接近的特征值 lik ；15北京航空航天大学数值分析 A计算实习报告图 3 最终结果输出第五章 关于迭代初始向量的选取讨论5.1 迭代初始向量的
25、选取对结果的影响为了更好地客观研究初始迭代向量对计算结果的影响，任取不同初始迭代向量求解 l1 和 l501 ，结果分别如下：（1）取初始迭代向量 u0 = 1,1,1,1 ，用两次幂法求得 l1 和 l501 的数值分别为-10.70011361487，9.724634101479，迭代次数分别为 594 和 1292；用反幂法求得 lS 的数值为-0.00555791079423，迭代次数为 70；（2）取初始迭代向量 u0 = 10,10,10,10，用两次幂法求得 l1 和 l501 的数值分别为-10.7001136149，9.72463410148，迭代次数分别为 594 和 12
26、92；用反幂法求得 lS 的数值为-0.005592461437513，迭代次数为 70；（3）取初始迭代向量 u0 = 100,100,100,100 ，用两次幂法求得 l1 和 l501的数值分别为-10.7001136149，9.72463410148，迭代次数分别为 594 和 1292；用反幂法求得 lS 的数值为-0.01454217452462，迭代次数为 76；（4）取初始迭代向量 u0 = 139,139,139,139 ，用两次幂法求得 l1 和 l501的数值分别为-10.7001136149，9.72463410148，迭代次数分别为 594 和 1292；用反幂法求得
27、 lS 的数值为-0.02869846741679，迭代次数为 82；（5）取初始迭代向量 u0 = 277,277,277,277 ，用两次幂法求得 l1 和l501 的数值分别为-10.7001136149，9.72463410148，迭代次数分别为 594 和1292；用反幂法求得 lS 的数值为-0.003610108303251，迭代次数为 26；5.2 讨论（1）对比分析关于不同初始迭代向量对 l1 和 l501 的求解影响，可以发现，16北京航空航天大学数值分析 A计算实习报告初始迭代向量的不同，不会影响用幂法对 l1 和 l501 数值的求解精度和迭代次数，因此验证了算法的初始
28、迭代向量可以任取的正确性。（2）对比分析不同初始迭代向量对 lS 的求解影响，可以发现，初始迭代向量的不同，影响着用反幂法对按模最小的特征值 lS 数值的求解精度和迭代次数，这与反幂法迭代算法中指出的非零初始迭代向量 U 可以任取的前提条件是相矛盾的，理论上不同的初始迭代向量 U 都不会影响着特征值数值的求解。如上情况出现的可能原因分析如下：（a）由于计算机的字长有限，所以参加运算的数据及其运算结果在计算机上存放会产生误差，反幂法主要是针对按模最小的特征值的求解，也就是求距离 0 最近的特征值，反幂法求得的特征值取倒数才是矩阵所求的按模最小的特征值，而在取倒数的过程中很有可能产生舍入误差，特殊地，矩阵 B = A - mk I 中按模最小的特征