数值分析之幂法及反幂法C语言程序实例.docx
- 文档编号:3151656
- 上传时间:2022-11-18
- 格式:DOCX
- 页数:14
- 大小:155.71KB
数值分析之幂法及反幂法C语言程序实例.docx
《数值分析之幂法及反幂法C语言程序实例.docx》由会员分享,可在线阅读,更多相关《数值分析之幂法及反幂法C语言程序实例.docx(14页珍藏版)》请在冰豆网上搜索。
数值分析之幂法及反幂法C语言程序实例
数值分析之幂法及反幂法C语言程序实例
1、算法设计方案:
①求
、
和
的值:
:
表示矩阵的按模最小特征值,为求得
直接对待求矩阵A应用反幂法即可。
、
:
已知矩阵A的特征值满足关系
要求
、及
时,可按如下方法求解:
a.对矩阵A用幂法,求得按模最大的特征值
。
b.按平移量
对矩阵A进行原点平移得矩阵
,对矩阵B用反幂法求得B的按模最小特征值
.
c.
则:
,
即为所求。
②求和A的与数
最接近的特征值
(k=0,1,…39):
求矩阵A的特征值中与
最接近的特征值的大小,采用原点平移的方法:
先求矩阵B=A-
I对应的按模最小特征值
则
+
即为矩阵A与
最接近的特征值。
重复以上过程39次即可求得
(k=0,1,…39)的值。
③求A的(谱范数)条件数
和行列式
:
在
(1)中用反幂法求矩阵A的按模最小特征值时,要用到Doolittle分解方法,在Doolittle分解完成后得到的两个矩阵分别为L和U,则A的行列式可由U阵求出,即:
det(A)=det(U)。
求得det(A)不为0,因此A为非奇异的实对称矩阵,则:
,
和
分别为模最大特征值与模最小特征值。
2、程序源代码:
#include #include〈stdio。 h〉 #include〈math.h〉 #defineN501//列 #defineM5//行 #defineR2//下带宽 #defineS2//上带宽 #defineK39 #definee1.0e—12//误差限 floatA[M][N];//初始矩阵 floatu[N];//初始向量 floaty[N],yy[N]; floatmaximum,value1,value2,value_1,value_N,value_s,value_abs_max; constfloatb=0.16f,c=-0。 064f; intmax_sign,max_position; voidInit_matrix_A()//初始化矩阵A { inti; for(i=2;i { A[0][i]=c; } for(i=1;i〈N;i++) { A[1][i]=b; } for(i=0;i〈N;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〈N—2;i++) { A[4][i]=c; } } voidInit_u()//初始化迭代向量 { inti; for(i=0;i〈N;i++) u[i]=1.0; } voidGet_max()//获得绝对值最大的数值的模 { inti; max_position=0; maximum=fabs(u[0]); for(i=1;i { if(maximum〈fabs(u[i])) { max_position=i; maximum=fabs(u[i]); } } if(u[max_position]<0) max_sign=—1; elsemax_sign=1; } voidGet_y()//单位化迭代向量 { inti; for(i=0;i y[i]=u[i]/maximum; } 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〈N—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_value()//获得迭代后特征值 { value2=value1; value1=max_sign*u[max_position]; } voidCheck_value()//幂法第二迭代格迭代 { Init_u(); Get_max(); Get_y(); Get_u(); Get_value(); while (1) { Get_max(); Get_y(); Get_u(); Get_value(); if(fabs((value2-value1)/value1) break; } } voidThe_value()//获取绝对值最大的特征值λ_501 { Check_value(); value_abs_max=value1; } voidThe_Other_value()//获取特征值λ_1 { inti; floatvalue_temp=value1; for(i=0;i〈N;i++) { A[2][i]-=value_temp; } Check_value(); value1+=value_temp; if(value1〈value_temp) { value_1=value1; value_N=value_temp; } else { value_N=value1; value_1=value_temp; } } intmin(inta,intb)//两值中取最小 { if(a〈b) returna; else returnb; } intmax(inta,intb)//两值中取最大 { if(a returnb; else returna; } voidResolve_LU() { intk,i,j,t; floattemp; for(k=1;k〈=N;k++) { for(j=k;j<=min(k+S,N);j++) { temp=0; for(t=max(max(1,k—R),j-S);t<=k-1;t++) temp+=A[k—t+S][t-1]*A[t—j+S][j-1]; A[k—j+S][j-1]=A[k—j+S][j-1]—temp; } for(i=k+1;i〈=min(k+R,N);i++) { temp=0; for(t=max(max(1,i-R),k—S);t<=k—1;t++) temp+=A[i—t+S][t—1]*A[t-k+S][k—1]; A[i—k+S][k-1]=(A[i—k+S][k-1]—temp)/A[S][k-1]; } } } voidBack_substitution()//方程组回代过程 { inti,t; floattemp=0; for(i=2;i〈N+1;i++) { for(t=max(1,i—R);t y[i—1]-=A[i—t+S][t—1]*y[t—1]; } u[N—1]=y[N-1]/A[S][N-1]; for(i=N—1;i>0;i—-) { temp=0; for(t=i+1;t<=min(i+S,N);t++) temp+=A[i-t+S][t—1]*u[t-1]; u[i-1]=(y[i—1]—temp)/A[S][i-1]; } } doubleDet_matrix()//求矩阵行列式值 { inti; doubledet=1; Init_matrix_A(); Resolve_LU(); for(i=0;i det=det*A[2][i]; returndet; } floatGet_norm()//获得迭代向量模 { inti; floatnormal=0; for(i=0;i〈N;i++) normal+=u[i]*u[i]; normal=sqrt(normal); returnnormal; } voidGet_yy(floatnormal)//迭代向量单位化 { inti; for(i=0;i〈N;i++) { y[i]=u[i]/normal; yy[i]=y[i]; } } voidGet_value_s()//获得绝对值最小的特征值 { inti; value2=value1; value1=0; for(i=0;i value1+=yy[i]*u[i]; value1=1/value1; } voidValue_min()//反幂法求绝对值最小的特征值 { floatnorm=0; intcount=0; value1=0,value2=0; Init_u(); norm=Get_norm(); Get_yy(norm); Back_substitution(); Get_value_s(); while(count〈10000) { count++; norm=Get_norm(); Get_yy(norm); Back_substitution(); Get_value_s(); if(fabs((value2-value1)/value1) break; } value_s=value1; } floatGet_cond_A()//求矩阵条件数 { floatcond1; cond1=fabs(value_abs_max/value_s); returncond1; } voidValue_translation_min()//偏移条件下反幂法求特征值 { inti,k; floattr; for(k=1;k { tr=value_1+k*(value_N-value_1)/40; Init_matrix_A(); for(i=0;i〈N;i++) A[2][i]-=tr; Resolve_LU(); Value_min(); value_s+=tr; printf(”k=%d=>>〉λi%d=%.13e\n",k,k,value_s); } } voidmain()
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 分析 反幂法 语言 程序 实例