数值分析实习第二题.docx
- 文档编号:11251691
- 上传时间:2023-02-26
- 格式:DOCX
- 页数:31
- 大小:697.69KB
数值分析实习第二题.docx
《数值分析实习第二题.docx》由会员分享,可在线阅读,更多相关《数值分析实习第二题.docx(31页珍藏版)》请在冰豆网上搜索。
数值分析实习第二题
数值分析实习第二题
一、算法设计方案
1.采用C语言进行算法设计输出。
2.由于矩阵A的为10×10的矩阵,所以采用双步位移的QR分解法求矩阵A的全部特征值。
3.先对矩阵A进行初始化赋值。
4.先采用Householder矩阵将矩阵A拟上三角化,得到矩阵A(n-1)。
5.为了完成题目要求,减少运算次数,此时运行独立函数QR_UP_HE()另外声明3个10×10局部的二维数组变量Q,R,RQ对A(n-1)进行一次QR分解,得到相对应的矩阵Q,R和乘积矩阵RQ,用于输出。
6.之后开始用双步位移法求解矩阵A的特征值。
在求解过程中,遇到解一元二次方程组,采用数组S_Re[2]来存根的实数部分,用数组S_Im[2]来存根的虚数部分,用公式
来求解方程
。
QR分解时候选用精度ε=10-12,不限制最大迭代次数。
7.将双步位移的判断逻辑进行合理设定如下:
1)变量m从9(采用10维数组)开始循环运算,当m小于零的时候即认为QR分解完成,结束循环。
2)先判断m是否为零,若为零,那么A[0][0]即为最后一个特征值,结束循环。
否则继续。
3)之后判断A[m][m-1]的绝对值是否小于精度要求,若小于,则认为A[m][m]即为其中一个特征值,运算维数降为m–1,跳过这次循环。
否则继续。
4)此时求二阶子阵
的两个特征值s1和s2。
5)此时如果m=1,那么这就是A剩下的最后两个特征值,结束循环。
否则继续。
6)如果A[m-1][m-2]的绝对值小于精度要求,那么即可判断s1和s2是A的其中两个特征值,运算维数降为m–2,跳过这次循环。
否则继续。
7)上述判断完之后,即可开始用双步位移的QR分解法循环求A得全部特征值了。
8.全部特征值求得后,将其打印出来,之后开始求A的相应于实特征值的特征向量。
特征向量采用解方程组(λI-A)X=0的方法进行,先用Gauss列主元的方法,将A上三角化。
由于b为0,所以不用记录A对角化过程中的行交换的次序。
9.在反解X的过程中,需先判断A的行是否全为0(即小于一个精度值,程序中精度值定为ε=10-9,取太小会导致程序判断出错,从而得出错误的特征值)。
若最后一行不为零,则对特征向量X的最后一个值赋1,即X=[0,0,0,…,0,1],若倒数第二行也不为零,则对特征向量X的倒数第二行的值也赋1,以此类推。
10.之后回代X,得到其中的一个特征向量X。
然后将其归一化,即可输出。
二、运算结果
1.矩阵经过拟上三角化后得到的矩阵
2.进行QR分解后所得到的矩阵Q、R和乘积矩阵RQ
3.矩阵的全部特征值
4.
相应于特征值的特征向量
三、源程序
#include
#include
#include
#defineEPS1.0e-12
#defineN10
doubleA[N][N];
doubleLambda_Re[N];
doubleLambda_Im[N];
voidINI_A(void)
//矩阵A进行赋值
{
chari,j;
for(i=0;i { for(j=0;j { if(j==i) { A[i][j]=1.5*cos(i+1+1.2*(j+1)); } else { A[i][j]=sin(0.5*(i+1)+0.2*(j+1)); } } } } voidVECT_MULTI(charMax_N,doubleVector[],doubleVect_Ret[],doublek) //向量数乘 { charn; for(n=0;n { Vect_Ret[n]=Vector[n]*k; } } doubleVECT_NORM(charRow,charStr,charEnd,doubleMat[][N]) //求矩阵A中列向量的模 { charn; doublenorm; for(n=Str,norm=0;n { norm+=Mat[n][Row]*Mat[n][Row]; } return(sqrt(norm)); } doubleVECT_NORM2(charMax_N,doubleVector_2[]) //向量2范数 { intn; doublenorm_2; for(n=0,norm_2=0;n { norm_2+=(Vector_2[n]*Vector_2[n]); } return(sqrt(norm_2)); } charSP_SGN(doubles) //取符号 { if(s>0) { return (1); } else { return(-1); } } doubleVECT_T_VEC(charMax_N,doubleVector_T[],doubleVector[]) //向量转置乘向量 { charn; doublepower_norm; for(n=0,power_norm=0;n { power_norm+=Vector_T[n]*Vector[n]; } return(power_norm); } voidVECT_PLUS(charMax_N,doubleVector_1[],doubleVector_2[],doubleVector_Ret[]) //向量加法 { chari; for(i=0;i { Vector_Ret[i]=Vector_1[i]+Vector_2[i]; } } voidMAT_VEC(charMax_N,doubleMatrix[][N],doubleVect[],doubleVect_Ret[]) //矩阵乘向量 { chari,j; doublecache[N]={0}; for(i=0;i { for(j=0;j { cache[i]+=Matrix[i][j]*Vect[j]; } } for(i=0;i { Vect_Ret[i]=cache[i]; } } voidMAT_T_VEC(charMax_N,doubleMatrix[][N],doubleVect[],doubleVect_Ret[]) //矩阵转置乘向量 { chari,j; doublecache[N]={0}; for(i=0;i { for(j=0;j { cache[i]+=Matrix[j][i]*Vect[j]; } } for(i=0;i { Vect_Ret[i]=cache[i]; } } voidVECT_EVA(charMax_N,doubleVect_new[],doubleVect_bef[]) //向量赋值 { chari; for(i=0;i { Vect_new[i]=Vect_bef[i]; } } voidMAT_MULTI(charMax_N,doubleMat_A[][N],doubleMat_B[][N],doubleRet_Mat[][N]) //矩阵乘法 { chari,j,k; for(i=0;i { for(j=0;j { Ret_Mat[i][j]=0; for(k=0;k { Ret_Mat[i][j]+=Mat_A[i][k]*Mat_B[k][j]; } } } } voidSOLVE_EQU(doubleCoe[2],doubleSol_Re[2],doubleSol_Im[2]) //解一元二次方程 { doubleDelta; Delta=Coe[0]*Coe[0]/4-Coe[1]; if(Delta>=0) { Sol_Re[0]=(-1)*Coe[0]/2+sqrt(Delta); Sol_Re[1]=(-1)*Coe[0]/2-sqrt(Delta); Sol_Im[0]=0; Sol_Im[1]=0; } else { Sol_Re[0]=(-1)*Coe[0]/2; Sol_Re[1]=(-1)*Coe[0]/2; Sol_Im[0]=sqrt((-1)*Delta); Sol_Im[1]=(-1)*sqrt((-1)*Delta); } } voidEVA_MAT(charMax_N,doubleBef[][N],doubleRet[][N]) //矩阵相等赋值 { chari,j; for(i=0;i { for(j=0;j { Ret[i][j]=Bef[i][j]; } } } voidEXCHANGE(double*No_1,double*No_2) //数据交换 { doubletemp; temp=*No_1; *No_1=*No_2; *No_2=temp; } voidUPPER_HESSENBERG(void) //矩阵拟上三角化 { charr,i,k,j; doubled,c,h; doublep[N],q[N],w[N],u[N]; for(r=0;r { for(i=r+2;i { if(A[i][r]! =0) { d=VECT_NORM(r,r+1,N,A); c=(-1)*SP_SGN(A[r+1][r])*d; h=c*c-c*A[r+1][r]; for(k=0;k { if(k>r+1) { u[k]=A[k][r]; } elseif(k==r+1) { u[k]=A[k][r]-c; } else { u[k]=0; } } VECT_MULTI(N,u,w,1/h); MAT_T_VEC(N,A,w,p); MAT_VEC(N,A,w,q); VECT_MULTI(N,u,w,-VECT_T_VEC(N,p,w)); VECT_PLUS(N,q,w,w); for(k=0;k { for(j=0;j { A[k][j]-=(w[k]*u[j]+u[k]*p[j]); } } break; } } } } voidQR_DISSOLUTION(charMax_N,doubleMat[][N]) //矩阵双步位移QR分解 { chari,j,r,k; doubled,c,h; doubleu[Max_N],w[Max_N],p[Max_N],q[Max_N]; for(r=0;r { for(i=r+1;i { if(Mat[i][r]! =0) { d=VECT_NORM(r,r,Max_N,Mat); c=(-1)*SP_SGN(Mat[r][r])*d; h=c*c-c*Mat[r][r]; for(k=0;k { if(k>r) { u[k]=Mat[k][r]; } elseif(k==r) { u[k]=Mat[k][r]-c; } else { u[k]=0; } } VECT_MULTI(Max_N,u,w,1/h); MAT_T_VEC(Max_N,Mat,w,p); for(k=0;k { for(j=0;j { Mat[k][j]-=(u[k]*p[j]); } } MAT_T_VEC(Max_N,A,w,p); MAT_VEC(Max_N,A,w,q); VECT_MULTI(Max_N,u,w,-VECT_T_VEC(Max_N,p,w)); VECT_PLUS(Max_N,q,w,w); for(k=0;k { for(j=0;j { A[k][j]-=(w[k]*u[j]+u[k]*p[j]); } } break; } } } } voidDOUBLE_DISPLACE(void) //双步位移法求特征值 { chari,j,m; doubleS_Re[2],S_Im[2],D_Coe[2]; doubleMk[N][N]; for(m=N-1;m>=0;) { if(m==0) { Lambda_Re[0]=A[0][0]; Lambda_Im[m]=0; break; } if(fabs(A[m][m-1])<=EPS) { Lambda_Re[m]=A[m][m]; Lambda_Im[m]=0; m--; continue; } D_Coe[0]=(-1)*(A[m-1][m-1]+A[m][m]); D_Coe[1]=A[m-1][m-1]*A[m][m]-A[m-1][m]*A[m][m-1]; SOLVE_EQU(D_Coe,S_Re,S_Im); if(m==1) { Lambda_Re[0]=S_Re[0]; Lambda_Re[1]=S_Re[1]; Lambda_Im[0]=S_Im[0]; Lambda_Im[1]=S_Im[1]; break; } if(fabs(A[m-1][m-2])<=EPS) { Lambda_Re[m-1]=S_Re[0]; Lambda_Re[m]=S_Re[1]; Lambda_Im[m-1]=S_Im[0]; Lambda_Im[m]=S_Im[1]; m=m-2; continue; } MAT_MULTI((m+1),A,A,Mk); for(i=0;i<=m;i++) { for(j=0;j<=m;j++) { Mk[i][j]=Mk[i][j]+D_Coe[0]*A[i][j]; if(i==j) { Mk[i][j]=Mk[i][j]+D_Coe[1]; } } } QR_DISSOLUTION((m+1),Mk); } } voidQR_UP_HE(void) //拟上三角矩阵QR分解 { doubleR[N][N]={0},Q[N][N]={0},RQ[N][N]; chari,j,r,k; doublep[N],w[N],u[N]; doubled,c,h; FILE*fp; EVA_MAT(N,A,R); for(i=0;i { Q[i][i]=1; } for(r=0;r { for(i=r+1;i { if(R[i][r]! =0) { d=VECT_NORM(r,r,N,R); c=(-1)*SP_SGN(R[r][r])*d; h=c*c-c*R[r][r]; for(k=0;k { if(k>r) { u[k]=R[k][r]; } elseif(k==r) { u[k]=R[k][r]-c; } else { u[k]=0; } } VECT_MULTI(N,u,p,1/h); MAT_VEC(N,Q,u,w); for(k=0;k { for(j=0;j { Q[k][j]-=(w[k]*p[j]); } } MAT_T_VEC(N,R,p,w); for(k=0;k { for(j=0;j { R[k][j]-=(u[k]*w[j]); } } break; } } } if((fp=fopen("A(n-1)-QR","w"))==NULL) //输出QR分解结果 { printf("cannotopenfile\n"); exit(0); } fprintf(fp,"SY1103104胡延勍数值分析第二题\n"); fprintf(fp,"拟上三角化得到矩阵A(n-1)进行QR分解后得到: \n"); fprintf(fp,"Q矩阵为: \n"); for(i=0;i { for(j=0;j { fprintf(fp,"%.13e",Q[i][j]); } fprintf(fp,"\n\n"); } fprintf(fp,"R矩阵为: \n"); for(i=0;i { for(j=0;j { fprintf(fp,"%.13e",R[i][j]); } fprintf(fp,"\n\n"); } MAT_MULTI(N,R,Q,RQ); fprintf(fp,"乘积矩阵RQ: \n"); for(i=0;i { for(j=0;j { fprintf(fp,"%.13e",RQ[i][j]); } fprintf(fp,"\n\n"); } fclose(fp); } voidGAUSS_PIVOT(doubleLamb,doublex[]) //Gauss选主元求特征向量 { chari,k,j,line; doublem; for(i=0;i { for(k=0;k { if(i==k) { A[i][k]=Lamb-A[i][k]; } else { A[i][k]=(-1)*A[i][k]; } } } for(k=0;k { for(i=k,line=k;i
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 分析 实习 第二
![提示](https://static.bdocx.com/images/bang_tan.gif)