Kiedy Newton-Krylov nie jest odpowiednim rozwiązaniem?

Ostatnio porównywałem różne nieliniowe solwery z scipy i byłem pod szczególnym wrażeniem przykładu Newtona-Kryłowa w książce kucharskiej Scipy, w którym rozwiązują równanie różniczkowe drugiego rzędu z nieliniowym terminem reakcji w około 20 wierszach kodu.

Zmodyfikowałem przykładowy kod, aby rozwiązać nieliniowe równanie Poissona ( zwane również równaniem Poissona-Boltzmanna , patrz strona 17 w tych uwagach) dla heterostruktur półprzewodnikowych, które mają postać:

\frac{{re}^{2)} ϕ}{re x^{2)}} - k (x) (p (x, ϕ) - n (x, ϕ) + {N.}^{+} (x)) = 0

$\frac{d^2\phi}{dx^2} - k(x) \left(p(x,\phi) - n(x,\phi) + N^{+}(x)\right) = 0$

(Jest to funkcja resztkowa przekazywana do solvera).

Jest to problem elektrostatyczny, w którym i są funkcjami nieliniowymi dla postaci . Szczegóły tutaj nie są ważne, ale chodzi o to, że funkcja nieliniowa zmienia się wykładniczo z więc funkcja resztkowa może zmieniać się w szerokim zakresie ( z niewielką zmianą w . $n(x,\phi)$ $p(x,\phi)$ $n_i(x) e^{-(E_i(x,\phi) - E_f)}$ $\phi$ $10^{-6} - 10^{16})$ $\phi$

Numerycznie rozwiązuję to równanie za pomocą Newtona-Kryłowa Scipy'ego, ale nigdy by się nie zbiegło (w rzeczywistości zawsze zgłaszałoby błąd przy obliczaniu jakobianu). Przeszedłem z Newton-Kryłowa solver do fsolve (który jest oparty na MINPACK hybrd) i pracował po raz pierwszy!

Czy istnieją ogólne powody, dla których Newton-Krylov nie jest odpowiedni do niektórych problemów? Czy równania wejściowe muszą być w jakiś sposób uwarunkowane?

Być może potrzeba więcej informacji, aby skomentować, ale dlaczego według ciebie fsolve działało w tej sprawie?

— boyfarrell
źródło

Miałem ten sam problem z niepowodzeniem Newtona-Kryłowa z Jakubem i stwierdziłem, że zmiana metody z „lgmres” na „gmres” ( sol = newton_krylov(func, guess, method='gmres')) rozwiązała problem. Nie jestem pewien, dlaczego, ale każdy z tym problemem może rozważyć zrobienie tego samego.

— Arthur Dent

Istnieją dwa problemy, które możesz napotkać.

Źle warunkowane

Po pierwsze, problem jest źle uwarunkowany, ale jeśli podasz resztkę, Newton-Kryłow wyrzuca połowę twoich znaczących cyfr, różnicując skończoną różnicę, aby uzyskać działanie jakobianów:

jot [x] y \approx \frac{fa (x + ϵ y) - fa (x)}{ϵ}

$J[x] y \approx \frac{F(x+\epsilon y) - F(x)}{\epsilon}$

Jeśli podasz analityczny jakobian, zachowasz wszystkie cyfry (np. 16 w podwójnej precyzji). Metody Kryłowa również opierają się na produktach wewnętrznych, więc jeśli twój Jakub jest źle uwarunkowany melodią , to jest efektywny w liczbie pojedynczej, a Kryłow może stagnować lub zwracać błędne rozwiązania. Może to również zapobiec zbieżności bezpośrednich solverów. Czasami można zastosować metody wielosiatkowe do homogenizacji do grubej siatki z kondycjonowalnym układem. Gdy problemu nie da się sformułować z lepszym uwarunkowaniem, warto pracować z poczwórną precyzją. (Jest to obsługiwane na przykład przez PETSc.) $10^{16}$

Zauważ, że te same problemy dotyczą metod quasi-Newtonowych, jednak bez różnicowania skończonego. Wszystkie skalowalne metody rozwiązywania problemów z operatorami niedokładnymi (np. Równania różniczkowe) muszą wykorzystywać informacje jakobskie do wstępnego kondycjonowania.

Jest prawdopodobne, że fsolvealbo nie używał informacji jakobskiej, albo używał metody dogleg lub zmiany, aby robić postępy metodą „gradientu opadania”, pomimo zasadniczo pojedynczego jakobiańskiego (tj. Skończone różnicowanie miałoby dużo „szumu” z arytmetyka o skończonej precyzji). Nie jest to skalowalne i fsolveprawdopodobnie zwiększa się wraz ze wzrostem rozmiaru problemu.

Globalizacja

Jeśli problemy liniowe zostaną rozwiązane dokładnie, możemy wykluczyć problemy związane z problemem liniowym (Kryłow) i skupić się na nich z powodu nieliniowości. Lokalne minima i nierówności cechują się powolną konwergencją lub powodują stagnację. Poisson-Boltzmann jest gładkim modelem, więc jeśli zaczniesz od wystarczająco dobrego wstępnego odgadnięcia, Newton zbiegnie się kwadratowo. Większość strategii globalizacji wymaga pewnego rodzaju kontynuacji w celu uzyskania wysokiej jakości wstępnego odgadnięcia ostatecznych iteracji. Przykłady obejmują kontynuację siatki (np. Full Multigrid), kontynuację parametrów i pseudotransient kontynuację. Ten ostatni ma ogólne zastosowanie do problemów w stanie ustalonym i oferuje pewną teorię globalnej konwergencji, patrz Coffey, Kelley i Keyes (2003) . Wyszukiwanie pokazuje ten artykuł, który może ci się przydać:Shestakov, Milovich i Noy (2002) Rozwiązanie nieliniowego równania Poissona-Boltzmanna za pomocą pseudo-przejściowej kontynuacji i metody elementów skończonych . Pseudotransjentowa kontynuacja jest ściśle związana z algorytmem Levenberga-Marquardta.

Dalsza lektura

— Jed Brown
źródło