解线性代数方程组的直接法之Gauss主元消去法及其C++编程代码

解线性代数方程组的直接法—Gauss主元消去法最近学习数值分析，为了能够边学习边训练编程，笔者决定对于所学算法一一分析记录，也便于自己以后复习时候的查阅，在这里顺便和大家共享一下我的学习经验，共勉今天主要写的是求解线性方程组的直接法之一Gauss主元消去法的基本算法步骤和用C++编程实现的具体过程，内容比较简单，代码已经经过测试，可以直接使用。首先就是Gauss顺序消去法的基本原理：主元消去法只是在顺序消元法进行前采取选主元措施。对于Ax=b这样一非齐次线性方程组，要求解x列向量，利用线性代数方面的知识，我们可以对(A|b)这样一个增广矩阵做初等行变换，化为矩阵A下三角为零的形式。这个时候再进行回代操作就可以计算出来x列向量啦。例如：对矩阵Ax=b;111213121222321,11,21,31,,1,2,3,Annnnnnnnnnnnaaaaaaaaaaaaaaaa　　　...　　...=...　　　...　　　...,121...nnbbbbb则对应A的增广矩阵如下：11112131212223221,11,21,31,1,1,2,3,.........(A|b).........nnnnnnnnnnnnnnbaaaaaaaabaaaabaaaab第一步：以(A|b)的第一行为基，对的第2行到第n行进行初等行变换；假设110a，将第1行的111iaa倍加到第i行(i=2,3,…,n)，通式为(2)(1)1111(j1,2,...,n,n1)iijijjaaaaa　注意这里的j是从1到n+1取值的；aij的上标（1）,（2）只是用来表示数据aij前一时刻和当前时刻的值经过消元，增广矩阵(A|b)变成下面的形式：(1)(1)(1)(1)(1)11112131(2)(2)(2)(2)222232(2)(2)(2)(2)1,21,31,1(2)(2)(2)(2),2,3,...0...(A|b)......0...0...nnnnnnnnnnnnbaaaabaaaaaabaaab到这里，第一步就进行完啦，下来的第2步一直到第n-1步都和上面道理一样，那么进行完n-1步以后，我们得到这样一个增广矩阵：(1)(1)(1)(1)(1)111121,11(2)(2)(2)(2)2222,2(n1)(n1)(n1)1,11,1(n1)(n1),...0...(A|b)......00...00...0nnnnnnnnnnnnbaaaabaaaaabab嘿嘿，到这里基本就可以大功搞成了。只要再进行回代操作很显然就可以求出来x列向量鸟。对于回代的公式是这样滴：(n)(n),(k)(k),1(k)11[a](kn1,n2,...,1)nnnnnkknkjjjkkkbxaxaxa　　有了上面的基础，那么下面咱们写C++实现就可谓是平步青云啦，磨刀不误砍柴工，我强烈建议爱学习的你能够在看明白上面的理论之后再进行C++编程过程。当然如果你已经对Guass顺序消元法烂熟于心，那就可以之间看看我的程序设计是否恰当了，如果我有的地方可能写的不够明白，如果你有什么不明白的地方可以给我留言，大家可以交流学习相互促进。好！下面计入C++设计层次，我首先做了一个对于5元线性非齐次方程的求解的C++程序;(1)顺序消元法的算法：我采用了base,i，j三个变量进行循环来完成顺序消元。base代表当前的基行，i代表需要消元的行，j代表对应与i行增广矩阵的列，注意哦，这里是增广矩阵的列(A|b)，是包含方程有段向量的。建立这样一个嵌套循环for(base=0;base=n-2;base++){doubleCurRatio=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//----------------------------------//--高斯主元消去法voidGauss_Run(doublea[][5],doubleb[5],doublex[5]);#endifCpp文件Gauss.cpp#includeGauss.h//--------------------------------------//寻找base列主元//返回值为该列主元最大值所在行//--------------------------------------intfindmaxbase(intbase,doublea[][5]){intmaxmum;//代表最大行inti;maxmum=base;for(i=base;i5;i++)if(a[i][base]a[maxmum][base])maxmum=i;returnmaxmum;}//---------------------------------------//将数组第base行和第row行换行//结果存储在array数组中//---------------------------------------voidChangeArrayRow(introw,intbase,doublearray[][6]){intj;for(j=0;j6;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]数组//---------------------------------------voidGauss_Run(doublea[][5],doubleb[5],doublex[5]){intbase,i,j;//i代表行，j代表列，base代表基行introw,k;//row代表行，k用于求解doublearray[5][6];doubleProcData;for(i=0;i5;i++)for(j=0;j6;j++){if(j==5)array[i][j]=b[i];elsearray[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++){doubleCurRatio=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文件#includestdio.h#includeGauss.hvoidmain(){inti,j;doubleA[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};doubleb[5]={6,7,8,9,10};doublex[5]={0};Gauss_Run(A,b,x);for(i=0;i4;i++){printf(%.3f\n,x[i]);}}main.cpp文件主要是用来验证采用Gauss主元消去法直接计算5元线性方程组的解。当然，你看到上面的代码会觉得很不满意，对啊，只能够求解5元，线性方程组，那怎么能够满足我们的需要，我们所需求的是求解n元线性方程组的才可以啊。不然好好的算法不就白白浪费了么？……但是要解决n元线性方程组的求解，问题就来啦：首先，你得给我提供一个可以解决n元线性方程组的函数借口吧，上面的函数接口是voidGauss_Run(doublea[][5],doubleb[5],doublex[5]);那么n元的呢？再者，n是任意数啊，那岂不是要求我的A方阵的大小，b向量的大小是任意数啊，这个岂不是要求数组A[n][n]和b[n]是任意大小的?任意分配数组大小，当时学C语言谭浩强可没告诉我怎么做啊！嘿嘿，还好，办法总比困难多，函数接口的问题比较好办，我们只要再写一个接口就行了呗，我这里采用的是对Gauss_Run()重载，voidGauss_Run(double*a,double*b,double*x,intn);上面这么一个重载函数中使用了a指针，b指针，x指针，和一个n?这是用来做什么？嘿嘿，代码里面有解释，我也要上课去了，这里就不详细说了。对于数组任意分配大小的问题我们可以采用C++动态分配数组的方法。在上代码之前我再说一个问题，这也是当时最恶心我的一个问题，就是动态分配数组。我采用的方法是double**a=newdouble*[n];for(inti=0;in;i++){a[i]=new[n+1]();}注意使用new创建的动态数组我们是需要释放滴,所以会有这样一段代码for(i=0;in;i++)delete[]a[i];delete[]a;问题就出在了动态数组上，我采用的这种方法产生的动态数组内存是不连续的，那么你可想而知。如果我现在动态a[n][n+1]数组，在我的另外一个函数中，例如voidfunc1(double*arr,intn)想使用func1(a[0],n)进行对数组的操作，如果你对指针很熟悉的话，你肯定知道我如果想改变数组a[n][n+1]中a[i][j]的值我就可以采用*(arr+i*(n+1)+j)=2；假设是想让a[i][j]=2;那么咋一看这是没问题的，但是！！！兄弟或者姐们儿，现在我们开辟的是动态数组，这货把地址分配的时候是多么的万恶，它是不连续的！！所以上面这种美好的作法在这样万恶的情况下是不能实现的！！对于具体问题的阐述在代码中我有注释。嗯，说了这么多我觉得也就这些问题吧如果你有你的编程思想想也可以写出来我们大家共享哦，那么，上代码！头文件Guass.h：#ifndef_GAUSS_H#define_GAUSS_H//----------------------------------//--高斯主元消去法voidGauss_Run(double*a,double*b,double*x,intn);#endifCpp文件Gauss.cpp#includeGauss.h//---------------------------------------//重载函数//a--A[n][n