mmn_test_func

#include <cmath>
#include <iomanip>
#include <iostream>
#include <vector>

double f(double x, double y) {
  double T = M_PI * x * y;
  double sin2T = std::pow(std::sin(T), 2);
  double cos2T = std::pow(std::cos(T), 2);
  double y2_x2 = y * y + x * x;
  return -2 * M_PI * std::exp(std::pow(std::sin(T), 2)) *
         (2 * M_PI * y2_x2 * cos2T * sin2T - M_PI * y2_x2 * sin2T +
          M_PI * y2_x2 * cos2T);
}

double u(double x, double y) {
  return std::exp(std::pow(std::sin(M_PI * x * y), 2));
}

std::vector<std::vector<double>> MinimalResiduals(const int nmax = 1000,
                                                  const double _eps = 0.0000001,
                                                  const int _n = 3,
                                                  const int _m = 3) {
  int Nmax = nmax; // максимальное число итераций (не менее 1)
  std::vector<std::vector<double>> v; // сеточная функция ѵ (x, y)
  v.resize(_n + 1);
  for (size_t i = 0; i < _n + 1; i++) {
    v[i].resize(_m + 1);
    for (size_t j = 0; j < _m + 1; j++) {
      v[i][j] = 0;
    }
  }
  std::vector<std::vector<double>> r_vec; // вектор невязки
  r_vec.resize(_n + 1);
  for (size_t i = 0; i < _n + 1; i++) {
    r_vec[i].resize(_m + 1);
    for (size_t j = 0; j < _m + 1; j++) {
      r_vec[i][j] = 0;
    }
  }
  std::vector<std::vector<double>> ar; // вектор невязки
  ar.resize(_n + 1);
  for (size_t i = 0; i < _n + 1; i++) {
    ar[i].resize(_m + 1);
    for (size_t j = 0; j < _m + 1; j++) {
      ar[i][j] = 0;
    }
  }
  double r1 = 0;
  int S = 0;         // счетчик итераций
  double eps = _eps; // заданная точность
  double eps_max = 0; // точность, достигнутая на текущей итерации
  double eps_cur = 0; // для подсчета точности на текущей итерации
  double a2, k2, h2; // ненулевые элементы матрицы (-А)
  const int n = _n, m = _m; // размерность сетки
  double r_max = 0;
  double a = 0;
  double b = 1;
  double c = 0;
  double d = 1; // границы области определения уравнения
  int i, j; // индексы
  double r = 0;
  double v_old; // старое значение преобразуемой компоненты вектора
  double v_new; // новое значение преобразуемой компоненты вектора и
  bool flag = false; // условие остановки
  double h = ((b - a) / n);
  double k = ((d - c) / m);
  double z = 0;
  double z_max = 0;
  double tau_s = 0;
  double num = 0; // tau=num/denominator
  double denominator = 1;
  h2 = -std::pow((n / (b - a)), 2);
  k2 = -std::pow((m / (d - c)), 2);
  a2 = -2 * (h2 + k2);
  for (int i = 0; i < n + 1; i++) {
    v[i][0] = u(a + i * h, a);
    v[i][m] = u(a + i * h, b);
  }
  for (int j = 0; j < m + 1; j++) {
    v[0][j] = u(c, c + j * k);
    v[n][j] = u(d, c + j * k);
  }
  do {
    z = 0;
    /* std::cout << "вектор невязки:" << std::endl; */
    /* for (int i = n - 1; i >= 0; i--) { */
    /*   for (int j = 0; j < m; j++) */
    /*     std::cout << std::left << std::setw(10) << r_vec[i][j]; */
    /*   std::cout << "\n"; */
    /* } */
    /* std::cout << "промежуточное решение:" << std::endl; */
    /* for (int i = n - 1; i >= 0; i--) { */
    /*   for (int j = 0; j < m; j++) */
    /*     std::cout << std::left << std::setw(10) << v[i][j]; */
    /*   std::cout << "\n"; */
    /* } */
    z_max = 0;
    eps_max = 0;
    r_max = -20;
    num = 0;
    denominator = 0;
    for (j = 1; j < m; j++) {
      for (i = n - 1; i > 0; i--) {
        r_vec[i][j] = f(a + i * h, c + j * k) -
                      (h2 * (v[i + 1][j] + v[i - 1][j]) +
                       k2 * (v[i][j + 1] + v[i][j - 1]) + v[i][j] * a2);
      }
    }
    for (j = 1; j < m; j++) {
      for (i = n - 1; i > 0; i--) {
        ar[i][j] =
            (h2 * (r_vec[i + 1][j] + r_vec[i - 1][j]) +
             k2 * (r_vec[i][j + 1] + r_vec[i][j - 1]) + r_vec[i][j] * a2);
        num += ar[i][j] * r_vec[i][j];
        denominator += ar[i][j] * ar[i][j];
      }
    }
    tau_s = num / denominator;
    for (j = 1; j < m; j++) {
      for (i = n - 1; i > 0; i--) {
        v_old = v[i][j];
        v_new = v_old + tau_s * r_vec[i][j];
        eps_cur = std::fabs(v_old - v_new);
        if (eps_cur > eps_max) {
          eps_max = eps_cur;
        }
        v[i][j] = v_new;
        r = fabs(f(a + i * h, c + j * k) -
                 (h2 * (v[i + 1][j] + v[i - 1][j]) +
                  k2 * (v[i][j + 1] + v[i][j - 1]) + v_old * a2));
        if (r > r_max)
          r_max = r;
        z = u(a + i * h, c + j * k) - v[i][j];
        if (z > z_max)
          z_max = z;
      }
    }

    S = S + 1;
    if ((eps_max <= eps) or (S >= Nmax)) {
      flag = true;
    }
  } while (!flag);

  std::cout << "\nСПРАВКА:" << std::endl;

  std::cout << "Количество выполненных итераций = " << S << std::endl;
  std::cout << "Точность на выходе = " << eps_max << std::endl;
  std::cout << "Невязка = " << r_max << std::endl;
  std::cout << "Норма общей погрешности = " << z_max << std::endl;
  return v;
}

int main() {
  int n = 100;
  int m = 100;
  std::vector<std::vector<double>> v(n, std::vector<double>(m, 0));
  v = MinimalResiduals(100000, 0.000001, n, m);
  n++;
  m++;
  /* std::cout << "Результат:" << std::endl; */
  /* for (int i = n - 1; i >= 0; i--) { */
  /*   for (int j = 0; j < m; j++) */
  /*     std::cout << std::left << std::setw(10) << v[i][j]; */
  /*   std::cout << "\n"; */
  /* } */
  return 0;
}