WWW.DISS.SELUK.RU

БЕСПЛАТНАЯ ЭЛЕКТРОННАЯ БИБЛИОТЕКА
(Авторефераты, диссертации, методички, учебные программы, монографии)

 

Pages:   || 2 |

«1 октября 2010 Н.Ю.Золотых Nikolai.Zolotykh Министерство образования и науки Российской Федерации Федеральное агентство по образованию Нижегородский государственный университет ...»

-- [ Страница 1 ] --

Некоторые страницы книги

С. А. Белов, Н. Ю. Золотых

Численные методы линейной алгебры

Если вы хотите приобрести всю книгу, напишите мне письмо.

1 октября 2010

Н.Ю.Золотых

Nikolai.Zolotykh@gmail.com

Министерство образования и науки Российской Федерации

Федеральное агентство по образованию Нижегородский государственный университет им. Н.И. Лобачевского С.А. Белов, Н.Ю. Золотых

ЧИСЛЕННЫЕ МЕТОДЫ

ЛИНЕЙНОЙ АЛГЕБРЫ

Лабораторный практикум Нижний Новгород Издательство Нижегородского госуниверситета 2005 УДК 519. ББК 22.193. Б Белов С.А., Золотых Н.Ю. Численные методы линейное алгебры. Лабораторный практикум. Нижний Новгород: Изд-во Нижегородского госуниверситета им. Н.И. Лобачевского, 2005.

264 с.

ISBN 5–85746–837–X Учебное пособие посвящено численным методам линейной алгебры и затрагивает следующие темы: матричные нормы и число обусловленности, прямые и итерационные методы решения линейных систем, задачу наименьших квадратов, технологию разреженных матриц. По каждой теме даются необходимый теоретический материал и задания к лабораторным работам. Помимо этого, книга содержит руководство по учебнопрактической библиотеке численных методов NL, доступный путеводитель по пакету LAPACK и его версии на языке C CLAPACK и краткое руководство по использованию системы Matlab.

Подготовлен в учебно-исследовательской лаборатории Математические и программные технологии для современных компьютерных систем (Информационные технологии) факультета ВМК ННГУ при поддержке Фонда содействия развитию малых форм предприятий в научно-технической сфере.

ISBN 5–85746–837–X ББК 22.193. c С.А. Белов, Н.Ю. Золотых, Оглавление Предисловие Обозначения и соглашения Глава 1. Машинная арифметика и анализ ошибок 1.1. Нормы векторов и матриц................. 1.2. Число обусловленности матрицы............. 1.2.1. Оценки для числа обусловленности........ 1.2.2. Решение возмущенной системы......... 1.2.3. Апостериорные оценки............... 1.2.4. Покомпонентная относительная обратная ошибка 1.3. Машинная арифметика................... 1.3.1. Числа с плавающей точкой............ 1.3.2. IEEE-арифметика.................. 1.4. Прямой и обратный анализ ошибок.





........... 1.4.1. Прямой анализ ошибок............... 1.4.2. Анализ чувствительности............. 1.4.3. Обратный анализ ошибок............. 1.5. Производительность алгоритмов............. 1.6. Библиотека NL........................ Глава 2. Прямые методы решения линейных систем 2.1. Введение........................... 2.2. Метод Гаусса......................... 2.3. LU -разложение....................... 2.3.1. Исключение по столбцу.............. 2.3.2. Исключение по строке............... 2.3.3. Компактная схема................. 2.3.4. Замечания о производительности......... 2.4. Схемы выбора ведущего элемента............. 2.4.1. Стратегия частичного выбора........... 4 Оглавление 2.4.2. Стратегия полного выбора............. 2.5. Анализ ошибок....................... Оглавление 6.1. Краткое введение в Matlab................ 6.2. Основные функции для работы с матрицами...... 6.3. Cистемы линейных уравнений............... 6.3.1. Решатели систем линейных уравнений...... 6.3.5. Сингулярное разложение.............. 6.3.7. Число обусловленности............... 6.3.8. Оценка числа обусловленности.......... 6.3.9. Обратная и псевдообратная матрицы...... 6.4. Разреженные матрицы................... 6.4.1. Создание и преобразование разреженных матриц 6.4.2. Алгоритмы переупорядочения........... 6.4.3. Неполная факторизация.............. 6.4.4. Решение систем линейных уравнений...... 8.1. Рекомендации по установке библиотеки......... 8.2. Особенности вызова функций............... Учебное пособие, предлагаемое вниманию читателя, посвящено численным методам линейной алгебры. В нем затрагиваются следующие темы: матричные нормы и число обусловленности, прямые и итерационные методы решения линейных систем, задача наименьших квадратов, технология разреженных матриц. По каждой теме мы даем необходимый теоретический материал и задания к лабораторным работам. Помимо этого книга содержит руководство по учебно-практической библиотеке численных методов NL, доступный путеводитель по пакету LAPACK и его версии на языке C CLAPACK и краткое руководство по использованию системы Matlab. Мы полагаем, что при выполнении задания студент воспользуется указанным программным обеспечением или, возможно, напишет свою программу с нуля.

На сегодняшний день объем информации в области методов вычислений практически необозрим. Существуют большие библиотеки программ, такие, как библиотека Netlib [32] (включающая пакеты LINPACK, EISPACK, LAPACK, BLAS и др.), библиотеки численного анализа Numerical Algorithms Group NAG [31], научная библиотека GSL ассоциации GNU [28], библиотека НИВЦ МГУ [24] и др. Указанные библиотеки широко применяются современными вычислителями, однако в учебных курсах и практикумах используются достаточно редко. К причинам последнего (среди прочих) следует отнести некоторую их громоздкость. Достаточно сложны соглашения о параметрах функций. Исходный код, как правило, написан на Фортране. Эти свойства часто затрудняют возможность использования исходного кода этих пакетов в рамках студенческих лабораторных работ. Поэтому мы не справились с искушением написать свою собственную библиотеку. Предлагаемая библиотека NL написана А.М. Тагуновым и Н.Ю. Золотых и протестирована стажерами лаборатории Информационные технологии А.Н. Половинкиным, К.В. Корняковым и В.В. Рябовым.





Она, конечно, уступает по производительности таким современным пакетам, как LAPACK, однако ее исходный код намного более прозрачен и студент всегда может обратиться к нему для уяснения деталей алгоритма или с целью внесения изменений.

Для работы с пособием необходимы знания линейной алгебры и математического анализа в объеме соответствующих университетских курсов для студентов естественно-научных специальностей. Предполагается, что лабораторные работы выполняются параллельно с изучением теоретического материала по численным методам. В качестве учебников по современным численным методам линейной алгебры мы рекомендуем [2, 4, 5], в качестве справочника [3]. При написании пособия мы сами активно пользовались этими книгами. Описание функций и подпрограмм в главах 6, 7 следует руководствам [18, 23]. По духу наше пособие близко учебникам [1, 10, 17, 21, 22].

В книге приводятся доказательства только некоторых ключевых утверждений, что вполне, как мы считаем, оправдывается поставленными задачами. Опущенные доказательства во многих случаях могут быть восстановлены: это, по мнению авторов, являлось бы хорошим упражнением. В тех случаях, когда известные нам доказательства сложны, мы всегда приводим ссылку на соответствующую литературу.

Мы благодарим за сотрудничество директора по развитию Нижегородской лаборатории корпорации Интел Л.В. Нестеренко. Во время работы над пособием авторы получали финансовую поддержку Фонда содействия развитию малых форм предприятий в научно-технической сфере. Мы искренне приПредисловие знательны его генеральному директору И.М. Бортнику и начальнику отдела развития инфраструктуры В.Н. Оганесяну.

Мы надеемся, что настоящее пособие будет полезно студентам и всем тем, кто на практике применяет численные методы линейной алгебры.

Обозначения и соглашения В пособии рассматриваются задачи линейной алгебры над полем вещественных чисел, поэтому все скаляры, векторы и матрицы считаются вещественными (не комплексными). Заметим, что некоторые приведенные результаты и описанные алгоритмы без труда переносятся и на комплексный случай, тогда как для некоторых алгоритмов и фактов такое перенесение невозможно.

Под арифметическим вещественным пространством понимается линейное пространство вещественных столбцов заданной высоты со стандартно определенными операциями сложения и умножения вектора на вещественный скаляр.

Как правило, векторы мы обозначаем малыми латинскими буквами (a, b, c,... ), матрицы большими латинскими буквами (A, B, C... ), а их элементы соответствующими малыми латинскими буквами с индексами (aj, bj, cj,..., aij, bij, cij,... ).

Если матрица задана выражением, например A + B, то для ссылки на ее элемент i-й строки j-го столбца будем писать (A + B)ij. Аналогичные соглашения примем и для векторов.

Под неравенством A B понимается система покомпонентных неравенств aij bij для всех i и j.

rank A det A Обозначения и соглашения diag(D1,..., Dk ) cond A, condp A condLS A gpp, gcp коэффициенты роста элементов в алгоритме LU -разложения с частичным и соответственно полным выбором ведущего Машинная арифметика Вещественное или комплексное линейное пространство называется нормированным, если каждому его вектору x поставлено в соответствие число, называемое нормой и обозначаемое x, так что выполнены следующие аксиомы нормированного пространства:

1) x 0, причем x = 0, только если x = 0;

для любых векторов x, y и любого числа. Примером нормы в арифметическом пространстве может служить функция где p 1. При p = 1, 2 и при p получаем соответственно манхеттенскую x 1, евклидову x 2 и чебышеву x нормы:

1.1. Нормы векторов и матриц Обобщением евклидовой нормы является функция где A симметричная положительно определенная матрица.

Обычно от нормы, введенной в пространстве квадратных матриц, требуют, чтобы помимо обычных аксиом 1)–3) она удовлетворяла свойству 4) A · B A · B для любых матриц A, B.

Такую норму будем называть матричной, или кольцевой. Далее рассматриваются только такие нормы.

Пусть в пространстве n-мерных арифметических векторов введена норма x, а в пространстве квадратных матриц порядка n введена норма A. Эти нормы называются согласованными, если для любого вектора x и любой матрицы A выполнено неравенство Матричную норму A, определенную соотношением называют подчиненной векторной норме x. Подчиненная норма всегда согласована со своей векторной нормой.

Матричная норма, подчиненная норме x p, обозначается A p. Нормы A 1, A 2, A называются манхеттенской, спектральной и чебышевой соответственно. Можно доказать, что 14 Глава 1. Машинная арифметика и анализ ошибок Для симметричной матрицы A имеем Матричная норма согласованы.

1.2. Число обусловленности матрицы Числом обусловленности невырожденной матрицы A называется величина Из определения следует, что если матричная норма подчинена векторной, то Если матрица A вырождена, то полагают cond A =.

Пусть в (1) в качестве матричной нормы A используется A p. Тогда число обусловленности обозначим condp A. Чаще других используется матхеттенское, спектральное, чебышево и фробениусово числа обусловленности: cond1 A, cond2 A, cond A, condF A соответственно.

Для спектрального числа обусловленности справедливо 1.2. Число обусловленности матрицы откуда для симметричной матрицы A получаем 1.2.1. Оценки для числа обусловленности Приведем несколько результатов, дающих возможность получать оценки числа обусловленности. Алгоритмы, позволяющие экспериментально оценить число обусловленности, будут даны в разделе 2.6. Пусть A невырожденная матрица, а A+B вырожденная. Тогда (см. задачу 15 из раздела 9.1 на стр.

239) Пример 1.1. Рассмотрим матрицу где 0 1. Пусть Для числа обусловленности произведения двух матриц справедливо двустороннее неравенство:

16 Глава 1. Машинная арифметика и анализ ошибок Для спектрального числа обусловленности имеем:

Как правило, оценка (5) слишком занижена.

Пример 1.2. Пусть Прямые вычисления дают:

Рассмотрим матрицу Оценка (5) дает cond2 B 49.71. Заметим, что Воспользовавшись (4), получаем, например, для спектрального числа обусловленности:

Заметим, что прямые вычисления числа обусловленности приводят нас к значению cond2 B = 2.20 · 103.

1.2. Число обусловленности матрицы Пусть вместо системы Ax = b решается некоторая близкая возмущенная система (A + A)x = b + b. Пусть x точное решение исходной системы, x + x точное решение возмущенной системы. Величина x называется абсолютной1, а Можно показать [2, 5] (см. также задачу 16 из раздела 9.1 на стр. 240), что если матричная и векторная нормы согласованы, В частности, Справедливо также неравенство Таким образом, если A / A мал, то относительная ошибо ка x / x при решении системы линейных уравнений не превосходит некоторой величины, пропорциональной относительной ошибке в коэффициентах системы, с коэффициентом пропорциональности cond A.

Чем больше число обусловленности матрицы, тем более чувствительно решение системы по отношению к возмущению ее Величину x также часто называют абсолютной ошибкой.

Иногда вместо x / x рассматривают x / x + x.

18 Глава 1. Машинная арифметика и анализ ошибок коэффициентов. Матрицы с большим числом обусловленности называются плохо обусловленными. Если матрица системы плохо обусловлена, то сама система также называется плохо обусловленной.

Пусть x точное, а x+x вычисленное решения системы Ax = b. Невязкой вектора x + x называется вектор r = b A(x + x).

При решении реальных задач точное решение x, как правило, неизвестно. Однако легко может быть вычислена невязка r и найдены апостериорные оценки точности найденного решения:

Вторая из приведенных оценок (как правило, более точная) часто обозначается FERR (forward error).

1.2.4. Покомпонентная относительная обратная ошибка В разделе 1.2.2 мы дали оценку относительной ошибки решения системы Ax = b при условии, что найденное решение x + x удовлетворяет возмущенной системе (A + A)(x + x) = b + b. Эта оценка пропорциональна числу обусловленности cond A и относительному возмущению.

Иногда эта оценка весьма завышена. Лучшие оценки дает альтернативная теория.

Пусть x точное, а x + x вычисленное решения системы Ax = b. Наименьшее число 0, для которого существуют возмущения A и b, удовлетворяющие оценкам |A| |A|, |b| |b| и уравнению (A + A)(x + x) = b + b, называется покомпонентной относительной обратной ошибкой (BERR, 1.3. Машинная арифметика backward error). Величина может быть вычислена по найденному приближенному решению x + x и вычисленной для него невязке r на основании формулы [5]:

Можно показать [5], что Величина |A1 | · |A| называется относительным покомпонентным числом обусловленности.

Для представления (аппроксимации) действительных чисел в компьютере наиболее часто используются числа с плавающей точкой. Число с плавающей точкой это запись числа в виде где (d0.d1 d2... dt ) представление числа в -ичной системе счисления, 0 di, e u, причем называется базой, di разрядами, (d0.d1 d2... dt ) мантиссой, e показателем.

Величины, t,, u зависят от компьютера3. Компьютеры, для которых = 2, называются двоичными.

Заметим, что здесь представлена несколько упрощенная математическая модель машинной арифметики с плавающей точкой (ср. с упомянутым ниже стандартом IEEE 754). Однако этой модели вполне достаточно для анализа численных методов.

20 Глава 1. Машинная арифметика и анализ ошибок Если d0 = 0, то число называется нормализованным, если d0 = 0, то число называется ненормализованным. Некоторые компьютеры работают только с нормализованными числами4, другие как с нормализованными, так и с ненормализованными. Преимущество нормализованных чисел единственность представления.

Обозначим через M максимальное положительное, а через m минимальное положительное число, представимое в компьютере. Если используются только нормализованные числа, то Пусть x некоторое вещественное число. Обозначим через (x) его представление в компьютере. Если x можно точно представить в виде (8), то следует ожидать, что (x) = x. Если x нельзя точно представить в виде (8), то либо |x| M, либо |x| M. В первом случае говорят, что произошло так называемое переполнение (overow): некоторые компьютеры генерируют ошибку, некоторые заменяют x числами Inf или Inf, представляющими + или соответственно. Во втором случае для представления x используется аппроксимация (x) x. Существует несколько способов аппроксимации. Один из них отбрасывание разрядов мантиссы числа x начиная с (t+1)-го (арифметика с отбрасыванием). Другой способ округление x до ближайшего числа, представимого в компьютере (арифметика с округлением)5.

Наряду с переполнением выделяют также случай появления машинного нуля, или потери значимости (underow). Это происходит, когда |x| m. Обычно это не приводит к катастрофическим проблемам, которые могут возникнуть при переполЧисло 0, не допускающее нормализованной записи, представляется отдельно.

Неоднозначности разрешаются дополнительными правилами. В двоичных компьютерах (например, поддерживающих описанный ниже стандарт IEEE 754) часто используют правило округления до ближайшего четного, т. е. до числа с dt = 0.

1.3. Машинная арифметика нении. Для большинства современных компьютеров в этом случае (x) = 0.

Назовем машинным эпсилон величину Можно сделать вывод, что в любой из двух упомянутых машинных арифметик, если нет переполнения и потери значимости, для относительной ошибки ((x) x)/x справедливо неравенство При анализе ошибок округления разбираемых ниже численных алгоритмов мы предполагаем это условие выполненным.

При выполнении арифметических операций над числами с плавающей точкой результат, как правило, также не может быть представлен в компьютере точно. Со сложением и вычитанием это происходит реже, с умножением и делением чаще. Пусть обозначает одну из операций +,,, /. Обозначим через x y результат выполнения в компьютере операции над числами x и y. Операция называется машинным сложением, машинным вычитанием, машинным умножением и машинным делением соответственно. Наиболее разумным способом выполнения в компьютере операции является способ, при котором если x, y заданы точно и в результате выполнения операции не произошло переполнения или потери значимости, то x y = (x y). При анализе ошибок округления разбираемых ниже численных алгоритмов мы предполагаем, что Иногда для простоты машинным эпсилон называют величину t для обеих арифметик.

22 Глава 1. Машинная арифметика и анализ ошибок это свойство выполнено:

Аналогичное свойство будем предполагать у операции извлечения квадратного корня.

Еще раз отметим, что условие (10) при реальных вычислениях на компьютере выполняется только при условии, что в результате выполнения операции не произошло переполнения или потери значимости, т. е. если x y = 0 и x y = ±Inf.

На большинстве современных компьютеров вместе с обычными числами, а также числами Inf, Inf используется так называемое нечисло NaN. Оно возникает как результат некорректных операций 0/0, Inf/Inf, 0 · Inf и т. п.

В 1985 г. на арифметику двоичных компьютеров Ассоциацией IEEE был утвержден стандарт № 754 [20]. Машины, сконструированные согласно этому стандарту, должны использовать арифметику с округлением. Машинные операции сложения, вычитания, умножения, деления должны удовлетворять (10). Аналогичное свойство должно быть выполнено для операции извлечения квадратного корня. Для представления вещественных чисел было разработано несколько форматов. Основные их них: числа двойной точности и числа одинарной точности.

Под числа одинарной точности (single) выделяется 32 бита:

1 бит под знак числа, 8 бит под показатель, 23 бита под мантиссу (см. таблицу 1.1). Стандарт предусматривает как нормализованные, так и ненормализованные числа. У нормализованных чисел разряд d0 (который всегда равен 1) не хранится явно. Ввиду использования таких ненормализованных чисел никогда при сложении и вычитании двух различных нормализованных чисел потеря значимости не наступит (однако она Положительные 0 11...10 11... нормализованные

ненормализованные

ненормализованные

нормализованные

Положительные 0 11...10 11... нормализованные

ненормализованные

ненормализованные

нормализованные

1.m 2e1023 1 11...10 11... 1.4. Прямой и обратный анализ ошибок может возникнуть при умножении или делении). Для нормализованных чисел M = 224 6 108. Кроме нормализованных и ненормализованных чисел есть также 0, +Inf, Inf, NaN.

Под числа двойной точности (double) выделяется 64 бита:

1 бит под знак числа, 11 бит под показатель, 52 бита под мантиссу (см. таблицу 1.2). Стандарт предусматривает как нормализованные, так и ненормализованные числа. У нормализованных чисел разряд d0 = 1 не хранится явно. Для нормализованных чисел M = 253 1 1016. Кроме нормализованных и ненормализованных чисел есть также 0, +Inf, Inf, NaN.

1.4. Прямой и обратный анализ ошибок При реализации алгоритмов численных методов на компьютере неизбежно появляются ошибки округления, поэтому возникает задача оценки точности вычисленного решения.

Обозначим через A входные данные задачи, через X = (A) результат их обработки по точному алгоритму. Будем считать A, X элементами некоторых нормированных пространств.

При реализации на компьютере алгоритм в силу отличий машинной арифметики от точной будет заменен другим алгоритмом M. Этот алгоритм по начальным данным A получает результат XM = M (A). Величины XM X, XM X / X называются абсолютной и относительной ошибкой соответственно. Иногда вместо XM X / X рассматривают XM X / XM. Исследование этих величин называется прямым анализом ошибок.

Итак, прямой анализ ошибок должен дать ответ на вопрос, насколько точное решение близко к вычисленному на компьютере. Относительно машинной арифметики делается, как правило, единственное предположение о выполнении ограничений (9), (10). Заметим, что попытки использования этих неравенств 26 Глава 1. Машинная арифметика и анализ ошибок в лоб для всех выполняемых в алгоритме операций, как правило, приводит к чрезмерно завышенным оценкам ошибок.

В настоящее время распространен более тонкий подход к прямому анализу ошибок. Анализ разбивается на два этапа:

1) анализ чувствительности задачи к возмущениям в данных;

2) обратный анализ ошибок.

Рассмотрим каждый из этих этапов.

Предположим, что вместо задачи с входными данными A решается близкая задача с данными A + A. Точным решением этой возмущенной задачи будет X + X = (A + A).

Анализ чувствительности заключается в исследовании погрешностей X и X / X.

Как правило, удается получить оценки вида где abs, rel некоторые константы, не зависящие от A, называемые абсолютным и относительным числами обусловленности соответственно. Эти величины показывают, как сильно ошибка во входных данных может влиять на выходной результат. Для краткости относительное число обусловленности называют просто числом обусловленности.

Анализ чувствительности задачи иногда называют теорией возмущений.

Для задачи решения системы линейных уравнений анализ чувствительности был проведен в разделах 1.2.2, 1.2.4. Далее Неравенства (11) могут быть справедливы лишь для достаточно малых значений A и A / A или быть асимптотическими с точностью до слагаемых o( A ) или o ( A / A ).

1.4. Прямой и обратный анализ ошибок в пособии будет проведен анализ чувствительности других задач. Например, чувствительность задачи наименьших квадратов рассматривается в разделе 3.2.

Перейдем ко второму этапу анализа ошибок. Пусть X = (A), XM = M (A). Предположим, что результат работы XM приближенного алгоритма можно рассматривать как результат работы точного алгоритма с возмущенными входными данными A + A, т. е. XM = (A + A). Если это возможно и A мал, то алгоритм называют обратно устойчио вым, при этом A называется эквивалентным возмущением, или обратной ошибкой. Исследование эквивалентного возмущения называется обратным анализом ошибок. В ходе этого исследования необходимо оценить A, A / A функциями, зависящими от M (и, возможно, зависящими от некоторых параметров задачи).

В разделе 2.5 будет подробно разобран обратный анализ ошибок алгоритма LU -разложения для решения системы линейных уравнений.

Подставляя полученные при обратном анализе оценки для A, A / A в (11), можно получить границы для абсолютной и относительной ошибок. Заметим, что для того, чтобы полученные оценки можно было подставить в (11), необходимо, чтобы A, A / A удовлетворяли требованиям, при которых неравенства (11) выполняются (см. предыдущее примечание). Отсюда и вытекает необходимость условия в определении обратно устойчивого алгоритма, чтобы A было небольшим.

Итак, в рассмотренном подходе обратный анализ ошибок возникает как подэтап в прямом анализе.

28 Глава 1. Машинная арифметика и анализ ошибок 1.5. Производительность алгоритмов Кроме точности получаемого решения важной характеристикой алгоритма является его производительность, или трудоемкость.

В настоящее время распространен следующий подход к анализу производительности алгоритмов, реализующих численные методы. Производительность алгоритма измеряют количеством флопов (op oating point operation), т. е. выполняемых алгоритмом операций над числами с плавающей запятой: машинных сложений, вычитаний, умножений, делений и операций извлечения квадратного корня.

Заметим, что такой подход не учитывает времени на выполнение операций с целыми числами (например, с индексами векторов и матриц), времени на пересылку данных между разными устройствами компьютера и т. д. Действительно, часто можно считать, что это время пренебрежимо мало по сравнению со временем, затрачиваемым программой на арифметику с плавающей запятой, или пропорционально ему, и поэтому в обоих случаях число флопов с точностью до постоянного множителя примерно совпадает с реальным временем работы программы. Тем не менее, информация о количестве флопов не всегда приводит к правильным выводам, особенно на современных компьютерных архитектурах. См. по этому поводу, например, раздел 2.3.4.

Учебно-исследовательская библиотека NL [33] представляет собой пакет функций, написанных на языке C, реализующих алгоритмы решения систем линейных уравнений:

• решение систем с помощью LU -, QR-разложений и разложения Холецкого;

1.6. Библиотека NL • сингулярное разложение; задача наименьших квадратов;

• прямые методы решения ленточных систем;

• операции с разреженными матрицами;

• прямые и итерационные методы решения разреженных систем.

Подробная инструкция по установке и эксплуатации библиотеки и справочник по всем функциям находятся в [33].

Приведем список всех заголовочных файлов библиотеки:

• nl.h подключает все остальные заголовочные файлы библиотеки;

• lu.h содержит объявления функций, реализующих и использующих LU -разложение квадратной матрицы;

• chol.h содержит объявления функций, реализующих и использующих разложение Холецкого;

• qr.h содержит объявления функций, реализующих и использующих QR-разложение матрицы;

• svd.h содержит объявления функций, строящих и использующих сингулярное разложение матрицы;

• sparse.h содержит объявления функций, поддерживающих работу с разреженными матрицами;

• band.h содержит объявления функций для работы с ленточными матрицами;

• gallery.h содержит объявления функций, генерирующих часто встречающиеся в приложениях матрицы;

• util.h содержит объявления вспомогательных функций;

• qsort.h содержит объявления функций, реализующих алгоритм быстрой сортировки (используется при работе с разреженными матрицами).

30 Глава 1. Машинная арифметика и анализ ошибок Библиотека поддерживает вычисления только с вещественными векторами и матрицами с элементами типа double (вещественное число двойной точности). Так как исходный код доступен, то все функции легко переписать, чтобы они работали с вещественными матрицами, в которых элементами являются числа одинарной точности. Некоторые алгоритмы могут быть распространены и на комплексный случай.

Для представления плотных векторов и матриц в библиотеке действуют следующие соглашения. В библиотеке объявлен тип данных nl index :

typedef size t nl index ;

использующийся для представления индексов векторов и матриц.

Вектор представляется в виде одномерного массива:

Для работы с вектором одного такого объявления недостаточно. Необходимо сперва выделить память. Например, создадим вектор длины n = 100:

a = nl dvector create(n);

Библиотечная функция nl dvector create резервирует память под заданное число элементов типа double (с помощью функции malloc) и возвращает указатель на вектор. Теперь a[i ] i-й элемент вектора a. Нумерация элементов начинается с нуля.

Значения элементов после вызова функции nl dvector create не определены. Как только вектор становится не нужен, необходимо освободить занимаемую им память:

nl dvector free(a);

1.6. Библиотека NL В библиотеке есть функции для копирования векторов, создания векторов с инициализированными значениями элементов, ввода из файла и с консоли, вывода в файл и на консоль и др.

Поддерживаются также векторы с элементами типа nl index.

Информация о плотной матрице хранится в библиотеке в виде набора массивов, соответствующих строкам матрицы, и массива указателей на эти массивы:

double **A;

Прежде чем работать с матрицей, необходимо выделить память. Создадим, например, матрицу размера 100 200:

double **A;

nl index m = 100;

nl index n = 200;

A = nl dmatrix create(m, n);

Библиотечная функция nl dmatrix create резервирует память под заданное число элементов типа double (с помощью функции malloc) и возвращает указатель на матрицу. Теперь A[i ][j ] элемент матрицы, стоящий в i-й строке и j-м столбце, причем нумерация строк и столбцов начинается с нуля. Значения элементов после вызова функции nl dmatrix create не определены. Для освобождения занимаемой памяти необходимо вызвать функцию nl dmatrix free:

nl dmatrix free(A, m);

вторым аргументом функции nl dmatrix free указывается число строк матрицы A.

В библиотеке есть функции для копирования матриц, создания матриц с инициализированными значениями элементов, ввода/вывода и др. Подробную документацию см. в [33].

Заметим, что существует другой широко распространенный способ хранения плотных матриц, в котором все ее элементы 32 Глава 1. Машинная арифметика и анализ ошибок расположены с последовательных ячейках памяти, например, по строкам:

A = malloc(m * n * sizeof (double));

Доступ к элементу следующий: A[i * n + j ]. При таком способе представления матрицы, если необходима последовательная обработка ее элементов, при переходе от одной строки к следующей не возникает дополнительных издержек. Особенно этот эффект заметен на современных вычислительных системах с многоуровневой системой памяти. Мы, тем не менее, остановились на первом способе в силу его большей наглядности, хотя он и менее эффективен и требует дополнительной памяти. Заметим, что библиотека содержит функции для конвертирования из первого способа во второй и обратно.

Еще раз отметим, что нумерация элементов векторов и нумерация строк и столбцов матриц начинается с нуля (как принято в C). Это надо иметь в виду, так как в книге при описании алгоритмов используется нумерация, начинающаяся с 1.

Прямые методы решения Метод решения систем линейных уравнений называется прямым (точным), если при отсутствии ошибок округления он находит точное решение за конечное число шагов.

В данной главе рассматриваются прямые методы решения квадратной определенной системы линейных уравнений Ax = b, т. е. системы с невырожденной матрицей A порядка n. Все методы данной главы основаны на алгоритме Гаусса. Переопределенные системы рассматриваются в главе 3, итерационные методы в главе 5.

Подробно прямые методы решения линейных систем освещаются, например, в [2, 4, 5].

Опишем метод Гаусса (а точнее: метод Гаусса без выбора ведущего элемента, или c диагональным выбором ведущего элемента) решения системы линейных уравнений Ax = b с квадратной невырожденной матрицей A. Обозначим 34 Глава 2. Прямые методы решения линейных систем Прямой ход метода состоит из n 1 шагов. На первом шаге при условии a11 = 0 для исходной системы

вычисляются n 1 множителей Гаусса а затем для каждого i = 2,..., n из i-го уравнения системы (12) вычитается 1-е уравнение, умноженное на µi1. После этих преобразований система принимает вид Опишем вычисления на k-м шаге (k = 1, 2,..., n 1). К его началу исходная система преобразована к виду a(0) x1 + a(0) x2 +... + a(0) xk +... + a(0) xn = b(0), 2.2. Метод Гаусса или, в матричных обозначениях, На k-м шаге при условии вычисляются множители Гаусса а затем для каждого i = k + 1,..., n из i-го уравнения системы (14) вычитается k-е уравнение, умноженное на µik. В результате система преобразуется к эквивалентному виду Элемент akk называется ведущим, или главным, элементом k-го шага. По окончании (n 1)-го шага получим верхнетреугольную систему 36 Глава 2. Прямые методы решения линейных систем или, в матричных обозначениях, U x = f, которую легко решить с помощью обратного хода. Обратный ход заключается в нахождении неизвестных по формулам Нетрудно проверить, что прямой ход метода Гаусса использует (2/3)n3 + O(n2 ) операций с плавающей точкой, в то время как обратный ход требует n2 + O(n) операций с плавающей точкой.

Описанный вариант метода Гаусса осуществим лишь при условии, что ведущие элементы не равны нулю.

LU -разложением матрицы A называется равенство где L нижнетреугольная матрица с единичной диагональю, верхнетреугольная матрица.

Пусть A = LU. Решение исходной системы Ax = b находится посредством последовательного решения треугольных систем Ly = b (прямая подстановка), U x = y (обратная подстановка). Решение верхнетреугольной системы рассматривалось в предыдущем разделе. Нижнетреугольная система Ly = b решается аналогично. Для ее решения также требуется n2 + O(n) операций с плавающей точкой.

Пусть требуется решить системы Ax = b(1),..., Ax = b(s), с одинаковыми левыми частями, но с разными правыми. В этом случае матрицу A достаточно разложить один раз. После того как разложение найдено, для решения каждой системы достаточно выполнить 2n2 + O(n) операций с плавающей точкой.

2.3. LU -разложение Один из способов нахождения LU -разложения основан на алгоритме Гаусса. Преобразования k-го шага алгоритма Гаусса эквивалентны домножению системы (14) слева на матрицу где µij множители Гаусса, вычисляемые по формуле (16).

Таким образом, прямой ход метода Гаусса заключается в преобразованиях, переводящих систему Ax = b в эквивалентную систему где U = M (n1) M (n2)... M (1) A, f = M (n1) M (n2)... M (1) b.

Обозначим L1 = M (n1) M (n2)... M (1). Легко проверить, что 38 Глава 2. Прямые методы решения линейных систем Таким образом, A = LU.

Заметим, что элементы матрицы L, расположенные ниже главной диагонали, и элементы матрицы U на главной диагонали и выше нее можно записывать в ячейки, в которых хранятся соответствующие элементы матрицы A.

Мы приходим к следующему алгоритму LU -разложения с помощью исключения по столбцу:

Заметим, что для осуществимости данной процедуры необходимо выполнение неравенства (15) для всех k = 1, 2,..., n1.

Процедура использует (2/3)n3 + O(n2 ) операций с числами с плавающей точкой.

Рассмотрим другую схему нахождения LU -разложения.

Основная идея алгоритма LU -разложения с помощью исключения по столбцу заключается в том, что на каждой (k-й) итерации ведущая (k-я) строка с подходящими множителями вычитается из строк, лежащих ниже, чтобы занулить все элементы матрицы, расположенные в k-м столбце ниже диагонали. Между тем возможно исключение по строке: на каждой (i-й) итерации можно вычитать из i-й строки все строки, расположенные выше, умноженные на подходящие коэффициенты, 2.3. LU -разложение так, чтобы занулить все элементы i-й строки левее диагонали.

Как и в исключении по столбцу, элементы ниже главной диагонали матрицы L и элементы матрицы U на главной диагонали и выше нее можно вычислять на месте матрицы A. Это сразу приводит нас к следующему алгоритму LU -разложения с помощью исключения по строке:

Как и в исключении по столбцу, приведенная схема использует (2/3)n3 + O(n2 ) операций с числами с плавающей точкой.

Приведем еще один способ нахождения LU -разложения матрицы. Этот способ называют компактной схемой, или методом Краута–Дулитла.

Предположим, что матрица A допускает разложение вида A = LU. Здесь, как и ранее, L = (lij ) нижнетреугольная матрица с единичной диагональю, U = (uij ) верхнетреугольная матрица, поэтому 40 Глава 2. Прямые методы решения линейных систем Матричное равенство A = LU эквивалентно системе n2 равенств В силу (19) из этих равенств получаем откуда легко получить расчетные формулы:

Как и для обычной схемы, элементы ниже главной диагонали матрицы L и элементы на главной диагонали матрицы U и выше нее можно вычислять на месте матрицы A.

Получаем следующий алгоритм построения LU -разложения матрицы с помощью компактной схемы (метод Краута– Дулитла) 2.3. LU -разложение Как и другие описанные здесь алгоритмы построения LU разложения, приведенная схема использует (2/3)n3 +O(n2 ) операций с числами с плавающей точкой.

2.3.4. Замечания о производительности Каждый из трех приведенных ранее алгоритмов построения LU -разложения: исключение по столбцу, исключение по строке, компактная схема заключается в выполнении трех вложенных циклов. Вычисление элемента aij в самом внутреннем цикле производится по одной и той же формуле Это позволяет назвать метод исключения по столбцу kijалгоритмом (по порядку следования переменных цикла), метод исключения по строке ijk-алгоритмом, компактную схему jik-алгоритмом. Заметим, что возможны также kji-, ikj- и jkiварианты алгоритма LU -разложения. Например, kji-вариант получается из kij-метода простой перестановкой двух внутренних циклов: в алгоритме на стр. 35 нужно переставить строки, отмеченные () и (). Чуть более сложно получаются ikj- и jki-модификации.

Заметим, что во всех шести вариантах алгоритма LU -разложения производятся одни и те же вычисления, но порядок этих вычислений различен. Понятно, что все методы выполняют одинаковое число арифметических операций с плавающей точГлава 2. Прямые методы решения линейных систем кой. Что же касается производительности реальных программных реализаций этих методов, то оказывается, что особенности архитектур современных компьютеров вносят сюда свои коррективы. Из-за наличия в современных компьютерах многоуровневого кэша последовательное извлечение и сохранение чисел, расположенных в соседних ячейках памяти, производится в несколько раз быстрее, чем последовательное выполнение тех же операций с элементами, расположенными в памяти далеко друг от друга. Предположим, что элементы матрицы хранятся в памяти по строкам: соседние элементы одной строки размещены в последовательных ячейках памяти. Ввиду вышеизложенного представляется, что наиболее быстрыми будут программные реализации алгоритмов, в которых в самом внутреннем цикле происходит изменение соседних элементов одной строки, т. е. алгоритмы с внутренним j-циклом: kij- и ikj-методы.

Обратим внимание, что все шесть вариантов алгоритма приводят к одинаковым LU -разложениям даже при наличии ошибок округления. Это замечание относится только к обычному режиму вычислений, когда все арифметические операции производятся с одной и той же точностью. Иногда в компьютере имеется возможность накапливать сумму произведений в регистре с длиной мантиссы, в два раза большей обычной (режим накопления). Это, например, позволяет вычислять выражения вида (20), (21) с большей точностью. Таким образом, точность вычисления LU -разложения с помощью компактной схемы можно повысить [2].

2.4. Схемы выбора ведущего элемента 2.4.1. Стратегия частичного выбора Алгоритм Гаусса и алгоритм нахождения LU -разложения без выбора ведущего элемента, описанные в двух предыдуСхемы выбора ведущего элемента щих разделах, осуществимы лишь при выполнении неравенства (15), т. е. тогда, когда ведущие элементы на каждом шаге ненулевые.

Рассмотрим k-й шаг алгоритма Гаусса. Пусть Переставим s-ю и k-ю строки и выполним те же преобразования, что и в алгоритме Гаусса, но с новым ведущим элементом. Приведенная схема называется стратегией выбора ведущего элемента по столбцу, или стратегией частичного выбора.

Для ее осуществления нужно выполнить (1/2)n2 + O(n) дополнительных операций сравнения.

Матрицу P, полученную из единичной перестановкой ее строк, назовем матрицей-перестановкой.

Перестановка s-й и k-й строк матрицы эквивалентна домножению ее слева на матрицу-перестановку Psk, которую можно получить из единичной матрицы перестановкой ее s-й и k-й строк. Таким образом, прямой ход метода Гаусса с выбором ведущего элемента по столбцу заключается в преобразованиях, приводящих систему Ax = b в эквивалентную систему M (n1) P (n1) M (n2) P (n2)... M (1) P (1) Ax = где M (k) матрицы вида (18), а P (k) матрицы-перестановки.

Таким образом, Можно показать, что из равенства (22) следует 44 Глава 2. Прямые методы решения линейных систем где P = P (n1) P (n2)... P (1) матрица-перестановка, L нижнетреугольная матрица с единичной диагональю. Разложение (23) также называется LU -разложением матрицы A. Очевидно, оно существует для любой невырожденной матрицы.

Если разложение P A = LU найдено, то для решения системы Ax = b достаточно вычислить произведение P b (которое получается перестановкой компонент вектора b), а затем решить две треугольные системы Ly = P b, U x = y.

Другая схема выбора ведущего элемента называется стратегией выбора ведущего элемента по всей матрице, или стратегией полного выбора.

Рассмотрим k-й шаг алгоритма Гаусса. Пусть Переставим местами s-ю и k-ю строки, а также t-й и k-й столбцы и выполним те же преобразования, что и в алгоритме Гаусса, но с новым ведущим элементом. Разумеется, при перестановке столбцов мы должны переобозначить переменные. Это и есть стратегия полного выбора. Заметим, что для ее осуществления нужно выполнить уже (1/3)n3 + O(n2 ) дополнительных операций сравнения, что уже сравнимо с количеством операций с плавающей точкой, осуществляемых при работе алгоритма.

Легко показать, что стратегия полного выбора приводит к разложению P1 AP2 = LU, где P1, P2 матрицы-перестановки, L нижнетреугольная матрица с единичной диагональю, U верхнетреугольная матрица.

Исследуем ошибки, возникающие при решении системы линейных уравнений Ax = b на основе LU -разложения. Алгоритм 2.5. Анализ ошибок состоит из двух этапов. На первом мы ищем LU -разложение матрицы A. На втором решаем две треугольные системы Ly = b, U x = y. Будем использовать обратный анализ ошибок.

2.5.1. Aнализ ошибок в LU -разложении Для определенности рассмотрим компактную схему построения LU -разложения. Анализ других вариантов алгоритма такой же. Также для удобства будем считать, что выбор ведущих элементов произведен заранее.

Пусть при вычислениях на компьютере были найдены нижнетреугольная с единичной диагональю матрица L и верхнетреугольная матрица U, для которых выполнено точное равенство A + A = LU, где A эквивалентное возмущение.

Обозначим aij = (A + A)ij. Оценим сверху A.

Cогласно (10) формула (20) при i j приводит к равенству где |q | M, |q | M. Преобразуем его правую часть:

Из последнего равенства выразим aij :

Для оценивания встречающихся произведений воспользуемся легко проверяемым двусторонним неравенством 46 Глава 2. Прямые методы решения линейных систем справедливым, если |j | M (j = 1, 2,..., d) и dM 1. Для IEEE-арифметики одинарной точности последнее ограничение означает d 107, что абсолютно не обременительно, так как при реальных вычислениях обычно d Воспользовавшись неравенством (24), теперь получаем откуда для всех i, j при i j с точностью до слагаемого O(2 ) Таким образом, |(A)ij | nM (|L||U |)ij + O(2 ).

Из формулы (21) аналогично получается точно такое же неравенство и для i j.

Итак, |A| nM |L|·|U |, откуда |A| nM |L| · |U | + O(2 ). Если норма матрицы зависит только от абсолютных значений ее элементов (это справедливо для манхеттенской, чебышевой и фробениусовой норм, но не для спектральной нормы), то получаем 2.5.2. Анализ ошибок при решении треугольных систем Найденное LU -разложение матрицы A (напомним, что оно является точным для матрицы A + A = LU ) используется при решении треугольных систем Ly = b, U x = y. Из-за ошибок округления при решении этих систем вместо точных решений 2.5. Анализ ошибок y, x будут найдены приближенные y +y, x+x соответственно. Имеем где L, U, A эквивалентные возмущения. Оценим норму каждого из них.

Для решения треугольной системы Ly = b воспользуемся формулами Cогласно (10) имеем где |q | M, |q | M (q = 1, 2,..., i 1). Преобразуем правую часть полученного равенства:

преобразований j= 48 Глава 2. Прямые методы решения линейных систем или согласно (24) где |j | jM +O(2 ) (j = 1, 2,..., i). Итак, вычисленное решеM ние y+y является точным решением возмущенной системы (L + L)(y + y) = b, где (L)ij jM |lij | + O(2 ), т. е.

Проводя аналогичные выкладки для системы U x = y, можно получить где A + A = LU, A = LU + LU + LU. Теперь из (25), (26), (27), получаем Для манхеттенской, чебышевой и фробениусовой (но не спектральной) норм эквивалентного возмущения теперь имеем 2.5. Анализ ошибок Неравенство (28) показывает, что для того, чтобы найденное приближенное решение было решением близкой системы, необходимо, чтобы элементы в L и U не росли слишком сильно.

Если используется алгоритм разложения без выбора ведущего элемента, то, очевидно, элементы в L и U могут сильно расти, что приводит к численной неустойчивости.

Теперь рассмотрим алгоритм с частичным выбором ведущего элемента. Экспериментальные данные показывают, что в этом случае почти всегда A L · U.

Легко видеть, что при использовании этой схемы выбора ведущего элемента |lij | 1 и поэтому Остается исследовать возможность роста коэффициентов матрицы U. Обозначим тогда Из экспериментов следует [5], что в среднем gpp, по-видимому, растет примерно как n2/3 или даже как n1/2. Однако существуют примеры матриц, для которых gpp = 2n1, и можно показать, что gpp 2n1. Например, для матриц A вида 50 Глава 2. Прямые методы решения линейных систем справедливо gpp = 2n1.

Собирая вместе (28), (29), (30), получаем Теперь можно получить оценки нормы невязки и точности решения системы с помошью LU -разложения со схемой частичного выбора ведущего элемента:

Полученные оценки слишком пессимистичны. На практике обычно множитель 3gpp n3 можно заменить на cn, где c некоторая умеренная константа (в пределах 10) [10]. Поэтому Из полученных приближенных равенств сделаем два важных практических вывода. Во-первых, метод LU -разложения с частичным выбором ведущего элемента дает приближенное решение x + x с малой относительной невязкой не зависимо от числа обусловленности матрицы A. Во-вторых, если M 10p, cond A 10q, то метод LU -разложения с частичным выбором ведущего элемента дает приближенное решение x + x с приблизительно p q верными десятичными знаками.

Теперь рассмотрим стратегию полного выбора ведущего элемента. Дж.Х. Уилкинсону принадлежит оценка по-видимому, тоже слишком завышенная. На практике в среднем gcp ведет себя примерно как n [5].

2.6. Оценка числа обусловленности 2.6. Оценка числа обусловленности Предположим, что система Ax = b решена посредством разложения P A = LU и мы хотим установить количество правильных цифр в вычисленном решении x. Для этого нам нужно оценить число обусловленности. Вычисление A 1 или A не вызывает проблем. Сложной задачей является вычисление A1 1 или A1. Для этого можно явно найти A1, однако даже если разложение P A = LU известно, это требует O(n3 ) операций с плавающей точкой.

Главная проблема оценки числа обусловленности состоит в том, чтобы при наличии разложения P A = LU оценить обусловленность за O(n2 ) операций. Алгоритмы, рещающие эту задачу, называют оценщиками. Известно много разных оценщиков. Рассмотрим один из них.

Оценщик Хэйджера [5] строит нижнюю оценку для A1 1.

На практике построенная оценка отличается от истинного значения нормы обратной матрицы в умеренное (порядка 10) число раз.

Обозначим B = A1 и пусть f (x) = Bx 1. Задачу нахождения величины A1 1 можно сформулировать как задачу максимизации функции f (x) при ограничении x 1 1.

Заметим, что упомянутый максимум достигается на одном из векторов ej, однако простой перебор всех векторов ej (вместе с подсчетом Bej 1 ) требует O(n3 ) операций. Легко показать, что функция f (x) выпукла и область, описываемая данным ограничением, также выпукла.

Идея оценщика Хэйджера заключается в том, что для максимизации используется градиентный алгоритм. Гарантированно этот алгоритм найдет лишь точку локального максимума, что дает нижнюю оценку числа обусловленности. Алгоритм перемещается от одного вектора ej к другому в направлении градиента f функции f (x).

52 Глава 2. Прямые методы решения линейных систем Пусть w = Bx. Обозначим sign w вектор, в котором j-я компонента равна 1, 1 или 0, если wj положительно, отрицательно или равно нулю соответственно. Можно показать [5], что градиент f (x) равен B T sign w.

Приведем формальное описание алгоритма Хэйджера, строящего нижнюю оценку для B 1. В алгоритме встречаются выражения вида Bx и B T y. Так как B = A1, то речь идет о решении линейных систем. Мы предполагаем, что нам известно LU -разложение матрицы A, поэтому решение этих систем требует O(n2 ) операций с плавающей точкой.

Можно показать [5], что на выходе алгоритма мы получим локальный максимум функции Bx 1.

Оценщик Хэйджера можно использовать при вычислении оценок FERR и BERR [5].

О других оценщиках см. в [4, 10].

Пусть найдено решение x системы Ax = b в арифметике одинарной точности. Нас может не удовлетворить его точность.

Если это так, то мы можем заново решить систему в арифметике двойной точности, что потребует не меньших затрат времеИтерационное уточнение ни и памяти. Есть более экономные способы получения более точного решения.

С помощью следущего алгоритма в арифметике двойной точности можно получить уточнение решения:

repeat решить систему Ad = r относительно вектора d Для решения системы Ad = r нужно воспользоваться уже имеющимся LU -разложением.

Итерации в алгоритме продолжаются до достижения требуемой точности. Если ошибок округления нет, то за одну итерацию мы получим d = 0, итерации завершатся и x является точным решением системы Ax = b.

Можно показать [5], что если cond(A)M 1, то в присутствии ошибок округления итерации приведенного алгоритма сойдутся к вектору x, для которого Таким образом, если cond(A)M значительно меньше 1, то решение можно вычислить с высокой точностью независимо от значения числа обусловленности.

Иногда возможность вычисления с удвоенной точностью отсутствует. Например, если вы изначально проводите вычисления с двойной точностью, то невязку r нужно вычислять с учетверенной точностью, что возможно далеко не всегда. Оказывается, к хорошим результатам приводит уже выполнение одного шага итерационного уточнения в арифметике обычной точности [5].

наименьших квадратов Пусть A матрица размера m n, а b столбец высоты m. Линейная задача наименьших квадратов заключается в нахождении такого вектора x, называемого псевдорешением системы Ax = b, на котором достигается минимум Таким образом, псевдорешение доставляет минимум евклидовой нормы невязки заданной системы Ax = b. Легко видеть, что если система Ax = b совместна, то любое ее решение является ее псевдорешением. Если m n, то система называется переопределенной.

Псевдорешение единственно тогда и только тогда, когда rank A = n (столбцы матрицы A линейно независимы).

Задача наименьших квадратов возникает при аппроксимации данных. Предположим, что имеются пары значений Пусть f1 (x), f2 (x),..., fn (x) заданные функции. Требуется найти функцию вида 3.1. Постановка задачи для которой f (xi ) yi (i = 1, 2,..., m), т. е. функцию, приближающую данные (33). Более точно: требуется найти значения параметров 1,..., n, на которых достигается минимум Легко видеть, что набор этих параметров является псевдорешением системы

И наоборот, если 1,..., n псевдорешение указанной системы, то функция f (x), определяемая по формуле (34), является решением задачи аппроксимации.

Часто в качестве функций f1 (x), f2 (x),..., fn (x) рассматриваются многочлены 1, x,..., xn1. Таким образом, речь идет об аппроксимации данных полиномами.

Мы рассмотрим следующие методы решения линейной задачи наименьших квадратов:

1) нормальные уравнения;

2) QR-разложение:

а) ортогонализация Грама–Шмидта, б) отражения Хаусхолдера, в) вращения Гивенса, 3) сингулярное разложение (SVD ).

66 Глава 3. Линейная задача наименьших квадратов Далее мы, как правило, сосредоточимся на случае полного ранга (rank A = n). В разделе 3.8 рассматривается (более трудная) задача наименьших квадратов неполного ранга.

Подробно линейная задача наименьших квадратов рассматривается, например, в [4, 5, 13].

3.2. Обусловленность задачи наименьших квадратов В разделе 1.2 было введено число обусловленности квадратной матрицы. Дадим определение числа обусловленности прямоугольной матрицы. Пусть матрица A имеет размеры m n, причем m n. Спектральным числом обусловленности матрицы A называется величина Это определение согласуется с определением cond2 A для квадратной матрицы (см. (2) на стр. 11). Если m n, то по определению cond2 A = cond2 AT.

Предположим, что rank A = n. Пусть x решение задачи наименьших квадратов min Ax b 2, r невязка, т. е.

r = b Ax, x + x решение возмущенной задачи наименьших квадратов min (A + A)x (b + b) 2. Предположим, что возмущение не слишком велико, а именно, 1/ cond2 A, где Можно показать [4] (см. также задачу 3 из раздела 9.3 на стр.

248), что 3.3. Нормальные уравнения Величина называется числом обусловленности для задачи наименьших квадратов. Получаем, что Итак, condLS (A, b) является мерой чувствительности задачи наименьших квадратов, показывающей, как сильно решение этой задачи зависит от возмущения входных данных.

Формулу (35) можно интерпретировать следующим образом. Если угол мал, то condLS (A, b) 2 cond2 A, т. е. задача наименьших квадратов обусловлена примерно так же, как и задача решения систем линейных уравнений.

Если угол не мал, но и не близок к /2, то condLS (A, b) cond2 A.

Если /2, т. е. задача близка к вырожденной, то величина condLS (A, b) растет неограничено, даже если cond2 A имеет умеренные значения.

Покажем, что множество псевдорешений системы Ax = b совпадает с множеством решений системы называемой системой нормальных уравнений.

Имеем 68 Глава 3. Линейная задача наименьших квадратов Если x псевдорешение, то Ax 2 + 2xT (AT Ax AT b) для любого x, откуда AT AxAT b = 0. Обратно, если AT Ax = AT b, то A(x + x) b 2 Ax b 2 для любого x, т. е. x псевдорешение.

Заметим, что матрица AT A квадратная и неотрицательно определенная. Система (36) всегда совместна. Она имеет единственное решение тогда и только тогда, когда псевдорешение исходной системы единственно. Последнее имеет место в том и только в том случае, когда столбцы матрицы A линейно независимы (т. е. rank A = n). В этом (и только в этом) случае матрица AT A положительно определена и для решения системы (36) мы можем применить разложение Холецкого. Тогда на вычисление матричного произведения AT A требуется mn2 + O(n2 ) операций с плавающей точкой (достаточно вычислить только элементы над диагональю), а на решение системы (36) требуется n3 /3 + O(n2 ) операций с плавающей точкой. Так как m n, то основная масса времени тратится на формирование матрицы AT A.

Недостатком метода нормальных уравнений является плохая обусловленность получаемой системы. Действительно, справедливо равенство cond2 (AT A) = cond2 A. Используя результаты анализа ошибок для метода Холецкого, приходим к выводу, что относительная ошибка получаемого решения задачи наименьших квадратов не превосходит некоторого умеренного кратного величины M cond2 A. Следует ожидать, что при решении нормальных уравнений с небольшой невязкой исходной системы будет потеряно вдвое больше разрядов, чем при использовании QR-разложения или сингулярного разложения.

Тем не менее использование нормальных уравнений заманчиво при решении линейной задачи наименьших квадратов при малых значениях n: например, при n = 2 система (36) имеет размеры 2 2.

3.4. QR-разложение Можно доказать, что для произвольной (m n)-матрицы A ранга n существуют и единственны матрицы Q размера m n с ортонормированными столбцами (т. е. QT Q = E) и верхнетреугольная матрица R размера n n с положительными диагональными элементами, такие, что Разложение (37) называется QR-разложением1 матрицы A.

Сперва покажем, как решать линейную задачу наименьших квадратов min Ax b 2, если QR-разложение матрицы A известно. Методы построения самого QR-разложения будут рассмотрены позднее.

Итак, QT Q = E. Дополним столбцы матрицы Q до ортонормированной системы m векторов2. Из построенных столбцов составим матрицу S. Имеем S T Q = 0. Матрица (Q, S) ортогональная, т. е. (Q, S)1 = (Q, S)T, поэтому для любого вектора x справедливо (Q, S)T x 2 = x 2. Теперь получаем:

Иногда QR-разложением называют представление матрицы A в виде произведения A = QR, в котором матрица Q квадратная ортогональная порядка m, а матрица R верхнетреугольная размера mn. Чтобы такое разложение получить по разложению (37), достаточно столбцы матрицы Q дополнить до ортонормированной системы m векторов, из построенных столбцов составить матрицу S и положить Обратный переход осуществляется отбрасыванием в Q последних m n столбцов, а в R последних m n строк.

Ср. предыдущее примечание.

70 Глава 3. Линейная задача наименьших квадратов Минимум достигается тогда и только тогда, когда первое слагаемое равно нулю. Итак, решением линейной задачи наименьших квадратов является решение треугольной системы при этом Заметим, что матрица R невырождена тогда и только тогда, когда rank A = n. В этом случае x = R1 QT b.

Рассмотрим несколько методов построения QR-разложения.

3.4.1. Ортогонализация Грама–Шмидта ческий метод построения QR-разложения. Он заключается в построении ортонормированного базиса q1,..., qn подпространства, натянутого на столбцы a1,..., an, причем Матрица, составленная из столбцов q1,..., qn, это матрица Q. Матрица перехода от базиса q1,..., qn к a1,..., an это матрица R, т. е. aj = i=1 rij qi Опишем этот алгоритм подробнее.

3.4. QR-разложение Алгоритм использует 2mn2 + O(mn) операций с плавающей точкой.

Заметим, что присваивание rij = qi aj в строке () привеT денного алгоритма можно заменить на rij = qi qj. Полученный таким образом новый алгоритм (в точной арифметике эквивалентный исходному) называется модифицированным процессом Грама–Шмидта.

При реальных вычислениях на компьютере (при наличии ошибок округления) процесс ортогонализации Грама–Шмидта (в том числе модифицированный) вместо векторов q1, q2,..., qn возвратит векторы q1, q2 + q2,..., qn + qn. Например, в режиме накопления где ak 2 2M ak 2, т. е. линейные оболочки полученных векторов совпадают с линейными оболочками слабо возмущенных исходных векторов [2].

Итак, в отношении выполнения равенства (38) процесс ортогонализации Грама–Шмидта ведет себя крайне устойчиво.

К сожалению, этого нельзя сказать о выполнении условия ортогональности полученной системы. В большинстве случаев процесс Грама–Шмидта приводит к значительному нарушению ортогональности. Модифицированный процесс Грама–Шмидта чуть лучше, но и он не свободен от указанного недостатка.

Пусть Q + Q реально вычисленная матрица. Можно показать [4], что 72 Глава 3. Линейная задача наименьших квадратов где E 2 M cond2 A.

Для исправления отмеченного недостатка используют переортогонализацию [2].

Тем не менее, как оказалось, ортогональность столбцов матрицы Q не играет решающей роли в задаче наименьших квадратов. Можно показать, что алгоритм решения задачи наименьших квадратов с помощью модифицированного процесса ортогонализации Грама–Шмидта обратно устойчив и, следовательно может быть использован на практике [4]. В любом случае, затраты в методах типа Грама–Шмидта чуть выше, чем в методе Хаусхолдера.

Метод отражений Хаусхолдера один из самых распространенных методов нахождения QR-разложения.

Пусть v n-мерный ненулевой вектор-столбец. Квадратная матрица P порядка n вида называется отражением Хаусхолдера или просто отражением. Вектор v называется вектором Хаусхолдера. Умножению матрицы P на вектор x можно дать следующую геометрическую интерпретацию: вектор P x получается отражением вектора x относительно гиперплоскости, ортогональной вектору v. Легко проверить, что матрица Хаусхолдера симметрична и ортогональна.

Пусть x = 0. Найдем отражение P, такое, что P x = e1 для какого-либо. Имеем и поэтому v коллинеарен xP x = xe1. Так как v определен с точностью до ненулевого множителя, то можно положить v = 3.4. QR-разложение x e1. Ввиду ортогональности матрицы P имеем P x = x 2, откуда = ± x 2. Поэтому Итак, в качестве вектора Хаусхолдера v можно взять вектор, вычисляемый по формуле (39), или любой ненулевой, ему коллинеарный.

С алгоритмом построения вектора Хаусхолдера связаны некоторые важные детали. Одна из них выбор знака в формуле (39). Если x почти коллинеарен e1, то вектор v = x sign(x1 ) x 2 e1 имеет малую норму. Вследствие этого возможно появление большой относительной ошибки при вычислении множителя 2/v T v. Эту трудность можно обойти, взяв с тем же знаком, что и первая компонента вектора x:

Это обеспечивает выполнение неравенства v 2 x 2 и гарантирует почти точную ортогональность вычисленной матрицы Другая деталь связана с выбором множителя для вектора, вычисляемого по формуле (39). Мы будем использовать такой множитель, что v1 = 1. Это несколько упрощает способ хранения ортогональных матриц, вычисленных по методу отражений Хаусхолдера.

Суть алгоритма QR-разложения на основе отражений сначала поясним на примере. Пусть m = 5, n = 4 и матрицы P1 и P2 таковы, что 74 Глава 3. Линейная задача наименьших квадратов Покажем, как на месте элементов f43, f53 получить нули. Найдем матрицу отражения P3, такую, что Если P3 = diag(1, 1, P3 ), то Выполнив n таких шагов, получим верхнюю треугольную матрицу R = Pn Pn1... P1 A. Откуда, так как матрицы Pj ортогональные и симметричные, имеем A = QR, где Q = P1 P2... Pn. Пусть Q состоит из первых n столбцов матрицы Q, а R состоит из первых n строк матрицы R, тогда A = QR. Это и есть искомое разложение.

Матрицы Pj не требуется формировать в явном виде: достаточно хранить соответствующий вектор vj, который можно записывать в j-й столбец матрицы A.

Для решения задачи наименьших квадратов необходимо вычислить произведение QT b = Pn Pn1... P1 b.

Трудоемкость алгоритма построения QR-разложения на основе отражений составляет 2n2 m (2/3)n3 + O(mn) арифметических операций с плавающей точкой. Для получения решения задачи наименьших квадратов при известном QR-разложении требуется только O(mn) арифметических операций.

Пусть при реальных вычислениях на компьютере в результате ошибок округления вместо Q получена матрица Q + Q.

3.4. QR-разложение Можно показать [4], что где E 2 M.

Можно показать [13], что при при наличии ошибок округления вычисленное по методу Хаусхолдера решение x+x задачи min Ax b 2 является точным решением близкой задачи где Как мы видим, относительная ошибка эквивалентного возмущения пропорциональна умеренному кратному числа M. Используя теорию возмущений из раздела 3.2, получаем, что при малой невязке относительная ошибка решения задачи наименьших квадратов примерно пропорциональна M cond2 A, а при большой невязке M cond2 A.

На выходе диагональ и верхняя треугольная часть матрицы A содержат элементы матрицы R. Матрица Q представлена элементами ниже диагонали и вектором t. Вектор t имеет длину n. Матрица Q определяется произведением матриц отражений Хаусхолдера.

Вращением Гивенса, или просто вращением, называется матрица 76 Глава 3. Линейная задача наименьших квадратов Умножению матрицы R на вектор x можно дать следующую геометрическую интерпретацию: вектор Rx получается вращением вектора x на угол. Вращение n-мерного пространства в плоскости xi, xj определяет матрица Очевидно, что матрица вращения ортогональна.

Если Таким образом, с помощью домножения слева вектора x на матрицу вращения R можно занулить одну (i-ю) компоненту.

Заметим, что при этом все компоненты, кроме i-й и j-й, не изменяются.

Алгоритм QR-разложения на основе вращений аналогичен алгоритму, использующему отражения, с тем отличием, что 3.4. QR-разложение элементы матрицы A при обработке j-го столбца в нем зануляются по одному (для определенности, например, сверху вниз).

Чтобы перейти, например, от матрицы производим умножения:

Для хранения элементов найденного разложения дополнительная память не нужна. Информация о матрице вращения, с помощью которой мы занулили соответствующий элемент, хранится на его месте.

Для этого применяют следующий метод [5]. Пусть s = sin, c = cos. Если |s| |c|, то в память записывается число s · sign(c), в противном случае записывается sign(s)/c.

По хранимому значению p восстановление чисел s и c осуществляется следующим образом. Если |p| 1, то s = p и c = 1 s2, в противном случае c = 1/p и s = 1 c2. Заметим, что вместо исходных s, c мы можем восстановить s, c, но для наших целей приемлема и такая пара.

78 Глава 3. Линейная задача наименьших квадратов Метод вращений в два полтора раза медленнее метода отражений. Существует метод быстрых вращений [4], который, все-таки, медленнее метода отражений.

Обратный анализ ошибок для метода Гивенса решения задачи наименьших квадратов приводит к оценкам, аналогичным (41) [4]. Относительная ошибка эквивалентного возмущения пропорциональна умеренному кратному числа M.

3.5. QR-разложение в библиотеке NL Метод QR-разложения на основе отражений Хаусхолдера реализован в функции qr decomp. После того как разложение найдено, для решения квадратной совместной системы линейных уравнений необходимо воспользоваться функцией qr solve.

Приведем пример.

Листинг xqr.c Пример использования функций из модуля qr.h Решение системы линейных уравнений Ax = b с помощью QR-разложения, #include “nl.h” int main() double *A, *b, *t;

3.5. QR-разложение в библиотеке NL A = nl dmatrix create(n, n);

b = nl dvector create(n);

t = nl dvector create(n);

A[0] =.1; A[1] =.5; A[2] =.6; b[0] = 1.2;

A[3] =.2; A[4] =.7; A[5] =.9; b[1] = 1.8;

A[6] =.3; A[7] = 1.1; A[8] = 1.3; b[2] = 2.7;

printf (“Матрица A:\n”);

nl dmatrix print(A, n, n, NULL);

qr decomp(A, n, n, t);

printf (“\nQR-разложение:\n”);

nl dmatrix print(A, n, n, NULL);

printf (“\nВектор t:\n”);

nl dvector print(t, n, NULL);

printf (“\nВектор b:\n”);

nl dvector print(b, n, NULL);

qr solve (A, n, t, b);

printf (“\nРешение системы Ax = b:\n”);

nl dvector print(b, n, NULL);

nl dmatrix free(A);

nl dvector free(b);

nl dvector free(t);

return 0;

80 Глава 3. Линейная задача наименьших квадратов Результаты работы программы:

Матрица A:

QRразложение:

0.374 1.390 1. 0.422 0.136 0. 0.633 0.410 0. Вектор t:

Вектор b:

Решение системы Ax = b:

1.000 1.000 1. Для решения задачи наименьших квадратов необходимо воспользоваться функцией qr least squares. Приведем пример.

Листинг xqrls.c Пример использования функций из модуля qr.h Решение задачи наименьших квадратов с помощью QR-разложения 3.5. QR-разложение в библиотеке NL #include “nl.h” int main() double *A, *b, *t, *r, *work ;

A = nl dmatrix create(m, n);

b = nl dvector create(m);

t = nl dvector create(n);

r = nl dvector create(m);

work = nl dvector create(m);

A[0] A[3] A[6] A[9] printf (“Матрица A:\n”);

nl dmatrix print(A, m, n, NULL);

qr decomp(A, m, n, t);

printf (“\nQR-разложение:\n”);

nl dmatrix print(A, m, n, NULL);

printf (“\nВектор t:\n”);

nl dvector print(t, n, NULL);

printf (“\nВектор b:\n”);

82 Глава 3. Линейная задача наименьших квадратов nl dvector print(b, m, NULL);

qr least squares(A, m, n, t, b, r, work );

printf (“\nПсевдорешение системы Ax = b:\n”);

nl dvector print(b, n, NULL);

printf (“\nНевязки:\n”);

nl dvector print(r, m, NULL);

nl dmatrix free(A);

nl dvector free(b);

nl dvector free(t);

return 0;

Результаты работы программы:

QRразложение:

5.477 12.780 18. Вектор t:

3.6. Сингулярное разложение (SVD ) Вектор b:

Псевдорешение системы Ax = b:

0.458 0.125 0. Невязки:

Подробное описание используемых функций см. в документации [33].

3.6. Сингулярное разложение (SVD ) Рассмотрим матрицу A размера m n, где m n. Можно доказать, что существуют матрицы U размера mn, V размера n n и размера n n, такие, что n 0. Разложение (42) называется сингулярным разложением3 (SVD ). Столбцы u1,..., un матрицы U называются левыми сингулярными векторами, cтолбцы v1,..., vn матрицы V правыми сингулярными векторами, величины 1,..., n Иногда сингулярным разложением называют представление матрицы A в виде произведения A = U V T, в котором обе матрицы U и V квадратные ортогональные порядка m и n соответственно. Чтобы такое разложение получить по разложению (42), достаточно столбцы матрицы U дополнить до ортонормированной системы m векторов, из построенных столбцов составить матрицу S и положить U = (U, S), =, V = V.

Обратный переход осуществляется отбрасыванием в U последних m n столбцов.

84 Глава 3. Линейная задача наименьших квадратов сингулярными числами. Заметим, что каждая из двух систем:

левых сингулярных векторов и правых сингулярных векторов ортонормирована. Из (42) следует A = n i ui vi.

При m n необходимо рассматривать сингулярное разложение матрицы AT.

Нетрудно доказать следующие важные свойства сингулярного разложения:

1) Справедливо cond2 A = 1 /n. Если матрица A квадратная, то A 2 = 1. Если A квадратная и невырожденная, то A1 2 = 1/n.

3) j (AT A) = j, причем столбцы матрицы V образуют ортонормированную систему из собственных векторов матрицы AT A.

4) Если rank A = n, то (единственное) псевдорешение системы Ax = b равно x = V 1 U T b.

Докажем свойство 4). Так как U T U = E, то столбцы матрицы U можно дополнить до ортонормированной системы m векторов4. Из построенных столбцов составим матрицу S. Имеем S T U = 0. Матрица (U, S) ортогональная, т. е. (U, S)1 = T, поэтому для любого вектора x имеем (U, S)T x x 2. Теперь получаем:

Минимум достигается тогда и только тогда, когда первое слагаемое равно нулю, т. е. при x = V 1 U T b. Заметим, что при Ср. предыдущее примечание.

3.7. Сравнение методов этом Подведем небольшой итог и сравним методы решения задачи наименьших квадратов в случае rank A = n. Напомним несколько важных моментов.

• Чувствительность задачи наименьших квадратов к возмущениям входных данных при малой невязке примерно пропорциональна величине cond2 A, при большой невязке пропорциональна величине cond2 A (раздел 3.2).

• Метод нормальных уравнений дает решение, относительная ошибка которого линейно зависит от M cond2 A (раздел 3.3).

• Подходы, использующие преобразования Хаусхолдера и Гивенса, сингулярное разложение и модифицированный процесс Грама–Шмидта, решают близкую задачу: относительная погрешность эквивалентного возмущения пропорциональна умеренному кратному числа M. Поэтому при малой невязке относительная ошибка примерно пропорциональна величине M cond2 A, при большой невязке M cond2 A (раздел 3.4.2).

Можно заключить, что при малой невязке и умеренном числе обусловленности метод нормальных уравнений уверенно находит решение задачи наименьших квадратов. В этом случае его нужно предпочесть другим методам, так как он превосходит их по скорости. С ростом числа обусловленности решение, полученное методом нормальных уравнений, будет все менее и менее точным. Если cond2 A по порядку близко к 1/ M (или больше), то этот метод неприменим.

86 Глава 3. Линейная задача наименьших квадратов решающих задачу наименьших квадратов Нормальные уравнения Модифицированный Грама–Шмидта Ортогонализация Хаусхолдера Ортогонализация Гивенса Сингулярное разложение (Голуб–Райнш) Сингулярное разложение (Шан) При малой невязке и большом числе обусловленности нужно использовать методы, основанные на QR-разложении. Однако при cond2 A 1/M и они, конечно, не могут дать удовлетворительного решения.

При большой невязке и методы, основанные на QR-разложении, и метод, основанный на нормальных уравнениях, дают решения с примерно равными погрешностями, зависящими (линейно) от квадрата числа обусловленности матрицы A.

В заключение приведем таблицу 3.1, сравнивающую описанные методы по числу выполняемых в них арифметических операций с плавающей точкой [4, 13]. Для сингулярного разложения приводится информация о двух вариантах алгоритмов его нахождения.

Методы, основанные на сингулярном разложении, ведут себя наиболее стабильно, однако они наиболее медленные. Их разумно использовать, если задача неполного или почти неполного ранга.

Использование системы История развития системы Matlab (сокращение от Matrix Laboratory ) насчитывает более двух десятков лет. Классический Matlab был написан Кливом Моулером (университет Нью-Мехико) в 1977 г. Он представлял собой интерактивную матричную лабораторию, позволяющую вызывать подпрограммы из пакетов LINPACK и EISPACK. До 1984 г. появлялись новые некоммерческие версии системы. В 1984 г. К. Моулер, С. Бангерт и Дж. Литтл образовали фирму MathWorks. С этого момента начинают выходить коммерческие версии системы. К настоящему моменту (март 2005 г.) последней является версия Matlab 7.0 R 14 SP 2. Сейчас Matlab представляет собой мощный математический пакет со своим языком программирования, гибкими графическими возможностями, средствами сопряжения с другими языками и несколькими десятками пакетов приложений.

Начиная с версии Matlab 6.0, для матричных вычислений вместо функций из библиотек LINPACK и EISPACK используются подпрограммы из пакета LAPACK. На некоторых классах задач это увеличило производительность в 2–3 раза.

В настоящем пособии мы, разумеется, описываем далеко не все возможности, предоставляемые системой Matlab для решения задач линейной алгебры. Для рассматриваемых далее функций не всегда приводятся все варианты их вызова. Подробную информацию всегда можно найти в фирменной доКраткое введение в MATLAB кументации [23]. Кроме того, по системе Matlab существует много хороших учебников, например, [11, 14] и др.

После запуска системы Matlab на экране появляется основное окно, содержащее несколько подокон. Одно из них имеет заголовок Command window это командное окно, в котором пользователь набирает команды, а Matlab выдает результаты. Результаты выполнения команд, содержащие графический вывод, выдаются в одно или несколько графических окон. Команды пользователя нужно вводить после приглашения системы, которое выглядит следующим образом:

Например, (ввод заканчивается клавишей Enter). Matlab выдаст ответ:

Теперь наберем:

(1 + sqrt(5))/ получим:

В этом примере мы использовали функцию sqrt для нахождения квадратного корня; ans это специальная переменная, в которую всегда засылается результат последней команды. Эту переменную можно использовать в следующей команде. Например, Пользователь может создавать свои переменные. Например, команда создает переменную с именем e и значением 2.7181. Теперь переменную e можно использовать. Например, Получим:

Функция exp вычисляет экспоненту ex. Запись 2.2627e это представление числа в форме с плавающей точкой. Его нужно понимать следующим образом:

Перед тем как использовать переменную, ее нужно инициализировать, т. е. присвоить ей некоторое значение (так и было в двух предыдущих примерах с переменными e и err ). Использовать неинициализированные переменные запрещено. Например, если p и/или q ранее не встречались в левой части присваивания, то следующая команда 6.1. Краткое введение в MATLAB приведет к сообщению об ошибке.

До сих пор знаком окончания команды являлся символ конца строки: ввод команды заканчивался клавишей Enter. В результате мы всегда получали эхо (отклик). В конце команды, перед тем как ввести Enter, можно поставить знак ; (точка с запятой). В этом случае отклика системы мы не получим.

Например, после того как мы введем e = 2 + 1/2 + 1/6 + 1/24 + 1/120 + 1/720;

переменной e будет присвоено вычисленное значение 2.7181 и мы сразу получим новое приглашение:

Далее в примерах знак приглашения мы печатать не будем.

В одной строке можно набирать несколько команд. Их нужно отделять либо символом, (запятая), либо символом ;

(точка с запятой). В первом случае отклик будет, во втором нет.

Именами переменных могут быть любые последовательности латинских букв (в любом регистре), цифр и знаков подчеркивания. Первым символом в имени может быть либо буква, либо символ подчеркивания. Ограничений на длину имени нет, но при сравнении имен роль играют только первые символы (31 символ). Кроме того, важен регистр: например, переменные Err и err разные.

На машинах, поддерживающих IEEE-арифметику, в системе Matlab под действительный скаляр отводится число двойной точности с плавающей точкой.

Кроме ans в системе Matlab есть другие встроенные переменные. Перечислим некоторые из них.

Переменная eps хранит расстояние между 1 и следующим за ним действительным числом, представимым в компьютере.



Pages:   || 2 |
 
Похожие работы:

«Федеральное агентство по образованию Сочинский государственный университет туризма и курортного дела Филиал Сочинского государственного университета туризма и курортного дела в г.Н.Новгород ТУВАТОВА В.Е. Маркетинг гостиниц Учебно-методическое пособие для студентов всех форм обучения Нижний Новгород 2009 1 ББК 65.432я73 Т 81 Туватова В. Е. Маркетинг гостиниц: учебно-методическое пособие для студентов всех форм обучения. – Н. Новгород: типография., 2009. с. В учебно-методическом пособии...»

«Федеральное государственное бюджетное образовательное учреждение высшего профессионального образования Саратовский государственный университет имени Н.Г. Чернышевского Усанов Д.А., Скрипаль А.В., Феклистов В.Б., Вениг С.Б. ИЗМЕРЕНИЕ ПАРАМЕТРОВ ПОЛУПРОВОДНИКОВ, МИКРОИ НАНОСТРУКТУР НА СВЧ САРАТОВ 2012 УДК 620.179.18 Усанов Д.А., Скрипаль А.В., Феклистов В.Б., Вениг С.Б. У74 Измерение параметров полупроводников, микро- и наноструктур на СВЧ (учебное пособие)– Саратов: Электронное издание Сарат....»

«МИНИСТЕРСТВО ОБРАЗОВАНИЯ РЕСПУБЛИКИ БЕЛАРУСЬ УЧРЕЖДЕНИЕ ОБРАЗОВАНИЯ ВИТЕБСКИЙ ГОСУДАРСТВЕННЫЙ ТЕХНОЛОГИЧЕСКИЙ УНИВЕРСИТЕТ ТЕХНОЛОГИЯ И ОБОРУДОВАНИЕ ТКАЦКОГО ПРОИЗВОДСТВА методические указания по преддипломной практике для студентов специализации 1-50 01 01 04 Технология тканей заочной формы обучения Витебск 2006 ОГЛАВЛЕНИЕ стр. Цели и задачи практики 1. Рабочее место и баланс времени студента 2. Методические указания 3. Содержание практики 4. Ознакомление с предприятием 4.1. Склад пряжи 4.2....»

«КАРЕЛЬСКИЙ НАУЧНЫЙ ЦЕНТР РОССИЙСКОЙ АКАДЕМИИ НАУК ИНСТИТУТ ЛЕСА Б. В. Раевский, А. А. Мордась СЕЛЕКЦИОННО-ГЕНЕТИЧЕСКАЯ ОЦЕНКА КЛОНОВ СОСНЫ ОБЫКНОВЕННОЙ НА ЛЕСОСЕМЕННЫХ ПЛАНТАЦИЯХ ПЕРВОГО ПОРЯДКА Учебно-методическое пособие Петрозаводск 2006 УДК 630*165: 582.475.4 СЕЛЕКЦИОННО-ГЕНЕТИЧЕСКАЯ ОЦЕНКА КЛОНОВ СОСНЫ ОБЫКНОВЕННОЙ НА ЛЕСОСЕМЕННЫХ ПЛАНТАЦИЯХ ПЕРВОГО ПОРЯДКА / Б. В. Раевский, А. А. Мордась. Петрозаводск: Карельский научный центр РАН, 2006. 91 с.: ил. 16, табл. 24. Библиогр. 41 назв....»

«МЕТОДИЧЕСКИЕ РЕКОМЕНДАЦИИ Диспансерное наблюдение, лечение и профилактика вирусных гепатитов у подростков и взрослых больных ВИЧинфекцией Москва - 2007 Данные рекомендации разработаны Федеральным научно-методическим центром по профилактике и борьбе со СПИДом с учетом протоколов ВОЗ для европейского региона Гепатит С и ВИЧ-инфекция: тактика ведения пациентов с сочетанной инфекцией, Гепатит В и ВИЧ-инфекция: тактика ведения пациентов с сочетанной инфекцией, Профилактика гепатитов A, B, C и...»

«Министерство образования и науки Российской Федерации Федеральное государственное автономное образовательное учреждение высшего профессионального образования Уральский федеральный университет имени первого Президента России Б.Н. Ельцина Институт государственного управления и предпринимательства МЕТОДИЧЕСКИЕ РЕКОМЕНДАЦИИ ПО НАПИСАНИЮ И ОФОРМЛЕНИЮ КОНТРОЛЬНЫХ И КУРСОВЫХ РАБОТ, ВЫПУСКНОЙ КВАЛИФИКАЦИОННОЙ РАБОТЫ БАКАЛАВРА, ДИПЛОМНОГО ПРОЕКТА СПЕЦИАЛИСТА, МАГИСТЕРСКОЙ ДИССЕРТАЦИИ Екатеринбург 2012...»

«САРАНСКИЙ КООПЕРАТИВНЫЙ ИНСТИТУТ АВТОНОМНОЙ НЕКОММЕРЧЕСКОЙ ОРГАНИЗАЦИИ ВЫСШЕГО ПРОФЕССИОНАЛЬНОГО ОБРАЗОВАНИЯ ЦЕНТРОСОЮЗА РОССИЙСКОЙ ФЕДЕРАЦИИ РОССИЙСКИЙ УНИВЕРСИТЕТ КООПЕРАЦИИ АРБИТРАЖНЫЙ ПРОЦЕСС Методические указания по организации самостоятельной работы САРАНСК 2008 УДК 347.918(076) Составитель В.А. Круглов Арбитражный процесс : метод. указания по организации самостоят. работы / [сост. В.А. Круглов] ; Саран. кооп. ин-т РУК. – Саранск, 2008. – 20 с. Методические указания содержат перечень тем...»

«Станислав Макиевский ШКОЛА 11 20 постановки рук барабанщика • рг бу ер ет Учебно-методическое П пособие тнк Са • р то зи по ом К о тв ьс ел ат зд И Издательство Композитор • Санкт-Петербург • ББК 85.315. М • рг Макиевский, С. А. бу ШКОЛА постановки рук барабанщика. — СПб. : Композитор • СанктМ Петербург, 2011. — 104 с., ил., нот. ер ISMN 979-0-66000-508- ет П Книга Станислава Макиевского Школа постановки рук барабанщика является прекрасным тучебным материалом для тех, кто решил посвятить себя...»

«База нормативной документации: www.complexdoc.ru МИНИСТЕРСТВО РОССИЙСКОЙ ФЕДЕРАЦИИ ПО ДЕЛАМ ГРАЖДАНСКОЙ ОБОРОНЫ, ЧРЕЗВЫЧАЙНЫМ СИТУАЦИЯМ И ЛИКВИДАЦИИ ПОСЛЕДСТВИЙ СТИХИЙНЫХ БЕДСТВИЙ ФЕДЕРАЛЬНОЕ ГОСУДАРСТВЕННОЕ УЧРЕЖДЕНИЕ ВСЕРОССИЙСКИЙ ОРДЕНА ЗНАК ПОЧЕТА НАУЧНОИССЛЕДОВАТЕЛЬСКИЙ ИНСТИТУТ ПРОТИВОПОЖАРНОЙ ОБОРОНЫ (ФГУ ВНИИПО МЧС РОССИИ) Л.М. Мешман, С.Г. Цариченко, В.А. Былинкин, В.В. Алешин, Р.Ю. Губин ПРОЕКТИРОВАНИЕ ВОДЯНЫХ И ПЕННЫХ АВТОМАТИЧЕСКИХ УСТАНОВОК ПОЖАРОТУШЕНИЯ Учебно-методическое пособие...»

«БЕЛГОРОДСКИЙ ГОСУДАРСТВЕННЫЙ НАЦИОНАЛЬНЫЙ ИССЛЕДОВАТЕЛЬСКИЙ УНИВЕРСИТЕТ Факультет управления и предпринимательства Кафедра менеджмента организации А.Н. Алимов, Е.В. Качурова СТРАТЕГИЧЕСКИЙ МЕНЕДЖМЕНТ Методические рекомендации по выполнению курсовой работы Белгород 2013 2 УДК 005.21(075.8) ББК 95.291.213.я73 А 50 Рекомендовано учебно-методической комиссией факультета управления и предпринимательства Белгородского государственного национального исследовательского университета в качестве учебного...»

«Министерство образования Российской Федерации Новгородский государственный университет имени Ярослава Мудрого _ _ М.П. ПОПОВА Великий Новгород 2002 Министерство образования Российской Федерации Новгородский государственный университет имени Ярослава Мудрого _ _ М.П. Попова Теория и методика воспитания Учебно-методическое пособие Великий Новгород 2002 1 ББК 74.00 Печатается по решению П 58 РИС НовГУ Рецензенты: Ключникова Г.А., канд. психол. наук, доцент; Никитина Н.И., канд. пед. наук, доцент...»

«1 Олимпиада для студентов и выпускников вузов – 2013 г. Демонстрационный вариант и методические рекомендации по направлению Международные отношения Профиль Международные отношения: европейские и азиатские исследования ДЕМОНСТРАЦИОННЫЙ ВАРИАНТ Время выполнения задания – 150 мин. I. ОБЩАЯ ЧАСТЬ Прочтите текст, изложите основные идеи автора и дайте их оценку (на русском языке). Purpose and Analysis in Political Science Var. I The science of international politics has, then, come into being in...»

«В. С. Березовский, И. В. Стеценко Создание электронных учебных ресурсов и онлайновое обучение Киев Издательская группа BHV 2013 УДК 37.091.64:004 ББК 74.202.4 Б48 Березовский В. С., Стеценко И. В. Б48 Создание электронных учебных ресурсов и онлайновое обучение: [Учебн. пособ.] / В. С. Березовский, И. В. Стеценко. — К.: Изд. группа BHV, 2013. — 176 с.: ил. ISBN 978-966-552-266-9 Изложены основные принципы разработки и создания учебного контента с помощью Adobe Captivate 6, а также организации и...»

«Периферийные устройства вычислительной техники Методические указания и контрольные задания для студентов-заочников 2013 Введение Предлагаемые программой разделы учебной дисциплины позволят студентам изучить: организацию системы ввода – вывода информации, классификацию периферийных устройств; аппаратную и программную поддержку работы периферийных устройств: контроллеры, адаптеры, мосты, прямой доступ к памяти, приостановки, прерывания, драйверы; современные и перспективные интерфейсы и шины...»

«Министерство образования и науки Российской Федерации Федеральное государственное автономное образовательное учреждение высшего профессионального образования Северо-Кавказский федеральный университет ПРИНЯТО Ученым советом СКФУ Протокол № 6 от 28 февраля 2013 г. ПОЛОЖЕНИЕ ОБ УЧЕБНО-МЕТОДИЧЕСКОМ КОМПЛЕКСЕ ДИСЦИПЛИНЫ ФЕДЕРАЛЬНОГО ГОСУДАРСТВЕННОГО АВТОНОМНОГО ОБРАЗОВАТЕЛЬНОГО УЧРЕЖДЕНИЯ ВЫСШЕГО ПРОФЕССИОНАЛЬНОГО ОБРАЗОВАНИЯ СЕВЕРО-КАВКАЗСКИЙ ФЕДЕРАЛЬНЫЙ УНИВЕРСИТЕТ Ставрополь, 2013 СОДЕРЖАНИЕ 1....»

«МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ РОССИЙСКОЙ ФЕДЕРАЦИИ Дальневосточный федеральный университет Инженерная школа Л.С. Иванова, В.И. Кондратьева АНАЛИЗ ФОРМИРОВАНИЯ ИЗДЕРЖЕК И ПРИБЫЛИ ПРЕДПРИЯТИЯ Методические указания к практическим занятиям Режим доступа: http://dvfu.ru/web/is/publikacii1 Формат pdf Объем 1 MB © Иванова Л.С., Кондратьева В.И., 2012 © Дальневосточный федеральный университет © Издательский дом Дальневосточного федерального университета, оформление, Владивосток Издательский дом...»

«Банковский Ю.В. ОСНОВЫ СПОРТИВНОЙ ТРЕНИРОВКИ В ГОРНЫХ ВИДАХ СПОРТА (Альпинизм, спортивное скалолазание, горный туризм) Москва, 1996 ББК 75.82 Б183 Методическое пособие подготовлено Мастером спорта по альпинизму, Заслуженным тренером Республики Таджикистан, Доцентом БАЙКОВСКИМ Ю.В. Рецензенты: Главный тренер сборной команды России по альпинизму, Заслуженный тренер России и СССР, Мастер Спорта, Снежный Барс ШАТАЕВ В.Н. Зав.кафедрой научных основ управления школой ФППК ОНО ДГПИ, профессор РЕПА...»

«МИНОБРНАУКИ РОССИИ Федеральное государственное бюджетное образовательное учреждение высшего профессионального образования Юго-Западный государственный университет (ЮЗГУ) Кафедра бухгалтерского учета, анализа и аудита УТВЕРЖДАЮ: Первый проректор – проректор по учебной работе Е.А. Кудряшов _ 2011 г. ПРЕДДИПЛОМНАЯ ПРОИЗВОДСТВЕННАЯ ПРАКТИКА Методические указания для студентов (слушателей), обучающихся по специальности 080109.65 Бухгалтерский учет, анализ и аудит Курск 201 УДК Составители: В.В....»

«К. В. Григорьева Методические указания Часть 1. Бескоалиционные игры в нормальной форме. Факультет ПМ-ПУ СПбГУ 2007 г. 1 ОГЛАВЛЕНИЕ ЗАНЯТИЕ №1 1.1. СОДЕРЖАНИЕ ТЕОРИИ ИГР 1.2. КЛАССИФИКАЦИЯ ИГР 1.3. ИГРА В НОРМАЛЬНОЙ ФОРМЕ 1.4. РАВНОВЕСИЕ ПО НЭШУ 1.5. ОПТИМАЛЬНОСТЬ ПО ПАРЕТО. ЗАНЯТИЕ №2 2.1. АНТАГОНИСТИЧЕСКИЕ ИГРЫ. СЕДЛОВАЯ ТОЧКА 2.2. ПРИНЦИП МАКСМИНА И МИНИМАКСА ЗАНЯТИЕ №3 3.1. СМЕШАННЫЕ СТРАТЕГИИ МАТРИЧНЫХ ИГР 3.2. СИТУАЦИЯ РАВНОВЕСИЯ В СМЕШАННЫХ СТРАТЕГИЯХ 3.3. СВОЙСТВА ОПТИМАЛЬНЫХ СМЕШАННЫХ...»

«УДК 658; 65.01(075.8) ББК -80*65.2/4–65.9я73 У 91 МИНОБРНАУКИ РОССИИ ФЕДЕРАЛЬНОЕ ГОСУДАРСТВЕННОЕ БЮДЖЕТНОЕ ОБРАЗОВАТЕЛЬНОЕ УЧРЕЖДЕНИЕ ВЫСШЕГО ПРОФЕССИОНАЛЬНОГО ОБРАЗОВАНИЯ ПОВОЛЖСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ СЕРВИСА (ФГБОУ ВПО ПВГУС) Рецензент к.э.н., доц. Маркова О. В. Кафедра Менеджмент УЧЕБНО-МЕТОДИЧЕСКОЕ ПОСОБИЕ Учебно-методическое пособие по дисциплине Организация У 91 производства и менеджмент / сост. А. Е. Бугаев, П. С. Мутaлапо дисциплине Организация производства и менеджмент пова. –...»






 
© 2013 www.diss.seluk.ru - «Бесплатная электронная библиотека - Авторефераты, Диссертации, Монографии, Методички, учебные программы»

Материалы этого сайта размещены для ознакомления, все права принадлежат их авторам.
Если Вы не согласны с тем, что Ваш материал размещён на этом сайте, пожалуйста, напишите нам, мы в течении 1-2 рабочих дней удалим его.