Форсайт Дж., Малькольм М., Моулер К. - Машинные методы математических вычислений (1040536), страница 9
Текст из файла (страница 9)
Матрица, показанная выше, имеет ширину ,пенты 3 и называется трехдиагональной матрицей. Некоторые численные методы, применяемые для хранимых матриц, весьма отличаются от методов, приспособленных для разреженных матриц. Методы для хранимых матриц более универсальны, и им будет уделено основное внимание. Эти методы можно модифицировать таким образом, чтобы обрабатывать матрицы, находящиеся во внешней памяти, ленточные матрицы и другие типы больших нли умеренно разреженных матриц. 3.1. Линейные системы с хранимыми матрицами В этом параграфе мы обсудим решение линейной системы алгебраических уравнений АХ=Ь с хранимой пхп-матрицей А и векторами Ь и х порядка п.
Будем предполагать, что А — невырожденная матрица. Если А вырохдена, то с точки зрения теории это обнаружится в процессе вычислений, хотя на практике установить вырожденность матрицы может быть нелегко. 3. ЛИНЕЙНЫЕ СИСТЕМЫ УРАВНЕНИЙ Почти всегда применяемый алгоритм, являющийся одним из старейших численных методов,— это метод последовательного исключения неизвестных, называемый обычно именем Гаусса.. Исследования, проведенные в период с 1955 по 1965 год, выявили важность двух аспектов гауссова исключения, не оцененных в более ранних работах: выбора ведущих элементов и надлежащей интерпретации влияния ошибок округления, Гауссово исключение и другие аспекты матричных вычислений подробно разбираются в книгах Форсайта, Молера (1967) и Стьюарта (1973).
Читатель, желающий иметь больше информации, чем мы сообщаем в этой главе, должен воспользоваться этими ссылками. Чтобы проиллюстрировать алгоритм, рассмотрим следующий пример порядка 3: — 3 26 х, =- 4 Зто, в сущности запись системы из трех уравнений 10х,— 7х, =7, — Зх, + 2х, + бх, = — 4, 5х, — х,+5х, =6. На первом шаге с помощью первого уравнения х, исключается из других уравнений. Зто достигается прибавлением первого уравнения, умноженного на 0.3, ко второму уравнению и при- бавлением первого уравнения, умноженного на — 0.5, к третьему уравнению. Величины 0.3 и — 0.5 называются множителями.
( О -0.1 6 ~(*,)-.(б.1). На втором шаге можно было бы использовать второе уравнение, чтобы исключить х, из третьего уравнения. Однако коэффициент при х., во втором уравнении — малое число — О.1. Вследствие этого последние два уравнения переставляются. В данном при- мере это делать не нужно, поскольку здесь нет ошибок округле- ния, но в общем случае такие перестановки очень важны. 0 2.5 5 х, == 2.5 Теперь с помощью второго уравнения можно исключить х, из третьего уравнения. Зто достигается прибавлением второго Е1. линеЙные системы с хРлиимыми млтРиикми 47 уравнения, умноженного на 0.04, к третьему уравнению. (Каков был бы множитель, если бы уравнения не переставлялисьу) 0 255 х, = 25 Последнее уравнение имеет теперь вид 6.2х,=6.2.
Решая его, получим х,=1. Зто значение подставляется во вто- рое уравнение: 2.5х,+(5) (1) =2.5. Следовательно, х,= — 1. Наконец, значения х, и х, подставляются в первое уравнение: 10х,+( — 7) ( — 1) =7. Следовательно, х„=О. Зто решение легко проверить, используя исходные уравнения: — 3 26 — ! =- 4 В общем случае гауссово исключение состоит из двух этапов; прямого хода и обратной подстановки. Прямой ход состоит из и — 1 шагов.
На я-м шаге кратные й-го уравнения вычитаются из оставшихся уравнений с целью исключить я-е неизвестное. Если коэффициент при хе <мал», то рекомендуется переставить уравнения перед исключением. Обратная подстановка заключается в решении последнего уравнения относительно х„, затем предпоследнего уравнения относительно х„ „ и т. д., пока из первого уравнения не будет вычислено х;. Алгоритм в целом может быть компактно записан в матричных обозначениях. Для рассмотренного выше примера положим М,= 03 10 Тогда МА= 0 — 01 6, М1Ь= 61 3, линейные системы уРАВнениЙ 48 Пусть Р,= 001, М,= 01 0 Тогда 10 — 7 О 7 М,Р,М,А = 0 2.5 6, М,Р,М,Ь = 2.6 0 0 6.2 6.2 Главное теперь состоит в том, что вследствие способа выбора матриц М,, Р, и М, произведение (7=М,Р,М,А является верхней треугольной матрицей, т.
е. все ненулевые элементы находятся в верхней правой половине этой матрицы. Для такой матрицы система уравнений (ух=с легко решается обратной подстановкой. Если положить с= М,Р,М,Ь, то система (ух=с имеет то же решение, что и исходная система Ах=Ь. Аналогичные соотношения справедливы и в общем случае. Пусть Рд, Ь=1,..., н — 1, обозначает матрицу, полученную из единичной матрицы 1 той же перестановкой строк, какая применялась к строкам А на й-м шаге исключения. Пусть МА обозначает матрицу, полученную из единичной матрицы записью в поддиагональные позиции Ь-го столбца множителей, использованных на Ь-м шаге. Каждая из матриц РА является матрицвй перестановки, каждая из матриц М, является простой нижней треугольной матрицей.
Обозначим через М и (А' произведения М = М„,Р„,... М,Р,М,Р„ и=МА. Тогда 17 — верхняя треугольная матрица, система (ух=МЬ легко решается относительно х и х является также решением исходной системы Ах=Ь. Матрица М вЂ” не обязательно нижняя треугольная, но она является произведением перестановок простых нижних треугольных матриц. Вследствие этого равенство (У=МА называется иногда треугольным разложением матрицы А, Следует подчеркнуть, что ничего нового здесь не было введено.
Треугольная факторизация есть просто гауссово исключение, выраженное в матричных обозначениях. Треугольные множители могут быть также вычислены другими алгоритмами; см. Форсайт, Малер (1967). Диагональные элементы матрицы (У называются ведущими влгментами; Ь-й ведущий элемент есть коэффициент при й-м з.ь линеяныв системы с хеАнимыми мзтгиплми 4д неизвестном в й-м уравнении на я-м шаге исключения. В нашем примере ведущими элементами являются 10, 2,5 и 6,2, И вычисление множителей, и обратная подстановка требуют деления на ведущие элементы. Поэтому алгоритм не может быть проведен, если какой-либо из ведущих элементов равен нулю.
Интуиция должна подсказать нам, что продолжать вычисления, в ходе которых некоторые из ведущих элементов оказались близки к нулю,— не слишком удачная идея. Чтобы убедиться в этом, изменим слегка наш пример так, чтобы получилась система — 3 2.099 6 к., = 3.901 Значение 2.000 элемента (2,2) матрицы изменено на 2.099, и правая часть также была изменена„чтобы вектор (О, — 1, 1) по- прежнему был точным ответом. Предположим, что решение вычисляется на гипотетической машине, которая имеет десятичную плавающую арифметику с пятью значащими цифрами.
Первый шаг исключения дает 0 — 1.Ох!О-' 6 х, = 6.001 Элемент (2,2) теперь очень мал в сравнении с прочими элементами матрицы. Тем не менее закончим исключение, не прибегая к перестановкам. Следующий шаг требует прибавления второго уравнения, умноженного на 2.5х10', к третьему. Для правых частей это означает, что нужно умножить 6.001 на 2.5х10'. Результат равен 1.50025х!0', и он не может быть представлен точно в нашей гипотетической плавающей системе. Он должен быть или усечен до 1.5002х10* или округлен до 1.5003х10'. Полученное число прибавляется затем к 2.5. Предположим, что используется арифметика с усечением. Тогда последнее уравнение принимает вид 1.5005 ~: 10'х, =1.5004 х! 0', откуда 1.0004 Х 10~ Поскольку точный ответ — это х,=!, кажется, что мы не сделали серьезных ошибок.
К сожалению, однако, х, приходится определять из уравнении — 1.0 х 10 'х,+ (6) (0.99993) = — 6.001, 3. линейные системы уР»ененип 50 что дает ~лв К 1О-» Наконец, х, определяется из первого уравнения 10х,+ ( — 7) ( — 1.5) = — 7, и зто приводит к х,= — 0.35. Вместо (О, — 1,1) мы получили ( — 0.35, — 1,50, 0.99993), Где же была допущена ошибка? Здесь не было «накопления ошибок округления», вызываемого выполнением тысяч арифметических операций. Матрица не близка к вырожденной. Трудности возникли в связи с выбором малого ведущего элемента на втором шаге исключения. Как результат этого множитель равен 2.5х10», и последнее уравнение имеет коэффициенты, превосходящие в 1О' раз по величине коэффициенты исходной задачи.