<8/7>集训日记
2017-08-07 22:40
232 查看
今天看了有关高斯消元的一些内容。
高斯消元法是线性代数中的一个算法,可用来求解线性方程组、求出矩阵的秩、求出可逆方阵的逆矩阵。
下面是一篇博客中讲的高斯消元法在编程中的应用。
首先,先介绍程序中高斯消元法的步骤:
(我们设方程组中方程的个数为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 !=
16ccb
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*
;
for ( int i =0; i < n; ++i )
matrix[i] =new double
;
}
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;
}
还有几道复杂些的高斯消元题,太晚了明天整理吧。
以上~
高斯消元法是线性代数中的一个算法,可用来求解线性方程组、求出矩阵的秩、求出可逆方阵的逆矩阵。
下面是一篇博客中讲的高斯消元法在编程中的应用。
首先,先介绍程序中高斯消元法的步骤:
(我们设方程组中方程的个数为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 !=
16ccb
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*
;
for ( int i =0; i < n; ++i )
matrix[i] =new double
;
}
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;
}
还有几道复杂些的高斯消元题,太晚了明天整理吧。
以上~