<8/7>集训日记

今天看了有关高斯消元的一些内容。

高斯消元法是线性代数中的一个算法，可用来求解线性方程组、求出矩阵的秩、求出可逆方阵的逆矩阵。

 

下面是一篇博客中讲的高斯消元法在编程中的应用。

首先，先介绍程序中高斯消元法的步骤：

 (我们设方程组中方程的个数为equ，变元的个数为var，注意：一般情况下是n个方程，n个变元，但是有些题目就故意让方程数与变元数不同)

1. 把方程组转换成增广矩阵。

2. 利用初等行变换来把增广矩阵转换成行阶梯阵。

 枚举k从0到equ– 1，当前处理的列为col(初始为0) ，每次找第k行以下(包括第k行)，col列中元素绝对值最大的列与第k行交换。如果col列中的元素全为0，那么则处理col + 1列，k不变。

3. 转换为行阶梯阵，判断解的情况。

① 无解

 当方程中出现(0, 0, …, 0, a)的形式，且a != 0时，说明是无解的。

② 唯一解

 条件是k = equ，即行阶梯阵形成了严格的上三角阵。利用回代逐一求出解集。

③ 无穷解。

 条件是k < equ，即不能形成严格的上三角形，自由变元的个数即为equ – k，即不确定的变元至少有var- k个。我们先把所有的变元视为不确定的。在每个方程中判断不确定变元的个数，如果大于1个，则该方程无法求解。如果只有1个变元，那么该变元即可求出，即为确定变元。

 

实例代码：

有多组测试数据。每组测试数据先输入一个整数n，表示方阵的阶。然后下面输入n阶方阵。输出其逆矩阵。若无逆矩阵，则输出No inverse matrix。

 

#include <iostream>

#include <cmath>

#include <algorithm>

using namespace std;

const double eps = 1e-6;

bool is_zero( const double num )//用于判断有无逆矩阵

{

    return fabs(num)< eps;

}

 

void create( double ** & matrix, const int n )

{

    matrix = newdouble* [n];

    for ( int i =0; i < n; ++i )

        matrix[i] =new double[n];

}

 

void input ( double ** matrix, const int n )

{

    for ( int i =0; i < n; ++i )

    {

        for ( intj  = 0; j < n; ++ j )

            cin>> matrix[i][j];

    }

}

 

bool inverse ( double ** matrix1, double ** matrix2,const int n )

{

    int i, j;

    for ( i = 0; i< n; ++ i )//初始化一个单位矩阵

    {

        for ( j =0; j < n; ++ j )

        {

            if ( i== j )

               matrix2[i][j] = 1;

            else

               matrix2[i][j] = 0;

        }

    }

    for ( i = 0; i< n; ++i )

    {

        introwmaxpos = i;

        for ( j = i+ 1; j < n; ++j )

        {

            if (matrix1[i][j] > matrix1[i][rowmaxpos] )

               rowmaxpos = j;

        }

        for ( j =i; j < n; ++ j )//按从大到小的顺序排列矩阵

        {

            swap(matrix1[j][rowmaxpos], matrix1[j][i]);

            swap(matrix2[j][rowmaxpos], matrix2[j][i]);

        }

        if (!is_zero(matrix1[i][i]) )

        {

            intdivisor = matrix1[i][i];

            for ( j= i; j < n; ++ j )//归一化矩阵

            {

                matrix1[i][j] /= divisor;

               matrix2[i][j] /= divisor;

            }

            for ( j= i + 1; j < n; ++ j )//高斯消元法处理行列式

            {

                intmultiple = matrix1[j][i];

                for( int k = i; k < n; ++ k )

                {

                   matrix1[i][j] -= matrix1[i][k] * multiple;

                   matrix2[i][j] -= matrix2[i][k] * multiple;

                }

            }

        }

        else

            returnfalse;

    }

    return true;

}

void output( double ** matrix, const int n )

{

    for ( int i =0; i < n; ++i )

    {

        for ( int j= 0; j < n; ++ j )

            cout<< matrix[i][j] << ' ';

       cout<<endl;

    }

}

void destroy( double ** matrix, const int n )

{

    for ( int i =0; i < n; ++ i )

        delete []matrix[i];

    delete []matrix;

}

int main()

{

    int n;

    double **matrix1;

    double **matrix2;

    while ( cin>> n )

    {

        create(matrix1, n );

        create(matrix2, n );

        input(matrix1, n);

        if (inverse(matrix1, matrix2, n) )

            output(matrix2, n );

        else

            cout<< "No inverse matrix" << endl;

        destroy(matrix1, n );

        destroy(matrix2, n );

    }

    return 0;

}

 

[POJ 1222 Extended Lights Out]是很多博客中提到的典型的求解线性方程组的题。

题意：有一个5*6的矩阵，每个元素代表1盏灯，0表示灯关，1表示灯开。若你改变了第a[i][j]盏灯的状态，相应地还会改变它四周a[i-1][j]、a[i][j-1]、a[i+1][j]、a[i][j+1]这些灯的状态。因为同一盏灯连续按两次又回到初始状态，所以每个灯只需按一次，若选择按第a[i][j]盏灯，则x[i][j] = 1，否则为0。现在给你这些灯的初始状态，让你找出一种方案x[][]使得所有灯全熄灭。

思路：对每盏灯，列一个异或方程，把可能影响它的灯和自己置1，其余灯置0。

比如说灯初始状态为a[i][j]，全灭的方案为x[i][j]。

则有——

a[i][j] ^ x[i-1][j]*1 ^ x[i][j-1]*1 ^ x[i+1][j]*1 ^x[i][j+1]*1 ^ x[i][j]*1 ^ (其余不相关的灯*0) =0。

两边同时异或a[i][j]得

x[i-1][j]*1 ^ x[i][j-1]*1 ^ x[i+1][j]*1 ^ x[i][j+1]*1 ^x[i][j]*1 ^ (其余不相关的灯*0) = a[i][j]。

30个方程组，30个变元，然后利用模板。

代码如下：（含模板）

 

#include<cstdio>  

#include<cstring>  

#include<algorithm>  

using namespace std;  

int a[31][31];//增广矩阵  

int equ, var;//equ个方程 var个变元  

int x[31];//解  

void init_A()  

{

 memset(a, 0, sizeof(a));  

 equ = 30, var = 30;//30个方程 30个变元  

 for(int i = 0; i < 5; i++)  

 {  

  for(int j = 0; j < 6; j++)  

  {  

   a[i*6+j][i*6+j] = 1;  

   if(i > 0) a[i*6+j][(i-1)*6+j] = 1;  

    if(j > 0) a[i*6+j][i*6+j-1] = 1;  

   if(i < 4) a[i*6+j][(i+1)*6+j] = 1;  

   if(j < 5) a[i*6+j][i*6+j+1] = 1;  

   }  

 }  

  for(int i = 0; i < 30; i++)  

  scanf("%d", &a[i][30]);  

}  

int Gauss()  

{  

  int k, max_r;  

  int col = 0;  

  for(k = 0; k < equ && col < var; k++, col++)  

  { //列主元法//找第col元素绝对值最大的行i(i > k) 与 当前行交换  

  max_r = k;  

  for(int i = k + 1; i < equ; i++)  

  if(a[i][col] > a[k][col]) max_r = i;  

  if(max_r != k)//找到 交换  

  {  

   for(int i = col; i < var+1; i++)  

   swap(a[k][i], a[max_r][i]);  

  }  

  if(a[k][col] == 0)//k行下面的元素全为0  

  {  

   k--;  

   continue;  

  }//加减消元  

  for(int i = k + 1; i < equ; i++)  

 {  

  if(a[i][col] != 0)  

  {  

  for(int j = col; j < var+1; j++)a[i][j] ^= a[k][j];  

   }  

  }  

 }//根据上三角矩阵 回代求解  

 for(int i = var-1; i >= 0; i--)  

 {  

   x[i] = a[i][var];  

  for(int j = i + 1; j < var; j++)  

  x[i] ^= (a[i][j] * x[j]);  

  }  

 return 0;  

}  

int main()  

{  

 int t, k = 1;  

 scanf("%d", &t);  

 while(t--)  

 {  

   init_A();  

   Gauss();  

    printf("PUZZLE #%d\n", k++);  

    for(int i = 0; i < 30; i++)  

    {  

    if((i+1) % 6 == 0)printf("%d\n", x[i]);  

    else printf("%d ", x[i]);  

    }  

   }  

 return 0;  

}  

还有几道复杂些的高斯消元题，太晚了明天整理吧。

以上~