Übungsblatt 4
Nichtlineare finite Elemente
Anmerkung: Die Abgabe von Arbeitsblatt 1 bis 4 ist verpflichtend und konstituiert die Studienleistung der Veranstaltung Simulationstechniken. Die Arbeitsblätter führen von der mathematischen Formulierung eines Modellproblems hin zur numerischen Lösung dieses Problems und bauen aufeinander auf. Zum Bestehen der Veranstaltung müssen auf jedem Blatt mindestens 50% der erzielbaren Punkte erreicht werden.
Geben Sie bei allen Aufgaben die Lösungswege und Zwischenergebnisse mit an. Das Endergebnis alleine ist nicht ausreichend! Wir empfehlen Ihnen die Nutzung von Python und Jupyter-Notebooks. Sollten Sie ein Jupyter-Notebook verwenden, dann können Sie dieses einfach direkt als Lösung bei uns einreichen. In allen anderen Fällen erzeugen Sie bitte ein PDF und legen die numerischen Codes als separate Datei dazu.
Sie werden durch die einzelnen Schritte der Modellimplementierung geleitet, und wir geben Hinweise zur Implementierung. Es ist nicht zwingend notwendig, diese 1-zu-1 zu verfolgen. Im Rahmen dieser Hinweise finden Sie Codeabschnitte, die Sie verwenden können. Sie dürfen natürlich auch die Codebeispiele aus dem Vorlesungsmaterial hier verwenden.
4.1 Newton-Raphson-Verfahren
4.1.1 Skalarwertige Funktionen
Implementierung eines Newton-Lösers
6 Punkte
Finden Sie eine Nullstelle für beliebige \(f(x)\) numerisch. Schreiben Sie dazu einen Löser, der \(f(x) = 0\) mittels eines einfachen Newton-Iterationsverfahren lokal löst. Nutzen Sie hierfür folgende Schnittstelle:
2 callback=None):
3 """
4 Newton solver expects scalar function fun and scalar
5 initial value x0.
6
7 Parameters
8 ----------
9 fun : callable(x)
10 Scalar function of a scalar.
11 x0 : float
12 Scalar initial value.
13 jac : callable(f, x)
14 Scalar derivative of f at position x.
15 tol : float, optional
16 Tolerance (with respect to the zero value) for
17 convergence.
18 maxiter : int, optional
19 Maximum number of Newton iterations.
20 callback : callable(x, **kwargs), optional
21 Callback function to handle logging and convergence
22 stats.
23
24 Returns
25 -------
26 x : float
27 Approximate local root.
28 """
Abbruch der Newton-Iteration
Das Newton-Verfahren muss zu einem bestimmten Punkt abgebrochen werden. Die Funktionssignatur sieht hierfür zwei Mechanismen vor: Der Parameter tol gibt eine Toleranz \(\varepsilon \) vor, um den sich der Funktionswert von Null unterscheiden darf. Das heißt die Schleife wird abgebrochen, wenn zum Zeitschritt \(i\) der Funktionswert \(|f(x_i)|<\varepsilon \) erfüllt ist. Um Situationen zu vermeiden, in denen keine Konvergenz erreicht werden kann, sieht die Funktion noch eine maximale Anzahl an Schritten (Parameter maxiter) vor. Die Iteration wird abgebrochen, wenn eines dieser beiden Kriterien erfüllt ist. Da das Newton-Verfahren bei Erfüllung des zweiten Kriterium nicht konvergiert ist, sollte hier ein Fehler gemeldet werden. Dies kann z.B. über eine Ausnahme mit dem raise-Befehl geschehen.
Callback
Die o.g. Schnittstelle enthält eine sogenannte Callback-Funktion. Diese sollte der Schnittstelle
2 """Solver callback for logging.
3
4 Parameters
5 ----------
6 x : float or np.ndarray
7 Current approximate solution
8 """
9 ...
gehorchen und, falls angegeben, einmal je Newton-Iteration mit aktueller Näherungslösung \(x\) aufgerufen werden. Dieser Callback kann dann beispielsweise über alle Zwischenlösungen Buch führen, aktuelle Konvergenzkriterien ausgeben und Ihnen so bei der Fehlersuche helfen.
Implementieren Sie einen solchen Callback. Was genau angezeigt wird, bleibt Ihnen überlassen.
Veranschaulichung an einem einfachen Polynom
Testen Sie Ihren Newton-Löser zunächst an dem einfachen Polynom \(f(x)=x^3 + x^2 - x + 1\). Wählen Sie hierfür verschiedene Startwerte. Für \(x_0 = 0.9\) und einem relativen Konvergenzkriterium von \(\varepsilon =0.01\) könnnte die Ausgabe Ihres Callbacks, die besuchten Näherungslösungen und der Konvergenzverlauf z.B. so aussehen:
X FUN JAC = === === 9.000e-01 1.639e+00 3.230e+00 3.926e-01 8.220e-01 2.475e-01 -2.929e+00 -1.262e+01 1.888e+01 -2.261e+00 -3.182e+00 9.810e+00 -1.936e+00 -5.741e-01 6.375e+00 -1.846e+00 -3.827e-02 5.533e+00 -1.839e+00 -2.168e-04 5.471e+00 -1.839e+00 -7.094e-09 5.470e+00 -1.839e+00 -2.220e-16 5.470e+00
Finden Sie nun die Nullstelle im selben Polynom, indem Sie bei \(x_0 = -1.1\) starten. Zeigen Sie die Schritte des Newton-Solvers für die Startpunkte \(x_0=0.9\) und \(x_0=-1.1\), in dem Sie die die aktuelle Position \(x_i\) gegen die Nummer der Iteration \(i\) graphisch plotten. Was beobachten Sie?
Anmerkung: Ihr Newton-Löser sollte für beliebige Funktionen funktionieren. Sie müssen hierzu die Ableitung der Funktion von Hand implementieren. Um z.B. die Nullstelle der Funktion \(f(x)=x^2-1\) zu finden benötigen Sie eine Implementierung der Ableitung \(\dif f/\dif x=2x\), z.B. in der folgenden Form:
Wir verlangen von Ihnen nicht, dass Ihre Implementierung in irgendeiner Form die Ableitung selbst ausrechnet. Ein häufiger Implementierungsfehler ist allerdings, dass die jac übergebene Funktion nicht der Ableitung der Funktion fun entspricht.
4.1.2 Vektorwertige Funktionen
7 Punkte
Verallgemeinern Sie die Implementierung Ihres Newton-Lösers auf vektorwertige Funktionen \(\v {f}(\v {x})\). Die Funktion fun in der obigen Funktionssignatur muss nun einen Vektor (also einen numpy Array) zurückliefern. Die Funktion jac liefert eine Matrix. Es kann Sinn machen, für die vektorwertige Implementierung eine separate Funktion zu implementieren. Nutzen Sie Ihren Solver um das Minimum der Funktion \(g(x,y)=2\exp (-10(x^2+y^2))+x^2+y^2+x\) zu finden.
Bitte beantworten Sie im Rahmen der Lösung dieser Aufgabe folgende Fragen:
- Wir haben bislang über das Newton-Verfahren für die Lösung von gekoppelten nichtlinearen Gleichungen gesprochen. Hier fragen wir nun nach der Minimierung einer Funktion und fordern damit die Lösung eines Optimierungsproblems. Wie passt dies zusammen?
- Welche spezielle Struktur hat die Jacobi-Matrix für ein Optimierungsproblem? Wie nennt man diese Matrix in diesem Fall auch?
- Sie müssen sich überlegen, wann die Newton-Iteration abgebrochen wird. Was könnte hier ein vernünftiges Kriterium sein und warum?
- Wie sieht die Iterationen Ihres Lösers in einem zweidimensionalen Plot in der \(x\)-\(y\)-Ebene, ausgehend von den Startpunkten \((1/2,1/2)\) und \((5,5)\), aus? Zeigen Sie Hintergrund die Funktionswerte \(g(x,y)\) farbkodiert. Dies geht beispielsweise mit der matplotlib-Funktion pcolormesh oder contour.
4.2 Poisson-Boltzmann-Gleichung
Wir wenden uns nun der Lösung der eindimensionalen Poisson-Boltzmann-Gleichung zu. In ihrer entdimensionalisierten Form lautet diese \begin {equation} \frac {\dif ^2 \Phi }{\dif x^2} = \sinh \Phi . \end {equation} Ziel dieser Aufgabe ist es, einen Löser für diese nichtlineare Gleichung zu implementieren.
4.2.1 Diskretisierung
3 Punkte
Zunächst muss die nichtlineare Poisson-Gleichung diskretisiert werden. Führen Sie diese Diskretisierung eigenständig durch und leiten Sie alle Schritte her. Unterscheiden Sie dabei zwischen der Formulierung der schwachen Form (inklusive der Reduktion der Differenzierbarkeitsanforderungen), des Galerkin-Ansatzes und des finite-Element Ansatzes. Nutzen Sie lineare Elemente und reguläre Gitter. Dies entspricht exakt dem Vorgehen aus den vorherigen Übungsblättern.
Anmerkung: Wenn Sie finite Elemente ansetzen, dann reicht es die Basisfunktionen mit Hilfe der Formfunktionen auszudrücken. Die exakten Ausdrücke für die Formfunktionen spielen hier noch keine Rolle.
4.2.2 Numerische Quadratur des Residuums
2 Punkte
Die finalen Gleichungen für das (diskrete) Residuum enthalten Terme der Form \((\varphi _k, \sinh \Phi _N)\), wobei \(\varphi _k\) die Basisfunktionen und \(\Phi _N\) die approximierte Lösung der DGL ist. Diese Terme sollen mit Hilfe der Gauß-Quadratur gelöst werden. Schreiben Sie die Skalarprodukte/Integrale die nicht analytisch gelöst werden können mit Hilfe der Art numerischer Quadratur. Fixieren Sie hier noch nicht die Anzahl der Quadraturpunkte, sondern leiten Sie diese Gleichungen für eine beliebige Anzahl Quadraturpunkte her.
4.2.3 Tangentenmatrix
3 Punkte
Leiten Sie die Tangentenmatrix her. Zeigen Sie, welche Form die Tangentenmatrix in der linearisierten Form annimmt. Vergleichen Sie diese linearisierte Tangentenmatrix mit der Systemmatrix des linearen Problems.
4.2.4 Implementation des Residuumvektors und der Tangentenmatrix
4 Punkte
Anmerkung: Die linearisierte Gleichung haben Sie im Rahmen von Übungsblatt 2 numerisch gelöst. Evtl. können Sie auf der Implementierung dieser Lösung hier aufbauen.
Implementieren Sie den Residuumsvektor und die Tangentenmatrix. Wir schlagen folgende Signaturen für die Funktionen vor, welche die Berechnung implementieren:
2 potential_left=0, potential_right=0,
3 dx=1, nb_quad=2, linear=False):
4 """
5 Assemble global residual vector for a specific potential.
6
7 Parameters
8 ----------
9 potential_g : np.ndarray
10 Current potential on the nodes (the expansion
11 coefficients); the length of the array is the number
12 of nodes.
13 potential_left : float
14 Left Dirichlet boundary condition.
15 potential_right : float
16 Right Dirichlet boundary condition.
17 dx : float, optional
18 Grid spacing. (Default: 1)
19 nb_quad : int, optional
20 Number of quadrature points. (Default: 2)
21 linear : bool, optional
22 Linearize mass matrix. (Default: False)
23
24 Returns
25 -------
26 residual_g : np.ndarray
27 Residual vector (same shape as ‘potential_g‘)
28 """
29 ...
30
31def tangent(potential_g,
32 potential_left=0, potential_right=0,
33 dx=1, nb_quad=2, linear=False):
34 """
35 Assemble global tangent matrix for a specific potential.
36
37 Parameters
38 ----------
39 potential_g : np.ndarray
40 Current potential on the nodes (the expansion
41 coefficients); the length of the array is the number
42 of nodes
43 potential_left : float
44 Left Dirichlet boundary condition.
45 potential_right : float
46 Right Dirichlet boundary condition.
47 dx : float, optional
48 Grid spacing. (Default: 1)
49 nb_quad : int, optional
50 Number of quadrature points. (Default: 2)
51 linear : bool, optional
52 Linearize mass matrix. (Default: False)
53
54 Returns
55 -------
56 tangent_gg : np.ndarray
57 Tangent matrix (quadratic, number of rows and columns
58 equal number of nodes)
59 """
60 ...
Bei der Implementierung sollten Sie schrittweise vorgehen: Ignorieren Sie zunächst die Dirichlet-Bedingungen und implementieren Sie lediglich den Laplace-Operator. Schreiben Sie den Code für die Massematrix (der Term \(\sinh \Phi \) der Differentialgleichung) erst dann, wenn der Laplace-Operator funktioniert. Die Option linear ist nützlich, um mit der linearisierten Lösung vergleichen zu können.
Anmerkung: Es ist wichtig, dass die Funktion tangent die korrekte(n) Ableitung(en) von residual liefert. Dies kann numerisch getestet werden. So können Sie z.B. numerisch über den Differenzenquotienten die Ableitung aus der Funktion residual berechnen und mit der analytischen Berechnung von tangent überprüfen. Im folgenden finden Sie einen Codeblock, der diese numerische Berechnung der Ableitung für Sie macht:
2 eps=1e-6):
3 """
4 Check that tangent_fun is gives the derivative_fun
5 using finite differences.
6
7 Parameters
8 ----------
9 value_g : numpy.nd_array
10 Nodal value for which to check the derivative
11 residual_fun : callable
12 Function that takes the values and returns an array
13 of residual values
14 tangent_fun : callable
15 Function that takes the values and return the
16 tangent/jacobian matrix
17 eps : float, optional
18 Finite difference used for numeric computation of
19 the derivative (Default: 1e-6)
20 """
21 nb_nodes = len(value_g)
22 tangent_gg = tangent_fun(value_g)
23 numeric_tangent_gg = np.zeros_like(tangent_gg)
24 for i in range(nb_nodes):
25 _value_g = value_g.copy()
26 _value_g[i] += eps
27 residual_plus_g = residual_fun(_value_g)
28 _value_g[i] -= 2 * eps
29 residual_minus_g = residual_fun(_value_g)
30 numeric_tangent_gg[:, i] = (
31 residual_plus_g - residual_minus_g
32 ) / (2 * eps)
33 np.testing.assert_array_almost_equal(
34 tangent_gg, numeric_tangent_gg)
35
36
37# Check if tangent is implemented correctly
38for i in range(10):
39 check_tangent(np.random.random(21)-0.5, residual, tangent)
4.2.5 Anwendung des Newton-Lösers
7 Punkte
Anmerkung: Falls ihr Newton-Löser für multidimensionale Probleme nicht funktioniert, dürfen Sie hier auch auf den Newton-Löser im Paket scipy zurückgreifen. Diesen finden Sie unter scipy.optimize.newton.
Implementierung Sie die Lösung der nichtlinearen Poisson-Boltzmann-Gleichung mit zwei Dirichlet Randbedingungen, jeweils am linken und rechten Rand. Zeigen Sie die Lösungen der Gleichung für \(1\), \(2\) und \(3\) Gaußsche Quadraturpunkte für ein System der Länge \(L=5\lambda \) und einem Potential von \(-1 k_B T/|e|\) auf der linken Elektrode und \(5 k_B T/|e|\) auf der rechten Elektrode, wobei \(\lambda \) die Debye-Länge ist. Zeigen Sie eine Lösung mit \(\approx 10\) und \(\approx 100\) Knoten. Vergleichen Sie diese Lösung der Lösung der linearisierten Gleichung.
Anmerkung:
- Sie müssen die Newton-Iteration von einem bestimmten Potential \(\Phi \) starten. Was ist hier ein guter Startpunkt?
-
Wenn Sie obige Funktionssignaturen verwandt haben, sollte der Aufruf des Newton-Lösers ungefähr so aussehen:
1nb_nodes = 21 # Number of notes
2potential0_g = ... # Initial condition
3potential_g = newton(residual, potential0_g, tangent)Unter Umständen wollen Sie noch eine Callback-Funktion übergeben, um die Konvergenz des Lösers zu verfolgen.
Zusätzlich sollten Sie noch die Potentialrandbedingungen und zusätzliche Parameter, wie z.B. die Anzahl der Quadraturpunkte, an die Funktionen residual und tangent übergeben. Dies können Sie z.B. mit lambda-Ausdrücken realisieren.
Punkte
Sie können auf diesem Übungsblatt insgesamt 32 Punkte erzielen.