Розробка алгоритму та програми чисельного розв'язку систем лінійних алгебраїчних рівнянь з розрідженою матрицею

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

Основна ідея методу скінчених елементів полягає в тому, що будь-яку неперервну величину, таку як температура, тиск і переміщення, можна апроксимувати дискретною моделлю, яка будується на множині кусочно-неперервних функцій, визначених на кінцевім числі підобластей. Кусочно-неперервні функції визначаються за допомогою значень неперервної величини в кінцевім числі точок розглядаємої області [1,2,3].

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

а) у розглядаємій області фіксується кінцеве число точок. Ці точки називаються вузловими точками або просто вузлами;

б) значення неперервної величини в кожній вузловій точці вважається змінною, яка повинна бути визначена;

в) область визначення неперервної величини розбивається на кінцеве число підобластей, називаних елементами. Ці елементи мають загальні вузлові точки і у сукупності апроксимують форму області;

г) неперервна величина апроксимується на кожному елементі функцією, яка визначається за допомогою вузлових значень цієї величини. Для кожного елемента визначається своя функція, але функції підбираються таким чином, щоб зберігалася неперервність величини уздовж границь елемента.

Способи розв'язку СЛАР, застосовувані разом із МСЕ, можна умовно розбити на дві групи: прямі (точні) і ітераційні (наближені). Прийнято вважати, що прямі методи найбільш ефективні при рішенні СЛАР невисоких порядків (не більше 10⁵ невідомих), а ітераційні - для розв'язку СЛАР великих і надвеликих систем (понад 10⁵ невідомих).

Якщо всю матрицю коефіцієнтів СЛАР вдається розмістити в пам'яті, то розв'язок такої системи не представляє якої-небудь складності. Однак на практиці часто доводиться мати справу з матрицями, розміри яких суттєво перевищують об'єм оперативної пам'яті комп'ютера. У таких випадках і доводиться використовувати різні модифікації як прямих, так і ітераційних методів.

Остаточне рішення про застосування прямих або ітераційних методів розв'язку СЛАР необхідно ухвалювати на основі аналізу структури досліджуваної математичної задачі. Прямі методи розв'язку СЛАР вигідніше використовувати, якщо необхідно розв'язувати багато однакових систем з різними правими частинами, або якщо матриця А не є додатньо-визначеною. Крім того, існують задачі з такою структурою матриці, для якої прямі методи завжди кращі, ніж ітераційні.

1.1 Точні методи розв'язку СЛАР

Розглянемо ряд точних методів розв'язку СЛАР [4,5].

1.1.1 Метод Гауса

Нехай дана система , де A - матриця розмірності m?m (квадратна). У припущенні, що a₁₁?0, перше рівняння системи

ділимо на коефіцієнт a₁₁, у результаті одержуємо рівняння:

Потім з кожного з інших рівнянь віднімається перше рівняння, помножене на відповідний коефіцієнт a_i1. У результаті ці рівняння перетворяться до виду:

Перше невідоме виявилося виключеним із усіх рівнянь, крім першого. Далі в припущенні, що , ділимо друге рівняння на коефіцієнт і виключаємо невідоме з усіх рівнянь, починаючи з третього і т.д. У результаті послідовного виключення невідомих, матриця системи рівнянь стане трикутною:

Сукупність проведених обчислень називається прямим ходом методу Гауса.

З m-го рівняння системи визначаємо x_m, з (m-1)-го рівняння визначаємо x_m-1 і т.д. до x₁. Сукупність таких обчислень називають зворотнім ходом методу Гауса.

Реалізація прямого методу Гауса вимагає N~2m²/3 арифметичних операцій, а зворотнього - N~m² арифметичних операцій.

1.1.2 Метод Крамера

Нехай дана система лінійних рівнянь, в якій число рівнянь дорівнює числу невідомих:

Припустимо, що визначник системи d?0. Якщо тепер замінити послідовно у визначнику стовпці коефіцієнтів при невідомих x_j стовпцем вільних членів b_i, то вийдуть відповідно n визначників d₁,...,d_n.

Теорема Крамера. Система n лінійних рівнянь з n невідомими, визначник якої відмінний від нуля, завжди сумісна і має єдиний розв'язок, який обчислюється по формулах:

Розв'язок довільних систем лінійних рівнянь. Нехай

довільна система лінійних рівнянь, де число m рівнянь системи менше числа n невідомих. Тоді в матриці СЛАР знайдуться r?m лінійно незалежних рядків, а інші m-r рядків виявляться їхніми лінійними комбінаціями. Перестановкою рівнянь можна добитися того, що ці r лінійно незалежних рядків займуть перші r місць.

Звідси випливає, що кожне з останніх m-r рівнянь системи (1.7) можна представити як суму перших r рівнянь (які називаються лінійно незалежними або базисними), узятих з деякими коефіцієнтами. Тоді система еквівалентна наступній системі r рівнянь з n невідомими:

Припустимо, що мінор r-ї розмірності, складений з коефіцієнтів при перших r невідомих, відмінний від нуля: М_r?0, тобто є базисним мінором. У цьому випадку невідомі, коефіцієнти при яких складають базисний мінор, називаються базисними невідомими, а інші n-r - вільними невідомими.

У кожному з рівнянь системи (1.8) перенесемо в праву частину всі члени з вільними невідомими x_r+1,...,x_n. Тоді одержимо систему, яка містить r рівнянь з r базисними невідомими. Оскільки визначник цієї системи є базисний мінор M_r, то система має єдиний розв'язок щодо базисних невідомих, який можна знайти по формулах Крамера. Даючи вільним невідомим довільні числові значення, одержимо загальний розв'язок вихідної системи.

Однорідна система лінійних рівнянь. Нехай дана однорідна система лінійних рівнянь з n невідомими:

Оскільки додавання стовпця з нулів не змінює рангу матриці системи, то на підставі теореми Кронекера-Kaпeллі ця система завжди сумісна і має, принаймні, нульовий розв'язок. Якщо визначник системи відмінний від нуля і число рівнянь системи дорівнює числу невідомих, то по теоремі Крамера нульовий розв'язок є єдиним.

У тому випадку, коли ранг матриці системи менше числа невідомих, дана система крім нульового розв'язку буде мати і ненульові розв'язки. Для знаходження цих розв'язків у системі (1.9) виділяємо r<n лінійно незалежних рівнянь, інші відкидаємо. У виділених рівняннях у лівій частині залишаємо r базисних невідомих, а інші n-r вільних невідомих переносимо в праву частину. Тоді приходимо до системи, розв'язуючи яку по формулах Крамера, виразимо r базисних невідомих x₁,...,х_r через n-r вільних невідомих.

1.1.3 Метод головних елементів

Нехай дана система n лінійних рівнянь з n невідомими:

де елементи a_ij (i,j=1,…,n) утворюють розширену матрицю системи . Виберемо найбільший по модулю і не належачий стовпцю вільних членів елемент a_pq матриці , який називається головним елементом, і обчислимо множники m_i=-a_iq /a_pq для всіх рядків з номерами i?p (р-й рядок, що містить головний елемент, називається головним рядком).

Далі до кожного другорядного i-го рядку додамо головний рядок, помножений на відповідний множник m_i.

У результаті одержимо нову матрицю, усі елементи q-го стовпця якої, крім a_pq, складаються з нулів.

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

Таким чином, побудуємо послідовність матриць. Для визначення невідомих x_j поєднуємо в систему всі головні рядки, починаючи з останнього.

Викладений метод розв'язку систем лінійних рівнянь називається методом головних елементів. Необхідна умова його застосування полягає в тому, що визначник матриці не дорівнює нулю [6,7].

1.1.4 Схема Халецького

Нехай система лінійних алгебраїчних рівнянь дана в матричному вигляді: , де А - квадратна матриця розмірності n; , - вектори-стовпці.

Представимо матрицю А у вигляді добутку нижньої трикутної матриці С і верхньої трикутної матриці В з одиничною діагоналлю, тобто А=СВ, де

Причому елементи с_ij і b_ij визначаються по формулах:

Рівняння можна записати в наступному вигляді:

Добуток матриці B на вектор-стовпець є вектором-стовпцем, який позначимо через :

Тоді рівняння (1.13) перепишемо у вигляді:

Тут елементи с_ij відомі, тому що матриця А системи вважається вже розкладеною на добуток двох трикутних матриць С і В.

Перемноживши матриці в лівій частині рівності (1.15), одержуємо систему рівнянь, з якої одержимо наступні формули для визначення невідомих:

Невідомі y_i зручно обчислювати разом з елементами b_ij.

Після того, як усі y_i визначені по формулах (1.16), підставляємо їх у рівняння (1.14) [8].

Оскільки коефіцієнти b_ij визначені в (1.12), то значення невідомих, починаючи з останнього, обчислюємо по наступних формулах:

1.1.5 Метод квадратного кореня

Наближені методи розв'язку систем лінійних рівнянь дозволяють одержувати значення коренів системи із заданою точністю у вигляді границі послідовності деяких векторів. Процес побудови такої послідовності називається ітераційним (повторюваним).

Ефективність застосування наближених методів залежить від вибору початкового вектора і швидкості збіжності процесу.

1.2.1 Метод простих ітерацій

Нехай дана система n лінійних рівнянь з n невідомими: . Припускаючи, що діагональні елементи a_ii?0 (i=1,...,n), виразимо x₁ через перше рівняння системи, x₂ - через друге рівняння і т.д. У результаті одержимо систему:

Позначимо b_i /a_ij=?_i-a_ij /a_ii=?_ij, де i,j=1,...,n. Тоді система запишеться, таким чином, у матричній формі: x=?+?x. Розв'яжемо цю систему методом простих ітерацій. За нульове наближення приймемо стовпець вільних членів. Кожне (k+1)-е наближення обчислюють по формулі:

Якщо послідовність наближень x⁽⁰⁾,...,x^(k) має границю , то ця границя є розв'язком системи, оскільки в силу властивості границі , тобто x=?+?x [4,6].

1.2.2 Метод Зейделя

Метод Зейделя являє собою модифікацію методу простих ітерацій.

Нехай дана СЛАР, приведена до нормального вигляду: x=?+?x. Вибираємо довільно початкові наближення невідомих і підставляємо в перше рівняння системи. Отримане перше наближення підставляємо в друге рівняння системи і так далі до останнього рівняння. Аналогічно будуємо другі, треті і т.д. ітерації.

Таким чином, припускаючи, що k-і наближення відомі, методом Зейделя будуємо (k+1)-і наближення по наступних формулах:

1.2.3 Метод релаксації

Одним з важливих достоїнств МСЕ є те, що в результаті його застосування, як правило, виникають позитивно визначені, симетричні і добре обумовлені матриці коефіцієнтів. Крім того, вони звичайно мають таку важливу властивість, як низька щільність, що при правильній нумерації вузлів кінцево-елементної розбивки тіла приводить до групування ненульових елементів стрічкою певної ширини уздовж головної діагоналі. Це і пояснює той факт, що більшість методів розв'язку СЛАР, застосовуваних разом із МСЕ, орієнтована на розв'язок систем зі стрічковою структурою матриці коефіцієнтів.

Симетричність матриці коефіцієнтів дозволяє не запам'ятовувати майже половину її елементів, а при рішенні СЛАР зменшити об'єм обчислень і, як наслідок, знизити кількість помилок округлення.

На жаль, існують цілі класи задач, для яких ширина стрічки прагне до її порядку. Наприклад, при дослідженні за допомогою МСЕ контактних задач механіки деформованого твердого тіла, у яких у взаємодії беруть участь більше двох тіл, ширина стрічки дорівнює розмірності глобальної матриці жорсткості всіх контактуючих тіл. Тому при більших порядках СЛАР ітераційні методи розв'язку часто виявляються кращими.

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

Усі ці факти, а також простота програмної реалізації більшості ітераційних методів розв'язку СЛАР і обумовлюють їхню популярність.

Застосування ітераційних методів дозволяє суттєво скоротити витрати оперативної пам'яті комп'ютера за рахунок того, що при формуванні глобальної матриці жорсткості досить фіксувати тільки її ненульові елементи. У зв'язку з цим при використанні кінцево-елементної технології виникає проблема розробки ефективних алгоритмів формування, зберігання і використання розріджених матриць.

Пам'ять, використовувана для зберігання розріджених матриць, складається з двох частин: основної пам'яті, у якій утримуються числові значення елементів матриць, і додаткової пам'яті, де зберігаються вказівники, індекси і інша інформація, необхідна для формування структури матриць і яка забезпечує доступ до числових значень їх елементів при виконанні процедур формування і розв'язку СЛАР, тобто так звані списки зв'язності. Способи зберігання і використання даних, що зберігаються в основній і додатковій пам'яті, досить різноманітні і визначаються, головним чином, обраним методом розв'язку СЛАР. Для структурної побудови інформації, що зберігається в додатковій пам'яті - списків зв'язності - широко використовується теорія графів.

Особливо складні алгоритми роботи з розрідженими матрицями виникають при використанні прямих методів розв'язку. У першу чергу це пов'язане з тим, що при розкладанні розрідженої матриці A=LL^T вона звичайно зазнає деякого заповнення, тобто нижній трикутний множник L має ненульові елементи в тих позиціях, де у вихідній матриці A стояли нулі (не виключено, що при цьому можуть з'явитися і нові нульові елементи). Крім того, прямі методи вимагають спеціальної нумерації вузлів кінцево-елементної сітки, яка забезпечує мінімальну ширину стрічки. Подібних проблем при використанні ітераційних методів не виникає [10,11].

У рамках даної роботи стосовно ітераційних методів розв'язку СЛАР викладені 2 алгоритми зберігання і використання розріджених матриць, які відрізняються простотою реалізації і високою обчислювальною ефективністю.

2.1 Перша схема

Матриця жорсткості, що виходить при застосуванні МСЕ, має симетричну структуру, що дозволяє в загальному випадку зберігати тільки верхню трикутну частину матриці. Однак для задач із великою кількістю невідомих це також приводить до проблеми нестачі пам'яті. Викладений у даній роботі метод дозволяє зберігати тільки ненульові члени матриці жорсткості. Суть його полягає в наступному.

49

Спочатку, з метою виявлення зв'язків кожного вузла з іншими, проводиться аналіз структури дискретизації області на КЕ. Наприклад, для КЕ-сітки, зображеної на рис. 2.1, відповідна структура зв'язків буде мати вигляд, зображений у табл. 2.1.

Таблиця 2.1 - Структура зв'язків КЕ-сітки

Тоді для зберігання матриці жорсткості необхідно построчно запам'ятовувати інформацію про коефіцієнти, відповідні до вузлів, з якими зв'язаний даний вузол. На рис. 2.2 наведені матриця жорсткості і її компактний вигляд для сітки, зображеної на рис. 2.1 [9].

Текст програми, що реалізує викладений алгоритм аналізу структури КЕ-розбивки тіла, наведений у додатку А [12,13].

Даний спосіб компактного зберігання матриці жорсткості дозволяє легко його використовувати разом з яким-небудь чисельним методом. Найбільш зручним для цієї мети здається використання вищевикладеного ітераційного методу Ланцоша, тому що на кожній ітерації потрібно тільки перемножувати матрицю коефіцієнтів СЛАР і заданий вектор. Отже, для використання запропонованого методу компактного зберігання матриці СЛАР необхідно побудувати пряме і зворотнє перетворення в первісну квадратну матрицю.

Нехай a_ij - елемент первісної квадратної матриці розмірністю n?n, а - її компактний вигляд. Тоді для зворотнього перетворення будуть слушні наступні співвідношення:

де m - кількість ступенів волі (m=1,2,3).

Для прямого перетворення будуть слушні співвідношення, зворотні до співвідношень (2.1) - (2.4).

2.2 Друга схема

При реалізації ітераційного розв'язку системи рівнянь досить частою є ситуація, коли необхідно виконати множення матриці системи A на який-небудь вектор, наприклад, на вектор вузлових невідомих або на вектор нев'язки , де [20]. Побудова добутків типу або є одним з вузьких місць ефективної реалізації всього ітераційного процесу розв'язку системи рівнянь, оскільки вимагає найбільших витрат процесорного часу. Нижче представлений алгоритм побудови добутку і безпосередньо пов'язаний з ним спосіб компактного зберігання матриці A [14-16].

Рисунок 2.3 - Кінцево-елементна модель

Розглянемо як приклад кінцево-елементну модель, зображену на рис. 2.3. З кожним вузлом сітки зв'язане деяке число вузлів, які можуть бути розбиті на дві множини: першу множину становлять вузли, номери яких менше i - номера розглядаємого вузла, а друга множина - вузли, номери яких більше i. У табл. 2.2 зображена структура матриці СЛАР, відповідна до вузлових зв'язків кінцево-елементної моделі, представленої на рис. 2.3.

Таблиця 2.2 - Структура матриці СЛАР

Матриця A зберігається у вигляді двох векторів (рис. 2.4): - члени головної діагоналі, - ненульові елементи, розташовані над головною діагоналлю. Для формування і використання елементів вектора створюються додаткові чотири вектори, що містять інформацію про зв'язки вузлів кінцево-елементної моделі.

Рисунок 2.4 - Структура векторного зберігання матриці A

Розглянемо i-е рівняння СЛАР . У загальному випадку воно має вигляд:

де a_ij - елементи матриці A; f_i - елемент вектора правої частини ; x_j - елемент вектора невідомих ; N - число невідомих.

У рівнянні (2.5) присутні суми добутків двох типів. Перша сума складається з добутків елементів a_ij - матриці A і елементів вектора невідомих , у яких j<i, а для добутків другої суми маємо j>i. Для побудови добутків першої суми використовуються вектори і , приклад структури яких показаний на рис. 2.5. У векторі зберігаються номери тих вузлів, які пов'язані з даним розглядаємим вузлом, і їх номера менше номера даного вузла. Наприклад, розглянемо вузол з номером 5 (див. рис. 2.3). Цей вузол входить в елементи з номерами 2 і 3. У силу цього він має зв'язок з вузлами 4, 2 і 1, номера яких менше 5. Номера цих вузлів записуються у вектор у послідовності, визначеній номерами елементів і правилом обходу їх вузлів. На рис. 2.3 стрілками показаний початок і напрямок обходу вузлів елементів. Цим пояснюється порядок появи номерів 4, 2 і 1 у векторі . У векторі зберігаються номери індексів, що координують у векторі розташування номера першого вузла всього масиву номерів вузлів, пов'язаних з розглядаємим вузлом. Наприклад, як видно на рис. 2.3, масив вузлів, пов'язаних з вузлом 5, починається з індексу 6. Останній індекс масиву визначається початковим індексом вузла, наступного за розглянутим мінус одиницю. Наприклад, для вузла 5 останній індекс масиву у векторі рівний 8=9-1, де номер 9 визначає початок масиву вузлів для вузла 6 (див. рис. 2.3).

Зовсім аналогічно будуються ще два вектори і , що містять інформацію про номери вузлів, пов'язаних з даним розглядаємим вузлом, але номера яких більше номера даного вузла. На рис. 2.5 зображено приклад конструкції зазначених векторів стосовно до кінцево-елементної моделі, зображеної на рис. 2.3. У векторі вказується початковий індекс, що фіксує початок списку вузлів, розміщених у векторі і пов'язаних з розглядаємим вузлом.

Рисунок 2.5 - Структура векторів-вказівників

У тих рідких випадках, коли той або інший розглядаємий вузол кінцево-елементної моделі не має зв'язку з вузлами, номера яких менше або більше номера даного вузла, тоді у відповідні гнізда векторів або заносяться нулі. Наприклад, вузол 4 (див. рис. 2.3) не має зв'язку з вузлами, номера яких менше 4, це враховується записом нуля в гніздо 4 вектора (рис. 2.5). Вузол 3 не має зв'язку з вузлами, номера яких більше 3 - цій обставині відповідає запис нуля в гніздо 3 вектора . Виключення становить запис в 1-е гніздо вектора , куди заноситься номер останнього гнізда вектора . Для всіх інших вузлів ці параметри визначаються тривіально. Для кожної кінцево-елементної моделі інформаційні вектори будуються один раз - відразу ж після того, як закінчена побудова масиву, що описує кінцеві елементи сітки.

Структура вектора (див. рис. 2.4) у значній мірі повторює структуру вектора . Ненульові наддіагональні елементи a_ij матриці A зберігаються у векторі порядково, а порядок їх розміщення в межах рядка визначається порядком розміщення індексів у векторі .

Тепер розглянемо побудову сум добутків, що входять у рівняння (2.5). Розгляд почнемо з другої суми, а саме:

При фіксованому індексі i визначаються номери першого і останнього гнізд вектора , у яких розміщаються номери вузлів, пов'язаних з j-м вузлом, і для яких j>i. Наприклад, для вузла 5 (див. рис. 2.3) маємо: номер 1-го гнізда 13, а номер останнього гнізда 14=15-1. Номера першого і останнього гнізд визначають у векторі діапазон розміщення ненульових елементів матриці A, які приймають участь у добутках суми (2.6). Крім того, ці ж номера гнізд визначають у векторі діапазон розміщення номерів відповідних елементів x_i вектора . Таким чином, наприклад, для вузла 5 сума (2.6) складається з двох додатків:

Побудова першої суми з рівняння (2.5)

алгоритмічно виконується складніше. При фіксованому номері i визначаються номери першого n і останнього m (n<m) гнізд вектора , у якому зберігаються номери j вузлів, пов'язаних з вузлом i, і для яких j<i. Цим самим задаються компоненти x_i вектора . Тепер необхідно визначити елементи a_ij=a_ji (i?j) матриці A, що зберігаються у векторі . Ця процедура виконується в такий спосіб. Послідовно з гнізд із n по m у векторі беруться номери вузлів j (j<i). По цих номерах j з вектора визначаються номери першого k і останнього l (k<l) гнізд вектора . Для кожного вузла j у діапазоні гнізд (l-k) вектора перебуває номер i, оскільки, по-перше, вузол i пов'язаний з кожним вузлом j, а, по-друге, i>j. Переглядаючи гнізда вектора і порівнюючи номера вузлів, що перебувають у цих гніздах, з номером i, можна визначити номер гнізда, у якому для даного вузла j і його діапазону номерів гнізд (l-k) перебуває номер вузла i. По номеру гнізда вузла i у векторі можна знайти у векторі елемент a_ij=a_ji (i?j, j<i).

Розглянемо приклад. Побудуємо суму (2.8) для вузла 5 кінцево-елементної моделі, зображеної на рис. 2.3. Із цим вузлом, як ми вже відзначали вище, зв'язані вузли 4, 2 і 1. Номера цих вузлів зберігаються у векторі у гніздах з 4 по 6, причому їх номера менше 5. Таким чином, компоненти вектора відомі: х₄, х₂ і x₁. Визначимо елементи матриці A. Вузлу 4 у векторі відповідають гнізда з 8 по 12. У цих гніздах зберігаються номери вузлів, пов'язаних з номером 4, причому всі ці номери більше 4. Один із цих номерів - 5. Цей номер перебуває в гнізді вектора з номером 12. У гнізді з цим же номером 12 у векторі зберігається елемент a₄₅=a₅₄ (j=4<i=5) матриці А. Аналогічно відшукуються елементи a₂₅ (j=2<i=5) і a₁₅ (j=1<i=5). Таким чином, одержуємо вираз для :

Остаточно з урахуванням (2.7) уся сума S₅ має такий вигляд:

Із представленого вище алгоритму видно, що обчислювальні витрати на побудову сум добутків і різні. Друга сума будується значно швидше. Це пов'язано з тим, що при побудові першої суми необхідно проводити додатковий порівняльний аналіз номерів вузлів, розміщених у векторі . При розгляді процедур алгоритму передбачалося, що всі необхідні дані розташовуються в оперативній пам'яті комп'ютера. Довжина одномірних масивів, у яких розміщаються вектори , і , залежить від структури розглядаємої кінцево-елементної моделі. У табл. 2.3 представлені дані по розмірності цих векторів для трьох типів кінцево-елементних моделей, побудованих у результаті розбивки одиничного куба на 20-вузлові квадратичні елементи. Куб розбивався в наступних пропорціях: 5:5:5 (КЕМ № 1), 10:10:10 (КЕМ № 2) і 20:20:20 (КЕМ № 3).

Таблиця 2.3 - Розмірність векторів КЕМ

Порівняння даного алгоритму з алгоритмами, розглянутими в роботах [17-19], показує наступне. Для реалізації запропонованого алгоритму потрібно більше оперативної пам'яті, оскільки для координування ненульових елементів використовуються не два масиви (вектори номерів рядків і стовпців), а чотири. Додаткові вектори і , розмірності яких не перевищують числа діагональних елементів, дозволяють локалізувати зони пошуку елементів-співмножників при побудові сум (2.6) і (2.8). За рахунок цього досягається економія часових витрат на розв'язок СЛАР.

На закінчення відзначу, що викладений вище алгоритм формування сум у рівнянні (2.5) може бути повністю застосований і для побудови інших аналогічних матричних добутків, необхідних для реалізації ітераційних процесів.

3 Оптимізація обчислень