北航研究生数值分析编程大作业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)
{
returna
a:
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次以上)时,得到的结果会趋于稳定,不同的初始向量和选定的精度水平,决定着程序是否出现以及何时出现假收敛。
当所选取初始向量的非零元素越多,以及非零元素的位置越靠后时,收敛会更加迅速、准确。