高斯列主元消元法求解线性方程组

一、高斯消去法的基本思想
    例1. 解方程组：
                

 

解 方程组矩阵形式为:AX=b,其中: 
            

 

第一步，消元过程：对增广矩阵进行消元

   

 

即得方程组

第二步, 回代过程:
             


此方法就是高斯消去法。

二、改进版算法

由高斯消去法知道，在消元过程中可能出现

=0的情况,
这时消去法将无法进行；即使主元素

≠0，但很小时，用其作除数,会导致其他元素数量级的严重增长和舍入误差的扩散,
最后也使得计算解不可靠。 
    引例 求解方程组
             

 

用4位浮点数进行计算。
   解： 方法1 用高斯消去法求解。

             


                                  
其中
         


计算解为:

         


显然计算解 

是一个很坏的结果，不能作为方程的近似解。
    方法2 交换行，避免绝对值小的主元做除数。
               

 

得计算解为：
         x=(-0.4900,-0.05113,0.3678)T ≈x*.
这个例子告诉我们，在采用高斯消去法解方程组时,小主元可能产生麻烦,故应避免采用绝对值小的主元素a 

。对一般矩阵来说，最好每一步选取系数矩阵(或消元后的低价矩阵)中绝对值最大的元素作为主元素，以使高斯消去法具有较好的数值稳定性,
这就是全主元素消去法, 在选主元时要花费较多机器时间，目前主要使用的是列主元消去法。 

本节主要介绍列主元消去法,并假定(2.1)的A∈Rn×n为非奇异的。
   1. 列主元素消去法
设方程组（2.1）的增广矩阵为：

                  


 
首先在A的第一列中选取绝对值最大的元素作为主元素，例如: 
         |ai1,1|= max |ai1|≠0， 1≤i≤n 
然后交换B的第一行与第i1 行，经第一次消元计算得
         (A|b)→(A(2)|b(2) ) 
 重复上述过程,设已完成第k-1步的选主元素,交换两行及消元计算，(A|b)约化为:

         

    (2.2)

其中A(k) 的元素仍记为aij ，b(k) 的元素仍记为bi 。

第k步选主元素(在A(k) 右下角方阵的第一列内选)，即确定ik ，使
              

 
交换(A(k) |b(k) ) 第k行与ik 列的元素，再进行消元计算,最后将原方
程组化为(k=1，2…，n-1):

                  

 
回代求解

         


    2. 高斯-若当消去法
高斯消去法始终是消去对角线下方的元素,现考察高斯消去法的一种修正，即消去对角线下方和上方的元素，这种方法称为高斯-若当(Gauss—Jordan)消去法。通过选主元，消元等过程最终化为： 

               

 
说明：用高斯-若当方法将A约化为单位矩阵,计算解就在常数位置得到,因此用不着回代求解,用高斯-若当方法解方程组其计算量要比高斯消去法大,但用高斯—若当方法求一个矩阵的逆矩阵还是比较合适的。
    定理4（高斯-若当法求逆矩阵）设A为非奇异矩阵，C=(A|In ),
如果对C应用高斯—若当方法化为（In|T），
则A-1=T 。

    例4 用高斯-若当方法求

 的逆矩阵以及

的解。

解：
         
         


c++代码实现如下：

/*
* test.cpp
*
*  Created on: 2014-1-5
*      Author: zhijian
*/

#include <stdio.h>
#include <math.h>

#define MAXN 100
#define ZERO 0.000000001	//定义浮点0
//#define DEBUG

//交换两个变量的值
void swap(double *a,double *b){
double temp = *a;
*a = *b;
*b = temp;
}

//B行的每一个元素加上A行的每一个元素的值*k
void lineOper(double A[],double B[],int n,double k){
for(int i = 0;i<n;i++)
B[i] += A[i] * k;
}
/*
* 解线性方程组Ax = b,
* 其中A为n*n的方阵，x，b为n维向量
* 返回是否有解
*/
bool solveEquations(double A[MAXN][MAXN],int n,double x[MAXN],double b[MAXN]){
double temp[MAXN][MAXN + 1];	//增广矩阵
//初始化增广矩阵
for(int i = 0;i<n;i++){
for(int j = 0;j<n;j++)
temp[i][j] = A[i][j];
temp[i]
= b[i];
}
//化为上三角矩阵
for(int j = 0;j<n;j++){
//寻找该列绝对值最大的行
int mmaxi = j;
for(int i = j+1;i<n;i++)
if(fabs(temp[i][j])>fabs(temp[mmaxi][j]))mmaxi = i;
if(fabs(temp[mmaxi][j])<ZERO)
return false;//无解
if(mmaxi != j){//交换两行
for(int i = 0;i<n+1;i++)
swap(&temp[mmaxi][i],&temp[j][i]);
}
//消元
for(int i = j + 1;i<n;i++){
double k = temp[i][j] / temp[j][j];
lineOper(temp[j],temp[i],n+1,-k);
}
}
//迭代求值
for(int j = n - 1;j>=0;j--){
double sum = 0;
for(int i = j+1;i<n;i++)
sum += x[i] * temp[j][i];
x[j] = (temp[j]
- sum) / temp[j][j];
}
return true;
}
void test(){
double A[MAXN][MAXN];
double x[MAXN];
double b[MAXN];
int n;
printf("请输入矩阵A的规模:\n");
scanf("%d",&n);
printf("请按行输入矩阵A:\n");
for(int i = 0;i<n;i++)
for(int j = 0;j<n;j++)
scanf("%lf",&A[i][j]);
printf("请输入向量b:\n");
for(int i = 0;i<n;i++)
scanf("%lf",&b[i]);
if(!solveEquations(A,n,x,b))
printf("无解\n");
else{
for(int i = 0;i<n;i++)
printf("%lf\n",x[i]);
}
}
int main(){
test();
return 0;
}

求逆的实现：

/*
* test.cpp
*
*  Created on: 2014-1-5
*      Author: zhijian
*/

#include <stdio.h>
#include <math.h>

#define MAXN 100
#define ZERO 0.000000001	//定义浮点0
#define DEBUG

//交换两个变量
void swap(double *a,double *b){
double temp = *a;
*a = *b;
*b = temp;
}

//交换两行变量
void swapLine(double A[],double B[],int n){
for(int i = 0;i<n;i++)
swap(&A[i],&B[i]);
}
//将一行（B）的每一个元素加上另外一行（A）的每一个元素*k
void lineOper(double A[],double B[],int n,double k){
for(int i = 0;i<n;i++)
B[i] += k * A[i];
}

/*
*
*
* 矩阵求逆
* 返回是否有逆
*/
bool inverse(double A[MAXN][MAXN],int n,double A_[MAXN][MAXN]){
double temp[MAXN][MAXN];	//辅助矩阵
//单位矩阵初始化
for(int i = 0;i<n;i++){
for(int j = 0;j<n;j++)
A_[i][j] = 0,temp[i][j] = A[i][j];
A_[i][i] = 1;
}
//化成上三角矩阵
for(int j = 0;j<n;j++){
int mmaxi = j;
for(int i = j + 1;i<n;i++)
if(fabs(temp[i][j])>fabs(temp[mmaxi][j]))mmaxi = i;
if(fabs(temp[mmaxi][j])<ZERO)return false;	//无逆
if(mmaxi != j){
swapLine(temp[mmaxi],temp[j],n);
swapLine(A_[mmaxi],A_[j],n);
}
for(int i = j+1;i<n;i++){
double k = temp[i][j] / temp[j][j];
lineOper(temp[j],temp[i],n,-k);
lineOper(A_[j],A_[i],n,-k);
}
}
//变成对角矩阵
for(int j = n-1;j>=0;j--){
for(int i = 0;i<j;i++){
double k = temp[i][j] / temp[j][j];
lineOper(temp[j],temp[i],n,-k);
lineOper(A_[j],A_[i],n,-k);
}
}
//单位化
for(int i = 0;i<n;i++){
double k = 1.0 / temp[i][i];
for(int j = 0;j<n;j++)
A_[i][j] *= k;
}
return true;
}

void test(){
double A[MAXN][MAXN];
double A_[MAXN][MAXN];
int n;
printf("输入方阵规模\n");
scanf("%d",&n);
for(int i = 0;i<n;i++)
for(int j = 0;j<n;j++)
scanf("%lf",&A[i][j]);
if(!inverse(A,n,A_)){
printf("无解\n");
}
else{
for(int i = 0;i<n;i++){
for(int j  = 0;j<n;j++)
printf("%lf ",A_[i][j]);
printf("\n");
}
printf("\n");
}
}
int main(){
test();
return 0;
}

参考链接：http://sxyd.sdut.edu.cn/zhanshi/shuzhifenxi/shuzhifenxi/2.1/szfx021.htm

http://sxyd.sdut.edu.cn/zhanshi/shuzhifenxi/shuzhifenxi/2.2/szfx022.htm