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