📄 homework2.c
字号:
//作业2//
#include<stdio.h>
#include<math.h>
#define N 10
//全局变量定义
double a[N][N],t[N][N],s[2][2];
double eps=1e-12;
int m,k;
//子函数定义
void hessenberg();
void DK();
void QR();
void two_step_move_QR();
void guass(double y);
//主函数
main()
{
double lamda[N][2];
double L=100;
int i,j;
//矩阵A初始化
for(i=0;i<N;i++)
{
for(j=0;j<N;j++)
a[i][j]=sin(0.5*(i+1)+0.2*(j+1));
a[i][i]=1.5*cos((i+1)+1.2*(i+1));
}
//存储矩阵A的初始值
for(i=0;i<N;i++)
for(j=0;j<N;j++)
t[i][j]=a[i][j];
//1.对矩阵进行上三角化
hessenberg();
//对拟上三角化后的矩阵A(n-1)进行QR分解
QR();
//2.设k，m初值
k=1;
m=N-1;
//3.判断a[m][m-1]是否为0，如果是则得到一个特征值，降阶计算，转4，否则转5
lable3:
if(fabs(a[m][m-1])<eps)
{
lamda[m][0]=a[m][m];
lamda[m][1]=0;
m=m-1;
goto lable4;
}
else
goto lable5;
//4.如果m为0，则得到一个特征值转11，如果m>0,则转3，否则直接转11
lable4:
if(m==0)
{
lamda[0][0]=a[0][0];
lamda[0][1]=0;
goto lable11;
}
else if (m>0)
goto lable3;
else
goto lable11;
//5.求二阶子阵的两个特征根
lable5:
DK();
//6.如果m=1，则得到A的两个特征值，转11，否则转7
if(m==1)
{
lamda[m-1][0]=s[0][0];
lamda[m-1][1]=s[0][1];
lamda[m][0]=s[1][0];
lamda[m][1]=s[1][1];
goto lable11;
}
else goto lable7;
//7.如果a[m-1][m-2]等于0，则得到A的两个特征值，降阶，转4，否则转5
lable7:
if(fabs(a[m-1][m-2])<eps)
{
lamda[m-1][0]=s[0][0];
lamda[m-1][1]=s[0][1];
lamda[m][0]=s[1][0];
lamda[m][1]=s[1][1];
m=m-2;
goto lable4;
}
else goto lable8;
//8.如果k到了最大迭代次数L，计算终止，未得到全部特征值，否则转9
lable8:
if(k==L)
printf("未求到全部特征值！\n\n");
else goto lable9;
//9.对A进行QR分解，得上三角形式
lable9:
two_step_move_QR();
//10.置k=k+1，转3
k=k+1;
goto lable3;
//11.A的全部特征值计算完毕，停止计算
lable11:
printf("迭代次数k=%d\n\n",k+1);
printf("全部特征值：\n");
for(i=0;i<N;i++)
printf("lamda[%d]=%1.11e+i(%1.11e)\n",i,lamda[i][0],lamda[i][1]);
printf("全部特征值计算完毕!\n\n");
//求全部实特征值对应的特征向量
for(i=0;i<N;i++)
if(lamda[i][1]==0)
guass(lamda[i][0]);
}
//对矩阵进行上三角化
void hessenberg()
{
double p[N],u[N],q[N],w[N];
double sum,h,c,d,t;
int i,j,r;
for(r=0;r<N-2;r++)
{
sum=0;
for(i=r+2;i<N;i++)
sum+=a[i][r]*a[i][r];
if(sum==0)
continue;
d=sqrt(sum+a[r+1][r]*a[r+1][r]);
if(a[r+1][r]<=0)
c=d;
else
c=-d;
h=c*(c-a[r+1][r]);
for(i=0;i<r+1;i++)
u[i]=0;
u[r+1]=a[r+1][r]-c;
for(i=r+2;i<N;i++)
u[i]=a[i][r];
for(i=0;i<N;i++)
{
p[i]=0;
for(j=r+1;j<N;j++)
p[i]+=a[j][i]*u[j];
p[i]=p[i]/h;
}
for(i=0;i<N;i++)
{
q[i]=0;
for(j=r+1;j<N;j++)
q[i]+=a[i][j]*u[j];
q[i]=q[i]/h;
}
t=0;
for(i=r+1;i<N;i++)
t+=p[i]*u[i];
t=t/h;
for(i=0;i<N;i++)
w[i]=q[i]-t*u[i];
for(i=0;i<N;i++)
for(j=0;j<N;j++)
a[i][j]=a[i][j]-w[i]*u[j]-u[i]*p[j];
for(i=0;i<N;i++)
for(j=0;j<N;j++)
{
if(fabs(a[i][j])<eps)
a[i][j]=0;
}
}
printf("拟上三角化矩阵A(n-1)：\n");
for(i=0;i<N;i++)
{
for(j=0;j<N;j++)
printf("%1.11e ",a[i][j]);
printf("\n");
}
printf("\n");
}
//求二阶子阵的两个特征根
void DK()
{
double det_D,x;
det_D=a[m-1][m-1]*a[m][m]-a[m-1][m]*a[m][m-1];
x=(a[m-1][m-1]+a[m][m])*(a[m-1][m-1]+a[m][m])-4*det_D;
if(x>0)
{
s[0][0]=((a[m-1][m-1]+a[m][m])+sqrt(x))/2;
s[0][1]=0;
s[1][0]=((a[m-1][m-1]+a[m][m])-sqrt(x))/2;
s[1][1]=0;
}
else
{
s[0][0]=(a[m-1][m-1]+a[m][m])/2;
s[0][1]=sqrt(-x)/2;
s[1][0]=(a[m-1][m-1]+a[m][m])/2;
s[1][1]=-sqrt(-x)/2;
}
}
//对A进行QR分解，得上三角形式
void two_step_move_QR()
{
double M[N][N],p[N],u[N],q[N],w[N],v[N];
double sum,h,c,d,t,s;
int i,j,r,l;
s=a[m-1][m-1]+a[m][m];
t=a[m-1][m-1]*a[m][m]-a[m][m-1]*a[m-1][m];
for(i=0;i<=m;i++)
for(j=0;j<=m;j++)
{
M[i][j]=0;
for(l=0;l<=m;l++)
M[i][j]+=a[i][l]*a[l][j];
M[i][j]-=s*a[i][j];
if(i==j)
M[i][j]+=t;
}
for(r=0;r<=m-1;r++)
{
sum=0;
for(i=r+1;i<=m;i++)
sum+=M[i][r]*M[i][r];
if(sum==0)
continue;
d=sqrt(sum+M[r][r]*M[r][r]);
if(M[r][r]<=0)
c=d;
else
c=-d;
h=c*(c-M[r][r]);
for(i=0;i<r;i++)
u[i]=0;
u[r]=M[r][r]-c;
for(i=r+1;i<=m;i++)
u[i]=M[i][r];
for(i=0;i<=m;i++)
{
v[i]=0;
for(j=r;j<=m;j++)
v[i]+=M[j][i]*u[j]/h;
}
for(i=0;i<=m;i++)
for(j=r;j<=m;j++)
M[i][j]-=u[i]*v[j];
for(i=0;i<=m;i++)
{
p[i]=0;
for(j=r;j<=m;j++)
p[i]+=a[j][i]*u[j]/h;
}
for(i=0;i<=m;i++)
{
q[i]=0;
for(j=r;j<=m;j++)
q[i]+=a[i][j]*u[j]/h;
}
t=0;
for(i=r;i<=m;i++)
t+=p[i]*u[i]/h;
for(i=0;i<=m;i++)
w[i]=q[i]-t*u[i];
for(i=0;i<=m;i++)
for(j=0;j<=m;j++)
a[i][j]=a[i][j]-w[i]*u[j]-u[i]*p[j];
for(i=0;i<N;i++)
for(j=0;j<N;j++)
{
if(fabs(a[i][j])<eps)
a[i][j]=0;
}
}
}
//用Guass法求全部实特征值对应的特征向量
void guass(double y)
{
double c[N][N],m[N][N];
double x[N],b[N];
int l,i,j;
double sum;
for(i=0;i<N;i++)
for(j=0;j<N;j++)
{
if(i==j)
c[i][j]=y-t[i][j];
else
c[i][j]=-t[i][j];
}
for(i=0;i<N;i++)
b[i]=0;
for(l=0;l<N-1;l++)
{
if(c[l][l]==0)
break;
else
for(i=l+1;i<N;i++)
{
m[i][l]=c[i][l]/c[l][l];
for(j=l+1;j<N;j++)
c[i][j]-=m[i][l]*c[l][j];
b[i]-=m[i][l]*b[l];
}
}
x[N-1]=1;
for(l=N-2;l>=0;l--)
{
sum=0;
for(j=l+1;j<N;j++)
sum+=c[l][j]*x[j];
x[l]=(b[l]-sum)/c[l][l];
}
printf("特征值lamda:%1.11e\n特征向量：\n",y);
for(i=0;i<N;i++)
printf("%1.11e\n",x[i]);
printf("\n");
}
//对拟上三角化后的矩阵A(n-1)进行QR分解
void QR()
{
double Q[N][N];
double b[N][N];
double p[N];
double u[N];
double w[N];
double sum,h,c,d;
int i,j,r;
for(i=0;i<N;i++)
{
for(j=0;j<N;j++)
Q[i][j]=0;
Q[i][i]=1;
}
for(i=0;i<N;i++)
for(j=0;j<N;j++)
b[i][j]=a[i][j];
for(r=0;r<N-1;r++)
{
sum=0;
for(i=r+1;i<N;i++)
sum=sum+b[i][r]*b[i][r];
if(sum==0)
break;
sum=0;
for(i=r;i<N;i++)
sum+=b[i][r]*b[i][r];
d=sqrt(sum);
if(b[r][r]<=0)
c=d;
else
c=-d;
h=c*(c-b[r][r]);
for(i=0;i<r;i++)
u[i]=0;
u[r]=b[r][r]-c;
for(i=r+1;i<N;i++)
u[i]=b[i][r];
for(i=0;i<N;i++)
{
sum=0;
for(j=r;j<N;j++)
sum+=b[j][i]*u[j];
p[i]=sum/h;
}
for(i=0;i<N;i++)
{
w[i]=0;
for(j=r;j<N;j++)
w[i]+=Q[i][j]*u[j];
}
for(i=0;i<N;i++)
for(j=r;j<N;j++)
Q[i][j]-=w[i]*u[j]/h;
for(i=r;i<N;i++)
for(j=0;j<N;j++)
b[i][j]-=u[i]*p[j];
for(i=r;i<N;i++)
for(j=0;j<N;j++)
{
if(fabs(b[i][j])<eps)
b[i][j]=0;
}
}
printf("QR分解后的R矩阵：\n");
for(i=0;i<N;i++)
{
for(j=0;j<N;j++)
printf("%1.11e ",b[i][j]);
printf("\n");
}
printf("\n");
printf("QR分解后的Q矩阵：\n");
for(i=0;i<N;i++)
{
for(j=0;j<N;j++)
printf("%1.11e ",Q[i][j]);
printf("\n");
}
printf("\n");
}⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -