constk

For_Ilyusha

Nov 18th, 2020
593
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. #include "SlaeSolver.h"
  2.  
  3. SlaeSolver::SlaeSolver(const ExtendedMatrix & other)
  4. {
  5.     matrix = new ExtendedMatrix(other);
  6.     roots = other.getColumn();
  7.     N = matrix->getRows();
  8. }
  9.  
  10. // -------------------- Метод Гаусса с выбором главного элемента --------------------
  11. vector<double> SlaeSolver::gaussMethod()
  12. {
  13.     ExtendedMatrix * matrixCopy = new ExtendedMatrix(*matrix);
  14.     double tmp = 0.0;
  15.  
  16.     // Прямой ход метода Гаусса
  17.     for (int n = 0; n != N; ++n)
  18.     {
  19.         tmp = matrix->at(n,n);
  20.         if (fabs(tmp) < 0.001)
  21.         {
  22.             double max = fabs(matrix->at(n + 1, n));
  23.             int iMax = n + 1;
  24.             for (int i = n + 2; i != N; ++i)
  25.             {
  26.                 if (fabs(matrix->at(i, n)) > max)
  27.                 {
  28.                     max = fabs(matrix->at(i, n));
  29.                     iMax = i;
  30.                 }
  31.             }
  32.             if (fabs(max) < 0.001)
  33.                 throw runtime_error("Error: Cant solve this system!");
  34.             else
  35.                 matrix->swap(n, iMax);
  36.         }
  37.         for (int i = N; i >= n; --i)
  38.             matrix->at(n,i) /= tmp;
  39.  
  40.         for (int i = n + 1; i != N; ++i)
  41.         {
  42.             tmp = matrix->at(i,n);
  43.             for (int j = N; j >= n; --j)
  44.                 matrix->at(i, j) -= tmp * matrix->at(n, j);
  45.         }
  46.     }
  47.  
  48.     // Обратный ход метода Гаусса
  49.     roots.at(N - 1) = matrix->at(N - 1, N);
  50.     for (int i = N - 2; i >= 0; --i)
  51.     {
  52.         roots.at(i) = matrix->at(i, N);
  53.         for (int j = i + 1; j != N; ++j)
  54.             roots.at(i) -= matrix->at(i, j) * roots.at(j);
  55.     }
  56.  
  57.     matrix = matrixCopy;
  58.     return roots;
  59. }
  60.  
  61. // -------------------- L-U разложение --------------------
  62. vector<double> SlaeSolver::LuDecompositionMethod()
  63. {
  64.     if (L == NULL && U == NULL) // Если L и U матрицы ещё не рассчитаны
  65.     {
  66.         ExtendedMatrix * matrixCopy = new ExtendedMatrix(*matrix);
  67.         double tmp = 0.0;
  68.  
  69.         L = new SquareMatrix(N);
  70.  
  71.         // Прямой ход метода Гаусса
  72.         for (int n = 0; n != N; ++n)
  73.         {
  74.             tmp = matrix->at(n, n);
  75.             if (fabs(tmp) < 0.001)
  76.             {
  77.                 double max = fabs(matrix->at(n + 1, n));
  78.                 int iMax = n + 1;
  79.                 for (int i = n + 2; i != N; ++i)
  80.                 {
  81.                     if (fabs(matrix->at(i, n)) > max)
  82.                     {
  83.                         max = fabs(matrix->at(i, n));
  84.                         iMax = i;
  85.                     }
  86.                 }
  87.                 if (fabs(max) < 0.001)
  88.                     throw runtime_error("Error: Cant solve this system!");
  89.                 else
  90.                     matrix->swap(n, iMax);
  91.             }
  92.  
  93.             L->at(n, n) = tmp; // Заполняем матрицу L (элементы главной диагонали)
  94.             for (int i = N; i >= n; --i)
  95.                 matrix->at(n, i) /= tmp;
  96.  
  97.             for (int i = n + 1; i != N; ++i)
  98.             {
  99.                 tmp = matrix->at(i, n);
  100.                 L->at(i, n) = tmp; // Заполняем матрицу L (элменты под главной диагональю)
  101.                 for (int j = N; j >= n; --j)
  102.                     matrix->at(i, j) -= tmp * matrix->at(n, j);
  103.             }
  104.         }
  105.  
  106.         // Берём матрицу U как матрицу с единицами на главной диагонали, получившуюся после
  107.         // выполнения прямого хода метода Гаусса
  108.         U = new SquareMatrix(matrix->getMatrix());
  109.         for (int i = 0; i != N; ++i)
  110.             U->at(i, i) = 1.0;
  111.  
  112.         matrix = matrixCopy;
  113.     }
  114.  
  115.     // Находем матрицу-столбец Y по формуле
  116.     vector<double> Y(N);
  117.     Y.at(0) = matrix->getColumn().at(0) / L->at(0, 0);
  118.     for (int i = 1; i != N; ++i)
  119.     {
  120.         Y.at(i) = matrix->getColumn().at(i);
  121.         for (int j = 0; j != i; ++j)
  122.             Y.at(i) -= L->at(i, j) * Y.at(j);
  123.         Y.at(i) /= L->at(i, i);
  124.     }
  125.    
  126.     // Находим корни по формуле
  127.     roots.at(N - 1) = Y.at(N - 1);
  128.     for (int i = N - 2; i != -1; i--)
  129.     {
  130.         roots.at(i) = Y.at(i);
  131.         for (int j = i + 1; j != N; ++j)
  132.             roots.at(i) -= U->at(i, j) * roots.at(j);
  133.     }
  134.  
  135.     return roots;
  136. }
  137.  
  138. // Задать столбец свободных членов
  139. void SlaeSolver::setColumn(const vector<double> column)
  140. {
  141.     this->matrix->setColumn(column);
  142.     roots = column;
  143. }
  144.  
  145. // Задать СЛАУ целиком
  146. void SlaeSolver::setSlae(const ExtendedMatrix & extendedMatrix)
  147. {
  148.     *(this->matrix) = extendedMatrix;
  149.     L->clear();
  150.     U->clear();
  151.     roots.clear();
  152.     roots = extendedMatrix.getColumn();
  153. }
  154.  
  155. // -------------------- Метод прогонки --------------------
  156. vector<double> SlaeSolver::sweepMethod()
  157. {
  158.     TridiagonalMatrix * m = new TridiagonalMatrix(this->matrix->getMatrix());
  159.  
  160.     bool notEnoughPresicionWarning = false;
  161.  
  162.     L->clear();
  163.     U->clear();
  164.  
  165.     int n = this->matrix->getRows();
  166.     vector<double> alphaArr;
  167.     vector<double> betaArr;
  168.     double alpha = 0.0, beta = 0.0, gamma = 0.0, a = 0.0, b = 0.0, c = 0.0, d = 0.0;
  169.  
  170.  
  171.     // Первая итерация прямого хода
  172.     b = m->getMainDiagonal().at(0);
  173.     c = m->getAboveDiagonal().at(0);
  174.     d = matrix->getColumn().at(0);
  175.  
  176.     if (!(fabs(b) >= fabs(a) + fabs(c) && fabs(b) > fabs(a)))
  177.         notEnoughPresicionWarning = true;
  178.  
  179.     gamma = b;
  180.     beta = d / gamma;
  181.     alpha = -c / gamma;
  182.     if (gamma == 0)
  183.         throw runtime_error("Can't solve this! Gamma is 0!");
  184.  
  185.     alphaArr.push_back(alpha);
  186.     betaArr.push_back(beta);
  187.  
  188.     // Прямой ход
  189.     for (int i = 1; i < n - 1; ++i)
  190.     {
  191.         a = m->getBelowDiagonal().at(i - 1);
  192.         b = m->getMainDiagonal().at(i);
  193.         c = m->getAboveDiagonal().at(i);
  194.         d = matrix->getColumn().at(i);
  195.  
  196.         if (!(fabs(b) >= fabs(a) + fabs(c) && fabs(b) > fabs(a)))
  197.             notEnoughPresicionWarning = true;
  198.  
  199.         gamma = b + a * alphaArr.at(i - 1);
  200.         beta = (d - a * betaArr.at(i - 1)) / gamma;
  201.         alpha = -c / gamma;
  202.         if (gamma == 0)
  203.             throw runtime_error("Can't solve this! Gamma is 0!");
  204.  
  205.         alphaArr.push_back(alpha);
  206.         betaArr.push_back(beta);
  207.     }
  208.  
  209.     // Последняя итерация прямого хода
  210.     a = m->getBelowDiagonal().at(n - 2);
  211.     b = m->getMainDiagonal().at(n - 1);
  212.     d = matrix->getColumn().at(n - 1);
  213.  
  214.     if (!(fabs(b) >= fabs(a) + fabs(c) && fabs(b) > fabs(a)))
  215.         notEnoughPresicionWarning = true;
  216.  
  217.     gamma = b + a * alphaArr.at(n - 2);
  218.     beta = (d - a * betaArr.at(n - 2)) / gamma;
  219.     if (gamma == 0)
  220.         throw runtime_error("Can't solve this! Gamma is 0!");
  221.  
  222.     roots.at(n - 1) = beta;
  223.  
  224.  
  225.     // Обратный ход
  226.     for (int i = n - 2; i > -1; --i)
  227.     {
  228.         alpha = alphaArr.back();
  229.         alphaArr.pop_back();
  230.  
  231.         beta = betaArr.back();
  232.         betaArr.pop_back();
  233.  
  234.         roots.at(i) = alpha * roots.at(i + 1) + beta;
  235.     }
  236.  
  237.     if (notEnoughPresicionWarning)
  238.         cout << "WARNING! Loss of accuracy is possible using sweep method!" << endl << endl;
  239.  
  240.     return roots;
  241. }
  242.  
  243. // -------------------- Тестирование --------------------
  244. void SlaeSolver::showRoots()
  245. {
  246.     cout << "Roots:" << endl;
  247.     for (int i = 0; i != roots.size(); ++i)
  248.         cout << "\t" << roots.at(i) << endl;
  249. }
  250.  
  251. void SlaeSolver::test()
  252. {
  253.     cout << "Slae matrix:" << endl;
  254.     matrix->show();
  255.     cout << endl;
  256.     showRoots();
  257.     cout << endl << "Checking roots:" << endl;
  258.     Tester * tester = new Tester(this->matrix, roots);
  259.     tester->test();
  260. }
  261.  
RAW Paste Data