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