数值分析一.docx

1、数值分析一数值分析大作业题目一DY1103103 郑国昆一、 算法设计1. 矩阵A的压缩存储对于矩阵A，由于A是501阶的带状方阵，上半带宽s和下半带宽r均为2，有nr+s+1，故可认为A为稀疏矩阵，为此可设置一个矩阵C(r+s+1,501)存储A的带内元素。C中第j列存放A的第j列带内元素，同时A的主对角线元素存放在A的第s+1行。2. 和计算a) 对矩阵A用幂法求出按模最大的特征值max1，即1、501其中之一。（若max1为正，则为501；反之则为1）b) 用max对矩阵A作带原点的平移，再利用幂法计算平移后矩阵按模最大的特征值max2，于是另一个待求的特征值即为max1+max23.

2、的计算为矩阵A按模最小的特征值，可以通过反幂法求得。4. 的计算可以对矩阵A进行平移后，再用反幂法求出按模最小的特征值， =+。5. 条件数（谱范数）的计算，其中、分别为A的按模最大和最小特征值。6. det(A)的计算矩阵A的行列式可先对矩阵A进行LU分解后，det（A）等于U所有对角线上元素的乘积。二、 程序开发环境及函数说明1. 开发环境程序源代码采用MATLAB的M文件编写，在WIN7 32位机上运行。MATLAB版本号为MATLAB7.12(R2011a)。2. 函数说明 matixM：矩阵A的压缩存储函数 resetM(R)：重置矩阵A的压缩存储矩阵 detaaa(M)：求矩阵的行

3、列式 myLU(M)：对矩阵进行LU分解 max3(a,b,c)：求a，b，c的最大值 powermethod(M,eps)：幂法求矩阵A的按模最大特征值 move_pm(M,beta,eps)：带原点平移的幂法求矩阵A的按模最大特征值 invpowermethod(M,eps)：反幂法求矩阵A的按模最小特征值 move_invpm(M,lam1,lam501,eps)：带原点平移的幂法求矩阵A的按模最小特征值三、 程序源代码1. main函数：%定义精度水平%eps = 1.0e-12;% 对矩阵A进行压缩存储存入矩阵M %M=matixM;%定义复位矩阵R%R=M;lam=zeros(1,

4、39);%幂法求按模最大特征值%beta=powermethod(M,eps);%原点平移幂法求按模最大特征值%beta1=move_pm(M,beta,eps);%比较上述两值得到lamda1和lam501%lam1=min(beta,beta1);disp(矩阵的最小特征值为： num2str(lam1,%10.16e) );lam501=max(beta,beta1);disp(矩阵的最大特征值为： num2str(lam501,%10.16e) );%反幂法求矩阵按模最小特征值%M=resetM(R);betas=invpowermethod(M,eps);disp(矩阵的按模最小特征

5、值为： num2str(betas,%10.16e) );%求矩阵A的条件数%condA=abs(beta/betas);disp(矩阵A的（谱范数）条件数为： num2str(condA,%10.16e) );%求矩阵的行列式%M=resetM(R);deta=detaaa(M);disp(矩阵A的行列式为： num2str(deta,%10.16e) );M=resetM(R);lam=move_invpm(M,lam1,lam501,eps);2. matixM.m：矩阵A的压缩存储函数function M=matixMb=0.16;c=-0.064;%定义压缩矩阵的行列数%r=5;s=

6、501;M=zeros(5,501);for i=1:1:s if i=2 M(1,i)=0; else M(1,i)=c; end; if i=1 M(2,i)=0; else M(2,i)=b; end; if i=(s-1) M(4,i)=b; else M(4,i)=0; end; if i=(s-2) M(5,i)=c; else M(5,i)=0; end; M(3,i)=(1.64-0.024*i)*sin(0.2*i)-0.64*exp(0.1/i);end;3. resetM(R)：重置矩阵A的压缩存储矩阵%矩阵复位函数%function M=reset(R)l.c=size

7、(R);M=zeros(l.c);M=R;4detaaa(M)：求矩阵的行列式 function det=deta(M)r,c=size(M);M=myLU(M);det=1; for j=1:c det=det*M(3,j); end5. myLU(M)：对矩阵进行LU分解 function C=myLU(C)%实现对矩阵A的LU分解，L为下三角矩阵,C为A的压缩存储矩阵%r,n=size(C);for k=1:n for j=k:min(k+2,n) a=0; for t=max3(1,k-2,j-2):k-1 a=a+C(k-t+3,t)*C(t-j+3,j); end C(k-j+3,

8、j)=C(k-j+3,j)-a; end if keps sum1 = 0; sum1=sum1+u*u; sum1=sqrt(sum1); y=u/sum1; for j=1:c sum=0; for k=1:c if abs(j-k)eps sum1 = 0; sum1=sum1+u*u; sum1=sqrt(sum1); y=u/sum1; temp=y; N=myLU(M); for i=2:c a=0; for t=max(1,i-2):i-1 a=a+N(i-t+3,t)*y(1,t); end y(1,i)=y(1,i)-a; end u(1,c)=y(1,c)/N(3,c);

9、for i=(c-1):-1:1 a1=0; for t2=i+1:min(i+2,c); a1=a1+N(i-t2+3,t2)*u(1,t2); end u(1,i)=(y(1,i)-a1)/N(3,i); end beta2=beta; beta=0; beta=beta+temp*u; sub=abs(beta-beta2)/abs(beta);endbeta=1/beta;10. move_invpm(M,lam1,lam501,eps)：带原点平移的幂法求矩阵A的按模最小特征值function deta=move_invpm(M,deta1,deta501,eps)r,c=size(

10、M);N=M;u=zeros(1,39);for k=1:39; M=resetM(N); u(k)=deta1+k*(deta501-deta1)/40; j=1:c; M(3,j)=M(3,j)-u(k); temp=invpowermethod(M,eps); deta(k)=temp+u(k); disp(与u num2str(k) num2str(u(k),%10.16e) 最接近的特征值lamdai num2str(k) 为： num2str(deta(k),%10.16e) );end四、 程序运行结果五、 分析初始向量选择对计算结果的影响矩阵的初始向量选择，对结果的影响很大，选

11、择不同的初始向量可能会得到不同阶的特征值。以幂法为例（反幂法原理相同），常见的初始向量选择有两种：1) ui=1;(i=1,2,3.n)其运行结果截图如下：2) u=0,ui=1(i为1-n任意整数)其运行结果截图如下：试验结果发现只有当i=62且i=201时，得到的结果才与第1种初始向量相同。参看对应的特征向量发现，相对较大的数正是集中于中间位置，两端的分类近似等于0。也就是说，对于第2种初始向量的选择越接近于特征向量所能得到的结果越准确，即能够保证该特征向量在参与产生初始向量时的系数a1越不可能为0。而第1种初始向量，肯定能够保证结果的准确，但迭代数次会较第2种方法多很多。实际操作中可以选择不同形式的初始向量2，利用求得的特征值，找到其中最大的即为需求的按模最大特征值。