Choliy.NumericaMethods.2011Dec01

З формальної точки зору, метод Гауса полягає у зведеннi квадратної матрицi системи рiвнянь до верхньої трикутної. Отримана таким способом система еквiвалентна початковiй, однак її розв’язати значно простiше i проводиться це починаючи з останнього рiвняння:

xn = a(nnn−1) .

Пiдставивши знайдене значення у передостаннє рiвняння

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

LU-розбиття. Будь-яку квадратну невироджену матрицю A можна розкласти у добуток нижньої трикутної L та верхньої трикутної U:

A = L · U, lij = 0, i < j, uij = 0, i > j.

Такий розклад неоднозначний, тому зазвичай кладуть uii = 1. Для знаходження елементiв матриць L, U запишемо

n

X

likukj = aij,

k=1

але оскiльки матрицi трикутнi, то n = min(i, j) i тодi:

Тепер у кожнiй сумi вiдокремимо останнiй доданок:

31

i дещо переписавши, маємо формули для обчислень:

=1

якими користуємося так: першою при j = 1, i ≥ 1, другою при i = 1, j ≥ 1, потiм знову першою при j = 2, i ≥ 2, тодi другою при i = 2, j ≥ 2 i так до n.

Щоб розв’язати систему рiвнянь з використанням LU-розбиття введемо замiну U~x = ~y:

тобто, одна система рiвнянь з заповненою матрицею A, замiнена на двi системи рiвнянь з трикутними матрицями. Розписавши цi двi системи покомпонентно i прийнявши до уваги трикутну структуру матриць, маємо:

XX

Вiдокремивши останнiй доданок у першому рiвняннi i перший доданок у другому, отримуємо формули для обчислень:

У формулi (4.1) для uij в знаменнику стоїть lii. Чим бiльший модуль цього числа, тим точнiшими будуть результати обчислень. На цьому будується правило вибору головного елемента i робиться перестановка рядочкiв та стовбчикiв так, щоб на дiагоналi був найбiльший за модулем елемент. Це дає модифiкацiю методу: LU розбиття з вибором головного елемента.

Приклад:

32

LU-розбиття можна отримати i на основi iншого пiдходу. Нехай d(i) - одинична матриця, у якої елемент dii замiнений не рiвним одиницi числом, а l(i) - нижня трикутна матриця з одиничною дiагоналлю, у якої лише i-й стовбчик пiд дiагоналлю ненульовий, а всi iншi елементи - нулi. Множення довiльної матрицi A (для прикладу вiзьмемо матрицю 3 × 3) на d(1) а тодi на l(1) дає такий результат:

A, то результат такого множення такий же, як пiсля першого кроку методу Гауса. Продовжуючи множення, готуючи кожен раз матрицi l(i), d(i) на основi елементiв матрицi A, пiсля послiдовностi крокiв, отримуємо:

маємо

ˆ

LA = U,

ˆ

де U - верхня трикутна матриця з одиничною дiагоналлю (за побудовою), а L - нижня трикутна, (оскiльки добуток нижнiх трикутних матриць є нижньою

ˆ−1

трикутною матрицею). Якщо покласти L = L,

L · U = A.

Таким чином, LU-розбиття можна виконати множенням на послiдовнiсть вiдповiдно пiдiбраних елементарних матриць. Знаходження матрицi L не складає труднощiв, оскiльки матрицi, оберненi до l(i) та d(i), знаходяться дуже просто, наприклад:

33

Обернена матриця знаходиться таким способом. Очевидно, що

A · A−1 = E.

Поклавши A−1 = X маємо

A · X = E.

Оскiльки стовбчики матрицi з ненульовим визначником незалежнi, то останнє рiвняння можна розбити на n лiнiйних систем:

A · ~x(i) = ~e(i), i = 1, . . . , n,

де ~e(i) - орт, нульовий вектор з одиницею у i-й позицiї. Матрицю досить розбити один раз, а тодi розв’язати n лiнiйних систем. Знайденi вектори ~x(i) є стовбчиками оберненої матрицi.

Знаходження визначника грунтується на тому, що

n

Y det A = det L · det U = lii.

i=1

Метод Холецького. Якщо матриця A симетрична. LU-розбиття можна за-

мiнити на

L · LT = A.

Тепер матриця L, а значить i LT - трикутна з неодиничною дiагоналлю. Мiркуючи зовсiм аналогiчно, запишемо

n

X

likukj = aij,

k=1

але U = LT , тому

n

X

likljk = aij.

k=1

Тут обов’язково i ≥ j, що дає

j

X

likljk = aij,

k=1

а пiсля вiдокремлення останнього доданка в сумi:

j−1

X

lik · ljk + lijljj = aij.

k=1

34

Це дає формули для обчислень:

Формули повиннi використовуватися у циклi по j = 1, . . . , n всерединi якого цикл по i = j +1, . . . , n. Iнша назва цього методу - метод квадратного кореня. QR-розбиття. Це ще одне розбиття матрицi на двi

A = QR,

але тепер матриця R - верхня трикутна, а матриця Q - ортогональна, для якої QQT = E, Q−1 = QT , тобто її транспонована спiвпадає з оберненою. Саме ця остання властивiсть вiдiграє головну роль у методi, що розглядається далi, а сам пiдхiд має дуже важливе значення також i при вiдшуканнi власних чисел та векторiв.

Перетворення Хаусхолдера. Знайдемо послiдовнiсть ортогональних матриць, якi утворять Q. Нехай w~ - довiльний вектор, для якого |w~| = 1. Розглянемо

матрицю

P = E − 2w~w~T .

Перевiримо її ортогональнiсть:

P 2 = (E − 2w~w~T )2 = E2 − 2 · E · 2w~w~T + 2w~w~T · 2w~w~T = = E2 − 4w~w~T + 4w~(w~T w~)w~T = E − 4w~w~T + 4w~w~T = P.

Отже P · P = E, або P = P −1. Далi розглянемо

P T = ET − (2w~w~T )T = E − 2w~w~T = E.

Отже P = P T , а разом з P = P −1 це означає, що матриця P ортогональна. З довiльного вектора ~u матрицю P утворюємо так:

Тепер вiзьмемо ~x - перший стовбчик матрицi A i побудуємо вектор:

~u = ~x |x| · ~e1,

де ~e1 - орт з одиницею в першому елементi i нулями у всiх iнших.

35

Запишемо чому рiвне H:

H = 12 (~x |x|~e1) (~x |x|~e1) = 12 |x|2 2|x|x1 + |x|2 = |x|2 |x|x1.

Знайдемо добуток:

Тобто, множення на матрицю P прибирає у векторi ~x всi елементи, окрiм першого. Максимiзацiя модуля вектора H покращує точнiсть результату (H в знаменнику), тому при побудовi вектора ~u знак перед |x| слiд вибирати так, щоб вiн спiвпадав зi знаком x1.

Приклад. Знайти QR-розбиття для матрицi

36

Другий крок. Працюємо з субматрицею 3 ×3, що утвориться, якщо з резуль-

Допишемо до отриманої матрицi зверху та злiва по додатковому стовбчику:

Третiй крок. Працюємо з субматрицею 2 × 2, вiдкинувши два першi рядочки

37

Щоб обчислити P3 до знайденої субматрицi 2 × 2 дописано два рядки зверху та два стовбчики злiва, як i в попередньому кроцi:

матриця R - трикутна. Визначник знаходиться аналогiчно:

n

Y det A = det Q · det R = (−1)n−1 rii.

i=1

Тут використано той факт, що i det Pi = −1.

Уточнення розв’язку. Розв’язок системи лiнiйних рiвнянь можна уточнити. Позначимо через ~x(0) отриманий розв’язок. Пiдставимо його в систему, маємо:

розв’язком. Процес можна продовжити до збiжностi, доки | ~x(i)| не зменшиться до бажаної точностi.

38

Дiючи згiдно з означенням норми, i вiднявши вiд другої системи першу,

Величина kA−1k·kAk називається числом зумовленостi, цей множник показує як вiдносна похибка правих частин переходить у вiдносну похибку результату. Теоретичний мiнiмум для нього - одиниця. Чим ближче значення числа зумовленостi до одиницi, тим ”кращий“ розв’язок.

Тридiагональнi матрицi. Загальний вигляд системи з такою матрицею:

39

Для розв’язання такої системи, комбiнуючи рядочки, прибираємо елементи ai, при цьому будуть змiнюватися також елементи bi та елементи вектора yi:

Новi значення записано злiва вiд знака рiвностi, попереднi - справа. Якщо зупинитися за два рядочки до кiнця, то залишкову систему, яка матиме вигляд

b0 c0 x y0

n−1 n−1 · n−1 = n−1 ,

можна розв’язати навiть методом Крамера. Маючи два значення xi−1 та xi можна знайти xi−2:

xi−2 = (yi0−1 − b0i−1xi−1 − c0i−1xi) .

Скориставшись цим рiвнянням для i −2, i −3, . . . , 1 знайдемо усi невiдомi xi. Цей спосiб називається методом прогонки.

Звiдси формула для знаходження X1:

X1 = (B1 − A12A−221B2) · (A11 − A12A−221A21)−1.

Аналогiчно поступивши з першим рiвнянням: (використавши A−111 та A21)

отримаємо:

X2 = (B2 − A21A−111B1) · (A22 − A21A−111A12)−1.

Єдиною складною задачею є знаходження A12A−221, але, оскiльки розмiр цiєї матрицi менший, нiж матрицi A, це може зменшити загальну кiлькiсть обчислювальної роботи.

40