北航数值分析大作业一.docx
- 文档编号:5405115
- 上传时间:2022-12-16
- 格式:DOCX
- 页数:14
- 大小:488.37KB
北航数值分析大作业一.docx
《北航数值分析大作业一.docx》由会员分享,可在线阅读,更多相关《北航数值分析大作业一.docx(14页珍藏版)》请在冰豆网上搜索。
北航数值分析大作业一
《数值分析A》
计算实习题目一
姓名:
学号:
学院:
2014年11月
一、算法设计方案:
1.A矩阵的存储与检索:
根据题目的设定,A矩阵为上半带宽与下半带宽均为2的带状矩阵,且除主对角线元素为变量ai外,其余元素均为常数b或c,故可以将A矩阵转存为新矩阵A[5][501],在存入时原A矩阵的主对角元素存入新A矩阵的第三行,且各元素的列号保持不变。
2.求解λ1、λ501、λs逻辑关系图:
3.求解A与数
最接近的特征值
:
可对A矩阵进行平移变换,对矩阵
采用反幂法求解按模最小的特征值
,再经过反平移得到
即为矩阵A与
最接近的特征值。
4.求矩阵A的(谱范数)条件数
和行列式
:
在采用反幂法求解
的过程中要对矩阵A进行Doolittle分解(LU分解),
为U矩阵对角线元素的乘积。
而根据公式
可以得到A的条件数,其中
和
分别为矩阵A按模最大特征值和按模最小特征值。
二、程序源代码
#include
#include
#include
#defineN501
#defineM5
#definee1.0e-12
doublea[M][N]={0.0},u[N]={0.0};
voidread_A()/*设定A矩阵函数*/
{
doubleb=0.16,c=-0.064;
inti;
doublea_0[N];
for(i=1;i<=N;i++)
a_0[i-1]=(1.64-0.024*i)*sin(0.2*i)-0.64*exp(0.1/i);
for(i=2;i a[0][i]=c; for(i=1;i a[1][i]=b; for(i=0;i a[2][i]=a_0[i]; for(i=0;i a[3][i]=b; for(i=0;i a[4][i]=c; } voidread_u()/*读取初始迭代向量*/ { inti; for(i=0;i u[i]=double(1231.0*rand()/1991.0); } double*unitization(doubleu[N],doubletotal)/*二范数向量单位化函数*/ {inti; doubley[N]; for(i=0;i y[i]=u[i]/total; return(y); } doublesecond_num(doubleu[N])/*求矩阵二范数函数*/ { inti; doublesum=0.0,total=0.0; for(i=0;i sum+=u[i]*u[i]; total=sqrt(sum); return(total); } doublemin(doublea,doubleb)/*最小值函数*/ { if(a returna; else returnb; } doublemax(doublea,doubleb)/*最大值函数*/ { if(a returnb; else returna; } double*Get_u(doublea[M][N],doubley[N])/*幂法中获得新迭代向量*/ { inti; doubleu[N]={0.0}; u[0]=a[2][0]*y[0]+a[1][1]*y[1]+a[0][2]*y[2]; u[1]=a[3][0]*y[0]+a[2][1]*y[1]+a[1][2]*y[2]+a[0][3]*y[3]; u[N-2]=a[4][N-4]*y[N-4]+a[3][N-3]*y[N-3]+a[2][N-2]*y[N-2]+a[1][N-1]*y[N-1]; u[N-1]=a[4][N-3]*y[N-3]+a[3][N-2]*y[N-2]+a[2][N-1]*y[N-1]; for(i=2;i u[i]=a[4][i-2]*y[i-2]+a[3][i-1]*y[i-1]+a[2][i]*y[i]+a[1][i+1]*y[i+1]+a[0][i+2]*y[i+2]; return(u); } voidResolve_LU(doublea[M][N])/*LU分解函数*/ { intk,i,j,t; doubletemp; for(k=1;k<=N;k++) { for(j=k;j<=min(k+2,N);j++) { temp=0; for(t=max(max(1,k-2),j-2);t<=k-1;t++) temp+=a[k-t+2][t-1]*a[t-j+2][j-1]; a[k-j+2][j-1]=a[k-j+2][j-1]-temp; } for(i=k+1;i<=min(k+2,N);i++) { temp=0; for(t=max(max(1,i-2),k-2);t<=k-1;t++) temp+=a[i-t+2][t-1]*a[t-k+2][k-1]; a[i-k+2][k-1]=(a[i-k+2][k-1]-temp)/a[2][k-1]; } } } double*back(doublea[M][N],doubleu[N])/*方程组回代函数*/ { inti,t; doubletemp=0.0; for(i=2;i<=N;i++) {temp=0.0; for(t=max(1,i-2);t temp+=a[i-t+2][t-1]*u[t-1]; u[i-1]-=temp; } u[N-1]=u[N-1]/a[2][N-1]; for(i=N-1;i>0;i--) { temp=0.0; for(t=i+1;t<=min(i+2,N);t++) temp+=a[i-t+2][t-1]*u[t-1]; u[i-1]=(u[i-1]-temp)/a[2][i-1]; } return(u); } doublefanmifa(doublea[M][N],doubleu[N])/*反幂法求特征值*/ {inti,j,k; doubley[N]={0.0},Y[N]={0.0},x,b,B,total,E; doubleea[M][N]={0.0}; doubleex[N]={0.0}; double*p,*q,*n; B=0.0; for(i=0;;i++) {b=1/B;B=0.0; total=second_num(u); p=unitization(u,total); for(j=0;j y[j]=*(p+j); for(j=0;j ex[j]=y[j]; for(j=0;j for(k=0;k ea[j][k]=a[j][k]; Resolve_LU(ea); q=back(ea,ex); for(j=0;j u[j]=*(q+j); total=second_num(u); n=unitization(u,total); for(j=0;j Y[j]=*(n+j); for(j=0;j B+=y[j]*u[j]; x=1/B; E=abs((1/B)-b)/abs(1/B); if(E<=e) break; } x=1/B; return(x); } doublemifa(doublea[M][N],doubleu[N])/*幂法求特征值*/ {inti,j; doubley[N]={0.0},Y[N]={0.0},x,b,B,total,E; double*p,*q,*n; B=0.0; for(i=0;;i++) {b=B;B=0.0; total=second_num(u); p=unitization(u,total); for(j=0;j y[j]=*(p+j); q=Get_u(a,y); for(j=0;j u[j]=*(q+j); total=second_num(u); n=unitization(u,total); for(j=0;j Y[j]=*(n+j); for(j=0;j B+=y[j]*u[j]; x=B; E=abs(B-b)/abs(B); if(E<=e) break; } x=B; return(x); } voiddisplace(doublea[M][N],doublex) { inti; for(i=0;i a[2][i]=a[2][i]+x; } voidmain()/*主程序*/ {/*第一部分求解三个特征值*/ inti; doubleEig_1,Eig_s,Eig_501; doubleEig_01,Eig_02,Eig_03; doubleA[M][N],U[39],Eig[39],Eig_temp[39]; doublecond_A=0.0,DetA=1.0; read_A(); read_u(); Eig_s=fanmifa(a,u);/*采用反幂法直接得到按模最小的特征值*/ read_u(); Eig_01=mifa(a,u);/*采用幂法获得按模最大的特征值作为平移量*/ displace(a,Eig_01); read_u(); Eig_02=fanmifa(a,u); Eig_03=Eig_02-Eig_01; Eig_1=min(Eig_01,Eig_03); Eig_501=max(Eig_01,Eig_03); printf("******************数值分析作业一**********************\n"); printf("结果为: \n"); printf(" (1)\n"); printf("特征值λ1=%19.11e\n",Eig_1); printf("特征值λ501=%19.11e\n",Eig_501); printf("特征值λs=%19.10e\n",Eig_s); /*第二部分求解指定的最接近特征值*/ printf(" (2)\n"); for(i=0;i<39;i++) { U[i]=Eig_1+((i+1)*(Eig_501-Eig_1)/40); read_A(); displace(a,-U[i]); read_u(); Eig_temp[i]=fanmifa(a,u); Eig[i]=Eig_temp[i]+U[i]; printf("k=%d,μ%d=%19.11e,最接近的特征值为: λ%d=%19.11e\n",i+1,i+1,U[i],i+1,Eig[i]); } /*第三部分求A的谱范数条件数cond(A)和行列式detA*/ printf("(3)\n"); read_A(); Resolve_LU(a); for(i=0;i DetA=DetA*a[2][i]; cond_A=abs(Eig_01/Eig_s); printf("cond(A)=%19.11e\n",cond_A); printf("DetA=%19.11e\n",DetA); system("pause"); } 三、程序运行结果: 结果截图如下 图1程序结果图 四、关于计算实习题目的讨论和思考: 1.迭代初始向量对结果的影响: 在求解矩阵特征值时使用的方法为幂法和反幂法,这两种方法都是迭代方法的一种形式,所以初始向量的选择对计算结果产生一定的影响。 当初始迭代向量中的0元素个数增大时(即 的无关数小于n),迭代收敛性趋向于不稳定,且收敛速度减慢。 在本程序中每一次使用幂法或者反幂法的初始迭代向量均为随机生成的。 为了进一步说明初始迭代向量的影响,这里单独提取出幂法迭代函数部分进行测试。 在这部分测试中,两次迭代循环的初始迭代向量为随机生成(完全不同),迭代结果如下图显示: 图2第一次迭代结果 从两个结果图的对比中可以发现,在精度要求(e-12)的要求下,因为结果显示即为12位有效数字,故虽两次迭代循环的初始迭代向量不同,但最终收敛的结果一致。 但在循环的过程中,可以发现第二次循环的初始向量中元素值为0的比例大于第一次循环,故第二次迭代的循环 图3第二次迭代结果 次数比第一次迭代多,耗时也较久。 接下来进一步讨论初始迭代向量u中各个分量的值对迭代的影响。 从迭代的结果和步数来看,当u[i]=u0(i=1,2,···,501)且u0的绝对值相对较大时,结果收敛但速度较慢,原因是u0与A矩阵元素数量级差距较大,导致了运算量增大;当u[i]=u0(i=1,2,···,501)且u0的绝对值相对较小时,结果收敛性不稳定,原因是u0较小会比较明显的受到计算舍入误差的影响,进而影响计算结果;当u[i]=u0(i=1,2,···,501)且各个u0之间数量级差距较大时,结果收敛性也不稳定,且收敛的速度减慢,原因是认为使迭代过程中的权重发生改变,使迭代更加复杂。 综合来看,迭代初始向量的选取最好应与矩阵A元素的数量级接近,且值为0的元素越少越好。 2.指数型格式控制“%m.ne”的设定: 根据题目的要求,输出格式应为有效数字超过12位的指数型格式,配合double浮点型数据的储存格式和“%m.ne”型指数输出格式的要求,在m位中有5位应分配给指数部分(E+000),并保证在小数点前可以保留符号,故整个输出位数设置为17,其中小数点后为11位数字,正好符合题目要求的输出格式。 3.循环语句使用过程中常见问题: 在编写程序过程中必须多次使用循环结构,有时还需要循环嵌套结构。 在C/C++语法中,for语句中循环变量的初值一般从0开始设置,这与实际矩阵或向量的行列数存在1的差值,且对循环终止条件的设定也需要相应的减小1。 若不注意这一点,很容易在程序循环过程中,数值的传递出现问题,最终导致结果错误。 在循环语句中,除了循环初值和终值问题以外,还有循环嵌套中一些用作存储累加值的变量初值设定。 若在某一次累加循环开始时未及时将累加值清零,则本次累加的终值就会出现错误,整个循环也就无法得到预计的值。 以上两点是在编写程序中出现的比较明显的错误。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 北航 数值 分析 作业