Алгебраїчна проблема власних значень

Размещено на http://www.allbest.ru/

Міністерство освіти і науки, молоді та спорту України

Національний університет «Львівська Політехніка»

Факультет прикладної математики

Курсова робота

З курсу чисельні методи

На тему: «Алгебраїчна проблема власних значень»

Львів 2014

Зміст

Вступ

1. Алгебраїчна проблема власних значень

1.1 Метод степенів

1.2 QR-метод

2. Чисельні експерименти

Висновок

Список літератури

Додаток

Вступ

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

, (*)

де -задана комплексна (дійсна чи недійсна) матриця розміру . Розв'язки рівняння (*) є коренями характеристичного рівняння . Ліва частина цього рівння є поліномом степеня по , і, відповідно, характеристичне рівняння має рівно коренів з урахуванням їх кратності. Якщо власне значення знайдено, то відповідний власний вектор можна, в принципі, виначити як розв'язок однорідної системи рівнянь

. (**)

Підкреслимо, що навіть якщо матриця дійсна, її власні значення (всі чи деякі), а відповідно, і власні вектори, можуть бути недійсними.

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

За означенням (за допомогою характеристичного рівняння) можна знаходити тільки власні значення матриць розмірності менше п'яти. Характеристичне рівняння має степінь рівну степені матриці.

1. Алгебраїчна проблема власних значень

Алгебраїчна проблема власних значень формулюється таким чином: знайти числа та вектори , для яких

(1)

де -- задана матриця з множини матриць з комплексними елементами;

-- множина комплексних чисел. Числа називаються власними числами (значеннями), а відповідні вектори -- правими власними векторами матриці .

Множина

утворює підпростір векторів простору , і цей підпростір має розмірність

Число є тоді і лише тоді власним числом матриці , коли , тобто коли >, і

Многочлен

називається характеристичним многочленом матриці А, і його корені є власними значеннями А. Якщо є різними коренями , то

Число називається кратністю власного значення, точніше алгебраїчною кратністю.

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

Приклад 1. Діагональна матриця

має характеристичний многочлен і є єдиним власним значенням, якому відповідає власний вектор . У даному випадку

Приклад 2. Матриця

має характеристичний многочлен і є єдиним власним значенням з алгебраїчною кратністю . Але . Отже, геометрична кратність

-- перший координатний вектор.

В силу рівностей

маємо: якщо є власним значенням матриці , то є власним значенням і матриці , а є власним значенням . Із співвідношень

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

де зберігає власні значення, а відповідні власні вектори змінює за формулою . Не змінюється при цьому характеристичний поліном, а також числа та : для це випливає з інваріантності характеристичного многочлена, а для -- з того, що для вектори є лінійно незалежними тоді і тільки тоді, коли лінійно незалежними є На еквівалентних перетвореннях грунтується багато чисельних методів розв'язування проблем власних значень і власних векторів.

Наведемо без доведення теореми про деякі форми представлення матриць.

Теорема 1. Нехай А --довільна -матриця і -- її різні власні значення з геометричними та алгебраїчними кратностями. Тоді для кожного власного значення існують натуральних чисел , таких, що

а також невироджена-матриця Т така, що матриця , яка називається нормальною формою Жордана матриці А, має вигляд



називаються жордановими клітинками. При цьому числа (а з ними і матриця J), з точністю до перестановки однозначно визначені, а матриця Т, взагалі кажучи, визначається неоднозначно.

Розщепимо матрицю за стовпчиками відповідно до нормальної форми Жордана:

Тоді із співвідношення маємо

(2)

Позначимо стовпчики матриці через тобто

Тоді із (2) та з означення матриць дістанемо

Зокрема помічаємо, що (перший стовпчик матриці ) є власним вектором матриці А, що відповідає власному значенню . Інші вектори називаються головними векторами, відповідними значенню і ми бачимо, що кожній клітинці Жордана відповідає один власний вектор і набір головних векторів. В цілому для кожної матриці А можна знайти базис простору (а саме стовпчики матриці ), який складається з власних та головних векторів матриці А.

Характеристичні поліноми окремих клітинок Жордана називаються елементарними дільниками матриці А. Таким чином, матриця А має лінійні елементарні дільники тоді і лише тоді, коли тобто жорданова нормальна форма матриці є діагональною матрицею. У цьому разі називається матрицею, що зводиться до діагонального вигляду, або матрицею, що припускає нормалізацію, тобто в є базис, який складається лише з власних векторів матриці , а головні вектори не з'являються. Отже, кожну матрицю , яка має різні власні значення, можна звести до діагнонального вигляду за допомогою перетворення подібності.

Далі розглянемо деякі класи матриць, що зводяться до діагонального вигляду. Власні вектори таких матриць утворюють базис в .

Якщо в перетворенні подібності припускати, що не довільні несингулярні матриці, то у загальному випадку не може бути зведена до форми Жордана. Для унітарних матриць Т, тобто має місце така теорема.

Теорема 2 (теорема Шура). Для довільної матриці існує унітарна матриця U така, що

де -- власні значення матриці А (не обов'язково різні).

Якщо А = А*, тобто є ермітовою матрицею, то (U*AU)* = U*A*U** = U*AU є також ермітовою і з теореми 4 випливає така теорема.

Теорема 3. Для довільної ермітової матриці А існує унітарна матриця U така, що

При цьому власні значення матриці А є дійсними, а і-й стовпчик матриці U є власним вектором, відповідним тобто А має п лінійно незалежних ортогональних власних векторів.

Узагальненням ермітових є нормальні матриці, для яких А*А = АА*.

Теорема 4. Матриця є нормальною (А*А = А А*) тоді і лише тоді, коли існує унітарна матриця U така, що

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

Для довільної прямокутної матриці A матриця А* А є невід'ємно визначеною з власними значеннями . Числа називаються сингулярними числами матриці А.

Теорема 5. Нехай . Тоді:

1) існує унітарна матриця U та унітарна матриця V такі, що матриця U*AV = є діагональною матрицею вигляду

>

де -- відмінні від нуля сингулярні числа матриці A, r -- ранг матриці А;

2) відмінними від нуля сингулярними числами матриці А* є також . Розвинення називається сингулярнозначним зображенням матриці А (розвинення за сингулярними значеннями).

Унітарні матриці U та V можна інтерпретувати таким чином: стовпчики матриці U являють собою m ортонормальних власних векторів ермітової матриці АА*, а стовпчики V є ортонормальними власними векторами ермітової матриці А*А, бо U*AA*U = , V*A*AV = . Діагональна матриця

є псевдооберненою до , або оберненою матрицею Мура -- Пенроуза, тобто матрицею, що задовольняє співвідношення: а; б)*; в) . Тому, як легко помітити, матриця

є в силу єдиності псевдооберненою до А.

Нехай власний вектор є спряженим до

.

Заради простоти вважатимемо, що -- просте власне значення, тобто воно є простим коренем многочлена .

Лема 1. Нехай в простими власними значеннями матриць A, , що відповідають власним векторам та . Тоді існує неперервно диференційовне відображення



деякого околу V з простору матриць з центром в А таке, що і є простим власним значенням матриці В для всіх причому

Д о в е д е н н я. Оскільки є простим власним значенням характеристичного многочлена , то

.

З теореми про неявну функцію випливає, що в деякому околі нуля існує неперервно диференційовне відображення

таке, що і є простим власним значенням матриці . Крім того, існує неперервно диференційовна функція

Така, що і х (t) є власним вектором матриці , який відповідає власному значенню , тобто

Помноживши цю рівність, а також рівність скалярно на , знайдемо



=

Після ділення на t і переходу до границі при дістанемо

Враховуючи, що

що і треба було довести.

Щоб оцінити обумовленість задачі маємо обчислити норму як лінійного відображення

Виберемо в просторі норму, а в просторі за норму вважатимемо модуль. Тоді

причому якщо то в першій нерівності досягається знак рівності. Це матиме місце, зокрема, для . Для матриці , очевидно, і в другій нерівності справжується знак рівності. Тому з урахуванням співвідношення маємо

де через позначено кут між власним вектором та спряженим вектором . Для нормальних матриць , для яких за означенням , та збігаються, тобто . Тому для таких матриць .

Сформулюємо здобуті результати у вигляді теореми.

Теорема 6. Абсолютне та відносне числа обумовленості задачі обчислення простого власного значення матриці стосовно матричної норми виражаються формулами

де -- власний вектор матриці А, що відповідає власному значенню , -- спряжений до нього власний вектор матриці *, тобто , . Зокрема, проблема власних значень для нормальних матриць є добре обумовленою з числом обумовленості .

Приклад матриці

з власними значеннями та збуреної матриці

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

тобто задача обчислення власного значення матриці стосовно абсолютної похибки взагалі не є коректно поставленою задачею. Інший приклад, розглянутий у вступі (п. 2), показує, що обчислення власних значень як коренів характеристичного многочлена (навіть добре розділених і однократних) хоча теоретично і можливе, але потребує великої обережності. Тому на практиці, як правило, застосовуються інші методи, деякі з яких ми розглянемо далі.

1.1. Метод степенів

Цей метод застосовується для ітераційного обчислення найбільшого за модулем власного значення і відповідного власного вектора, а також власних значень, для яких відомі вже досить хороші наближення.

Нехай -- матриця, яку можна звести до діагонального вигляду, з власними значеннями , для яких

Припустимо додатково, що не існує відмінних від власних значень, для яких , тобто існує r > 0 таке, що

Оскільки зводиться до діагонального вигляду, то має лінійно незалежних власних векторів які утворюють базис в . Нехай -- довільний вектор, а послідовність утворюється за правилом

Вектор можна подати у вигляді

і тому в силу маємо

Припустимо додатково, що для має місце , тобто має «досить загальний вигляд». Нормуючи яким-небудь способом, наприклад,



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

Маємо, що для ( -- просте власне число) граничний вектор є незалежним від вибору (якщо лише ). Якщо (--многократне власне значення), то залежить від співвідношення коефіцієнтів , тобто від початкового наближення . Швидкість збіжності оцінюється величиною і є тим більшою, чим менше це число. З викладеного вище також випливає, що метод буде збігатися не до , а до і відповідного власного вектора, якщо вибрати в розвиненні для (і якщо додатково не існує відмінних від власних значень з тим самим модулем). Проте це справедливо лише теоретично, бо внаслідок неминучих похибок заокруглення ми вже після першої ітерації для дістанемо

з деяким і метод, таким чином, збігатиметься до .

Якщо не зводиться до діагонального вигляду і має єдине власне значення з найбільшим модулем, то ми можемо подати через власні і головні вектори матриці А і аналогічно довести, що для досить загального вигляду метод степенів збігається до та відповідного власного вектора.

Практична обмеженість описаного вище методу полягає в тому, що він повільно збігається, якщо абсолютні значення власних чисел близькі, а також у тому, що за ним обчислюють лише одне (максимальне за модулем) власне значення і власний вектор. Цих недоліків, проте, можна позбутися за допомогою методу зворотних ітерацій, запропонованого в 1945 р. Віландтом. Цей метод застосовується тоді, коли для якогось із власних значень наприклад , відоме «досить хороше наближення» , тобто

Тоді для «досить довільного» початкового вектора будується послідовність

Якщо то існує і остання рівність еквівалентна

тобто цей метод є звичайним методом степенів для матриці та власних чисел , причому

Припустивши знову, що зводиться до діагонального вигляду і має власні вектори , і поклавши дістанемо для простого

Швидкість збіжності тим більша, чим меншим є відношення для тобто чим кращим є наближення . Зауважимо, що матриця для досить хороших наближень є майже виродженою, проте можна довести, що оскільки ми шукаємо лише напрям власного вектора, то ця задача є добре обумовленою і не виникає жодних труднощів. Зазначимо також, що для обчислення досить один раз обчислити -розвинення матриці і потім з його допомогою обчислювати розв'язки систем з різними правими частинами.

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

Найпростіша ідея, яка здавалось би веде до розв'язування повної проблеми власних значень, полягає в тому, щоб за допомогою перетворень подібності, скажімо, за допомогою унітарних матриць, звести матрицю до діагонального вигляду. Перша ж спроба зробити це за допомогою матриць Хаусгольдера показує, що це неможливо, бо нулі, які утворюються після множення на зліва, «руйнуються» після множення на * справа:

Проте легко помітити, що таким чином можна матрицю звести до так званої форми Хессенберга:

-матриця Хаусгольдера, побудована за обведеним тривимірним вектором. Якщо -- ермітова матриця, то замість згаданих вище ситуацій матимемо такі:



причому в силу ермітовості матриці матриця також буде ермітовою:

Сформулюємо ці міркування у вигляді такої леми.

Лема 2. Нехай . Тоді існує унітарна матриця Р, яка є добутком (п -- 2) матриць Хаусгольдера і така, що

-- форма Хессенберга для неермітових матриць та

-- тридіагональна ермітова матриця для ермітоеих матриць (, тобто -- дійсні числа).

Д о в е д е н н я. Продовжуючи описаний вище процес за допомогою матриць відображення Хаусгольдера дістанемо

1.2. QR-метод

На ідеї зведення матриці до «простішого» вигляду за допомогою перетворень подібності будується ряд чисельних методів. Проте розглянемо спочатку загальніший -метод, який ефективно використовується на практиці. Історично цей метод, запропонований Ф. Френсісом у 1959 р. і В. М. Кублановською в 1961 p., можна розглядати як розвинення -методу Рутісхаузера (1958 p.). Алгоритмічно вони дуже схожі, лише в останньому -- замість -розвинення матриць застосовується -розвинення, що є недоліком, бо -розвинення існує не для всіх матриць (навіть невироджених). -poзвинення матриці існує завжди, але з точністю до так званої фазової матриці

Дійсно, якщо А = QR, де Q -- унітарна матриця (), то разом з Q є унітарною матриця QS і .

Суть -алгоритму полягає в побудові послідовності матриць за правилом: (-розвинення), Ми бачили, що -розвинення матриць можна виконати, наприклад, за допомогою чисельно стійких перетворень Хаусгольдера

де є верхньою трикутною матрицею. Оскільки для перетворень Хаусгольдера маємо , то

Лема 3. Матриці а також матриці

мають такі властивості:

1) є подібною до тобто

2)

3)

Д о в е д е н н я. 1) Оскільки, то

2) Це твердження є елементарним наслідком 1).

3) Із властивості 2) випливає, що

Наступна теорема, доведення якої проведемо за Вілкінсоном, доводить збіжність -алгоритму.

Теорема 7. Нехай матриця має такі власні значення, що (3)

і нормальну жорданову форму в якій матриця Y має LR-розвинення

Тоді матриці визначені в QR-алгоритмі, мають такі властивості збіжності: існують фазові матриці

1)

2)

де -- елементи матриці

Д о в е д е н н я. Оскільки і існує , то

(4)

де є -розвиненням невиродженої матриці в добуток унітарної матриці та невиродженої верхньої трикутної матриці . Неважко помітити, що матриця є нижньою трикутною матрицею, причому

Оскільки для , то для , а також

Тому із (11) маємо

(5)

де

Розглянемо далі додатно визначену матрицю

в якій і скористаємося тим фактом, що будь-яка додатно визначена матриця має розвинення

де є верхньою трикутною матрицею з додатними діагональними елементами. При цьому зазначимо, що множник залежить неперервно від матриці , тому

Матриця

Це означає, що матриця має -розвинення

Тому з (12) маємо

де матриця є унітарною, а матриця -- верхньою трикутною. З іншого боку, з леми 4 випливає; що матриця має також -розвинення



Оскільки -розвинення матриць єдине з точністю до деякої фазової матриці, то існують матриці

і оскільки , то

Матриця є верхньою трикутною вигляду

Збіжність -алгоритму можна довести і при слабкіших обмеженнях, ніж в теоремі 9. Якщо, наприклад,

 (6)

то для матриць можна довести таку теорему.

Теорема 8. Нехай із умов теореми 7 замість (3) виконується (6). Тоді

а)

б)

в) хоча матриці

при не збігаються, а їхні власні значення збігаються до та тобто

де збіжність має місце для елементів, які позначені та , а також для власних значень -матриць, позначених , які збігаються до та .

Можна також довести, що коли матриця зводиться до діагонального вигляду

але Y не мае -розвинення, то -метод все одно збігається, лише в граничній матриці діагональні елементи, які є власними значеннями матриці , не обов'язково будуть впорядкованими за значеннями модуля.

Звернемо увагу на такі недоліки -алгоритму в описаній вище формі: а) якщо матриця повністю заповнена, то на кожному кроці потрібно виконати операцій; б) збіжність алгоритму невисока, якщо власні значення погано розділені, тобто .

Щоб усунути перший недолік, тобто зменшити кількість арифметичних операцій, матрицю за допомогою перетворень подібності зводять спочатку до форми Хессенберга (для матриць загального вигляду):

або до тридіагонального вигляду (для ермітових, зокрема симетричних, матриць):

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

Лема 4. Матриці мають такі властивості:

1) подібні до , тобто ;

2) якщо є симетричною, то симетричними є і матриці ;

3) якщо є симетричною тридіагональною матрицею, то такими самими є і , якщо реалізуються за допомогою поворотів Гівенса.

Д о в е д е н н я. 1) Це твердження доведене в лемі 3.

2) Це твердження є наслідком властивості 1) і того, що для

3) Нехай є симетричною тридіагональною матрицею. Реалізуємо за допомогою поворотів Гівенса так, що Позначаючи через елементи, які виключаються, а через нові ненульові елементи, дістаємо:



За властивістю 2) матриця має бути знову симетричною, тому всі елементи, які позначені , в насправді дорівнюють нулю. Отшже, є тридіагональною.

Вказаний вище недолік б) -алгоритму, а саме повільна збіжність при можна усунути за допомогою так званої стратегії зсувів (Shift-strategy). Для цього на кожному ітераційному кроці і застосовується параметр зсуву і послідовність визначається таким чином: (-розвинення), Неважко довести, що для тридіагональних симетричних матриць (див. леми 3, 4)

Послідовність збігається до зі швидкістю

Параметри мають лежати якомога ближче до . Вілкінсон запропонував такий метод вибору параметрів зсуву для симетричної тридіагональної матриці : якщо

то за вибирається те власне значення -матриці

яке лежить ближче до . Крім методу явних зсувів для прискорення збіжності застосовується метод неявних зсувів, який, проте, виходить за рамки цієї книги.

Для симетричних тридіагональних матриць потрібно виконати операцій, щоб відшукати кожне власне значення, отже операцій, щоб відшукати всі власні значення.

Крім власних значень, нас цікавлять також власні вектори, які для симетричних матриць можна знайти таким чином: якщо , то стовпчики матриці апроксимують власні вектори матриці , тобто .

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

1. Звести до тридіагонального вигляду

-- симетрична і тридіагональна матриця;

2. Знайти власні наближені значення за допомогою -алгоритму з поворотами Гівенса: добутком поворотів Гівенса

3. Стовпчики є апроксимаціями власних векторів :

Цей алгоритм потребує ~ множень для перетворення до тридіагонального вигляду і для виконання -алгоритму. Для загальних матриць спочатку виконується перетворення до форми Хессенберга, а потім за допомогою -алгоритму -- до нормальної форми Щура (комплексної верхньої трикутної матриці).

2. Чисельні експерименти

Нехай необхідно знайти власні значення та власні функції деякої матриці. Розклад цієї матриці за допомогою програмної реалізації QR-алгоритму наведеної в додатку буде мати наступний вигляд:

вектор матриця програмний алгоритм

Покладемо і знайдемо QR-розклад матриці .

,

Матрицю визначимо перемноженням отриманих в результаті QR-розкладу матриць в зворотньому порядку :

.

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

, .

Перемноживши отримані матриці в зворотньому порядку знаходимо матрицю

.

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

, .

Видно, що під діагональні елементи першого стовпця стають достатньо малими, і, відповідно, діагональний елемент може бути прийнятий власним значенням. В той же час прослідковується комплексно-спряжена пара власних значень, яка відповідає блоку, утвореному елементами другого і третього стовпців , , , . Незважаючи на значні зміни в ході ітерацій самих цих елементів, власні значення, які відповідають даному блоку і визначаються з розв'язку квадратного рівняння , змінюються несуттєво (в межах допустимої похибки). Таким чином, кінцевий роз'язок задачі можна записати у вигляді:

Висновок

Вибір найбільш ефективного методу визначення власних значень або власних векторів залежить від ряду факторів (тип рівнянь, число і характер власних значень, вид матриці і т.д.). Алгоритми визначення власних значень можна розподілити на 3 групи:

1) прямі, які ґрунтуються на розкритті характеристичних визначників і розв'язанні характеристичних рівнянь;

2) ітераційні, які ґрунтуються на багатократному застосуванні ітераційного алгоритму, який наближає власний вектор, що одержується в кожному циклі, до точного розв'язку;

3) перетворення подібності, які використовують властивості подібних матриць, що мають однакові власні значення і ортогональні власні вектори.

Список літератури

1. Гаврилюк І.П., Макаров В.Л. Методи обчислень: Підручник: У2ч.-- К. і Вища шк., 1995,

2. Дж. Ортега, У. Пул Введение в численные методы. Решения дифференциальных уравнений М.: Наука 1986.

3. Квєтний Р.Н., Богач І.В., Бойо О.Р., Софина О. Ю., Шушура О. М. Комп'ютерне моделювання систем та процесів. Методи обчислень. Частина 1. В.: ВНТУ 2013.

Додаток

#include<iostream>

#include<vector>

#include<ostream>

#include <algorithm>

#include<iomanip>

#include<math.h>

using namespace std;

class Vect_op;

typedef vector<double>Vect;

typedef vector<Vect>Matr;

typedef vector<Vect_op>Matr_op;//a matrix which contain vectors Vect_op

Matr EnterMatrix();

void PrintMatrix(Matr a);

int fact(int a);

Matr InvertMatrix(Matr a);

class Vect_op// a class for overloading the operators so that I can add multiply etc.. these vectors

{

public:

Vect_op()

{

Vect a;

for(int i = 0; i < 3; i++)

a.push_back(0);

itsVect = a;

}

Vect_op(Vect& a){ itsVect = a;}//copyctor

~Vect_op(){};

void setVect(const Vect& a){ itsVect = a;}

Vect getVect(){ return itsVect;}

Vect_op operator+(const Vect_op& a)

{

Vect b = a.itsVect;

Vect c = a.itsVect;

for(int i = 0; i < b.size(); i++)

c[i] = itsVect[i] + b[i];

Vect_op z(c);

return z;

}

Vect_op operator-(const Vect_op& a)

{

Vect b = a.itsVect;

Vect c = a.itsVect;

for(int i = 0; i < b.size(); i++)

c[i] = itsVect[i] - b[i];

Vect_op z(c);

return z;

}

double operator*(const Vect_op& a)

{

Vect b = a.itsVect;

double res = 0;

for(int i = 0; i < b.size(); i++)

res += (itsVect[i] * b[i]);

return res;

}

Vect_op operator*(const double& x)

{

Vect b = itsVect;

for(int i = 0; i < b.size(); i++)

b[i] *= x;

Vect_op w(b);

return w;

}

double operator[](int x)

{

Vect b = itsVect;

return b[x];

}

double sq_norm()

{

double res = 0;

for(int i = 0; i < itsVect.size(); i++)

res += (itsVect[i] * itsVect[i]);

return res;

}

friend

ostream& operator<<( ostream& xout, const Vect_op& v)

{

Vect a = v.itsVect;

for(int i = 0; i < a.size(); i++)

{

xout << " element [" << i << "] = " << a[i] << endl;

}

return xout;

}

Matr_op setMatr_op( Matr a)//a function for turning an ordinary matrix into a matrix filled with vectors Vect_op

{

Matr_op b;

for(int i = 0; i < a.size(); i++)

{

Vect_op v(a[i]);

b.push_back(v);

}

return b;

}

Matr getMatr( Matr_op a)

{

Matr b;

Vect v;

for(int i = 0; i < a.size(); i++)

{

v = a[i].getVect();

b.push_back(v);

}

return b;

}

private:

Vect itsVect;

};

int main()

{

Vect a1;

Vect b;

Vect c;

double x;

a1.push_back(3);

a1.push_back(7);

a1.push_back(-2);

b.push_back(5);

b.push_back(-4);

b.push_back(10);

Vect_op v(a1);

Vect_op w(b);

x = v * w;

cout << x << endl;

Vect_op y(a1);

y = v * x;

c = y.getVect();

for(int i = 0; i < c.size(); i++)

{

cout << " element " << i << " = " << c[i] << endl;

}

double z = w.sq_norm();

cout << " squared norm of w = " << z << endl;

Vect_op u(a1);

u = v + w;

c = u.getVect();

for(int i = 0; i < c.size(); i++)

{

cout << " element " << i << " = " << c[i] << endl;

}

cout << v[1] << endl;

cout << w << endl;

cout << u << endl;

int i;

int counter = 0;

int nb = 0;

cout << " number of rows or columns:\n";

cin >> nb;

cout << " Enter the square matrix:\n";

Vect row(nb);

Vect determinant;

Matr a;

int product = 1;

int determ = 0;

int ve1 = 0;

for(i = 0; i < nb; i++)

{

for(int j = 0; j < nb; j++)

{

if(j<(nb-1))

cout << " enter element\n";

else cout << " enter last element of this row\n";

cin >> row[j];

}

a.push_back(row);

}

for(i = 0; i < nb; i++)

{

cout << endl;

for(int j = 0; j < nb; j++)

{

cout << a[i][j] << " ";

}

}

cout << endl;

vector<int>vec;

for(i = 1; i <= nb; i++)

vec.push_back(i);

vector<int>::iterator iter1 = vec.begin();

vector<int>::iterator iter2 = vec.end();

for(i = 0; i < fact(nb); i++)

{

/* copy(vec.begin(), vec.end(), ostream_iterator<int>(cout, " "));;*/

vector<int>vec1(nb);

copy(vec.begin(), vec.end(), vec1.begin());

for(int j = 0; j < (nb); j++)

cout << vec1[j];

for(int i = 0; i < vec1.size(); i++)

{

for(int j = i; j < vec1.size(); j++)

{

if(vec1[j]<vec1[i])

counter++;

}

}

if(counter%2==0)

cout << " even ";

else

cout << " odd ";

for(int w = 0; w < nb; w++)

{

ve1 = (vec1[w])-1;

product *= (a[w][ve1]);

}

if(counter%2!=0)

product = (- product);

determinant.push_back(product);

counter = 0;

product = 1;

next_permutation(iter1, iter2);

cout << endl;

}

for(int x = 0; x < determinant.size(); x++)

{

determ += determinant[x];

}

for(i = 0; i < nb; i++)//print the matrix

{

cout << endl;

for(int j = 0; j < nb; j++)

{

cout << a[i][j] << " ";

}

cout << endl;

cout << " The determinant is:\n";

cout << determ << endl;

if(determ!=0) cout << " determinant is not zero \n";

else cout << " determinant equal zero: no QR decomposition possible\n";

Matr e;

Vect d;

for(int j = 0; j < (a[0].size()); j++)

{

for(int i = 0; i < a.size(); i++)

{

d.push_back(a[i][j]);

}

e.push_back(d);

d.clear();

}

Vect row1;

Matr q,f;

for(int i = 0; i < a.size(); i++)

{

for(int j = 0; j < a.size(); j++)

{

row1.push_back(0);

}

q.push_back(row1);

f.push_back(row1);

}

Matr_op qz, fz, ez;

qz = v.setMatr_op(q);

fz = v.setMatr_op(f);

ez = v.setMatr_op(e);

for(int i = 0; i < q.size(); i++)

{

for(int j =0; j < q.size(); j++)

{

if(i==0)

qz[i] = ez[i];

if(i!=0)

{

if(j==0)

fz[j].setVect(f[0]);

if(j!=0)

{

fz[j] = fz[j-1] + ( qz[j-1] * ( ez[i] * qz[j-1] ))* ( 1 / qz[j-1].sq_norm() );

if(j==i)

{

qz[i] = ez[i] - fz[j];

j = q.size();

}

fz = v.setMatr_op(f);

}

Matr qw;

qw = v.getMatr(qz);

//turn the vectors of qw into an orthonormal basis ( they are only orthogonals...)

for(int i = 0; i < qz.size(); i++)

{

double n = qz[i].sq_norm();

n = sqrt(n);

Vect ort; double w;

for(int j =0; j < (qz[0].getVect()).size(); j++)

{

w = (qz[i].getVect())[j] /= n;

ort.push_back(w);

}

qz[i].setVect(ort);

}

qw = v.getMatr(qz);

//finding r the upper triangular matrix:

Matr r = f;

Matr_op rz;

rz = v.setMatr_op(r);

for(int i = 0; i < ez.size(); i++)

{

Vect tria; double w;

for(int j =0; j < ez.size(); j++)

{

/*if(j<i)

{

w = ((rz[i]).getVect())[j] = 0;

tria.push_back(w);

}

else if(j>=i) */

{

w = ((rz[i]).getVect())[j] = ez[i] * qz[j];

tria.push_back(w);

}

rz[i].setVect(tria);

}

r = v.getMatr(rz);

// put the Q and R matrices upside down...( qw and r ):

Matr Q, R;

Q = InvertMatrix(qw);

R = InvertMatrix(r);

cout << " \n The Q matrix is:\n";

PrintMatrix(Q);

cout << " \n The R matrix is:\n";

PrintMatrix(R);

return 0;

}

Matr EnterMatrix()

{

int nb1 = 0;

int nb2 = 0;

Matr a;

cout << " Number of rows\n";

cin >> nb1;

cout << " Number of columns\n";

cin >> nb2;

cout << " Enter the matrix:\n";

Vect row(nb2);

for(int i = 0; i < nb1; i++)

{

for(int j = 0; j < nb2; j++)

{

if(j<(nb2-1))

cout << " enter element\n";

else cout << " enter last element of this row\n";

cin >> row[j];

}

a.push_back(row);

}

return a;

}

void PrintMatrix(Matr a)

{

int nb = a.size();

int nb1 = (a[0].size());

for(int i = 0; i < nb; i++)

{

cout << endl;

cout.width(1);

for(int j = 0; j < nb1; j++)

{

if((a[i][j])>0)

cout << setw(1) ;

//cout.width(5);<< " "<< " ";

cout << setw(14) << a[i][j] ;

//<< " "

}

cout << endl;

}

int fact(int a)

{

if(a<=1)

return a;

return (a * fact(a - 1));

}

Matr InvertMatrix(Matr a)

{

Matr e;

Vect d;

for(int j = 0; j < (a[0].size()); j++)

{

for(int i = 0; i < a.size(); i++)

{

d.push_back(a[i][j]);

}

e.push_back(d);

d.clear();

}

return e;

}

Размещено на Allbest.ru

...