lu算法.cpp

<计算方法>ＬＵ算法过程及程序

/*******************************************************************
*程序名称：LU算法 
*学生姓名：吴龙武  学号：031150003
*完成时间：05年12月
*备注：平台/工具：WIN98+VC6.0
*********************************************************************/

#include<iostream.h>
#define max_n      20

void main()
{
	int i, j, k, n, m;
	double sum1=0, sum2=0;
	double a[max_n][max_n], L[max_n][max_n], U[max_n][max_n], 
			b[max_n], x[max_n], y[max_n];
	cout<< "未知数的个数：";
	cin>>n;
	cout<<"请输入a[i][j]和b[i]:\n";
	for (i=0; i<n; i++)
	{
		for (j=0; j<n; j++)
		{
			cout<<"请输入a["<<i<<"]["<<j<<"]:  ";
			cin>>a[i][j];
		}
	}
	for (i=0; i<n; i++)
	{
		cout<<"请输入b["<<i<<"]:  ";
		cin>>b[i];
	}
	
	//L,U矩阵
	for(i=0; i<n; i++)
	{	L[i][i] = 1;
	   for(j=0; j<n;j++)
	   {	if(j>i)
				L[i][j]=0;	   
	        if(j<i)
			    U[i][j] =0;
	   }

	}
	for (j=0; j<n; j++)
	{	
		U[0][j] = a[0][j];	         //U的第一行等于a的第一行
	}

	for (i=0; i<n; i++)
	{
		L[i][0] = a[i][0] / U[0][0];       //L的第一列的值
	}
	for(k=2; k<=n; k++)
	{
		for(j=k; j<=n; j++)
		{
			for(m=1; m<=k-1; m++)
		    	 sum1 = L[k-1][m-1] * U[m-1][j-1] +sum1 ;
			U[k-1][j-1] = a[k-1][j-1] - sum1; 	
		}                     //U矩阵已全部知道
	}
	for(k=2; k<=n; k++)
	{	for(i=k+1; i<=n; i++)
		{	
			for(m=1; m<=k-1; m++)
			 	sum2 = L[i-1][m-1] * U[m-1][k-1] +sum2 ;
			L[i-1][k-1] = (a[i-1][k-1]- sum2) / U[k-1][k-1]; 		
		}                                  //L矩阵已全部知道

	}                    /*至此，L,U两个矩阵已知道*/
	cout<<endl;
	cout<<"a矩阵：\n";
	for(i=0; i<n; i++)
	{
		for (j=0; j<n; j++)
		{	cout <<   a[i][j]<<"        ";}
			cout<<endl;
		}
		cout<<endl;
		cout<<"L矩阵：\n";
		for(i=0; i<n; i++)
		{	for (j=0; j<n; j++)
			{cout<<L[i][j]<<"       ";}
			cout<<endl;
		}
		cout<<endl;
		cout<<"U矩阵：\n";
		for(i=0; i<n; i++)
		{	for (j=0; j<n; j++)
			{cout<<U[i][j]<<"       ";}
			cout<<endl;
		}


	/*求LY=b*/
	double sum3=0;
	for (i=0; i<n; i++)
	{
		if (i != 0)
		{
			for (j=i-1; j>=0; j--)
				sum3=L[i][j] * y[j];		
		}
		y[i] = b[i] - sum3;
	}                   /*现已知道y[i]*/

	/*求UX=Y*/
	double sum4=0; 
	for (i=n-1; i>=0; i--)
	{
		if(i != n-1)
		{
			for (j=n-1; j>i; j--)
				sum4= U[i][j] * x[j];
		}
		x[i] = (y[i] - sum4)/U[i][i];
		
		cout << "方程组的解是："<<endl;
		cout<<"x["<<i<<"]= "<<x[i]<<endl;
	}
}

/*
未知数的个数：2
请输入a[i][j]和b[i]:
请输入a[0][0]:  2
请输入a[0][1]:  3
请输入a[1][0]:  7
请输入a[1][1]:  1
请输入b[0]:  5
请输入b[1]:  8

a矩阵：
2        3
7        1

L矩阵：
1       0
3.5       1

U矩阵：
2       3
0       -9.5
方程的解是：
x[1]= 1
方程的解是：
x[0]= 1
Press any key to continue  
*/