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