數值計算---求希爾伯特矩陣的條件數
阿新 • • 發佈：2019-02-08
這幾天數值計算老師交給我們一個課程設計，計算希爾伯特矩陣的條件數，觀察其隨維數的變化情況。
下面是程式，主要用到冪法和反冪法。
#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