Основні прямі та ітераційні методи рішення СЛАР в MathCAD
Як відомо, рішення систем лінійних алгебраїчних рівнянь (СЛАР) - дуже поширений на практиці тип завдань. Теорію можна почитати за посиланням, а тут наведемо основні розрахунки як для прямих (аналітичних), так і для ітераційних (наближених) методів рішення СЛАР.
Почнемо з прямих. Класичний метод оберненої матриці в MathCAD легко реалізувати за допомогою стандартної функції lsolve або ж за допомогою операції звернення матриці, код наводити не будемо через його тривіальності. Вже двоє запитали "а як все ж знайти рішення методом зворотної матриці?" :)
Ввівши дані, як на картинці нижче, під даними написати одну з формул x: = A -1 * b або x: = lsolve (A, b)
Потім ще нижче зробити x =
А ось метод Крамера запрограмуємо. Елемент вектора рішення xi в ньому виходить у вигляді дробу, знаменником якої є визначник матриці системи, а чисельником - визначник матриці Ai. отриманої з вихідної заміною i-го стовпця стовпцем вільних членів b. Для зручності будемо в усьому документі нумерувати рядки і стовпці матриць з одиниці, тобто, встановимо значення системної змінної ORIGIN: = 1. Також визначимо загальні для всіх методів матрицю і вектор правої частини системи:
Визначення тестової СЛАР
Умовою існування і єдиності рішення СЛАР у всіх випадках є умова det A ≠ 0. тобто визначник матриці A не дорівнює нулю. Також має сенс зробити перевірку отриманого рішення, вважаючи значення нев'язки. рівне нормі різниці векторів A * x (x - знайдене рішення) і b. В ідеалі невязка повинна бути рівною нулю, але через неминуче накопичення похибок операцій над числами вона виявиться дорівнює малому числу ε. відповідному похибки методу. У MathCAD скалярний оператор "модуль" з панелі інструментів калькулятора в застосуванні до різниці векторів дасть якраз значення нев'язки, перевіримо це твердження на невеликому тесті:
З урахуванням усього сказаного, реалізуємо метод Крамера і перевірку отриманого рішення:
Метод Крамера рішення СЛАР
У тілі функції det - НЕ модуль, а схожа зовні кнопка "Визначник" з панелі "Матриці"!
Класичний метод Гаусса з приведенням матриці до верхнього трикутного вигляду детально вивчається в базовому курсі вищої математики. Реалізуємо найпростішу його різновид, що вибирає провідний елемент на головній діагоналі матриці, тобто, що працює в припущенні, що значення A1,1, ≠ 0. Так як ця підпрограма "нульового" рівня, назвемо її Gauss0. а більш складну Gauss напишемо окремо.
Просте програмування методу Гаусса в MathCAD
Для зручності вектор правої частини b записаний як (n + 1) -й стовпець матриці A. таку матрицю системи називають розширеною.
Реалізація більш "повноцінного" методу Гаусса з вибором ведучого елемента (і перестановкою при необхідності рядків матриці) виконана в доданому до статті документі MathCAD, по крайней мере, систему з нулями на головній діагоналі матриці підпрограма Gauss вирішила. Її додатковий параметр - похибка ε. починаючи з якої значення | Ai, j |<ε считается равным нулю. В случае ошибки (нет решения) подпрограмма возвращает вектор из n значений "бесконечность".
На практиці неважко побачити загальні для всіх прямих методів недоліки підходу - трудомісткість обчислень, що вимагає брати зворотні матриці або вважати визначники, що випливає з неї досить швидке накопичення похибки, нарешті, неможливість знайти рішення з заздалегідь заданої. а не закладеної в алгоритм точністю.
До певної міри уникнути цих недоліків дозволяють ітераційні методи, послідовно наближають рішення формулами виду xi (k + 1) = f (xi (k)). де k = 0,1. - номер кроку, до тих пір, поки обрана міра різниці між двома сусідніми векторами прібліженіямй | x (k + 1) -x (k) | не стане менше заданого малого значення ε. У найпростішому випадку рішення СЛАР з матрицею розмірності 2 * 2 методом Якобі (він же - метод простих ітерацій) буде виглядати так:
Всі невідомі значення xi присутні і в лівій, і в правій частинах нових рівнянь. Вибравши певний вектор початкового наближення x (0). порахуємо по ньому нове наближення x (1). потім підставимо його в праві частини рівнянь і порахуємо x (2) і т.д. до виконання умови збіжності. А воно, до речі, досить просто - метод Якобі сходиться, якщо матриця системи має діагональне переважання. тобто, на головній діагоналі знаходяться найбільші в своїх рядках елементи. Наша тестова матриця вже має діагональне переважання, а в більшості інших випадків цього можна домогтися, виконуючи перетворення над рівняннями системи, подібні до тих, що робить розширена процедура Gauss.
Вибір вектора початкового наближення x (0) на практиці також зазвичай простий, приймають x (0) = b. тобто, вектору правої частини системи. Можна і просто "занулити" вектор x (0).
Приходимо до наступної процедури вирішення:
Рішення СЛАР методом простих ітерацій (Якобі)
Зверніть увагу, що нам довелося "схитрувати" при розрахунку сум s1 і s2 - MathCAD просто не зможе обчислити суму з нижньою межею підсумовування = 1 і верхнім = 0 (або нижнім n і верхнім n-1). З тієї ж причини додаткові перевірки зроблені і в процедурі Gauss.
В основній "нескінченний" цикл підпрограми має сенс додати аварійний вихід оператором break. наприклад, по виконанні 10000 кроків.
Також, в цьому і наступному методі в рядку з break точніше був би критерій виходу | max (x1-x0) | ≤ε. де | | - значок модуля числа з панелі калькулятора.
Ітераційний метод Гаусса-Зейделя відрізняється від методу простих ітерацій лише тим, що для підрахунку i -й компоненти (k + 1) -го наближення до шуканого вектора рішення використовуються вже обчислені на цьому. тобто (K + 1) -му кроці нові значення перших i-1 компонент. а не просто береться вектор x0 цілком з попереднього кроку. В нашій подпрограмме досить замінити x0 на x1 в операторі розрахунку суми s1 :) У мене точність рішення на використаному тесті виросла при цьому вчетверо.
P.S. І ще про прямі "гауссоподобние" методи розв'язування СЛАР. Якщо Вам потрібно не покрокове програмування, а достатньо застосування стандартних функцій, є спосіб простіше. За допомогою стандартної функції augment можна отримати розширену матрицю системи (поставивши "поруч" матрицю A і вектор b), а за допомогою rref привести матрицю до ступінчастого вигляду з одиничним базисним мінор. Потім залишиться витягти рішення за допомогою методу submatrix (останній стовпець матриці, яку повернув метод rref).
Ще один спосіб вирішення СЛАР методом Гауса в MathCAD