Формирование матрицы

Учитывая сказанное, создадим программный модуль, который позволяет проверить наши возможности управления последовательностью valarray на примере задачи, близкой к реальности. Самым сложным моментом в реализуемом плане является создание функции f (), с помощью которой генерируется матрица заданной структуры, но произвольной размерности п. При генерации она помещается в последовательность типа 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 меньше. Я думаю, что вы, дорогой читатель, сами найдете достойное применение сечениям.