Guest User

Untitled

a guest
Feb 28th, 2012
3,268
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. // Метод Якоби (метод простой итерации) решения СЛАУ
  2. // СЛАУ задается в файле input.txt, ее решение выдается в output.txt
  3. // Допустим, надо решить такую систему:
  4. // 4*x1 -   x2 + 0*x3 = 2
  5. //  -x1 + 4*x2 -   x3 = 6
  6. // 0*x1 -   x2 + 4*x3 = 2
  7. // Тогда файл input.txt должен иметь следующий вид:
  8. //  3
  9. //  4 -1  0  2
  10. // -1  4 -1  6
  11. //  0 -1  4  2
  12.  
  13. #include <iostream>
  14. #include <fstream>
  15. #include <vector>
  16. using namespace std;
  17.  
  18. const double epsilon = 0.0001; // требуемая точность
  19.  
  20. void ReadLinearSystem(string file_name,
  21.                       vector< vector<double> >& linear_system)
  22. {
  23.   linear_system.resize(0);
  24.  
  25.     ifstream file_linear_system(file_name);
  26.    
  27.   if (!file_linear_system.is_open()) {
  28.     return;
  29.     }
  30.  
  31.   int size_linear_system = 0;
  32.   file_linear_system >> size_linear_system;
  33.  
  34.   vector< double > equation(size_linear_system + 1);
  35.   for (int i = 0; i < size_linear_system; i++) {
  36.     for (int j = 0; j < size_linear_system + 1; j++) {
  37.       file_linear_system >> equation[j];      
  38.     }  
  39.     linear_system.push_back(equation);
  40.   }
  41.  
  42.   file_linear_system.close();
  43. }                                                                              
  44.  
  45. // Считается норма матрицы размера N,  M
  46. // Предполагается, что матрица непуста
  47. double MatrixNormRow(vector< vector<double> > matrix)
  48. {
  49.   int n = matrix.size();
  50.   int m = matrix[0].size();
  51.   double norm = 0;
  52.   for (int i = 0; i < n; i++) {
  53.     double sum = 0;
  54.     for (int j = 0; j < m - 1; j++) { // Заменить m - 1 на m при копипасте
  55.       sum += fabs(matrix[i][j]);
  56.     }
  57.     if (norm < sum) {
  58.       norm = sum;
  59.     }
  60.   }
  61.  
  62.   return norm;
  63. }
  64.  
  65. // Считается норма матрицы размера N,  M
  66. // Предполагается, что матрица непуста
  67. double MatrixNormColumn(vector< vector<double> > matrix)
  68. {
  69.   int n = matrix.size();
  70.   int m = matrix[0].size();
  71.   double norm = 0;
  72.   for (int j = 0; j < m - 1; j++) { // Заменить m - 1 на m при копипасте
  73.     double sum = 0;
  74.     for (int i = 0; i < n; i++) {
  75.       sum += fabs(matrix[i][j]);
  76.     }
  77.     if (norm < sum) {
  78.       norm = sum;
  79.     }
  80.   }
  81.  
  82.   return norm;
  83. }
  84.  
  85. double VectorNorm(vector<double> vec)
  86. {
  87.   double norm = 0;
  88.   for (int i = 0; i < vec.size(); i++) {
  89.     double abs_vec_i = fabs(vec[i]);
  90.     if (norm < abs_vec_i) {
  91.       norm = abs_vec_i;
  92.     }
  93.   }
  94.   return norm;
  95. }
  96.  
  97. bool MissStopCriterion(vector<double> a, vector<double> b)
  98. {
  99.   vector<double> c = a;
  100.   for (int i = 0; i < a.size(); i++) {
  101.     c[i] -= b[i];
  102.   }
  103.   if (VectorNorm(c) < epsilon) {
  104.     return false;
  105.   } else {
  106.     return true;
  107.   }
  108. }
  109.  
  110. // http://evlm.stuba.sk/~partner2/DBfiles/Numericalanalysis/Linear%20systems%20of%20equations/Jacobi_method_for_systems_of_linear_equations_EN.pdf
  111. // На вход подается система линейных уравнений в виде матрицы размером N, N + 1.
  112. // В части матрицы размером N, N заданы коеффициенты системы. В поледнем
  113. // столбце - свободные члены. Функция возвращает решение системы с помощью
  114. // метода Якоби.
  115. vector<double> JacobiMethod(vector< vector<double> > linear_system)
  116. {
  117.   int linear_system_size = linear_system.size();
  118.   vector<double> solution(linear_system_size);
  119.  
  120.   // Построение матрицы С и вектора d
  121.   for (int i = 0; i < linear_system_size; i++) {
  122.     int a_ii = linear_system[i][i];
  123.     if (a_ii != 0) {
  124.       linear_system[i][i] = 0;
  125.       linear_system[i][linear_system_size] /= a_ii;
  126.       for (int j = 0; j < linear_system_size; j++) {
  127.         if (i != j) {
  128.           linear_system[i][j] /= -a_ii;
  129.         }
  130.       }
  131.     } else {
  132.       linear_system[i][i] = 1;
  133.       for (int j = 0; j < linear_system_size; j++) {
  134.         if (i != j) {
  135.           linear_system[i][j] *= -1;
  136.         }
  137.       }
  138.     }
  139.   }
  140.  
  141.   // Проверка сходимости
  142.   bool convergence = true;
  143.   if (MatrixNormRow(linear_system) > 1) {
  144.     if (MatrixNormColumn(linear_system) > 1) {
  145.       convergence = false;
  146.     }
  147.   }
  148.  
  149.   if (!convergence) {
  150.     return solution;
  151.   }
  152.  
  153.   // Есть сходимость, запуск итерационного процеса
  154.   vector<double> cur_iteration(linear_system_size);
  155.   for (int i = 0; i < linear_system_size; i++) {
  156.     cur_iteration[i] = linear_system[i][linear_system_size];
  157.   }
  158.   vector<double> prev_iteration(linear_system_size);
  159.  
  160.   while (MissStopCriterion(cur_iteration, prev_iteration)) {
  161.     prev_iteration = cur_iteration;
  162.     cur_iteration.assign(linear_system_size, 0);
  163.     for (int i = 0; i < linear_system_size; i++) {
  164.       for (int j = 0; j < linear_system_size; j++) {
  165.         cur_iteration[i] += linear_system[i][j] * prev_iteration[j];
  166.       }
  167.       cur_iteration[i] += linear_system[i][linear_system_size];
  168.     }
  169.   }
  170.   solution = cur_iteration;
  171.  
  172.   return solution;
  173. }
  174.  
  175. void WriteLinearSystemSolution(string file_name,
  176.                                const vector<double>& linear_system_solution)
  177. {
  178.   ofstream file_linear_system_solution(file_name);
  179.   for (int i = 0; i < linear_system_solution.size(); i++) {
  180.     file_linear_system_solution << "x" << i + 1 << " = "
  181.                                 << linear_system_solution[i] << endl;
  182.   }
  183.   file_linear_system_solution.close();
  184. }
  185.  
  186. int main (void)
  187. {
  188.   vector< vector<double> > linear_system;
  189.     ReadLinearSystem("input.txt", linear_system);
  190.   vector<double> linear_system_solution = JacobiMethod(linear_system);
  191.   WriteLinearSystemSolution("output.txt", linear_system_solution);
  192.     return 0;
  193. }
RAW Paste Data