Joyce' Blog
0%

完全主元高斯消去法解线性方程组

发表于 更新于

计算方法，老师很好，思路清晰，娓娓道来，然而，忘得比学的快，实在无奈…嗯还是完成今天课上说的用完全主元高斯消去法解线性方程组吧。先记一下自己踩过的几个坑：

  • 记录完全主元时，不小心定义成了int类型，导致进行取舍后系数为0，返回无解的提示。应定义为double；
  • 上课没理解老师说的“换回去”是什么意思，结果跑出来才大彻大悟——调整完全主元的位置时，行交换不会影响解的顺序，但是第i.j列进行交换时，xi,xj的值也会被交换。最后在修改代码时增加了数组order，用于记录xi的位置变化；
  • 最麻烦的可能还是三层循环里的i.j.k的初始条件吧，调整了好几次……
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
/*
-----完全主元高斯消去法解线性方程组-----
输入：未知数的个数N，N*N的系数矩阵A，N*1的向量Y
输出：N*1的向量X

数据结构：矩阵、线性方程组
*/
#include"LinearEqu.h"
int main()
{
	double a[] = {//方程系数矩阵
		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 b[] = { 1.8471, 1.7471, 1.6471, 1.5471 };//方程右端项
	LinearEqu equ(4);//定义一个四元方程组对象
	equ.setLinearEqu(a, b);//设置方程组
	if (!equ.solve())
		cout << "Fail" << endl;
	system("pause");
	return 0;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
#pragma once
#include<iostream>
#include<iomanip>
using namespace std;
class Matrix
{
public:
	double* elements;
	Matrix(int size = 2);
	~Matrix();
	void setMatrix(double* values);
	double element(int i, int j) {return elements[i*size + j];}
	void printMatrix(double* p,int s);
	int size;
};
1
2
3
4
5
6
7
8
9
10
11
12
13
14
#pragma once
#include"Matrix.h"
class LinearEqu : public Matrix
{
public:
	LinearEqu(int size = 2);
	~LinearEqu();
	void setLinearEqu(double* a, double* b);
	bool solve();
	void printLinearEqu();
	void printSolution();
	double* sums;
	double* solution;
};
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
#include"Matrix.h"
Matrix::Matrix(int size):size(size) {
	elements = new double[size*size];
}
Matrix::~Matrix() { delete[] elements; }
void Matrix::setMatrix(double* values)
{
	for (int i = 0; i < size*size; i++)
		elements[i] = values[i];
}
void Matrix::printMatrix(double* p, int s)
{
	cout << "The Matrix is:" << endl;
	for (int i = 0; i <size; i++)
	{
		for (int j = 0; j < s; j++)
			cout <<setw(12)<<p[i*s + j] ;
		cout << endl;
	}
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
#include"LinearEqu.h"
void swap(double &a, double &b)
{
	double tmp = a;
	a = b;
	b = tmp;
}
LinearEqu::LinearEqu(int size):Matrix(size) {
	sums = new double[size];
	solution = new double[size];
}
LinearEqu::~LinearEqu() {
	delete[]sums;
	delete[] solution;
}
void LinearEqu::setLinearEqu(double* a, double* b) {
	setMatrix(a);
	for (int i = 0; i < size; i++)
		sums[i] = b[i];
}
bool LinearEqu::solve() {
	printMatrix(elements, size);
	int* order = new int[size];
	for (int i = 0; i < size; i++)
		order[i] = i + 1;
	for (int i = 0; i < size; i++)
	{
		int im =i,jm = i;//最大元所在的行号列号
		double max = elements[i*size + i];
		for (int j = i; j < size; j++)
		{			
			for (int k = j; k < size; k++)
			{
				if (max < elements[j*size + k])
				{
					max = elements[j*size + k];
					jm = k;
					
					im = j;
				}
			}
			
		}
		if (max == 0)
			return false;
		else
		{

			if (im != i)//行交换
			{
				for (int k = i; k < size; k++)
				{
					swap(elements[im*size + k], elements[i*size + k]);
					swap(sums[im], sums[i]);
				}
			}
			if (jm != i)//列交换
			{
				swap(order[jm], order[i]);
				for (int k = 0; k < size; k++)
				{
					swap(elements[k*size + i], elements[k*size + jm]);
				}
			}
		}
		//消去
		double m = elements[i*size + i];
		for (int j = i + 1; j < size; j++)
		{
			double n = elements[j*size + i] / m;
			for (int k = i ; k < size; k++)
			{
				elements[j*size + k] = elements[j*size + k]- elements[i*size + k] * n;
			}
			sums[j] = sums[j] - sums[i] * n;
		}
		//判断剩下的一个元素是否等于0（等于0则矩阵的秩不等于矩阵的行数，则无解）
		//这里没有想好，待补充

		printMatrix(elements, size);
	}
	

	//回代
	for (int i = size - 1; i >= 0; i--)
	{
		double m = elements[i*size + i];
		for (int j = i - 1; j >= 0; j--)
		{
			double n = elements[j*size + i] / m;
			elements[j*size + i] = 0;
			sums[j] = sums[j] - sums[i] * n;
		}
	}
	
	for (int i = 0; i < size; i++)
		cout << "x[" << order[i] << "] = " << (solution[i] = sums[i] / elements[i*size + i]) << endl;
	return true;
}
void LinearEqu::printLinearEqu() {
	printMatrix(elements,size);
	printMatrix(sums, 1);
}
void LinearEqu::printSolution() {
	printMatrix(solution, 1);
}