lu.cpp

数值与符号计算LU分解法

#include <iostream.h>
#include <math.h>

bool Luf_cp(double *a,int *cpp,int n,double eps=1.0e-8);      //作LU分解
void solve_lu_cpp(double *b,double *lu,int *cpp,int n,double eps=1.0e-8); //求AX=b的解

int main(int argc, char* argv[])
{
	//变量的数目
	int n; 
	//系数矩阵
	double *a; 
	//当前行主元变量序号
	int *cpp;
	//常数列
	double *b;  
	int i;

	cout<<"请输入变量值"<<endl;
	cin>>n;
	if(n<1)  //非法变量数
	{
		cout<<"请输入不小于1的数";
		return 0;
	}
	
	//分配空间
	a=new double[n*n];
	b=new double[n];
	cpp=new int[n];
	
	//读入系数矩阵
	cout<<"请输入系数矩阵"<<endl;
	for(i=0;i<n*n;i++)
		cin>>a[i];

	//读入常数列
	cout<<"请输入常量列"<<endl;
	for(i=0;i<n;i++)
		cin>>b[i];
	
	//初始化变量位置向量
	for(i=0;i<n;i++)
		cpp[i]=i;
	
	/*如果LU分解成功
	 *求解Ax=b*/
	if(Luf_cp(a,cpp,n))   
	{
		solve_lu_cpp(b,a,cpp,n); 

		cout<<"解答如下:"<<endl;
		for(i=0;i<n;i++)
			cout<<b[i]<<' ';
		cout<<endl;
	}
	else
		cout<<"sorry，不能用高斯方法解决!"<<endl;

	//输出LU矩阵
	cout<<endl<<"LU矩阵如下:"<<endl;
	for(i=0;i<n*n;i++)
	{
		if(i%n==0)
			cout<<endl;
		cout<<a[i]<<'\t';
	}
	
	cout<<endl;
	//释放存储空间
	delete []a;
	delete []b;
	delete []cpp;

	return 0;
}

/*LU分解，返回职位是否成功，
 *若矩阵非奇异返回true，否则false
 *a输入为系数矩阵，返回时为LU矩阵（左下角为响应的变换向量）
 *cpp为变换后各变量位置
 *n为变量数，eps为精度要求
 *绝对值小于eps的数近似为零*/
bool Luf_cp(double *a,int *cpp,int n,double eps)      
{
	int row,colume;

	//初始化变量位置向量
	for(row=0;row<n;row++)
		cpp[row]=row;

	for(colume=0;colume<n-1;colume++)
	{
		int maxrow=colume;

		//选择当前列的最大值
		for(row=colume+1;row<n;row++)
			if(fabs(a[row*n+colume])>fabs(a[maxrow*n+colume]))
				maxrow=row;
		
		//如果整列为0，返回失败
		if(a[maxrow*n+colume]==0)
			return false;

		//如果需要，交换最大值所在行与当前行
		if(maxrow!=colume)
		{
			double temp;
			int temp2;

			temp2=cpp[maxrow];
			cpp[maxrow]=cpp[colume];
			cpp[colume]=temp2;

			for(int col=0;col<n;col++) //交换
			{
				
				temp=a[maxrow*n+col];
				a[maxrow*n+col]=a[colume*n+col];
				a[colume*n+col]=temp;
			}
		}

		//变换一下各行
		for(row=colume+1;row<n;row++)
		{
			double percent;
			
			percent=a[row*n+colume]/a[colume*n+colume];
			a[row*n+colume]=percent;
			for(int col=colume+1;col<n;col++)
			{
				a[row*n+col]-=percent*a[colume*n+col];
				if(fabs(a[row*n+col])<eps)
					a[row*n+col]=0.0;
			}
		}
	}
	return true;
}

/*求解Ax=b
 *a输入为系数矩阵，返回时为LU矩阵（左下角为响应的变换向量）
 *cpp为变换后各变量位置
 *n为变量数，eps为精度要求
 *绝对值小于eps的数近似为零*/
void solve_lu_cpp(double *b,double *lu,int *cpp,int n,double eps) 
{
	int row,colume;
	double *tb=new double[n];

	//重排常数列,tb为重排后数组
	for(row=0;row<n;row++)
		tb[row]=b[cpp[row]];

	//计算常数列
	for(row=1;row<n;row++)
	for(colume=0;colume<row;colume++)
		tb[row]-=lu[row*n+colume]*tb[colume];

	//从第n行到第1行回溯求解
	for(row=n-1;row>=0;row--)
	{
		for(colume=row+1;colume<n;colume++)
			tb[row]-=lu[row*n+colume]*tb[colume];

		tb[row]/=lu[row*n+row];
		if(fabs(tb[row])<eps)
			tb[row]=0;
	}

	//按x1--xn顺序重排解向量,b中是结果
	for(row=0;row<n;row++)
		b[row]=tb[row];

	delete []tb;
}