bigpo.ru
добавить свой файл
1 2 3 4
Содержание


Вступление…………………………………………………...3

Раздел 1. Теоретические основы метода отражений (Хаусхолдера) для решения СЛАУ………………...5

1.1. Суть метода отражений для решения СЛАУ…..5

1.2. Модификации метода Хаусхолдера……………..10

1.3. Рамки использования метода отражений для решения СЛАУ………………………..…………18

Раздел 2. Программная реализация метода отражений

на языке программирования Pascal………...……24

2.1. Программный код решения СЛАУ методом отражений………………………………………24

2.2. Документация программы и ее характеристика………………………………...29

Раздел 3. Заключительные положения работы………..30

3.1. Сравнение метода отражений с остальными методами решения СЛАУ……………………...30

3.2. Перспективы программы и возможность ее модификации…………………………………….37

Заключение…………………………………………………39

Список использованной литературы…………………...41

Приложения………………………………………………...48


Вступление


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

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

Несомненно, что каждый метод имеет свои преимущества и свои недостатки, например, мало операций обеспечивают скорость работы программы и получения результата, но проигрыш в точности этого самого результата. И наоборот, если ставить перед собой цель достижение точности результата, то алгоритм решения такой задачи будет не только громоздким, но и нерациональным во временных затратах. Для анализа и определения лучшей модификации для решения СЛАУ в данной работе был выбран метод отражений (Хаусхолдера). Именно при помощи этого метода и его модификаций будет проведен сравнительный анализ алгоритмов решения СЛАУ.

Целью работы выступает описание и глубокий анализ самого метода отражений и его модификаций на примере решения нескольких видов СЛАУ.

На основании поставленной цели формируются и задачи курсового проекта:

  1. Провести теоретическое описание метода отражений и его модификаций (Раздел 1);

  2. Построить алгоритмы программы, составить программный код метода отражений и провести анализ работы этой программы (Раздел 2);

  3. Провести сравнение нескольких модификаций метода отражений (Хаусхолдера) и выбрать оптимальный метод решения СЛАУ, сформировать заключение работы;

  4. Провести сравнительный анализ всех методов решения СЛАУ и очертить границы их использования (Раздел 3).

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


Раздел 1. Теоретические основы метода отражений (Хаусхолдера) для решения СЛАУ.


1.1. Суть метода отражений для решения СЛАУ (метод Хаусхолдера).


Алгоритм Хаусхолдера приводит симметричную матрицу A размером (n x n) к трехдиагональной форме за (n - 2) ортогональных преобразований. Каждое преобразование уничтожает требуемую часть целой строки и целого соответствующего столбца. Базовой матрицей является матрица Хаусхолдера P, имеющая форму

P = 1 - 2wwT,

где w -- действительный вектор такой что |w|2 = 1. (В обозначениях, принятых здесь, матричное, или диадное, произведение векторов a и b записывается как abT, в то время как скалярное произведение этих векторов aTb). Матрица P является ортогональной, поскольку

P2 = (1 - 2wwT)(1 - 2wwT) = 1 - 4wwT + 4w(wTw)wT = 1.

Таким образом, P = P-1. Но PT = P, следовательно PT = P-1, что и доказывает ортогональность.

Перепишем P как

P = 1 - (uuT)/H,   где  H = |u|2/2,

и теперь u может быть любым вектором. Пусть x -- вектор, составленный из первого столбца матрицы A. Выберем

u = x - s|x|e1,

где e1 -- единичный вектор [1, 0,...,0]T, s -- либо 1, либо -1, а выбор знака s будет сделан позже. Тогда

Px = x - (u(x - s|x|e1)Tx)/H = x - (2u(|x|2 - s|x|x1))/(2|x|2 - 2s|x|x1) = x - u = s|x|e1.

Это показывает, что матрица Хаусхолдера P действует на заданный вектор x, уничтожая все его элементы, кроме первого. [15]

Для приведения симметричной матрицы A к трехдиагональной форме, вектор x, который будет образовывать первую матрицу Хаусхолдера, составим из нижних n - 1 элементов первого столбца. Тогда будут обнулены нижние n - 2 элемента:

 

Здесь матрицы записаны в блочной форме, причем (n-1)P1 обозначает матрицу Хаусхолдера размерами ((n-1) x (n-1)). Величина k это просто плюс или минус модуль вектора [a21, ..., an1]T.

Полная ортогональная трансформация теперь будет выглядеть как:
 

Здесь был использован факт, что PT = P.

Для второй трансформации матрица Хаусхолдера составляется на основе нижних (n - 2) элементов второго столбца, и из нее конструируется матрица преобразования:



Единичный блок в левом верхнем углу обеспечивает сохранение трехдиагональности, полученной на предыдущей стадии, в то время как матрица Хаусхолдера (n-2)P2 размера (n-2) добавляет еще один столбец и строку к трехдиагональному результату. Ясно, что цепочка (n-2) таких преобразований приведет матрицу A к трехдиагональному виду.

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

p = (Au)/H.

Тогда

AP = A(1 - (uuT)/H) = A - puT,   A' = PAP = A - puT - upT + 2KuuT,

где скаляр K определяется как:

K = (uTp)/(2H).

Если мы обозначим

q = p - Ku,

то будем иметь

A' = A - quT - uqT.

Последняя формула удобна для вычислений.

Программа редукции Хаусхолдера, основанная на алгоритме из [15] и помещенная ниже, реально начинает действие с n-ого столбца матрицы A, а не с первого, как было написано выше. В деталях алгоритм выглядит следующим образом. На стадии m (m = 1, 2, ..., n-2) вектор u имеет вид:

u = [ai1, ai2, ..., ai,i-2, ai,i-1 + , 0, ..., 0].

Здесь

i = n - m + 1 = n, n-1, ..., 3.

Величина  (что соответствует s|x| в прежних обозначениях) определяется как

2 = (ai1)2 + ... + (ai,i-1)2.

Знак  выбирается так, чтобы совпадать со знаком ai,i-1 для уменьшения ошибки округления.

Переменные вычисляются в следующем порядке: , u, H, p, K, q, A'. На каждой m-ой стадии матрица A становится трехдиагональной в последних m - 1 строках и столбцах.

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

Q = P1P2...Pn-2

на матрицу этих векторов. Произведение трансформирующих матриц Q формируется рекурсией после нахождения всех Pi:

Qn-2 = Pn-2, ..., Qj = PjQj+1, ..., Q = Q1.

На входе программы задается действительная симметричная матрица в массиве a[1...n][1...n]. На выходе этот массив содержит элементы матрицы Q. Вектор d[1...n] содержит множество диагональных элементов A', вектор e[1...n] содержит внедиагональные элементы A', начиная с e[2]. Заметим, что поскольку исходное содержимое массива a уничтожается, то при производстве серии последовательных вычислений с A он должен быть скопирован перед передачей в программу.

Промежуточные результаты не требуют дополнительной памяти. На стадии m вектора p и q не равны нулю только для элементов с индексами 1 ... i (вспомним, что i = n - m + 1), а вектор u -- только для элементов с индексами 1 ... i - 1. Элементы вектора e определяются в порядке n, n - 1, ..., так что вектор p можно держать на месте неопределенных элементов e. Вектор q может записываться поверх p, так как после определения q вектор p уже не требуется. Вектор u размещается в i-ой строке массива a, а u/H -- в i-ом столбце. После проведения редукции, матрицы Qj вычисляются с использованием записанных в массив a векторов u и u/H. Поскольку матрица Qj является единичной для последних n - j + 1 строк и столбцов, то требуется вычислить ее элементы лишь до строки и столбца n - j. При этом можно перезаписать u и u/H в соответствующих строке и столбце, поскольку они более не нужны.

Наша программа содержит еще одну тонкость. Если величина 2 равна нулю или "мала" на какой-либо из стадий, то соответствующее преобразование может быть пропущено. Простой критерий типа

2 < (наименьшее представимое число)/(машинная точность)

в большинстве случаев является слишком сильным. В действительности берется более аккуратный критерий. Определим число

 = k=1i-1(|aik|).

При  от нуля до машинной точности преобразование пропускается. В противном случае мы переопределяем

aik --> aik/

и проводим преобразование с масштабированными величинами. (Преобразование Хаусхолдера зависит только от отношения элементов).

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


1.2. Модификации метода Хаусхолдера.





Интервальный метод Хаусхолдера для комплексных линейных систем.

Рассматривается метод Хаусхолдера для внешнего оценивания мно­жеств решений систем линейных алгебраических уравнений с комплексными интер­вальными параметрами. Для этого используется нормализованное приведение мат­рицы системы к верхней "треугольной" форме с помощью последовательности орто­гональных преобразований.

Введение

Исследуется интервальная система линейных алгебраических уравнений (ИСЛАУ)



где элементами nхn-матрицы A и n-вектора правой части b являются ком­плексные интервалы.

Множество решений ИСЛАУ в общем случае не имеет простого описания. Поэтому мы ограничимся нахождением интервального вектора, целиком содер­жащего это множество решений. Как известно, полученный вектор называется внешней оценкой множества решений ИСЛАУ.

Одним из способов нахождения внешней оценки множества решений ИС­ЛАУ является метод Хаусхолдера [1], его вещественный вариант известен так­же в отечественной литературе как метод отражений [2, 3]. Суть этого метода заключается в том, что матрица системы A приводится к верхней треуголь­ной форме с помощью последовательности ортогональных преобразований. В отличие от метода Гаусса решения ИСЛАУ [4], метод Хаусхолдера дает более узкие внешние оценки множеств решений, если интервалы в матрице и правой части "не слишком широки". Кроме того, интервальный метод Хаусхолдера поз­воляет получать ответ в тех задачах, к которым интервальный метод Гаусса неприменим.

Отметим, что метод Хаусхолдера ранее применялся для внешнего оцени­вания множеств решений ИСЛАУ с вещественными интервальными парамет­рами. Целью работы является обобщение метода Хаусхолдера на случай ком­плексных ИСЛАУ.

Чтобы сделать изложение нашей работы самодостаточным, мы подробно описываем множество комплексных интервалов и основные свойства арифме­тик комплексных интервалов.

Обозначения:

  • – множество вещественных чисел,

– множество комплексных чисел,

  • – множество замкнутых интервалов на ,

  • — множество n-мерных векторов с элементами из

– множество -матриц с элементами из

  1. Комплексные интервальные арифметики

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

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

1.1.1. Прямоугольники в качестве комплексных интервалов Определение 1. Множество



называется прямоугольным комплексным интервалом.

Множество всех таких комплексных интервалов обозначим через

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

где

а каждый элемент — как сумму



откуда видно, что .

Элементы и из считаются равными, если



Определенное здесь отношение равенства рефлексивно, симметрично и транзитивно.

Теперь обобщим вещественную интервальную арифметику на комплексный случай.

Определение 2. Пусть - арифметическая операция над элементами из Тогда если



то мы полагаем



где



  • определяется аналогично.

Считается, что в случае деления Это требование выполняется при условии . иначе деление не определено.

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

1.2.2. Круги в качестве комплексных интервалов

Определение 3. Пусть — комплексное число и Множество



называется круговым комплексным интервалом или кругом.

Множество всех кругов обозначим через . Круг z с центром a и радиусом r будем записывать в виде .

Комплексные числа можно рассматривать как специальные элементы из , имеющие вид . Ясно, что .

Два круга и считаются равными, если и .

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

Определение 4. Пусть - арифметическая операция над комплексными числами. Тогда если и , то



при условии, что ,

при условии, что .

Здесь обозначает модуль комплексного числа а - сопряженное к

Соберем теперь вместе наиболее важные свойства операций на и Если не оговорено противное, то можно понимать и как обозначение множества с операциями из определения 2, и как обозначение множества с операциями из определения 4.

Теорема 1. Пусть Тогда

I. (коммутативность).

II. для из . (ассоциативность) .

III. для из

IV. (соответственно ) и (соответственно ) ) — определенные единственным образом нейтральные элементы сложения (нуль) и умножения (единица).

V. не имеет делителей нуля.

VI. Элемент z множества имеет противоположный и обратный элементы, только если и в случае обратного элемента Однако

VII. (субдистрибутивность), дляa [6].

Особо отметим, что ассоциативный закон в общем случае не выполняется при перемножении элементов



      1. Метод Хаусхолдера

Дана система интервальных уравнений

(1)

с невырожденной матрицей и вектором

Напомним, что интервальная матрица A называется невырожденной, если невырождены все матрицы A из A.

Множество решений интервальной системы уравнений (1) определяется в виде:



Как было уже отмечено, множество решений ИСЛАУ в общем случае не име­ет простого описания, поэтому мы ограничимся нахождением интервального вектора , целиком содержащего это множество решений.

Пусть интервальный вектор из . Интервальное расширение нормы вектора , обозначаемое как

, берется следующим образом:



Где — сопряженное с .

Если , то обозначим через интервальный вектор, полученный нормированием интервального вектора :



Где - знак середины действительной части комплексного интервала для всех

Определим интервальную -матрицу H как

(2)

Н является интервальным расширением эрмитовой унитарной матрицы Ха­усхолдера, порожденный интервальным вектором v. Кроме того, H содержит эрмитовы и унитарные матрицы отражений H, порожденные комплексными векторами . В дальнейшем, интервальное расширение нормы вектора и интервальное расширение матрицы Хаусхолдера мы будем называть просто нормой интервального вектора и интервальной матрицей Хаусхолдера соот­ветственно. [7,8]

Пусть где — вектор-столбцы интервальной матрицы. Мы собираемся вычислить интервальную матрицу Q и верхнюю треугольную интервальную nхn-матрицу R с помощью последова­тельности трансформаций Хаусхолдера. Здесь следует уточнить, что мы полу­чим не строго треугольную интервальную матрицу R, так как из-за некоторых особенностей классической интервальной арифметики, мы не можем получить 0 при выполнении арифметических операций с интервалами ненулевой ширины, но обязательным условием является включение нуля.

Обозначим через действительный p-вектор, первый компонент которого равен 1, а остальные - 0. Предположим, что = причем условимся, что, если вещественный интервал является уравновешанным (т.е. ), то

Тогда H — интервальная матрица Хаусхолдера, порождаемая вектором

Мы имеем



Где для , .

Первая строка матрицы R:



Положим Мы получим интервальную матрицу .



Обозначим элементы как

Предположим, что и пусть

Тогда H — интервальная -матрица Хаусхолдера, порождаемая вектором

Имеем



Где для

Вторая строка матрицы R:



Положим



Далее, повторяем процедуру построения матриц R и Q, в данном случае



Интервальный метод Хаусхолдера позволяет вычислить внешнюю оценку x множества решений решением верхней треугольной интервальной системы.

где

Данный вывод следует из того, что множество решений В следствии монотонности включения матриц

Далее, воспользуемся формулами обратной подстановки для нахождения ком­понентов интервального вектора x:



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


1.3. Рамки использования метода отражений для решения СЛАУ.


1.3.1. Численные эксперименты интервального метода для комплексных СЛАУ.

Рассмотрим пример применения метода Хаусхолдера для нахождения внеш­ней оценки множества решений ИСЛАУ (1) с матрицей A и правой частью b [4]:



Одним из условий существования решения является невырожденность ин­тервальной матрицы A:

,

В данном случае, следовательно, матрица A – невырождена.

Далее, определяем интервальные векторы



Следовательно,



Вычислим векторv, порождающий интервальную матрицу Хаусхолдера:



Определим вектор w нормированием вектора:



Далее, вычисляем интервальную матрицу Хаусхолдера по формуле (2):



Построим первую строку интервальной матрицы R:



В итоге получим матрицу



cледовательно,



Вычислим последний диагональный элемент матрицы R:



На этом построение матриц R и Q закончено, причем мы получили матрицу Q = H. [9]

Теперь нам необходимо определить интервальный вектор c:



Далее, воспользуемся формулами обратной подстановки Гаусса для нахожде­ния компонентов интервального вектора x:



В итоге мы получим внешнюю оценку



целиком содержащую множество решений нений.

1.3.2. Преобразование Хаусхолдера с использованием QR - разложение матрицы.

Используя преобразование Хаусхолдера, построить QR - разложение матрицы



Решение.

1. Положим A0 = A и найдем ортогональную матрицу Хаусхолдера H1, такую что в

матрице A1 = H1A0 все поддиагональные элементы первого столбца равны нулю. С этой целью компоненты вектора v определим, используя элементы первого столбца матрицы

A0:



В результате получен вектор

Найдем соответствующую этому вектору матрицу Хаусхолдера:



В заключение первого шага вычислим матрицу A1:



Таким образом, после первого шага получена матрица с нулевыми поддиагональными элементами в первом столбце. [10]

2. На втором шаге проделаем аналогичную процедуру, обнуляя поддиагональный элемент второго столбца.



Т.е. искомый вектор

Далее найдем соответствующую ему матрицу Хаусхолдера:



и вычислим матрицу A2:




Таким образом, исходная матрица А приведена к верхнему треугольному виду, т.е. получена матрица R = A2 искомого разложения.

Результирующая ортогональная матрица преобразования Q получается в результате перемножения матриц Hi, i = 1,2:



В заключение выпишем окончательный результат А QR= в явном виде:







1.3.3. Применение метода отражений с использованием алгоритма QR и указанной точностью.

С помощью QR - алгоритма вычислить собственные значения матрицы A из предыдущего примера с точностью= 0,01 ,

Решение.

1. Положим А(0) = А и найдем QR - разложение этой матрицы

Эта процедура подробно рассмотрена в предыдущем примере. Получены следующие



Матрицу A(1) определим перемножением полученных в



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

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


2. а). Находим QR - разложение аналогично примеру (используя преобразование Хаусхолдера):




б). Перемножая полученные выше матрицы в обратном порядке находим матрицу



Продолжая итерационный процесс, получим соответственно на 6-ой и 7-ой итерациях следующие матрицы:



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

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

окончательное решение задачи можно записать в виде:






следующая страница >>