数值分析一.docx
《数值分析一.docx》由会员分享，可在线阅读，更多相关《数值分析一.docx（12页珍藏版）》请在冰豆网上搜索。
数值分析一
数值分析大作业题目一
DY1103103郑国昆
一、算法设计
1.矩阵A的压缩存储
对于矩阵A，由于A是501阶的带状方阵，上半带宽s和下半带宽r均为2，有n>>r+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+
max2
3.
的计算
为矩阵A按模最小的特征值，可以通过反幂法求得。
4.
的计算
可以对矩阵A进行
平移后，再用反幂法求出按模最小的特征值
，
=
+
。
5.条件数（谱范数）
的计算
，其中
、
分别为A的按模最大和最小特征值。
6.det（A）的计算
矩阵A的行列式可先对矩阵A进行LU分解后，det（A）等于U所有对角线上元素的乘积。
二、程序开发环境及函数说明
1.开发环境
程序源代码采用MATLAB的M文件编写，在WIN732位机上运行。
MATLAB版本号为MATLAB7.12（R2011a）。
2.函数说明
●matixM：
矩阵A的压缩存储函数
●resetM（R）：
重置矩阵A的压缩存储矩阵
●detaaa（M）：
求矩阵的行列式
●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,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（['矩阵的按模最小特征值为：
'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的压缩存储函数
functionM=matixM
b=0.16;
c=-0.064;
%定义压缩矩阵的行列数%
r=5;
s=501;
M=zeros（5,501）;
fori=1:
1:
s
ifi<=2
M（1,i）=0;
else
M（1,i）=c;
end;
ifi==1
M（2,i）=0;
else
M（2,i）=b;
end;
ifi<=（s-1）
M（4,i）=b;
else
M（4,i）=0;
end;
ifi<=（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的压缩存储矩阵
%矩阵复位函数%
functionM=reset（R）
[l.c]=size（R）;
M=zeros（l.c）;
M=R;
4．detaaa（M）：
求矩阵的行列式
functiondet=deta（M）
[r,c]=size（M）;
M=myLU（M）;
det=1;
forj=1:
c
det=det*M（3,j）;
end
5.myLU（M）：
对矩阵进行LU分解
functionC=myLU（C）
%实现对矩阵A的LU分解，L为下三角矩阵,C为A的压缩存储矩阵%
[r,n]=size（C）;
fork=1:
n
forj=k:
min（k+2,n）
a=0;
fort=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,j）=C（k-j+3,j）-a;
end
ifkfori=k+1:
min（k+2,n）
a=0;
fort=max3（1,i-2,k-2）:
k-1
a=a+C（i-t+3,t）*C（t-k+3,k）;
end
C（i-k+3,k）=（C（i-k+3,k）-a）/C（3,k）;
end
end
end
6.max3（a,b,c）：
求a，b，c的最大值
functiont=max3（a,b,c）
temp=max（a,b）;
t=max（temp,c）;
7.powermethod（M,eps）：
幂法求矩阵A的按模最大特征值
%幂法球按模最大特征值源代码%
functionbeta=powermethod（M,eps）
[r,c]=size（M）;
i=1:
501;
u（i）=1;
sub=1;%计算误差时的中间变量%
beta=0;
beta2=0;
whilesub>eps
sum1=0;
sum1=sum1+u*u';
sum1=sqrt（sum1）;
y=u/sum1;
forj=1:
c
sum=0;
fork=1:
c
ifabs（j-k）<=2
sum=sum+M（j-k+3,k）*y（k）;
u（1,j）=sum;
end
end
end
beta2=beta;
beta=0;
fori=1:
c
beta=beta+y（1,i）*u（1,i）;
end
sub=abs（beta-beta2）/abs（beta）;
end
8.move_pm（M,beta,eps）：
带原点平移的幂法求矩阵A的按模最大特征值
functiondeta=move_pm（M,detam,eps）
[r,c]=size（M）;
forj=1:
c
M（3,j）=M（3,j）-detam;
end
deta=powermethod（M,eps）;
deta=deta+detam;
9.invpowermethod（M,eps）：
反幂法求矩阵A的按模最小特征值
%反幂法球按模最小特征值源代码%
functionbeta=invpowermethod（M,eps）
[r,c]=size（M）;
%初始化迭代向量%
i=1:
501;
u（i）=1;
sub=1;%计算误差时的中间变量%
%中间变量%
beta=0;
beta2=0;
whilesub>eps
sum1=0;
sum1=sum1+u*u';
sum1=sqrt（sum1）;
y=u/sum1;
temp=y;
N=myLU（M）;
fori=2:
c
a=0;
fort=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）;
fori=（c-1）:
-1:
1
a1=0;
fort2=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）;
end
beta=1/beta;
10.move_invpm（M,lam1,lam501,eps）：
带原点平移的幂法求矩阵A的按模最小特征值
functiondeta=move_invpm（M,deta1,deta501,eps）
[r,c]=size（M）;
N=M;
u=zeros（1,39）;
fork=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
四、程序运行结果
五、分析初始向量选择对计算结果的影响
矩阵的初始向量选择，对结果的影响很大，选择不同的初始向量可能会得到不同阶的特征值。
以幂法为例（反幂法原理相同），常见的初始向量选择有两种：
1）u[i]=1;（i=1,2,3...n）
其运行结果截图如下：
2）u=0,u[i]=1（i为1-n任意整数）
其运行结果截图如下：
试验结果发现只有当i>=62且i<=201时，得到的结果才与第1种初始向量相同。
参看对应的特征向量发现，相对较大的数正是集中于中间位置，两端的分类近似等于0。
也就是说，对于第2种初始向量的选择越接近于特征向量所能得到的结果越准确，即能够保证该特征向量在参与产生初始向量时的系数a1越不可能为0。
而第1种初始向量，肯定能够保证结果的准确，但迭代数次会较第2种方法多很多。
实际操作中可以选择不同形式的初始向量2，利用求得的特征值，找到其中最大的即为需求的按模最大特征值。
升级会员 