|
Учитывая сказанное, создадим программный модуль,
который позволяет проверить наши возможности
управления последовательностью valarray на
примере задачи, близкой к реальности. Самым
сложным моментом в реализуемом плане является
создание функции f (), с помощью которой
генерируется матрица заданной структуры, но
произвольной размерности п. При генерации она
помещается в последовательность типа valarray.
Вторая функция (f и) проста. С ее помощью
вычисляются коэффициенты уже известного вектора
решений1:
//====== Глобально заданная размерность системы
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";
}
Так как мы знаем структуру вектора правых
частей, то, изменяя граничные условия и порядок
системы, можем проверить достигнутую технику
работы с сечениями.
Тестирование обнаруживает появление численных
погрешностей (в пределах Ю"15),
обусловленных ограниченностью разрядной сетки, в
случаях когда диапазон изменения искомой
величины не кратен размеру расчетной области.
Стоит отметить, что сечения хоть и являются
непривычным инструментом, для которого хочется
найти наилучшее применение, но в рамках нашей
задачи вполне можно обойтись и без него.
Например, алгоритм умножения матрицы на вектор
можно выполнить таким образом:
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 меньше. Я думаю,
что вы, дорогой читатель, сами найдете достойное
применение сечениям. |
