北航研究生数值分析编程大作业1.docx
- 文档编号:23986550
- 上传时间:2023-05-23
- 格式:DOCX
- 页数:20
- 大小:343.60KB
北航研究生数值分析编程大作业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;j detA*=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;j TC[i][j]=C[i][j]; } for(i=0;i ty[i]=y[i]; Solve(C,y,u); for(i=0;i { for(j=0;j C[i][j]=TC[i][j]; } for(i=0;i y[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 选取 时(i 通过以上类似的实验可以大致看出这样的规律: 的值趋近于 有两种情况: (1)当 的元素中,1的个数较多时; (2)在1的个数相同的条件下,1的分布越靠中后段, 观察 对应的特征向量可以发现: (1)随着i的增加,特征向量元素的绝对值不断增大,即绝对值较大的数集中于中后位置。 因此,如果初始向量的非零元素集中在中后段,该初始向量会更容易逼近对应的特征向量,得到的结果也越准确。 对于,初始向量的非零元素集中在前半段的情况进行实验,会发现当算法中不考虑给定的精度水平,强制性执行足够高次数(大约在300多次以上)的迭代,运算结果也会趋近于 。 这就说明,程序之前没有得到准确结果的原因,是因为迭代次数不够。 当迭代次数在100到200次左右时,每一次迭代所造成的相对误差小于给定的精度水平,因此,如果由精度水平来控制循环迭代的次数,程序将错误地判断已经收敛,但实际上,当继续迭代到300次以上时,运算结果会突然变化,直至最终稳定在 。 由此,可以得出结论,当迭代次数足够高(300次以上)时,得到的结果会趋于稳定,不同的初始向量和选定的精度水平,决定着程序是否出现以及何时出现假收敛。 当所选取初始向量的非零元素越多,以及非零元素的位置越靠后时,收敛会更加迅速、准确。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 北航 研究生 数值 分析 编程 作业