ur[i]=0;
for(inti=r+1;i<=MAXLEN;i++)
ur[i]=R[i][r];
ur[r]=R[r][r]-cr;
for(inti=1;i<=MAXLEN;i++)
{
wr[i]=0;
for(intj=1;j<=MAXLEN;j++)
wr[i]+=Q[i][j]*ur[j];
}
for(inti=1;i<=MAXLEN;i++)
for(intj=1;j<=MAXLEN;j++)
Q[i][j]=Q[i][j]-wr[i]*ur[j]/hr;
for(inti=1;i<=MAXLEN;i++)
{
pr[i]=0;
for(intj=1;j<=MAXLEN;j++)
pr[i]+=R[j][i]*ur[j];
pr[i]=pr[i]/hr;
}
for(inti=1;i<=MAXLEN;i++)
for(intj=1;j<=MAXLEN;j++)
R[i][j]=R[i][j]-ur[i]*pr[j];
}
}
//双步位移的QR方法
booltwoStepQR(doubleA[][MAXLEN+1],intL,Compflags[])
{
intk=1,m=MAXLEN;//第二步
intnum=1;
while(true)
{
if(fabs(A[m][m-1])<=THTA)//第三步
{
flags[num].real=A[m][m];
flags[num++].imag=0;
m=m-1;
if(m==1)//第四步
{
flags[num].real=A[m][m];
flags[num++].imag=0;
returntrue;
}
if(m<=0)
returntrue;
if(m>1)
continue;
}
else
{
doublea=1;//第五步
doubleb=-(A[m-1][m-1]+A[m][m]);
doublec=A[m-1][m-1]*A[m][m]-A[m-1][m]*A[m][m-1];
doubledelta=b*b-4*a*c;
Comps1,s2;
if(delta>=0)
{
s1.real=(-b-sqrt(delta))/(2*a);
s1.imag=0;
s2.real=(-b+sqrt(delta))/(2*a);
s2.imag=0;
}
else
{
s1.real=(-b)/(2*a);
s1.imag=-sqrt(-delta)/(2*a);
s2.real=(-b)/(2*a);
s2.imag=sqrt(-delta)/(2*a);
}
if(m==2)//第六步
{
flags[num++]=s1;
flags[num++]=s2;
returntrue;
}
else
{
if(fabs(A[m-1][m-2])<=THTA)//第七步
{
flags[num++]=s1;
flags[num++]=s2;
m=m-2;
if(m==1)//转第四步
{
flags[num].real=A[m][m];
flags[num++].imag=0;
returntrue;
}
if(m==0)
returntrue;
if(m>1)
continue;
}
else
{
if(k>=L)//第八步
returnfalse;
else
{
doubles=A[m-1][m-1]+A[m][m];//第九步
doublet=A[m-1][m-1]*A[m][m]-A[m-1][m]*A[m][m-1];
double(*MK)[MAXLEN+1]=newdouble[MAXLEN+1][MAXLEN+1];
double(*AK)[MAXLEN+1]=newdouble[MAXLEN+1][MAXLEN+1];
for(inti=1;i<=MAXLEN;i++)
{
for(intj=1;j<=MAXLEN;j++)
{
if((i>m)||(j>m))
AK[i][j]=0;
else
AK[i][j]=A[i][j];
}
}
double(*AK2)[MAXLEN+1]=matrixMut(AK,AK);
for(inti=1;i<=MAXLEN;i++)
{
for(intj=1;j<=MAXLEN;j++)
{
if((i>m)||(j>m)){
MK[i][j]=0;
continue;
}
MK[i][j]=AK2[i][j]-s*AK[i][j];
if(i==j)
MK[i][j]+=t;
}
}
double(*Q)[MAXLEN+1],(*R)[MAXLEN+1],(*QT)[MAXLEN+1],(*temp)[MAXLEN+1];
QRMethod(MK,Q,R,MAXLEN-1);
QT=matrixTrans(Q);
temp=matrixMut(QT,AK);
delete[]A;
A=matrixMut(temp,Q);
delete[]AK2;
delete[]MK;
delete[]Q;
delete[]R;
delete[]QT;
delete[]temp;
delete[]AK;
k++;
}
}
}
}
}
}
//解线性方程组
voidsolveEqua(double(*mat)[MAXLEN+1],doublex[])
{
for(intk=1;k{
inti=k;
for(intj=k;j<=MAXLEN;j++)
{
if(mat[j][k]>mat[i][k])
i=j;
}
for(intj=k;j<=MAXLEN;j++)
{
doubletemp=mat[k][j];
mat[k][j]=mat[i][j];
mat[i][j]=temp;
}
for(inti=k+1;i<=MAXLEN;i++)
{
doublemik=mat[i][k]/mat[k][k];
for(intj=k+1;j<=MAXLEN;j++)
mat[i][j]=mat[i][j]-mik*mat[k][j];
}
}
x[MAXLEN]=1;
for(intk=MAXLEN-1;k>=1;k--)
{
x[k]=0;
for(intj=k+1;j<=MAXLEN;j++)
x[k]-=mat[k][j]*x[j];
x[k]=x[k]/mat[k][k];
}
}
intmain()
{
doubleA[MAXLEN+1][MAXLEN+1],(*ssjResult)[MAXLEN+1];
createMatrix(A);
cout.precision(11);
cout.setf(ios:
:
scientific|ios:
:
uppercase);
ssjResult=triangulation(A);
cout<<"矩阵A的拟上三角化"<printMatrix(ssjResult);
double(*Q)[MAXLEN+1]=NULL,(*R)[MAXLEN+1]=NULL;
QRMethod(ssjResult,Q,R,1);
cout<cout<<"拟上三角矩阵A的QR分解第一步"<cout<<"R*Q"<double(*QR)[MAXLEN+1];
QR=matrixMut(R,Q);
printMatrix(QR);
Compflags[MAXLEN+1];
cout<cout<<"A^(n-1)和R*Q的差值"<for(inti=1;i<=MAXLEN;i++)
{
for(intj=1;j<=MAXLEN;j++)
cout<cout<}
cout<if(twoStepQR(matrixCopy(ssjResult),ITIME,flags))
{
for(inti=1;i<=MAXLEN;i++)
{
cout<<"λ"<
if(flags[i].imag>0)
cout<<"+"<elseif(flags[i].imag<0)
cout<cout<}
cout<for(inti=1;i<=MAXLEN;i++)
{
if(flags[i].imag==0)
{
doubletemp[MAXLEN+1][MAXLEN+1];
doublex[MAXLEN+1];
for(intj=1;j<=MAXLEN;j++)
{
for(intk=0;k<=MAXLEN;k++)
{
temp[j][k]=-A[j][k];
if(j==k)
temp[j][k]+=flags[i].real;
}
}
solveEqua(temp,x);
cout<<"实特征值λ"<
for(intj=1;j<=MAXLEN;j++)
cout<cout<}
}
}
else
{
cout<<"Fail!
"<}
return0;
}
程序输出
1.矩阵A的拟上三角化
2.对
执行QR方法第一步,并打印第一步做完后的矩阵RQ
3.用双步位移的QR方法计算矩阵A的全部特征值
4.A的相应于实特征值的特征向量
计算结果
1.矩阵A的拟上三角化
矩阵A得到的拟上三角化矩阵为:
-8.94521698228E-001-9.93313649183E-002-1.09983175888E+000
-7.66503870908E-0011.70760114146E-001-1.93488255889E+000
-8.39020870525E-0029.13256511314E-001-6.40797700919E-001
1.94673367868E-001
-2.34787836242E+0002.37205792160E+0001.82799855232E+000
3.26655688471E-0012.08236058364E-0012.08898700994E+000
1.84786191029E-001-1.26301526608E+0006.79069466850E-001
-4.67215088650E-001
0.00000000000E+0001.73595446995E+000-1.16502336748E+000
-1.24674444352E+000-6.29822548908E-001-1.98482018099E+000
2.97575006080E-0016.33930059659E-001-1.30851892877E-001
3.04030103610E-001
0.00000000000E+0000.00000000000E+000-1.29293756392E+000
-1.12623922590E+0001.19078291192E+000-1.30877298390E+000
1.86015166267E-0014.23673393688E-001-1.01960082655E-001
1.94366091451E-001
0.00000000000E+0000.00000000000E+0000.00000000000E+000
1.57771115303E+0008.16935832816E-0014.46153172383E-001
-4.36509254161E-002-4.66597916719E-0012.94123156618E-001
-1.03442111366E-001
0.00000000000E+0000.00000000000E+0000.00000000000E+000
0.00000000000E+000-7.72897513499E-001-1.60102824405E+000
-2.91268547483E-001-2.43433785832E-0016.73628608451E-001
2.62477290494E-001
0.00000000000E+0000.00000000000E+0000.00000000000E+000
0.00000000000E+0000.00000000000E+000-7.29677394636E-001
-7.96545627983E-0039.71073910201E-001-1.29896736857E-001
2.78024208124E-002
0.00000000000E+0000.00000000000E+0000.00000000000E+000
0.00000000000E+0000.00000000000E+0000.00000000000E+000