北航数值分析大作业一Word格式文档下载.doc
《北航数值分析大作业一Word格式文档下载.doc》由会员分享，可在线阅读，更多相关《北航数值分析大作业一Word格式文档下载.doc（10页珍藏版）》请在冰点文库上搜索。
使用位移的方式求得另一特征值即可。
2求A的与数最接近的特征值。
题目可看成求以为偏移量后，按模最小的特征值。
即以为偏移量做位移，使用反幂法求出按模最小特征值后，加上，即为所求。
3求A的（谱范数）条件数和行列式detA。
矩阵A为非奇异对称矩阵，可知，
（1-1）
其中为按模最大特征值，为按模最小特征值。
detA可由LU分解得到。
因LU均为三角阵，那么其主对角线乘积即为A的行列式。
二算法实现
1幂法
使用如下迭代格式：
（2-1）
终止迭代的控制理论使用，
实际使用
（2-2）
由于不保存A矩阵中的零元素，只保存主对角元素a[501]及b,c值。
那么上式中简化为：
（2-3）
2反幂法
（2-4）
其中，解方程求出。
求解过程中使用LU分解，由于A为5对角矩阵，选择追赶法求取LU分解。
求解过程如下：
追赶法求LU分解的实现：
（2-5）
由上式推出分解公式如下：
（2-6）
推导出回代求解公式如下：
（2-7）
（2-8）
3及A行列式求解
（2-9）
由式（2-5）可得：
（2-10）
三源程序
#include<
stdio.h>
math.h>
doubleep=1e-12,b=0.16,c=-0.064;
intj=0;
doublepower（doublea[501]）;
//幂法
doubleinv_power（doublea[501]）;
//反幂法
doubledet（doublea[501]）;
//求det
intmain（） //主程序
{
inti,k;
doubleA[501],B[501],beta_1,beta_501,beta_s,beta_k;
doublemu;
for（i=0;
i<
501;
i++）
A[i]=（1.64-0.024*（i+1））*sin（0.2*（i+1））-0.64*exp（0.1/（i+1））;
beta_1=power（A）;
//第一问
printf（"
λ1\t=%.12e\t迭代次数：
%d\n"
beta_1,j）;
i++）//位移
B[i]=A[i]-beta_1;
beta_501=power（B）+beta_1;
printf（"
λ501\t=%.12e\t迭代次数：
beta_501,j）;
beta_s=inv_power（A）;
λs\t=%.12e\t迭代次数：
beta_s,j）;
for（k=1;
k<
=39;
k++）//第二问
{
mu=beta_1+k*（beta_501-beta_1）/40;
for（i=0;
B[i]=A[i]-mu;
beta_k=inv_power（B）+mu;
printf（"
λi%d\t=%.12e\t迭代次数：
k,beta_k,j）;
}
cond（A）2=%.12e\n"
beta_1/beta_s）;
//第三问
detA\t=%.12e\n"
det（A））;
}
doublepower（doublea[501]）//幂法
inti=0,N=5000;
doubleb=0.16,c=-0.064;
doubleu[501],y[501];
doublem=1,beta;
u[i]=1;
j=0;
while（j<
N）
{
for（i=0;
{
y[i]=u[i]/fabs（m）;
}
u[0]=a[0]*y[0]+b*y[1]+c*y[2];
u[1]=b*y[0]+a[1]*y[1]+b*y[2]+c*y[3];
u[499]=c*y[497]+b*y[498]+a[499]*y[499]+b*y[500];
u[500]=c*y[498]+b*y[499]+a[500]*y[500];
for（i=2;
499;
{u[i]=c*y[i-2]+b*y[i-1]+a[i]*y[i]+b*y[i+1]+c*y[i+2];
}
beta=0;
for（i=0;
{
if（fabs（u[i]）>
=fabs（beta））
beta=u[i];
}
if（beta<
0）
if（fabs（fabs（beta）-fabs（m））/fabs（beta）<
ep）
break;
if（fabs（beta-m）/fabs（beta）<
break;
m=beta;
j++;
}
returnbeta;
doubleinv_power（doublea[501]）//反幂法
doublep[501],r[501],t[501],q[501],u[501],y[501];
doublebeta,m=1;
inti,N=1000;
p[0]=a[0];
t[0]=b/p[0];
r[1]=b;
p[1]=a[1]-r[1]*t[0];
q[0]=c/p[0];
q[1]=c/p[1];
t[1]=（b-r[1]*q[0]）/p[1];
for（i=2;
r[i]=b-c*t[i-2];
p[i]=a[i]-c*q[i-2]-r[i]*t[i-1];
q[i]=c/p[i];
t[i]=（b-r[i]*q[i-1]）/p[i];
}
u[i]=1;
j=0;
y[i]=u[i]/fabs（m）;
u[0]=y[0]/p[0];
u[1]=（y[1]-r[1]*u[0]）/p[1];
for（i=2;
u[i]=（y[i]-c*u[i-2]-r[i]*u[i-1]）/p[i];
u[499]=u[499]-t[499]*u[500];
for（i=498;
i>
=0;
i--）
u[i]=u[i]-t[i]*u[i+1]-q[i]*u[i+2];
beta=0;
if（beta<
return1/beta;
doubledet（doublea[501]）//求det
doubledet_A=1;
doublep[501],r[501],t[501],q[501];
inti;
}
det_A=det_A*p[i];
returndet_A;
四程序结果
五计算过程中的现象
使用作为终止迭代条件时，出现迭代无法终止的情况，通过调试发现按模最大特征值为负时，当k充分大后，迭代向量各分量不断变号，使得与异号，判别式不收敛。
因此将终止迭代条件修改为，程序实现如下：
if（beta<
if（fabs（beta-m）/fabs（beta）<
break;
从迭代次数可以看出与收敛较慢，由按模最大特征值与按模次大特征值的比值越小，收敛速度越慢，可知存在与和的模相近的特征值。
【本文档内容可以自由复制内容或自由编辑修改内容期待你的好评和关注，我们将会做得更好】
页脚下载后可删除，如有侵权请告知删除！
升级会员 