资源描述
一、算法的设计方案：
（一）各所求值得计算方法
1、最大特征值λ501，最小特征值λ1，按模最小特征值λs的计算方法
首先使用一次幂法运算可以得到矩阵的按模最大的特征值λ，λ必为矩阵A的最大或最小特征值，先不做判断。对原矩阵A进行一次移项，即（A-λI），在进行一次幂法运算，可以得到另一个按模最大特征值λ0。比较λ和λ0的大小，较大的即为λ501，较小的即为λ1。
对矩阵A进行一次反幂法运算，即可得到按模最小特征值λs。
2、A与μk值最接近的特征值λik的计算方法
首先计算出k所对应的μk值，对原矩阵A进行一次移项，即（A-μkI），得到一个新的矩阵，对新矩阵进行一次反幂法运算，即可得到一个按模最小特征值λi。则原矩阵A与μk值最接近的特征值 λik=λi+μk。
3、A的（谱范数）条件数cond（A）2的计算方法
其中矩阵A的按模最大和按模最小特征值。
（二）程序编写思路。
由于算法要求A的零元素不存储，矩阵A本身为带状矩阵，所以本题的赋值，LU分解，反幂法运算过程中，均应采用Doolittle分解法求解带状线性方程组的算法思路。
幂法、反幂法和LU分解均是多次使用，应编写子程序进行反复调用。
二、源程序：
#include<stdio.h>
#include<iostream>
#include<stdlib.h>
#include<math.h>
#include<float.h>
#include<iomanip> /*头文件*/
/*定义全局变量*/
#define N 502 /*取N为502，可实现从1到501的存储，省去角标变换的麻烦*/
#define epsilon 1.0e-12 /*定义精度*/
#define r 2 /*r，s为带状矩阵的半带宽，本题所给矩阵二者都是2*/
#define s 2
double c[6][N]; /*定义矩阵存储压缩后的带状矩阵*/
double fuzhi(); /*赋值函数*/
void LUfenjie(); /*LU分解程序*/
int max(int a,int b); /*求两个数字中较大值*/
int min(int a,int b); /*求两个数字中较小值*/
double mifa(); /*幂法计算程序*/
double fanmifa(); /*反幂法计算程序*/
double fuzhi() /*赋值程序，按行赋值，行从1到5，列从1到501，存储空间 的第一行第一列不使用，角标可以与矩阵一一对应，方便书写程序*/
{
int i,j;
i=1;
for(j=3;j<N;j++)
{c[i][j]=-0.064;}
i=2;
for(j=2;j<N;j++)
{c[i][j]=0.16;}
i=3;
for(j=1;j<N;j++)
{c[i][j]=(1.64-0.024*j)*sin(0.2*j)-0.64*exp(0.1/j);}
i=4;
for(j=1;j<N-1;j++)
{c[i][j]=0.16;}
i=5;
for(j=1;j<N-2;j++)
{c[i][j]=-0.064;}
return(c[i][j]);
}
int max(int a,int b)
{ return((a>b)?a:b);
}
int min(int a,int b)
{ return((a<b)?a:b);
}
void LUfenjie() /*LU分解程序，采用的是带状矩阵压缩存储后的LU分解法*/
{
double temp;
int i,j,k,t;
for(k=1;k<N;k++)
{ for(j=k;j<=min(k+s,N-1);j++)
{
temp=0;
for(t=max(1,max(k-r,j-s));t<=(k-1);t++)
{temp=temp+c[k-t+s+1][t]*c[t-j+s+1][j];}
c[k-j+s+1][j]=c[k-j+s+1][j]-temp;
}
for(i=k+1;i<=min(k+r,N-1);i++)
{
temp=0;
for(t=max(1,max(i-r,k-s));t<=(k-1);t++)
{temp=temp+c[i-t+s+1][t]*c[t-k+s+1][k];}
c[i-k+s+1][k]=(c[i-k+s+1][k]-temp)/c[s+1][k];
}
}
}
double mifa() /*幂法计算程序*/
{
double u0[N],u1[N];
double temp,Lu,beta=0,beta0;
int i,j;
for(i=1;i<N;i++) /*选取迭代初始向量*/
{u0[i]=1;}
do
{
beta0=beta;
temp=0;
for(i=1;i<N;i++)
{temp=temp+u0[i]*u0[i]; }
Lu=sqrt(temp);
for(i=1;i<N;i++)
{
u1[i]=u0[i]/Lu;
}
for(i=1;i<N;i++)
{
temp=0;
for(j=max(i-1,1);j<=min(i+2,N-1);j++)
{temp=temp+c[i-j+s+1][j]*u1[j]; }
u0[i]=temp;
}
temp=0;
for(i=1;i<N;i++)
{temp=temp+u1[i]*u0[i]; }
beta=temp;
}while(fabs(beta-beta0)/fabs(beta)>=epsilon); /*迭代运行条件判断*/
return(beta);
}
double fanmifa() /*反幂法计算程序*/
{
double u0[N],u1[N],u2[N];
double temp,Lu,beta=0,beta0;
int i,t;
for(i=1;i<N;i++) /*选取迭代初始向量*/
{u0[i]=1;}
do
{
beta0=beta;
temp=0;
for(i=1;i<N;i++)
{temp=temp+u0[i]*u0[i]; }
Lu=sqrt(temp);
for(i=1;i<N;i++)
{
u1[i]=u0[i]/Lu;
u2[i]=u1[i];
}
fuzhi();
LUfenjie();
/*带状矩阵压缩存储并进行LU分解后，求解线性方程组得到迭代向量uk，即程序中的u0*/
for(i=2;i<N;i++)
{ temp=0;
for(t=max(1,i-r);t<=(i-1);t++)
{temp=temp+c[i-t+s+1][t]*u2[t];}
u2[i]=u2[i]-temp;
}
u0[N-1]=u2[N-1]/c[s+1][N-1];
for(i=N-2;i>=1;i--)
{ temp=0;
for(t=i+1;t<=min(i+s,N-1);t++)
{temp=temp+c[i-t+s+1][t]*u0[t];}
u0[i]=(u2[i]-temp)/c[s+1][i];
}
temp=0;
for(i=1;i<N;i++)
{temp=temp+u1[i]*u0[i]; }
beta=temp;
beta=1/beta; /*beta即为所求特征值，可直接返回*/
}while(fabs(beta-beta0)/fabs(beta)>=epsilon); /*迭代运行条件判断*/
return(beta);
}
void main()
{
double u[40]; /*定义数组，存放k值运算得到的μk值*/
double lambda1,lambda501,lambdak,a,b,d,cond,det;
int i,j,k;
fuzhi();
a=mifa(); /*幂法计算按模最大值*/
fuzhi();
d=fanmifa(); /*反幂法计算按模最小值*/
fuzhi();
for(j=1;j<N;j++)
{
c[3][j]=c[3][j]-a;
}
b=mifa()+a; /*移项后幂法计算按模最大值*/
if(a>b) /*比较两个按模最大值大小，并相应输出最大特征值λ501和最小特征值λ1*/
{
lambda1=b;
lambda501=a;
printf("矩阵A最小的特征值lambda1=%13.11e\n",lambda1);
printf("矩阵A最大的特征值lambda501=%13.11e\n",lambda501);
}
else
{
lambda1=a;
lambda501=b;
printf("矩阵A最小的特征值lambda1=%13.11e\n",lambda1);
printf("矩阵A最大的特征值lambda501=%13.11e\n",lambda501);
}
printf("矩阵A按模最小特征值lambdas=%13.11e\n",d); /*输出按模最小特征值λs*/
for(k=1;k<40;k++) /*对每一个进行移项反幂法运算，求出最接近μk的特征值并输出*/
{
u[k]=(lambda501-lambda1)*k/40+lambda1;
fuzhi();
for(j=1;j<N;j++)
{c[3][j]=c[3][j]-u[k];
}
lambdak=fanmifa()+u[k];
i=k;
printf("矩阵A最接近uk特征值lambdak%d=%13.11e\n",i,lambdak);
}
cond=fabs(a/d);
printf("A的条件数=%13.11e\n",cond); /*计算A条件数并输出*/
fuzhi(); /*计算A的行列式值并输出*/
LUfenjie();
det=1;
for(i=1;i<N;i++)
{
det=det*c[3][i];
}
printf("行列式的值detA=%13.11e\n",det);
}
三、程序的运行结果：
四、初始向量的选取对计算结果的影响：
（一）选取形式不变，数值变换
1、取u0为[0.5,0.5………..0.5]，运行结果如下：
2、取u0为[50,50………..50]，运行结果如下：
从运行结果来看，此类初始向量的选取对结果不会产生影响，即使选成0，结果也不变化。
（二）改变选取形式
取u0为[1,0………..0]，运行结果如下：
可以看出，改变初始向量的类型之后，运行结果出现了明显的变化，说明迭代初始向量的选取对迭代运算的结果是有影响的。为了观察初始向量与运行结果之间的联系，下面以幂法为例进行试验。
1、取u0为[1,0………..0]，运行幂法程序，运行结果为：
对特征值对应的特征向量uk进行输出，找出数据中最接近1的数据，在数组第22个位置，且第1位周围的数值与1比较接近。其他位置的数值对比1均为小量。如下图：
2、取u0为[0……1…..0]，运行幂法程序，1位于第200个位置运行结果为：
对特征值对应的特征向量uk进行输出，找出数据中最接近1的数据，在数组第198个位置，且第200位周围的数值与1比较接近。其他位置的数值对比1均为小量。如下图：
3、取u0为[0……1…..0]，运行幂法程序，1位于第300个位置运行结果为：
对特征值对应的特征向量uk进行输出，找出数据中最接近1的数据，在数组第288个位置，且第300位周围的数值与1比较接近。其他位置的数值对比1均为小量。如下图：
4、取u0为[0……1…..0]，运行幂法程序，1位于第400个位置运行结果为：
对特征值对应的特征向量uk进行输出，找出数据中最接近1的数据，在数组第413个位置，且第400位周围的数值与1比较接近。其他位置的数值对比1均为小量。如下图：
5、取u0为[0……1…..0]，运行幂法程序，1位于第500个位置运行结果为：
对特征值对应的特征向量uk进行输出，找出数据中最接近1的数据，在数组第477个位置，且第500位周围的数值与1比较接近。其他位置的数值对比1均为小量。如下图：
综上可以看出，迭代初始向量的选取对迭代结果有明显的影响。输出的最终迭代向量（即矩阵计算出特征值所对应的特征向量）与所选的迭代初始向量形式上很接近。因此可以总结如下：
当所选取的迭代初始向量与矩阵的特征向量uk最为接近时，那么，幂法迭代计算结果输出的应该就是特征向量uk对应所对应的特征值。
以上试验只是考虑了初始向量u0为[0……1…..0]的情况，如果已知一个矩阵的特征向量，那么以其中任一特征向量uk为迭代初始向量，进行幂法运算，输出结果必为特征向量uk对应的特征值λk,而不一定是按模最大特征值。例如：矩阵 A= 的特征值分别为-1,2,5，三者对应的特征向量分别为（-2，-2，-1）T，（-2,1,2）T，（1，-2,2）T，如果以（-2，-2，-1）T为初始迭代向量进行幂法运算，得到的按模最大特征值为-1，而不是5.
同理，在反幂法运算中也可以得到相同的结论。
展开阅读全文
浙ICP备2021020529号-1 | 浙B2-20240490 
