Advertisement
cd62131

Jacobi method

May 25th, 2019
572
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 1.63 KB | None | 0 0
  1. #include <array>
  2. #include <iomanip>
  3. #include <iostream>
  4. template <std::size_t Size>
  5. static int jacobi(std::array<double, Size> &x,
  6.                   const std::array<std::array<double, Size>, Size> &a,
  7.                   const std::array<double, Size> &b, const double e,
  8.                   int count) {
  9.   ++count;
  10.   std::array<double, Size> next;
  11.   bool last = true;
  12.   for (std::size_t i = 0; i < Size; ++i) {
  13.     next[i] = b[i];
  14.     for (std::size_t j = 0; j < Size; ++j) {
  15.       if (i == j) { continue; }
  16.       next[i] -= a[i][j] * x[j];
  17.     }
  18.     next[i] /= a[i][i];
  19.     if (std::abs(next[i] - x[i]) >= e) { last = false; }
  20.   }
  21.   if (last) { return count; }
  22.   x = next;
  23.   return jacobi(x, a, b, e, count);
  24. }
  25. template <class T, std::size_t Size>
  26. std::ostream &operator<<(std::ostream &out, const std::array<T, Size> &a) {
  27.   out << '[';
  28.   for (auto &&i = a.cbegin(); i != a.cend(); ++i) {
  29.     out << (i == a.cbegin() ? "" : ", ") << *i;
  30.   }
  31.   out << ']';
  32.   return out;
  33. }
  34. int main() {
  35.   const int n = 5;
  36.   std::array<std::array<double, n>, n> a = {{{9., 3., 1., 1., 1.},
  37.                                              {1., 8., 4., 1., 1.},
  38.                                              {1., 1., 4., 1., 1.},
  39.                                              {1., 1., 4., 7., 1.},
  40.                                              {1., 1., 4., 4., 8.}}};
  41.   std::array<double, n> b = {-8., 18., -12., -15., 24.};
  42.   // initialize zero vector
  43.   std::array<double, n> x = {0., 0., 0., 0., 0.};
  44.   const double e = 1.e-6;
  45.   int count = jacobi(x, a, b, e, 0);
  46.   std::cout << std::setprecision(6) << x << std::endl << count << std::endl;
  47. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement