解线性方程组迭代法之Guass-Seidel迭代法及其算法实现

上篇给大家讲解了迭代法中的Jacobi迭代法,这篇将给大家讲解一个收敛更快的方法:Guass-Seidel迭代法,在方法上,Guass-Seidel迭代法和Jacobi迭代法大同小异,它的迭代式与Jacobi迭代法的差不多,唯一的差别就是,在Jacobi迭代法中,我们是按照顺序依次求
X1^k+1,X2^k+1,…,Xn^k+1,即每个分量都由前一次迭代的Xi^k计算得到,而在这个过程中,我们求Xn^k+1时,前面的n-1项都已经求出来了,我们就可以将其替换成Xi^k+1(i<n)而后面i>n的部分,就按照前一次迭代的结果来计算,由于有了Jacobi迭代法的基础,这里就直接看代码实现
初始化

double** init_Matrix(int r, int c)
{
	double** p = new double* [r];
	int d = c + 1;
	for (int i = 0; i < r; i++)
	{
		p[i] = new double[d];
		memset(p[i], 0, sizeof(double) * d);
	}
	cout << "请输入线性方程组对应的增广矩阵:" << endl;
	for (int i = 0; i < r; i++)
	{
		for (int j = 0; j < d; j++)
		{
			cin >> p[i][j];
		}
	}
	return p;
}

判断精度

bool isRight(double** p, int r, double* x)
{
	double sum1 = 0, flag = 0, sum2 = 0;
	for (int i = 0; i < r; i++)
	{
		sum1 = 0, flag = 0;
		for (int j = 0; j < r; j++)
		{
			sum1 += x[j] * p[i][j];
		}
		flag = fabs(p[i][r] - sum1);
		if (flag > one_Precision)//解代入单个方程式的误差过大
		{
			return false;
		}
		else
		{
			sum2 += flag;
		}
	}
	if (sum2 > total_Precision)//整体误差过大
	{
		return false;
	}
	return true;
}

作用与Jacobi迭代法相同
迭代计算

void Iteration(double** p, int r, double* x)
{
	int k = 0;//最大迭代次数
	double sum = 0;
	while (true)
	{
		for (int i = 0; i < r; i++)
		{
			sum = 0;
			for (int j = 0; j < r; j++)
			{
				if (j == i)
				{
					continue;
				}
				else
				{
					sum -= p[i][j] * x[j];
				}
			}
			x[i] = (p[i][r] + sum) / p[i][i];
		}
		printf("第%d次迭代结果为:", ++k);
		for (int i = 0; i < r; i++)
		{
			printf("%f\t", x[i]);
		}
		cout << endl;
		if (k >= MAX_time)
		{
			cout << "超出迭代次数上限!停止迭代" << endl;
			return;
		}
		if (isRight(p, r, x))//精度符合要求
		{
			cout << "精度符合要求,停止迭代,共迭代:" << k << "次" << endl;
			return;
		}
	}
}

这里与Jacobi迭代法不同的地方就是,不需要存储上一次的运算结果,直接用每个解分量之前的解分量去求,未知的部分就用上一次的迭代结果
综合运用

void Gauss_seidel_main()
{
	int i = 0, j = 0;
	cout << "请输入线性方程组对应系数矩阵的行和列:" << endl;
	cin >> i >> j;
	double** p = init_Matrix(i, j);
	double* X = new double[i];//第n+1次跌代
	memset(X, 0, sizeof(double) * i);
	Iteration(p, i, X);
	for (int i = 0; i < j; i++)
	{
		delete[]p[i];
	}
	delete[]p;
	delete[]X;
}

同样别忘记了,释放内存,防止内存泄漏
完整代码及测试数据

#include<iostream>
#include<Windows.h>
#include<string>
#define one_Precision   1e-5//单个方程误差
#define total_Precision 3e-5//整体方程误差
#define MAX_time  1000//最大迭代次数
using namespace std;
/*
测试数据
3 3
10 -1  0 9
-1 10 -2 7
0  -2 10 6


3 3
20 -1 2 74
2   8 1 -4
1  -2 4 56

3 3
8 -3  2  20
4 11 -1  33
2  1  4  12

5 5
28 -3   0 0  0  10
-3 38 -10 0 -5  0
0 -10 25 -15 0  0
0  0 -15  45 0  0
0 -5   0  0  30 0
*/
double** init_Matrix(int r, int c)
{
	double** p = new double* [r];
	int d = c + 1;
	for (int i = 0; i < r; i++)
	{
		p[i] = new double[d];
		memset(p[i], 0, sizeof(double) * d);
	}
	cout << "请输入线性方程组对应的增广矩阵:" << endl;
	for (int i = 0; i < r; i++)
	{
		for (int j = 0; j < d; j++)
		{
			cin >> p[i][j];
		}
	}
	return p;
}
bool isRight(double** p, int r, double* x)
{
	double sum1 = 0, flag = 0, sum2 = 0;
	for (int i = 0; i < r; i++)
	{
		sum1 = 0, flag = 0;
		for (int j = 0; j < r; j++)
		{
			sum1 += x[j] * p[i][j];
		}
		flag = fabs(p[i][r] - sum1);
		if (flag > one_Precision)//解代入单个方程式的误差过大
		{
			return false;
		}
		else
		{
			sum2 += flag;
		}
	}
	if (sum2 >total_Precision)//整体误差过大
	{
		return false;
	}
	return true;
}
//用一维数组x存储每次迭代的解向量
void Iteration(double** p, int r, double* x)
{
	int k = 0;//迭代次数
	double sum = 0;
	while (true)
	{
		for (int i = 0; i < r; i++)
		{
			sum = 0;
			for (int j = 0; j < r; j++)
			{
				if (j == i)
				{
					continue;
				}
				else
				{
					sum -= p[i][j] * x[j];
				}
			}
			x[i] = (p[i][r] + sum) / p[i][i];
		}
		printf("第%d次迭代结果为:", ++k);
		for (int i = 0; i < r; i++)
		{
			printf("%f\t", x[i]);
		}
		cout << endl;
		if (k >= MAX_time)
		{
			cout << "超出迭代次数上限!停止迭代" << endl;
			return;
		}
		if (isRight(p, r, x))//精度符合要求
		{
			cout << "精度符合要求,停止迭代,共迭代:" << k << "次" << endl;
			return;
		}
	}
}
void Gauss_seidel_main()
{
	int i = 0, j = 0;
	cout << "请输入线性方程组对应系数矩阵的行和列:" << endl;
	cin >> i >> j;
	double** p = init_Matrix(i, j);
	double* X = new double[i];//第n+1次跌代
	memset(X, 0, sizeof(double) * i);
	Iteration(p, i, X);
	for (int i = 0; i < j; i++)
	{
		delete[]p[i];
	}
	delete[]p;
	delete[]X;
}
int main(void)
{
	Gauss_seidel_main();
	system("pause");
	return 0;
}