Очевидно, что мы должны протестировать нашу программу. К счастью, это сделать несложно.
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);
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. Случайные числа
Если вы попросите любого человека назвать случайное число, то они назовут 7 или 17, потому что эти числа считаются самыми случайными. Люди практически никогда не называют число нуль, так как оно кажется таким идеально круглым числом, что не воспринимается как случайное, и поэтому его считают наименее случайным числом. С математической точки зрения это полная бессмыслица: ни одно отдельно взятое число
нельзя назвать случайным. То, что мы часто называем случайными числами — это последовательность чисел, которые подчиняются определенному закону распределения и которые невозможно предсказать, зная предыдущие числа. Такие числа очень полезны при тестировании программ (они позволяют генерировать множество тестов), в играх (это один из способов гарантировать, что следующий шаг в игре не совпадет с предыдущим) и в моделировании (мы можем моделировать сущность, которая ведет себя случайно в пределах изменения своих параметров).
Как практический инструмент и математическая проблема случайные числа в настоящее время достигли настолько высокой степени сложности, что стали широко использоваться в реальных приложениях. Здесь мы лишь коснемся основ теории случайных чисел, необходимых для осуществления простого тестирования и моделирования. В заголовке
<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;
return v;
}
Обратите внимание на использование числа
1.0
, гарантирующего, что все вычисления будут выполнены в арифметике с плавающей точкой. Иначе при каждом делении целого числа на
RAND_MAX
мы получали бы
0
.
Сложнее получить целое число из заданного диапазона, например
[0:max]
. Большинство людей сразу предлагают следующее решение:
int val = rand%max;
Долгое время такой код считался совершенно неудовлетворительным, поскольку он просто отбрасывает младшие разряды случайного числа, а они, как правило, не обладают свойствами, которыми должны обладать числа, генерируемые традиционными датчиками случайных чисел. Однако в настоящее время во многих операционных системах эта проблема решена достаточно успешно, но для обеспечения переносимости своих программ мы рекомендуем все же скрывать вычисления случайных чисел в функциях.