矩阵变换---高斯全主元消去法---解线性非齐次方程组

写在前面，本文暂时是针对于有唯一解的非齐次线性方程组。

代码比较复杂，不是难，在文末，我把 “ 矩阵变换成三角矩阵 ” 的功能封装成了一个函数，不想看过程的可以直接使用。是非面向对象的。求逆矩阵也是类似方法。

---

1. 简介

下图是初始时的增广矩阵，解方程组的关键就是将矩阵变换成三角矩阵，于是此方程组的解为 [ 1, 2, 3, 4 ] ，具体的变换方法就是下面要介绍的高斯全主元消去法。

               

---

2. 算法思路：

对每一行依次处理

主元：第 i 行的主元是a [ i ] [ i ]，主元不能为零，如果为零就要找到合适的行与本行交换，产生新的主元。

• 判断该行的主元是否为0，如果是0，需要向下找主元所在列的第一个非零元素，主元所在行与该元素所在行交换，该元素代替成为新的主元。
• 把主元化为1（对该行的所有元素分别除以主元素）
• 用主元素将主元素所在列的其他元素消成0（假设要用a[0][0]消掉a[1][0],就让第0行元素乘以负的a[1][0]再加到第1行）
• 再对下一行用同样的方法处理，处理完最后一行之后，系数矩阵变成E，此时的b就是解。

---

3. 两个实例的运行过程：

画圈的是处理每一行时的主元

以下是运行过程：（第一个是增广矩阵）

0  0  2  1  10
4  1  1  1  13
1  1  1  2  14
1  0  1  0  4
--------------------------
4  1  1  1  13
0  0  2  1  10
1  1  1  2  14
1  0  1  0  4
--------------------------
1  0.25  0.25  0.25  3.25
0  0  2  1  10
1  1  1  2  14
1  0  1  0  4
--------------------------
1  0.25  0.25  0.25  3.25
0  0  2  1  10
0  0.75  0.75  1.75  10.75
0  -0.25  0.75  -0.25  0.75
--------------------------
1  0.25  0.25  0.25  3.25
0  0.75  0.75  1.75  10.75
0  0  2  1  10
0  -0.25  0.75  -0.25  0.75
--------------------------
1  0.25  0.25  0.25  3.25
0  1  1  2.33333  14.3333
0  0  2  1  10
0  -0.25  0.75  -0.25  0.75
--------------------------
1  0  0  -0.333333  -0.333333
0  1  1  2.33333  14.3333
0  0  2  1  10
0  0  1  0.333333  4.33333
--------------------------
1  0  0  -0.333333  -0.333333
0  1  1  2.33333  14.3333
0  0  2  1  10
0  0  1  0.333333  4.33333
--------------------------
1  0  0  -0.333333  -0.333333
0  1  1  2.33333  14.3333
0  0  1  0.5  5
0  0  1  0.333333  4.33333
--------------------------
1  0  0  -0.333333  -0.333333
0  1  0  1.83333  9.33333
0  0  1  0.5  5
0  0  0  -0.166667  -0.666667
--------------------------
1  0  0  -0.333333  -0.333333
0  1  0  1.83333  9.33333
0  0  1  0.5  5
0  0  0  -0.166667  -0.666667
--------------------------
1  0  0  -0.333333  -0.333333
0  1  0  1.83333  9.33333
0  0  1  0.5  5
0  0  0  1  4
--------------------------
1  0  0  0  1
0  1  0  0  2
0  0  1  0  3
0  0  0  1  4
--------------------------
1  0  0  0  1
0  1  0  0  2
0  0  1  0  3
0  0  0  1  4
x0=1
x1=2
x2=3
x3=4

---

4. 代码的具体实现：

#include<iostream>
using namespace std;
//本程序中所有用到的数组都是动态创建的
//二维矩阵用一维数组存放，所以必须用指针访问*(p+i)



class Matrix{ //矩阵类
	public:
		Matrix(int dims=2);
		void setMatrix(double * rmatr);//设置矩阵的值,必须要接收参数
		void printM();
		~Matrix();

	protected://用protected不用private保证子类依旧可以访问，为了使接下来的子类继续可以访问，用公有继承
		double * MatrixA; //存放一维数组首地址,为什么用double
		int index;//矩阵的维数
};

void Matrix::printM(){
	for(int i=0;i<index;i++){//i代表二维数组的行数
		for(int j=0;j<index;j++){
			cout<<*(MatrixA+i*index+j)<<" ";
		}
		cout<<endl;
	}

}
Matrix::Matrix(int dims){//初始化成员数据,声明时已经给了参数默认值，声明时就不能再写了
	MatrixA=new double[index*index];//分配成一维的形式
	index=dims;
}

Matrix::~Matrix(){
	delete[] MatrixA;
}

void Matrix::setMatrix(double * rmatr){
	//MatrixA = rmatr;
	for(int i=0;i<index*index;i++){
		MatrixA[i]=rmatr[i];
	}
}


class Linequ:public Matrix{ //线性方程组类
	private:
		double * solu;//存放解向量x的数组的首地址
		double * sums;//存放b的数组的首地址
	public:
		void setLinequ(double *a,double *b);//设置方程组,参数是干嘛的，A和b吗？？？
		void printL();//显示
		void Solve();//求解方程组，为什么返回值是int
		void showX();//输出方程的解x
		~Linequ();
		Linequ(int dims=2);//参数为阶数
};

Linequ::~Linequ(){ //清除动态分配的数组
	delete[] sums;
	delete[] solu;
}

//调用子类构造函数时，需要自动调用父类的构造函数
//构造函数的作用就是初始化成员数据
Linequ::Linequ(int dims):Matrix(dims){
	index=dims;
	solu = new double[index];
	sums = new double[index];
}


//线性方程组的系数矩阵是Matrix类所需要的
void Linequ::setLinequ(double *a,double *b){
	setMatrix(a);//初始化A矩阵，因为已经继承过来了，不用写Matrix::
	//sums=b;//只需要指针赋值就可以了
	//不知道为什么得循环赋值
	for(int i=0;i<index;i++){
		sums[i]=b[i];
	}
	
}
void Linequ::printL(){//输出增广矩阵
	 
	for(int i=0;i<index;i++){
		for(int j=0;j<index;j++){ 
			cout<<*(MatrixA+i*index+j)<<"  ";
			
		}
		cout<<sums[i]<<endl;
	
	}

}
void Linequ::Solve(){//MatrixA是系数矩阵，[MatrixA|sums]是增广矩阵
	double zhuyuan,a,c;int times;
	double *b = new double[index];

	for(int i=0;i<index;i++){//对每一行处理
		
		if(*(MatrixA+i*index+i)==0){//判断主元素MatrixA[i][i]是否为零
			//等于零要与下面的非零行交换
			times=i+1;//times记录下来下一行行号
			//判断具体要与哪一行交换
			while(i<index-1    &&   *(MatrixA+(i+1)*index+i)==0        ){//i不是最后一行
				times++;
			}
			for(int j=0;j<index;j++){//交换i行和times行，b为辅助数组
				b[j]=*(MatrixA+i*index+j);
				*(MatrixA+i*index+j)=*(MatrixA+times*index+j);
				*(MatrixA+times*index+j)=b[j];
			}
			c=sums[i];
			sums[i]=sums[times];sums[times]=c;
		}printL();
	cout<<"--------------------------"<<endl;
		//else{
			//1.  将主元化为1（对i行所有元素都除以主元的值）
			zhuyuan = *(MatrixA+i*index+i);//先把本轮主元的初始值保存下来，因为后面主元的值会变为1，该行其他元素的操作就变成除以1了
			for(int j=0;j<index;j++){
				*(MatrixA+i*index+j)=(*(MatrixA+i*index+j))/zhuyuan;
			}
			sums[i]=sums[i] /zhuyuan; 
			printL();
			cout<<"--------------------------"<<endl;
			//2.  将主元所在列的其他元素都用主元消成0，主元当前是i行i列，也就是接下来被处理的元素列为i,行不能为i
			for(int k=0;k<index;k++){//k代表行号,消去系数矩阵中位于[k][i]位置的元素
				if(k!=i){
					a = (*(MatrixA+k*index+i));//用变量a记录一下[k][i]位置的值
					//将i行元素乘以负的[k][i]再加到k行
					for(j=0;j<index;j++){
						(*(MatrixA+k*index+j))=(*(MatrixA+k*index+j)) + ( (*(MatrixA+i*index+j)) * (-1) * a );
					}
					sums[k]=sums[k] + ( sums[i] * (-1) * a );
				}
			}
	//	}
	printL();
	cout<<"--------------------------"<<endl;
	}
	for(int m=0;m<index;m++){
		solu[m]=sums[m];
	} 
}
void Linequ::showX(){
	for(int i=0;i<index;i++){
		cout<<"x"<<i<<"="<<solu[i]<<endl;
	}

}

void main(){//针对性处理四阶矩阵

	//系数矩阵
/*	double a[16]={0.2368,0.2471,0.2568,1.2671,
		0.1968,0.2071,1.2168,0.2271,
		0.1581,1.1675,0.1768,0.1871,
		1.1161,0.1254,0.1397,0.1490};
		double a[16]={2,1,1,1,
		1,2,1,1,
		1,1,2,1,
		1,1,1,2};*/
	double a[16]={0,0,2,1,
		4,1,1,1,
		1,1,1,2,
		1,0,1,0};

	//b向量
	//double b[4]={1.8471,1.7471,1.6471,1.5471};
		//double b[4]={11,12,13,14};
		double b[4]={10,13,14,4};
	Linequ equl(4);//新建一个线性方程组对象，4阶,在构造函数中分配了空间
	
	equl.setLinequ(a,b);//设置方程组，即将线性方程租对象具体赋值（用上面的数组）
	
	equl.printL();//将刚刚初始化的这个增广矩阵输出
	cout<<"--------------------------"<<endl;
	equl.Solve();
	equl.printL();
	equl.showX();
}

---

5. 封装的函数

#include<iostream>
using namespace std;

 void printL(int index,double *A,double*b){
	for(int i=0;i<index;i++){
		for(int j=0;j<index;j++){
			cout<<*(A+i*index+j)<<"  ";
		}
		cout<<b[i]<<endl;
	}
}



void Solve(int index,double*A,double*b){//index是维数
	double zhuyuan,a,c;int times;
	double *t = new double[index];

	for(int i=0;i<index;i++){//对每一行处理
		
		if(*(A+i*index+i)==0){//判断主元素MatrixA[i][i]是否为零
			//等于零要与下面的非零行交换
			times=i+1;//times记录下来下一行行号
			//判断具体要与哪一行交换
			while(i<index-1    &&   *(A+(i+1)*index+i)==0        ){//i不是最后一行
				times++;
			}
			for(int j=0;j<index;j++){//交换i行和times行，b为辅助数组
				t[j]=*(A+i*index+j);
				*(A+i*index+j)=*(A+times*index+j);
				*(A+times*index+j)=t[j];
			}
			c=b[i];
			b[i]=b[times];b[times]=c;
		}
		printL(index,A,b);
		cout<<"--------------------------"<<endl;
		//1.  将主元化为1（对i行所有元素都除以主元的值）
		zhuyuan = *(A+i*index+i);//先把本轮主元的初始值保存下来，因为后面主元的值会变为1，该行其他元素的操作就变成除以1了
		for(int j=0;j<index;j++){
			*(A+i*index+j)=(*(A+i*index+j))/zhuyuan;
		}
		b[i]=b[i] /zhuyuan; 
		printL(index,A,b);
		cout<<"--------------------------"<<endl;
		//2.  将主元所在列的其他元素都用主元消成0，主元当前是i行i列，也就是接下来被处理的元素列为i,行不能为i
		for(int k=0;k<index;k++){//k代表行号,消去系数矩阵中位于[k][i]位置的元素
			if(k!=i){
				a = (*(A+k*index+i));//用变量a记录一下[k][i]位置的值
				//将i行元素乘以负的[k][i]再加到k行
				for(j=0;j<index;j++){
					(*(A+k*index+j))=(*(A+k*index+j)) + ( (*(A+i*index+j)) * (-1) * a );
				}
				b[k]=b[k] + ( b[i] * (-1) * a );
			}
		}

	printL(index,A,b);
	cout<<"--------------------------"<<endl;
	}
	 
}


void main(){//实验数据：A[index*index]={0,0,2,1,4,1,1,1,1,1,1,2,1,0,1,0};b ={10,13,14,4};结果为{1,2,3,4}

	//动态定义数组
	int index = 4;
	double *A=new double[index*index];
	double *b=new double[index];

	//初始化数组
	cout<<"请输入系数矩阵A："<<endl;
	for(int i=0;i<index*index;i++){
		cin>>A[i];
	}
	cout<<"请输入b："<<endl;
	for(i=0;i<index;i++){
		cin>>b[i];
	}
	
	printL(index,A,b);//将刚刚初始化的这个增广矩阵输出
	cout<<"--------------------------"<<endl;
	
	//调用矩阵变换函数
	Solve(index,A,b); //最后b就是结果

	//输出结果
	for(i=0;i<index;i++){
		cout<<"x"<<i<<" = "<<b[i]<<endl;
	}
}