Программирование

  // обнуляя элементы, стоящие ниже диагонали:

  for (Index j = 0; j

    const double pivot = A(j, j);

    if (pivot == 0) throw Elim_failure(j);

    // обнуляем элементы, стоящие ниже диагонали в строке i

    for (Index i = j+1; i

      const double mult = A(i, j) / pivot;

      A[i].slice(j) = scale_and_add(A[j].slice(j),

      –mult, A[i].slice(j));

      b(i) –= mult * b(j); // изменяем вектор b

    }

  }

}

Опорным называется элемент, лежащий на диагонали в строке, которую мы в данный момент обрабатываем. Он должен быть ненулевым, потому что нам придется на него делить; если он равен нулю, то генерируется исключение.

Vector back_substitution(const Matrix& A, const Vector& b)

{

  const Index n = A.dim1();

  Vector x(n);

  for (Index i = n – 1; i >= 0; ––i) {

    double s = b(i)–dot_product(A[i].slice(i+1),x.slice(i+1));

    if (double m = A(i, i))

      x(i) = s / m;

    else

      throw Back_subst_failure(i);

  }

  return x;

} 

24.6.2. Выбор ведущего элемента

Для того чтобы избежать проблем с нулевыми диагональными элементами и повысить устойчивость алгоритма, можно переставить строки так, чтобы нули и малые величины на диагонали не стояли. Говоря “повысить устойчивость”, мы имеем в виду понижение чувствительности к ошибкам округления. Однако по мере выполнения алгоритма элементы матрицы будут изменяться, поэтому перестановку строк приходится делать постоянно (иначе говоря, мы не можем лишь один раз переупорядочить матрицу, а затем применить классический алгоритм).

void elim_with_partial_pivot(Matrix& A, Vector& b)

{

  const Index n = A.dim1();

  for (Index j = 0; j < n; ++j) {

    Index pivot_row = j;

    // ищем подходящий опорный элемент:

    for (Index k = j + 1; k < n; ++k)

      if (abs(A(k, j)) > abs(A(pivot_row, j))) pivot_row = k;

    // переставляем строки, если найдется лучший опорный

    // элемент

    if (pivot_row != j) {

      A.swap_rows(j, pivot_row);

      std::swap(b(j), b(pivot_row));

    }

    // исключение:

    for (Index i = j + 1; i < n; ++i) {

      const double pivot = A(j, j);

      if (pivot==0) error("Решения нет: pivot==0");

        onst double mult = A(i, j)/pivot;

      A[i].slice(j) = scale_and_add(A[j].slice(j),

      –mult, A[i].slice(j));

      b(i) –= mult * b(j);

    }

  }

}

Для того чтобы не писать циклы явно и привести код в более традиционный вид, мы используем функции swap_rows() и scale_and_multiply().

24.6.3. Тестирование

Очевидно, что мы должны протестировать нашу программу. К счастью, это сделать несложно.

void solve_random_system(Index n)

{

  Matrix A = random_matrix(n); // см. раздел 24.7

  Vector b = random_vector(n);

  cout << "A = " << A << endl;

  cout << "b = " << b << endl;

  try {

    Vector x = classical_gaussian_elimination(A, b);

    cout << "Решение методом Гаусса x = " << x << endl;

    Vector v = A * x;

    cout << " A * x = " << v << endl;

  }

  catch(const exception& e) {

    cerr << e.what() << std::endl;

  }

}

Существуют три причины, из-за которых можно попасть в раздел catch.

• Ошибка в программе (однако, будучи оптимистами, будем считать, что этого никогда не произойдет).

• Входные данные, приводящие к краху алгоритма classical_elimination (целесообразно использовать функцию elim_with_partial_pivot).

• Ошибки округления.

Тем не менее наш тест не настолько реалистичен, как мы думали, поскольку случайные матрицы вряд ли вызовут проблемы с алгоритмом classical_elimination.