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

一、高斯消去法的基本思想
例1. 解方程组：
高斯列主元消元法求解线性方程组

解方程组矩阵形式为: AX=b ,其中:
高斯列主元消元法求解线性方程组

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

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

即得方程组

第二步, 回代过程:
高斯列主元消元法求解线性方程组

此方法就是高斯消去法。

二、改进版算法

由高斯消去法知道，在消元过程中可能出现高斯列主元消元法求解线性方程组 =0的情况, 这时消去法将无法进行；即使主元素≠0，但很小时，用其作除数,会导致其他元素数量级的严重增长和舍入误差的扩散, 最后也使得计算解不可靠。
    引例求解方程组
              高斯列主元消元法求解线性方程组

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



其中

计算解为:



显然计算解高斯列主元消元法求解线性方程组是一个很坏的结果，不能作为方程的近似解。
    方法2 交换行，避免绝对值小的主元做除数。

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

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

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

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

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

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

(2.2)

其中 A^(k) 的元素仍记为 a_ij ，b^(k) 的元素仍记为 b_i 。

第 k 步选主元素 ( 在 A^(k) 右下角方阵的第一列内选 ) ，即确定 i_k ，使
高斯列主元消元法求解线性方程组

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

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

回代求解

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

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

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

秒客网

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

相关文章