北航数值分析大作业2学硕Word文档下载推荐.docx
- 文档编号:16975267
- 上传时间:2022-11-27
- 格式:DOCX
- 页数:26
- 大小:330.79KB
北航数值分析大作业2学硕Word文档下载推荐.docx
《北航数值分析大作业2学硕Word文档下载推荐.docx》由会员分享,可在线阅读,更多相关《北航数值分析大作业2学硕Word文档下载推荐.docx(26页珍藏版)》请在冰豆网上搜索。
//简单QR分解函数
voidqr_decompose(doubledec_matrix[][10],doubleeps);
//求解特征值函数
voideigenvalue(doubleeg_matrix[][10],doubleeps,inttm);
//反幂法求特征向量函数
voidip_method(doubleipm_mat[][10],doublesub_epsilon,doubleip);
//矩阵LU分解函数
voidlua(doublelua_mat[][10]);
//定义两个全局变量
doubleglobala[10][10];
intM[10];
//***************************************************************************
voidmain()//主函数
{
inti,j;
intn=10;
doublea[10][10]={0};
doubleepsilon=1e-12;
for(i=0;
i<
n;
i++)//矩阵赋初值
for(j=0;
j<
j++)
{
if(i!
=j)
a[i][j]=sin(0.5*(i+1)+0.2*(j+1));
else
a[i][j]=1.5*cos((i+1)+1.2*(j+1));
}
i++)//存储全局变量,以供调用
globala[i][j]=a[i][j];
triangulate(a,epsilon);
//矩阵拟上三角化并输出
qr_decompose(a,epsilon);
//对拟上三角矩阵进行QR分解
intl=100;
eigenvalue(a,epsilon,l);
//求解特征值及其特征向量
}
*****************************************************************
voidtriangulate(doubletri_mat[][10],doubleeps)//矩阵拟上三角化
intr,i,j,t;
doubled=0,c=0,h=0;
doublev[10]={0};
doublep[10]={0};
doubleq[10]={0};
doublew[10]={0};
for(r=0;
r<
=n-3;
r++)
{
t=0;
for(i=r+2;
i++)//判断是否全为零
if(tri_mat[i][r]!
=0)
{
t=1;
break;
}
if(t==0)continue;
else
doubled_sum=0;
for(i=r+1;
i++)//计算dr
d_sum=d_sum+tri_mat[i][r]*tri_mat[i][r];
d=sqrt(d_sum);
if(tri_mat[r+1][r]<
=0)
c=d;
c=-d;
h=c*(c-tri_mat[r+1][r]);
doubleu[10]={0};
//创建向量u
u[r+1]=tri_mat[r+1][r]-c;
for(i=r+2;
i++)
u[i]=tri_mat[i][r];
for(i=0;
v[i]=u[i]/h;
i++)//求pr
doublepsum=0;
for(t=0;
t<
t++)
{
psum=psum+tri_mat[t][i]*v[t];
}
p[i]=psum;
}
i++)//求qr
doubleqsum=0;
qsum=qsum+tri_mat[i][t]*v[t];
q[i]=qsum;
doubletr=0;
//求tr
doubletsum=0;
for(t=0;
tsum=tsum+p[t]*v[t];
tr=tsum;
t++)//计算wr
w[t]=q[t]-tr*u[t];
i++)//计算A(r+1)
for(j=0;
tri_mat[i][j]=tri_mat[i][j]-w[i]*u[j]-u[i]*p[j];
}
printf("
\n****************矩阵A拟上三角化所得矩阵A(n-1)******************\n\n"
);
doubletp;
if(fabs(tri_mat[i][j])<
eps)
tp=0;
printf("
%20.11e"
tp);
tri_mat[i][j]);
printf("
\n\n"
*************************************************************************************
voiddouble_qr_decompose(doubledec_mat[][10],doubleak_mat[][10],intmt,doubleeps)//双步QR分解
inti,j,r,t;
doubled,c,h;
doubleq_mat[10][10]={0};
for(i=0,j=0;
=mt,i<
=mt;
i++,j++)
q_mat[i][j]=1;
doublew0[10]={0};
doubletr=0;
=mt-1;
for(i=r+1;
if(dec_mat[i][r]!
for(i=r;
d_sum=d_sum+dec_mat[i][r]*dec_mat[i][r];
if(dec_mat[r][r]<
h=c*(c-dec_mat[r][r]);
u[r]=dec_mat[r][r]-c;
u[i]=dec_mat[i][r];
doublev0[10]={0};
i++)//求v0r
doublev0sum=0;
v0sum=v0sum+dec_mat[t][i]*v[t];
v0[i]=v0sum;
i++)//求下一个B
dec_mat[i][j]=dec_mat[i][j]-u[i]*v0[j];
psum=psum+ak_mat[t][i]*v[t];
qsum=qsum+ak_mat[i][t]*v[t];
tr=0;
i++)//求解tr
tr=tr+p[i]*v[i];
i++)//求wr
w[i]=q[i]-tr*u[i];
//求下一个A
ak_mat[i][j]=ak_mat[i][j]-u[i]*p[j]-w[i]*u[j];
}
******************************************************************
voidqr_decompose(doubledec_matrix[][10],doubleeps)//简单QR分解
j==i&
&
doubledec_mat[10][10];
//为了保护数值,另存矩阵
dec_mat[i][j]=dec_matrix[i][j];
n-1;
i++)//求w0r
doublew0sum=0;
w0sum=w0sum+q_mat[i][t]*u[t];
w0[i]=w0sum;
i++)//求下一个Q
q_mat[i][j]=q_mat[i][j]-w0[i]*v[j];
psum=psum+dec_mat[t][i]*v[t];
i++)//求下一个A
dec_mat[i][j]=dec_mat[i][j]-u[i]*p[j];
doubleqr[10][10]={0};
i++)//计算RQ
doubleqrsum=0;
qrsum=qrsum+dec_mat[i][t]*q_mat[t][j];
qr[i][j]=qrsum;
*************矩阵A(n-1)进行简单QR分解所得矩阵Q******************\n\n"
i++)//输出矩阵Q
if(fabs(q_mat[i][j])<
eps)
printf("
q_mat[i][j]);
*************矩阵A(n-1)进行简单QR分解所得矩阵R******************\n\n"
i++)//输出矩阵R
if(fabs(dec_mat[i][j])<
dec_mat[i][j]);
**********************乘积矩阵RQ***********************\n\n"
i++)//输出矩阵RQ
if(fabs(qr[i][j])<
qr[i][j]);
***************************************************************************
voideigenvalue(doubleeg_matrix[][10],doubleeps,inttm)//求取特征值与特征向量
inti,j,k,m,t;
doubledet;
doubleeg_mat[10][10];
//存储矩阵
doubleegv[10][2]={0};
//存储矩阵特征值
doublesk=0,tk=0;
doubles[2][2];
doubledelta;
eg_mat[i][j]=eg_matrix[i][j];
doublematrix[10][10]={0};
doublesum;
m=n-1;
for(k=1;
;
k++)
if(fabs(eg_mat[m][m-1])<
=eps)
egv[m][0]=eg_mat[m][m];
egv[m][1]=0;
m=m-1;
if(m==0)
egv[0][0]=eg_mat[0][0];
egv[0][1]=0;
else
if(m>
0)continue;
det=eg_mat[m-1][m-1]*eg_mat[m][m]-eg_mat[m-1][m]*eg_mat[m][m-1];
delta=(eg_mat[m][m]+eg_mat[m-1][m-1])*(eg_mat[m][m]+eg_mat[m-1][m-1])-4*det;
2;
i++)//存储二阶子阵特征值
s[i][j]=0;
if(delta>
s[0][0]=((eg_mat[m][m]+eg_mat[m-1][m-1])+sqrt(delta))/2;
s[1][0]=((eg_mat[m][m]+eg_mat[m-1][m-1])-sqrt(delta))/2;
delta=-delta;
s[0][0]=(eg_mat[m][m]+eg_mat[m-1][m-1])/2;
s[0][1]=sqrt(delta)/2;
s[1][0]=(eg_mat[m][m]+eg_mat[m-1][m-1])/2;
s[1][1]=-sqrt(delta)/2;
if(m==1)
egv[i][j]=s[i][j];
break;
if(fabs(eg_mat[m-1][m-2])<
for(i=0;
for(j=0;
{
egv[m-1+i][j]=s[i][j];
}
m=m-2;
if(m==0)
egv[0][0]=eg_mat[0][0];
egv[0][1]=0;
break;
else
if(m>
if(k==tm)
printf("
计算中止,未得到全部特征值!
\n"
else
sk=eg_mat[m][m]+eg_mat[m-1][m-1];
tk=eg_mat[m][m]*eg_mat[m-1][m-1]-eg_mat[m][m-1]*eg_mat[m-1][m];
for(i=0;
=m;
i++)//Mk
for(j=0;
{
sum=0;
for(t=0;
t++)
sum=sum+eg_mat[i][t]*eg_mat[t][j];
matrix[i][j]=sum-sk*eg_mat[i][j];
if(i==j)
matrix[i][j]=matrix[i][j]+tk;
}
double_qr_decompose(matrix,eg_mat,m,eps);
//双步位移QR分解
if(k<
tm)
特征值为:
*******************实部*****************************虚部************\n"
lambda[%d][%d]=%17.11e"
i+1,j+1,egv[i][j]);
if(egv[i][1]==0)
属于特征值λ[%d]=%17.11e的特征向量为:
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 北航 数值 分析 作业