Основные прямые и итерационные методы решения СЛАУ в 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
оператор |A|
- это не модуль числа с панели "Калькулятор", а похожая внешне кнопка "Определитель" с панели "Матрицы"!
Классический метод Гаусса с приведением матрицы к верхнему треугольному виду подробно изучается в базовом курсе высшей математики. Реализуем самую простую его разновидность, выбирающую ведущий элемент на главной диагонали матрицы, то есть, работающую в предположении, что значение 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 методом простых итераций будет выглядеть так:
Из первого уравнения выразим неизвестное x1
, а из второго - x2
, получаем систему:
x1 = (b1 - A1,2*x2) / A1,1
x2 = (b2 - A2,1*x1) / A2,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
:) У меня точность решения на использованном тесте выросла при этом вчетверо.
Метод Якоби является вариантом метода простых итераций, в котором используемые в итерационной процедуре матрица Cij = -Aij / Aii, если i≠j
C
и вектор β
определяются по формулам
Cij = 0, если i=j
βi = bi / Aii, i =1,2,...,n
Вот расчёт методом Якоби с критерием выхода "максимальный модуль разности между проекциями
x1(k)i
и x0(k)i
стал меньше либо равен заданной точности ε
:
Метод Якоби
При расчёте r
применены операторы "Векторизовать" с панели "Матрицы" и "Модуль" с Калькулятора, а при расчёте
элементов С
и β
деление выполнено "в строчку" для экономии места.
Прямые и итерационные методы решения СЛАУ: скачать документ .xmcd Mathcad 14/15 в архиве .zip (78 Кб)
P.S. И ещё про прямые "гауссоподобные" методы решения СЛАУ. Если Вам нужно не пошаговое программирование, а достаточно применения стандартных функций, есть способ проще.
С помощью стандартной функции augment
можно получить расширенную матрицу системы (поставив "рядом" матрицу A
и вектор b
), а с помощью rref
привести матрицу к ступенчатому виду с единичным базисным минором. Потом останется извлечь решение с помощью метода submatrix
(последний столбец матрицы, которую вернул метод rref
).
Норма вектора |A*x-b|
, как и другие нормы в статье, берётся кнопкой |x|
с Калькулятора, а не похожей на неё кнопкой с панели "Матрицы".
Еще один способ решения СЛАУ методом Гаусса в MathCAD
21.09.2013, 11:24 [64374 просмотра]