数值计算---求希尔伯特矩阵的条件数

最新推荐文章于 2023-11-25 16:43:15 发布

春天来啦_花开啦

最新推荐文章于 2023-11-25 16:43:15 发布

阅读量1.3w

点赞数 14

CC 4.0 BY-SA版权

分类专栏：数值计算文章标签： c语言编程

本文链接：https://blog.youkuaiyun.com/zhangminflying/article/details/23399201

数值计算专栏收录该内容

1 篇文章

订阅专栏

本文详细介绍了如何使用幂法和反幂法计算希尔伯特矩阵的条件数，并通过编程实验证明了条件数随矩阵维数的增加而显著增长的趋势。研究发现，随着矩阵维数的增加，希尔伯特矩阵的条件数急剧增加，揭示了数值稳定性问题的重要性。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

这几天数值计算老师交给我们一个课程设计，计算希尔伯特矩阵的条件数，观察其随维数的变化情况。

下面是程序，主要用到幂法和反幂法。

#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