这几天数值计算老师交给我们一个课程设计,计算希尔伯特矩阵的条件数,观察其随维数的变化情况。
下面是程序,主要用到幂法和反幂法。
#include <iostream> #include <cmath> using namespace std; #define MAX_N 100 // 最大迭代次数 //函数声明 void Hn(double *h,int n); //希尔伯特矩阵 double getLambda_1(double *h,double e,int n); //求最大特征值 double getLambda_2(double **L,double **U,double e,int n); //求最小特征值 void LU(double **L ,double **U,double *A,int n); //LU分解法 double getCond(int n); //求对阵矩阵条件数 //主函数 int main(){ for(int i=0;i<=20;i++) { cout<<"维数 n= "<<i<<"条件数 cond="<<getCond(i)<<"\n"; } return 0; } /************************************************************************/ /* 求对阵矩阵条件数 /* int n 矩阵维数 /************************************************************************/ double getCond(int n) { double cond=0; // 条件数 double *hn=new double[n*n]; Hn(hn,n); //初始化希尔伯特矩阵 double lamda1=getLambda_1(hn,0.0005,n); //最大特征值 //cout<<"lamda1:"<<lamda1<<"\n"; double **L ; //L矩阵 double **U; //U矩阵 //该矩阵采用下三角 存储 节约内存 L=new double *[n]; for(int k=1;k<=n;k++) { L[k-1]=new double[k]; } //该矩阵采用上三角 存储 节约内存 U=new double*[n]; for(int k=n;k>0;k--) { U[n-k]=new double[n]; } LU(L,U,hn,n); //LU分解法求LU矩阵 double lamda2=getLambda_2(L,U,0.0005,n); //最小特征值 //cout<<"lamda2:"<<lamda2<<"\n"; cond=fabs(lamda1)/fabs(lamda2); //计算条件数 //cout<<"维数为:"<<n<<" 条件数为: "<<cond<<endl; //释放矩阵内存空间 for(int i=0;i<n;i++) { delete []U[i]; delete []L[i]; } delete []L; delete []U; delete []hn; //释放内存空间 return cond; //返回条件数 } /************************************************************************/ /* 希尔伯特矩阵 /*double *h 存放矩阵指针 /*int n 维数 /************************************************************************/ void Hn(double *h,int n) { for(int i=0;i<n;i++) for(int j=0;j<n;j++) { h[j+n*i]=1/((i+j+1)*1.0); } } /************************************************************************/ /* 得到按模最大特征值 幂法 /* double *h 希尔伯特矩阵 /* double e 容许误差限 /* int n 维数 /************************************************************************/ double getLambda_1(double *h,double e,int n) { double *v0=new double[n]; //初始向量v0 for(int i=0;i<n;i++) if(i==n-1) v0[i]=1; else v0[i]=0; double m0=100, m=0; int a=0; //标识迭代次数 while (fabs(m-m0)>=e&&(++a<MAX_N)) { m0=m; double *u=new double[n]; for(int i=1;i<=n;i++) { u[i-1]=0; for (int j=1;j<=n;j++) { u[i-1]=u[i-1]+h[j-1+n*(i-1)]*v0[j-1]; } } m=u[0]; int l=1; for(int i=2;i<=n;i++) { if(fabs(m)<fabs(u[i])) { m=u[i-1]; l=i; } if(v0[l-1]<0) m=-m; for(int i=1;i<=n;i++) { v0[i-1]=u[i-1]/m; } } delete []u; } delete []v0; return m; } /************************************************************************/ /* LU分解法 /* double *L 分解后的矩阵L /* double *U 分解后的矩阵U /* double *A 原矩阵 /* int n 矩阵维数 /************************************************************************/ void LU(double **L ,double **U,double *A,int n) { for(int k=1;k<=n;k++) { L[k-1][k-1]=1; for(int j=k;j<=n;j++) { double sum1=0; for(int s=1;s<=k-1;s++) { sum1+=L[k-1][s-1]*U[s-1][j-1]; } U[k-1][j-1]=A[j-1+n*(k-1)]-sum1; } if(k!=n) { for(int i=k+1;i<=n;i++) { double sum2=0; for(int s=1;s<=k-1;s++) { sum2+=L[i-1][s-1]*U[s-1][k-1]; } L[i-1][k-1]=(A[k-1+n*(i-1)]-sum2)/U[k-1][k-1]; } } } } /************************************************************************/ /* 求解按模最小特征值 反幂法 先对矩阵A 进行LU分解 再解方程组 Lw=v 和Uu=w 求出u /* double **L L矩阵 /* double **U U矩阵 /* double e 容许误差限 /* int n 向量维数 /************************************************************************/ double getLambda_2(double **L,double **U,double e,int n) { double *v=new double[n]; for(int i=0;i<n;i++) if(i==n-1) v[i]=1; else v[i]=0; double m0=5,m=0; int a=0 ; //标识迭代次数 while (fabs(m-m0)>=e&&(++a)<=MAX_N) { m0=m; double *w=new double[n]; double *u=new double[n]; //LU分解法 for (int k=1;k<=n;k++) { double sum1=0; for(int s=1;s<=k-1;s++) { sum1+=L[k-1][s-1]*w[s-1]; } w[k-1]=v[k-1]-sum1; //解 Uy=b; } for(int k=n;k>=1;k--) { double sum2=0; for(int s=k+1;s<=n;s++) { sum2+=U[k-1][s-1]*u[s-1]; } u[k-1]=(w[k-1]-sum2)/U[k-1][k-1]; //解Lx=y } m=u[0]; int l=1; for(int i=2;i<=n;i++) { if(fabs(m)<fabs(u[i-1])) { m=u[i-1],l=i; } if(v[l-1]<0) { m=-m; } } for(int i=1;i<=n;i++) { v[i-1]=u[i-1]/m; } delete u; delete w; //cout<<" m= "<<m; } delete v; return 1/m; }
以下是输出结果:维数 n= 0条件数 cond=0
维数 n= 1条件数 cond=1
维数 n= 2条件数 cond=19.2814
维数 n= 3条件数 cond=524.06
维数 n= 4条件数 cond=15513.6
维数 n= 5条件数 cond=476611
维数 n= 6条件数 cond=1.49513e+007
维数 n= 7条件数 cond=4.75381e+008
维数 n= 8条件数 cond=1.52582e+010
维数 n= 9条件数 cond=4.9314e+011
维数 n= 10条件数 cond=1.6024e+013
维数 n= 11条件数 cond=5.21367e+014
维数 n= 12条件数 cond=1.6469e+016
维数 n= 13条件数 cond=5.34184e+017
维数 n= 14条件数 cond=5.44566e+017
维数 n= 15条件数 cond=4.03223e+017
维数 n= 16条件数 cond=5.06755e+017
维数 n= 17条件数 cond=5.13573e+017
维数 n= 18条件数 cond=3.72509e+017
维数 n= 19条件数 cond=5.5365e+017
维数 n= 20条件数 cond=2.47509e+017