Zagadnienia

9.1 Szybkość zbieżności
9.2 Metoda Banacha
9.3 Metoda Newtona
- 9.3.1 Analiza zbieżności metody Newtona
  - 9.3.1.1 Lematy techniczne
- 9.3.2 Implementacja metody Newtona

9. Układy równań nieliniowych. Metoda Newtona

Zajmiemy się teraz nową, trudniejszą klasą zadań: metodami rozwiązywania układów równań nieliniowych

F⁢x=0,

gdzie F:RN⊃D→RN.

W matematyce stosowanej bardzo wiele modeli zyskuje na realizmie, gdy uwzględnić w nich zjawiska nieliniowe. Takie modele są jednak znacznie trudniejsze — zarówno w analizie, jak i w komputerowych symulacjach. Nieustanny wzrost mocy obliczeniowej komputerów, jaki obserwujemy w ciągu ostatnich kilkudziesięciu lat powoduje, że rozwiązywanie dużych układów nawet silnie nieliniowych równań stało się możliwe nawet dla posiadaczy ,,zwykłych” komputerów osobistych.

Od razu zauważmy, że już na poziomie matematycznej poprawności tego zadania możemy napotkać kłopoty, albowiem to zadanie może nie mieć rozwiązań, może mieć dokładnie jedno rozwiązanie, albo skończenie wiele rozwiązań, albo nieskończenie wiele izolowanych rozwiązań, albo… i tak dalej, i tak dalej.

Aby więc szukać jakiegokolwiek rozwiązania powyższego równania, musimy najpierw mieć pewność, że ono istnieje. Nie będziemy się tą kwestią tu zajmować, traktując ją jako domenę matematyki czystej, bez numerycznego zacięcia.

Z kursu Matematyki obliczeniowej I znamy kilka metod rozwiązywania nieliniowego równania skalarnego, f⁢x=0: metodę bisekcji, stycznych, siecznych, …. Okazuje się, że niektóre z nich da się w mniej lub bardziej naturalny sposób uogólnić na przypadek wielowymiarowy. W niniejszym rozdziale wyjaśnimy, jak metodę stycznych (Newtona) uogólnić na przypadek wielowymiarowy.

9.1. Szybkość zbieżności

Rozważmy ciąg xk⊂RN taki, że xk→x*∈RN gdy k→∞.

Definicja 9.1 (Zbieżność wykładnicza)

Powiemy, że xk→x* wykładniczo z rzędem co najmniej q>1, jeśli

∃C≥0⁢⁢∃K∈N⁢⁢xk+1-x*≤C⁢xk-x*q⁢⁢∀k≥K.

(9.1)

Jeśli wykładnik zbieżności q=2, mówimy wtedy o zbieżności kwadratowej; gdy q=3, zbieżność nazywamy sześcienną (albo: kubiczną). W praktyce, zbieżność kwadratową uznaje się za satysfakcjonująco szybką, a zbieżność kubiczną — za bardzo szybką. Jednak nie zawsze taką szybkość zbieżności można uzyskać; rozsądnym praktycznym minimum jest zbieżność liniowa, w której błąd w każdym kroku maleje o stały czynnik:

Definicja 9.2 (Zbieżność liniowa)

Powiemy, że xk→x* przynajmniej liniowo, jeśli

∃0≤α≤1⁢⁢∃K∈N⁢⁢xk+1-x*≤α⁢xk-x*⁢⁢∀k≥K.

Jednak gdy α≈1, zbieżność liniowa może być w praktyce nieakceptowalnie wolna.

Oczywiście powyższe nie wyczerpuje wszystkich możliwych szybkości zbieżności. Mianowicie, ciąg może być zbieżny szybciej niż liniowo (choć niekoniecznie zaraz kwadratowo):

Definicja 9.3 (Zbieżność superliniowa)

Powiemy, że xk→x* superliniowo, jeśli

lim supk→∞⁡xk+1-x*xk-x*=0.

Czasem trudno o samym ciągu powiedzieć, jaki ma wykładnik zbieżności, natomiast łatwiej jest wskazać ciąg o znanym wykładniku zbieżności, który majoryzuje błąd popełniany na każdej iteracji.

Definicja 9.4 (r–zbieżność)

Powiemy, że xk→x* r–wykładniczo/r–liniowo/r–superliniowo, jeśli istnieje ciąg ξk⊂R taki, że

xk-x*≤ξk

oraz ξk→0, odpowiednio, wykładniczo/liniowo/superliniowo.

Aby móc porównywać (jakościowo) różne metody iteracyjne rozwiązywania równań nieliniowych — o różnych kosztach jednej iteracji i o różnych wykładnikach zbieżności — wprowadza się pojęcie wskaźnika efektywności metody.

Definicja 9.5 (Efektywność metody iteracyjnej)

Wskaźnik efektywności E metody iteracyjnej, zbieżnej wykładniczo z wykładnikiem q, której jedna iteracja kosztuje m flopów, definiujemy jako [17]

E=q1/m.

Wskaźnik efektywności metody służy szybkiemu jakościowemu porównywaniu metod. Motywacją dla takiej definicji jest wymaganie, by przy zadanym (tym samym) całkowitym koszcie iteracji W obu metod, porównać uzyskane wykładniki redukcji błędu. Jeśli bowiem metoda jest rzędu q>1, to błąd początkowy po k iteracjach zostanie zredukowany do poziomu x0-x*qk. Ponieważ całkowity koszt k iteracji wynosi W=k⁢m, to qk=qW/m=q1/mW=EW. Definiując w ten sposób efektywność metody iteracyjnej, czynimy jednak pewne oproszczenie, polegające na tym, że pomijamy wpływ stałych na redukcję błędu (patrz oszacowanie (9.1)). Z tego powodu na przykład (formalnie) każda metoda zbieżna liniowo formalnie miałaby taką samą efektywność — równą 1 — choć oczywiście realna sprawność różnych metod liniowych może być różna, m.in. właśnie w zależności od stałej redukcji błędu.

Rozważmy dwa przykłady dużych układów równań nieliniowych, motywowanych — podobnie jak już kilka razy wcześniej w trakcie tego przedmiotu — zadaniami z równań różniczkowych cząstkowych.

Przykład 9.1 (stacjonarne równanie Allena–Cahna)

Rozważmy równanie

		-Δ⁢u+δ⁢u⁢1-u2=g⁢x,y		w ⁢Ω,
		u=0		na ⁢∂⁡Ω.

Dla uproszczenia przyjmiemy, że Ω jest kwadratem jednostkowym. Wprowadzając dyskretyzację tego równania najprostszą¹¹Po raz kolejny przypomnijmy, że (z wyjątkiem trywialnych przypadków) najprostsze metody dyskretyzacji, choć bardzo ułatwiają programowanie, to jednak są zwykle najmniej efektywnym sposobem rozwiązywania równań różniczkowych, jeśli wziąć pod uwagę koszt obliczeniowy potrzebny do uzyskania satysfakcjonującego przybliżenia rozwiązania. metodą — metodą różnic skończonych na jednorodnej siatce — patrz przykład 5.2 — otrzymujemy, po uwzględnieniu warunku brzegowego, nieliniowy układ równań algebraicznych,

PN⁢UN+f⁢UN=0,

gdzie PN jest operatorem liniowym zdefiniowanym w przykładzie 5.2, natomiast dla zadanej macierzy U∈RN×N

f⁢Ui⁢j=δ⁢ui⁢j⁢1-ui⁢j2-gi⁢j.

(Przypomnijmy, że dokonujemy tutaj utożsamienia macierzy N×N z wektorem długości N2, ktorego elementami są kolejne wiersze macierzy.)

Przykład 9.2 (ewolucyjne równanie Allena–Cahna)

Rozważmy równanie

		ut-Δ⁢u+δ⁢u⁢1-u2=g⁢x,y		w ⁢Ω,
		u=0		na ⁢∂⁡Ω.

Dla uproszczenia przyjmiemy, że Ω jest kwadratem jednostkowym. Wprowadzając dyskretyzację tego równania najprostszą metodą — metodą różnic skończonych na jednorodnej siatce po zmiennej przestrzennej oraz niejawnym schematem Eulera po zmiennej czasowej (patrz wykład z Numerycznych równań różniczkowych, otrzymujemy, po uwzględnieniu warunku brzegowego, nieliniowy układ równań algebraicznych z niewiadomą UNk,

UNk-UNk-1τ+PN⁢UNk+f⁢UNk=0,

gdzie PN jest operatorem liniowym zdefiniowanym w przykładzie 5.2, a f jest zadane jak w poprzednim przykładzie. (Przypomnijmy, że dokonujemy tutaj utożsamienia macierzy N×N z wektorem długości N2, ktorego elementami są kolejne wiersze macierzy.) Zostawiając po lewej stronie tylko wyrazy z niewiadomymi, dostajemy równanie

I+τ⁢PN⁢UNk+τ⁢f⁢UNk=UNk-1,

a więc takie, w którym — dla dostatecznie małych τ — część nieliniowa staje się mała w porównaniu z UN.

9.2. Metoda Banacha

Wiele równań nieliniowych można w naturalny sposób sformułować jako zagadnienie punktu stałego,

x=Φ⁢x,

dla którego najprostszą metodą iteracyjną jest doskonale nam znana metoda Banacha (iteracji prostej),

xk+1=Φ⁢xk.

Przypomnijmy (bez dowodu) twierdzenie o zbieżności tej metody:

Twierdzenie 9.1 (o zbieżności metody Banacha)

Niech D⊂RN będzie niepusty, ograniczony i domknięty, i niech Φ:D→D będzie kontrakcją, tzn.

∃γ<1⁢⁢Φ⁢x-Φ⁢y≤γ⁢x-y⁢⁢∀x,y∈D.

Wtedy Φ ma dokładnie jeden punkt stały x*=Φ⁢x*∈D oraz dla dowolnego x0∈D iteracja

xk+1=Φ⁢xk

jest zbieżna do x*, przy czym dla każdego k=0,1,…,

xk+1-x*≤γ⁢xk-x*.

Z powyższego twierdzenia wynika więc, że iteracja Banacha jest zbieżna przynajmniej liniowo ze współczynnikiem redukcji błędu γ<1.

Przykład 9.3

Stosując niejawny schemat Eulera¹²Por. np. wykład z Numerycznych równań różniczkowych. Podobna metoda — pod nazwą metody łamanych Eulera — bywa wykorzystywana w dowodzie twierdzenia o istnieniu rozwiązania równania różniczkowego. Nas tutaj interesuje ona wyłącznie jako metoda numeryczna.

yn+1=yn+τ⁢F⁢tn+1,yn+1

do przybliżonego rozwiązania układu równań różniczkowych zwyczajnych y′⁢t=F⁢t,y∈RN musimy dla znanego yn znaleźć yn+1 jako punkt stały odwzorowania

Φ⁢x=yn+τ⁢F⁢tn+1,x.

Jeśli F jest lipschitzowska ze stałą L w normie ∥⋅∥ ze względu na drugi argument (por. twierdzenie Picarda–Lindelöfa), to

Φ⁢x-Φ⁢y=τ⁢F⁢x-F⁢y≤τ⁢L⁢x-y⁢⁢∀x,y∈RN,

zatem Φ jet kontrakcją na RN, gdy τ⁢L<1, czyli gdy krok czasowy τ schematu Eulera jest dostatecznie mały względem h.

Czasami wygodniej jest skorzystać z ,,lokalnej” wersji twierdzenia o zbieżności metody Banacha:

Twierdzenie 9.2 (o lokalnej zbieżności metody Banacha)

Niech Φ:RN⊃D→RN, gdzie D jest otwarty i niepusty, i niech dla zadanego x0∈D istnieje kula B=x∈RN:x-x0≤r taka, że B⊂D oraz Φ jest na B kontrakcją ze stałą γ.

Jeśli dodatkowo x0-Φ⁢x0≤1-γ⁢r, to w B istnieje dokładnie jeden punkt stały x*∈B odwzorowania Φ, a metoda iteracji prostej jest doń zbieżna przynajmniej liniowo ze stałą redukcji błędu γ.

Dowód

Wystarczy wykazać, że Φ⁢B⊂B i skorzystać z poprzedniego twierdzenia. Na mocy warunku kontrakcji w B, dla x∈B mamy

Φ⁢x-x0≤Φ⁢x-Φ⁢x0+Φ⁢x0-x0≤γ⁢x-x0+1-γ⁢r≤γ⁢r+1-γ⁢r=r,

zatem faktycznie Φ⁢x∈B.

∎

9.3. Metoda Newtona

W przypadku równania skalarnego, f:R→R, metoda stycznych (zwana też metodą Newtona) rozwiązywania równania f⁢x=0 jest zadana wzorem

xk+1=xk-f⁢xkf′⁢xk.

Przez analogię, gdy F:RN⊃D→RN można więc byłoby zdefiniować wielowymiarową metodę Newtona wzorem

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

(9.2)

gdzie F′⁢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:

Zbiór D jest otwarty i niepusty, a F jest różniczkowalna w D.
Istnieje rozwiązanie x*∈D:

F⁢x*=0.
Pochodna F′:D→L⁢RN jest lipschitzowska ze stałą L:

∃L>0⁢⁢F′⁢x-F′⁢y≤L⁢x-y⁢⁢∀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,

A=supx≠0⁡A⁢xx.
Macierz pochodnej w rozwiązaniu, F′⁢x*, jest nieosobliwa.

9.3.1.1. Lematy techniczne

Przez K⁢x*,δ będziemy oznaczali kulę otwartą o środku w x* i promieniu δ,

K⁢x*,δ=x∈RN:x-x*<δ.

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=∫01F′⁢x+t⁢y-x⁢y-x⁢⁢d⁢t.

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 F′(⋅) jest ciągła na K¯, zatem także g′⁢t=F′⁢x+t⁢y-x⁢y-x jest ciągła na 0,1 (i całkowalna). Teza lematu wynika więc z podstawowego twierdzenia rachunku różniczkowego, że

∫01g′⁢t⁢⁢d⁢t=g⁢1-g⁢0.

∎

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:

F′⁢x≤2⁢F′⁢x*,
F′⁢x-1≤2⁢F′⁢x*-1,
12⁢F′⁢x*-1⁢x-x*≤F⁢x≤2⁢F′⁢x*⁢x-x*,

Ćwiczenie 9.1

Wykaż, że przy założeniach standardowych x* musi być izolowanym rozwiązaniem równania F⁢x=0.

Rozwiązanie:

Z lematu 9.2 o oszacowaniu funkcji, istnieje otoczenie U∋x* takie, że

12⁢F′⁢x*-1⁢x-x*≤F⁢x.

Jeśli więc dla x~ leżącego w tym otoczeniu F⁢x~=0, to znaczyłoby to, że x~-x*≤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:

xk+1-x*≤C⁢xk-x*2⁢⁢∀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-F′⁢xk-1⁢F⁢xk.

Na mocy założeń standardowych i lematu o wartości średniej,

F⁢xk=F⁢xk-F⁢x*=∫01F′⁢x*+t⁢xk-x*⁢xk-x*⁢⁢d⁢t=∫01F′⁢x*+t⁢ek⁢ek⁢⁢d⁢t,

zatem

ek+1=ek-F′⁢xk-1⁢∫01F′⁢x*+t⁢ek⁢ek⁢⁢d⁢t=F′⁢xk-1⁢∫01F′⁢xk-F′⁢x*+t⁢ek⁢ek⁢⁢d⁢t.

Korzystając z lipschitzowskości pochodnej i raz jeszcze z lematu 9.2 o oszacowaniu pochodnej, dostajemy

	ek+1	≤F′⁢xk-1⁢∫01F′⁢xk-F′⁢x*+t⁢ek⁢⁢ek⁢⁢d⁢t
		≤2∥F′(x*)-1∥L∥ek∥2∫01(1-t)dt=:C∥ek∥2.

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

ek+1≤C⁢ek2=C⁢ek⁢⁢ek≤C⁢δ⁢ek,

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 lipschitzowskości pochodnej warunkiem

∃α∈0,1⁢⁢∃L>0⁢⁢F′⁢x-F′⁢y≤L⁢x-yα⁢⁢∀x,y∈D,

(czyli założyć, że F′ jest w swej dziedzinie hölderowska 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, modyfikując oszacowanie F′⁢xk-F′⁢x*+t⁢ek.

Ćwiczenie 9.3

Jak, przy założeniach standardowych, oszacować normę błędu, xk-x*, przez normę residuum, F⁢xk?

9.3.2. Implementacja metody Newtona

Implementując metodę Newtona, nie będziemy rzecz jasna nigdy explicite wyznaczali macierzy odwrotnej F′⁢xk-1, tylko rozwiązywali układ równań z macierzą F′⁢xk — tę jednak będziemy już musieli wyznaczyć. Pamiętając także o tym, że metoda Newtona jest metodą iteracyjną, 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 F′⁢x;

rozwiąż F′⁢x⁢s=F⁢x;

x = x - s;

end

return(x);

Ćwiczenie 9.4

Macierzowe zadanie własne dla symetrycznej macierzy A∈RN×N, znajdowania pary x,λ∈RN×R spełniającej

A⁢x=λ⁢x,⁢x≠0

można potraktować jako kwadratowe równanie nieliniowe dla F:RN+1→RN+1 danej na przykład wzorem

F⁢x,λ=A⁢x-λ⁢x1-12⁢xT⁢x.

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

F′⁢x,λ=A-λ⁢I-x-xT0.

Zatem metodę Newtona można — jeśli tylko F′⁢x,λ nieosobliwa — zaimplementować, korzystając na przykład ze wzoru Shermana–Morrisona (por. ćwiczenie 5.23).

Zauważmy, że F jest funkcją gładką, jej pochodna jest lipschitzowska ze stałą Lipschitza równą… no właśnie, ile? Mamy

F′⁢x,λ-F′⁢y,μ=μ-λ⁢Iy-xy-xT0,

zatem

	F′⁢x,λ-F′⁢y,μ1	=max⁡μ-λ+y-x1,…,y-x1
		≤∥y-x∥1+\|μ-λ\|=∥xλ-yμ∥1.

A więc — w normie indukowanej normą ∥⋅∥1 — F′ 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

QT1⁢F′⁢x,λ⁢Q1=Λ-λ⁢IQT⁢xxT⁢Q0.

Niech λ*=λi. Wtedy x*=C⁢qi 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 QT⁢x=C⁢ei, 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 macierzą nieosobliwą i rozważmy iterację

Xk+1=Xk+Xk⁢I-A⁢Xk,

startującą z macierzy kwadratowej X0. Wykaż, że jeśli I-A⁢X0<1, to Xk→A-1 gdy k→∞.

Wskazówka:

Wykaż, że I-A⁢Xk+1=I-A⁢Xk2, 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. przykład 9.1), że równanie nieliniowe, które nas interesuje w przypadku stacjonarnym ma postać

F⁢UN=PN⁢UN+f⁢UN,

a w przypadku ewolucyjnym —

F⁢UN=UN+τ⁢PN⁢UN+f⁢UN.

Oba przypadki możemy zatem ogarnąć, definiując funkcję

Fτ,α⁢UN=α⁢UN+τ⁢PN⁢UN+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τ,α′⁢UN. Ponieważ (traktując f jako funkcję na RN2)

f⁢UNj=δ⁢uj⁢1-uj2+gj,

∂⁡fj∂⁡uk⁢UN=1-3⁢uj2,k=j,0,k≠j,

zatem ostatecznie

Fτ,α′⁢UN⁢⁢V=α⁢I+τ⁢PN+diag⁢1-3⁢uj2⁢⁢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 rozwiązywania równania Allena–Cahna opisywanego w przykładzie 9.1, będziemy rozważać dyskretyzacje jednowymiarowego stacjonarnego równania Allena–Cahna,

TN⁢uj+δ⁢uj⁢1-uj2=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)
%
% 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 = [25, 250]; % zestaw wartosci testowych parametru p
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);

tic; [u, resid, info, output] = newton('acfun', 'acdfun', u0, atol, rtol, 1, maxit);
	time(k) = toc;

opts.TolFun = atol;
	opts.MaxIter = maxit;
	tic; [U, fU, Info, Output] = fsolve('acfun',u0, opts); timeg(k) = toc;  
	
	% rozwiazanie
	plot(x,[u,U], ['b-;u_{Newton};';'r-;u_{fsolve};']); pause(1);
	
	% residuum
	semilogy(x,abs([acfun(u),acfun(U)]), ['b-;|F(u_{Newton})|;';'r-;|F(u_{fsolve})|;']); pause(1);
	
	% historia zbieznosci metody Newtona
	semilogy(resid,'b-o;|F(u_{Newton})|;'); pause(1);
	
	[norm(acfun(u)),norm(acfun(U))]
end

Ćwiczenie 9.6

Przeprowadź podobne eksperymenty numeryczne dla równania Allena–Cahna w kwadracie, jak to opisywano we wcześniejszych przykładach.

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