北航数值分析大作业一.docx
- 文档编号:23664982
- 上传时间:2023-05-19
- 格式:DOCX
- 页数:19
- 大小:217.35KB
北航数值分析大作业一.docx
《北航数值分析大作业一.docx》由会员分享,可在线阅读,更多相关《北航数值分析大作业一.docx(19页珍藏版)》请在冰豆网上搜索。
北航数值分析大作业一
《数值分析》大作业一
1、算法设计方案:
矩阵A的存储与检索:
将矩阵A[501][501],转存为新矩阵A[5][501]。
在存入时,按每一行依次存入,存入时原矩阵与新矩阵的列号保持不变,其中原矩阵的主对角线上的元素存为新矩阵的第三行,最后所有元素存入新矩阵中。
求λ1,λ501,λs
λs:
λs表示矩阵的按模最小特征值,为求得λs直接对待求矩阵A应用反幂法即可。
λ1、λ501:
已知矩阵A的特征值满足关系λ1<…<…λ501,要求λ1、及λ501时,可按如下方法求解:
a.对矩阵A用幂法,求得按模最大的特征值λn1。
b.按平移量λn1对矩阵A进行原点平移得矩阵,对矩阵B=A+λn1I用反幂法求得B的按模最小特征值λn2。
c.λn3=λn2一λn1
d.则:
λ1=min(λn1,λn3),λ501=max(λn1,λn3)
求A中与数μk=λ1+k(λ501-λ1)/40最接近的特征值λik(k=1,…39):
先求矩阵B=A-
I对应的按模最小特征值
,则
+
即为矩阵A与
最接近的特征值。
重复以上过程39次即可求得
(k=0,1,…39)的值。
求A的(谱范数)条件数
和行列式
:
在
(1)中用反幂法求矩阵A的按模最小特征值时,要用到Doolittle分解方法,在Doolittle分解完成后得到的两个矩阵分别为L和U,detA等于U所有对角线上元素的乘积。
,λmax和λmin分别为模最大特征值与模最小特征值。
2、程序源代码:
#include
#include
#include
#defineN501/*定义列数*/
#defineM5/*定义行数*/
#definee1.0e-12/*误差限*/
doubleA[M][N];/*初始矩阵*/
doubleu[N];/*初始向量*/
doubley[N],yy[N];
doublemaximum,TZ1,TZ2,TZ_1,TZ_N,TZ_s,TZ_abs_max;
intmax_sign,max_position;
/*对矩阵A进行初始化程序,存为A[5][501]*/
voidCSH_A()
{
doubleb=0.16,c=-0.064;
inti;
for(i=2;i A[0][i]=c; for(i=1;i A[1][i]=b; for(i=0;i A[2][i]=(1.64-0.024*(i+1))*sin(0.2*(i+1))-0.64*exp(0.1/(i+1)); for(i=0;i A[3][i]=b; for(i=0;i A[4][i]=c; } /*初始化迭代向量,且随机赋值*/ voidCSH_u() { inti; for(i=0;i u[i]=double(15*rand()/32767); } /*幂法中对迭代向量进行单位化程序*/ voidGet_y() { inti; for(i=0;i y[i]=u[i]/maximum; } /*幂法中求解迭代向量的无穷范数*/ voidGet_max() { inti; max_position=0; maximum=fabs(u[0]); for(i=1;i { if(maximum { max_position=i; maximum=fabs(u[i]); } } if(u[max_position]<0) max_sign=-1; else max_sign=1; } /*获得新迭代向量*/ voidGet_u() { inti; 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]; } /*获得迭代后特征值*/ voidGet_TZ() { TZ2=TZ1; TZ1=max_sign*u[max_position]; } /*求解迭代向量无穷范数的幂法*/ voidCheck_TZ() { inti; CSH_u(); Get_max(); Get_y(); Get_u(); Get_TZ(); for(i=0;;i++) { Get_max(); Get_y(); Get_u(); Get_TZ(); if(fabs((TZ2-TZ1)/TZ1) break; } } /*获取绝对值最大的特征值λ_501*/ voidThe_TZ() { Check_TZ(); TZ_abs_max=TZ1; } /*获取特征值λ_1*/ voidThe_Other_TZ() { inti; doubleTZ_temp=TZ1; for(i=0;i { A[2][i]-=TZ_temp; } Check_TZ(); TZ1+=TZ_temp; if(TZ1 { TZ_1=TZ1; TZ_N=TZ_temp; } else { TZ_N=TZ1; TZ_1=TZ_temp; } } /*两值中取最小*/ intmin(inta,intb) { if(a returna; else returnb; } /*两值中取最大*/ intmax(inta,intb) { if(a returnb; else returna; } /*把分解矩阵为LU*/ voidResolve_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]; } } } /*方程组回代过程*/ voidGO_Back_substitution() { inti,t; doubletemp=0; for(i=2;i { for(t=max(1,i-2);t y[i-1]-=A[i-t+2][t-1]*y[t-1]; } u[N-1]=y[N-1]/A[2][N-1]; for(i=N-1;i>0;i--) { temp=0; for(t=i+1;t<=min(i+2,N);t++) temp+=A[i-t+2][t-1]*u[t-1]; u[i-1]=(y[i-1]-temp)/A[2][i-1]; } } /*求矩阵行列式值*/ doubleDet_matrix() { inti; doubledet=1; CSH_A(); Resolve_LU(); for(i=0;i det=det*A[2][i]; returndet; } /*反幂法中获得迭代向量2一范数*/ doubleGet_norm() { inti; doublenormal=0; for(i=0;i normal+=u[i]*u[i]; normal=sqrt(normal); returnnormal; } /*反幂法中对迭代向量单位化*/ voidGet_yy(doublenormal) { inti; for(i=0;i { y[i]=u[i]/normal; yy[i]=y[i]; } } /*获得绝对值最小的特征值*/ voidGet_TZ_s() { inti; TZ2=TZ1; TZ1=0; for(i=0;i TZ1+=yy[i]*u[i]; TZ1=1/TZ1; } /*反幂法求绝对值最小的特征值*/ voidTZ_min() { doublenorm=0; intcount=0; TZ1=0,TZ2=0; CSH_u(); norm=Get_norm(); Get_yy(norm); GO_Back_substitution(); Get_TZ_s(); while(count<10000) { count++; norm=Get_norm(); Get_yy(norm); GO_Back_substitution(); Get_TZ_s(); if(fabs((TZ2-TZ1)/TZ1) break; } TZ_s=TZ1; } /*求矩阵条件数*/ doubleGet_cond_A()/*求矩阵条件数*/ { doublecond1; cond1=fabs(TZ_abs_max/TZ_s); returncond1; } voidTZ_trans_min()/*偏移条件下反幂法求特征值*/ { inti,k; doubletr,uu[40]; for(k=1;k<40;k++) { tr=TZ_1+k*(TZ_N-TZ_1)/40; CSH_A(); for(i=0;i A[2][i]-=tr; Resolve_LU(); TZ_min(); TZ_s+=tr; uu[k]=TZ_s+k*((TZ_N-TZ_1))/k; printf("与数μ%2d=%.12e最近的特征值λ%d=%.12e\n",k,uu[k],k,TZ_s); printf("\n"); } } /*主程序*/ voidmain() { doublecond; doubleTZ_det; CSH_A();/*初始化矩阵A*/ The_TZ();/*获取绝对值最大的特征值λ_501*/ The_Other_TZ();/*获取特征值λ_1*/ printf("输出结果\n"); printf("_____________________________________\n"); printf("输出特征值λ1=%.12e\n",TZ_1); printf("_____________________________________\n"); printf("\n"); printf("_____________________________________\n"); printf("输出特征值λ501=%.12e\n",TZ_N); printf("_____________________________________\n"); printf("\n"); TZ_det=Det_matrix();/*求矩阵行列式值*/ TZ_min();/*反幂法求绝对值最小的特征*/ printf("_____________________________________\n"); printf("输出特征值λs=%.12e\n",TZ_s); printf("_____________________________________\n"); printf("\n"); cond=Get_cond_A();/*求矩阵条件数*/ TZ_trans_min();/*偏移条件下反幂法求特征值*/ printf("条件数cond(A)=%.12e\n",cond); printf("\n"); printf("行列式det(A)=%.12e\n",TZ_det); } 3、程序运行结果: 第一次输出结果: 第二次输出结果: 4、迭代初始向量的选取对计算结果的影响 计算实习中求矩阵A的具有某些特征的特征值,主要用到的方法是幂法和反幂法,这两种方法从原理上看都是迭代法,因此迭代初始向量的选择对计算结果会产生一定影响,主要表现在收敛速度上。 初始迭代向量中元素等于0的个数越少,收敛结果越稳定;初始迭代向量中元素等于0的个数越多,收敛结果越不稳定; 迭代初始值u[i]=s(i=1,2,…,501)且s的绝对值值极大,收敛结果可以稳定但收敛速度减慢,其原因为s的数量级与矩阵A中元素数量级差距过大,导致迭代次数以及运算量增大; 迭代初始值u[i]=s(i=1,2,…,501)且s的绝对值值极小,收敛结果并不稳定,且收敛速度减慢,其原因是计算机舍入误差将会影响计算结果; 迭代初始值u[i]=s(i=1,2,…,501)之间数量级偏差很大,收敛结果亦不稳定,且收敛速度减慢,其原因是人为使迭代过程中的权重发生较大区别,使迭代复杂化 结论,对于迭代初始值的选取应尽量与矩阵A中元素数量级保持相近,应保证相近的数量级,且元素等于0的个数越少越好。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 北航 数值 分析 作业