Wann ist Newton-Krylov kein geeigneter Löser?

Kürzlich habe ich verschiedene nichtlineare Löser von scipy verglichen und war besonders beeindruckt vom Newton-Krylov-Beispiel im Scipy-Kochbuch, in dem sie eine Differentialgleichungsgleichung zweiter Ordnung mit nichtlinearem Reaktionsterm in etwa 20 Codezeilen lösen.

Ich habe den Beispielcode geändert, um die nichtlineare Poisson-Gleichung ( auch Poisson-Boltzmann-Gleichung genannt , siehe Seite 17 in diesen Hinweisen) für Halbleiter-Heterostrukturen zu lösen , die die folgende Form haben:

\frac{d^{2} ϕ}{d 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$

(Dies ist die Restfunktion, die an den Solver übergeben wird.)

Dies ist ein elektrostatisches Problem, bei dem und nichtlineare Funktionen für die Form . Die Details hier sind nicht wichtig, aber der Punkt ist, dass die nichtlineare Funktion exponentiell mit variiert, so dass die Restfunktion über einen großen Bereich variieren kann ( $n(x,\phi)$ $p(x,\phi)$ $n_i(x) e^{-(E_i(x,\phi) - E_f)}$ $\phi$ $10^{-6} - 10^{16})$ mit einer leichten änderung in . $\phi$

Ich löse diese Gleichung numerisch mit Newton-Krylov von scipy, aber sie würde niemals konvergieren (tatsächlich würde sie immer einen Fehler bei der Berechnung des Jacobian melden). Ich habe von einem Newton-Krylov- Solver auf fsolve (basierend auf MINPACK hybrd) gewechselt und es hat zum ersten Mal funktioniert!

Gibt es allgemeine Gründe, warum Newton-Krylov für bestimmte Probleme nicht geeignet ist? Müssen die Eingangsgleichungen irgendwie konditioniert werden?

Möglicherweise sind weitere Informationen erforderlich, um einen Kommentar abzugeben, aber warum hat fsolve Ihrer Meinung nach in diesem Fall funktioniert?

python nonlinear-equations poisson scipy Boyfarrell
quelle

Ich hatte das gleiche Problem mit Newton-Krylov, der mit dem Jacobian versagte, und stellte fest, dass das Ändern der Methode von "lgmres" zu "gmres" ( sol = newton_krylov(func, guess, method='gmres')) das Problem behebt . Ich weiß nicht genau warum, aber jeder, der dieses Problem hat, könnte das Gleiche in Betracht ziehen.

Arthur Dent

Antworten:

Es gibt zwei Probleme, auf die Sie wahrscheinlich stoßen.

Schlechte Konditionierung

Erstens ist das Problem schlecht konditioniert, aber wenn Sie nur ein Residuum bereitstellen, wirft Newton-Krylov die Hälfte Ihrer signifikanten Ziffern weg, indem er das Residuum endlich differenziert, um die Aktion des Jakobianers zu erhalten:

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

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

$10^{16}$

Beachten Sie, dass die gleichen Probleme für Quasi-Newton-Methoden gelten, jedoch ohne endliche Differenzierung. Alle skalierbaren Methoden für Probleme mit nicht kompakten Operatoren (z. B. Differentialgleichungen) müssen Jacobi-Informationen zur Vorkonditionierung verwenden.

Es ist wahrscheinlich, dass fsolveentweder keine Jacobi-Informationen verwendet wurden oder dass eine Dogleg-Methode oder eine Verschiebung verwendet wurde, um Fortschritte mit einer "Gradientenabstieg" -Methode zu erzielen, obwohl eine im Wesentlichen singuläre Jacobi-Methode (dh eine endliche Differenzierung würde viel "Rauschen" verursachen) endliche Präzisionsarithmetik). Dies ist nicht skalierbar und wird fsolvewahrscheinlich langsamer, wenn Sie das Problem vergrößern.

Globalisierung

Wenn die linearen Probleme genau gelöst sind, können wir Probleme im Zusammenhang mit dem linearen Problem (Krylov) ausschließen und uns auf solche konzentrieren, die auf Nichtlinearität zurückzuführen sind. Lokale Minima und nicht glatte Merkmale verlangsamen die Konvergenz oder verursachen eine Stagnation. Poisson-Boltzmann ist ein glattes Modell. Wenn Sie also mit einer ausreichenden anfänglichen Schätzung beginnen, konvergiert Newton quadratisch. Die meisten Globalisierungsstrategien beinhalten eine Art Fortsetzung, um eine qualitativ hochwertige erste Schätzung für die endgültigen Iterationen zu erhalten. Beispiele hierfür sind die Netzfortführung (z. B. Full Multigrid), die Parameterfortführung und die pseudotransiente Fortführung. Letzteres ist allgemein auf stationäre Probleme anwendbar und bietet einige globale Konvergenztheorien (siehe Coffey, Kelley und Keyes (2003)) . Bei einer Suche wird dieses Dokument angezeigt, das für Sie nützlich sein kann:Shestakov, Milovich und Noy (2002) Lösung der nichtlinearen Poisson-Boltzmann-Gleichung unter Verwendung der pseudotransienten Fortsetzung und der Finite-Elemente-Methode . Die pseudotransiente Fortsetzung ist eng mit dem Levenberg-Marquardt-Algorithmus verwandt.

Weitere Lektüre

Jed Brown
quelle