fork download
copy 
  1. #include <iostream> 
  2. #include <locale>
  3. #include <cmath>
  4. #include <vector>
  5. #include <algorithm>
  6. #include <fstream>
  7. #include <iomanip>
  8. #include <cstdlib> // Для вызова system()
  9.  
  10. #define EPSILON 0.001
  11.  
  12. using namespace std;
  13.  
  14. void printVector(const std::vector<double>& vec) {
  15.     for (size_t i = 0; i < vec.size(); ++i) {
  16.         std::wcout << L"x" << (i + 1) << L" = " << std::fixed << std::setprecision(10) << vec[i];
  17.         if (i < vec.size() - 1) {
  18.             std::wcout << L", ";
  19.         }
  20.     }
  21.     std::wcout << std::endl;
  22. }
  23.  
  24. bool isConvergent(const vector<vector<double>>& matrix) {
  25.     // Проверяем, удовлетворяет ли матрица условию диагонального преобладания для сходимости.
  26.     int n = matrix.size();
  27.     for (int i = 0; i < n; i++) {
  28.         double sum = 0;
  29.         for (int j = 0; j < n; j++) {
  30.             if (i != j) {
  31.                 sum += abs(matrix[i][j]);
  32.             }
  33.         }
  34.         if (abs(matrix[i][i]) <= sum) {
  35.             return false;
  36.         }
  37.     }
  38.     return true;
  39. }
  40.  
  41. double calculateResidual(const vector<vector<double>>& A, const vector<double>& x, const vector<double>& b) {
  42.     double residual = 0.0;
  43.     int n = A.size();
  44.     for (int i = 0; i < n; ++i) {
  45.         double Ax_i = 0.0;
  46.         for (int j = 0; j < n; ++j) {
  47.             Ax_i += A[i][j] * x[j];
  48.         }
  49.         residual += pow(Ax_i - b[i], 2);
  50.     }
  51.     return sqrt(residual);
  52. }
  53.  
  54. vector<double> gaussSeidel(vector<vector<double>>& A, vector<double>& b, double relaxation, int& iterations) {
  55.     if (relaxation == 0.0) {
  56.         // При релаксации = 0 решаем стандартным методом Гаусса-Зейделя
  57.         relaxation = 0.01;
  58.     }
  59.     int n = A.size();
  60.     vector<double> x(n, 0.0); // Начальное приближение устанавливаем в 0.0
  61.     iterations = 0;
  62.  
  63.     double maxDiff;
  64.     do {
  65.         maxDiff = 0.0;
  66.         for (int i = 0; i < n; i++) {
  67.             double sigma = 0.0;
  68.             for (int j = 0; j < n; j++) {
  69.                 if (j != i) {
  70.                     sigma += A[i][j] * x[j];
  71.                 }
  72.             }
  73.             double newValue = (1 - relaxation) * x[i] + relaxation * (b[i] - sigma) / A[i][i];
  74.             maxDiff = max(maxDiff, fabs(newValue - x[i]));
  75.             x[i] = newValue;
  76.         }
  77.         iterations++;
  78.     } while (maxDiff > EPSILON); // Ограничение на количество итераций
  79.  
  80.     return x;
  81. }
  82.  
  83. void plotResults(const vector<double>& relaxation_values, const vector<int>& iteration_counts) {
  84.     ofstream dataFile("data.dat");
  85.     for (size_t i = 0; i < relaxation_values.size(); ++i) {
  86.         dataFile << relaxation_values[i] << " " << iteration_counts[i] << "\n";
  87.     }
  88.     dataFile.close();
  89.  
  90.  
  91. }
  92.  
  93. int main() {
  94.     setlocale(LC_ALL, "Russian");
  95.     vector<vector<double>> A = {
  96.         {6, 0, 3, -1},
  97.         {3, 13, -1, 2},
  98.         {5, 7, 15, 1},
  99.         {2, 7, -1, 11}
  100.     };
  101.     vector<double> b = { 16, 8, 85, 52 };
  102.  
  103.     if (!isConvergent(A)) {
  104.         wcout << L"Данная матрица не удовлетворяет условиям сходимости." << endl;
  105.         return -1;
  106.     }
  107.     int k = 0;
  108.     vector<double> relaxation_values;
  109.     vector<int> iteration_counts;
  110.     for (double relaxation = 0.1; relaxation <= 1.8; relaxation += 0.1) {
  111.         int iterations = 0;
  112.  
  113.         vector<double> solution = gaussSeidel(A, b, relaxation, iterations);
  114.         if (relaxation >= 0.9 && k == 0) {
  115.             wcout << fixed << L"Параметр релаксации: " << relaxation << L", Число итераций: " << iterations << endl;
  116.             wcout << L"Решение: ";
  117.             wcout.precision(10);
  118.             printVector(solution);
  119.             // Вычисляем и выводим конечную невязку
  120.             double residual = calculateResidual(A, solution, b);
  121.             wcout << L"Конечная невязка (||Ax - b||): " << residual << endl;
  122.             k = 1;
  123.  
  124.         }
  125.  
  126.  
  127.  
  128.         // Сохраняем результаты для построения графика
  129.         relaxation_values.push_back(relaxation);
  130.         iteration_counts.push_back(iterations);
  131.     }
  132.  
  133.  
  134.  
  135.     // Построение графика с помощью Gnuplot
  136.     plotResults(relaxation_values, iteration_counts);
  137.  
  138.     return 0;
  139. }
  140.  
Success #stdin #stdout 0.01s 5272KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
???????? ??????????: 1.000000, ????? ????????: 5
???????: x1 = 0.9999295622, x2 = 0.0000498321, x3 = 5.0000120494, x4 = 4.9999821910
???????? ??????? (||Ax - b||): 0.0005590819
https://ideone.com/bmy2NG
language:
C++ (gcc 8.3)
created:
10 days ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!