Równania nieliniowe

Jak wiadomo z podstawowego kursu analizy numerycznej, zazwyczaj nie istnieje bezpośredni algorytm rozwiązania równania nieliniowego postaci

,

gdzie . Do rozwiązywania równań powyższej postaci używa się metod iteracyjnych, takich jak metoda iteracji prostej Banacha:

,

czy też metoda Newtona

gdzie jest macierzą pochodnej (zwanej także jakobianem) odwzorowania F w punkcie x_n. Czytelnika zainteresowanego, jakie założenia na funkcję F oraz przybliżenie początkowe x₀ gwarantują zbieżność powyższych procesów iteracyjnych, odsyłamy do podręczników analizy numerycznej [6], [18], oraz licznych monografii, np. [15], [11].

Zauważmy, że w przypadku iteracji Newtona nie ma potrzeby (ani sensu) odwracania macierzy F'(x_n) dla obliczenia x_n+1. Kolejne przybliżenie rozwiązania x^* uzyskamy najpierw wyznaczając poprawkę na x_n , rozwiązując względem układ z macierzą jakobianu, a następnie kładąc . Oczywiście, ponieważ macierz pochodnej zmienia się w każdym kroku, rozwiązanie jednego układu liniowego jest dużo tańsze niż najpierw wyznaczenie macierzy odwrotnej, a dopiero potem pomnożenie jej przez wektor F(x_n).

Reasumując, aby znaleźć nowe przybliżenie, musimy skonstruować macierz pochodnej, a następnie rozwiązać układ liniowy z tą macierzą. Jeśli do rozwiązywania tego układu użyjemy metody iteracyjnej (co może być uzasadnione w przypadku, gdy mamy do czynienia z macierzą rzadką dużego rozmiaru), taka metoda będzie wtedy tzw. niedokładną (ang. inexact) metodą Newtona. Możliwe są też inne modyfikacje metody Newtona, o mniejszym koszcie jednego kroku (często właśnie konstrukcja macierzy jakobianu jest zaskakująco kosztownym etapem metody Newtona). W takich przypadkach czasem stosuje się wariant, w którym macierz jakobianu jest uaktualniana dopiero po wykonaniu kilku kroków metody. Innym rozwiązaniem jest aproksymacja macierzy jakobianu przy pomocy różnic skończonych (wielowymiarowy analog metody siecznych). Zainteresowanych przeglądem rozmaitych metod rozwiązywania układów równań nieliniowych odsyłamy do zwięzłej książki [11] oraz do monografii [15].

Solvery nieliniowe dostępne w PETSc dotyczą wyłącznie metod typu Newtona. Nic dziwnego: implementacja metody iteracji (nomen omen) prostej jest co prawda nieskomplikowana, ale jej szybkość zbieżności pozostawia dużo do życzenia. Natomiast metoda Newtona jest zbieżna kwadratowo, a ponieważ w każdym kroku algorytmu musimy rozwiązać układ równań liniowych, będziemy mieli doskonałą okazję wykorzystania pełnej mocy i elastyczności SLESu.

Na niekorzyść algorytmu Newtona zdaje się działać fakt, że wymaga on dobrego przybliżenia początkowego. Na szczęście, tę wadę metody Newtona często można obejść w przypadku dyskretyzacji równań różniczkowych cząstkowych, gdy "dobre'' przybliżenie początkowe można zgadnąć, lub wyliczyć w mniej lub bardziej wyrafinowany sposób. Ponadto, w PETSc zastosowano heurystyczne techniki "globalizujące'' zbieżność metody Newtona. Bazują one na spostrzeżeniu, że w przypadku, gdy startujemy ze złego przybliżenia, poprawka jest co prawda wektorem skierowanym w dobrym kierunku (tzn. w stronę rozwiązania x*), ale za dużej długości. Taka sytuacja prowadzi do "przestrzelenia'' rozwiązania i w efekcie do niezbieżności całego procesu. Dlatego wprowadza się dodatkową operację poszukiwania optymalnej długości wektora poprawki (tzw. line search methods). Więcej informacji na ten temat znajduje się w [11], rozdział 8.