Раздел: Документация
0 ... 91 92 93 94 95 96 97 ... 177 Пример 8.25. Влияние обусловленности системы на точность результата Классической плохо обусловленной системой линейных уравнений является система, матрица коэффициентов которой является матрицей Ги.гьберта. Матрица Гильберта А — это матрица, значение элементов которой определяется <рормулой Aj—1/(1iviei wj — индексы элемента (отсчет строк и столбцов матрицы немегсм с нуля). Создадим функцию, которая будет форм»фо-вать матрицы Гилберта произвольной размерности. A(N):= for i 6 0.. N - I for j с 0.. N - 1 1 M. A(3) j + i + 1 M
Очевидно, что чем больше будет матрица Гильберта, тем в меньшей степени будут отличаться нижние рядки. Соответственно, тем хуже будет матрица обусловлена. Посмотрим, как зависят опреде.гитель матрицы Гильберта и ее число обусловленности от ее размерности. л-4 А(2) =0.083А(3) =4.63x10 conde(A(2)) = 19.333 conde(A{3» = 526.159 а<4) = 1,653x10 conde(A(4))= 1.561x10 Наше предположение подтвердилось. Стоит ожидать, что с увеличением размерности матрицы Гильберта погрешность определешгя корней будет возрастать. Проверим это. В первую очередь создадим функцию, которая будет формировать вектор правых частей, соразмерный матрице Гильберта. Так как нас интересует только погрешность, не имеет значения, какие у системы будут корни. Поэтому мы будем заполнять вектор правых частей одними единицами. В(М):= for icO.N - 1 V.<- I Чтобы можно было оценить погрешность численного метода, необходимо решить систему точно. В нашем случае это возможно, так как значения элементов матрицы Гильберта могут быть представлены простыми дробями, а в векторе правых частей имеются только целые числа. Результат, возвращенный функцией l.iolve, присвоим переменной Тг. N:=3 Тс:= lsolve(A(N),B(N)) 3 % -24 Далее находим среднюю ошибку численного метода. Для этого от вектора с определенными им корнями отнимаем вектор Тг. При этом будет получен вектор абсолютных ошибок. Чтобы найти ошибки относительные, данный вектор делим па вектор Т.. применив оператор векторизации. Затем умножаем вектор относительных ошибок на 100. чтобы получить вектор прииеитных ошибок. И, наконец, находим среднее арифметическое ошибок, задействовав функцию mean. Все описанные действия легко совместить в одной формуле. егг% := mean! (Isolve* А( N). В( N)) - Т, )00 етг% = 2,536x10 13 Изменяя значения N, находим среднюю ошибку н число обусловлейпасти для матриц Гильберта разной размерности. Полученные данные заносим в табл. 8,1. Таблица 8.1. Зависимость средней ошибки и числа обусловленности от размерности матрицы Гильберта
Как видно из табл. 8,1, с увеличением числа обусловленности погрешность работы численного метода резко возрастает. Ошибка определения корня может составлять сотни процентов, делая результат абсолютно бесполезным. В нашем случае это наблюдается, когда размерность матрицы Гильберта достшает 13. Пример 8.25 показывает, сколь важно оценивать степень обусловленности матрицы коэффициентов. Если система обусловлена плохо, то численному решению нужно Предпочесть аналитическое. Если же это невозможно, то к результату следует относиться крайне осторожно. В отдельных случаях точность численного решения можно повысить, используя специальные техники вроде сингулярного разложения (о нем МЫ поговорим ниже). Традиционно основным методом решения СЛАУ считается метод Гаусса. Его изучают на занятиях по математике во всех вузах, его описывают во всех книгах по численным методам. Однако на практике метод Гаусса не применяется, так как су шествуют более эффективные методы. В основе функции Isolve лежит алгоритм, основанный иа LU-разложенин матрицы. Мы, следуя установленному а этой главе принципу обучения, реализуем данный метод сампстоятельио. LU разложением квадратной матрицы А называется ее представление в виде произведения нижней треугольной матрицы L (то есть матрицы, все ненулевые элементы которой лежат только на главной диагонали и ниже ее) и верхней треугольной матрицы U (матрицы, элементы которой располагаются на главной диагонали и выше нее). К примеру, LU-разложение матрицы размерности 3x3 можно записать следующим образом: 1,1 4.1 О 1 U„ „ ил , ил А (А 0,0 ii 0.1 ~0,2 и 1,1 о и 1.2 10.0 "0.1 1,0 1,1 А0.2 1.2 12.0 "2,1 "2.2. Зная разложение матрицы коэффициентов А на треугольные матрицы L и U, решить систему очень просто. Действительно, если справедливо А Х-В (здесь X — вектор неизвестных, В - вектор правых частей), то исполняться будет и L-UX-B, Произведение UX даст вектор. Обозначим его Y. Найти вектор V можно, решив систему LY~B. Так как матрица L является нижней треугольной, то jto несложно сделать методом прямой подстановки. Данный метод заключается в следующих действиях. □Переменная Ya равна В/Ц, так как в верхней (.троке матрицы L все остальные зде-мекты — вто нули. □Перемножение второй строки матрицы коэффициентов L и вектора неизвестных Y даст уравнение LYL,yY,-!*,. Так как значение Y0 было получено иа прошлом щаге. то из данного уравнения легко найти У",: Yl-(BI-L1fl-Y„)/L11. □В общем случае умножение n-й строки матрицы L на вектор Y даст следующее урав- ,+L -Y-В Так так неизвестные от Y. ДО Y , были найдены ранее, то в данном уравнении будет только одно неизвестное — Y.. Выразив его, получим рабочую формулу метода прямой подстановки; п-1 L .-Y. П,1 I i = 0 П,[1 Найдя вектор Y, решаем систему UX-Y, определяя корни исходной системы АХ-В. Так как матрица U — верхняя треугольная, то это легко сделать с помощью метода обратной подстановки. Идея данного метода точно такая же, как у метода прямой подстановки. Отличие состоит в том, что неиавестн ые определяются в противоположном порядке. Рабочая формула метода обратной подстановки следующая: last(Y) U n.n Итак, польза от представления матрицы коэффициентов а виде произведения треугольных матриц заключается в том, что при этом становится возможным нахождение корней посредством подстановок. Написать же программы прямой и обратной подстановки, зная рабочие формулы данных методов, проще простого: direct(L,B) for о е 1„1ая(В) В - V L ,Y. П £t D.l 1 Y- nL и, и invert(U.Y) ; N <- lastCY) XN *~ YN + UN.N for rcN-1.0 N Y - У U -Х- U Как же можно на практике произвести LU-разложение? В Mathcad для этого используется весьма эффективный алгоритм Краута (Crout). описать который можно следующими шагами. 0 ... 91 92 93 94 95 96 97 ... 177
|
