hilbert.cpp

来自「验证hilbert矩阵是病态矩阵」· C++ 代码 · 共 54 行

#include<iostream.h>

const int n=20;

void main()
{
	int i,j,k;
	double a[n][n],b[n],deta;
	for(i=0;i<n;i++){
		b[i]=0.0;
		for(j=0;j<n;j++){
			a[i][j]=1.0/(i+1+j);
			b[i]=b[i]+a[i][j];
		}
		b[i]=-b[i];
	}
	deta=1.0;
	for(k=0;k<n;k++){
		deta=a[k][k]*deta;
		for(j=0;j<n;j++){
			if(j!=k)
				a[k][j]=-a[k][j]/a[k][k];
		}
		b[k]=-b[k]/a[k][k];
		a[k][k]=1.0/a[k][k];
		
		for(i=0;i<n;i++){
			if(i!=k){
				for(j=0;j<n;j++){
					if(j!=k)
						a[i][j]=a[i][j]+a[i][k]*a[k][j];
				}
				b[i]=b[i]+a[i][k]*b[k];
				a[i][k]=a[i][k]*a[k][k];
			}
		}
	}
	cout<<"H-1="<<endl;
	for(i=0;i<n;i++){
		for(j=0;j<n;j++){
			cout<<a[i][j]<<"\t";
		}
		cout<<endl;
	}
	cout<<"x="<<endl;
	for(i=0;i<n;i++){
		cout<<b[i]<<"\t";
	}
	cout<<endl;
	cout<<"detH= "<<deta<<endl;
}