高斯消元法一般用来计算线性方程的解，同时也可计算矩阵的秩，可逆矩阵的逆矩阵，以及行列式的值。
一般情况下，高斯消元是对一个有 n 个方程以及 n 个变量的方程组进行求解，形式如下：
A11* x1+ A 12* x2+ ... + A 1n* xn= B1
A21* x1+ A 22* x2+ ... + A 2n* xn= B2
.
.
A n1* x1+ A n2* x2+ ... + Ann* xn= Bn
其中 A ij 已知， Bi 已知，求 xi。
有关线性方程求解，先了解如下两个定理
定理一 线性方程组存在唯一解的充要条件是它的系数矩阵 A 的秩与其增广矩阵 A'= （A, B) 的秩相等。
定理二 设线性方程组 AX= B 的增广矩阵 A'= ( A, B ) 经初等行变换后所得到的矩阵为 C'= (C, D), 则矩阵 C' 所对应的线性方程组 CX= D 与原线性方程组 AX= B 同解，即它们有相同的解集合。
高斯消元法是基于定理二，经过一系列初等变换后，化为形如以下的方程
A11* X1+ A12* X2+ A13* X3+ ... + A1n* Xn= B1
A22* X2+ A23* X3+ ... + A2n* Xn= B2
A33* X3+ ... + A3n* Xn= B3
.
.
An* Xn= Bn
求解这样形式的方程组可以先求出 Xn，再求 Xn-1 ，依次即可求出所有 Xi 的值。
化为如上形式的方程组步骤如下：
1. 得到原线性方程组的增广矩阵。
2. 在第 c( c 从 1 开始 ) 列中找出一个行号大于 c 的绝对值最大的值 m，设该最大值在第 r 行。
3. 如果 r!= c，则将第 c 行与第 r 行的所有元素交换。
4. 对于所有第 k( k> c ) 行的值，我们对这行的第 j( j> c)个值执行以下变换，Akj= Akj- Akc/ m* Acj。
5. c 加上 1，转 2。
如下例，增广矩阵如下：
2 1 -2 10
3 2 2 1
5 4 3 4
先找出第 1 列行号大于1的绝对值最大的 5，与第 1 行交换，得到
5 4 3 4 L1
3 2 2 1 L2
2 1 -2 10 L3
将 L2 执行如下变换， L2= L2- 3/ 5* L1
将 L3 执行如下变换， L3= L3- 2/ 5* L1
得到
5 4 3 4 L1
0 -2/5 1/5 -7/5 L2
0 -3/5 -16/5 42/5 L3
在第 2 中找出行号大于 2 绝对值最大的 -3/5
L2 与 L3 交换，得到
5 4 3 4
0 -3/5 -16/5 42/5 L2
0 -2/5 1/5 -7/ 5 L3
执行如下变换 L3= L3- (-2/5)/ (-3/5)* L2
得到
5 4 3 4
0 -3/5 -16/5 42/5
0 0 7/ 3 -7
得到 X3= -7/ (7/ 3)= -3
X2= ( 42/5- (-16/5)* X3 )/ (-3/5)= 2
X1= ( 4- X3* 3- X2* 4 )/ 5= 1
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
int const N= 100;
double const eps= 1e-10;
double A[N][N], B[N];
int n;
void print(){
for( int i= 0; i< n; ++i ){
for( int j= 0; j< n+ 1; j++ )
printf("%.3lf ", A[i][j] );
puts("");
}
puts("");
}
int gause(){
for( int c= 0; c< n; ++c ){
double m= A[c][c];
int row= c;
for( int i= c+ 1; i< n; ++i ) // 找到绝对值最大的
if( fabs( A[i][c] )> fabs( m ) ){
m= A[i][c], row= i; }
if( fabs(m)< eps ) return 0;
if( c ){ // 交换
for( int i= 0; i< n+ 1; ++i ){
double t= A[c][i];
A[c][i]= A[row][i]; A[row][i]= t;
}
}
for( int i= c+ 1; i< n; ++i ){ // 变换
A[i][c]/= m;
for( int j= c+ 1; j< n+ 1; j++ )
A[i][j]-= A[i][c]* A[c][j];
A[i][c]= 0;
}
}
for( int i= n- 1; i>= 0; i-- ){
double t= 0;
for( int j= i+ 1; j< n; ++j )
t+= B[j]* A[i][j];
B[i]= ( A[i][n]- t )/ A[i][i];
}
return 1;
}
/*
3
2 1 -2 10
3 2 2 1
5 4 3 4
4
1 -1 5 -1 0
1 1 -2 3 0
3 -1 8 1 0
1 3 -9 7 0
*/
int main(){
while( scanf("%d",&n)!= EOF ){
for( int i= 0; i< n; ++i )
for( int j= 0; j< n+ 1; j++ )
scanf("%lf", &A[i][j] );
if( !gause() ){
puts("没有唯一解或者无解"); continue; }
print();
for( int i= 0; i< n; ++i ) printf("X%d: %.6lf\n", i, B[i] );
}
return 0;
}
|
阅读(1278) | 评论(0) | 转发(0) |
