Гурский Д., Турбина Е. - Вычисления в MathCad 12 (1077322), страница 70
Текст из файла (страница 70)
Проще это сделать, не переписывая код функции lu_brain, а создав новую функцию, использующую lu_brainкак подпрограмму. Назовем мы эту функцию LU. Ее алгоритм будет заключаться в следующих действиях.•Пытаемся произвести LU-разложение, вызвав функцию lu_brain. Если при этомвозникает ошибка («отловить» ее можно, используя оператор on error), то генерируем матрицу перестановок и умножаем на нее разлагаемую матрицу. Полученнуюматрицу передаем функции lu_brain. Описанные действия повторяются до тех пор,пока разложение не удастся осуществить.302 •Глава 8.
Решение уравнений и систем уравнений•Вполне может оказаться так, что матрица вырождена и, соответственно, ее нельзяпредставить в виде произведения треугольных матриц. Это означает, что ошибкаделения на нуль будет возникать, какая бы матрица перестановок ни была использована. Поэтому важно ограничить количество итераций, в течение которых разложение должно быть произведено, например величиной в 2-N, где N — размерностьматрицы.
Если лимит на количество итераций окажется превышенным, нужно возвратить сообщение об ошибке, информирующее пользователя о сингулярностиматрицы.• Если у матрицы коэффициентов А при произведении LU-разложения были переставлены строки, то при решении системы уравнений такая же перестановка должна быть осуществлена и над вектором правых частей В. Чтобы это можно былосделать, необходимо знать, какая матрица перестановок была использована алгоритмом LU-разложения.
Поэтому соответствующая матрица должна возвращатьсяфункцией LU как третий элемент вектора ответа. Если же матрица перестановок неприменялась, то возвращена должна быть единичная матрица. Кстати, именно такработает встроенная функция Mathcad lu, отвечающая за LU-разложение. Она возвращает матрицу перестановок, L и U матрицы, слитые в одну матрицу.Код функции LU следующий:(п <- О N <- rows(A))LU(A) :=while Ion error lu_brain(A)error("Error! Probably, matrix is singular!') if n > 2-Nperm_matr<— per(N)A <— perm_matr-A(n <— n + 1 continue )(answer <— lu_brain(A)break)answer. - <— perm_matr if n > 0identity(N) otherwiseanswer 0,2answerПроверим эффективность функции LU, разложив с ее помощью матрицу, с которой несправилась функция lu_brain:res =_v-21 О4О 112 6М:=48-11261 ООО 717О 0 1OO-25V.0 1res := LU(M)ОЛreso,oreso,i12 648-1V-235 ,8.2.
Решение систем уравнений * 3 0 3Проверка показывает, что функция LU справляется со своей задачей хорошо. Значит,нам осталось объединить возможности функций LU, direct и invert в программе, непосредственно производящей решение СЛАУ. Ее код будет очень прост. Единственнаятонкость состоит в том, что нужно не забыть умножить вектор правых частей на матрицу перестановок, использованную при вычислении LU-разложения. Также стоитпроверить матрицу коэффициентов на сингулярность, вычислив ее ранг (для этогоможно использовать встроенную функцию rank).lu_solver(A,B):=errorC'Matrix is singular") if rank(A) Ф rows (A)(M<- LU(A) L <- 1VLU <- M.peВ <— pernvB(Zl <- direct(L.B) Z2 <-invert(U,Zl))Z2Решаем посредством функции lusolver систему линейных уравнений и сравниваемполученный ответ с точным аналитическим решением:(1А:= 42 2>1 -1_2 4 5 ;(ЛВ:= 1 K:=lsolve(A,B)->2912-16lu_solver(A,B)-K= 6.661x10-16.-16V.-6.106x10V9Итак, написанный нами алгоритм работает очень неплохо.
Точность решения системыиз трех уравнений находится на уровне предельной точности численных расчетов.Причем ответ рассчитывается моментально.Почему мы разбираем особенности решения систем линейных уравнений столь подробно? Дело в том, что это, пожалуй, наиболее важный вопрос в теории численных методов.Алгоритмы решения СЛАУ используются доброй половиной всех практически важных численных методов. К примеру, на них основывается нахождение обратной матрицы, они применяются при решении систем нелинейных уравнений, краевых задач,дифференциальных уравнений в частных производных; без них невозможны многиеалгоритмы интерполяции и регрессии, методы отыскания собственных значений матрици многое-многое другое. По причине основополагающей роли методов решения СЛАУ,в курсах численных методов именно с них принято начинать изложение материала.Продемонстрируем, как в Mathcad, на основании алгоритма решения СЛАУ, вычисляется обратная матрица.Пусть дана квадратная матрица М и нужно найти матрицу, обратную данной.
По определению обратной матрицы, М-М~'=1, где I — единичная матрица соответствующей размерности. Исходя из того, как производится матричное умножение, можно записатьМ-(М~').=1., где (М~'). — j-й столбец обратной матрицы, I. — j-й столбец единичнойматрицы. По сути, М-(М~')=1. является системой линейных уравнений, где М — матрица коэффициентов, (М~') — вектор неизвестных, I — вектор правых частей. Решивданную систему, можно найти j-й столбец обратной матрицы. Таким образом, чтобынайти обратную матрицу для матрицы размерности N, нужно N раз решить систему изN линейных уравнений.304 •Глава 8.
Решение уравнений и систем уравненийНа языке программирования Mathcad реализовать функцию, вычисляющую обратнуюматрицу, можно следующим кодом:2inverter(M):= (N <- rows(M) В <г- identity(N))21error("Matrix must be squarte!") if rows(M) *• cols(M) A:= 4 1 -1"2 4error("Matrix is singular!") if rank(M) *• Nfor i e 0.. N - 11 -0.222 -0.444Alu solverinverter(A) = -2112 -0.889 -0.778,/LU-разложение квадратной матрицы, которому мы посвятили столько внимания, полезно не только при решении СЛАУ. В Mathcad оно используется для нахожденияопределителя матрицы. В линейной алгебре доказано, что определитель равен произведению диагональных элементов верхней треугольной матрицы U.
Попробуем, воспользовавшись данным фактом, самостоятельно создать функцию, вычисляющуюопределитель матрицы.det(M) :=error("Matrix must be squarte!" ) if rows(M) *• cols(M)U <- LU(ML ,rows(M)-lПi=0U.A:=124-14-25det(A) = 9=9Функция Lsolve использует для решения СЛАУ LU-разложение матрицы коэффициентов.
Однако существуют и другие матричные разложения, полезные при решениисистем линейных уравнений: сингулярное разложение, QR-разложение, разложениеХолецкого. Каждое из них предназначено для эффективного решения отдельного типаСЛАУ.Наиболее востребованным на практике матричным разложением является так называемое сингулярное разложение. Дело в том, что такой подход дает возможность решать системы уравнений, матрицы коэффициентов которых практически вырождены.Во многих случаях сингулярное разложение позволяет найти решение там, где методLU-разложения или метод Гаусса оказываются бессильными.Сингулярное разложение заключается в представлении матрицы коэффициентовв виде произведения трех квадратных матриц той же размерности: A=U-diag(S)-VT.Здесь diag(S) — диагональная матрица, непулевыми элементами которой являются такназываемые сингулярные числа разлагаемой матрицы.
В Mathcad вектор сингулярныхчисел возвращает встроенная функция svds. Матрицы U и V можно найти, используяфункцию svd (они возвращаются слитыми в одну матрицу).Чтобы найти на основании сингулярного разложения решение СЛАУ, нужно последовательно решить три системы линейных уравнений, в качестве матриц коэффициентов которых выступают полученные в результате проведения разложения матрицыU, V и diag(S).
Так как обычно данные матрицы обусловлены куда лучше исходнойматрицы, то для решения соответствующих систем можно применять стандартные методы вроде метода LU-разложения. В случае матрицы diag(S) решение можно полу-8.2. Решение систем уравнений * 3 0 5чить и проще, поделив элементы вектора правых частей на соответствующие им элементы главной диагонали.sing_solver(A,B) := (m <- svd(A) п <- rows(m) k <- cols(m))U <r- submatrix m,0,Il,0,k- 12V<- submatrix m, —,n - l , 0 , k - 1T(zi<-lu_solver(U,B) Z2<-(Z1 -=- svds(A)) Z3<-lu_solver(y ,Z2))Z3Проверьте эффективность программы sing_solver самостоятельно, решив плохо обусловленную систему (например, матрица коэффициентов которой является матрицейГильберта).Интересной особенностью сингулярного разложения является то, что оно может бытьиспользовано по отношению к неквадратной матрице.
Это дает возможность решатьсистемы уравнений, в которых количество уравнений и количество неизвестных несовпадает. Также на основании сингулярного разложения может быть реализован алгоритм расчета линейной регрессии.В Mathcad имеются функции и других матричных разложений, которые могут бытьиспользованы для повышения эффективности решения специфических видов СЛАУ.Так, если матрица коэффициентов системы линейных уравнений является симметричной (то есть А-А т ) и положительно определенной (все элементы больше или равнынулю), то вместо LU-разложения гораздо техничнее применять так называемое разложение Холецкого.
Данное разложение заключается в представлении матрицы коэффициентов А в виде произведения нижней треугольной матрицы L и ее транспонированной формы: A-LL T . Чтобы найти решение СЛАУ на основании разложения Холецкого,нужно решить две системы линейных уравнений, матрицами коэффициентов в котоTрых являются L и L . Так как данные матрицы являются соответственно нижней треугольной и верхней треугольной, то для этого можно использовать методы прямойи обратной подстановки.Как это ни удивительно, но симметричные матрицы встречаются в приложениях довольно часто (к примеру, такой матрицей является матрица Гильберта).
Поэтому разложение Холецкого имеет заметное практическое значение.В Mathcad разложение Холецкого можно провести, обратившись к встроенной "функции chotecky. Результатом работы данной функции является матрица L.chol(A,B):= m <- cholesky(A)zl-direct(m, B)z2«Tt(n ,zl)invert\z211n2 31 1 1A:= - - 2 3 4113 415хB:= 2f 27chol(A,B) = -192v210y306 •Глава 8. Решение уравнений и систем уравненийВ отдельных случаях может быть полезно так называемое QR-разложение. Данное разложение заключается в представлении матрицы в виде произведения ортогональнойматрицы Q (то есть QT=Q~') и верхней треугольной матрицы R. Впрочем, на практикеQR-разложение для решения СЛАУ применяется нечасто, так как оно требует почтив два раза больше операций по сравнению с LU-разложением. Большее значениеQR-разложение имеет при решении специфических задач (например, определениисобственных значений и векторов матрицы).Все способы решения СЛАУ, которые мы рассмотрели выше, являются прямыми.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
