Поскольку мы видим трехмерный мир, при моделировании чаще используются трехмерные объекты класса
Matrix
(например, в физическом моделировании в декартовой системе координат).
int grid_nx; // разрешение сетки; задается вначале
int grid_ny;
int grid_nz;
Matrix<double,3> cube(grid_nx, grid_ny, grid_nz);
Если добавить время в качестве четвертого измерения, то получим четырехмерное пространство, в котором необходимы четырехмерные объекты класса
Matrix
. И так далее.
24.6. Пример: решение систем линейных уравнений
Если вы знаете, какие математические вычисления выражает программа для численных расчетов, то она имеет смысл, а если нет, то код кажется бессмысленным. Если вы знаете основы линейной алгебры, то приведенный ниже пример покажется вам простым; если же нет, то просто полюбуйтесь, как решение из учебника воплощается в программе с минимальной перефразировкой.
Данный пример выбран для того, чтобы продемонстрировать реалистичное и важное использование класса
Matrix
. Мы решим систему линейных уравнений следующего вида:
a1,1x1 + ... + a1,nxn = b1
...
an,1x
1
+ ... + an,nxn = bnгде буквы x обозначают n неизвестных, а буквы a и b — константы. Для простоты предполагаем, что неизвестные и константы являются числами с плавающей точкой.
Наша цель — найти неизвестные, которые одновременно удовлетворяют указанные n уравнений. Эти уравнения можно компактно выразить с помощью матрицы и двух векторов.
Ax = b
где A — квадратная матрица n×n коэффициентов:
Векторы x и b векторы неизвестных и константа соответственно.
В зависимости от матрицы A и вектора b эта система может не иметь ни одного решения, одно решение или бесконечно много решений. Существует много разных методов решения линейных систем. Мы используем классическую схему, которая называется исключением Гаусса. Сначала мы преобразовываем матрицу A и вектор b, так что матрица А становится верхней треугольной, т.е. все элементы ниже диагонали равны нулю. Иначе говоря, система выглядит так.
Алгоритм несложен. Для того чтобы элемент в позиции (i,j) стал равным нулю, необходимо умножить строку i на константу, чтобы элемент в позиции (i,j) стал равным другому элементу в столбце j, например a(k, j). После этого просто вычтем одно уравнение из другого и получим a(i,j)==0. При этом все остальные значения в строке i изменятся соответственно.
Если все диагональные элементы окажутся ненулевыми, то система имеет единственное решение, которое можно найти в ходе обратной подстановки. Сначала решим последнее уравнение (это просто).
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
}
}
}
Опорным называется элемент, лежащий на диагонали в строке, которую мы в данный момент обрабатываем. Он должен быть ненулевым, потому что нам придется на него делить; если он равен нулю, то генерируется исключение.
