Программирование. Принципы и практика использования C++ Исправленное издание

Алгоритм несложен. Для того чтобы элемент в позиции (i,j) стал равным нулю, необходимо умножить строку i на константу, чтобы элемент в позиции (i,j) стал равным другому элементу в столбце j, например a(k, j). После этого просто вычтем одно уравнение из другого и получим a(i,j)==0. При этом все остальные значения в строке i изменятся соответственно.

Если все диагональные элементы окажутся ненулевыми, то система имеет единственное решение, которое можно найти в ходе обратной подстановки. Сначала решим последнее уравнение (это просто).

a_n,nx_n = b_n

Очевидно, что x[n] равен b[n]/a(n,n). Теперь исключим строку n из системы, найдем значение x[n–1] и будем продолжать процесс, пока не вычислим значение x[1].

При каждом значении n выполняем деление на a(n,n), поэтому диагональные значения должны быть ненулевыми. Если это условие не выполняется, то обратная подстановка завершится неудачей. Это значит, что система либо не имеет решения, либо имеет бесконечно много решений.

24.6.1. Классическое исключение Гаусса

Посмотрим теперь, как этот алгоритм выражается в виде кода на языке С++. Во-первых, упростим обозначения, введя удобные имена для двух типов матриц, которые собираемся использовать.

 typedef Numeric_lib::Matrix<double, 2> Matrix;
 typedef Numeric_lib::Matrix<double, 1> Vector;

Затем выразим сам алгоритм.

 Vector classical_gaussian_elimination(Matrix A,Vector b)
 {
  classical_elimination(A, b);
  return back_substitution(A, b);
 }

Иначе говоря, мы создаем копии входных матрицы A и вектора b (используя механизм передачи аргументов по значению), вызываем функцию для решения системы, а затем вычисляем результат с помощью обратной подстановки. Такое разделение задачи на части и система обозначений приняты во всех учебниках. Для того чтобы завершить программу, мы должны реализовать функции classical_elimination и back_substitution. Решение также можно найти в учебнике.

 void classical_elimination(Matrix& A,Vector& b)
 {
  const Index n = A.dim1;
  // проходим от первого столбца до последнего,
  // обнуляя элементы, стоящие ниже диагонали:
  for (Index j = 0; j<n–1; ++j) {
  const double pivot = A(j, j);
  if (pivot == 0) throw Elim_failure(j);

  // обнуляем элементы, стоящие ниже диагонали в строке i
  for (Index i = j+1; i<n; ++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. Тестирование