内容发布更新时间 : 2024/5/22 5:52:01星期一 下面是文章的全部内容请认真阅读。
解线性代数方程组的直接法
—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
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(\ } }