13.7. Операції з матрицями

Матричні обчислення відіграють значну роль при розв’язанні прикладних задач. Вони покладені в основу математичного апарата квантової й статистичної механіки й фізики, радіоелектроніки, теорії коливань, де фундаментальне значення мають характеристичні рівняння, власні числа й вектори. Алгоритми лінійної алгебри є складовою частиною алгоритмів розв’язування дуже великої кількості лінійних і нелінійних задач моделювання і обробки даних, задач оптимізації, багатьох інших задач у наукових дослідженнях, технічних розробках, розв’язуванні економічних проблем. У даному підрозділі розглядаються декілька алгоритмів лінійної алгебри, які використовуються найбільш часто.

Звертаємо увагу на те, що у загально-математичних формулюваннях нумерація елементів масивів починається з 1, в той час як у програмних реалізаціях, які наводяться, нумерація, згідно з традиціями мови С, починається з 0.

Задача обчислення сідлової точки матриці

Задана прямокутна матриця А розміром NM. Потрібно обчислити її мінімаксну сідлову точку за формулою :

, (1)

або максимінну сідлову точку:

. (2)

Алгоритм

Алгоритм випливає прямо з формул (1) і (2).

Нижче наведені дві функції, які повертають значення сідлових точок матриці з елементами типу double.

// Приклад 1

double MinMax (double** A, int N, int M )

{ int i,j; double minmax;

double* D = new double[N]; // створення допоміжного масиву

for (i=0; i<N; i++)

{ D[i] = A[i][0];

for (j=1; j<M; j++) if (A[i][j]>D[i]) D[i] = A[i][j];

}

// знайти найменший серед найбільших

minmax = D[0];

for (i=1; i<N; i++) if (D[i] < minmax) minmax = D[i];

delete[] D;

return minmax;

}

// Приклад 2

double MaxMin (double** A, int N, int M )

{ int i,j; double maxmin;

double* D = new double[N];

for (i=0; i<N; i++)

{ D[i] = A[i][0];

for (j=1; j<M; j++) if (A[i][j]<D[i]) D[i] = A[i][j];

}

maxmin = D[0];

for (i=1; i<N; i++) if (D[i] > maxmin) maxmin = D[i];

delete [] D;

return maxmin;

}

Задача обчислення визначника матриці

Задана квадратна числова матриця:

. (3)

Обчислити її визначник det(A).

Алгоритм

До еквівалентних перетворень відносять такі перетворення, які не змінюють визначник матриці. Використовуючи певну послідовність еквівалентних перетворень, вихідну матрицю можна звести до трикутного вигляду:

. (4)

Визначник такої матриці легко обчислюється шляхом отримання добутку діагональних елементів:

det (A) = det (A') = . (5)

Якщо при виконанні перетворень матриці виконувались операції перестановки рядків, то формула для обчислення визначника дещо змінюється:

, (6)

де p - кількість операцій парної перестановки рядків.

Для того, щоб перетворити елемент матриці a_kj, що лежить нижче головної діагоналі, у нуль, треба отримати нові елементи k-того рядка за такою формулою:

,. (7)

З останньої формули видно, що діагональний елемент a_ii повинен бути відмінним від нуля. Щоб цього уникнути, і, одночасно, збільшити точність і надійність обчислення визначника, застосовується прийом, названий «вибором головного елемента». Цей прийом полягає в наступному. Припустимо в i-м стовпці необхідно перетворити в нуль елементи, що лежать нижче головної діагоналі. Серед елементів цього стовпця, починаючи з a_ii і до останнього елемента a_ni, знаходимо максимальний по модулю елемент. Потім міняємо місцями i-тий рядок і рядок, що містить максимальний по модулю елемент, місцями. Таким чином, на місці діагонального елемента завжди буде максимальний (по абсолютній величині) з можливих елементів. Якщо знайдений максимальний елемент дорівнює нулю, то обчислення можна завершити: визначник у цьому випадку буде дорівнювати нулю.

Нижче наводиться приклад функції DetMatr(A,N), яка обчислює визначник матриці A розміром NN. Обчислення виконуються відповідно до дійсного типу double.

// Приклад 3

double DetMatr(double** A, unsigned N )

{ int i, j, k=0;

double D, R;

double** C= new double*[N]; // створення робочої

for (i=0; i<N; i++) C[i]=new double[N]; // матриці С

for (i=0; i<N; i++)

for (j=0; j<N; j++) C[i][j] = A[i][j];

// копіювання матриці

for (i=0; i<N-1; i++) // основний цикл

{ // вибираємо головний елемент

D=fabs(C[i][i]);

k=i;

for (j=i+1; j<N; j++)

if (fabs(C[j][i])>D) { k = j; D = fabs(C[j][i]); }

// міняємо, якщо треба, рядки місцями

if (k!=i) swp(C[k],C[i]);

// перетворення в нуль елементів стовпця R = 1./C[i][i];

for (k=i+1; k<N; k++)

{ D = C[k][i]*R;

for (j=i; j<N; j++) C[k][j] = C[k][j]- D*C[i][j];

}

}

// отримаємо визначник як добуток елементів // головної діагоналі

D= 1.0;

for (i=0; i<N; i++) D *= C[i][i];

for (i=0; i<N; i++) delete[] C[i]; // звільнення пам’яті

delete[] C;

return D;

}

Задача обчислення оберненої матриці

Надана квадратна числова матриця А з розмірами nn. Отримати матрицю A^-1, яка відповідає такій умові:

A A^-1 = A^-1 A = I , (8)

де I - одинична діагональна матриця.

Алгоритм

Створимо робочу прямокутну матрицю С розміром n2n, приписавши праворуч до вихідної матриці A одиничну діагональну матрицю E такого ж розміру:

. (9)

Виконаємо над цією розширеною матрицею перетворення, мета яких - перетворити у нулі всі елементи, що лежать нижче головної діагоналі лівої частини матриці. Операції перетворення рядків виконуємо над рядками розширеної матриці в цілому, включаючи одиничну підматрицю. Перетворення проводяться в три етапи:

1. Перетворюємо в нулі всі елементи підматриці А, які лежать нижче головної діагоналі.

2. Перетворюємо в нулі всі елементи підматриці А, що лежать вище головної діагоналі.

3. Перетворюємо всі діагональні елементи підматриці А в одиниці, шляхом ділення кожного рядка розширеної матриці на відповідний діагональний елемент.

Після виконання цих кроків, підматриця А стане одиничною, а підматриця Е перетвориться в шукану обернену матрицю А^-1. У наведеному алгоритмі можна також застосувати вибір головного елемента, відповідна операція виконується на етапі 1.

Нижче наведено функцію ObrMatr, яка реалізує поданий вище алгоритм. Вхідна матриця зберігається в масиві А, зворотна записується в масив В, N - розмір матриці. Двовимірний масив B відповідних розмірів повинен бути створений заздалегідь.

// Приклад 4