Шрифт:
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);
return v;
}
И вновь простая операция над объектом класса
Matrix
делает за нас большую часть работы. Как указывалось в разделе 24.5.3, операции вывода объектов класса Matrix
описаны в заголовке MatrixIO.h
. Функции random_matrix
и random_vector
просто используют случайные числа (раздел 24.7). Читатели могут написать эти функции в качестве упражнения. Имя Index
является синонимом типа индекса, используемого в библиотеке Matrix
, и определено с помощью оператора typedef
(раздел A.15). Этот тип включается в программу с помощью объявления using
.
using Numeric_lib::Index;
24.7. Случайные числа
<cstdlib>
из стандартной библиотеки есть такой код:
int rand; // возвращает числа из диапазона
// [0:RAND_MAX]
RAND_MAX // наибольшее число, которое генерирует
// датчик rand
void srand(unsigned int); // начальное значение датчика
// случайных чисел
Повторяющиеся вызовы функции
rand
генерируют последовательность чисел типа int
, равномерно распределенных в диапазоне [0:RAND_MAX]
. Эта последовательность чисел называется псевдослучайной, потому что она генерируется с помощью математической формулы и с определенного места начинает повторяться (т.е. становится предсказуемой и не случайной). В частности, если мы много раз вызовем функцию rand
в программе, то при каждом запуске программы получим одинаковые последовательности. Это чрезвычайно полезно для отладки. Если же мы хотим получать разные последовательности, то должны вызывать функцию srand
с разными значениями. При каждом новом аргументе функции srand
функция rand
будет порождать разные последовательности. Например, рассмотрим функцию
random_vector
, упомянутую в разделе 24.6.3. Вызов функции random_vector(n)
порождает объект класса Matrix<double,1>
, содержащий n
элементов, представляющих собой случайные числа в диапазоне от [0:n]
:
Vector random_vector(Index n)
{
Vector v(n);
for (Index i = 0; i < n; ++i)
v(i) = 1.0 * n * rand / RAND_MAX;