线性方程组的求解

ONE PIECE2

已于 2024-01-11 22:44:56 修改

阅读量6.6k

点赞数 11

分类专栏：数学文章标签：线性代数

于 2022-05-08 17:59:14 首次发布

本文链接：https://blog.csdn.net/qq_63788960/article/details/124643992

版权

数学专栏收录该内容

2 篇文章 1 订阅

订阅专栏

克莱姆法则

求解线性方程组有一种比较简单易行的方法就是用克莱姆法则通过行列式的计算以解出方程，下面给出行列式解方程的代码并分析优缺点；

对于一个n元一次方程组，如果可以将其化为n阶行列式就能使用克莱姆法则；例如：

$a_{11}x_{1}+a_{12}x_{2}+a_{13}x_3+...+a_{1n}x_{n} = b_1 \\ a_{21}x_{1}+a_{22}x_{2}+a_{23}x_3+...+a_{2n}x_{n} = b_2\\ ... ~~~~~~~~~~... ~~~~~~~~~... ~~~~~~~~~~~~~~... ~~~~~~...\\ ... ~~~~~~~~~~... ~~~~~~~~~... ~~~~~~~~~~~~~~... ~~~~~~...\\ a_{n1}x_{1}+a_{n2}x_{2}+a_{n3}x_3+...+a_{nn}x_{n} = b_n$

有 D=

用（b1,b2,...bn)T替换D的第一列、第二列、。。第n列得到D1D2。。Dn

行列式计算的代码：

double rec(int n, vector<double>D) {
	if (n == 1) {
		return D.back();
	}
	int e = 1;//设置符号
	double sum = 0;
	for (int k = 0; k < n; k++) {
		vector<double>d;
		for (int i = 1; i < n; i++) {//i==1表示从第一行开始，因为第0行是被展开的量；
			for (int j = 0; j < n; j++) {
				if (k == j)continue;     //划去同列元素；
				d.push_back(D[i*n + j]);
 
			}
		}
		sum += D[k] * e*rec(n - 1, d);
		e *= -1;
	}
	return sum;
}

分别计算出各个行列式的值后

X1=D1/D,X2=D2/D..Xn=Dn/D;非常的简单；

#include<iostream>
#include<vector>
using std::vector;
double rec(int n, vector<double>D) {
	if (n == 1) {
		return D.back();
	}
	int e = 1;//设置符号
	double sum = 0;
	for (int k = 0; k < n; k++) {
		vector<double>d;
		for (int i = 1; i < n; i++) {//i==1表示从第一行开始，因为第0行是被展开的量；
			for (int j = 0; j < n; j++) {
				if (k == j)continue;     //划去同列元素；
				d.push_back(D[i*n + j]);

			}
		}
		sum += D[k] * e*rec(n - 1, d);
		e *= -1;
	}
	return sum;
}
int main() {
	int n;
	std::cin >> n;
	vector<double>D[20], B;
	//得到D的行列式与B
	for (int i = 0; i < n; i++) {
		for (int j = 0; j < n; j++) {
			double d;
			std::cin >> d;
			D[0].push_back(d);
		}
		double d;
		std::cin >> d;
		B.push_back(d);
	}
	//利用B得到Di
	for (int k = 0; k < n; k++) {
		D[k + 1] = D[0];
		for(int i=0;i<n;i++)
		D[k + 1][i*n + k] = B[i];
	}
	//求解；
	double d = rec(n, D[0]);
	for (int i = 1; i <= n; i++) {
		double x = rec(n, D[i]) / d;
		std::cout << "x" << i << "=" << x << "\n";
	}
	return 0;
}

克莱姆法则计算的优点：

1：不会在中间计算中丢失精度，对于整数系数的方程可比较简单地通过gcd直接消除精度差

2：行列式计算的代码写好后，整个的代码就很好写了，没多少细节处理；

缺点：

时间复杂度为阶乘，无法求出10组以上的方程组

高斯消元

先看简单的情况，n个方程组：

对于上面的方程组可以写成一个增广矩阵；

$\begin{pmatrix} a11 , a12 ,......... ,b1 \\ a21, a22,.........,b2\\ a31,a32.......,....b3\\ ....\\ an1,.an2.........,bn \end{pmatrix}$

若将其化为如下的上三角矩阵，就能很容易的进行求解了：

$\begin{pmatrix} 1 , a12 ,... ,b1 \\ 0, 1,......,b2\\ 0,0,1,....b3\\ ....\\ 0,.0..,1,br\\ 0,0,0,..,0 \end{pmatrix}$

从第一行第一列开始，设当前行为r，列为c，

找行：找出r及r下面所有行中，第c列元素绝对值最大的一行，与当前行交换；

消元：对于该行的第一个非0元素，即第c列，将其化为1，利用该行将下面所有与这个1所在的列的元素消为0；（若该行第c列为0，也就是找出的绝对值最大值为0，说明有一组无效方程，则将c进行+1）；

开一个二维数组存储一个n行n+1列的矩阵：

算法流程：为方便计算将行列从0开始

for(r=0,c=0;c<n;c++)://r表示当前行，c表示当前列，r同时计算矩阵的有效方程数；

1：k=r

2：for(i=r+1;i<n;i++):if(|akc|<|aic|)k=i//得到第c列元素绝对值最大的行；

3：if(|akc|==0)continue;//第c列最大值为0，则该列全为0直接进行下一层循环；

4：swap(ar,ak)//交换这两行；

5：arj/=arc//使该行第一个非0元素，即第c列化为1

6：for(i=r+1,i<n;i++):aij-=aic*arj//利用第r行将下面的第c列消为0；

7：r++//有效方程数加一

矩阵变换代码：

for (r = 0, c = 0; c < n; c++) {
		int k = r;
		for (int i = r + 1; i < n; i++)
			if (fabs(a[i][c]) > fabs(a[k][c]))
				k = i;
		if (fabs(a[k][c]) < eps)continue;//如果c列最大的一行趋于0，之后的行同样趋于0，则该行视为无效方程组
		for (int i = c; i <= n; i++)std::swap(a[r][i], a[k][i]);//交换第r行与第k行
        for (int i = n; i >= c; i--)a[r][i] /= a[r][c];//c列元素归一：
		for (int i = r + 1; i < n; i++) 
			for (int j = n; j >= c; j--) 
				a[i][j] -= a[i][c] * a[r][j];//消元
			 
		
		r++;//秩+1
	}

$\begin{pmatrix} 1 , a12 ,... ,b1 \\ 0, 1,......,b2\\ 0,0,1,....b3\\ ....\\ 0,.0..,1,br\\ 0,0,0,..,0 \end{pmatrix}$

得到该矩阵后，

若r==n,显然xn=bn;解得xn后可以代到上一行解出xn-1。。。。最后解出x1

因此可以从下往上依次算出所有x

回代代码块：

for (int i = n - 1; i >= 0; i--) {
			for (int j = i + 1; j < n; j++) {
				a[i][n] -= a[i][j] * a[j][n];
			}
		}

解得xi=ain

若r<n,如果存在第r行后的b不为0，则左右两式矛盾（等式左边全部系数消为0,则等式左边为0）

否则有无穷多个解

完整代码:

#include<iostream>
#include<vector>
#include<cstring>
#include<iomanip>
const double eps = 1e-5;
using std::vector;
double a[107][107];
int n;
void read() {
	std::cin >> n;
	for (int i = 0; i < n; i++) {
		for (int j = 0; j <= n; j++) {
			std::cin >> a[i][j];
		}
	}
}
void print() {
	for (int i = 0; i < n; i++) {
		std::cout << "x" << i + 1<<" = " << a[i][n] << "\n";
	}
}

void solve() {
	read();//读入数据
	int r, c;
	for (r = 0, c = 0; c < n; c++) {
		int k = r;
		for (int i = r + 1; i < n; i++)
			if (fabs(a[i][c]) > fabs(a[k][c]))
				k = i;
		if (fabs(a[k][c]) < eps)continue;//如果c列最大的一行趋于0，之后的行同样趋于0，则该行视为无效方程组
		for (int i = c; i <= n; i++)std::swap(a[r][i], a[k][i]);//交换第r行与第k行
        for (int i = n; i >= c; i--)a[r][i] /= a[r][c];//c列元素归一：
		for (int i = r + 1; i < n; i++) 
			for (int j = n; j >= c; j--) 
				a[i][j] -= a[i][c] * a[r][j];//消元
			 
		
		r++;//秩+1
	}
    //得到一个上三角的矩阵：

	if (r < n) {//无解或无穷解
		for (int i = r ; i < n; i++) {
			if (fabs(a[i][n]) > eps) {
				std::cout << "方程无解\n";
				return;
			}
		}
		//不存在矛盾式子则有无穷解
		std::cout << "方程有无穷解\n";
	}
	else {
		//利用上三角的性质，解出最后一项从下到上往回代
		for (int i = n - 1; i >= 0; i--) {
			for (int j = i + 1; j < n; j++) {
				a[i][n] -= a[i][j] * a[j][n];
			}
		}
		print();//输出解
	}
}
int main() {
	solve();
	return 0;
}

现在看n元，m个方程组的情况；

若m<n;那么就将上面的代码行数改为m就行，最后就会得出有解或无解的答案；

若m>n; r<n依然得到正解，r==n，因为下面还会有方程组即：

$\begin{pmatrix} 1 , a12 ,... ,b1 \\ 0, 1,......,b2\\ 0,0,1,....b3\\ ....\\ 0,.0..,1,br\\ 0,0,0,..,br+1\\ 0,0,0,...,br+2\\ 0,0,0,...,bm \end{pmatrix}$

还需要判断是否矛盾；因此将判断无解的代码提出来；

for (int i = r; i < m; i++) {
		if (fabs(a[i][n]) > eps) {
			std::cout << "无解";
			return;
		}
	}
	if (r < n) {
		std::cout << "无穷解";
	}
	else {
		for (int i = n - 1; i >= 0; i--) {
			for (int j = i + 1; j < n; j++) {
				a[i][n] -= a[i][j] * a[j][n];
			}
		}
		print();
	}

实际上最开始n个方程组的代码也可这样写，这是因为将数组定义再全局中自动初始化为0了；

最终代码：

#include<iostream>
#include<vector>
#include<cstring>
#include<iomanip>
#include<cmath>
const double eps = 1e-9;
using std::vector;
double a[107][107];
int n, m;
void read() {
	std::cin >> m >> n;//m个方程组，n个元;
	for (int i = 0; i < m; i++) {
		for (int j = 0; j <= n; j++) {
			std::cin >> a[i][j];
		}
	}
}
void print() {
	for (int i = 0; i < n; i++) {
		std::cout << "x" << i + 1 << " = " << a[i][n] << "\n";
	}
}
void solve() {
	read();
	int r, c;
	for (r = 0, c = 0; c < n; c++) {
		int k = r;
		for (int i = r + 1; i < m; i++)
			if (fabs(a[i][c]) > fabs(a[k][c]))
				k = i;
		if (fabs(a[k][c]) < eps)continue;
		for (int i = c; i <= n; i++)std::swap(a[r][i], a[k][i]);
		for (int i = n; i >= c; i--)a[r][i] /= a[r][c];
		for (int i = r + 1; i < m; i++)
			for (int j = n; j >= c; j--)
				a[i][j] -= a[i][c] * a[r][j];
		r++;
	}
	for (int i = r; i < m; i++) {
		if (fabs(a[i][n]) > eps) {
			std::cout << "无解";
			return;
		}
	}
	if (r < n) {
		std::cout << "无穷解";
	}
	else {
		for (int i = n - 1; i >= 0; i--) {
			for (int j = i + 1; j < n; j++) {
				a[i][n] -= a[i][j] * a[j][n];
			}
		}
		print();
	}
}
int main() {
	solve();
	return 0;
}

高斯消元法的缺点：需要处理的细节较多，中间计算过程对精度的处理需要一定的技巧（例如不使用临时变量保存double型变量）

优点：1:时间复杂度为 $O(n^3)$

2:可以解n元m个方程组

ONE PIECE2

关注

11
点赞
踩
27

收藏

觉得还不错? 一键收藏
1
评论
线性方程组的求解

求解线性方程组有一种比较简单易行的方法就是用克莱姆法则通过行列式的计算以解出方程，下面给出行列式解方程的代码并分析优缺点；对于一个n元一次方程组，如果可以将其化为n阶行列式就能使用克莱姆法则；例如：a11x1+a12x2+a13x3+.....+a1nxn=b1a21x1+a22x2+a23x3+.....+a2nxn=b2. .ai1x1+ai2x2+ai3x3+......
复制链接

扫一扫