Формирование матрицы
Учитывая сказанное, создадим программный модуль, который позволяет проверить наши возможности управления последовательностью valarray на примере задачи, близкой к реальности. Самым сложным моментом в реализуемом плане является создание функции f (), с помощью которой генерируется матрица заданной структуры, но произвольной размерности n. При генерации она помещается в последовательность типа valarray. Вторая функция (f и) проста. С ее помощью вычисляются коэффициенты уже известного вектора решений:
//====== Глобально заданная размерность системы
int n;
//====== Граничные условия
double UO, UN;
//====== Функция вычисления коэффициентов
//====== трехдиагональной матрицы
double f ()
{
//====== Разовые начальные установки
static int raw = -1, k = -1, col = 0;
//====== Сдвигаемся по столбцам
col++;
//====== k считает все элементы
//====== Если начинается новая строка
if (++k % n == 0)
{
col =0; // Обнуляем столбец
raw++; // Сдвигаемся по строкам
}
//====== Выделяем три диагонали
return col==raw? -2.
: col == raw-1 И col==raw+l? 1.
: 0.;
}
double fu()
{
//==== Вычисления вектора правых частей по формуле (5)
static double
dU = (UN-UO)/(n+1),
d = U0; return d += dU;
}
В функции main создается valarray с трехдиагональной матрицей и производится умножение матрицы на вектор решений. Алгоритм умножения использует сечение, которое вырезает из valarray текущую строку матрицы:
void main()
{
//======= Размерность задачи и граничные условия
n =4;
UO = 100.;
UN = 0.;
//======= Размерность valarray (вся матрица)
int nn = n*n;
//======= Матрица и два вектора
valarray<double> a(nn), u(n), v(n);
//======= Генерируем их значения
generate (&а[0], &a[nn], f); generate (&u[0], &u[n], fu);
out ("Initial matrix", a); out ("Initial vector", u);
//======= Умножение матрицы на вектор
for (int i=0; i<n; i++) {
//======= Выбираем i-ю строку матрицы
valarray<double> s = a[slice (i*n, n, 1)];
//======= Умножаем на вектор решений
//======= Ответ помещаем в вектор v <
transform(&s[0], &s[n], &u[0], &v[0], multiplies<double>());
//======= Суммируем вектор, получая
//======= i-ый компонент вектора правых частей
cout << "\nb[" << i << "] = " << v.sum(); }
cout <<"\n\n";
}
Так как мы знаем структуру вектора правых частей, то, изменяя граничные условия и порядок системы, можем проверить достигнутую технику работы с сечениями.
Тестирование обнаруживает появление численных погрешностей, обусловленных ограниченностью разрядной сетки, в случаях когда диапазон изменения искомой величины не кратен размеру расчетной области. Стоит отметить, что сечения хоть и являются непривычным инструментом, для которого хочется найти наилучшее применение, но в рамках нашей задачи вполне можно обойтись и без него. Например, алгоритм умножения матрицы на вектор можно выполнить таким образом:
for (int i=0, id=0; i<n; i++, id+*n)
{
transform(&a[id], &a[id+n], &u[0], &v[0],
multiplies<dovible> ());
cout << "\nb[" << i << "] = " << v.sum();
}
Эффективность этой реализации, несомненно, выше, так как здесь на один valarray меньше. Я думаю, что вы, дорогой читатель, сами найдете достойное применение сечениям.