数值分析带双步位移的QR分解求特征值算法Word文件下载.docx
- 文档编号:21504756
- 上传时间:2023-01-30
- 格式:DOCX
- 页数:15
- 大小:66.83KB
数值分析带双步位移的QR分解求特征值算法Word文件下载.docx
《数值分析带双步位移的QR分解求特征值算法Word文件下载.docx》由会员分享,可在线阅读,更多相关《数值分析带双步位移的QR分解求特征值算法Word文件下载.docx(15页珍藏版)》请在冰豆网上搜索。
线性方程组的求解采用列主元素Gauss消去法。
#include<
stdio.h>
math.h>
#defineERR1.0e-12//误差限
#defineN10//矩阵行列数
#defineL1.0e5//最大迭代次数
doubleA[N][N]={0};
voidInit_A()//初始化矩阵
{
inti,j;
for(i=0;
i<
N;
i++)
for(j=0;
j<
j++)
{
if(i==j)
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));
}
}
/*
voidDisplay_A()
{
printf("
%8.4f"
A[i][j]);
printf("
\n"
);
}
*/
intSgn(doublea)
if(a>
=0)
return1;
else
return-1;
voidOn_To_The_Triangle()//矩阵拟上三角化
inti,j,r,flag=0;
doublecr,dr,hr,tr,temp;
doubleur[N],pr[N],qr[N],wr[N];
for(r=1;
r<
=N-2;
r++)
flag=0;
for(i=r+2;
=N;
if(A[i-1][r-1]!
{
flag=1;
break;
}
if(0==flag)
continue;
dr=0;
for(i=r+1;
dr+=A[i-1][r-1]*A[i-1][r-1];
dr=sqrt(dr);
if(0==A[r][r-1])
cr=dr;
elsecr=-Sgn(A[r][r-1])*dr;
hr=cr*cr-cr*A[r][r-1];
for(i=1;
=r;
ur[i-1]=0;
ur[r]=A[r][r-1]-cr;
ur[i-1]=A[i-1][r-1];
temp=0;
for(j=1;
temp+=A[j-1][i-1]*ur[j-1];
pr[i-1]=temp/hr;
temp+=A[i-1][j-1]*ur[j-1];
qr[i-1]=temp/hr;
temp=0;
temp+=pr[i-1]*ur[i-1];
tr=temp/hr;
wr[i-1]=qr[i-1]-tr*ur[i-1];
A[i-1][j-1]=A[i-1][j-1]-wr[i-1]*ur[j-1]-ur[i-1]*pr[j-1];
voidGet_Roots(doubleeigenvalue[][2],intm,doubless,doublett)//求一元二次方程的根
doublediscriminant=ss*ss-4*tt;
//
if(discriminant<
0)
*(*(eigenvalue+m-2))=0.5*ss;
*(*(eigenvalue+m-2)+1)=0.5*sqrt(-discriminant);
*(*(eigenvalue+m-1))=0.5*ss;
*(*(eigenvalue+m-1)+1)=-0.5*sqrt(-discriminant);
else
*(*(eigenvalue+m-2))=0.5*(ss+sqrt(discriminant));
*(*(eigenvalue+m-2)+1)=0;
*(*(eigenvalue+m-1))=0.5*(ss-sqrt(discriminant));
*(*(eigenvalue+m-1)+1)=0;
voidGet_Mk(doublemk[][N],intm,doubless,doublett)//获取Mk,用于带双步位移的QR分解
inti,j,k;
m;
*(*(mk+i)+j)=0;
for(k=0;
k<
k++)
*(*(mk+i)+j)+=A[i][k]*A[k][j];
*(*(mk+i)+j)-=ss*A[i][j];
if(j==i)
*(*(mk+i)+j)+=tt;
voidQR_Reslove(doublemk[][N],intm)//QR分解
doubleur[N],vr[N],pr[N],qr[N],wr[N];
doubleB[N][N],C[N][N];
{
B[i][j]=*(*(mk+i)+j);
C[i][j]=A[i][j];
=m-1;
=m;
if(B[i-1][r-1]!
for(i=r;
dr+=B[i-1][r-1]*B[i-1][r-1];
if(0==B[r-1][r-1])
elsecr=-Sgn(B[r-1][r-1])*dr;
hr=cr*cr-cr*B[r-1][r-1];
r;
ur[r-1]=B[r-1][r-1]-cr;
ur[i-1]=B[i-1][r-1];
temp+=B[j-1][i-1]*ur[j-1];
vr[i-1]=temp/hr;
for(i=0;
for(j=0;
B[i][j]-=ur[i]*vr[j];
temp+=C[j-1][i-1]*ur[j-1];
temp+=C[i-1][j-1]*ur[j-1];
C[i-1][j-1]=C[i-1][j-1]-wr[i-1]*ur[j-1]-ur[i-1]*pr[j-1];
A[i][j]=C[i][j];
voidDisplay_Eigenvalue(doublevalue[][2])//显示特征值
inti;
{
λ%d=%8.4f"
i+1,*(*(value+i)));
if(*(*(value+i)+1)>
+%8.4f"
*(*(value+i)+1));
elseif(*(*(value+i)+1)<
printf("
intQR_With_Double_Step_Displacement(doubleeigenvalue[][2])//带双步位移QR分解求特征值
inti,j,k=1,m=N;
doubles,t;
doubleMk[N][N];
2;
eigenvalue[i][j]=0;
do
k++;
if(m==1)
eigenvalue[m-1][0]=A[m-1][m-1];
m--;
elseif(m==2)
s=A[m-2][m-2]+A[m-1][m-1];
t=A[m-2][m-2]*A[m-1][m-1]-A[m-1][m-2]*A[m-1][m-2];
Get_Roots(eigenvalue,m,s,t);
//求一元二次方程的根
m=0;
elseif(m==0)
return0;
elseif(fabs(A[m-1][m-2])<
=ERR)
else
t=A[m-2][m-2]*A[m-1][m-1]-A[m-1][m-2]*A[m-2][m-1];
if(fabs(A[m-2][m-3])<
Get_Roots(eigenvalue,m,s,t);
m-=2;
continue;
else
Get_Mk(Mk,m,s,t);
//获取Mk,用于带双步位移的QR分解
QR_Reslove(Mk,m);
//QR分解
while(k<
L);
Errer\n"
//迭代过程不成功,报错
voidSelect_Principal_Element(intk)//Gauss消元法求解方程组之选主元
inti,max=k;
doubletrans[N];
for(i=k+1;
if(fabs(A[i][k])>
fabs(A[max][k]))
max=i;
if(k==max)
return;
for(i=k;
trans[i]=A[k][i];
A[k][i]=A[max][i];
A[max][i]=trans[i];
voidEliminant(intk)//Gauss消元法求解方程组之消元
doublerate;
rate=A[i][k]/A[k][k];
for(j=k;
A[i][j]=A[i][j]-rate*A[k][j];
voidBack_Substitution(doubleEigenvector[])//Gauss消元法求解方程组之回代
doubleTemp_M[N][N+1];
N-1;
for(j=i;
Temp_M[i][j]=A[i][j];
Temp_M[i][N]=0;
Temp_M[N-1][N-1]=1;
Temp_M[N-1][N]=1;
for(k=N-1;
k>
=0;
k--)
*(Eigenvector+k)=Temp_M[k][N];
*(Eigenvector+k)=*(Eigenvector+k)-*(Eigenvector+i+1)*Temp_M[k][i+1];
*(Eigenvector+k)=*(Eigenvector+k)/Temp_M[k][k];
voidDisplay_Eigenvector(doubleEigenvector[],doubleeigenvalue)//显示特征值对应特征向量
Whenλ=%8.4f\nEigenvector=\n"
eigenvalue);
*(Eigenvector+i));
voidGet_Eigenvector(doublevalue[][2])//利用Gauss消元法求解特征向量
doubleeigenvalue;
doubleEigenvector[N];
if(*(*(value+i)+1)==0)
Init_A();
//初始化矩阵
eigenvalue=*(*(value+i));
A[j][j]-=eigenvalue;
Select_Principal_Element(j);
//Gauss消元法求解方程组之选主元
Eliminant(j);
//Gauss消元法求解方程组之消元
Back_Substitution(Eigenvector);
//Gauss消元法求解方程组之回代
Display_Eigenvector(Eigenvector,eigenvalue);
//显示特征值对应特征向量
main()
doubleeigenvalue[N][2];
Init_A();
//初始化矩阵
On_To_The_Triangle();
//矩阵上三角化
QR_With_Double_Step_Displacement(eigenvalue);
//带双步位移QR分解求特征值
Contactme:
****************\n\n"
Display_Eigenvalue(eigenvalue);
//显示特征值
Get_Eigenvector(eigenvalue);
//利用Gauss消元法求解特征向量
return0;
PS:
Forfurtherdetails,pleasecontactme:
****************
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 分析 带双步 位移 QR 分解 特征值 算法