Как в mathcad решить систему линейных уравнений

Zoloto585CPA

Как в mathcad решить систему линейных уравнений

БлогNot. Основные прямые и итерационные методы решения СЛАУ в MathCAD

Основные прямые и итерационные методы решения СЛАУ в MathCAD

Как известно, решение систем линейных алгебраических уравнений (СЛАУ) — весьма распространённый на практике тип задач. Теорию можно почитать по ссылке, а здесь приведём основные расчёты как для прямых (аналитических), так и для итерационных (приближённых) методов решения СЛАУ.

Начнём с прямых. Классический метод обратной матрицы в MathCAD легко реализовать с помощью стандартной функции lsolve или же посредством операции обращения матрицы, код приводить не будем из-за его тривиальности.

А вот метод Крамера запрограммируем. Элемент вектора решения xi в нём получается в виде дроби, знаменателем которой является определитель матрицы системы, а числителем – определитель матрицы Ai , полученной из исходной заменой i-го столбца столбцом свободных членов b . Для удобства будем во всём документе нумеровать строки и столбцы матриц с единицы, то есть, установим значение системной переменной ORIGIN:=1 . Также определим общие для всех методов матрицу и вектор правой части системы:

Условием существования и единственности решения СЛАУ во всех случаях является условие det A≠0 , т.е., определитель матрицы A не равен нулю. Также имеет смысл сделать проверку полученного решения, посчитав значение невязки, равное норме разности векторов A*x ( x — найденное решение) и b . В идеале невязка должна быть равной нулю, но из-за неизбежного накопления погрешностей операций над вещественными числами она окажется равна малому числу ε , соответствующему погрешности метода. В MathCAD скалярный оператор "модуль" с панели инструментов калькулятора в применении к разности векторов даст как раз значение невязки, проверим это утверждение на небольшом тесте:

С учётом всего сказанного, реализуем метод Крамера и проверку полученного решения:

В теле функции det оператор |A| — это не модуль числа с панели "Калькулятор", а похожая внешне кнопка "Определитель" с панели "Матрицы"!

Классический метод Гаусса с приведением матрицы к верхнему треугольному виду подробно изучается в базовом курсе высшей математики. Реализуем самую простую его разновидность, выбирающую ведущий элемент на главной диагонали матрицы, то есть, работающую в предположении, что значение A1,1,≠0 . Так как эта подпрограмма "нулевого" уровня, назовём её Gauss0 , а более сложную Gauss напишем отдельно.

Для удобства вектор правой части 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 🙂 У меня точность решения на использованном тесте выросла при этом вчетверо.

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

Вот расчёт методом Якоби с критерием выхода "максимальный модуль разности между проекциями x1 (k) i и x0 (k) i стал меньше либо равен заданной точности ε :

Метод Якоби
Метод Якоби

При расчёте r применены операторы "Векторизовать" с панели "Матрицы" и "Модуль" с Калькулятора, а при расчёте элементов С и β деление выполнено "в строчку" для экономии места.

P.S. И ещё про прямые "гауссоподобные" методы решения СЛАУ. Если Вам нужно не пошаговое программирование, а достаточно применения стандартных функций, есть способ проще. С помощью стандартной функции augment можно получить расширенную матрицу системы (поставив "рядом" матрицу A и вектор b ), а с помощью rref привести матрицу к ступенчатому виду с единичным базисным минором. Потом останется извлечь решение с помощью метода submatrix (последний столбец матрицы, которую вернул метод rref ).

Норма вектора |A*x-b| , как и другие нормы в статье, берётся кнопкой |x| с Калькулятора, а не похожей на неё кнопкой с панели "Матрицы".

Решение нелинейных уравнений и систем уравнений в пакете MathCAD

Вычисление корней численными методами включает два основных этапа:

· отделение корней;

· уточнение корней до заданной точности.

Рассмотрим эти два этапа подробно.

Отделение корней нелинейного уравнения

Учитывая легкость построения графиков функций в MathCAD , в дальнейшем будет использоваться графический метод отделения корней.

Пример. Дано алгебраическое уравнение

Определить интервалы локализации корней этого уравнения.

Пример. Дано алгебраическое уравнение

Определить интервалы локализации корней этого уравнения.

На рисунке приведен график функции , построенный в MathCAD . Видно, что в качестве интервала изоляции можно принять интервал . Однако уравнение имеет три корня. Следовательно, можно сделать вывод о наличии еще двух комплексных корней. ¨

Уточнение корней нелинейного уравнения

Для уточнения корня используются специальные вычислительные методы такие, как метод деления отрезка пополам, метод хорд, метод касательных (метод Ньютона) и многие другие.

Функция root . В MathCAD для уточнения корней любого нелинейного уравнения (не обязательно только алгебраического) введена функция root , которая может иметь два или четыре аргумента, т.е. или , где – имя функции или арифметическое выражение, соответствующее решаемому нелинейному уравнению, – скалярная переменная, относительно которой решается уравнение, – границы интервала локализации корня.

Пример. Используя функцию , найти все три корня уравнения , включая и два комплексных.

Zoloto585CPA

Заметим, что для вычисления всех трех корней использовался прием понижения порядка алгебраического уравнения, рассмотренный в п. 8.1.1. ¨

Функция root с двумя аргументами требует задания (до обращения к функции) переменной начального значения корня из интервала локализации.

Пример 8.1.5. Используя функцию root , вычислить изменения корня нелинейного уравнения при изменении коэффициента а от 1 до 10 с шагом 1.

Функция polyroots . Для вычисления всех корней алгебраического уравнения порядка (не выше 5) рекомендуется использовать функцию polyroots . Обращение к этой функции имеет вид polyroots (v) , где v – вектор, состоящий из n +1 проекций, равных коэффициентам алгебраического уравнения, т.е. . Эта функция не требует проведения процедуры локализации корней.

Пример. Используя функцию polyroots , найти все три корня уравнения , включая и два комплексных

Блок Given . При уточнении корня нелинейного уравнения можно использовать специальный вычислительный блок Given , имеющий следующую структуру:

Решаемое уравнение задается в виде равенства, в котором используется «жирный» знак равно, вводимый с палитры Логичес­кий .

Ограничения содержат равенства или неравенства, которым должен удовлетворять искомый корень.

Функция Find уточняет корень уравнения, вызов этой функции имеет вид Find ( x ), где x – переменная, по которой уточняется корень. Если корня уравнения на заданном интервале не существует, то следует вызвать функцию Minerr ( x ), которая возвращает приближенное значение корня.

Для выбора алгоритма уточнения корня необходимо щелкнуть правой кнопкой мыши на имени функции Find ( x ) и в появившемся контекстном меню (см. рисунок) выбрать подходящий алгоритм.

Аналогично можно задать алгоритм решения и для функции Minerr ( x ).

Использование численных методов в функциях Find ( x ), Minerr ( x ) требует перед блоком Given задать начальные значения переменным, по которым осуществляется поиск корней уравнения.

Пример. Используя блок Given , вычислите корень уравнения в интервале отделения .

Решение систем уравнений

В зависимости от того, какие функции входят в систему уравнений, можно выделить два класса систем:

· алгебраические системы уравнений;

· трансцендентные системы уравнений.

Среди алгебраических систем уравнений особое место занимают системы линейных алгебраических уравнений (СЛАУ).

Системы линейных алгебраических уравнений

Системой линейных алгебраических уравнений (СЛАУ) называется система вида:

В матричном виде систему можно записать как

где – матрица размерности , – вектор с проекциями.

Для вычисления решения СЛАУ следует использовать функцию lsolve , обращение к которой имеет вид: lsolve (А, b ), где А – матрица системы, – вектор правой части.

Решение систем нелинейных уравнений

MathCAD дает возможность находить решение системы уравнений численными методами, при этом максимальное число уравнений в MathCAD 2001 i доведено до 200.

Для решения системы уравнений необходимо выполнить следующие этапы.

Задание начального приближения для всех неизвестных, входящих в систему уравнений. При небольшом числе неизвестных этот этап можно выполнить графически, как показано в примере.

Пример. Дана система уравнений:

Определить начальные приближения для решений этой системы.

Видно, что система имеет два решения: для первого решения в качестве начального приближения может быть принята точка (-2, 2), а для второго решения – точка (5, 20). ¨

Вычисление решения системы уравнений с заданной точностью . Для этого используется уже известный вычислительный блок Given .

Функция Find вычисляет решение системы уравнений с заданной точностью, и вызов этой функции имеет вид Find ( x ), где x – список переменных, по которым ищется решение. Начальные значения этим переменным задаются в блоке < Начальные условия >. Число аргументов функции должно быть равно числу неизвестных.

Следующие выражения недопустимы внутри блока решения:

· ограничения со знаком ¹ ;

· дискретная переменная или выражения, содержащие дискретную переменную в любой форме;

· блоки решения уравнений не могут быть вложены друг в друга, каждый блок может иметь только одно ключевое слово Given и имя функции Find (или Minerr ).

Пример. Используя блок Given , вычислить все решения системы предыдущего примера. Выполнить проверку найденных решений.

Решение системы линейных уравнений в MathCAD

Решение системы линейных уравнений методом Зейделя в Mathcad
Помогите, пожалуйста! Необходимо решить систему методом Зейделя с точностью 0,001 в mathcad.

Решение систем линейных уравнений в Mathcad
Решаю систему линейных уравнений в Mathcad. Но не могу почему-то применить функцию find. В чем.

Решение систем линейных уравнений в Mathcad
Почему не выполняется проверка?В результате проверки должен получиться вектор b, а получается.

Программирование в Mathcad. Решение систем линейных алгебраических уравнений
Решить систему уравнений с тремя неизвестными методом Гаусса-Жордана, методом простой итерации.

Zoloto585CPA

Добавить комментарий

Ваш адрес email не будет опубликован. Обязательные поля помечены *