Дж. Райе МАТРИЧНЫЕ ВЫЧИСЛЕНИЯ И МАТЕМАТИЧЕСКОЕ ОБЕСПЕЧЕНИЕ Перевод с английского О. Б. Арушаняна под редакцией В. В. Воеводина МОСКВА «МИР» 1984 MATRIX COMPUTATIONS AND MATHEMATICAL SOFTWARE John R. Rice Professor of Mathematics and Computer Science Purdue University McGraw-Hill Book Company New York St. Louis San Francisco Auckland Bogota Hamburg Johannesburg London Madrid Mexico Montreal New Delhi Panama Paris Sao Paulo Singapore Sydney Tokyo Toronto 1981 Предисловие редактора перевода Известно, что в настоящее время стоимость математического обеспечения вычислительных машин составляет существенную часть общей стоимости вычислительных систем. Это соотношение, по-видимому, не только сохранится в будущем, но, скорее всего, изменится в сторону еще большего увеличения относительной стоимости математического обеспечения. Очевидно, что в этих условиях высокое качество становится одной из важнейших ха- рактеристик математического обеспечения, особенно с точки зре- ния гарантии достоверности его работы. Математическое обеспечение создается специалистами, которых почему-то принято называть программистами, хотя, как правило, они имеют специальное математическое образование. Но если роль математического обеспечения постоянно возрастает, то этого никак нельзя сказать в целом о положении программиста в науч- ном обществе. К сожалению, до сих пор для специалиста-матема- тика считается престижнее доказать какую-нибудь теорему и построить какой-нибудь численный метод на бумаге, чем сделать и внедрить в практику хорошую программу. Причин подобного отношения к труду программистов довольно много, но главная заключается в устойчивой недооценке и даже непонимании тех трудностей, с которыми приходится сталкиваться при создании высококачественного математического обеспечения. Вниманию читателя предлагается книга Дж. Раиса «Матричные вычисления и математическое обеспечение». В ней автор попытался на примере достаточно простых задач, связанных с решением си- стем линейных алгебраических уравнений и проблемы среднеквад- ратичного приближения, рассказать и детально проанализировать особенности и трудности создания соответствующего математиче- ского обеспечения. Эту книгу не стоит читать «по диагонали» — пользы от такого чтения будет немного. Внешне материал выглядит достаточно простым, математический аппарат хорошо известен, хотя некото- рые аспекты, связанные с оценкой достоверности получаемых ББК 22.19 Р18 УДК 512.8+518.12 Райе Дж. Р 18 Матричные вычисления и математическое обеспечение: Пер. с англ.—М.: Мир, 1984.—264 с., ил. В книге рассматриваются этапы конструирования практических вычислительных алгоритмов на примере решения систем линейных уравнений. Материал изложен просто и понятно. Приводится описание нескольких пакетов и библиотек программ, созданных в последние годы в США н Великобритании и нашедших широкое прак- тическое применение- Для математиков-прикладников, программистов, студентов и аспирантов университетов, 1702070000-145 ББК 22.19 р 041 (01)-84 39"84 ч- 1- 518 Редакция литературы по математическим наукам Джон Райе МАТРИЧНЫЕ ВЫЧИСЛЕНИЯ И МАТЕМАТИЧЕСКОЕ ОБЕСПЕЧЕНИЕ Научи, ред. К. Г. Батаев, И. А. Маховая. Мл. научн. ред. Т. А. Денисова. Художник В. В. Кошмин. Художественный редактор В. И, Шаповалов. Технический редактор М. А. Страшнова. Корректор С. А. Денисова ИБ № 3865 Сдано в набор 25.07.83. Подписано к печати 29.12.83. Формат 60Х901/,,. Бумага типог- рафская № 2. Гарнитура литературная. Печать высокая. Усл. печ. л. 16,50. Усл. кр.-отт. 16,50. Уч.-изд. л. 16,63. Изд. № 1/2913. Тираж 14000 экз. Зак. 2002. Цена 1 р. 40 к. ИЗДАТЕЛЬСТВО «МИР» 129820, Москва, И-110, ГСП, 1-й Рижский пер., 2 Ордена Октябрьской Революции и ордена Трудового Красного Знамени Первая Образцовая типография имени А. А. Жданова Союзполиграфпрома при Государственном комитете СССР по делам издательств, полиграфии и книжной торговли. Москва, М-54, Валовая, 28 1981 by McGraw-Hill, Inc. Перевод на русский язык, «Мир», 1984 Предисловие Эта книга объединяет два очень разных, но дополняющих друг друга предмета: матричное исчисление и математическое обеспе- чение 1). Решение линейных уравнений представляет собой типич- ный образец численных расчетов, которыми занимались еще в древности. Это основной «строительный блок» для алгоритмов решения задач науки, техники, промышленности и любой другой области человеческой деятельности, в которой используются мате- матические модели. Математическое обеспечение стало наиболее распространенным средством проведения численных расчетов вслед- ствие значительного повышения надежности и существенного сни- жения себестоимости его компонентов. Матричное исчисление может послужить идеальным введением для изучающих математическое обеспечение, поскольку оно является ясной, простой (но в то же время не слишком простой) областью знаний, полезной для всех технических приложений. Цель той части книги, которая посвящена матричному исчис- лению, заключается в том, чтобы представить основные алгоритмы решения линейных систем и расчетов по методу наименьших квад- ратов (статистической регрессии). Более сложные задачи вычис- ления собственных значений и сингулярного разложения здесь опущены. Однако мы непосредственно сталкиваемся с наиболее трудной проблемой, возникающей в численных расчетах: как уз- нать, что получены верные ответы? И опять матричное исчисление оказывается идеальным средством для исследования этой централь- ной проблемы. Особое значение придается пяти аспектам математического обеспечения: учету человеческого фактора, оценке рабочих харак- теристик математического обеспечения, технологии создания его компонентов, конструированию алгоритмов и обучению студентов навыкам разработки программ. Эти аспекты относятся к общим вопросам проектирования математического обеспечения, и здесь также матричное исчисление оказывается очень полезным. Оценку рабочих характеристик математического обеспечения легче всего проводить на известных стандартных задачах; задачи 1) В данной книге под математическим обеспечением понимаются программные средства решения задач вычислительной математики, а также математической статистики.— Прим. перев, 6 Предисловие редактора перевода результатов, могут быть интересны и квалифицированному специа- листу. Главное заключается в другом — именно на этих хорошо известных в математике задачах показано, как трудно создать хорошую программу даже в той области, где кажется, что уже нет никаких проблем. В книге затрагивается очень много вопросов, связанных с раз- работкой высококачественного математического обеспечения. Мно- гие из них рассмотрены весьма схематично. Однако большим до- стоинством книги является то, что рассмотрен весь спектр основных вопросов. Эта книга, безусловно, будет интересна нашим читателям, особенно тем из них, кто хочет понять, почему до сих пор так мало хороших программ и как создавать хорошие программы. В. В. Воеводин Глава 1 ВВЕДЕНИЕ Цель изучения прикладной математики и вычислительных наук1) состоит главным образом в овладении средствами анализа матема- тических моделей, возникающих в хозяйственной и организацион- ной сферах деятельности. В этой книге рассматриваются вопросы анализа моделей простейшего вида, основанных на системах ли- нейных уравнений. Матричное исчисление изучается здесь не только с алгоритмической точки зрения, но и с точки зрения раз- работки математического обеспечения. Таким образом, мы будем исследовать как алгоритмические средства, так и вопрос о том, как сделать эти средства легко доступными для использования. Прежде чем перейти к изложению этого материала, кратко обрисуем роль теории, методов и математического обеспечения в решении рас- сматриваемых задач. Математические модели. Математические модели бывают двух разновидностей: статические (или аналитические) и динамические. Статическая модель представляет собой совокупность уравнений, формул, определений, таблиц, соотношений и данных, которые математически описывают ситуацию или явление в предположении об их завершенности. Динамическая модель является совокупностью тех же объектов, которые описывают, как изменяется ситуация от одного состояния к другому. В динамических моделях существен- ным параметром является протяженность процесса, и они приме- няются итерационно, начиная с некоторого начального состояния с последующим наблюдением за тем, что происходит на более позднем этапе. Не все, но большинство динамических моделей достаточно явно описывают переход от состояния к состоянию, и такой подход к построению динамических моделей часто легко сопоставить с анализом, требуемым для получения той же инфор- мации на основе статической модели (если она существует). Цель значительной части математического и научного анализа состоит в преобразовании заданной модели в явную модель: ответ == формула (данные, параметры, переменные). К сожалению, это часто невозможно сделать, и некоторые простей- шие статические модели определяют ответ А неявно. Простые 1) В оригинале — computer science.— Прим, перев. 8 Предисловие матричного исчисления как раз подходят для этой цели. Однако даже в этой области нелегко дать такого рода оценку, поскольку цели, преследуемые при создании математического обеспечения, часто оказываются противоречивыми, а различные виды затрат трудно сопоставлять ввиду отсутствия единого критерия. В. ис- тории развития математического обеспечения было гораздо больше неудач, чем успехов. Тем не менее мы должны проявлять настой- чивость — слишком высоки затраты кустарного производства. Кроме того, использование компонент математического обеспечения делает проверенные на практике идеи и опыт специалистов доступ- ными даже для тех, кто недостаточно знаком со специальными приемами численных расчетов. Учебные проекты, приведенные в гл. 12, преследуют три цели: (1) практика в использовании математического обеспечения, соз- данного специалистами, (2) приобретение опыта в практических вопросах проектирования математического обеспечения и кон- струирования алгоритмов, (3) ознакомление с более значительными программными комплексами. По каждому проекту учащимся сле- дует написать хорошо организованный и аргументированный от- чет. Рекомендуется большинство проектов давать на группу из двух-четырех студентов; это позволит им не только приобрести дополнительный опыт коллективной работы, но и участвовать в существенно более глубоких исследованиях. Предполагается, что читатели знакомы с линейной алгеброй или теорией матриц в объеме вводного курса и имеют солидную подготовку в программировании — предпочтительно на Фортране, так как большая часть существующего математического обеспечения написана на этом языке. Книга предназначена для студентов млад- ших и средних курсов, однако старшекурсникам и аспирантам она также может быть полезной. При более высокой подготовке учащиеся могут, например, более глубоко проанализировать пред- лагаемые проекты. Общий объем материала достаточен для курса длительностью в полсеместра или семестр. Книга пригодна также в качестве под- готовительного этапа для более широкого или более специализи- рованного курса. Она основана на курсе лекций, в котором часть материала читалась в течение 8 недель, а затем лектор переходил к линейному программированию и нелинейной оптимизации, Я благодарю Карла де Бура, Ричарда Хэнсона, Роберта Линча и Гарольда Стоуна за ценные предложения как по деталям, так и по организации книги в целом. Я также благодарю Джона Рейда, предоставившего примеры разреженных матриц для четвертой главы. Разрешения на использование материала, включенного в Приложение, получены от издательств Academic Press, ACM и ^1АМ. Джон Р. Райе 10 Глава 1_______________________________•____ примеры (В и С — исходные данные): А2 + A cos A+6.1= 4.2 ^В^ТТ—^-, ^=2A(t)sint-B+-c^)-, А(0)=1.0. Заметим, что второй пример определяет статическую модель, ко- торая в действительности является известным математическим эквивалентом динамической модели. С такого рода моделями поступают двояким образом. Первый путь состоит в решении задачи, т. е. проводятся различного рода манипуляции и применяются аппроксимации до тех пор, пока ответ не будет получен в явном виде. Такие манипуляции могут исполь- зовать конкретные значения данных и параметров, и поэтому их результатом не является построение явной математической модели; вместо нее строится частное решение. Это обычно происходит при применении численных методов. Второй путь заключается в оптимизации. Некоторые параметры или переменные являются свободными, и их следует подобрать таким образом, чтобы получить «наилучшее» значение для неко- торой другой переменной. Приведем простой пример: X^H^R—R^T^R-^—min (О, 1.0—1.5R*H). Здесь Х — жалование, выдаваемое в конечном счете служащему на руки. Предположим, что мы можем управлять параметром R, который представляет собой номинальную ставку заработной платы без учета различного рода добавок или вычетов. Можно применить дифференциальное исчисление, чтобы получить уравнение, опреде- ляющее оптимальное значение R: dX/dR=H—2R+T*R——[min(0, 1.0—1.5R»H)]=0. Значение R, которое является решением этого уравнения, должно определять наибольшее или наименьшее итоговое жалование X. В данном случае легко видеть, что R определяет наибольшее зна- чение X. Отметим пять обстоятельств, имеющих отношение к процессу оптимизации: (1) Оптимизация явной аналитической модели при- водит к другой математической модели. Вторая модель редко бы- вает явной. (2) Почти всякая оптимизация выполняется посредством приравнивания производной нулю. Исключения сводятся к опре- деленным тривиальным случаям или к тем случаям, когда каждый из них может быть исследован по отдельности. Это обстоятельство часто имеет скрытый характер, так как математические объекты, которые должны быть подвергнуты дифференцированию, недиффе- ренцируемы в обычном смысле (вспомним «min» в формуле вычис- Введение 11 ления итогового жалования). (3) При оптимизации наблюдается тенденция к существенному увеличению сложности, и поэтому зачастую оптимизируемые модели менее сложны, чем решаемые. (4) Большинство методов анализа статических моделей применимо также к динамическим моделям. (5) Легко построить модели, ре- шение или оптимизация которых невероятно трудны. Существуют примеры моделей исключительной важности, полный расчет кото- рых (известными на настоящий момент методами) потребовал бы использования всех потенциальных ресурсов страны (человеческих и других) в целом в течение многих лет. Например, что произойдет, если какие-либо компоненты атомной станции выйдут из строя, в результате чего тепловыделяющие стержни расплавятся? Или какое влияние на экономику будет иметь сокращение налогов на 10 миллиардов долларов? В действительности модели упрощаются таким образом, чтобы они стали доступными для обработки. Однако при этом упрощения не должны приводить к потере связи с реальностью. Часто бывает очень трудно обеспечить эту связь; трудно даже узнать, потеряна она или нет. Средства анализа. Данная книга не касается моделирования как такового. Мы интересуемся средствами решения и анализа моделей. Эти средства распадаются на три основные категории: Теория дает общее понимание модели и процесса решения за- дачи. Это один из путей достижения качественного представления о том, что происходит в действительности. Методы дают средства решения задачи. Метод есть совокуп- ность указаний или шагов решения задачи. Иногда (возможно, часто) можно решить задачу или исследовать модель, не достигнув полного понимания метода или самой задачи. Математическое обеспечение — это законченный метод, вопло- щенный в программе для ЭВМ. В самом лучшем случае достаточно просто нажать кнопку РЕШИТЬ, чтобы получить ответ. Глава 2 ОСНОВНЫЕ ПОНЯТИЯ ЛИНЕЙНОЙ АЛГЕБРЫ 2.А. СПИСОК ПОНЯТИЙ, КОТОРЫЕ ДОЛЖНЫ ЗНАТЬ СТУДЕНТЫ Существуют три распространенных источника основных поня- тий линейной алгебры, относящихся к матричному исчислению: элементарный курс линейной алгебры, вводный курс в численный анализ и введение в линейный анализ. Обычно курсы линейной алгебры наименее удовлетворительны для наших целей, поскольку они делают упор на алгебраические аспекты предмета, а не на его анализ, курсы линейного анализа обеспечивают наилучшие пред- посылки. Большинство студентов ощущают, что их подготовка в линейной алгебре ниже того уровня, который необходим для удов- летворительного (но неглубокого) понимания матричного исчис- ления. Материал этой главы изложен преднамеренно сжато в пред- положении, что читатель имеет возможность пользоваться справоч- никами. 2.А.1. Векторы Векторы представляют собой направленные отрезки (они имеют длину, направление и положение) в М-мерном пространстве. Они рассматриваются как векторы-столбцы, если не оговорено про- тивное: [У1\ У2 у= • . .Ум. Операция транспонирования, которая меняет столбцы на строки и наоборот, обозначается верхним индексом Т. Векторы обычно выражаются через базис: выбирается стандартная совокупность векторов bi, ba, . . ., bn, и все другие векторы выражаются через базис bi, i==l, 2, . . ., N, y=yibi+y2ba+. . .+УмЬк- Коэффициенты yi в этом представлении называются компонентами _______________ ____Основные понятия линейной алгебры 13 вектора у, а само представление записывается в компактной форме У=(У1> У2, • . ., Ум)7. Базисные векторы определяют систему координат, и компоненты У) являются координатами точки конца вектора в этой системе координат. Стандартные арифметические операции (для векторов х, у, z и скаляра а): Сложение: х+у=у+х= (х,+у,, х^+Уа, ..., Хм+Уы)"1', х—у=—(у—х); (x+y)+z=x+(y4-z). Умножение на скаляр: ax==(aXi, ax;,, ..., ахм), а(х+у)==ах+ау. Совокупность векторов Xi, Хг, .... Хщ линейно независима, т если ни одна их линейная комбинация ^ KiXi не равна нулю, за исключением нулевой комбинации, т. е. m ^ о1[Х)=0 означает, что ОС[=0 для всех i. Пространство натянуто на совокупность векторов Xi, Ха, . . ., Хщ, если каждый вектор пространства может быть выражен через линейную комбинацию векторов Xi, Xs, . . ., Xm. Размерность векторного пространства есть минимальное число из этих векторов, требуемых для получения пространства; каждый базис М-мерного пространства должен иметь в своем составе N векторов. N-мерные векторные пространства для краткости часто называются ^-про- странствами, и ^-векторы являются векторами в N-пространстве. Скалярное, или внутреннее, произведение двух векторов х и у имеет вид N х^у^х-у^^ XiYi, где Xi, yi — компоненты векторов х и у. Два вектора ортогональ- ны (перпендикулярны), если x^y=Q. Размер вектора может быть измерен нормой |)xJ, где [|х||2__„Т„ __ У у 2 || Х || — Х Х — ^J X, . В случае трех измерений это — обычная евклидова длина. Двой- ные вертикальные черточки обозначают знак нормы. Часто бывают удобны две другие нормы: ||х||„=тах|х,|, ! Mx-Sxi. 14 Глава 2________________________________________ Евклидова норма обозначается IMa^ES^T72' ^гол ^ между двумя векторами определяется соотношением x^v с0 1|х|И1у|Ь • Формат записи векторов должен быть согласован с форматом записи матриц, поэтому векторы обычно записываются как векторы- столбцы, т. е. вектор х==(2,1,4,—2)7 в действительности является вектором f 2\ х= 4 • _9 Это усложняет изложение текста, поэтому мы будем записывать векторы горизонтально со знаком транспонирования Т, если ука- зание формата столбца не является необходимым. Иногда мы будем использовать также векторы-строки, которые представляют собой векторы, чья запись в матричном формате имеет на самом деле горизонтальный вид, например строка матрицы. Мы ввели векторы как абстрактные объекты и затем перешли к их конкретному представлению с использованием базиса. Базис, который определяет систему координат, существен для вычислений (поскольку вычислительные машины не манипулируют с векторами непосредственно), однако абстрактное геометрическое рассмотре- ние часто более полезно для понимания того, что происходит в процессе вычислений. Задачи подраздела 2. A.I 1. Вычислить [|x||i, ЦхЦг и ||х||<„ для каждого из следующих векторов: (а) (1, 2, 3, 4)'1, (б) (О, —1, 0, -2, О, 1)', (в) (О, 1, —2, 3, —4)', (г) (4.1, —3.2, 8.6, —1.5, —2.5)1. 2. Вычислить угол в между следующими парами векторов: (а) (1, 1, 1, I)7, (—1, 1, -1, 1)Т, (б) (1, 0, 2, O)T, (0, 1, —1, 2)Т, (в) (1, 2, —1, 3)Т, (2, 4, —2, 6)Т, (г) (1.1, 2.3, —4.7, 2.0)Т, (3.2, 1.2, —2.3, —4.7)"1. 3. Найти евклидову норму вектора (2, —1, 4) и скалярное произведение век- торов (2, —3, 4, D'r и (4, —3, 2)1'. 4. Определить, принадлежит ли вектор (6, 1, —6, 2)T пространству, натяну- тому на векторы (1, 1, —1, I)'1, (—1, 0, 1, l)1 и (1, —1, —1, О)'1'. Какова размер- ность пространства, натянутого на эти три вектора? 5. Образуют ли векторы bi=(l, 0, 0, —l)1, Ьг=(0, —1, 1, 0)1', ba=(0, 0, —1, I)1', b^(l, 1, О, О)'1' ______________________Основные понятия линейной алгебры 15 базис для векторов размерности 4? 6. Доказать, что три вектора bi=(l, 2, 3, 4)Т, b,=(2, 1, О, 4)Т, Ьз=(0, 1, 1, 4)Т линейно независимы. Образуют ли они базис трехмерного пространства? 7. Предположим, что Хд, х^, . . . , Хщ линейно независимы, а х;, Хз, ... . . ., Хд, Хд+i линейно зависимы. Показать, что Xm+i является линейной комби- нацией Xi, Х2, . . . , Хщ. 8. Пусть Xi=(l, 2, 1), Х2==(1, 2, 3), х,э=(3, 6, 5). Показать, что эти векторы линейно зависимы, но на них может быть натянуто двумерное пространство. Найти базис этого пространства. 9. Показать, что для любых двух векторов х и у и любой из трех норм (1, 2 или оо) имеет место неравенство 11|х||-||у|||<||х-у||. Это неравенство выполняется для всех векторных норм. 10. Показать, что для двух ненулевых векторов х и у равенство ||х+у||а= ==(|х||2+|1у||2 выполняется тогда и только тогда, когда у=ах, где а — неотрица- тельная константа. 11. Показать, что для двух п-векторов х и у выполняется неравенство ||х+у|1<:]|х1]+||у|| для любой из трех норм (1, 2 и оо). Указание: для норм 1 и оо докажите неравенство для двумерных векторов и затем примените индукцию. 2. А. 2. Матрицы Абстрактные объекты, которые приводят к матрицам, представ- ляют собой линейные отображения, преобразования, или функции, определенные над векторами. Поскольку для векторов были вве- дены координаты, то эти линейные функции могут быть конкретно представлены посредством двумерной таблицы чисел, т. е. матри- цей. Например: /I 6-2N А== 4 17 -12 = (а,,). \0 42 6.1 / Если у есть линейная функция от х, то каждая компонента у^ вектора у есть линейная функция компонент Xj вектора х. Таким образом, для каждого k мы имеем yk=akiXi+ak2X2+. . ,+akNXN. Коэффициенты собираются в матрицу А, которая и представ- ляет линейную функцию (отображение, преобразование, или соот- ношение) между векторами х и у. Линейная функция обозна- чается Ах. Действия над матрицами определяются правилами, выполня- емыми для линейных отображений. Так, А+В есть представление суммы двух линейных функций, представленных в свою очередь 16 Глава 2______________________________________ матрицами А и В. Имеем А+В==С, где Ci)=aij+bi]. Произведение АВ представляет собой эффект применения отображения В и по- следующего применения отображения А. Покажем, что АВ==С, где Cij==2aikbki- Для этого обозначим у=Ах и z==By. Следова- тельно, нам необходимо определить матрицу С такую, что z==Cx. Выразим это соотношение на языке компонент: N N Yk== S akjX„ г,==2Ь,ьУь. j=l k=l Таким образом, N / N \ N / N \ 2i= 2 bik[ 2 3k,x, 1=2(2 b^k, )xj k=l \j=l /' )=1 \k=l ' ! ' и поэтому сц определяется указанной выше формулой. Элемент, стоящий на пересечении i-й строки и j-ro столбца матрицы С, есть скалярное произведение i-й строки матрицы А и j-ro столбца мат- рицы В. Справедливы арифметические правила А+В=В+А, АВ=т^=ВА, за исключением специальных случаев. Транспонированная матрица А'1 получается отражением мат- рицы А относительно ее диагонали (элементов а;,), т. е. a^==aji. Единичная матрица имеет единицы на диагонали, а остальные ее элементы равны нулю: /10 (Л 1= 0 1 0 . \0 О I/ Единичная матрица обязательно квадратная, т. е. имеет одина- • ковое число строк и столбцов. Легко видеть, что 1А==А1=А. Об- ратная матрица А"1 есть такая матрица, для которой выполняется равенство А^А^. Не всякая матрица имеет обратную. Дейст- вительно, произведение АВ может быть равно нулю, даже если А или В не были нулевыми матрицами: /I iy i iy/о ел =1 i 1=1 1-норма и оо-норма широко применяются, поскольку они легко вычисляются. Легко показать, что р (A)s?S|]A|| для любой нормы, поэтому 1-норма и оо-норма обеспечивают простую оценку спект- рального радиуса матрицы. Интерес к матрицам возникает из их связи с линейными функ- циями нескольких переменных (такие функции представляют со- бой естественный способ описания простых математических моде- лей объектов, зависящих от нескольких переменных). Линейные функции со многими переменными «существуют» как абстрактные функции точно так же, как векторы «существуют» как абстрактные Линейная (рункция со i i многими переменными .y=F(x) ^ "' ki, / -^ /^^ / > Векторы Векторы Рис. 2.1. Диаграмма эффекта функции со многими переменными. объекты. Для векторов мы определили понятия, с которыми мы можем манипулировать и делать вычисления с помощью введения систем координат; то же самое имеет место и для функций со мно- гими переменными, как показано на рис. 2.1. Если мы выразим х и у в обозначениях некоторой системы координат, то должна существовать формула для вычисления координат (yi, уз, ... • • i Ук)1 через координаты x=(Xi, Ха, . . ., х^)7. Эта формула включает в себя двумерную таблицу чисел, скажем A=(aij), и на самом деле имеет вид у=Ах. Таким образом, матрица есть конкретное представление линейной функции со многими переменными, полученное от введения системы координат для векторов. Различные выборы систем координат 20 Глава 2___________________________________ дают различные матричные представления для одной и той же функции. Если В^^АС, то А и В называются подобными матрицами; эта формула выражает эффект изменения систем координат (или другого выбора базисных векторов). Поэтому все свойства мат- рицы А, которые являются свойствами лежащей в ее основе ли- нейной функции F, не меняются после применения подобного пре- образования B^C'^AC. Эти свойства относятся к собственным значениям и собственным векторам. Как пример полезности собственных значений и собственных векторов рассмотрим, что произойдет, если в качестве базисных векторов выбраны собственные векторы (предположим, что их достаточно для этой цели). В этом случае матричное представление линейной функции есть просто диагональная матрица! Далее, компоненты вектора Ах являются компонентами вектора х, умно- женными на соответствующие собственные значения. Поэтому анализ матриц и матричное исчисление существенно упрощаются, если получить собственные векторы матрицы А и использовать их в качестве базисных векторов. Задачи подраздела 2. А. 2 /I 2\ 1. Показать, что матрица А=1 | вырождена. 2. Найти ненулевое решение системы Xi—2x2+ Хз=0, Xi+ Xg—2хз=0. 3. Пусть 123 11 2 101 А=0-12], В= -1 1-1 , С=012). ^2 0 2-7 \ 1 0 2^ ^2 О К (а) Вычислить АВ и ВА и показать, что АВ^&ВА. (б) Найти (А+В)+С и А+(В+С). (в) Показать, что (АВ)С==А(ВС). (г) Показать, что (АВ)Т=ВТАТ. 4. Показать, что матрица А вырождена: 3 1 О А= 2 -1 -1 . \4 3 \) 5. Для матрицы А задачи 4 вычислить РА и АР, где Р — матрица переста- новок: 010 fioo). \о о \' Основные понятия линейной алгебры 21 6. Найти матрицу перестановок Pi, которая для квадратной матрицы А чет- вертого порядка переставляет второй и четвертый столбцы. Найти Рд, которая переставляет первую и третью строки матрицы А. 7. Напишите следующую систему линейных уравнений в матричной форме Ах=Ь: 2xi+X2—Зхз+4х4—2=0, Х1+Х2—Хз=1, 2х2+Х4—7хз=8, ' l+Xl+X2+Xg+X4=5. 8. Показать, что определенная выше матричная норма может быть также вы- числена как || А |]= max ^. II х И Ф О II Х II 9. Матрица размеров NX 1 является вектором-столбцом- Показать, что ре- зультат транспонирования вектора-столбца есть вектор-строка. 10. Пусть В — диагональная матрица. Показать, что матрица ВА вычисляется умножением первой строки матрицы А на Ьц, второй строки — на Ьзг и т. д. Как вычислить АВ? 11. Проверить следующие системы уравнений на совместность (т. е. показать существование по крайней мере одного решения): (а) х+ y+2z+ w=5, (б) х+ У+ z+ w=0, 2х+3у— z—2w=2, x+3y+2z+4w=0, 4x+5y+3z =7, 2х + z— w=0. 12. Показать, что матрица -Q не вырождена, вычислив обратную к ней матрицу. 13. Вычислить обратные матрицы для /I 1\ /I 0\ Чо.)- ^(i.) и для АВ. Показать, что A^B-S^AB)-1, но B-'A-^AB)-!. Объяснить это, используя определение произведения матриц, представляющих композицию двух линейных функций. 14. Найти одно решение неоднородной системы: Х1+2хз=3, 2xi+4x2=6, и затем решение соответствующей однородной системы: Xi+2x2=0, 2xi+4X2=0. Показать, что еще одно решение неоднородной системы может быть получено при- бавлением к уже найденному ее решению решения однородной системы, умножен- ного на константу. 15. Составьте программу, выполняющую сложение или умножение матриц. Входными параметрами программы являются две квадратные матрицы А и В порядка N, а также переключатель MODE, определяющий вид работы: MODE=1 22 Глава 2________________________________________ для сложения и MODE=2 для умножения. Выходным параметром должна быть матрица С=А+В или С=АВ. 16. Найти число операций сложения и умножения, необходимых для умноже- ния NXN-матрицы на М-вектор и перемножения двух NxN-матриц. 17. Вычислить ||A||i и ]|А||«, для матрицы 1\ 0 2 —П _ 6 —4 3 О А= 4 0—4 2 " Л 5 1 6, 18. Найти ненулевую ЗХЗ-матрицу А и ненулевой вектор х, такие, что Ах==0. 19. Пусть А=1—2хх1, где х—вектор-столбец. Показать, что А—ортого- нальная матрица и А2^!. 20. Для невырожденной матрицы А показать, что х^^Ах^О тогда и только тогда, когда х1^1^. 21. Пусть х — вектор-строка. Показать, что xxт=x*x=(||x||a)2. Чему равно х^? 22. Пусть для любого вектора х функция F определена одним из следующих способов: (а) Р(х)=(хз, xi)1; (б) F (x)=(xi+X2, 0, Хз)Т; (В) F(X)=(X2, 0, Xi—X2, Хз, О)1. Показать, что каждый из этих способов определяет F как линейную функцию, 23. Найти для каждой из линейных функций, определенных в задаче 22, пред- ставляющую ее матрицу. 24. Рассмотрим матрицу /cos9 —sin6\ Л -1 I \sin 9 cos9/ Дать геометрическую интерпретацию функции, представленной матрицей А. По- казать, что А — ортогональна. 25. Описать геометрически эффект линейных функций, представленных сле- дующими матрицами: (а) / cos9 0 sin9\ (б) /5 0 0\ ( 0 1 0 ) , 0 5 О) , \—sin9 0 cos6/ \0 0 5/ (в) Л 0 0. (г) ,-\ О (К (о о о), (01 °)- V) О I/ \- О О —I-7 26. Пусть А — ортогональная матрица и Ах==Хх, где х=,&0. Показать, что 1^1=1. 27. Пусть вектор х имеет компоненты (xi, Хг)1 в обычной эвклидовой системе координат на плоскости [т. е. базисными векторами являются (1, О)1 и (О, I)7]. Какие компоненты имеет вектор х в базисе bi=(l, I)1, ba==(l, —l)1^? Дать геомет- рическую интерпретацию соотношения между системой координат, определенной базисом bi и Ьа, и обычной системой координат. Основные понятия линейной алгебры 23 28. Показать, что матрицы А и В равны тогда и только тогда, когда Ах=Вх для всех векторов х. 29. Пусть матрица А имеет самое большее один линейно независимый вектор- столбец. Показать, что тогда матрица А имеет только одну линейно независимую вектор-строку. Показать, далее, что А может быть записана в виде ху1 для неко- торых векторов-столбцов х и у. 30. Показать, что |]АВ||<:ДА|| ||В(| для матричной нормы, определенной в этой главе. Это соотношение имеет место не для всех матричных норм. 31. Доказать правильность данной формулы для нормы |[А||„. 32. Доказать правильность данной формулы для нормы ||A||i. 33. Показать для NXN-матрицы А, что тах | а» | ^ [| А |] <; N тах | а» I. 1,1 i, ] 34. Рассмотреть , 6 13 —17. /6—4 —1\ А=( 13 29 3S\, A-i=|-4 11 7 . \—17 —38 SO^ \—1 7 5/ (а) Описать геометрически множество S={Ax| ЦхД^!}, которое является об- разом единичной сферы при линейном отображении, представленном матрицей А. (б) Что представляют собой ||А||„, ЦА-Ч!.», ||А||] и ||A-4[i? 35. Пусть в правиле Крамера каждый определитель вычисляется как сумма (п—I)! произведений, каждое из которых состоит из п множителей. Оценить, сколь- ко тогда потребуется операций умножения, чтобы решить Ах=Ь для п=10, 100 и 1000. Сколько потребуется на это решение машинного времени, если каждая операция умножения (и соответствующие другие операции) выполняется за 1 микросекунду? 36. Повторить решение задачи 35 при условии, что каждый определитель вы- числяется за гт^/З операций умножения. 37. Показать, что если строки и столбцы матрицы А линейно независимы, то А — квадратная невырожденная матрица. 38. Пусть А — представляет линейную функцию F [т. е. y=F(x) может быть описана как у=Ах]. Пусть изменения координат в каждом из двух векторных пространств выполняются так, что координаты (wi, w^, . . . , Wn)1 вектора х в новой системе координат связаны с координатами (xi, Хз, . . . , Хп)7 соотношением w=Bx и прдобным же образом новые координаты (Zi, z^, . . . , Zn,)7 вектора у вычисляются посредством соотношения z==Cy. Показать, чтоСАВ"! есть матрица, которая представляет функцию F в двух новых координатных системах. 39. Пусть D — диагональная матрица. Показать, что DA=AD для всех мат- риц А тогда и только тогда, когда D=al для некоторой константы а. 40. Определим е* разложением I+A+A2/2!+AЗ/ЗI+A4/4!+. . . . Доказать следующие свойства матричной экспоненциальной функции'. (а) iHl^ellAll»; (б) еЛ+и^е^^ тогда и только тогда, когда АВ=ВА; (в) e^-^I; (г) B=eAt удовлетворяет уравнению —.—=АВ. Глава 3 КЛАССЫ И ПРОИСХОЖДЕНИЕ ЗАДАЧ МАТРИЧНОГО ИСЧИСЛЕНИЯ З.А. СИСТЕМЫ ЛИНЕЙНЫХ УРАВНЕНИЙ Ах=Ь З.А.1. Модели физических систем Существуют многочисленные физические задачи, которые ес- тественным образом сводятся к линейным уравнениям. Рассмот- рим, например, баланс сил, действующих на какую-либо конст- рукцию (мост, строение, каркас самолета), как это показано на Рис. 3.1. рис. 3.1. Силы веса моста и грузовика сбалансированы силами, приложенными в местах опоры моста. Эти силы распространяются вдоль железных балок, и в каждой узловой точке их равнодейст- вующая должна равняться нулю (в противном случае мост должен прийти в движение). Если разложить силы на горизонтальные (х) и вертикальные (у) составляющие, то в каждом узле будем иметь два уравнения: сумма горизонтальных (х) составляющих сил=0, сумма вертикальных (у) составляющих сил=0. Эти уравнения объединяются в одну большую линейную си- стему, которая может быть решена, если нужно вычислить неиз- вестные силы, действующие на различные балки. Некоторые из этих сил известны (веса балок и грузовика), и они образуют члены, которые переносятся в правую часть системы, так что система неоднородна. Для электрической цепи (рис. 3.2) на основе законов Кирх- гофа аналогично описывается баланс сил токов в узлах (точках разветвления) и напряжений по замкнутым контурам. Эти фор- мулы линейны для значений сопротивлений и источников энергии Классы задач матричного исчисления 25 и сводятся к системе линейных уравнений относительно всех не- известных величин. В решении такого типа задач выделяются три этапа: 1. Формулируется математическая модель системы (конструк- ции здания, электрической цепи и т. д.). Явно вычисляются а^, основанные на физических законах, которым подчиняется система. ^-—————vWV——————у\ У/^\\ fn -/ /^\^f ^С\ L 7 0\ ^-—w^—-^—MV—z-——и——-^ Рис. 3.2. 2. Описывается вынуждающая функция (внешние силы, дейст- вующие на конструкцию, или внешние нагрузки и подводимая мощность цепи). Явно вычисляются bj, основанные на анализе внешних условий, характеризующих конкретную ситуацию. 3. Решается линейная система Ах=Ь. Первые два этапа зависят от знания физического состояния исследуемого объекта (строительное дело, физика и т. д.). На третьем этапе используются алгоритмы численных расчетов; одна из целей математического обеспечения — предоставить надежные и эффективные средства для такого рода вычислений. Многие физические системы моделируются дифференциальными уравнениями, например: ^+cos выполнить: Решить для у(х) уравнение dгy/dx2+cos(x)yl-ly^log(x+4) и положить у'(х)=у(х). Если тах|у'(х)—у^^х)! достаточно мал, то выйти из цикла. Конец цикла. Согласно этому подходу, во внутреннем цикле вычислений исполь- зуется процедура решения линейного дифференциального урав- нения. Ясно, что процедура внутри цикла может быть выбрана более эффективной, чем многократное повторение решения линей- ной задачи. Здесь уместно отметить два момента: Классы задач матричного исчисления 27 1. Программное средство типа «решить Ах=Ь» может быть частью внутреннего цикла более общих вычислений и тем самым существенно, чтобы оно выполнялось эффективно. 2. То же самое программное средство может быть «скрыто» по отношению к решению общей проблемы, и важно, чтобы это сред- ство было надежным. Оно может быть скрыто так глубоко, что программист в действительности и не знает, что используется данное средство, и любая ошибка в нем очень трудно диагности- руема. З.А.2. Выравнивание данных методом наименьших квадратов или регрессия Пусть заданы 100 результатов наблюдений d(t|<) (данные в 100 точках) и предположим, что функция d(t) моделируется фор- мулой d(t)=cigi(t)+c,g2(t)+. . .+c,g,(t), где gj(t) — известные функции переменной t. Коэффициенты или параметры Ci, Сз, . . ., с, должны быть подобраны так, чтобы модель согласовывалась с данными в том или ином смысле. Теоретически имеем 7 d (4) = S c,g, (tJ для k=l, 2, ..., 100, но в силу ошибок в регистрации наблюдений или плохого предпо- ложения о характере модели или того и другого, вместе взятых, точное равенство в этих 100 уравнениях не достигается для любого выбора семи параметров Cj. Можно, однако, определить G] так, 100 чтобы сумма ^ ri квадратов невязок k= 1 rk=d(tk)-2c,g,(tk) i =1 была минимальна. Это означает, что параметры с; оптимизируются так, чтобы минимизировать среднеквадратическую ошибку. Набор параметров G], для которого 100 E(Ci, c„ .... с,)= ^rk2 k=! минимальна, определяется системой уравнений Ц-=0, i=l,2,..,7. 28 Глава 3_________________________ Имеем 100 Г 7 1 S d(tk)-Sc,g,(t„) gi(tk)=0, i=l, 2, .... 7 k=l [ j =1 J или, переписав эту систему, 2c,2g,(tk)g,(tk)=^d(tk)gi(t,), i==l,2, ...,7. 1=1 k=l k=l Эти выкладки приводят к системе уравнений Gc=d, где Q — матрица с элементами 100 gi,= s gi(tk) gj (tk)> d — вектор с компонентами 100 di== S d(t„)g,(tk), k=l с — вектор искомых параметров. Полученная система называется системой нормальных уравнений. З.А.З. Организационные модели Две рассмотренные физические системы (мосты и электрические цепи) представляют собой определенные формы организации фи- зических величин. Многие другие формы организации также при- водят к линейным моделям. Например, в экономике оплата счета имеет своим результатом увеличение суммы денег у одного лица и уменьшение — у другого. Деньги «делаются» различными пу- тями (открытием золотого рудника, рисованием картин, сооруже- нием письменного стола из пиломатериалов), и можно записать уравнения баланса потока денег (или эквивалентных форм богат- ства) для экономического региона (деревни, города, страны). Не- смотря на то что существует немало неопределенностей в реальном составлении балансов, экономисты создали много больших мате- матических моделей такого рода экономических процессов. В случае установившегося состояния (то есть в случае, когда с течением времени не происходит никаких изменений) эти модели часто пред- ставляют собой большие системы линейных уравнений. Подобная ситуация имеет место и для инвентарного баланса автомобильной компании. Компания владеет автомобилями, на- ходящимися на фабриках, складах, в пути и в розничной торговле. Существуют простые линейные соотношения, которые моделируют изменения в инвентаре. Несмотря на то что нет необходимости в Классы задач матричного исчисления 29 формулировке «решить задачу инвентаризации», существуют раз- нообразные манипуляции, которые нужно выполнять с этими линейными моделями. Например, компания нуждается в миними- зации своих затрат на складирование и может попытаться найти такое распределение автомобилей, внесенных в инвентарный спи- сок, которое обеспечивает это. Очевидное решение — поместить все автомобили в наиболее дешевое место для хранения продук- ции, но если таким местом окажется сама фабрика, то ни один автомобиль не будет продан. Поэтому определенные ограничения должны быть установлены для процесса минимизации. Это — при- мер задачи линейного программирования, предмета, близко примы- кающего к матричному исчислению. Формулирование, анализ и манипулирование с линейными системами уравнений представляют собой составные части процесса решения всех задач такого рода. З.Б. МАТРИЧНЫЕ УРАВНЕНИЯ АХ==В Матричное уравнение АХ=В является в действительности просто множеством линейных систем Axi=bi, Ах2=Ьз, . . ., Axni=bm для различных векторов решений xi и векторов правых частей bi. Матрица А в уравнении Ах=Ь обычно представляет модель (кон- струкцию моста, схему электрической цепи, регрессионную мо- дель), a b представляет вынуждающую функцию или внешние по отношению к модели объекты (нагрузки на мост и электрическую цепь, данные экспериментальных наблюдений). Поэтому для фик- сированной модели (А) может возникнуть необходимость в рас- смотрении многих различных воздействующих на модель объектов (b), то есть в решении матричного уравнения АХ=В. Существо дела состоит в том, что если решается линейная система с одной матрицей и многими правыми частями, то может иметь смысл пре- образовать матрицу А к специальной простой форме, с тем чтобы облегчить процесс решения этих многих линейных систем. Могут быть выделены следующие этапы решения этой задачи: 1. Вычислить матрицу А, моделирующую рассматриваемую систему. 2а. Вычислить вынуждающую функцию bi. За. Решить линейную систему Axi=bi. 26. Вычислить вынуждающую функцию Ьа. 36. Решить линейную систему Ах2=Ьг и т. д. Эти этапы могут быть реорганизованы следующим образом) 1. Вычислить матрицу А, моделирующую рассматриваемую систему. 30 Глава 3____ _____ 2. Вычислить вынуждающие функции bi, ba, . . . и сформировать матрицу В. 3. Решить матричное уравнение АХ=В, где Х — матрица со столб- цами Xi, Ха, . . . . Преимущество последнего подхода состоит в том, что на третьем этапе известно о предстоящих значительных вычислениях и они могут быть упорядочены для достижения большей эффективности. Первый подход может быть реализован так же эффективно, но это требует более тщательного программирования, чтобы избежать повторения работы на этапе 36, которая уже выполнена на этапе За. Хорошая программа решает матричное уравнение АХ==В эффективно, не вынуждая программиста вникать в детали про- цесса решения. З.В. ОБРАЩЕНИЕ МАТРИЦ А-1 В математических и технических текстах часто встречаются формулы вида: Х=А-Ч), y=B-l(I+2A)b. z^B-^A+IXC-^+A)^ которые, естественно, предполагают обращение матриц для вы- числения векторов х, у и z. Такой подход весьма неэффективен: почти никогда нет необходимости обращать матрицы, чтобы вы- числить другие объекты. В гл. 5 будет показано, что вычисление обратной матрицы требует почти в три раза больше операций и, быть может, в два раза больше памяти вычислительной машины, чем решение линейной системы уравнений. В частности, любое матричное выражение, примененное к вектору, может быть эф- фективно вычислено без обращения любой из входящих в это вы- ражение матрицы независимо от того, сколько обратных матриц имеется в выражении. В большинстве случаев явное обращение матрицы приводит к вычислительным потерям. Пусть решается матричное уравнение АХ=В. Вы можете ска- зать: «Ну, хорошо, непосредственное вычисление А~1 — дорого, но уж получив обратную матрицу, я быстро смогу найти все \i по формуле X)=A'"^bl». Даже эта аргументация неверна, поскольку существуют другие способы (например, LU-разложение матрицы А), которые гораздо дешевле для вычислений, чем А"1, и которые позволяют получить каждое х, из Ь) так же быстро, как и умноже- нием bi на А~1. Однако существуют некоторые задачи статистики и техники, в которых цель вычислений состоит в нахождении обратной мат- рицы как таковой без последующего ее использования для опре- деления чего-либо другого. В этих задачах элементы матрицы А"1 имеют смысл «влияющих параметров», которые показывают, как Классы задач матричного исчисления 31 вынуждающие объекты воздействуют на модель. Если нет необхо- димости в непосредственном исследовании элементов обратной матрицы, то весьма маловероятно, что следует вычислять обратную матрицу. З.Г. ОПРЕДЕЛИТЕЛИ det(A) Определители бесполезны для проведения вычислений в линейной алгебре. Исторически они были введены в математическую и на- учную терминологию в прошлом столетии, и сейчас они состав- ляют обязательную часть системы образования и математического языка. Определители нужны только в очень немногих разделах теории матричного исчисления, и редко возникает необходимость вычислить в явном виде значение определителя. Например, правило Крамера чрезвычайно неэффективно для решения линейных систем уравнений (см. задачу 35 гл. 2). Часто можно слышать: «Я не могу решить систему, поскольку определитель равен нулю». Однако может ли быть решена система или нет, обнаруживается почти для всех способов вычисления определителя прежде, чем мы получим значение определителя. Сама по себе величина определителя не имеет отношения к достоверности вычислений. Можно быть в полной безопасности, когда det (A)=10~48, и можно получить со- вершенно неправильные результаты, когда det (А)=6.2. Чтобы убедиться в этом, достаточно умножить каждое уравнение системы 20-го порядка на множитель, равный 100. В результате определи- тель системы увеличится в 1020 раз без какого-либо воздействия на вычислительные трудности или вырожденность системы. Типы матриц 33 Рис. 4.1. Пример расположения ненулевых элементов в двух больших разрежен- ных матрицах. Каждая точка указывает ненулевой элемент. На практике возни- кают гораздо большие разреженные матрицы, однако если их структуру умень- шить доданных размеров, то точки были бы не видны [пример любезно предостав- лен Джоном Рейдом (John Reid, A.E.R.E. Harwell)]. 4. А. 2. Упорядоченная разреженность Многие приложения приводят к матрицам, которые не только разрежены, но и имеют особую структуру расположения ненулевых элементов. Наиболее простые из них — ленточные матрицы, в которых ai,=0 для |i—Л>ширины ленты (рис. 4.2). Например, матрица конечных разностей в подразделе З.А.1 имеет ширину ленты, равную 1. Такие матрицы называются трехдиагональными. 2 № 2002 Типы матриц 35 4. Б. СПЕЦИАЛЬНЫЕ МАТРИЦЫ, ВОЗНИКАЮЩИЕ ИЗ АНАЛИЗА В методах для численных расчетов рассматриваются матрицы, чья специальная форма позволяет легче проводить вычисления. Так, например, цель метода исключения Гаусса — получить мат- рицу простой структуры, для которой уравнения могут быть легко решены. Рассмотрим некоторые специальные матрицы. 4.Б.1. Диагональные матрицы Для диагональных матриц ац==0 при всех i и j, кроме i==j. Собственные векторы интересны именно потому, что если они выбраны в качестве базиса (или координатных векторов), то пред- ставление соответствующей линейной функции со многими перемен- ными описывается диагональной матрицей. 4.Б.2. Треугольные матрицы Два вида треугольных матриц показаны на рис. 4.5. Они об- ладают тем свойством, что над ними легко выполнять различные действия. Линейные системы с треугольными матрицами могут быть \.0 а^О \ а,^0 \^ 9-ля к} Г~\\ ^ля '"У а Нижняя треугольная б Верхняя треугольная Рис. 4.5. решены прямой или обратной подстановкой. Собственные значения треугольной матрицы суть диагональные элементы. 4.Б.З. Ортогональные матрицы Если А—ортогональная матрица, то А'А=1 или t\\=h.~\. В этом случае легко решать линейные системы Ах=Ь, так как x=ATЬ. Ортогональные матрицц не меняют длину (ЦхЦг) любого вектора. Действительно, пусть у=Ах. Тогда |x|]^=xтx=x'1 Ix=x A^x^Ax^Ax^y'y^lyl. Это важно в численных расчетах, когда нежелательно, чтобы вы- полняемые действия увеличивали как ошибки, возникающие при задании параметров исходной задачи, так и ошибки округления, возникающие в процессе вычислений. 2» Типы матриц 37 7. Пусть А — симметричная матрица, а В — любая другая матрица. Пока- зать, что В^В — симметричная матрица. 8. Показать, что произведение двух симметричных матриц есть симметричная матрица тогда и только тогда, когда матрицы коммутативны, т. е. АВ==ВА. 9. Пусть А — нижняя треугольная ортогональная матрица. Показать, что А — диагональная матрица. Какие у нее диагональные элементы? 10. Показать, что если столбцы матрицы А линейно независимы, то А'А— положительно определенная матрица (т. е. х^^хХ) для всех х-т^О). Для мат- рицы .1000 1020\ A=(l000 1000 ) Ч 000 1000/ вычислить BsAU, используя арифметику с четырьмя десятичными цифрами. Показать, что В не является положительно определенной матрицей, рассмотрев вектор х==(1, —I)7. Объяснить это кажущееся противоречие. 11. Рассмотрим систему уравнений, полученную дискретизацией дифферен- циального уравнения в подразделе З.А.1. Показать, что при дискретизации любо- го дифференциального уравнения вида а (х) у"+Ь (х) у'+с (х) y=f (х), У ^^У0' У (Xl) =У1> соответствующая линейная система — трехдиагональна. 12. Показать, что нормальные уравнения, выведенные в подразделе З.А.2, имеют симметричную матрицу коэффициентов. 13. Дать интуитивные объяснения, почему организационные модели должны быть разреженными для каждой из следующих структур: (а) электрическая сеть для коммунальных услуг; (б) региональная система канализации и очистки воды; (в) министерство почт; (г) сборочный конвейер для автомобилей. 14. Показать, как вычислить вектор z=B-l(2A+I)(C-Ч-A)b без явного вычисления обратных матриц, используя только операции матричной и векторной арифметики и решение систем линейных уравнений. 15. Рассмотрим системы уравнений: (а) ixi=bi, i=l, 2, . . . , М; (б) xi/i=ci, i=l, 2, . . . , М; (в) Xj/10=di, i=l, 2, . . ., М. Выписать коэффициенты матриц и вычислить их определители для М=5, К), 20, 50. Предполагаете ли вы встретить какие-либо затруднения при решении этих систем уравнений? Метод исключения Гаусса и LU- разложение 39 Эти уравнения могут быть решены обратной подстановкой. Отметим, что вычисление МРЬ в точности соответствует прямому ходу. Пример 5.1: Решение вырожденной совместной системы. Форму- лировка теоремы 5.1 позволяет нам сделать попытку решения вы- рожденной системы, такой, как х+ У+ z+ w=4, 2х+3у+ z+ w=7, или Ах=Ь, х+ У+ z+ w=4, x+2y+2z+2w=7. После первого шага исключения Гаусса имеем x+y+z+w=4, множители строки 1, y—z—w=—l, —2, 0=0, -1, y+z+w=3, —1. После следующего шага имеем x+y+z+ w=4, множители строки 2, у—z— w=—1, 0=0, О, 2z+2w=4, —1. Теперь переставим местами третье и четвертое уравнения, чтобы получить x+y+z+ w=4, у—z— w==—1, 2z+2w=4, 0=0. Матрица перестановок, которая выполняет эту операцию, имеет вид (1000 р ^ о 1 о о 3 0001 0010 В этом частном примере мы не будем выполнять действия послед- него шага. Матрица М определяется множителями, использован- ными на последовательных шагах метода исключения Гаусса. Напомним, что основной шаг исключения Гаусса, такой, как «ум- ножить первую строку на т^ и прибавить ее к четвертой строке», выполняется умножением слева на матрицу __ Метод исключения Гаусса и LU-разложение 41 It 0 0 0 1 О О О /I 1 1 1\ /I 1 1 1' -2 100010023111=01-1-1 -1-110000111111 00 2 2 1 001 0010 \1222J\00 О О M P A = U Суммарный Суммарный Исходная Верхняя здлрект эуузект матрица треугольная множителей перестановок матрица Чтобы продолжить решение, мы вычислим произведение МРЬ, которое является вектором-столбцом (4, —1, 4, О)1. Отметим, что на практике мы не вычисляем M, а храним именно множители, по- скольку мы можем вычислить Mb на их основе так же быстро, как и на основе самой матрицы M. Теперь попытаемся выполнить обратную подстановку, чтобы найти решение. Мы никоим образом не можем решить последнее уравнение, однако оно автоматически удовлетворяется при любом значении w. Поэтому можно присвоить w произвольное значение. Продемонстрируем обратную подстановку для трех случаев: Случай w=0: Случай w==l: 2z=4, т. е. z=2, 2z=4—2, т. е. z=l, у=-1+2=1, у=_1+2=1, х=4—3=1; х=4—3=1; Случай w==6.7s 2z=4—13.4, т. е. z=—4.7, у=-1+2=1, х=4—3=1. Это, конечно, иллюстрирует тот факт, что вырожденная система Ах==Ь, которая является совместной (т. е. имеет хотя бы одно решение), имеет бесконечно много решений. Отметим, что если четвертая компонента вектора b не равна 7, то последнее уравнение становится уравнением «нуль = что-то ненулевое», которое невозможно решить, и таким образом вырож- денная система Ах=Ь совсем не будет иметь решения. Несложные рассуждения показывают, что: 1. А вырождена, если U имеет нуль на диагонали. 2. Для вырожденной матрицы А система Ах=Ь совместна, если при обратной подстановке нули на диагонали появляются только в уравнениях вида 0=0. Неизвестные, соответствующие диагональным нулям в этих уравнениях, могут быть определены произвольно. Доказательство теоремы 5.1. Рассмотренный выше простой пример фактически только иллюстрирует механизм доказательства, однако формальное математическое доказательство все же необ- Метод исключения Гаусса и LU- разложение 43 где L^_2 — нижняя треугольная матрица порядка k—2 и Q — мат- рица перестановок порядка п—k+1. Можно проверить, что мат- рица PkMk-iPic — нижняя треугольная (результатом этого про- изведения является воздействие Q на В). Тогда мы имеем PkMk-i = PkM^PkPk == Mk-iPk, где Mk_i — нижняя треугольная. Это дает нам возможность за- писать АК = MkPkM^Pk-iA == MkMk_iPkPk-iA = MkP„A, что и требовалось получить. Здесь P^=PkPk-r Правильность ги- потезы математической индукции для k=l устанавливается из тех же соображений, и доказательство завершено. Доказательство теоремы 5.2. Прежде всего необходимо пока- зать, что отличие от нуля главных миноров означает, что не требует- ся проводить перестановки строк. Доказательство проводится также по индукции относительно k. Согласно гипотезе математи- ческой индукции aii=/=0, так что матрица Pi не нужна. Напомним, что треугольная матрица вырождена тогда и только тогда, когда хотя бы один из ее диагональных элементов равен нулю. Тем самым верхняя треугольная подматрица в A^_i не вы- рождена (по предположению относительно главных миноров) и никаких перестановок строк не требуется. Если а^==0, то главный минор порядка k равен нулю, что является противоречием. Таким образом, a^k^O и доказательство по индукции протекает без пере- становок строк. Построение, используемое при доказательстве теоремы 5.1, дает однозначно определенное разложение, за исключением, воз- можно, различных выборов строк для перестановок. Далее, мат- рица М, а тем самым и обратная к ней матрица L имеют единицы на диагонали, т. е. /ц=1 для всех i. Поскольку никаких переста- новок не нужно производить, они и не требуются в теореме 5.2; полученное разложение единственно и доказательство закончено. 5.Б. ВЫБОР ВЕДУЩЕГО ЭЛЕМЕНТА В МЕТОДЕ ИСКЛЮЧЕНИЯ ГАУССА Процесс перестановки строк в методе исключения Гаусса, в результате которого получаются ненулевые диагональные эле- менты, называется выбором ведущих элементов, а элементы мат- рицы, используемые для исключения, называются ведущими эле- ментами. Отличие от нуля ведущего элемента достаточно для теоретической правильности метода исключения Гаусса, однако для получения надежного результата требуется большая осторож- ность. Это видно из следующего примера. Пример 5.2: Влияние малых ведущих элементов на процесс ___ __ Метод исключения Гаусса и LU-разложение 45 ошибки. Мало что можно сделать, чтобы обойти первую трудность. Можно лишь попытаться идентифицировать такого рода задачи. Выбор ведущего элемента представляет собой технический прием для устранения второй трудности, связанной с эффектом увеличе- ния ошибок округления. В примере 5.2 именно надлежащий выбор ведущего элемента устранил эту трудность. Не так легко ответить на общий вопрос: насколько хорошо выбор ведущего элемента регулирует увеличение ошибок? В одной из наиболее известных статей фон Неймана и Гольдстайна в 1947 г. показано, что можно ожидать увеличения ошибок в 1012 или боль- шее число раз при решении скромной, скажем 20х20, системы методом Гаусса. Казалось бы, безнадежно решать многие линейные системы такого или большего порядка. Однако отдельные исследо- ватели игнорировали этот факт и пытались так или иначе решать большие системы. Оказалось, что можно с достаточной точностью решать системы из 100, 200 и даже 400 уравнений, если позволяли затраты машинного времени. Это открытие привело к попыткам лучше понять метод Гаусса, и к концу 60-х годов были выяснены следующие четыре факта: Факт 1 (доказанный). Полный выбор ведущего элемента пред- ставляет собой надежную стратегию, поскольку ошибки никогда чрезмерно не возрастают. Коэффициент увеличения ошибок не превосходит fn=(n»2»3i^*4i/8»51/1*.. .^п1^"-!»)^2 для любой системы уравнений порядка п. Нижеследующая таблица показывает, что f„ растет довольно медленно с возрастанием п п 5 10 15 20 30 50 80 100 f„ 5.73 18.30 39.09 69.77 155.5 536.17 1915.51 3552.41 и любая ошибка (или ошибки в коэффициентах системы, или ошибки в процессе вычислений), влияющая на решение, умножается самое большее на коэффициент увеличения. Факт 2 (замеченный экспериментально). Вероятность того, что при использовании частичного выбора ведущего элемента возник- нут трудности, связанные с ростом ошибок, весьма мала. В дей- ствительности скорость увеличения ошибок для частичного выбора редко превосходит более чем в 2—4 раза скорость роста ошибок для полного выбора. Наблюдаемая на практике величина коэф- фициента увеличения fn обычно не превосходит 8 и часто близка к 1, особенно для плохо обусловленных задач. Факт 3 (на основе примера). Была построена специальная пхп- матрица, для которой при использовании частичного выбора веду- Метод исключения Гаусса и LU-разложение 47 4. При помощи калькулятора решить систему, используя только четыре пра- вильно округляемые цифры на каждом шаге. Решить систему три раза: два раза с ведущими элементами ац, a^i и один раз с полным выбором ведущего элемента. Сравнить полученные ответы с точным решением (1.000, 0.5000)' 0.002110xi+0.08204x2=0.084l5, 0.3370xi + 12.84ха=6.757. 5. Показать, что матрица А не вырождена, но не может быть представлена произведением LU нижней и верхней треугольных матриц. Объяснить, почему. 123 А= 241. ^—1 0 2/ 6. Рассмотреть теорему 5.1 и объяснить, (а) какова связь теоремы с обычным исключением Гаусса; (б) когда необходимо осуществлять выбор ведущего элемен- та; (в) почему теорема полезна для решения системы Ах=Ь. 7. Дать пример линейной системы Ах=Ь, у которой det (А) много меньше ошибок округления (по сравнению с величиной элементов матрицы А), но для которой нет вычислительных трудностей, имеющих место при решении систем Ах=Ь. 8. Выполнить исключение Гаусса для решения системы 2х+3у— z+ w=5, x+2y+3z— w=—6, 4x+5y—9z+6w=28, x+ у—4z—2w=7 и получить все ее решения. 9. Матрица А называется матрицей с диагональным преобладанием, если 1ац1> 2 I3!)! Р'ля всех i- Пусть А—симметричная матрица с диагональным пре- обладанием. После первого шага исключения Гаусса матрица примет вид (а у). \0 Ai/ (а) Показать, что Af — матрица с диагональным преобладанием. (б) Показать, что для симметричной матрицы с диагональным преобладанием процесс исключения Гаусса не зависит от того, производится ли выбор ведущего элемента или нет. 10. Показать, что (а) для матрицы с диагональным преобладанием процесс исключения Гаусса без выбора ведущего элемента не может потерпеть неудачу; (б) если А — матрица с диагональным преобладанием, то А — невырожденная матрица. Указание: использовать задачу 9. 11. Пусть А—симметричная положительно определенная матрица (т.е. х^Ах^ для всех х=^0). Показать, что (а) А-1—положительно определенная матрица; (б) А может быть представлена в виде LL'1', где L — нижняя треугольная матрица с положительными диагональными элементами; (в) для этой матрицы А процесс исключения Гаусса без выбора ведущего элемента не может потерпеть неудачу. Указание: использовать задачу 9. 12. Можно показать, что коэффициент увеличения ошибок округления мал для некоторых матриц специального вида. Доказать, что для следующих случаев имеют место указанные предельные значения коэффициентов увеличения ошибок: (а) для трехдиагональной матрицы коэффициент увеличения <3; Метод исключения Гаусса и LU- разложение 49 ВЫЧИСЛЕНИЕ МНОЖИТЕЛЕЙ И ЗАПОМИНАНИЕ ИХ В СТОЛБЦЕ К m = ац, == aik/akk ЦИКЛ ПО ЭЛЕМЕНТАМ СТРОКИ I For 3 = k + 1 to n do aij=aii—m*an) ВЫЧИСЛЕНИЕ ПРАВОЙ ЧАСТИ СТРОКИ I bi=bj—m*b^ конец цикла по i конец цикла по k ОБРАТНАЯ ПОДСТАНОВКА ЦИКЛ ПО СТРОКАМ В ОБРАТНОМ ПОРЯДКЕ ПОДСТАНОВКА ИЗВЕСТНЫХ КОМПОНЕНТ РЕШЕНИЯ В ПРАВУЮ ЧАСТЬ n г=Ь;— 2 a^jX, ПРОВЕРКА НА НУЛЬ ВЕДУЩЕГО ЭЛЕМЕНТА If (|ац|> BUFFER) then Xi=r/an ПРОВЕРКА НА НУЛЬ ПРАВОЙ ЧАСТИ ПРИСВОИТЬ НУЛЬ КОМПОНЕНТЕ РЕШЕНИЯ И ПЕ- ЧАТЬ СООБЩЕНИЯ else if (| г К BUFFER) then х;=0 ИЛИ ПЕЧАТЬ СООБЩЕНИЯ ОБ ОШИБКЕ else печать «несовместная система, А—вырождена» конец цикла по i Укажем, как матрицы М и U хранятся на месте матрицы А. Каждый множитель помещается на месте обнуляемого элемента А. Поскольку диагональные элементы М равны единице, то они не хранятся явно. Хотя запоминать матрицы М и U не нужно, когда решается лишь система Ах==Ь, существуют другие ситуации, когда весьма полезно сохранять обе матрицы М и U. 5.В.2. Масштабирование и анализ машинного нуля Алгоритм Гаусса требует анализа ситуаций, когда некоторые получаемые числа равны нулю. С математической точки зрения было бы правильным положить BUFFER =0, однако на практике этого сделать нельзя. Числа, которые при использовании точной арифметики должны быть нулями, почти наверняка не равны нулю при использовании арифметики с плавающей точкой из-за ошибок округления. Чтобы определить, является ли число нулем, необ- ________________Метод исключения Гаусса и LU-разложение 51 aij выполняется самое большее 2п2 арифметических операций, поэтому c=n2 — самое большее значение, которое следует рас- сматривать. Однако можно ожидать, что будет иметь место значи- тельное взаимное уничтожение ошибок округления и потому c=n — достаточно надежный выбор. Большинство численных расчетов проводится с относительно высокой точностью (от 10 до 15 десятич- ных цифр), и в этих случаях можно выбирать значения с доста- точно произвольно. Подведем итоги, отметив следующие три факта: Факт 1. Вычисления более устойчивы (менее чувствительны к изменениям), если все элементы матрицы приблизительно одина- ковы. Факт 2. Не известны способы масштабирования матриц в общем случае, и существуют матрицы, которые не могут быть удовле- творительно масштабированы. Факт 3. На практике обычно масштабируют строки (делением каждого уравнения на его наибольший коэффициент). Это доста- точно надежный способ для обычно встречающихся задач, и есть основания полагать, что апостериорные проверки, которые будут рассмотрены в гл. 8, дают предупреждающие диагностики, если возникают затруднения при вычислениях. 5.В.З. Модификации алгоритма Гаусса Существуют много модификаций основного алгоритма исклю- чения Гаусса. Выбор конкретной модификации представляет собой одну из задач проектирования математического обеспечения. Проблема выбора будет рассмотрена подробнее в гл. 6. Здесь мы перечислим некоторые более очевидные модификации: 1. Элементы матрицы А сохраняются; должны быть выделены дополнительные участки памяти в качестве рабочих вместо ис- пользования для этих целей частей матрицы А. 2. При выборе ведущих элементов строки матрицы А явно не переставляются; порядок исключения неизвестных определяется совокупностью указателей. 3. Матрица А факторизуется к виду LU; множители L и U ис- пользуются на следующем этапе для решения системы. 4. Модификация, позволяющая решать систему со многими правыми частями, т. е. матричное уравнение АХ==В. 5. Вычисление А"1 посредством решения матричного уравнения АХ=1. Здесь правые части суть последовательно столбцы единич- ной матрицы, а решения уравнений образуют столбцы матрицы А"1. В задаче 19 этой главы показано, что вычисление А"1 может Метод исключения Гаусса и LU- разложение 53 тогда, когда неудобно хранить промежуточные результаты, на- пример, в случае использования калькуляторов. Суммирование во внутреннем цикле можно выполнять без отдельной записи про- межуточных частичных сумм. Это позволяет существенно эконо- мить время и, кроме того, увеличить надежность вычислений, по- скольку запись чисел от руки приводит к появлению ошибок. Второе преимущество проявляется для вычислительных машин, обладающих сравнительно быстрой арифметикой с расширенной точностью или особенно быстрой операцией скалярного произ- ведения. Пусть вычисление сумм k -I k -1 2j BijSik' JL Ski^j с удвоенной точностью относительно дешево и требует только на 5 процентов больше времени, чем вычисление с обычной точностью, Было бы правдоподобным предположить, что использование ариф- метики с удвоенной точностью позволяет исключить влияние большинства ошибок округления в вычислениях, и детальный анализ, который здесь опускается, доказывает, что это имеет место. Раньше существовали вычислительные машины, которые обладали быстрой арифметикой с расширенной точностью для вычисления такого рода сумм, и алгоритм Краута был для них весьма привле- кателен. Весьма вероятно, что это свойство снова станет широко распространенным. Некоторые современные вычислительные ма- шины допускают особенно эффективное вычисление вышеупомяну- тых сумм в режиме арифметики с обычной точностью, и алгоритм Краута выгоден для таких машин. Наконец, отметим, что существует специальная модификация исключения Гаусса для положительно определенных симметричных матриц. Если А — симметричная матрица, то представляется ве- роятным, что при решении системы Ах=Ь можно избежать поло- вины трудозатрат по сравнению с решением системы с матрицей общего вида. Такой экономии достигает разложение Холесского L-lA=LT, для которого вычисление матрицы L может быть выпол- нено за половину операций, требуемых для LLJ-разложения. Существо метода видно из доказательства возможности раз- ложения методом математической индукции по размеру матрицы. Рассмотрим сначала случай 2х2-матрицы /ац а„\ ,/„ 0\ f\ —— I I • *—' ~~~ I • • ) » \а„ а^У V2i W Потребуем, чтобы A=LL1, или, что то же самое, ац==ч1> ^^'ii^u ^u^'n'm а^а == t-zi "т "ч.ч' Метод исключения Гаусса и LV-разложение 55 For k = 1 to n do ВЫЧИСЛЕНИЕ ОЧЕРЕДНОГО ВЕКТОРА-СТРОКИ, КРОМЕ ДИАГОНАЛЬНОГО ЭЛЕМЕНТА For i=l to k—1 do / i-^l \ , aki=l aki—.S aijakji/ai, конец цикла по i ВЫЧИСЛЕНИЕ ДИАГОНАЛЬНОГО ЭЛЕМЕНТА У—k^г- ^= |/ Skk—.S Sk) конец цикла по k ТЕПЕРЬ МНОЖИТЕЛЬ L РАСПОЛОЖЕН В НИЖНЕМ ТРЕУГОЛЬНИКЕ МАТРИЦЫ А Алгоритм 5.3 может быть сделан более устойчивым заменой оператора, содержащего квадратный корень, следующим опера- тором: k -1 t=3kk— 2 а^', i=i If t>0 then akk=^F, else а^=0. Это позволит избежать останова программы, когда элемент близок к нулю или точно равен нулю и ошибки округления случайно делают t отрицательным. 5.В.4. Подсчет числа операций Трудозатраты при решении системы Ах=Ь с помощью исклю- чения Гаусса могут быть оценены непосредственно путем подсчета арифметических операций. Это делается следующим образом: матрица А на k-м шаге имеет вид как на рис. 5.2. Г\————1 t- t 0\______ 1 [n > о . Рис. 5.2. I—.______________I ' • \' Метод исключения Гаусса и LV-разложение 57 дозатраты на вычисление чего-либо (в нашем случае на решение системы Ах==Ь) в терминах основных машинных операций (в нашем случае арифметических операций). Очевидный и важный вопрос: является ли исключение Гаусса наиболее быстрым способом реше- ния Ах=Ь? Для достаточно больших значений п применение алго- ритма Штрассена приводит к методу, имеющему меньшее число операций. Он требует порядка п2'7 сложений и умножений. В на- стоящее время проводится теоретическое изучение и эксперимен- тирование по вопросу о том, насколько выгоден этот новый алго- ритм на практике. Предварительные результаты показывают, что он становится более эффективным для п, больших 250. Задачи разд. 5.В 1. Реализовать алгоритм 5.1 в виде подпрограммы на языке Фортран. Дать обоснования выбора последовательности параметров обращения. Объяснить вы- бор значения машинно-зависимой константы BUFFER. 2. Реализовать алгоритм 5.1 в виде трех подпрограмм на языке Фортран. В подпрограмме FACTOR LU-разложение формируется на месте матрицы А без выбора ведущего элемента. Подпрограмма FSUB использует результаты работы подпрограммы FACTOR для соответствующих преобразований правой части си- стемы. Подпрограмма BSUB выполняет обратную подстановку. Дать полное описание последовательностей параметров обращений к подпрограммам. Составить набросок программы на Фортране, использующей эти три подпрограммы для решения следующих задач: — решить две системы с матрицами А и В размеров 10Х 10 и 20Х20; — решить семь систем Axi=bi с векторами правых частей, расположенными в массиве В (10, 7); — решить 16 систем Dyi=ei с векторами правых частей, расположенными в массиве Е (20, 16). 3. Рассмотреть две библиотечные подпрограммы FACTOR (А) и SOLVE (A, X, В), в которых для простоты различные аргументы, специфицирующие размеры А или системы, опущены. Подпрограмма FACTOR вычисляет треугольное разло- жение на месте матрицы А. Подпрограмма SOLVE решает систему Ах=Ь в пред- положении, что матрица А представлена своим треугольным разложением. Пусть требуется вычислить вектор у по формуле у=В (SC+D-^x-HI—D-'Dz, Составить набросок программы на Фортране, использующей FACTOR и SOLVE для вычисления у без явного обращения матриц. Для краткости можно использо- вать в операциях векторной и матричной арифметик операторы типа «прибавить вектор V к вектору S» или «умножить вектор V на матрицу М». Пояснить исполь- зование переменных и способы хранения информации. 4. На основе алгоритма 5.1 .составить эффективный алгоритм приведения трехдиагональной матрицы к верхней треугольной форме без выбора ведущего элемента. 5. Составить набросок подпрограммы на Фортране, реализующей алгоритм задачи 4 в модификации, рассчитанной на экономию машинной памяти. Описать в деталях структуры данных и последовательность параметров обращений. 6. Обобщить алгоритм и программы задач 4 и 5 на случай ленточных матриц с шириной ленты К (без выбора ведущего элемента). Метод исключения Гаусса и LU- разложение 59 могут быть определены решением систем Axi=ej для i=l, 2 , . . . , N. Пусть с помощью первой части алгоритма 5.1 матрица А уже представлена своим треу- гольным разложением. На основе той части алгоритма 5.1, которая осуществляет обратную подстановку, составить эффективный алгоритм вычисления х,. При составлении алгоритма учесть специальный вид правых частей е\. 20. Подсчитать число операций для алгоритма задачи 19 и оценить экономию трудозатрат по сравнению с решением систем Axi=bi для i==l, 2, . , , , N, где Ь, — произвольные векторы правых частей. ________________Метод исключения Гаусса и LU-разложение 61 aij выполняется самое большее 2п2 арифметических операций, поэтому c=n2 — самое большее значение, которое следует рас- сматривать. Однако можно ожидать, что будет иметь место значи- тельное взаимное уничтожение ошибок округления и потому c=n — достаточно надежный выбор. Большинство численных расчетов проводится с относительно высокой точностью (от 10 до 15 десятич- ных цифр), и в этих случаях можно выбирать значения с доста- точно произвольно. Подведем итоги, отметив следующие три факта: Факт 1. Вычисления более устойчивы (менее чувствительны к изменениям), если все элементы матрицы приблизительно одина- ковы. Факт 2. Не известны способы масштабирования матриц в общем случае, и существуют матрицы, которые не могут быть удовле- творительно масштабированы. Факт 3. На практике обычно масштабируют строки (делением каждого уравнения на его наибольший коэффициент). Это доста- точно надежный способ для обычно встречающихся задач, и есть основания полагать, что апостериорные проверки, которые будут рассмотрены в гл. 8, дают предупреждающие диагностики, если возникают затруднения при вычислениях. 5.В.З. Модификации алгоритма Гаусса Существуют много модификаций основного алгоритма исклю- чения Гаусса. Выбор конкретной модификации представляет собой одну из задач проектирования математического обеспечения. Проблема выбора будет рассмотрена подробнее в гл. 6. Здесь мы перечислим некоторые более очевидные модификации: 1. Элементы матрицы А сохраняются; должны быть выделены дополнительные участки памяти в качестве рабочих вместо ис- пользования для этих целей частей матрицы А. 2. При выборе ведущих элементов строки матрицы А явно не переставляются; порядок исключения неизвестных определяется совокупностью указателей. 3. Матрица А факторизуется к виду LU; множители L и U ис- пользуются на следующем этапе для решения системы. 4. Модификация, позволяющая решать систему со многими правыми частями, т. е. матричное уравнение АХ==В. 5. Вычисление А"1 посредством решения матричного уравнения АХ=1. Здесь правые части суть последовательно столбцы единич- ной матрицы, а решения уравнений образуют столбцы матрицы А"1. В задаче 19 этой главы показано, что вычисление А"1 может