Metody numeryczne w fizyce

(1)

Metody numeryczne w fizyce

FZP002934wcL

rok akademicki 2016/17 semestr letni

Wykład 10

Karol Tarnowski

karol.tarnowski@pwr.edu.pl A-1 p. 411B

(2)

• Ortogonalizacja

• Rozkład QR

• Metoda QR wyznaczania wartości własnych

Plan wykładu

(3)

• D. Kincaid, W. Cheney, Analiza

numeryczna, Wydawnictwa Naukowo- Techniczne, 2006, rozdział 5

(4)

• iloczyn skalarny wektorów z

pozwala określić pojęcie ortogonalności

• układ wektorów v₁, v₂, …, v_k jest ortogonalny jeśli

• dodatkowo jeśli

układ jest ortonormalny

Ortogonalizacja

, Cⁿ

, 0 dla

i j

v v  i  j

, 1 dla 1

i i

v v   i k

(5)

• niech wektory v₁, v₂, …, v_n będą bazą ortonormalną przestrzeni

• wtedy każdy wektor wyraża się jednoznacznie jako kombinacja liniowa wektorów bazy

Ortogonalizacja

Cn

x C n

1 n

i i i

x c v



 

1 1 1

, , ,

n n n

j i i j i i j i ij j

i i i

x v c v v c v v c  c

  

      

1

,

n

i i i

x x v v



 

(6)

• dla danego ciągu wektorów liniowo niezależnych x₁, x₂, … z przestrzeni unitarnej procedura

Grama-Schmidta pozwala zbudować ciąg

ortonormalny u₁, u₂, … złożony z kombinacji liniowych tamtych wektorów

• ciąg wektorów

jest bazą ortonormalną podprzestrzeni span{x₁, x₂, … , x_k}

Ortogonalizacja Grama- Schmidta

1

2

, ,

k k k i i k k i i

i k i k

u x x u u x x u u



 

 

 



 





(7)

• dowód indukcyjny względem k

• dla k = 1

• załóżmy, że układ u₁, u₂, …, u_k-1 jest bazą ortonormalną podprzestrzeni

span{x₁, x₂, … , x_k-1}

• rozważmy

Ortogonalizacja Grama-Schmidta

Dowód

1

1 1 2 1

u  x ^ x

k k , i i

i k

v x x u u



 



, , , ,

, , , , 0

j k j k i i j

i k

k j k i ij k j k j

i k

v u x u x u u u

x u x u  x u x u



  

    



1 k 2

u  v ^ v

(8)

Ortogonalizacja Grama-Schmidta stosowana do kolumn A₁, A₂, …, A_n macierzy A rozmiaru m na n może być interpretowana jako jej rozkład QR na czynniki B i T. A = BT

Ortogonalizacja Grama-Schmidta

Algorytm

for j = 1 to n do

for i = 1 to j-1 do t_ij = <A_j,B_i>

enddo

C_j = A_j – S_i<jt_ijB_i t_jj = ||C_j||

B_j = t_jj^-1C_j enddo

(9)

B – kolumny stanowią układ ortogonalny T – na przekątnej są same jedynki

Ortogonalizacja Grama-Schmidta

Algorytm zmodyfikowany

for k = 1 to n do d_k = ||A_k||² t_kk = 1

for j = k+1 to n do t_kj = d_k^-1<A_j,B_k>

A_j = A_j - t_kjA_k enddo

enddo

(10)

Rozważmy układ m równań z n niewiadomymi Ax = b (m>n)

Szukamy takiego x, które spełnia ten układ

„najlepiej”, czyli takiego, że norma wektora residualnego b - Ax jest minimalna. Jeśli

rząd macierzy A jest n to ten wektor jest określony jednoznacznie.

Zadanie najmniejszych

kwadratów

(11)

Rozwiązaniem zadania najmniejszych kwadratów jest wektor x taki, że A^H(Ax-b) = 0.

Dowód. Rozważmy wektor y należący do Cⁿ.

Ponieważ A^H(Ax-b) = 0, więc wektor Ax – b jest ortogonalny względem przestrzeni rozpiętej na kolumnach macierzy A.

Wektor A(x - y) należy do tej przestrzeni więc

Zadanie najmniejszych kwadratów

Lemat

 

, 0

b Ax A x y  

   

2 2

2 2 2 2

2 2

b Ay b Ax Ax Ay

b Ax A x y b Ax A x y b Ax

     

         

(12)

Jeśli macierz A rozłożono na czynniki B i T, to rozwiązanie x układu Ax = b w sensie

najmniejszych kwadratów jest dokładnym rozwiązaniem układu

Zadanie najmniejszych kwadratów

 

^H ^H

Tx  B B ^¹ B b,

 

H H H

A Ax b  0 A Ax  A b

 

^H

 

H H H H H

A Ax  BT BTx T B BTx T  B B Tx

 

^H

  

H H H H H H H

A b  BT b T B b T  B B B B ^¹ B b

 

^{B B}^H ^¹ ^ ^diag



^d¹^¹^,^d²^¹^, ^,^dⁿ^¹



(13)

Zadanie najmniejszych kwadratów

A A AH

 

 

  

 

       

2

1 1 1

0 0 , 1 1 1 .

0 0

1 1 1

0 0

(14)

Celem metody jest rozkład A = QR macierzy A rozmiaru m na n, gdzie Q ma być macierzą

unitarną stopnia m, a R macierzą trójkątną górną rozmiaru m na n.

Algorytm daje Q^H i R takie, że Q^HA = R.

Q^H powstaje jako iloczyn macierzy unitarnych postaci

Metoda Householdera

0 .

0

k

H m k

I

I _ vv

 

  

 

(15)

Metoda zaczyna się od wyboru takiego v, żeby macierz I – vv^H była unitarna i żeby pierwsza

kolumna iloczynu (I – vv^H)A miała, jak w R, postać (b, 0, …, 0).

Niech A₁ oznacza pierwszą kolumnę A.

Chcemy więc, żeby było (I – vv^H)A₁ = b e⁽¹⁾.

Można to osiągnąć wybierając liczbę zespoloną

b taką, że oraz jest rzeczywiste.

Ponadto , gdzie

Metoda Householdera

A1 2

b  _A₁_,b_e^{ }¹



¹  ¹



v  A  be   ₁  b ^{ }¹ 2 A e 2.

(16)

Metoda Householdera

1 2

A ei^

b  a₁₁  a e₁₁ ⁱ^

 ¹  

1, 11 11 1 2 ⁱ

A be  a b  a A e ^{ }^      0   



¹  ¹



v   A  be

 



^ ^



   

1 11 11

11 11

i i

i i i

v a a e e

a e e a e

  

  

 b  b

     

   

(17)

Dlatego definiujemy b wzorem

Ostatecznie algorytm konstrukcji jest następujący

Metoda Householdera

1 2 1 2 11 11

A ei^ A a a b    

 

1 2 11 11 1

1

2 2 H

A a a

y A e

y

v y

U I vv b

b



 



 

(18)

Kolejne kroki rozkładu QR są podobne.

Po k krokach macierz A jest pomnożona z lewej przez k macierzy unitarnych, a wynik tych mnożeń jest macierzą blokową

gdzie J jest trójkątną górną.

Metoda Householdera

 

     

1 2 1 ,

0

k k k n k

k k

m k k m k n k

J H

U U U U A

W

  



    

 

  

 

 

(19)

Istnieje wektor taki, że I – vv^H jest macierzą unitarną stopnia m – k, a iloczyn

(I – vv^H)W ma wyzerowaną pierwszą kolumnę (bez pierwszego wyrazu).

Metoda Householdera

v C m k^

   

 

     

 



   



^ ^{ } ^

        

 

         

1 1 1 1

0

0 0

0

k k k n k

k k

H

m k m k m k k m k n k

k k k n k

H

m k k m k m k m k n k

k k k n k

m k k m k n k

J H

I

I vv W

J H

I vv W

J H

W

  



       

  

       

      

        

 

  

  

   

   

 

 

 

  

 

 

 

  

(20)

• Twierdzenie Schura gwarantuje, że dowolna macierz kwadratowa jest

unitarnie podobna do macierzy trójkątnej UAU^H = T.

• Metoda QR jest procedurą iteracyjną, prowadzącą do macierzy T.

Metoda QR obliczania wartości

własnych

(21)

Do rozkładu macierzy A na czynniki unitarny Q i trójkątny górny R możemy wykorzystać metodę Householdera.

A = QR

Teraz dodatkowo chcemy, aby elementy przekątniowe macierzy R były nieujemne.

Jeśli rozkład nie ma tej własności konstruujemy unitarną macierz przekątniową D tak, że

Metoda QR obliczania wartości własnych

dla 0 1 w p.p.

ii ii ii ii

ii

d r r r

d

 



(22)

Metoda QR obliczania wartości własnych

ˆ ˆ A QDD R QR H 

A₁ = A

for k = 1 to M do

rozkład A_k = Q_k R_k A_k+1 = R_k Q_k

enddo

1

H H

k k k k k k k k k k

A  Q R  Q R Q Q  Q A Q_