lu.txt

通过LU分解解线形方程组

#include "iostream"
#include "math.h"
#define MAX_SIZE 100
using namespace std;
void findans(double*Xn,double L[][MAX_SIZE],double U[][MAX_SIZE],double* b,int n)
{
	double *Yn=new double[n];
	int i,j;
	double sum = 0;
	Yn[0]=b[0];
	for(i=2;i<=n;i++)
	{
		sum=0;
		for(j=1;j<=i-1;j++)

			sum=sum+L[i-1][j-1]*Yn[j-1];
		Yn[i-1]=b[i-1]-sum;
	}
	Xn[n-1]=Yn[n-1]/U[n-1][n-1];
	for(i=n-1;i>0;i--)
	{
		sum=0;
		for(j=i+1;j<=n;j++)
			sum=sum+U[i-1][j-1]*Xn[j-1];
		Xn[i-1]=(Yn[i-1]-sum)/U[i-1][i-1];
	}
}
void findLU(double ark[][MAX_SIZE],double L[][MAX_SIZE],double U[][MAX_SIZE],int n)
{
	//ark是方阵，L U 将要存的是上三角和下三角
	//n是方阵的阶数
	int k,j,i;
	double sum;

	for(k=1;k<=n;k++)
	{
		for(j=k;j<=n;j++)
		{
			sum=0;
			for(i=1;i<k;i++)
			{
				sum=sum+L[k-1][i-1]*U[i-1][j-1];
				if(L[k - 1][i - 1]*U[i-1][j-1] != 0 )//我加的
				{	if( L[k - 1][i - 1] == 1)
				printf("U[%d][%d]: U[%d][%d]\n",k-1,j-1,i,j-1);
				else if(U[i-1][j-1] == 1)
					printf("U[%d][%d]:L[%d][%d]\n",k-1,j-1,k-1,i-1);
				else
					printf("U[%d][%d]:L[%d][%d]*U[%d][%d]\n",k-1,j-1,k-1,i-1,i,j-1);
				}

			}

			U[k-1][j-1]=ark[k-1][j-1]-sum;
		}

		for(i=k+1;i<=n;i++)
		{    
			sum=0;
			for(j=1;j<k;j++)
			{
				sum=sum+L[i-1][j-1]*U[j-1][k-1];
				if( L[i-1][j-1]*U[j-1][k-1] != 0)//我加的
				{
					if( L[i-1][j-1] == 1 )
						printf("L[%d][%d]:U[%d][%d]\n",i-1,k-1,j-1,k-1);
					else if(U[j-1][k-1]== 1)
						printf("L[%d][%d]:L[%d][%d]\n",i-1,k-1,i-1,j-1);
					else 
						printf("L[%d][%d]:L[%d][%d]*U[%d][%d]\n",i-1,k-1,i-1,j-1,j-1,k-1);
				}
			}
			L[i-1][k-1]=(ark[i-1][k-1]-sum)/U[k-1][k-1];

		}
	}
}

void show(double*Xn,double L[][MAX_SIZE],double U[][MAX_SIZE],int n)
{
	int i,j;

	cout<<"======================================="<<endl;
	cout<<"///////////下面为输出LU///////////////"<<endl;
	cout<<"======================================="<<endl;
	cout<<"L的行列式如下："<<endl;
	for(i=0;i<n;i++)
	{
		for(j=0;j<n;j++)
		{
			if(j<i)
				printf("%-.9f\t",L[i][j]);
			if(j>i)
				cout<<"0\t\t"<<" ";
			if(j==i)
				cout<<"1\t\t"<<" " ;
		}
		cout<<endl;
	}
	cout<<endl<<endl;
	cout<<"U的行列式如下："<<endl;
	for(i=0;i<n;i++)
	{
		for(j=0;j<n;j++)
		{
			if(j>=i)
				printf("%.11f ",U[i][j]);
			else
				cout<<"0.00000000000"<<" ";			 
		}
		cout<<endl;
	}
	cout<<endl<<endl;
	cout<<"======================================"<<endl;
	cout<<"//////////下面为输出结果/////////////"<<endl;
	cout<<"======================================"<<endl;
	cout<<endl;
	cout<<"线性方程组的解如下：\n";
	for(int column=0;column<n;column++)
		cout<<"X"<<column<<"="<<Xn[column]<<" ,";
	cout<<endl<<endl;
}

int main(int argc, char* argv[])
{
	int n=0;
	char readChar=' ';
	int row=0;
	int column=0;
	int i=0;
	int j=0;
	int k=0;

	//freopen("newLUinput.txt","r",stdin);
	//freopen("newLUoutput.txt","w",stdout);

	cout<<"请输入方阵的阶数(0<n<101):\nn=";
	cin>>n;
	/*double **L=new double*[n];//动态分配空间
	double **U=new double*[n];
	double **ark=new double*[n]; //建立方阵
	for(int r=0;r<n;r++)
	{
	ark[r]=new double[n];
	L[r]=new double[n];
	U[r]=new double[n];
	}  
	double *b=new double[n];
	double *Xn=new double[n];*/
	double L[MAX_SIZE][MAX_SIZE],U[MAX_SIZE][MAX_SIZE],ark[MAX_SIZE][MAX_SIZE],Xn[MAX_SIZE],b[MAX_SIZE];


	//cout<<"输入方正A的元素，

	for(int i =0;i<n;i++)
		for(int j =0;j<n;j++)
			cin>>ark[i][j];
       for(int i=0;i<n;i++)
		   cin>>b[i];
	
	findLU(ark,L,U,n);
	findans(Xn,L,U,b,n);
	show(Xn,L,U,n);
	return 0;
}