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

просто).

an,nxn = bn

Очевидно, что 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. Тестирование

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

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.

Для того чтобы проверить наше решение, выводим на экране произведение A*x, которое должно быть равно вектору b (или достаточно близким к нему с учетом ошибок округления). Из-за вероятных ошибок округления мы не можем просто ограничиться инструкцией

if (A*x!=b) error("Неправильное решение");

Поскольку числа с десятичной точкой являются лишь приближением действительных чисел, получим лишь приближенный ответ. В принципе лучше не применять операторы == и != к результатам вычислений с десятичными точками: такие числа являются лишь приближениями.

В библиотеке Matrix нет операции умножения матрицы на вектор, поэтому эту функцию нам придется написать самостоятельно.

Vector operator*(const Matrix& m,const Vector& u)

{

  const Index n = m.dim1();

  Vector v(n);

  for (Index i = 0; i < n; ++i) v(i)