9. Układy równań nieliniowych. Metoda Newtona
9.3. Metoda Newtona
W przypadku równania skalarnego, f : R → R, metoda stycznych (zwana też metodą New-tona) rozwiązywania równania f (x) = 0 jest zadana wzorem
xk+1= xk− f (xk) f0(xk).
Przez analogię, gdy F : RN ⊃ D → RN można więc byłoby zdefiniować wielowymiarową
metodę Newtona wzorem
xk+1 = xk− F0(xk)−1F (xk), (9.2) gdzie F0(xk) oznaczałoby macierz pochodnej F w punkcie xk. Jak za chwilę się przekonamy, jest to rzeczywiście bardzo dobry pomysł, a tak określona metoda zachowuje (z zachowaniem właściwych proporcji) wszystkie cechy metody Newtona znane nam z przypadku skalarnego!
9.3.1. Analiza zbieżności metody Newtona
Wielowymiarową metodę Newtona i różne jej warianty będziemy analizowali przy pewnych dość ogólnych założeniach dotyczących funkcji F : RN ⊃ D → RN. W skrócie, będziemy je nazywali za [9] założeniami standardowymi :
1. Zbiór D jest otwarty i niepusty, a F jest różniczkowalna w D. 2. Istnieje rozwiązanie x∗ ∈ D:
F (x∗) = 0. 3. Pochodna F0 : D → L(RN) jest lipschitzowska ze stałą L:
∃L > 0 kF0(x) − F0(y)k ¬ Lkx − yk ∀x, y ∈ D,
przy czym norma po lewej stronie nierówności jest normą operatorową indukowaną przez normę wektorową w RN obecną po prawej stronie, tzn. dla liniowego operatora A : RN →
RN,
kAk = sup
x6=0
kAxk kxk .
4. Macierz pochodnej w rozwiązaniu, F0(x∗), jest nieosobliwa.
Lematy techniczne
Przez K(x∗, δ) będziemy oznaczali kulę otwartą o środku w x∗ i promieniu δ,
K(x∗, δ) = {x ∈ RN : kx − x∗k < δ}.
Lemat 9.1 (użyteczna wersja twierdzenia o wartości średniej). Niech będą spełnione założenia standardowe i niech K będzie otwartą kulą w D. Wtedy dla każdego x, y ∈ K,
F (y) − F (x) =
Z 1 0
F0(x + t(y − x)) (y − x) dt.
Dowód. Ustalmy x, y ∈ K. Ponieważ K ⊂ D jest wypukły, to x + t(y − x) dla t ∈ [0, 1] i dla t ∈ [0, 1] jest dobrze określona funkcja g(t) = F (x + t(y − x)).
Na mocy założeń standardowych F0(·) jest ciągła na ¯K, zatem także g0(t) = F0(x + t(y −
x))(y − x) jest ciągła na [0, 1] (i całkowalna). Teza lematu wynika więc z podstawowego
twier-dzenia rachunku różniczkowego, że Z 1
0
108 9. Układy równań nieliniowych. Metoda Newtona Lemat 9.2 (o lokalnym oszacowaniu funkcji i pochodnej). Przy założeniach standardowych, istnieje dostatecznie małe δ > 0 takie, że dla dowolnego x ∈ K(x∗, δ) zachodzi:
1. kF0(x)k ¬ 2kF0(x∗)k,
2. kF0(x)−1k ¬ 2kF0(x∗)−1k,
3. 1
2kF0(x∗)−1kkx − x
∗k ¬ kF (x)k ¬ 2kF0(x∗)kkx − x∗k,
Ćwiczenie 9.1. Wykaż, że przy założeniach standardowych x∗ musi być izolowanym rozwią-zaniem równania F (x) = 0.
Rozwiązanie. Z lematu9.2o oszacowaniu funkcji, istnieje otoczenie U 3 x∗ takie, że 1
2kF0(x∗)−1kkx − x
∗k ¬ kF (x)k.
Jeśli więc dla ˜x leżącego w tym otoczeniu F (˜x) = 0, to znaczyłoby to, że k˜x − x∗k ¬ 0, a więc
˜
x = x∗. Czyli w U nie ma innych rozwiązań niż x∗.
Twierdzenie 9.3 (o zbieżności metody Newtona). Przy standardowych założeniach, istnieją C > 0 (dostatecznie duże) i δ > 0 (dostatecznie małe) takie, że jeśli x0 ∈ K(x∗, δ), to ciąg (xk) zadany metodą Newtona (9.2) jest dobrze określony,
xk∈ K(x∗, δ), oraz xk→ x∗. Co więcej, ciąg ten jest zbieżny kwadratowo:
kxk+1− x∗k ¬ Ckxk− x∗k2 ∀k ∈ N. (9.3)
Dowód. Dowód będzie opierał się na wykazaniu oszacowania (9.3), pozostałe elementy tezy będą jego konsekwencją.
Na wstępie wybierzmy δ takie, by zachodził lemat 9.2 o oszacowaniu funkcji i pochodnej. Oznaczając ek = xk− x∗ mamy ze wzoru Newtona
ek+1 = ek− F0(xk)−1F (xk). Na mocy założeń standardowych i lematu o wartości średniej,
F (xk) = F (xk) − F (x∗) = Z 1 0 F0(x∗+ t(xk− x∗))(xk− x∗) dt = Z 1 0 F0(x∗+ tek)ekdt, zatem ek+1= ek− F0(xk)−1 Z 1 0 F0(x∗+ tek)ekdt = F0(xk)−1 Z 1 0 F0(xk) − F0(x∗+ tek) ekdt.
Korzystając z lipschitzowskości pochodnej i raz jeszcze z lematu 9.2 o oszacowaniu pochodnej, dostajemy kek+1k ¬ kF0(xk)−1k Z 1 0 kF0(xk) − F0(x∗+ tek)k kekk dt ¬ 2kF0(x∗)−1k Lkekk2 Z 1 0 (1 − t) dt =: Ckekk2.
Stąd wynika, że jeśli Cδ < 1 oraz xk∈ K(x∗, δ), to także xk+1 ∈ K(x∗, δ), a więc ciąg (xk) jest dobrze określony. Co więcej, gdy Cδ < 1 to
kek+1k ¬ Ckekk2 = Ckekk kekk ¬ Cδkekk,
a więc ciąg ek jest zbieżny (co najmniej liniowo — a w rzeczywistości co najmniej kwadratowo) do zera.
Ćwiczenie 9.2. Wykaż, że jeśli w założeniach standardowych zastąpić warunek
lipschitzow-skości pochodnej warunkiem
∃α ∈ (0, 1] ∃L > 0 kF0(x) − F0(y)k ¬ Lkx − ykα ∀x, y ∈ D,
(czyli założyć, że F0 jest w swej dziedzinie h¨olderowska z wykładnikiem α), to metoda Newtona będzie lokalnie zbieżna z wykładnikiem zbieżności co najmniej 1 + α.
Wskazówka. Wystarczy powielić dowód twierdzenia przy standardowych założeniach, modyfiku-jąc oszacowanie kF0(xk) − F0(x∗+ tek)k.
Ćwiczenie 9.3. Jak, przy założeniach standardowych, oszacować normę błędu, kxk−x∗k, przez
normę residuum, kF (xk)k?
9.3.2. Implementacja metody Newtona
Implementując metodę Newtona, nie będziemy rzecz jasna nigdy explicite wyznaczali ma-cierzy odwrotnej F0(xk)−1, tylko rozwiązywali układ równań z macierzą F0(xk) — tę jednak będziemy już musieli wyznaczyć. Pamiętając także o tym, że metoda Newtona jest metodą ite-racyjną, stosując ją będziemy musieli zadbać o postawienie sensownego kryterium zatrzymania metody. Odkładając na bok pytanie o konkretny warunek stopu (por. rozdział 7.5), możemy schematycznie zapisać algorytm realizujący metodę Newtona w następujący sposób:
Schemat metody Newtona
function Newton(x, F, stop) while not stop do begin
oblicz macierz pochodnej F0(x); rozwiąż F0(x)s = F (x);
x = x − s;
end return(x);
Ćwiczenie 9.4. Macierzowe zadanie własne dla symetrycznej macierzy A ∈ RN ×N, znajdowa-nia pary (x, λ) ∈ RN× R spełniającej
Ax = λx, x 6= 0
można potraktować jako kwadratowe równanie nieliniowe dla F : RN +1 → RN +1 danej na przykład wzorem
F (x, λ) = Ax − λx
1 −12xTx
!
.
Zdefiniuj metodę Newtona dla F . Wykaż, że jeśli λ jest jednokrotną wartością własną, to metoda Newtona dla F będzie lokalnie zbieżna kwadratowo.
Rozwiązanie. Mamy
F0(x, λ) = A − λI −x −xT 0
!
.
Zatem metodę Newtona można — jeśli tylko F0(x, λ) nieosobliwa — zaimplementować, korzy-stając na przykład ze wzoru Shermana–Morrisona (por. ćwiczenie5.23).
110 9. Układy równań nieliniowych. Metoda Newtona
Zauważmy, że F jest funkcją gładką, jej pochodna jest lipschitzowska ze stałą Lipschitza równą... no właśnie, ile? Mamy
F0(x, λ) − F0(y, µ) = (µ − λ)I y − x
(y − x)T 0 !
,
zatem
kF0(x, λ) − F0(y, µ)k1 = max{|µ − λ| + |(y − x)1|, . . . , ky − xk1} ¬ ky − xk1+ |µ − λ| = k x λ ! − y µ ! k1.
A więc — w normie indukowanej normą k · k1 — F0 jest lipschitzowska ze stałą równą 1. Pozostaje jeszcze sprawdzić, czy w rozwiązaniu — parze własnej x∗, λ∗— macierz pochodnej jest nieosobliwa. Ponieważ dla macierzy symetrycznej A = QΛQT, gdzie Λ jest macierzą diagonalną z wartościami własnymi macierzy A, a kolumny macierzy ortogonalnej Q są wektorami własnymi odpowiadającymi tym wartościom własnym, to
QT 1 ! F0(x, λ) Q 1 ! = Λ − λI QTx xTQ 0 ! .
Niech λ∗ = λi. Wtedy x∗ = Cqi dla pewnej stałej C. Mamy więc, że Λ − λ∗I jest macierzą
diagonalną, która na diagonali ma zero jedynie na i-tej pozycji (na mocy założenia, że λ∗ jest pojedynczą wartością własną) oraz QTx = Cei, gdzie ei oznacza i-ty wektor jednostkowy. Nietrudno sprawdzić wprost, że taka macierz jest pełnego rzędu.
Ćwiczenie 9.5 (Metoda Schultza wyznaczania macierzy odwrotnej, [15]). Niech A będzie ma-cierzą nieosobliwą i rozważmy iterację
Xk+1= Xk+ Xk(I − AXk),
startującą z macierzy kwadratowej X0. Wykaż, że jeśli kI − AX0k < 1, to Xk → A−1 gdy
k → ∞.
Wskazówka. Wykaż, że kI − AXk+1k = kI − AXkk2, a metoda jest w rzeczywistości metodą Newtona zastosowaną do równania macierzowego F (X) = A − X−1 = 0.
Przykład 9.4 (Metoda Newtona dla równania Allena–Cahna). Przypomnijmy (zob.
przy-kład 9.1), że równanie nieliniowe, które nas interesuje w przypadku stacjonarnym ma postać
F (UN) = PNUN + f (UN), a w przypadku ewolucyjnym —
F (UN) = UN + τ (PNUN+ f (UN)). Oba przypadki możemy zatem ogarnąć, definiując funkcję
Fτ,α(UN) = αUN + τ (PNUN + f (UN)),
która dla α = 1 i τ 1 daje nam przypadek ewolucyjny, a dla α = 0 i τ = 1 redukuje się do przypadku stacjonarnego równania Allena–Cahna.
Wystarczy więc tylko wyznaczyć pochodną Fτ,α0 (UN). Ponieważ (traktując f jako funkcję na RN2)
to ∂fj ∂uk(UN) = ( 1 − 3u2j, k = j, 0, k 6= j, zatem ostatecznie Fτ,α0 (UN) V = (αI + τ (PN+ diag(1 − 3u2j))) V.
Przykład 9.5 (Numeryczne eksperymenty z metodą Newtona dla równania Allena–Cahna).
Dla uproszczenia, w eksperymentach numerycznych dotyczących działania różnych metod roz-wiązywania równania Allena–Cahna opisywanego w przykładzie 9.1, będziemy rozważać dys-kretyzacje jednowymiarowego stacjonarnego równania Allena–Cahna,
(TNu)j+ δuj(1 − u2j) = gj, j = 1, . . . , N,
gdzie δ = 10 oraz gj = j/(N + 1). %
% rozwiazania zadania z macierza rozrzedzona %
disp(’Matematyka obliczeniowa II’);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [x, resid, info, output] = newton(nazwa f, nazwa df, x0, atol, rtol, step, maxit)
%
% [x, resid, info, output] = newton(nazwa f, x0, atol, rtol, step, maxit)
To tylko fragment skryptu Octave. Możesz go uruchomić na http: //mst.mimuw.edu.pl/lecture.php?lecture=mo2&part=Ch9.
Ćwiczenie 9.6. Przeprowadź podobne eksperymenty numeryczne dla równania Allena–Cahna
10. Wariacje na temat metody Newtona
Podstawowy element metody Newtona — rozwiązywanie równania zlinearyzowanego — w przypadku, gdy N jest duże, może stanowić wąskie gardło całego procesu iteracyjnego. Dla-tego w tym i następnym rozdziale poszukamy skutecznych metod ominięcia Dla-tego ograniczenia; jednak na początek przytoczymy inną wersję twierdzenia o zbieżności, która (przy silniejszych założeniach) zagwarantuje nam także istnienie rozwiązań.