Zagadnienia

10.1 Twierdzenie o istnieniu rozwiązań
10.2 Obniżanie kosztu iteracji

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. Dlatego w tym i następnym rozdziale poszukamy skutecznych metod ominięcia 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ń.

10.1. Twierdzenie o istnieniu rozwiązań

W dotychczasowej analizie zakładaliśmy, że rozwiązanie x* równania F⁢x=0 istnieje. Okazuje się, że nieco modyfikując założenia standardowe można udowodnić i istnienie rozwiązania, i lokalną zbieżność metody Newtona do tego rozwiązania.

Twierdzenie 10.1 (Kantorowicza, o istnieniu i zbieżności)

Niech F:RN⊃D→RN, przy czym D jest niepusty, otwarty i wypukły, i niech F będzie różniczkowalna w D oraz istnieje L>0 takie, że

F′⁢x-F′⁢y≤L⁢x-y⁢⁢∀x,y∈D.

Niech istnieje x0∈D takie, że F′⁢x0 jest nieosobliwa oraz dla pewnych dodatnich β,η zachodzi

F′⁢x0-1⁢F⁢x0≤β,⁢F′⁢x0-1≤η,

przy czym α=L⁢β⁢η≤1/2. Określmy δ=1-1-2⁢α/L⁢β i załóżmy, że K¯⁢x0,δ⊂D.

Wówczas ciąg xk zdefiniowany metodą Newtona jest dobrze określony, co więcej xk∈K¯⁢x0,δ oraz ma granicę x*, która jest jedynym miejscem zerowym F w K¯⁢x0,δ. Ponadto

xk-x*≤L⁢β⁢2k⁢2⁢α2k

dla k≥0.

Dowód

Dowód tego twierdzenia można znaleźć np. w [12, rozdział 12.6.1].

∎

10.2. Obniżanie kosztu iteracji

Z punktu widzenia praktycznej realizacji, każdy krok metody Newtona ma parę potencjalnie kosztownych momentów, które teraz przedyskutujemy.

Wyznaczenie wzoru na F′⁢x może być trudne: na przykład, gdy F jest zadana bardzo skomplikowaną procedurą (a nie jednolinijkowym wzorem), wtedy ani ręczne, ani automatyczne różniczkowanie F może nie prowadzić do poprawnych czy zadowalających rezultatów. Może też być bardzo żmudne.
Wyznaczenie wartości F⁢x, a zwłaszcza: F′⁢x, może być kosztowne. Rzeczywiście, F′⁢x to N2 elementów macierzy pochodnej, a dodatkowo — poza najprostszymi przypadkami — pochodna jest zwykle dana bardziej skomplikowanym wyrażeniem niż sama funkcja. Rozsądne jest więc przyjęcie¹³O ile F′⁢x nie jest macierzą rozrzedzoną., że w praktyce koszt wyznaczenia F′⁢x będzie rzędu c⋅N⋅(koszt obliczenia F(x)), z raczej dużą stałą c. Nawet wyznaczenie pojedynczej wartości F⁢x może być w niektórych zadaniach bardzo kosztowne. W [3] wspomina się o F:R50→R50, której jedna wartość kosztowała 100 godzin pracy ówczesnego szybkiego komputera.
Rozwiązanie układu z F′⁢x jest potencjalnie najbardziej kosztownym fragmentem iteracji. Rzeczywiście, gdy F′⁢x jest macierzą gęstą bez żadnych szczególnych własności, to koszt rozwiązania układu

F′⁢x⁢s=F⁢x

jest rzędu O⁢N3.

Te obserwacje prowadzą do kilku heurystycznych sposobów modyfikacji metody Newtona tak, by obniżyć koszt jednej jej iteracji. Czasem — choć nie zawsze, i dlatego te sposoby są ważne! — obniżenie kosztu iteracji będzie skutkowało pogorszeniem szybkości zbieżności metody, ale nawet i wówczas może okazać się, że ostatecznie metoda wolniej zbieżna, ale tańsza, będzie bardziej efektywna od metody Newtona.

10.2.1. Uproszczona metoda Newtona

W uproszczonej metodzie Newtona, najdroższy fragment iteracji Newtona — wyznaczenie rozkładu trójkątnego macierzy pochodnej — usuwamy poza pętlę:

xk+1=xk-F′⁢x0-1⁢F⁢xk.

Rzeczywiście, jeśli bowiem na samym początku iteracji dokonamy kosztem O⁢N3 flopów rozkładu LU (lub QR) macierzy F′⁢x0 i zapamiętamy czynniki rozkładu, to potem każdy z układów równań:

F′⁢x0⁢s=F⁢xk

będziemy mogli rozwiązać już tylko kosztem nie wyższym niż O⁢N2! Oznacza to, że koszt jednej iteracji spada nam do O⁢N2, jest więc N-krotnie mniejszy niż w przypadku oryginalnej metody Newtona. Ceną za przyspieszenie jest wyraźne zmniejszenie szybkości zbieżności metody:

Twierdzenie 10.2 (o zbieżności uproszczonej metody Newtona)

Przy standardowych założeniach, istnieją stałe dodatnie C i δ takie, że jeśli x0∈K⁢x*,δ, to uproszczona metoda Newtona jest zbieżna przynajmniej liniowo do x* oraz

xk+1-x*≤C⁢x0-x*⁢xk-x*⁢⁢∀k∈N.

(10.1)

Dowód

Zostanie podany później, gdy będziemy dysponowali wygodnym narzędziem do badania tego typu metod.

∎

Ćwiczenie 10.1

Zaimplementuj uproszczoną metodę Newtona.

Rozwiązanie:

Wystarczy zmodyfikować algorytm realizujący metodę Newtona:

Uproszczona metoda Newtona

function simplerNewton(x, F, stop)

oblicz macierz pochodnej F′⁢x;

wyznacz rozkład F′⁢x, np. Q⁢R=F′⁢x;

while not stop do begin

rozwiąż układ z macierzą trójkątną, R⁢s=QT⁢F⁢x;

x = x - s;

end

return(x);

Ćwiczenie 10.2

Metoda Szamańskiego to metoda, w której po m krokach uproszczonej metody Newtona dokonujemy wyznaczenia i rozkładu macierzy pochodnej (czyli restartujemy metodę uproszczoną, biorąc za przybliżenie początkowe ostatnio wyznaczone przybliżenie rozwiązania). Jest to więc coś pośredniego pomiędzy prawdziwą metodą Newtona (w której w każdej iteracji wyznacza się i rozkłada macierz pochodnej) a uproszczoną metodą Newtona (gdzie macierz pochodnej wyznacza się tylko jeden raz).

Formalnie możemy więc przyjąć, że m kroków uproszczonej metody Newtona to jest ,,jeden duży krok” i określić iterację Szamańskiego następująco:

	xk+1m	=xk-F′⁢xk-1⁢F⁢xk
	xk+j+1m	=xk+jm-F′(xk)-1F(xk+jm),j=1,…,m-1.

Znajdź wykładnik p taki, że

xk+1-x*≤C⁢xk-x*p.

Rozwiązanie:

Z definicji metody, xk+1 wyznaczamy jako m-tą iterację uproszczonej metody Newtona startującej z xk. Na mocy oszacowania błędu uproszczonej metody Newtona mamy więc

xk+j+1m-x*≤C⁢xk-x*⁢xk+jm-x*,

skąd

xk+1-x*≤C⁢xk-x*⁢xk+m-1m-x*≤C⁢xk-x*2⁢xk+m-2m-x*≤⋯≤Cm⁢xk-x*m+1.

A więc metoda ma wykładnik m+1.

Ćwiczenie 10.3

Porównaj efektywność metody Szamańskiego i Newtona w zależności od rozmiaru zadania N, przy następujących modelowych założeniach:

koszt wyznaczenia F⁢x jest równy c0⁢N
koszt wyznaczenia F′⁢x jest równy c1⁢N2
koszt wyznaczenia rozkładu F′⁢x jest równy cr⁢N3
koszt wyznaczenia rozwiązania układu z macierzą F′⁢x o danym rozkładzie jest równy cs⁢N2

Rozwiązanie:

Oczywiście musimy założyć m>1, w przeciwnym razie obie metody są tożsame. Dla metody Newtona mamy

ENewton=21/cr⁢N3+c1+cs⁢N2+c0⁢N,

a dla metody Szamańskiego

ESzamański=m+11/cr⁢N3+c1+m⁢cs⁢N2+m⁢c0⁢N.

Stąd ESzamański≥ENewton wtedy i tylko wtedy, gdy log2⁡ESzamański≥log2⁡ENewton. Dalej przekształcając, dostajemy warunek wiążący ze sobą m i N:

log2⁡m+1≥1+m⁢⁢O⁢1N1+O⁢1N,

co zachodzi dla N dostatecznie dużych względem m.

Ćwiczenie 10.4 (liniowa niezmienniczość iteracji Newtona)

Rozważmy równanie F⁢x=0 i jego liniową transformację

G⁢y=A⁢F⁢B-1⁢y=0,

gdzie A,B są nieosobliwymi macierzami rozmiaru N×N. Wykaż, że jeśli tylko y0=B⁢x0, to iteracja Newtona dla G,

yk+1=yk-G′⁢yk-1⁢G⁢yk

jest równoważna metodzie Newtona dla F w tym sensie, że zawsze zachodzi yk=B⁢xk.

Rozwiązanie:

Iteracja Newtona dla F to

xk+1=xk-F′⁢xk-1⁢F⁢xk.

Mnożąc lewostronnie przez B i podstawiając gdzie trzeba xk=B-1⁢yk mamy

yk+1=yk-B⁢F′⁢B-1⁢yk-1⁢F⁢B-1⁢yk=yk-B⁢F′⁢B-1⁢yk-1⁢A-1⁢⁢A⁢F⁢B-1⁢yk.

Na zakończenie wystarczy zauważyć, że B⁢F′⁢B-1⁢yk-1⁢A-1=A⁢F′⁢B-1⁢yk⁢B-1-1=G′⁢yk-1.

Ćwiczenie 10.5

Wykaż, że jeśli spełnione są założenia standardowe oraz metoda Newtona jest zbieżna, to dla dostatecznie dużych k zachodzi

F′⁢xk-1⁢F⁢xk+1≤14⁢F′⁢xk-1⁢F⁢xk.

Może to więc być dość prosty test (niezmienniczy ze względu na liniowe transformacje zmiennej niezależnej i zależnej!), pozwalający przypuszczać, że nie zachodzi zbieżność: na przykład, gdy na pewnym kroku metody stwierdzimy F′⁢xk-1⁢F⁢xk+1≥34⁢F′⁢xk-1⁢F⁢xk, możemy wtedy uznać prawdopodobny brak zbieżności i zakończyć iterację.

10.2.2. Wpływ niedokładności wyznaczenia pochodnej lub funkcji na zbieżność metody Newtona

Twierdzenie 10.3 (o lokalnej stabilności metody Newtona)

Przy standardowych założeniach, istnieją stałe dodatnie C, δ oraz δ~ takie, że jeśli xk∈K⁢x*,δ oraz Δk<δ~, to

xk+1=xk-F′⁢xk+Δk-1⁢F⁢xk+ϵk

(10.2)

jest dobrze określony oraz

ek+1≤C⁢ek2+Δk⁢ek+ϵk,

(10.3)

gdzie ek=xk-x*.

Dowód

Dowód znajdziemy np. w [9].

∎

Z twierdzenia o lokalnej stabilności metody Newtona można wywnioskować wiele twierdzeń o zbieżności tych modyfikacji metody Newtona, które sprowadzają się do zaburzenia oryginalnej metody Newtona, na przykład twierdzenia o zbieżności uproszczonej metody Newtona.

Dowód

(Twierdzenia 10.2, o zbieżności uproszczonej metody Newtona.) Jeden krok uproszczonej metody Newtona

xk+1=xk-F′⁢x0-1⁢F⁢xk

możemy zapisać w języku twierdzenia o lokalnej stabilności metody Newtona, gdzie

Δk=F′⁢x0-F′⁢xk,⁢ϵk=0.

Niech więc δ będzie dostatecznie małe tak, by dla xk zachodziło twierdzenie o lokalnej stabilności metody Newtona w kuli K⁢x*,δ. Dodatkowo załóżmy, że x0∈K⁢x*,δ. Wtedy

Δk=F′⁢x0-F′⁢xk≤L⁢x0-xk≤L⁢x0-x*+xk-x*

i na mocy właśnie tw. o lokalnej stabilności,

ek+1≤C~⁢ek2+L⁢e0+ek⁢ek≤C⁢e0⁢ek.

Ponieważ z założenia e0<δ, to z powyższego wynika, że

ek+1≤C⁢δ⁢ek

jeśli więc δ jest na tyle małe, by dodatkowo C⁢δ<1, to ciąg xk zawiera się w K⁢x*,δ oraz ek musi być zbieżny do zera przynajmniej liniowo.

∎

10.2.3. Metoda z przybliżoną pochodną

Jak już wspominaliśmy, innym kłopotliwym momentem w realizacji metody Newtona jest wyznaczanie macierzy pochodnej. Opracowanie dokładnej formuły obliczania F′⁢x może być żmudne, podatne na ludzkie pomyłki. Uzyskany wzór może być ciężki w implementacji i kosztowny w użyciu.

Ale, ponieważ metoda Newtona jest jedynie metodą przybliżoną, można zaryzykować użycie przybliżonej pochodnej — najlepiej: przybliżonej niskim kosztem. Jednym z pomysłów mogłoby być wyznaczenie przybliżenia pochodnej na podstawie ilorazów różnicowych, zgodnie z poniższym twierdzeniem.

Twierdzenie 10.4 (o różnicowej aproksymacji pochodnej kierunkowej)

Przyjmijmy, że są spełnione założenia standardowe. Niech x będzie ustalonym punktem w D i niech w∈RN. Niech h∈R będzie na tyle małe, by x+h⁢w∈D. Wtedy

F⁢x+h⁢w-F⁢xh-F′⁢x⁢w≤L2⁢w2⋅h.

Twierdzenie to gwarantuje więc, że pochodną w kierunku wektora w można aproksymować ilorazem różnicowym

F′⁢x⁢w≈F⁢x+h⁢w-F⁢xh

z błędem rzędu O⁢h.

Dowód

Twierdzenie jest prostą konsekwencją twierdzenia o wartości średniej i założenia lipschitzowskości pochodnej:

	1h⁢F⁢x+h⁢w-F⁢x	=1h∫01F′(x+thw)⋅hwdt=∫01F′(x)w+(F′(x+thw)-F′(x))wdt
		=F′⁢x⁢w+∫01F′⁢x+t⁢h⁢w-F′⁢x⁢w⁢⁢d⁢t.

Stąd

1h⁢F⁢x+h⁢w-F⁢x-F′⁢x⁢w≤∫01F′⁢x+t⁢h⁢w-F′⁢x⁢w⁢⁢d⁢t≤∫01t⁢h⁢w⁢w⁢⁢d⁢t=L2⁢w2⋅h.

∎

Z powyższego wynika pomysł na przybliżenie całej macierzy pochodnej ilorazami różnicowymi: wystarczy wziąć F′⁢xk≈Ak, gdzie j-ta kolumna Ak, Akj, jest wyznaczona jako iloraz różnicowy

Aj=1h⁢F⁢x+h⁢ej-F⁢x≈F′⁢x⁢ej.

Twierdzenie 10.5 (o zbieżności metody z pochodną przybliżoną ilorazami różnicowymi)

Przy standardowych założeniach, istnieją dodatnie δ i h (dostatecznie małe) takie, że jeśli hk jest ciągiem takim, że 0<hk≤h oraz x0∈K⁢x*,δ, to ciąg zdefiniowany wzorem

xk+1=xk-Ak-1⁢F⁢xk,

w którym kolumny macierzy Ak zadane są ilorazami różnicowymi:

Akj=1h⁢F⁢xk+hk⁢ej-F⁢xk,⁢j=1,…,N,

jest dobrze określony i zbieżny przynajmniej liniowo do x*. Co więcej,

jeśli hk→0, to xk→x* superliniowo;
jeśli hk≤M⁢xk-x* dla pewnej stałej M, to xk→x* przynajmniej kwadratowo;
jeśli |hk|≤M∥F(xk∥ dla pewnej stałej M, to xk→x* przynajmniej kwadratowo.

Dowód

Dla wygody, dowód poprowadzimy w normie ∥⋅∥1 — tzw. kolumnowej normie macierzowej (w przestrzeni skończenie wymiarowej i tak wszystkie normy są równoważne). Mamy więc

Ak1=maxj⁡Akj1,

oraz

Ak-F′⁢xk1=maxj⁡Akj-F′⁢xk⁢ej1.

Tymczasem, na mocy lematu o różnicowej aproksymacji pochodnej,

Akj-F′⁢xk⁢ej1≤L2⁢ej12⁢hk=L2⁢hk.

Oznaczmy Δk=Ak-F′⁢xk, wtedy Δk1≤L2⁢hk. Na mocy twierdzenia o lokalnej stabilności metody Newtona z (dla δ,h dostatecznie małych) mamy, że Ak jest nieosobliwa oraz zachodzi

ek+11≤C⁢ek12+hk⁢ek1.

Stąd już wynika teza twierdzenia. Ostatni warunek kwadratowej zbieżności wynika z oszacowania residuum przez błąd,

F⁢xk1≤2⁢F′⁢x*1⁢ek1,

gdy tylko xk jest dostatecznie blisko x*.

∎

10.2.4. Niedokładna metoda Newtona

Jak pamiętamy, najkosztowniejszą częścią jednej iteracji Newtona jest rozwiązywanie układu równań z macierzą pochodnej. Alternatywą dla metody bezpośredniej rozwiązywania równania poprawki

F′⁢xk⁢s=F⁢xk

mogłaby być, zwłaszcza w przypadku gdy N jest duże, metoda iteracyjna (mielibyśmy zatem do czynienia z iteracją wewnątrz iteracji). Na k-tym kroku metody Newtona moglibyśmy zatrzymywać wewnętrzną iterację stosując np. residualne kryterium stopu z parametrem wymuszającym ηk,

F⁢xk-F′⁢xk⁢s≤ηk⁢F⁢xk.

Taką modyfikację metody Newtona nazywa się niedokładną metodą Newtona.

Twierdzenie 10.6 (o błędzie niedokładnej metody Newtona)

Przy standardowych założeniach, istnieją dodatnie stałe C i δ takie, że jeśli xk∈K⁢x*,δ oraz s spełnia warunek

F′⁢xk⁢s-F⁢xk≤ηk⁢F⁢xk,

to następna iteracja niedokładnej metody Newtona dana wzorem

xk+1=xk-s

spełnia

ek+1≤C⁢ek+ηk⁢ek.

Dowód

Niech δ będzie na tyle małe, by zachodził lemat 9.2 o oszacowaniu funkcji i pochodnej, a także by zachodziło twierdzenie 9.3 o zbieżności metody Newtona z punktem startowym xk. Aby udowodnić twierdzenie, najpierw postaramy się wyłuskać związek pomiędzy naszą metodą, a prawdziwą metodą Newtona (o której już sporo wiemy). Mamy bowiem

xk+1=xk-s=xk-F′⁢xk-1⁢F⁢xk+F′⁢xk-1⁢F⁢xk-s=xk+1Newton+F′⁢xk-1⁢r,

gdzie r jest residuum rozwiązania równania na poprawkę,

r=F⁢xk-F′⁢xk⁢s.

Stąd wynika, że

xk+1-x*=xk+1Newton-x*+F′⁢xk-1⁢r

i konsekwentnie

xk+1-x*≤xk+1Newton-x*+F′⁢xk-1⁢r.

Ponieważ na mocy twierdzenia o zbieżności metody Newtona mamy xk+1Newton-x*≤C~⁢xk-x*2, to wystarczy oszacować ostatni człon nierówności. Z definicji r mamy

F′⁢xk-1⁢r≤F′⁢xk-1⁢F⁢xk-F′⁢xk⁢s≤2⁢F′⁢x*-1⋅ηk⁢F⁢xk≤4⁢cond⁡F′⁢x*⁡ηk⁢xk-x*.

(W ostatniej nierówności skorzystaliśmy z lematu o oszacowaniu funkcji i pochodnej.)

∎

W powyższym dowodzie dostaliśmy rezultat, w którym stała C w oszacowaniu błędu zależała od uwarunkowania F′⁢x*, skąd moglibyśmy wysnuć wniosek, że jeśli F′⁢x* jest źle uwarunkowana, to ηk należy brać patologicznie małe. W rzeczywistości nie jest aż tak źle i można osłabić założenia na ciąg ηk, ale pod warunkiem zmiany normy, w której badamy błąd, tak, by uwzględniać zachowanie się pochodnej: ∥⋅∥=∥F′(x*)⋅∥ (zob. [9]).

Wniosek 10.1 (twierdzenie o zbieżności niedokładnej metody Newtona)

Przy standardowych założeniach, istnieją δ i η (dostatecznie małe) takie, że jeśli x0∈K⁢x*,δ oraz 0≤ηk≤η to niedokładna metoda Newtona z parametrem wymuszającym ηk jest zbieżna przynajmniej liniowo do x*. Ponadto,

jeśli ηk→0, to zbieżność jest superliniowa;
jeśli ηk≤Mη⁢F⁢xk dla pewnego ustalonego Mη, to zbieżność jest przynajmniej kwadratowa.

Ćwiczenie 10.6

Udowodnij powyższy wniosek.

Niedokładna metoda Newtona jest bardzo popularną metodą w przypadku, gdy N jest bardzo duże: wówczas trudno myśleć o sposobach rozwiązywania równania poprawki innych niż jakaś metoda iteracyjna. Dodatkowym plusem tej metody jest to, że znakomicie można ją łączyć z przybliżaniem pochodnej kierunkowej ilorazami różnicowymi; rzeczywiście, głównym składnikiem metody iteracyjnej jest mnożenie F′⁢xk⁢w, co można tanio przybliżać w sposób opisany w rozdziale 10.2.3. W ten sposób oszczędzamy na kilku polach na raz:

nie musimy wykonywać, zazwyczaj bardzo kosztownego, rozkładu macierzy pochodnej;
nie musimy nawet wyznaczać macierzy pochodnej, lecz w zamian wystarczy, że w każdej wewnętrznej iteracji, jedną dodatkową wartość F⁢xk+h⁢w;
ponieważ po kilku iteracjach zbliżamy się do rozwiązania, że xk+1≈xk, w ten sposób od razu dysponujemy dobrym przybliżeniem startowym dla iteracji wewnętrznej: w sprzyjających okolicznościach metody Kryłowa mogą z tego zrobić dodatkowy pożytek;
możemy ograniczyć koszt rozwiązywania układu z macierzą pochodnej (liczbę iteracji wewnętrznych), stosując łagodne kryterium stopu, gdy jesteśmy daleko od rozwiązania.

Przykład 10.1 (Numeryczne eksperymenty z wariantami metody Newtona dla równania Allena–Cahna)

Porównamy jakość pracy kilku metod:

klasycznej metody Newtona
uproszczonej metody Newtona
metody Szamańskiego o m=4 krokach
metody Newtona z pochodną przybliżoną ilorazami różnicowymi (ze stałym, ale bardzo małym krokiem, h≈ϵ, gdzie ϵ to precyzja arytmetyki.

Do porównania wybierzemy:

wykresy uzyskanego rozwiązania, u⁢x, dla każdej z metod
wykresy residuum, F⁢u⁢x, dla każdej z metod
historię zbieżności metody, F⁢uj2
liczbę iteracji, czas pracy, końcową wartość residuum.

%
% 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)
%
% Rozwiazuje uklad rownan nieliniowych F(x) = 0, gdzie F jest zadana M-funkcja o 
% nazwie ”nazwa_f”.
% Parametry wejscia:
% 	nazwa_f - nazwa M-funkcji zadajacej funkcje F(). 
% 		Ta funkcja *musi* przyjmowac takze argumenty macierzowe!
% 	x0 - przyblizenie poczatkowe rozwiazania
%	atol, rtol - tolerancja bledu residualnego
%	step - co ile iteracji uaktualniac macierz pochodnej; step = 1 daje standardowa metode Newtona, step > maxit daje metode cieciw
%	maxit - maksymalna dopuszczalna liczba iteracji
% Wyjscie:
%	x - koncowe przyblizenie rozwiazania
%	resid - wektor zawierajacy historie zbieznosci, resid(i) - wartosc
% 		normy residuum w i-tej iteracji
%	info - kod zakonczenia: = 1 - sukces, = 0 - przekroczona maxit
%	output - struktura zawierajaca statystyki

% kilka pomocniczych definicji
epsq = sqrt(eps);	% pierwiastek z podwojonej precyzji maszyny
unit = ones(size(x0));	% wektor zlozony z samych jedynek

% poczatkowe rozwiazanie przyblizone
x = x0;
	
% wyznaczenie residuum poczatkowego
fx = feval(nazwa_f, x);
iter = 1;
resid(iter) = norm(fx);

while( (resid(iter) >= (atol + rtol*resid(1))) && (iter < maxit) )

if (mod(iter-1, step) == 0)
		% aproksymacja pochodnej - wszystkie kolumny jednoczesnie
		h = epsq*max(abs(x),unit);
		xh = x*unit' + diag(h);
		if exist(nazwa_df) > 0 
			Df = feval(nazwa_df, x);
		else
			Df = (feval(nazwa_f, xh) - fx*unit')./(unit*h');
		end
		% rozklad macierzy pochodnej do wykorzystania
		[L, U, p] = lu(Df, 'vector');
	end
	% poprawka przyblizonego rozwiazania, x = x - Df\fx;
	x = x - U\(L\fx(p));
	
	% wyznaczenie residuum
	fx = feval(nazwa_f, x);
	iter = iter + 1;
	resid(iter) = norm(fx);
end
info = ~(iter == maxit);
output.iterations = iter - 1;
end

function y = acfun(u) % lewa strona rownania Allena-Cahna
global T; % macierz laplasjanu
global delta;
global x; % wektor wezlow zmiennej przestrzennej
y = T*u + delta*u.*(1-u.^2) - repmat(delta*sqrt(x),[1, size(u,2)]);
end

function y = acdfun(u) % pochodna acfun
global T; % macierz laplasjanu
global delta;
y = T + delta*spdiag(1-3*u.^2);
end

global T;
global delta;
global x;

atol = 1e-8; rtol = 0; maxit = 300; delta = 10;
Nrange = [250]; % zestaw wartosci testowych rozdzielczosci dyskretyzacji
for k = 1:length(Nrange)
	N = Nrange(k); % rozmiar macierzy
	x = linspace(0,1,N+2)'; x = x(2:end-1);
	e = ((N+1)^2)*ones(N,1); T = spdiags([-e, 2*e, -e], [-1,0,1], N,N);
	u0 = ones(N,1);

% metoda Newtona
	tic; [un, residn] = newton('acfun', 'acdfun', u0, atol, rtol, 1, maxit);
	timen = toc;
	% historia zbieznosci metody 
	semilogy(residn,'b-o;|F(u_{Newton})|;'); pause(1);

% uproszczona metoda Newtona
	tic; [uun, residun] = newton('acfun', 'acdfun', u0, atol, rtol, maxit+1, maxit);
	timeun = toc;
	% historia zbieznosci metody 
	semilogy(residun,'r-o;|F(u_{uproszcz Newton})|;'); pause(1);
	
	% metoda Szamanskiego z odswiezeniem co 4 iteracje
	tic; [usz, residsz] = newton('acfun', 'acdfun', u0, atol, rtol, 4, maxit);
	timesz = toc;
	% historia zbieznosci metody 
	semilogy(residsz,'g-o;|F(u_{Szamanski})|;'); pause(1);
	
	% metoda Newtona z przyblizona pochodna
	tic; [unfd, residnfd] = newton('acfun', ”, u0, atol, rtol, 1, maxit);
	timenfd = toc;
	% historia zbieznosci metody 
	semilogy(residnfd,'m-o;|F(u_{Newton przybl poch})|;'); pause(1);

% rozwiazanie
	plot(x,[un,uun,usz,unfd], ['b-;u_{Newton};';'r-;u_{Newton upr};';'g-;u_{Szamanski};';'m-;u_{Newton przybl poch};']); pause(1);
	
	% residuum
	semilogy(x,abs(acfun([un,uun,usz,unfd])), ['b-;|F(u_{Newton})|;';'r-;|F(u_{Newton upr})|;';'g-;|F(u_{Szamanski})|;';'m-;|F(u_{Newton przybl poch})|;']); pause(1);
	
	
	[length(residn), length(residun), length(residsz), length(residnfd)]
	[timen,timeun,timesz,timenfd]
	[residn(end), residun(end), residsz(end), residnfd(end)]
end

Zwróć uwagę na to, że odpowiednio dobierając częstość uaktualniania macierzy pochodnej, można znacząco poprawić efektywność metody (w porównaniu do metody Newtona). Z drugiej strony, zbyt rzadka aktualizacja może przeszkodzić w uzyskaniu zbieżności, jak to dzieje się w naszym przykładzie w przypadku uproszczonej metody Newtona.

Ćwiczenie 10.7

Sprawdź, jak zmienią się wyniki (uzyskane rozwiązanie, jego dokładność, a także — koszt metody i szybkość jej zbieżności), gdy zmienisz jeden z poniższych parametrów

osłabisz nieliniowość, biorąc δ=1,
nieco wzmocnisz nieliniowość, biorąc δ=15,
zwiększysz rozmiar zadania N do 400,
zmniejszysz N do 20,
zmniejszysz tolerancję błędu do 10-12.

Wskazówka:

Zdaje się, że dla N=20 mamy niejednoznaczność rozwiązania? Jak to zweryfikować?

Ćwiczenie 10.8

Zaimplementuj niedokładną metodę Newtona i wykorzystaj ją w skrypcie rozwiązującym (najlepiej: dwuwymiarowe) równanie Allena–Cahna. Jak wpłynęło to na efektywność metody, w porównaniu do standardowej metody Newtona?

Przykład 10.2

Monografii Kelley'a towarzyszy zestaw skryptów, pozwalających przetestować działanie różnych metod iteracyjnych w kilku praktycznych zadaniach [9].

Poniższy skrypt, którego podstawowe kody źródłowe pochodzą ze strony http://www4.ncsu.edu/~ctk/newton, umożliwi Ci zapoznanie się z zadaniem rozwiązania równania Chandrasekhara [9, rozdział 5.6].

Aby w pełni wykorzystać możliwości skryptu, zachęcamy do uruchomienia go na własnym (podłączonym do internetu) komputerze i następnie do wnikliwej obserwacji zmiany wyników przy zmianie, czy to parametrów zadania, czy to parametrów pracy solvera.

Treść automatycznie generowana z plików źródłowych LaTeXa za pomocą oprogramowania wykorzystującego LaTeXML.