数值分析线性方程组迭代法之SOR迭代法详解及其C语言算法实现

SOR迭代法,又名逐次超松弛迭代法,与Jacobi迭代法和Guass-Seidel迭代法相比,收敛速度更快其原理如下(想详细了解,可以点击这里数值分析(东北大学)):
1.构造迭代式时,要加上一个大于0的松弛因子w,这样可以加快其收敛速度

2.根据上式进行分析:

3.得到迭代式:

得到迭代式以后,就可以选择合适的初始解进行计算了,由于迭代法的收敛性与初始向量无关,与系数矩阵的谱半径有关,所以在计算时的初始解向量不妨设为0向量即可

代码实现

1.初始化:

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];
		}
	}

	cout << "请输入SOR的松弛因子w(0<w):" << endl;
	while (!(cin >> w)||(cin.fail()))
	{
		if (w>0)
		{
			break;
		}
		else
		{
			cin.clear();
			cin.sync();
			cout << "输入不规范的w,请重新输入!" << endl;
		}
	}
	return p;
}

2.判断是否达到精度要求:

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;
}

3.迭代计算:

void Iteration(double** p, int r, double* x)
{
	int k = 0;//迭代次数
	double sum = 0,tmp=0;
	while (true)
	{
		for (int i = 0; i < r; i++)
		{
			sum = 0,tmp=0;
			for (int j = 0; j < r; j++)
			{
				if (j == i)
				{
					tmp -= x[j];
				}
				else
				{
					sum -= p[i][j] * x[j];
				}
			}
			x[i] = x[i]+w*((p[i][r] + sum) / p[i][i])+w*tmp;
		}
		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;
		}
	}
}

完整代码:

#include<iostream>
#include<Windows.h>
#include<math.h>
#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
3 2 -3 -2
1 1  1  6
1 2 -1  2

3 3
1  2 -3 1
2 -1  3 5
3 -2  2 1



4 4
4 -3  6  7 11
1  1  3  4 10
-2 9 -7  1 10
3  3 -4 20 25
*/
double w =0;//先赋初始值为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];
		}
	}

	cout << "请输入SOR的松弛因子w(0<w):" << endl;
	while (!(cin >> w)||(cin.fail()))
	{
		if (w>0)
		{
			break;
		}
		else
		{
			cin.clear();
			cin.sync();
			cout << "输入不规范的w,请重新输入!" << endl;
		}
	}
	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;
}
void Iteration(double** p, int r, double* x)
{
	int k = 0;//迭代次数
	double sum = 0,tmp=0;
	while (true)
	{
		for (int i = 0; i < r; i++)
		{
			sum = 0,tmp=0;
			for (int j = 0; j < r; j++)
			{
				if (j == i)
				{
					tmp -= x[j];
				}
				else
				{
					sum -= p[i][j] * x[j];
				}
			}
			x[i] = x[i]+w*((p[i][r] + sum) / p[i][i])+w*tmp;
		}
		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 SOR_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)
{
	SOR_main();
	system("pause");
	return 0;
}