西华数学与计算机学院上机实践报告
课程名称：计算方法A
指导教师：严常龙
上机实践名称：Jacobi法求实对称矩阵特征值
上机实践编号：7
年级：2010级 姓名： 学号： 上机实践成绩： 上机实践日期：yyyy.mm. dd 上机实践时间： 一、目的
1．通过本实验加深对Jacobi法的构造过程的理解；
2．能对Jacobi法提出正确的算法描述编程实现，得到计算结果。
二、内容与设计思想
自选方阵，用Jacobi法求其全部特征值及其特征向量。
可使用实例：
1A2
328939 8
三、使用环境
操作系统：
软件平台：
四、核心代码及调试过程
#include
#include
#include
#define n 3
#define e 0.000001
main()
{
int i,j,p,q,cj=0,k;
float a[n][n],sum=0.0,s,max,t,t1,t2,c,d,v[n],L[n][n],L1[n][n],L2[n][n];
a[0][0]=1;a[0][1]=2;a[0][2]=3;
a[1][0]=2;a[1][1]=8;a[1][2]=9;
a[2][0]=3;a[2][1]=9;a[2][2]=8;
for(i=0;i
{
for(j=0;j
{
if(j==i)
{
L[i][j]=1;
}
else
{
L[i][j]=0;
}
}
}
memcpy(L2,L,sizeof(L2));//复制数组 memcpy(L1,L,sizeof(L1));
printf("矩阵A为\n");
for(i=0;i
{
for(j=0;j
{
printf("%f\t",a[i][j]);
}
printf("\n");
}
sum=0.0;
for(i=0;i
{
for(j=0;j
{
if(i!=j)
{
sum=sum+a[i][j]*a[i][j];
}
}
}
while(sum>e)
{
max=fabs(a[0][1]);p=0;q=1;
for(i=0;i
for(j=0;j
{
if(i!=j)
{
if(max
{
max=fabs(a[i][j]); p=i;q=j;
}
}
}
}
s=(a[q][q]-a[p][p])/(a[p][q]*2.0);
if(s==0)
{
t=1;
}
else
{
t1=-s-sqrt(s*s+1);
t2=-s+sqrt(s*s+1);
if(fabs(t1)>fabs(t2))
{
t=t2;
}
else
{
t=t1;
}
}
c=1.0/sqrt(1+t*t);
d=t/sqrt(1+t*t);
for(i=0;i
{
if(i!=p && i!=q)
{
a[i][p]=c*a[p][i]-d*a[q][i]; a[i][q]=c*a[q][i]+d*a[p][i]; a[q][i]=a[i][q];
a[p][i]=a[i][p];
}
}
a[p][p]=a[p][p]-t*a[p][q];
a[q][q]=a[q][q]+t*a[p][q];
a[p][q]=0;
a[q][p]=0;
sum=0.0;
for(i=0;i
{
for(j=0;j
{
if(i!=j)
{
sum=sum+a[i][j]*a[i][j]; }
}
}
L2[p][p]=c;
L2[q][q]=c;
L2[p][q]=d;
L2[q][p]=-d;
//L=L1*L2;
//L1=L;
for(i=0;i
{
for(j=0;j
{
L[i][j]=0;
for(k=0;k
{
L[i][j]=L[i][j]+L1[i][k]*L2[k][j];
}
}
}
memcpy(L1,L,sizeof(L1));
}
printf("矩阵B为\n");
for(i=0;i
{
for(j=0;j
{
printf("%f\t",L[i][j]);
}
printf("\n");
}
printf("正交矩阵为\n"); for(i=0;i
{
for(j=0;j
{
printf("%f\t",a[i][j]);
}
printf("\n");
}
printf("特征值为\n"); for(i=0;i
{
v[i]=a[i][i];
printf("v[%d]=%f\n",i+1,v[i]);
}
printf("特征向量为\n");
for(i=0;i
{
printf("t[%d]=(\t",i+1); for(j=0;j
{
printf("%f\t",L[j][i]); }
printf(")T");
printf("\n");
}
}
五、总结
六、附录
西华数学与计算机学院上机实践报告
课程名称：计算方法A
指导教师：严常龙
上机实践名称：Jacobi法求实对称矩阵特征值
上机实践编号：7
年级：2010级 姓名： 学号： 上机实践成绩： 上机实践日期：yyyy.mm. dd 上机实践时间： 一、目的
1．通过本实验加深对Jacobi法的构造过程的理解；
2．能对Jacobi法提出正确的算法描述编程实现，得到计算结果。
二、内容与设计思想
自选方阵，用Jacobi法求其全部特征值及其特征向量。
可使用实例：
1A2
328939 8
三、使用环境
操作系统：
软件平台：
四、核心代码及调试过程
#include
#include
#include
#define n 3
#define e 0.000001
main()
{
int i,j,p,q,cj=0,k;
float a[n][n],sum=0.0,s,max,t,t1,t2,c,d,v[n],L[n][n],L1[n][n],L2[n][n];
a[0][0]=1;a[0][1]=2;a[0][2]=3;
a[1][0]=2;a[1][1]=8;a[1][2]=9;
a[2][0]=3;a[2][1]=9;a[2][2]=8;
for(i=0;i
{
for(j=0;j
{
if(j==i)
{
L[i][j]=1;
}
else
{
L[i][j]=0;
}
}
}
memcpy(L2,L,sizeof(L2));//复制数组 memcpy(L1,L,sizeof(L1));
printf("矩阵A为\n");
for(i=0;i
{
for(j=0;j
{
printf("%f\t",a[i][j]);
}
printf("\n");
}
sum=0.0;
for(i=0;i
{
for(j=0;j
{
if(i!=j)
{
sum=sum+a[i][j]*a[i][j];
}
}
}
while(sum>e)
{
max=fabs(a[0][1]);p=0;q=1;
for(i=0;i
for(j=0;j
{
if(i!=j)
{
if(max
{
max=fabs(a[i][j]); p=i;q=j;
}
}
}
}
s=(a[q][q]-a[p][p])/(a[p][q]*2.0);
if(s==0)
{
t=1;
}
else
{
t1=-s-sqrt(s*s+1);
t2=-s+sqrt(s*s+1);
if(fabs(t1)>fabs(t2))
{
t=t2;
}
else
{
t=t1;
}
}
c=1.0/sqrt(1+t*t);
d=t/sqrt(1+t*t);
for(i=0;i
{
if(i!=p && i!=q)
{
a[i][p]=c*a[p][i]-d*a[q][i]; a[i][q]=c*a[q][i]+d*a[p][i]; a[q][i]=a[i][q];
a[p][i]=a[i][p];
}
}
a[p][p]=a[p][p]-t*a[p][q];
a[q][q]=a[q][q]+t*a[p][q];
a[p][q]=0;
a[q][p]=0;
sum=0.0;
for(i=0;i
{
for(j=0;j
{
if(i!=j)
{
sum=sum+a[i][j]*a[i][j]; }
}
}
L2[p][p]=c;
L2[q][q]=c;
L2[p][q]=d;
L2[q][p]=-d;
//L=L1*L2;
//L1=L;
for(i=0;i
{
for(j=0;j
{
L[i][j]=0;
for(k=0;k
{
L[i][j]=L[i][j]+L1[i][k]*L2[k][j];
}
}
}
memcpy(L1,L,sizeof(L1));
}
printf("矩阵B为\n");
for(i=0;i
{
for(j=0;j
{
printf("%f\t",L[i][j]);
}
printf("\n");
}
printf("正交矩阵为\n"); for(i=0;i
{
for(j=0;j
{
printf("%f\t",a[i][j]);
}
printf("\n");
}
printf("特征值为\n"); for(i=0;i
{
v[i]=a[i][i];
printf("v[%d]=%f\n",i+1,v[i]);
}
printf("特征向量为\n");
for(i=0;i
{
printf("t[%d]=(\t",i+1); for(j=0;j
{
printf("%f\t",L[j][i]); }
printf(")T");
printf("\n");
}
}
五、总结
六、附录