Гурский Д., Турбина Е. - Вычисления в MathCad 12 (1077322), страница 68
Текст из файла (страница 68)
Впрочем, можно обойтись без всех этих операций, если использоватьдля решения системы оператор solve (аналитический расчет) или же блок Given-Find(численный расчет). Подробно о данном блоке мы поговорим ниже, а пока ограничимся лишь примером.Пример 8.24. Решение СЛАУ в случае неравенства количество уравненийи количества неизвестныхЕсли количество неизвестных превышает количество уравнений, то одни неизвестные можновыразить как функциональные зависимости от других неизвестных:294 *Глава 8. Решение уравнений и систем уравненийх+ у + z+ k = 4x-y+z-k=Osolve, х, у , z —>(1 - к + 2 1)- х + у + z+ k = 2Если количество уравнений превышает количество неизвестных, то найти корни аналитическиможно, использовав оператор solve.х+ у + z = 3х - у + z= 1solve,х,у, z-» (1 1 1)-х + у + z = 1V-x-y-z = -S,Численно решить систему, в которой количество уравнений превышает количество неизвестных, позволяет блок Given-Find.х:=0у:=0z:=0Givenx-y + z=lx+y + z =3-x+y +z = lfind(x,y,z) =-x-y-z = -31Говоря о численном решении СЛАУ, важно затронуть вопрос о точности используемых в Mathcad алгоритмов.
В случае обычных систем, вроде приводимых в задачниках, точность функции LsoLve или блока Given-Find очень высока. Ответ обычно веренвплоть до 14-15-го знака мантиссы, что близко к предельной теоретической точности численных методов на машинах с 64-битовыми числами с плавающей точкой. Однако с увеличением размеров системы точность обычно падает.
Также крайне важно, насколькоразнородные уравнения образуют систему. Если коэффициенты двух или несколькихуравнений близки, то при малых возмущениях в матрице коэффициентов А или векторе правых частей В произойдут большие изменения в векторе неизвестных Х=А~'-В.Это означает, что численный метод будет неустойчив и, следовательно, резко повысится вероятность получения ответа, содержащего большую погрешность.Системы, содержащие практически линейно зависимые уравнения, называются плохообусловленными. Количественной мерой обусловленности может являться определитель: чем он ближе к нулю, тем хуже обусловлена система.
Однако есть и более объективная характеристика, называемая числом обусловленности (condition number). Число обусловленности определяется как произведение нормы матрицы коэффициентовА на норму обратной ей матрицы. Чем больше число обусловленности, тем выше вероятность получения ответа со значительной погрешностью. В Mathcad есть функции,позволяющие определять число обусловленности. Данные функции различаются тем,на основании каких норм находится число обусловленности.
Перечислим их: condl(норма LI), cond2 (норма L2), conde (эвклидова норма), condi (бесконечная норма).Расчет данные функции могут вести как численно, так и аналитически.8.2. Решение систем уравнений * 2 9 5Пример 8.25. Влияние обусловленности системы на точность результатаКлассической плохо обусловленной системой линейных уравнений является система, матрицакоэффициентов которой является матрицей Гильберта. Матрица Гильберта А — это матрица,значение элементов которой определяется формулой А{ =l/(l+i+j), где i и j — индексы элемента(отсчет строк и столбцов матрицы ведется с нуля). Создадим функцию, которая будет формировать матрицы Гилберта произвольной размерности.A(N):=1for i e 0.. N - 1for j e 0..
N - 1lM. . < - — —•>Jj +i+1M112311123 4111345Очевидно, что чем больше будет матрица Гильберта, тем в меньшей степени будут отличатьсянижние рядки. Соответственно, тем хуже будет матрица обусловлена. Посмотрим, как зависятопределитель матрицы Гильберта и ее число обусловленности от ее размерности.|А(2)| =0.083—4|А(3)| = 4 . 6 3 X 1 0conde(A(2)) = 19.333conde(A(3)) = 526.159-7|А(4)| = 1.653x10conde(A(4)) = 1.561x10Наше предположение подтвердилось. Стоит ожидать, что с увеличением размерности матрицыГильберта погрешность определения корней будет возрастать. Проверим это.
В первую очередьсоздадим функцию, которая будет формировать вектор правых частей, соразмерный матрицеГильберта. Так как нас интересует только погрешность, не имеет значения, какие у системы будут корни. Поэтому мы будем заполнять вектор правых частей одними единицами.B(N):=for i e 0.. N - 1V. <r- 1Чтобы можно было оценить погреигаость численного метода, необходимо решить систему точно.
В нашем случае это возможно, так как значения элементов матрицы Гильберта могут бытьпредставлены простыми дробями, а в векторе правых частей имеются только целые числа. Результат, возвращенный функцией lsolve, присвоим переменной Т..3N:=3Т. :=lsolve(A(N),B(N))-»\-24Далее находим среднюю ошибку численного метода. Для этого от вектора с определенными имкорнями отнимаем вектор Т.. При этом будет получен вектор абсолютных ошибок. Чтобы найтиошибки относительные, данный вектор делим на вектор Тс, применив оператор векторизации.Затем умножаем вектор относительных ошибок на 100, чтобы получить вектор процентных ошибок.
И, наконец, находим среднее арифметическое ошибок, задействовав функцию mean. Все описанные действия легко совместить в одной формуле.2 9 6 •:• Глава 8. Решение уравнений и систем уравнений.- 13err% = 2.536x10егг% := meanf |"(lsolve(A(N),B(N)) - Т с )+ T J J - 1 0 0Изменяя значения N, находим среднюю ошибку и число обусловленности для матриц Гильбертаразной размерности. Полученные данные заносим в табл.
8.1.Таблица 8 . 1 . Зависимость средней ошибки и числа обусловленности от размерностиматрицы Гильберта5Размерность, N3Числообусловленности526.1Средняя ошибка2.5-10"'75981311114.8-104.8-105.0-105.3107.210"3.3-1070.00030.1733145.7-1О17732.877Как видно из табл. 8.1, с увеличением числа обусловленности погрешность работы численного метода резко возрастает. Ошибка определения корня может составлять сотнипроцентов, делая результат абсолютно бесполезным. В нашем случае это наблюдается,когда размерность матрицы Гильберта достигает 13.Пример 8.25 показывает, сколь важно оценивать степень обусловленности матрицыкоэффициентов.
Если система обусловлена плохо, то численному решению нужнопредпочесть аналитическое. Если же это невозможно, то к результату следует относиться крайне осторожно. В отдельных случаях точность численного решения можноповысить, используя специальные техники вроде сингулярного разложения (о нем мыпоговорим ниже).Традиционно основным методом решения СЛАУ считается метод Гаусса. Его изучаютна занятиях по математике во всех вузах, его описывают во всех книгах по численнымметодам.
Однако на практике метод Гаусса не применяется, так как существуют болееэффективные методы. В основе функции Isolve лежит алгоритм, основанный на LUразложении матрицы. Мы, следуя установленному в этой главе принципу обучения,реализуем данный метод самостоятельно.LU разложением квадратной матрицы А называется ее представление в виде произведения нижней треугольной матрицы L (то есть матрицы, все ненулевые элементы которой лежат только на главной диагонали и ниже ее) и верхней треугольной матрицы U (матрицы, элементы которой располагаются на главной диагонали и выше нее).К примеру, LU-разложение матрицы размерности 3x3 можно записать следующим образом:0),оL2,0L2,0J.2U0,0U0,lU0,2о и 1,1UоU2,2Оv0,01,2l,2iVДД2,0 A 2 , l А 2 , 2Зная разложение матрицы коэффициентов А на треугольные матрицы L и U, решить систему очень просто.
Действительно, если справедливо А-Х=В (здесь X — вектор неизвестных, В — вектор правых частей), то исполняться будет и L-U-X-B. Произведение U-Xдаст вектор. Обозначим его Y. Найти вектор Y можно, решив систему L-Y=B. Так какматрица L является нижней треугольной, то это несложно сделать методо'м прямойподстановки. Данный метод заключается в следующих действиях.8.2. Решение систем уравнений* 297Переменная Yo равна Bo/Loo, так как в верхней строке матрицы L все остальные элементы — это нули.Перемножение второй строки матрицы коэффициентов L и вектора неизвестных Yдаст уравнение L10-Y0+L, J-Y^BJ.
Так как значение Yo было получено на прошломшаге, то из данного уравнения легко найти Y(: Y 1 =(B 1 -L ]0 -Y 0 )/L ] ,.•В общем случае умножение n-й строки матрицы L на вектор Y даст следующее уравнение^, „-Y.+L, ,-Y,+...+L .-Y ,+L -Y =B . Так как неизвестные от Y. до Y .были1,001.11п,п-1п-1пп.ппип-1найдены ранее, то в данном уравнении будет только одно неизвестное — Yn. Выразив его, получим рабочую формулу метода прямой подстановки:п-1L .Y.Y =пi=0п,пНайдя вектор Y, решаем систему UX=Y, определяя корни исходной системы А'Х=В.Так как матрица U — верхняя треугольная, то это легко сделать с помощью метода обратной подстановки.
Идея данного метода точно такая же, как у метода прямой подстановки. Отличие состоит в том, что неизвестные определяются в противоположномпорядке. Рабочая формула метода обратной подстановки следующая:last(Y)Y - У U -х.X =о1 — П1Uп,пИтак, польза от представления матрицы коэффициентов в виде произведения треугольных матриц заключается в том, что при этом становится возможным нахождениекорней посредством подстановок. Написать же программы прямой и обратной подстановки, зная рабочие формулы данных методов, проще простого:invert(U.Y) :=direct(L,B) :=for n e 1.. last(B)n-1Вn - УL ,Y.Z_, n,i lN <- last(Y)XN ^ ~ YN+UN,Nfor n € N - 1.. 0Ni=0Yn - Z_(V Un , i.-X.ln,nXn<-i =nUn,nXКак же можно на практике произвести LU-разложение? В Mathcad для этого используется весьма эффективный алгоритм Краута (Crout), описать который можно следующими шагами.2 9 8 •:• Глава 8.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
