Нестационарный метод Зейделя
Пытаюсь реализовать т.н. нестационарный метод Зейделя, в котором используется перестановка уравнений в зависимости от вектора поправок. Без перестановки уравнений метод сходится за 8 итераций, с перестановкой - или растет число шагов или вовсе наблюдается расходимость. Подскажите, пожалуйста, что я делаю неверно. Проблемное место в коде закомментировал.
#include <iostream>
#include <algorithm>
#include <time.h>
#include <cmath>
#include <math.h>
#include <algorithm>
#include <iomanip>
using namespace std;
const double error = 0.00001; // число, близкое к нулю
void fillMatrixPredefined(double**& arr, int &size)
{
// Выделяем память и заполняем матрицу значениями
int N = 7,k=1;
size=5;
arr = new double*[size];
for(int i=0; i<size; i++)
{
arr[i] = new double[size+1];
for(int j=0; j<size; j++)
{
if(i<j)
arr[i][j] = (double)(i+1.0) + (double)(j+1.0) - (double)(N/3.0) - (double)k;
else if(i>j)
arr[i][j] = (double)(i+1.0) + (double)(j+1.0) - (double)(N/5.0) - (double)k;
else
arr[i][j] = (double)(i+1.0) + (double)(j+1.0) + (double)(N/4.0) + (double)k;
}
arr[i][size] = 3.0 * (double)(i+1.0) + (double)(N/2.0) + (double)k;
}
}
/*void fillMatrixPredefined(double**& arr, int &size)
{
// Выделяем память и заполняем матрицу значениями
size=3;
arr = new double*[size];
for(int i=0; i<size; i++)
{
arr[i] = new double[size+1];
}
arr[0][0]=12;
arr[0][1]=4;
arr[0][2]=7;
arr[0][3]=5;
arr[1][0]=5;
arr[1][1]=9;
arr[1][2]=3;
arr[1][3]=11;
arr[2][0]=3;
arr[2][1]=4;
arr[2][2]=11;
arr[2][3]=9;
}*/
void addRow(double**& arr, const int &size, const int &row1,
const int &row2, const double &coeff1, const double &coeff2)
{
for(int i=0; i<(size+1); i++)
{
arr[row1][i] = arr[row1][i] * coeff1 + arr[row2][i] * coeff2;
}
}
void swapRows(double*& row1, double*& row2)
{
// обмен двух строк матрицы
// путем обмена указателей на строки
double* temp = row1;
row1 = row2;
row2 = temp;
}
bool isDDM(double** arr, const int &size)
{
// проверим, является ли матрица
// диагонально доминирующей
double rowSum = 0.0;
for(int i=0; i<size; i++)
{
rowSum = -abs(arr[i][i]);
for(int j=0; j<size; j++)
{
rowSum += abs(arr[i][j]);
}
if(rowSum > abs(arr[i][i])) return false;
}
return true;
}
bool isBConvergent(double** arr, const int &size)
{
for(int i=0; i<size; i++)
{
double rowSum = 0.0;
for(int j=0; j<size; j++)
{
if(j!=i)
rowSum += abs(arr[i][j]);
}
if(rowSum >=1.0) return false;
}
return true;
}
void printMatrix(double** arr, const int& size)
{
// Выводим матрицу на экран
cout << "----------------------------------" << endl;
for(int i=0; i<size; i++)
{
for(int j=0; j<size+1;j++)
{
cout << setw(8) << setprecision(3) << arr[i][j] << " ";
}
cout << endl;
}
cout << "----------------------------------" << endl;
}
bool reachedAccuracy(const double* prev, const double* cur, const int& size)
{
for(int i=0; i<size; i++)
{
if(abs(prev[i]-cur[i]) > error)
return false;
}
return true;
}
int main()
{
//cout << fixed << showpoint;
int size = 0;
double** arr = nullptr;
fillMatrixPredefined(arr, size);
double* x_K_plus_1 = new double[size];
double* x_K = new double[size];
double* correction = new double[size];
std::fill(x_K_plus_1, x_K_plus_1+size,0);
std::fill(x_K, x_K+size, 0);
std::fill(correction, correction+size,0);
// Приводим исходную матрицу к дигональному преобладанию
printMatrix(arr, size);
cout << "Is diagonally dominating? " << isDDM(arr,size) << endl;
// A[5] = A[5] - 5*A[1]
addRow(arr, size, 4, 0,1.0, -5.0);
// A[4] = 2*A[4] - 10*A[2]
addRow(arr, size, 3, 1,2.0,-10.0);
// A[3] = 12*A[3] - 6*A[2]
addRow(arr, size, 2, 1,12.0,-6.0);
// A[2] = 9*A[2] + A[4] - A[1]
addRow(arr, size, 1, 3,9.0,1.0);
addRow(arr, size, 1, 0,1.0,-1.0);
// A[1] = 100*A[1] + 22*A[5] - A[3] + A[4]
addRow(arr, size, 0, 4,100.0,22.0);
addRow(arr, size, 0, 2,1.0,-1.0);
addRow(arr, size, 0, 3,1.0,1.0);
swapRows(arr[0],arr[4]); // Обмен строк 1 и 5
swapRows(arr[1],arr[3]); // Обмен строк 2 и 4
printMatrix(arr, size);
cout << "Is diagonally dominating? " << isDDM(arr,size) << endl;
// Преобразуем систему к виду x = bx + B
// Нормализуем матрицу (чтобы суммы по строкам не превышали единицы)
// Строить новую матрицу необходимости нет
// (просто не будем в дальнейшем использовать главную диагональ)
// Вектор свободных членов B является последним столбцом имеющейся расширенной матрицы
for(int i=0; i<size; i++)
{
double divisor = arr[i][i];
for(int j=0; j<(size+1); j++)
{
arr[i][j] /= divisor;
if(j!=size) arr[i][j] *= -1;
} arr[i][i] = 0.0;
}
printMatrix(arr, size);
if(isBConvergent(arr, size))
cout << "Method converges, we can proceed" << endl;
else {
cout << "Method does NOT converge, please check matrices" << endl;
return -1;
}
// Инициализируем вектор (x)^(0) начального приближения
for(int i=0; i<size; i++)
x_K[i] = arr[i][size];
// Простые итерации
int counter = 0;
int test = 0;
int order[] = { 0,1,2,3,4};
while(!reachedAccuracy(x_K, x_K_plus_1, size))
{
std::copy(x_K_plus_1, x_K_plus_1+size, x_K);
std::fill(x_K_plus_1, x_K_plus_1+size,0);
/*for(int i=0; i<size; i++)
{
if ((*std::max_element(correction, correction+size))<error) {
break;
}
int k = std::distance(correction,std::max_element(correction, correction+size));
std::swap(order[k], order[i]);
std::swap(correction[k], correction[i]);
correction[k]=0.0f;
}*/
for(int i=0; i<size; i++)
{
for(int j=0; j<(size+1); j++)
{
if(order[i]>j)
x_K_plus_1[order[i]] += arr[i][j] * (j==size ? 1 : x_K_plus_1[j]);
else
x_K_plus_1[order[i]] += arr[i][j] * (j==size ? 1 : x_K[j]);
}
}
for(int i=0; i<size; i++)
{
correction[i] = abs(x_K_plus_1[i] - x_K[i]);
}
counter++;
}
std::cout << "--------------------------------------------" << endl;
for(int i=0; i<size; i++)
{
cout << setprecision(8) << "X_" << i << "=" << x_K_plus_1[i] << endl;
}
cout << "--------------------------------------------" << endl;
cout << "Performed " << counter << " simple iterations" << endl;
// освобождаем память, выделенную ранее для матрицы
for(int i=0; i<size; i++)
{
delete[] arr[i];
}
delete[] arr;
delete[] x_K_plus_1;
delete[] x_K;
delete[] correction;
system("pause");
}