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

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);

}

<div class="fb2-code"><code> }</code></div>

}

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

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) = dot_product(m[i], u);