数值分析第二次作业Word文档下载推荐.docx
- 文档编号:17652974
- 上传时间:2022-12-07
- 格式:DOCX
- 页数:24
- 大小:129.46KB
数值分析第二次作业Word文档下载推荐.docx
《数值分析第二次作业Word文档下载推荐.docx》由会员分享,可在线阅读,更多相关《数值分析第二次作业Word文档下载推荐.docx(24页珍藏版)》请在冰豆网上搜索。
//初始化矩阵
voidcreateMatrix(doubleA[][MAXLEN+1])
for(inti=1;
i<
=MAXLEN;
i++)
for(intj=1;
j<
j++)
{
if(i!
=j)
A[i][j]=sin(0.5*i+0.2*j);
else
A[i][j]=1.52*cos(i+1.2*j);
}
}
//打印矩阵
voidprintMatrix(doubleA[][MAXLEN+1])
cout<
<
A[i][j]<
"
"
;
endl;
intsgn(doublevalue)
if(value>
0)
return1;
if(value<
=0)
return-1;
//矩阵的拟上三角化
double(*triangulation(doubleA[][MAXLEN+1]))[MAXLEN+1]
double(*tempA)[MAXLEN+1]=newdouble[MAXLEN+1][MAXLEN+1];
tempA[i][j]=A[i][j];
for(intr=1;
r<
=MAXLEN-2;
r++)
boolcanContinue=true;
for(inti=r+2;
if(tempA[i][r]!
=0){
canContinue=false;
break;
if(canContinue)
continue;
doubledr=0,cr=0,hr=0;
for(inti=r+1;
dr+=tempA[i][r]*tempA[i][r];
dr=sqrt(dr);
cr=-sgn(tempA[r+1][r])*dr;
hr=cr*cr-cr*tempA[r+1][r];
doubleur[MAXLEN+1],pr[MAXLEN+1],qr[MAXLEN+1],wr[MAXLEN+1],tr;
r+1;
ur[i]=0;
ur[i]=tempA[i][r];
ur[r+1]-=cr;
tr=0;
pr[i]=0;
qr[i]=0;
pr[i]+=tempA[j][i]*ur[j];
qr[i]+=tempA[i][j]*ur[j];
pr[i]=pr[i]/hr;
qr[i]=qr[i]/hr;
tr+=pr[i]*ur[i];
tr=tr/hr;
wr[i]=qr[i]-tr*ur[i];
tempA[i][j]=tempA[i][j]-wr[i]*ur[j]-ur[i]*pr[j];
if(fabs(tempA[i][j])<
THTA)
tempA[i][j]=0;
returntempA;
//矩阵相乘
double(*matrixMut(doubleA[][MAXLEN+1],doubleB[][MAXLEN+1]))[MAXLEN+1]
double(*result)[MAXLEN+1]=newdouble[MAXLEN+1][MAXLEN+1];
result[i][j]=0;
for(intk=1;
k<
k++)
result[i][j]+=A[i][k]*B[k][j];
returnresult;
//矩阵转置
double(*matrixTrans(doubleA[][MAXLEN+1]))[MAXLEN+1]
result[j][i]=A[i][j];
//矩阵复制
double(*matrixCopy(doubleA[][MAXLEN+1]))[MAXLEN+1]
result[i][j]=A[i][j];
//QR方法
voidQRMethod(doubleA[][MAXLEN+1],double(*&
Q)[MAXLEN+1],double(*&
R)[MAXLEN+1],intstep)
Q=newdouble[MAXLEN+1][MAXLEN+1];
R=newdouble[MAXLEN+1][MAXLEN+1];
if(i==j)
Q[i][j]=1;
Q[i][j]=0;
R[i][j]=A[i][j];
=step;
if(R[i][r]!
for(inti=r;
dr+=R[i][r]*R[i][r];
cr=-sgn(R[r][r])*dr;
hr=cr*cr-cr*R[r][r];
doubleur[MAXLEN+1],wr[MAXLEN+1],pr[MAXLEN+1],tr;
r;
ur[i]=R[i][r];
ur[r]=R[r][r]-cr;
wr[i]=0;
wr[i]+=Q[i][j]*ur[j];
Q[i][j]=Q[i][j]-wr[i]*ur[j]/hr;
pr[i]+=R[j][i]*ur[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)//第四步
returntrue;
if(m<
if(m>
1)
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>
s1.real=(-b-sqrt(delta))/(2*a);
s1.imag=0;
s2.real=(-b+sqrt(delta))/(2*a);
s2.imag=0;
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;
if(fabs(A[m-1][m-2])<
=THTA)//第七步
m=m-2;
if(m==1)//转第四步
if(m==0)
if(k>
=L)//第八步
returnfalse;
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];
if((i>
m)||(j>
m))
AK[i][j]=0;
AK[i][j]=A[i][j];
double(*AK2)[MAXLEN+1]=matrixMut(AK,AK);
m)){
MK[i][j]=0;
MK[i][j]=AK2[i][j]-s*AK[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[])
MAXLEN;
inti=k;
for(intj=k;
if(mat[j][k]>
mat[i][k])
i=j;
doubletemp=mat[k][j];
mat[k][j]=mat[i][j];
mat[i][j]=temp;
for(inti=k+1;
doublemik=mat[i][k]/mat[k][k];
for(intj=k+1;
mat[i][j]=mat[i][j]-mik*mat[k][j];
x[MAXLEN]=1;
for(intk=MAXLEN-1;
k>
=1;
k--)
x[k]=0;
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);
矩阵A的拟上三角化"
printMatrix(ssjResult);
double(*Q)[MAXLEN+1]=NULL,(*R)[MAXLEN+1]=NULL;
QRMethod(ssjResult,Q,R,1);
拟上三角矩阵A的QR分解第一步"
R*Q"
double(*QR)[MAXLEN+1];
QR=matrixMut(R,Q);
printMatrix(QR);
Compflags[MAXLEN+1];
A^(n-1)和R*Q的差值"
ssjResult[i][j]-QR[i][j]<
if(twoStepQR(matrixCopy(ssjResult),ITIME,flags))
λ"
="
flags[i].real<
if(flags[i].imag>
+"
flags[i].imag;
elseif(flags[i].imag<
if(flags[i].imag==0)
doubletemp[MAXLEN+1][MAXLEN+1];
doublex[MAXLEN+1];
for(intk=0;
temp[j][k]=-A[j][k];
if(j==k)
temp[j][k]+=flags[i].real;
solveEqua(temp,x);
实特征值λ"
的特征向量"
x[j]<
endl<
Fail!
return0;
程序输出
计算结果
矩阵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+000-7.72897513499E-001-1.60102824405E+000
-2.91268547483E-001-2.43433785832E-0016.73628608451E-001
2.62477290494E-001
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
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 分析 第二次 作业