Fakultät für Mathematik und Informatik 23. Januar 2014 TU Bergakademie Freiberg
Dr. M. Helm, Dr. A. Franke-Börner
Numerik linearer und nichtlinearer Parameterschätzprobleme Computerübung zur Inversion eines 1D-Wärmeleitungsproblems
Zum Abschluss unserer Veranstaltung soll das stationäre eindimensionale Wärmeleitungspro- blem
(c(x)u 0 (x)) 0 = 1 (x ∈ (0, 1)), u(0) = u(1) = 0 (1) bearbeitet werden. Aus Messwerten von u an verschiedenen Orten x i ∈ (0, 1) wollen wir dabei die Wärmeleitfähigkeit c(x) > 0 rekonstruieren.
Wir wählen dazu n + 1 = 101 auf [0, 1] äquidistante Gitterpunkte x i = ih (i = 0, . . . , n + 1), h = 1
n + 1 .
Die Lösung u(x) wird gemäß u i = u(x i ) (i = 1, . . . , n) diskretisiert, während wir die Leitfä- higkeiten mittig zwischen den Stützstellen schätzen:
c i = c(x i − h
2 ), i = 1, . . . , n + 1.
Die standardmäßige Finite-Differenzen-Diskretisierung von (1) lautet dann
A(c)u = f (2)
mit
A(c) =
c 1 + c 2 −c 2
−c 2 c 2 + c 3 −c 3
. .. . .. . ..
−c n−1 c n−1 + c n −c n
−c n c n + c n+1
, f = h 2
1 1 .. .
1 1
.
Als Lösung der diskreten Aufgabe (2) bei gegebenem Parameter c ergibt sich
u(c) = A −1 (c) f . (3)
Wir werden die Lösung nur an p = 19 ausgewählten Gitterpunkten messen, nämlich bei x 5 = 0.05, x 10 = 0.1, x 15 = 0.15, . . . x 95 = 0.95.
Als Daten verwenden wir somit einen Vektor der Form
d = QA −1 (c true ) f + r ε .
Inversion der stationären 1D-Wärmeleitungsgleichung 2
Hierbei ist r ε ein Vektor N (0, ε 2 )-verteilter Zufallszahlen (auch ε = 0 zulässig) und Q ∈ R p×n mit
Q(j, :) =
e j ∈ R p , falls x j Beobachtungspunkt;
0, sonst.
der sogenannte Beobachtungsoperator. Für unsere Tests werden wir für c true die Diskretisie- rung von
c(x) = 1 + 3 exp(−100(x − 0.18) 2 ) + 1.5 exp(−50(x − 0.6) 2 ) + 2 exp(−100(x − 0.8) 2 ) verwenden.
Die Inversion erfolgt durch Minimierung des regularisierten Zielfunktionals f (c) = 1
2 kQA −1 (c) f − dk 2 2 + β 2
2 kW (c − c ref )k 2 2 (4) mittels Gauß-Newton-Verfahren. Der Vektor c ref beinhaltet evtl. Vorkenntnisse über c und kann bei uns konstant gewählt werden. Die Regularisierungsmatrix W besitzt die Gestalt
W =
k 1 I k 2 D 1 k 3 D 2
∈ R [(n+1)+(n−1)+(n−1)]×(n+1) .
Die Matrizen D 1 , D 2 ∈ R (n−1)×(n+1) repräsentieren die diskretisierte erste bzw. zweite Ablei- tung von c(x), und haben die Gestalt
D 1 = 1 2h
−1 0 1 . .. ... ...
−1 0 1
bzw. D 2 = 1 h 2
1 −2 1 . .. ... ...
1 −2 1
.
In unseren Experimenten könnte man z. B. zunächst k 1 = k 2 = k 3 wählen.
Ein wichtiger Punkt bei der Verwendung des Gauß-Newton-Verfahrens ist die Kenntnis der ersten Ableitung J (c) = r 0 (c) von
r(c) =
Qu(c) − d βW (c − c ref )
.
mit u(c) = A −1 (c) f notwendig. Es gilt r 0 (c) =
Qu 0 (c) βW
∈ R (p + 3n − 1) × (n + 1),
so dass nur noch die Berechnung der Ableitung u 0 (c) ∈ R n×(n+1) Schwierigkeiten bereitet.
Die Spalten von bestehen u 0 (c) sind gerade die partiellen Ableitungen von u(c). Für diese gilt
∂u
∂c j = ∂A −1 (c)
∂c j f
Durch partielles Differenzieren der Gleichung A(c)A −1 (c) = I ergibt sich
∂A −1 (c)
∂c j = −A −1 (c) ∂A(c)
∂c j A −1 (c),
Inversion der stationären 1D-Wärmeleitungsgleichung 3
wobei es sich bei den Ableitungen A j = ∂A(c) ∂c
j
