北航研究生数值分析编程大作业1.docx
《北航研究生数值分析编程大作业1.docx》由会员分享，可在线阅读，更多相关《北航研究生数值分析编程大作业1.docx（20页珍藏版）》请在冰点文库上搜索。
北航研究生数值分析编程大作业1
数值分析大作业
1、算法设计方案
1、矩阵初始化
矩阵
的下半带宽r=2,上半带宽s=2，设置矩阵
，在矩阵C中检索矩阵A中的带内元素
的方法是：
。
这样所需要的存储单元数大大减少，从而极大提高了运算效率。
2、利用幂法求出
幂法迭代格式：
当
时，迭代终止。
首先对于矩阵A利用幂法迭代求出一个
，然后求出矩阵B，其中
（
为单位矩阵），对矩阵B进行幂法迭代，求出
，之后令
，比较
，大者为
，小者为
。
3、利用反幂法求出
反幂法迭代格式：
当
时，迭代终止，
。
每迭代一次都要求解一次线性方程组
，求解过程为：
（1）作分解
对于
执行
（2）求解
（数组b先是存放原方程组右端向量，后来存放中间向量y）
使用反幂法，直接可以求得矩阵按模最小的特征值
。
求与数
最接近的特征值
，对矩阵
实行反幂法，即可求出对应的
。
4、求出A的条件数和行列式
根据
，其中分子分母分别对应按模最大和最小的特征值。
的计算：
由于
其中
为下三角矩阵，且对角线元素为1，故
，所以有
，又
为上三角矩阵，故
为对其对角线上各元素的乘积，最后可得
。
2、程序源代码
（1）定义所需要的函数：
#include
#include
#include
#defineN501
#defineR2
#defineS2
intmin（inta,intb）;//求最小值
intmax（inta,intb,intc）;//求最大值
doubleFan_two（doublex[N]）;//计算二范数
voidFenjieLU（double（*C）[N]）;//解线性方程组的LU分解过程
voidSolve（double（*C）[N],double*b,double*x）;//解线性方程组的求解过程
doublePowerMethod（doubleC[][N],doubleu[N],doubley[N],doublebta,doubleD）;//幂法
doubleInversePowerMethod（doubleC[][N],doubleu[N],doubley[N],doublebta,doubleD）;//反幂法
};
（2）程序的主函数，Main.cpp代码如下：
voidmain（）
{
doubleC[R+S+1][N];
doubleu[N];
doubley[N];
doublemiu[39];
doubleC1[R+S+1][N];
doublebta=1.0;
doubleNamda1,Namda501,NamdaS;
doubleNamda[39];
doubleCondA2;
doubledetA=1.0;
doubleD=1.0e-12;
inti,j,k;
FILE*fp;
fp=fopen（"Namda.txt","w"）;
//对数组进行初始化//
inti,j;
for（i=0;i{
u[i]=1;
}
for（i=0;i{
for（j=0;j{
if（i==0||i==4）
{
C[i][j]=-0.064;
}
elseif（i==1||i==3）
{
C[i][j]=0.16;
}
elseif（i==2）
{
C[i][j]=（1.64-0.024*（j+1））*sin（0.2*（j+1））
-0.64*exp（0.1/（j+1））;
}
}
}
//幂法求Namda1//
Namda1=PowerMethod（C,u,y,bta,D）;
printf（"\n================================================\n"）;
printf（"Namda1=%12.11e",Namda1）;
printf（"\n================================================\n"）;
//幂法求Namda501//
bta=1.0;
for（i=0;i{
for（j=0;j{
if（i==2）
C1[i][j]=C[i][j]-Namda1;
else
C1[i][j]=C[i][j];
}
}
Namda501=algorism.PowerMethod（C1,u,y,bta,D）+Namda1;
printf（"\n================================================\n"）;
printf（"Namda501=%12.11e",Namda501）;
printf（"\n================================================\n"）;
//反幂法求NamdaS//
bta=1.0;
NamdaS=InversePowerMethod（C,u,y,bta,D）;
printf（"\n================================================\n"）;
printf（"NamdaS=%12.11e",NamdaS）;
printf（"\n================================================\n"）;
//反幂法求Namda[k]//
printf（"\n================================================\n"）;
for（k=0;k<39;k++）
{
miu[k]=Namda1+（k+1）*（Namda501-Namda1）/40.0;
bta=1.0;
for（i=0;i{
for（j=0;j{
if（i==2）
C1[i][j]=C[i][j]-miu[k];
else
C1[i][j]=C[i][j];
}
}
Namda[k]=InversePowerMethod（C1,u,y,bta,D）+miu[k];
fprintf（fp,"与%12.11e最接近的特征值为:
%12.11e\n",miu[k],Namda[k]）;
}
printf（"求与miu[k]最接近的Namda[k]的计算结果已经输出到文件Namda.txt中"）;
printf（"\n================================================\n"）;
//求A的谱范数//
printf（"\n================================================\n"）;
printf（"A的谱范数为:
%12.11e",sqrt（Namda501））;
printf（"\n================================================\n"）;
//求A的条件数//
CondA2=fabs（Namda1/NamdaS）;
printf（"\n================================================\n"）;
printf（"A的谱范数的条件数Cond（A）2为:
%12.11e",CondA2）;
printf（"\n================================================\n"）;
//求det（A）2的值//
for（j=0;jdetA*=C[2][j];
printf（"\n================================================\n"）;
printf（"行列式A的值为:
%12.11e",detA）;
printf（"\n================================================\n"）;
fclose（fp）;
_getch（）;
return;
}
（3）成员函数的实现
intmin（inta,intb）
{
returnaa:
b;
}
intmax（inta,intb,intc）
{
inttemp;
temp=a>b?
a:
b;
returntemp>c?
temp:
c;
}
doubleFan_two（doublex[N]）
{
doublesum=0.0;
inti;
for（i=0;i{
sum+=pow（x[i],2）;
}
returnsqrt（sum）;
}
voidFenjieLU（double（*C）[N]）
{
doublesum=0;
inti,j,k,t;
for（k=0;k{
j=k;
i=k+1;
while
（1）
{
if（j==min（k+S+1,N））
break;
for（t=max（0,k-R,j-S）;t<=k-1;t++）
{
sum+=C[k-t+S][t]*C[t-j+S][j];
}
C[k-j+S][j]=C[k-j+S][j]-sum;
sum=0.0;
j++;
if（k==N-1）
break;
if（i==min（k+R+1,N））
break;
for（t=max（0,i-R,k-S）;t<=k-1;t++）
{
sum+=C[i-t+S][t]*C[t-k+S][k];
}
C[i-k+S][k]=（C[i-k+S][k]-sum）/C[S][k];
sum=0;
i++;
}
}
}
voidSolve（double（*C）[N],double*b,double*x）
{
doublesum=0;
inti,t;
sum=0;
for（i=1;i{
for（t=max（0,i-R）;t<=i-1;t++）
{
sum+=C[i-t+S][t]*b[t];
}
b[i]=b[i]-sum;
sum=0;
}
x[N-1]=b[N-1]/C[S][N-1];
for（i=N-2;i>=0;i--）
{
for（t=i+1;t<=min（i+S,N-1）;t++）
{
sum+=C[i-t+S][t]*x[t];
}
x[i]=（b[i]-sum）/C[S][i];
sum=0;
}
}
doublePowerMethod（doubleC[][N],doubleu[N],doubley[N],doublebta,doubleD）
{
doubleita;
doublesum=0;
doubletemp=0.0;
inti,j,k=0;
while（fabs（bta-temp）/fabs（bta）>D）
{
temp=bta;
ita=Fan_two（u）;
for（i=0;i{
y[i]=u[i]/ita;
}
for（i=0;i{
for（j=max（0,i-R）;j{
sum+=C[i-j+S][j]*y[j];
}
u[i]=sum;
sum=0;
}
for（i=0;i{
sum+=y[i]*u[i];
}
bta=sum;
sum=0;
k++;
}
returnbta;
}
doubleInversePowerMethod（doubleC[][N],doubleu[N],doubley[N],doublebta,doubleD）
{
doubleTC[R+S+1][N];
doublety[N];
doubleita;
doublesum=0;
doubletemp=0.0;
inti,j,k=0;
FenjieLU（C）;
while（abs（1/bta-1/temp）/abs（1/bta）>D）
{
temp=bta;
ita=Fan_two（u）;
for（i=0;i{
y[i]=u[i]/ita;
}
//用到临时存储数组TC[][]和ty[][]是因为函数Solve执行过程中会改变A[][]和y[][]
for（i=0;i{
for（j=0;jTC[i][j]=C[i][j];
}
for（i=0;ity[i]=y[i];
Solve（C,y,u）;
for（i=0;i{
for（j=0;jC[i][j]=TC[i][j];
}
for（i=0;iy[i]=ty[i];
for（i=0;i{
sum+=y[i]*u[i];
}
bta=sum;
sum=0;
k++;
}
bta=1.0/bta;
returnbta;
}
3、程序运行结果
下图为主程序运行结果
其中
的结果输出在Namda.txt文件中，结果如下：
四、分析迭代初始向量对计算结果的影响
选择不同的初始向量
可能会得到不同的特征值。
选取
时，运行结果如下：
选取
时，运行结果如下：
选取
时（i=int（N/2）时为0），运行结果如下：
选取
时（i=int（N/2）时为1），运行结果如下：
通过以上类似的实验可以大致看出这样的规律：
的值趋近于
有两种情况：
（1）当
的元素中，1的个数较多时;
（2）在1的个数相同的条件下，1的分布越靠中后段，
观察
对应的特征向量可以发现：
（1）随着i的增加，特征向量元素的绝对值不断增大，即绝对值较大的数集中于中后位置。
因此，如果初始向量的非零元素集中在中后段，该初始向量会更容易逼近对应的特征向量，得到的结果也越准确。
对于，初始向量的非零元素集中在前半段的情况进行实验，会发现当算法中不考虑给定的精度水平，强制性执行足够高次数（大约在300多次以上）的迭代，运算结果也会趋近于
。
这就说明，程序之前没有得到准确结果的原因，是因为迭代次数不够。
当迭代次数在100到200次左右时，每一次迭代所造成的相对误差小于给定的精度水平，因此，如果由精度水平来控制循环迭代的次数，程序将错误地判断已经收敛，但实际上，当继续迭代到300次以上时，运算结果会突然变化，直至最终稳定在
。
由此，可以得出结论，当迭代次数足够高（300次以上）时，得到的结果会趋于稳定，不同的初始向量和选定的精度水平，决定着程序是否出现以及何时出现假收敛。
当所选取初始向量的非零元素越多，以及非零元素的位置越靠后时，收敛会更加迅速、准确。
升级会员 