Раздел: Документация
0 ... 93 94 95 96 97 98 99 ... 177 На языке программирования Mathcad описанный алгоритм можно реализовать следующим образом: рсп>0:я&г i б 0.. N - 1 cle*r Vj «- i white 1 v «— clear v for j«0.. N - I md index«- round(md(last(v))) a *- break if lasi(v) = 0 (v «- submairixfv. 1, lasi(v).O.O) continue) if mdjndcx=0 (v *— Bubmatrixfv.OJasv) - 1,0,0) continue ) if md index = last(v) jv«-ete£fubiranixfyAm break if M * identity(rO M В контексте написанных ранее программ функция per не содержит никаких новых приемов. Единственное, стоит пояснить, как производится случайный выбор элемента на вектора v. Для этого с помощью функции rnd генерируется случайное число от 0 до п, где п — индекс последнего в векторе элемента. Так как сгенерированное число может быть дробью, оно округляется до ближайшего целого, для чего используется функция round. Проверим, как работает функция per, сгенерировав три матрицы перестановки размерности 3: репЧ) Итак, функция per работает так, как было задумано. Теперь мы должны так модифицировать код, производящий LU-разложен и с, чтобы он при возникновении ошибки деления на нуль осуществлял перестановку строк матрицы. Проще это сделать, не переписывая код функции lu bratn, а создав новую функцию, использующую lu brain как подпрограмму. Назовем мы эту функцию LTJ. Ее алгоритм будет заключаться в следующих действиях. □ Пытаемся произвести LU-раэложенне, вызвав функцию lu brain. Если при этом возникает ошибка («отловить» ее можно, используя оператор on error), то генерируем матрицу перестановок и умножаем на нее разлагаемую матрицу. Полученную матрицу передаем функцни lu brain. Описанные действия повторяются до тех лор. пока разложение не удастся осуществить.
□В пол не может оказаться так, что матрица вырождена н, соответственно, ее нельзя прелставить в виде произведения треугольных матриц. Это означает, что ошибка деления на нуль будет возникать, какая бы матрица перестановок ни была использована. Поэтому важно ограничить количество итераций, в течение которых разложение должно быть произведено, например величиной в 2-N, где N — размерность матрицы. Если лимитна количество итераций окажется превышенным, нужно возвратить сообщение об ошибке, информирующее пользователя о сингулярности матрицы. □Если у матрицы коэффициентов А при произведении LU-разложения были переставлены строки, то при решении системы уравнений такая же перестановка должна быть осуществлена н над вектором правых частей В. Чтобы это можно была сделать, необходимо знать, какая матрица перестановок была использована алгоритмом LU-разложения. Поэтому соответствующая матрица должна возвращаться функцией LU как третий элемент вектора ответа. Если же матрица перестановок не применялась, то возвращена должна быть единичная матрица. Кстати, именно так работает встроенная функция Mathcad ш, отвечающая за LU-разложение. Она возвращает матрицу перестановок, L и U матрицы, слитые в одну матрицу, Код функции LU следующий? ЩА):= (Df-0 N <- rows(A)) while I on error lu brain(A) error("Error! Probably, matrix is singular!) if n > 2- N perm malr«— per(N) A 4— perm matr-A (»«-п+ 1 continue ) (answer <-lu bram(A) break) answer *- perm matr if п > 0 0.2 answer answern ,, <— identityfN) otherwise Проверим эффективность функции LU, разложив с ее помощью матрицу, с которой не справилась функция liii4.nn ( 1 2 бЛ М:= 4 8 -2 3 -1 5 У res :=ЩМ) res =
(reS0,2l,IW0,0rc,0,l = 1 1 6У 4 8-1 <-1 3 5 j Проверка показывает, что функция LU справляется со своей задачей хорошо. Значит, нам осталось объединить возможности функций LU. direct и invert в программе, непосредственно производящей решение СЛАУ. Ее код будет очень прост. Единственная тонкость состоит о том, что нужно не забыть умножить вектор правых частей на матрицу перестановок, использованную при вычислении LU-разложення. Также стоит проверить матрицу коэффициентов на сингулярность, вычислив ее ранг (для этого можно использовать встроенную функцию rank). lu solver(A,B) := error (""Matrix is singular") if rank(A) * rows(A) (м <- ЩА) Lt-Mjj й U «- MQ , рет ш <- Mfl J В 4- (Zl Z2 pcrnvB e— direct(L,B) 72 < invert(U,Zl)) Решаем посредством функции lu solver систему линейных уравнений и сравниваем полученный ответе точным аналитическим решением:
К.:=1зо1уе(\.В) 9 1 2 1и во1ует(А,В)-К = -2.22x10 - 16 6.661x10 16 1,-6.106*10 г 16 Итак, написанный нами алгоритм работает очень неплохо. Точность решения системы ИЗ трех уравнений находится на уровне предельной точности численных расчетов. Причем ответ рассчитывается моментально. Почему мы разбираем особенности решения систем линейных уравнений столь подробно? Дело в том, что это, пожалуй, наиболее важный вопрос в теории численных методов. Алгоритмы решения СЛАУ используются доброй половиной всех практически важных численных методов. К примеру, на них основывается нахождение обратной матрицы, они применяются при решении систем нелинейных уравнений, краевых задач, дифференциальных уравнений в частных производных; без них невозможны многие алгоритмы интерполяции и регрессии, методы отыскания собственных значений матриц и многое-многое другое. По причине основополагающей роли методов решения СЛАУ, в курсах численных методов именно с них принято начинать изложение материала Продемонстрируем, как в Mathcad, на основании алгоритма решения СЛАУ, вычисляется обратная матрица. Пусть дана квадратная матрица М и нужно найти матрицу, обрапгую данной. По определению обратной матрицы. ММ 1-1,где I - единичная матрица соответствующей размерности. Исходя га того, как производится матричное умножение, можно записать М(М"-1р, где (М")— j-й столбец обратной матрицы, I, — j-Й столбец единичной матрицы. По сути, М-(М-1)-является системой линейных уравнений, где М — матрица коэффициентов, (Мч). — вектор неизвестных, Ij — вектор правых частей. Решив данную систему, можно найти j-й столбец обратной матрицы. Таким образом, чтобы найти обратную матрицу для матрицы размерности N, нужно N раз решить систему из N линейных уравнений. 0 ... 93 94 95 96 97 98 99 ... 177
|
