shuzhifangfa.cpp

三角分解法中的lu分解法

#include <iostream>
#include <math.h>
using namespace std;
void daolite(double a[3][3])
{ //由a[][]分解为L[][](单位下三角矩阵),U[][](上三角矩阵)
	double u[3][3],double y[3],double x[3];
	double l[3][3]; 
	double b[3]={1.0,3.0,4.0};
    int i,r,j,k;
	for(r=0;r<3;r++)
		for(j=0;j<3;j++)
		{
			u[r][j]=0;
		}
		for(r=0;r<3;r++)
			for(j=0;j<3;j++)
			{
			if(r==j) l[r][j]=1;
			else l[r][j]=0;
			}
	for( r=0;r<3;r++)
	{
		for( j=r;j<3;j++)
		{
			double s=0.0;
			for( k=0;k<r;k++)
            s+=l[r][k]*u[k][j];//计算：∑l[r][k]*u[k][j]
            u[r][j]=a[r][j]-s;//计算：u[r][j]=a[r][j]-∑l[r][k]*u[k][j]
		} 
		for( i=r+1;i<3;i++)
		{
            double s=0.0;
            for( k=0;k<r;k++)
            s+=l[i][k]*u[k][r];//计算：∑l[i][k]*u[k][r]
            l[i][r]=(a[i][r]-s)/u[r][r];//计算：l[i][r]=(a[i][r]-∑l[i][k]*u[k][r])/u[r][r]
		}
	}
	cout<<"其中u矩阵为："<<endl;//输出u矩阵
	 for( r=0;r<3;r++)
		{
			for( j=0;j<3;j++)
			{
		cout<<u[r][j]<<"\t";
			}
                cout<<endl;
		} 
	cout<<"其中l矩阵为："<<endl;//输出l矩阵

	 for( i=0;i<3;i++)
		{
			for( r=0;r<3;r++)
			{
		cout<<l[i][r]<<"\t";
			}
                cout<<endl;
		} 
	 
	    for( i=0;i<3;i++)
		{double s=0.0;
			for(k=0;k<i;k++)
			{
			  s=s+l[i][k]*y[k];
			}
			y[i]=b[i]-s;
		}
	  
     for(i=2;i>=0;i--) 
	 {  double sum=0.0;

          for(k=i+1;k<2;k++)
		  {
                sum+=u[i][k]*x[k]; 
		  }
           x[i]=(y[i]-sum)/u[i][i]; 
	 } 
     cout<<"方程式组的解为：x=(";//输出x[]
     for (i=0;i<3;i++)
	 {
          cout<<x[i];
          if(i!=2)cout<<",";
	 }
     cout<<")"<<endl;


}
int main()
{ 
double a[3][3];
int i,j;
cout<<"请输入矩阵:"<<endl;
		for(i=0;i<3;i++)
			for(j=0;j<3;j++)
			cin>>a[i][j];		
	   cout<<"你输入的矩阵为:"<<endl;
               
		for( i=0;i<3;i++)
		{
			for(j=0;j<3;j++)
			{
		cout<<a[i][j]<<"\t";
			}
                cout<<endl;
		}  
    daolite(a);
	return 0;
}