Основные прямые и итерационные методы решения СЛАУ в 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 :) У меня точность решения на использованном тесте выросла при этом вчетверо.

Метод Якоби является вариантом метода простых итераций, в котором используемые в итерационной процедуре матрица C и вектор β определяются по формулам

Cij = -Aij / Aii, если i≠j
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 [63863 просмотра]