lu.cpp

LU分解算法 及 分解后矩阵解线性方程组

/*06401 060995 耿晓阳*/
/*LU原理： AX=B A=LU =>  LUX=B => L(UX)=B UX=Y  =>  LY=B
           由于L为下三角矩阵，对应的方程组的解可由迭代求出*/
/*对应关系   a=>A    q=>X      b=>B
           LU分解后   a=>上L下U（U的1省略）   q=>Y*/ 

#include<iostream>
#include<cstdlib>
#include<cmath>
#define N 3               //矩阵的维数，为了简便 先定成3 
using namespace std; 
bool LU(double* a,double *q,int n)   //a[i][j]=&(a+i*n+j)
{ 
    int i,j,k;
    for(i=0;i<n;i++)
         if(a[i*n+i]==0.0)
              return false;           //对角线元素若为零，则无法进行LU分解 
    for(k=0;k<n;k++)
    {
        for(j=k+1;j<n;j++)
            a[k*n+j]=a[k*n+j]/a[k*n+k];
        for(i=k+1;i<n;i++)
           for(j=k+1;j<n;j++)
              a[i*n+j]=a[i*n+j]-a[i*n+k]*a[k*n+j];
    }
    return true;                  //有两种情况会false,  1.不是n*n矩阵    2.对应行列式的值为0.....但是我没做 
} 

bool solve(double* b, double const* lu, double* q, int n) 
{
     int i,j,t;
     double m[N],k[N];
     for(t=0;t<n;t++)
     {
          m[t]=0.0;
          k[t]=0.0;
     }                                //初始化中间变量数组  
     q[0]=b[0]/lu[0];
     for(i=1;i<n;i++)
     {
        for(j=0;j<=i-1;j++)
           m[i]=m[i]+lu[i*n+j]*q[j];
        q[i]=(b[i]-m[i])/lu[i*n+i];
     }                                //求Y矩阵的值 
     b[n-1]=q[n-1];
     for(i=n-2;i>=0;i--)
     {
         for(j=i+1;j<=n-1;j++)
             k[i]=k[i]+lu[i*n+j]*b[j];
         b[i]=q[i]-k[i];
     }                               //求X矩阵的值 
     return true;
}
int main()
{
    int i,j;
    int n; 
    n=N;
    //double a[9]={1.00,2.34,2.1,4.50,6.00,8.66,76.34,7.7,0.98};   //先写个3*3的吧，太多编的太累 
    double a[9]={1,1,1,1,2,1,1,1,2};    //为了验证用的能手算的矩阵
    //double a[9]={1,1,1,1,1,1,1,1,1};    //错误的矩阵，无结果 
    double q[3]={1.0,2.0,3.0};                                   //和a矩阵维数一致 ...要改维数，这也得改 
    //double b[3]={3.3,6.7,4.00};                                //同上 
    double b[3]={2.00,3.00,3.00};
    if(LU(a,q,n)==true)                 //调用LU分解函数
    { 
      cout<<"LU分解之后的矩阵:\n";
      for(i=0;i<n*n;i++)
      {
         cout<<a[i]<<"  ";
         if(i%n==2)
            cout<<"\n";      //按行列打印矩阵 
      }
    }
    cout<<"\n\n\n";
    
    if(solve(b,a,q,n)==true)
    {
       cout<<"线性方程组的结果为：\n"; 
       for(i=0;i<n;i++)
          cout<<"  "<<"X["<<i<<"]="<<b[i]<<'\n';    //打印结果 
    }
    
    cout<<endl;
    system("pause");
    return 0;
}