Rozwiązywanie układów równań nieliniowych

Draft wersja 0.12

Artur Wilkowski

20 listopada 2013

Streszczenie

Dokument obejmuje zagadnienia związane z metodami rozwiązy-

wania układów równań nieliniowych (poszukiwania miejsc zerowych

funkcji wielu zmiennych). W dokumencie przedstawione są następują-

ce metody: metoda iteracji prostej oraz metoda Newtona

1

Opis problemu

Dana jest funkcja wektorowa y = f (x) i f : D →

n

n

R , gdzie D ⊂ R .

Należy wyznaczyć jej miejsca zerowe, czyli rozwiązać równanie f (x) = 0.

Odpowiednikiem problemu jest zagadnienie rozwiązania układu równań w postaci

 f



1( x 1 , . . . , xn) = 0





. . .





 fn( x 1 , . . . , xn) = 0

jeśli rozwiniemy oznaczenia x = [ x 1 , . . . , xn] T , f = [ f 1 , . . . , fn] T i 0 = [0 , . . . , 0] T .

2

Sposób rozwiązania

Poszukiwanie zerowych funkcji możemy podzielić na dwa etapy:

1. Lokalizacja miejsc zerowych - określenie obszarów w których znajduje się jedno miejsce zerowe

2. Określenie dokładnej wartości miejsc zerowych

1

W kolejnych sekcjach skupimy się na iteracyjnych metodach określenia do-kładnych wartości miejsc zerowych. W dalszych częściach opisane zostaną metody lokalizacji miejsc zerowych dla funkcji wielu zmiennych.

3

Metody iteracyjne

Do rozwiązania problemu stosujemy metody iteracyjne - kolejne przybliżenia miejsca zerowego obliczane są na podstawie przybliżeń poprzednich. Metody, którymi będziemy się posługiwać to metoda iteracji prostej oraz metoda Newtona dla funkcji wielowymiarowych.

3.1

Metoda iteracji prostej

Równanie w postaci f (x) = 0 przekształcamy do postaci x = g(x) i aż do osiągnięcia zbieżności prowadzimy proces iteracyjny

x( k+1) = g(x( k))

x(0) - przybliżenie początkowe

Jeżeli spełnione są założenia twierdzenia Banacha o kontrakcji, czyli D jest zbiorem domkniętym w

n

R i funkcja g : D → D spełnia warunek Lipschitza ze stałą 0 ¬ L < 1 taką że

kg(x 0) − g(x 00) k ¬ Lkx 0 − x 00k, dla każdego x 0, x 00 ∈ D

to istnieje jedno rozwiązanie równania f (x) = 0 w D i metoda iteracji prostej jest zbieżna.

Uwaga. Zbiór jest domknięty w przestrzeni metrycznej, jeżeli dla dowolnego ciągu zbieżnego o wyrazach należących do zbioru, jego granica również należy do zbioru.

Jeżeli D jest wypukły i g ∈ C 1( D), do zbieżności wystarczy aby max kg0(x) k < 1

x ∈D

Uwaga. Zbiór wypukły to podzbiór przestrzeni Euklidesowej, o tej własności że każdy odcinek którego końce znajdują sie w zbiorze, w całości zawiera się w zbiorze.

2

Należy zauważyć, że w przypadku gdy spełnione są założenia twierdzenia Banacha w zbiorze D możliwe jest wstępne określenie wymaganej liczby iteracji niezbędnych do osiągnięcia odpowiedniej dokładności. Korzystając wielokrotnie z warunku Lipschitza i wiedząc, że miejsce zerowe x ∗ jest punk-tem stałym przekształcenia g (zatem g(x ∗) = x ∗) otrzymujemy oszacowanie kx( k) − x ∗k ¬ Lkkx(0) − x ∗k Teraz, ażeby uzyskać wymaganą dokładność kx( k) − x ∗k ¬ ε, wystarczy aby Lkkx(0) − x ∗k ¬ ε

Rozwiązując tą nierówność ze względu na k otrzymujemy oszacowanie liczby iteracji metody, które należy wykonać w celu uzyskania żądanej dokładności ε

k ­ log L kx(0) − x ∗k

Odległość między przybliżeniem początkowym i pierwiastkiem kx(0) − x ∗k można ograniczyć z góry na przykład przez średnicę d obszaru (zbioru), w którym zlokalizowaliśmy miejsce zerowe, wtedy liczba iteracji może być okre-

ślona jako

ε

k ­ log L d

Uwaga. Średnicą zbioru nazywamy supremum odległości wszystkich par punktów tego zbioru

3.2

Metoda Newtona

Metoda Newtona dla przypadku funkcji wielowymiarowych określana jest analogicznie do przypadku jednowymiarowego. Przypomnijmy, dla liczby wymiarów n = 1 mieliśmy wzór iteracyjny

f ( x( k))

x( k+1) = x( k) − f0( x( k)) Dla n ­ 1 mamy natomiast

x( k+1) = x( k) − [ f 0(x( k))] − 1 f (x( k)), i x(0) - przybliżenie początkowe gdzie [ f 0(x k)] − 1 jest odwrotnością macierzy Jacobiego.

3

W praktyce nie odwracamy macierzy Jacobiego, zamiast tego w każdej iteracji metody rozwiązujemy następujący układ równań

[ f 0(x k)](x k+1 − x k) = f (x k) Pochodne wyznaczamy w sposób analityczny, jeżeli jest taka możliwość, w przeciwnym przypadku stosujemy różniczkowanie numeryczne.

4

Lokalizacja miejsc zerowych funkcji wielu

zmiennych

Celem zadania lokalizacji miejsc zerowych, jest określenie liczby miejsc zerowych oraz określenie obszarów w których te miejsca zerowe się znajdują, w ten sposób żeby w każdym obszarze było dokładnie jedno miejsce zerowe.

Podczas lokalizacji miejsc zerowych funkcji można się wspomagać następu-jącymi sposobami:

1. Tablicowanie funkcji. Tworzymy tablicę wartości funkcji, ewentual-nie rysujemy wykres funkcji i dokonujemy wizualnej analizy przebiegu funkcji. Metoda ta jest znacznie trudniejsza do realizacji niż w przypadku jednowymiarowym, gdyż teraz poszukujemy miejsc, w których

zmienia się jednocześnie znak kilku funkcji skalarnych (zakładamy funkcje ciągłe). Dobrym pomysłem na wstępne ograniczenie zakresu poszukiwań jest obliczenie normy wektora wartości funkcji kf (x) k dla siatki punktów i skoncentrowanie się na obszarach o bliskich zeru wartościach kf (x) k.

2. Metoda graficzna. Jest szczególnie użyteczna dla układów dwóch równań (o dwóch zmiennych). Każde z równań określa pewną krzywą

na płaszczyźnie. Po narysowaniu obu krzywych określamy obszary, w których krzywe mają punkty przecięcia.

5

Praktyczne kryteria stopu metod

Wyznaczając kolejne przybliżenia x( k) musimy określić moment, w którym przerwiemy obliczenia przyjmując, że osiągnęliśmy już wymaganą dokład-ność. W przypadku metody iteracji prostej, jeżeli jesteśmy w stanie określić 4

parametr L warunku Lipschitza, można na wstępie oszacować wymaganą liczbę iteracji prowadzących do uzyskania odpowiedniej dokładności. W ogólnym przypadku jednak odwołujemy się najczęściej do praktycznych kryteriów stopu, które nie nakładają silnych wymagań na postać analizowanej funkcji, czy metodę rozwiązania.

Praktyczne kryteria stopu w stosowane w rozwiązywaniu układów rów-nań liniowych są zbliżone do kryteriów stopu dla równań jednej zmiennej, jednak tutaj odległość między punktami jest charakteryzowana przez normy wektorowe, a nie wartości bezwzględne.

Dla wybranego ε praktyczne kryteria stopu mogą zostać określone jako 1. kx( k+1) − x( k) k ¬ ε

2. kx( k+1) − x( k) k ¬ εkx( k+1) k 3. kf (x( k)) k ¬ ε

6

Podstawy matematyczne

6.1

Macierz Jacobiego

Macierz Jacobiego jest macierzą pochodnych cząstkowych pierwszego rzędu funkcji. Dla funkcji f : D →

n

n

R , gdzie D ⊂ R

i x = [ x 1 , . . . , xn] T , f =

[ f 1 , . . . , fn] T macierz Jacobiego możemy określić jako

 ∂f 1

. . .

∂f 1 

( ∂f ) n

∂x 1

∂xn

J

i





f (x) =

=

. . .

. . .

. . .

∂x





j

i,j=1

∂fn

. . .

∂fn

∂x 1

∂xn

Macierz Jacobiego określa pochodną funkcji różniczkowalnej.

Przykład. Mamy daną następującą funkcję dwóch zmiennych f ( x 1 , x 2) =

[ f 1( x 1 , x 2) , f 2( x 1 , x 2)] T



 f 1( x 1 , x 2) = x 2 + x 3 + 4

1

2

 f 2( x 1 , x 2) = sin x 1 + cos x 2

Pochodna funkcji f jest określana przez macierz Jacobiego i ma postać

!

2 x

f 0( x

1

3 x 22

1 , x 2) = Jf ( x 1 , x 2) =

cos x 1 − sin x 2

5

6.2

Normy wektorów i macierzy

6.2.1

Normy wektorów

Normy wektorów umożliwiają określenie długości wektorów. Norma zastoso-wana do różnicy pomiędzy wektorami kx − y k umożliwia określenie stopnia podobieństwa pomiędzy nimi (w przypadku dwóch wektorów identycznych uzyskujemy najmniejszą możliwą wartość normy, czyli 0). Ważnymi normami są normy Schurra określone równaniem

n

!1 /p

kx k

X

p =

|xi|p

dla 1 ¬ p ¬ ∞

i=1

Wśród których w powszechnym użyciu jest norma Euklidesowa

v

u

n

kx k

uX

2 = t

x 2 i

i=1

jak również norma Manhattańska

n

kx k

X

1 =

|xi|

i=1

czy norma maksimum (nieskończoność)

kx k∞ = max |xi|

1 ¬i¬n

6.2.2

Normy macierzy kwadratowych

Wśród norm macierzy w użyciu są tzw. normy indukowane przez odpowiednie normy wektorów. Normy indukowaną macierzy A określa się jako kAk = max kAx k

kx k=1

gdzie kx k oraz kAx k jest wybraną normą wektorową. Wśród przykładów norm indukowanych macierzy są normy indukowane przez wektorowe normy Schurra ( p-normy):

n

kAk

X

1 = max

|aij|

1 ¬j¬n i=1

6

lub

n

kAk

X

∞ = max

|aij|

1 ¬i¬n j=1

oraz

q

kAk 2 =

λ max( AT A)

gdzie funkcja λ max oznacza maksymalną wartość własną macierzy.

Wśród norm, które nie są normami indukowanymi używana jest np. norma Frobeniusa, obliczana w sposób zbliżony do normy Euklidesowej dla wektorów

v

n

n

u

kAk

uX X

F =

a 2

t

ij

i=1 j=1

7

Zadania

Zadanie 1. Zlokalizować rozwiązania układu równań



 x 2 + x 2 = 1

1

2

 x 2 = x 3

1

następnie użyć metody Newtona oraz metody iteracji prostej do odnalezienia rozwiązań, tak by spełnione było dokładności w postaci kf (x) k 2 ¬ 10 − 9 , gdzie x = ( x 1 , x 2) jest odnalezionym rozwiązaniem Zadanie 2. Zlokalizować rozwiązania układu równań



 x 2 + x 2 = 1

1

2

 x 3 + x 3 = x

1

2

1 x 2

następnie użyć metody Newtona oraz metody iteracji prostej do odnalezienia rozwiązań, tak by spełnione było kryterium dokładności w postaci kf (x) k 2 ¬

ε, gdzie ε = 10 − 9 i x = [ x 1 , x 2] T jest odnalezionym rozwiązaniem Zadanie 3. Rozwiązać za pomocą metody iteracji prostej układ równań



 x 1 = 1 e sin( x 1+ x 2)

7

 x 2 = 3 cos x

5

1 sin x 2

rozpoczynając obliczenia od punku x = [0 , 0] T . Przyjąć praktyczne kryterium stopu kx( k+1) − x( k) k ¬ εkx( k+1) k dla ε = 10 − 6 . Zbadać teoretycznie zbieżność metody iteracji prostej.

7