Tauruser

3.2 Gauss消去法

3.2.1　顺序消去法

　　Gauss消去法就是将方程组(3.1.1)通过(n-1)步消元，将(3.1.1)转化为上三角方程组
　　　　　　　(3.2.1)
再回代求此方程组的解.
　　下面记增广矩阵，即
　　　
　　第1步　设，计算l，记为，若用乘第一行加到第i行，可消去,用Gauss变换矩阵表示
　　　　
令　　　
其中
　　一般地，假定已完成了(k-1)步消元，即已将转化为以下形式：
　　　
　　第k步，假定，计算
　　　　　　　(3.2.2)
记，，则
　　　
其中
　　　　　(3.2.3).
当k=1,2，…,n-1则可得到，即方程组(3.2.1).

　　直接回代解(3.2.1)得，
　　　　　(3.2.4)
并且有，由以上顺序消去过程可得如下定理.
　　定理2.1　设非奇异，则通过两行互换总可使，k=1,2，…,n-1.可将方程组(3.1.1)转化为(3.2.1)并求得方程组(3.1.1)的解为(3.2.4)，且有.
　　如果不做行交换，则使的条件如下.

　　定理2.2　非奇异，且各阶顺序主子式， 则，k=1,2，…,n-1.
　　证明　用归纳法，当，故.现假设(k-1)成立，即，对i=1，2，…,k-1已推出，故Gauss消去法能进行(k-1)步消元，A已约化为，即
　　　

　　　　　

故对k=1,2,…,n均成立，证毕.
　　在整个消去法消元过程中，k从1到(n-1)共需乘除法运算次数为
　　　　　
加减法次数为
　　　　　
回代过程中由公式(3.2.4)可知乘除法次数为，加减法次数为，于是Gauss消去法的乘除法总次数为，加减法次数为
　　例3.4　用Gauss消去法解方程组
　　　　　　
并求detA.
　　解　消元得
　　　
再由(3.2.4)回代，得解

讲解：

　　Gauss 消去法是将方程组AX＝b,通过消元转化为上三角方程组（3，2，1）求解，消元第一步做完后有
　　　　　　　　　
　　用矩阵表示
　　第K－1步完成后得到
　　当，可做K步，得到
　　得到，对应的方程组就是（3.2.1），利用公式（3.2.4）就可求得解。
　　定理2.2给出了进行顺序消去法的条件，即A的所有顺序生子式，而方程（3.1.1）解存在唯一的条件是

好了，原理讲完了，贴我的例程。

#include  < iostream >
#include  < vector >
#include  < cmath >
using   namespace  std;
class  CGAUSSSOLVEEQU
{
private :
    vector < vector < double >>  m_equset;
    vector < double >  m_answer;
    unsigned  int  m_n;
public :
     void  inputEquSet();
     void  solveEquSet();
     void  outputAnswer();
} ;
void  CGAUSSSOLVEEQU::inputEquSet()
{
     double  dtemp;
    vector < double >  vtemp;
    cout << " 请输入你的方程个数： " ;
    cin >> m_n;
    cout << " 请按照向量的形式输入各变量的系数。最后一位为b。每个方程一行: " << endl;
     for (unsigned  int  i( 0 );i < m_n; ++ i)
     {
        m_equset.push_back(vtemp);
         for (unsigned  int  j( 0 );j <= m_n; ++ j)
         {    
            cin >> dtemp;
            m_equset[i].push_back(dtemp);
        }
         if (m_equset[i].size() != m_n + 1 )
         {
            cout << " 输入有误，请重新输入上一个方程。 " << endl;
             -- i;
        }
    }
    
}

void  CGAUSSSOLVEEQU::solveEquSet()
{
    vector < vector < double >> ::iterator iter;
    iter = m_equset.begin();
     for (unsigned  int  m( 0 );m < m_n - 1 ; ++ m)
     {
         // 将绝对值最大的主元素移上去。此举是为了减少误差
         for (vector < vector < double >> ::iterator iter2 = iter + 1 ;iter2 != m_equset.end(); ++ iter2)
         {
             if (fabsl(iter -> front()) < fabsl(iter2 -> front()))
             {
                swap( * iter, * iter2);
            }
        }
         // 进行消元
         for (unsigned  int  i = m + 1 ;i < m_n; ++ i)
         {
             double  dm;
            dm = m_equset[i][m] / m_equset[m][m];
             for (unsigned  int  j = m;j < m_n + 1 ; ++ j)
             {
                m_equset[i][j] -= dm * m_equset[m][j];
            }
        }
         ++ iter;
    }
     // 初始化m_answer向量
     for (unsigned  int  i( 0 );i < m_n; ++ i) m_answer.push_back( 0 );
     // 求解答案
    m_answer[m_n - 1 ] = m_equset[m_n - 1 ][m_n] / m_equset[m_n - 1 ][m_n - 1 ];

     for ( int  i = m_n - 2 ;i >= 0 ; -- i)
     {
        m_answer[i] = m_equset[i][m_n];
         for ( int  j = m_n - 1 ;j > i; -- j)
            m_answer[i] -= m_answer[j] * m_equset[i][j];
        m_answer[i] /= m_equset[i][i];
    }

    
}

void  CGAUSSSOLVEEQU::outputAnswer()
{
     for (unsigned  int  i( 1 );i <= m_n; ++ i)
     {
        cout << " x( " << i << " )= " << m_answer[i - 1 ] << endl;
    }
}

int  main()
{
    CGAUSSSOLVEEQU myEqu;
    myEqu.inputEquSet();
    myEqu.solveEquSet();
    myEqu.outputAnswer();
     return   0 ;
}

// Power By Tauruser 2006.6.4