解线性代数方程组的直接法

—Gauss主元消去法

最近学习数值分析，为了能够边学习边训练编程，笔者决定对于所学算法一一分析记录，也便于自己以后复习时候的查阅，在这里顺便和大家共享一下我的学习经验，共勉

今天主要写的是求解线性方程组的直接法之一Gauss主元消去法 的基本算法步骤和用C++编程实现的具体过程，内容比较简单，代码已经经过测试，可以直接使用。

首先就是Gauss顺序消去法的基本原理：

主元消去法只是在顺序消元法进行前采取选主元措施。

对于Ax=b这样一非齐次线性方程组，要求解x列向量，利用线性代数方面的知识，我们可以对(A|b)这样一个增广矩阵做初等行变换，化为矩阵A下三角为零的形式。这个时候再进行回代操作就可以计算出来x列向量啦。

例如：

对矩阵Ax=b;

?a11　 a12　 a13　 ... a1n??b1????b?a　 a　 a ... a22232n??21?2?? ,b??...? A = ?...????a　a　a　... a?n?1,1n?1,2n?1,3n?1,n??bn?1????b?a　 a 　a 　... a?n?n,2n,3n,n?n,1?则对应A的增广矩阵如下：

?a11??a21(A|b)????an?1,1?a?n,1a12a22an?1,2an,2a13a23...a1n...a2n...an?1,3...an?1,nan,3...an,nb1??b2?...?

?bn?1?bn??第一步：以(A|b)的第一行为基，对的第2行到第n行进行初等行变换； 假设a11?0，将第1行的?ai1倍加到第i行(i=2,3,…,n)，通式为 a11(2)(1)aij?aij?ai1a1j　(j?1,2,...,n,n?1) 注意这里的j是从1到n+1取值的； a11aij的上标（1）,（2）只是用来表示数据aij前一时刻和当前时刻的值

经过消元，增广矩阵(A|b)变成下面的形式：

?a(1)?11?0(A|b)????0?0??(1)a12(2)a22(1)a13(2)a23...a1(1)n(2)...a2n...(2)an?1,2(2)an,2(2)(2)an...an?1,3?1,n(2)an,3(2)...an,nb1(1)??(2)b2?...?

?(2)bn?1??(2)bn??到这里，第一步就进行完啦，下来的第2步一直到第n-1步都和上面道理一

样，那么进行完n-1步以后，我们得到这样一个增广矩阵：

?a(1)?11?0(A|b)????0?0??(1)a12(2)a22(1)...a1,n?1(2)...a2,na1(1)n(2)a2n...00(n?1)...an?1,n?1(n?1)an?1,n(n?1)an,n...0b1(1)??(2)b2?...?

?(n?1)bn?1??(n?1)bn??嘿嘿，到这里基本就可以大功搞成了。

只要再进行回代操作很显然就可以求出来x列向量鸟。对于回代的公式是这样滴：

(n)?bn?xn?(n)an,n? ?n(k)?x?1[a(k)?a(k?n?1,n?2,...,1)?k,n?1kjxj]　　(k)?kakkj?k?1?有了上面的基础，那么下面咱们写C++实现就可谓是平步青云啦，磨刀不误砍柴工，我强烈建议爱学习的你能够在看明白上面的理论之后再进行C++编程过程。当然如果你已经对Guass顺序消元法烂熟于心，那就可以之间看看我的程序设计是否恰当了，如果我有的地方可能写的不够明白，如果你有什么不明白的地方可以给我留言，大家可以交流学习相互促进。

好！下面计入C++设计层次，我首先做了一个对于5元线性非齐次方程的求解的C++程序;

(1)顺序消元法的算法：

我采用了base, i，j 三个变量进行循环来完成顺序消元。

base代表当前的基行，i代表需要消元的行，j代表对应与i行增广矩阵的列，注意哦，这里是增广矩阵的列(A|b)，是包含方程有段向量的。 建立这样一个嵌套循环

for(base=0; base<=n-2; base++) {

double CurRatio = a[i][base]/a[base][base]; for(i=base+1; i<=n-1; i++) for(j=base; j<=n; j++) {

a[i][j] -= CurRatio * a[base][base];

}

}

注意base是从0-n-2 不是从1-n-1; a[i][j] -= CurRatio * a[base][base]; 这里为什么是写成这样而非

a[i][j] -= a[i][base]/a[base][base] * a[base][base]; 直截了当，嘿嘿，你想想。

在这之后，就可以进行回代操作了； 嗯，大致对于采用Gauss主元消去法直接计算5元线性方程组，大致要说的就是这些了

下面是具体做法：

头文件 Guass.h： #ifndef _GAUSS_H #define _GAUSS_H

//---------------------------------- //-- 高斯主元消去法

void Gauss_Run(double a[][5], double b[5], double x[5] ); #endif

Cpp文件Gauss.cpp #include \

//-------------------------------------- // 寻找base列主元

// 返回值为该列主元最大值所在行

//-------------------------------------- int findmaxbase(int base, double a[][5]) {

int maxmum; //代表最大行 int i;

maxmum = base;

for(i=base; i<5; i++)

if( a[i][base]>a[maxmum][base] ) maxmum = i;

return maxmum; }

//--------------------------------------- // 将数组第base行和第row行换行 // 结果存储在array数组中

//---------------------------------------

void ChangeArrayRow(int row, int base, double array[][6])

{

int j;

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

array[row][j] += array[base][j];

array[base][j] = array[row][j] - array[base][j]; array[row][j] = array[row][j] - array[base][j]; } }

//--------------------------------------- // 求解5元非齐次线性方程组

// a[5][5]是系数矩阵，b[5]是Ax=b等式右端向量; // 将计算结果存入x[5]数组

//---------------------------------------

void Gauss_Run(double a[][5], double b[5], double x[5] ) {

int base, i, j; //i代表行，j代表列，base代表基行 int row, k; //row代表行，k用于求解 double array[5][6]; double ProcData;

for(i=0; i<5; i++) for(j=0; j<6; j++) {

if(j==5) array[i][j] = b[i]; else array[i][j] = a[i][j]; }

for(base=0; base<=5-2; base++) {

// 寻找主元

row = findmaxbase(base, a); if(row!=base) {

ChangeArrayRow(row, base, array); }

//顺序消元

for(i=base+1; i<=5-1; i++) {

double CurRatio = array[i][base]/array[base][base]; for(j=base; j<=5; j++) {

array[i][j] -= CurRatio * array[base][j]; } } }

//回代求解

for(k=4; k>=0; k--) {

if(k==4) x[k] = array[k][5]/ array[k][k]; else {

ProcData = 0;

for(j = k+1; j<=4; j++) {

ProcData += array[k][j]*x[j]; }

x[k] = (array[k][5] - ProcData)/array[k][k]; } } }

main.cpp文件

#include #include \

void main() {

int i, j;

double A[5][5] = { 1, 2, 3, 4, 5, 2, 1, 3, 5, 4, 3, 2, 1, 4, 5, 5, 4, 2, 3, 1, 4, 5, 1, 2, 3 };

double b[5] = {6, 7, 8, 9, 10}; double x[5] = {0}; Gauss_Run(A, b, x); for(i=0; i<4; i++) {

printf(\ } }