wissenschaftliches rechnen an ausgewahlten beispielen¨ · er f¨urin der realit¨at gestellte...
TRANSCRIPT
Wissenschaftliches Rechnenan ausgewahlten Beispielen
Prof. Dr. Gabriel Wittum
Dr. Nicolas Neuß
Universitat HeidelbergInterdisziplinares Zentrum f¨ur Wissenschaftliches Rechnen
Sommersemester 2000
In LATEX gesetzt von Matthias Kleiser
Wissenschaftliches Rechnen Seite 2
Inhaltsverzeichnis
1 Modellierung und Simulation 3
1.1 Einfuhrung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.2 Differentialgleichungen 2. Ordnung . . . . . . . . . . . . . . . . 6
1.3 Approximation der Dgl.:”Diskretisieren“ . . . . . . . . . . . . . 9
1.3.1 Zeitdiskretisierung . . . . . . . . . . . . . . . . . . . . . 9
1.3.2 Ortsdiskretisierung - Differenzenverfahren . . . . . . . . 11
1.3.3 Finite Elemente . . . . . . . . . . . . . . . . . . . . . . . 12
1.3.4 Finite-Volumen-Verfahren . . . . . . . . . . . . . . . . . 15
2 Das HautproblemDiffusion von Arzneimitteln durch die menschliche Haut 21
2.1 Modellierung . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.2 Diskretisierung . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.3 Wahl des Gitters . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
2.4 Wahl eines geeigneten L¨osungsverfahrens . . . . . . . . . . . . . 25
2.5 Algorithmus des Mehrgitterverfahrens . . . . . . . . . . . . . . . 28
2.6 Charakterisierung der Membran durch Lag-Zeit . . . . . . . . . . 28
3 Numerik der Homogenisierung 32
3.1 Problemstellung . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
3.2 Homogenisierung eindimensionaler Diffusionsgleichungen . . . . 34
3.3 Der mehrdimensionale Fall . . . . . . . . . . . . . . . . . . . . . 37
3.4 Homogenisierung von zeitabh¨angigen Problemen . . . . . . . . . 41
4 Eigenschwingungen von Seen 47
4.1 Modell . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
4.2 Herleitung der Flachwassergleichungen . . . . . . . . . . . . . . 49
4.3 Dreiecksgitter . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
4.4 Das Algebraische Mehrgitterverfahren (AMG) . . . . . . . . . . . 55
4.5 Mehrgitterverfahren f¨ur Eigenwertprobleme . . . . . . . . . . . . 59
Wissenschaftliches Rechnen Seite 3
1 Modellierung und Simulation
1.1 Einfuhrung
Wissenschaftliches Rechnen: scientific computing, computational science (Rech-nergest¨utzte Wissenschaft).
Das Wissenschaftliche Rechnen besch¨aftigt sich grunds¨atzlich mit der Projektionreeller Probleme in die Welt der Mathematik, der sog.Modellierung, mit demdortigen Losen des Problems, d.h. Extrahieren der n¨otigen Situation, und deranschließenden R¨uckubersetzung in die Realit¨at (Simulation).
Als Erster beschritt diesen Weg Pythagoras von Samos etwa um 540 v.Chr., indemer fur in der Realitat gestellte Fragen mathematischeAquivalente suchte und daf¨urallgemeing¨ultige, auch auf verwandte Probleme anwendbare Aussagen fand.
realeObjekte mathematische
ObjekteModellierung
Simulation
Lösen
WissenschaftlichesRechnen
Realität Mathematik
Beispiel:
Rad, Muhlstein - KreisFeld - Rechteck
etc.
Das Ziel des Wissenschaftlichen Rechnens ist es, die Realitat durch Quantifizie-ren und Mathematisieren zu verstehen. Das mathematische Objekt enthalt diebenotigte Information, jedoch meist implizit .
Wissenschaftliches Rechnen Seite 4
Prozesse:
Prozesse sind dynamische, d.h. zeitabhangige Vorgange in der Realitat (Naturund Technik), welche meist durch Systeme von partiellen Differentialgleichungenmodelliert werden konnen.
Frage: Wie bekommt man Modelle ?
Modell:
(lat.: modellum; Abstraktion) zuerst benutzt in der Planung des Kirchenbaus um1100 n.Chr.; Beschrankung auf das Wesentliche; so einfach wie notig, so kompli-ziert wie notig.
In der Physik (Leibniz, Newton):
Die Erhaltungssatze von Masse, Energie, Impuls und Drehimpuls.
Beispiel: Diffusion
Vc
c(~x; t) ist die Konzentration einer Substanz im Volumen V
Annahmen:
� Es befinden sich keine Quellen und Senken in V,
� Die Substanz ist inkompressibel, d.h. keine Anderung der Dichte.
1.Schritt: Formulieren des Erhaltungssatzes
Nach dem Modell der Massenerhaltungkann sich die Konzentration c(~x; t) in Vnur andern, wenn Substanz uber den Rand von V fließt.
ZV
dc
dt=Z@V
~j � ~n = Fluß uber den Rand von V
~n ist der Normalenvektor zum Rand des Volumens V, der Strom- oder Flußdich-tevektor ~j wird auf den Rand projiziert. Einheit: [j] = [c]
s� [LE]
Wissenschaftliches Rechnen Seite 5
2.Schritt: Phanomenologischer Zusammenhang zwischen ~j und ~c, umdas Modell auf eine Differentialgleichung zu reduzieren�!
”konstitutiveBeziehung“
Fick’sches Gesetzfur die Diffusion bzw. Warmeleitung:
~j = D � rxc(~x; t)
d.h. die Stromdichte hangt ab vom Gradienten der Konzentration (Konzentrati-onsgefalle) und der Durchlassigkeit D (vgl. Maxwell-Gleichung) mit D > 0.
3.Schritt:”Einsetzen“
(1:1)ZV
dc
dt=Z@V
~j � ~n =Z@V
D � rxc(~x; t) � ~n
Der Gaußsche Integralsatzliefert :
ZV
dc
dt=ZV
divD � rc(~x; t)
und schließlich die differentielle Form :
(1:2)dc
dt= divD � rc
Wissenschaftliches Rechnen Seite 6
1.2 Differentialgleichungen 2. Ordnung
(1:3) k � u(~x) = f(~x) ; ~x 2 � Rd mit
k =dX
i;j=1
@
@xiaij
@
@xj+
dXi=1
bi@
@xi+ q
Definition 1.1.: (Typeinteilung)
Sei A = (aij)di;j=1 und ~b =
0BB@
b1...bd
1CCA
Dann heißt die Differentialgleichung (1.3) :
� elliptisch, wenn A f pos:
neg:g definit,
� parabolisch, wenn A f pos:
neg:g semidefinit und (Ajb) 2 Rd�d+1 vollen
Rang d hat,
� hyperbolisch, wenn d-1 Eigenwerte von A dasselbe Vorzeichen habenund einer das gegenseitige.
Notabene :
� A ist oBdA symmetrisch nach dem Satz von Schwarz.
� Die Begriffsbildung stammt aus der Aquivalenz zu den entsprechenden geo-metrischen Gleichungen.
Beispiel: A =
1 23 4
!; b = 0; q = 0
k =@2
@x21+2
@
@x2
@
@x1+3
@
@x1
@
@x2+4
@2
@x22
S:v:Schwarz=
@2
@x21+5
@
@x2
@
@x1+4
@2
@x22=
=@2
@x21+
5
2
@
@x2
@
@x1+
5
2
@
@x1
@
@x2+ 4
@2
@x22
Wissenschaftliches Rechnen Seite 7
) A ist symmetrisierbar in A0 =
1 5
252
4
!; b = 0; q = 0
Eigenschaften von Losungen der Dgln. aus Definition 1.1.:
elliptisch: Bei einem elliptischen Problem handelt es sich um ein ortsabhangigesRandwertproblem (RWP), dessen Losung global abhangig ist(siehe (1.5)).
Beispiel: Potentialgleichung (1.4) �ru = 0 in
K1(0) = f(x; y) 2 R2 : k(x; y)k < 1g = + Randbedingung uj@ = u0
Polarkoordinaten : x = r sin'; y = r cos'
! Sk := rk sin k'; ck = rk cos k'
genugen der Potentialgleichung (1.4) und sind beliebig oft differenzierbar (Fourier-Basis).Liegen die Randwerte u0 als Fourier-Reihe vor :
u0(1; ') =1Xj=1
(aj cos j' + bj sin j') + a0
dann ist:
(1:5) u(r; ') =1Xj=1
rj (aj cos j' + bj sin j') + a0:
parabolisch: zeitabhangiges Problem
K = ( ddt� Kx) ; Differentialoperator�Kx ist elliptisch im Ort
du
dt= Kx u + f ! Anfangsrandwertproblem(ARWP)
Beispiel:Diffusion, Warmeleitung, . . .
Beispiel (1.2) ist parabolisch, wenn die Matrix D > 0 und konstant ist.
Sei nun D = 1 und die Dimension = 1, d.h. d = 2 (Ort und Zeit),Ortsintervall = (0; �); t > 0.
Wissenschaftliches Rechnen Seite 8
Anfangs- und Randbedingung:
c(x; t = 0) = c0(x); c(0; t); c(�; t)
Anfangspunkt als Fourierreihe gegeben:
c0(x) =1Xk=0
ak cos kx
) ck(x; t) = e�k2t cos kx ist Losung von
dc
dt= D
d2c
dx2; D = 1:
denn:dck(x; t)
dt= �k2 ck(x; t) =
d2c
dx2
) allgemeine Losung:
c(x; t) =1Xk=0
ak e�k2t cos kx
Fur t ! 1 gilt:
c(x; t)t!1�! a0 = const:
”Dampfung“
Die allgemeine Zeitdampfung1Pi=0
ai e��it ui(x) (f = 0) ist abhangig von
den Eigenwerten �i von Kx; je großer die �i , desto schneller wird gedampft.Dies ist das generelle Verhalten parabolischer Differentialgleichungen.
Wissenschaftliches Rechnen Seite 9
1.3 Approximation der Dgl.: ”Diskretisieren“
� Ersetze das Gebiet G durch ein Gitter Gh
� Ersetze Differentialgleichungen und Funktionen durch endlichdimensionaleBeziehung
Gh
Gh
GG
1.3.1 Zeitdiskretisierung
Die Geometrie der Zeitdiskretisierung ist eindimensional:
ti = ti�1 + �i�1; Zeitschrittweite �i > 0Anfangszeitpunkt t0 = 0
Die einfachste Variante ist das explizite Euler-Verfahren (Polygonzug-Methode):Ersetze
dc
dt=
c(t+ � ) � c(t)
�+ O(� ) , d.h. Fehler geht mit � ! 0 zuruck
durch die explizite Form
c(x; t+ � ) � c(x; t)
�= kx c(x; t)
) c(x; t+ � ) = c(x; t) + � kx c(x; t) + f = (1 + � kx) c(x; t) + f
Beispiel:du(t)
dt= ��u(t) ; � > 0 (gewohnliche Dgl.)
Wissenschaftliches Rechnen Seite 10
hat die Losung:u(t) = e��t
Es sei t0 = 0, y(t0) = y0 sowie � > 0 konstant.Mit explizitem Euler gilt:
u(t+ � )� u(t)
�= ��u(t)
) u(t+ � ) = u(t)� �� u(t) = (1� �� )u(t)
) u(k� ) = (1� �� )k u0; �� > 0
Es mußalso folgende Stabilitatsbedingung(”CFL-Bedingung“ ) erfullt sein:
(1:6) 1� � � > 0 () � <1
�
Im Fall einer parabolischen Dgl. in x und t hangt � von h (Schrittweite im Ort)ab, wobei typischerweise gilt:
� = �(h) = O(h2)
Die Dgl. dudt
= Kx u ist parabolisch mit unendlicher Ausbreitungsgeschwindig-keit und elliptischem DifferentialoperatorKx ! Kx;h (n�n - Matrix mit n � 1
h2).
Damit ist � < 1O(h�2) = O(h2) und (1.6) beschrankt � nach oben.
Konsequenz: Explizite Verfahren sind bei parabolischen Gleichungen fur großeZeiten ungeeignet!
Alternative: Implizites Euler-Verfahren
u(t+ � )� u(t)
�= k u(t+ � )
Beispiel:du(t)
dt= ��u(t) ; � > 0
Mit implizitem Euler:
u(t+ � )� u(t)
�= ��u(t+ � ) ; � > 0
) u(t+ � ) = ��� u(t+ � ) + u(t) =1
1 + ��u(t)
Wissenschaftliches Rechnen Seite 11
) u(k� ) =1
(1 + �� )ku(0); �� > 0
Strategie zur Zeitdiskretisierung parabolischer Dgln.
� Verwende ein implizites Verfahren.
� Wahle den Zeitschritt � entsprechend dem Abklingverhalten der Losungzu Beginn: kleine Schritte ! hohe Genauigkeit,spater: immer großere Schritte! hohe Effizienz,etwa: Startzeitschritt �0 > 0
�i+1 = � �i ; � > 1:
1.3.2 Ortsdiskretisierung - Differenzenverfahren
parabolischer Ansatz: dudt
= div D � ru in � Rd (Ortskoordinate).
Nach der Zeitdiskretisierung uber ein implizites Verfahren zerlegen wir nun dasGebiet in ein kartesisches Gitter h mit Schrittweite h und diskretisieren dieDgl. auf den Teilintervallen (eindimensional):
�
�-h
h
Vorwarts:
(1:7)du
dx=
u(x+ h)� u(x)
h+O(h)
Taylor-Entwicklung:
u(x+ h) = u(x) + hu0(x) +h2
2u00(�)
) u(x+ h)� u(x)
h=
u(x) + hu0(x) + h2
2u00(�)� u(x)
h= u0(x) +
h
2u00(�)
Wissenschaftliches Rechnen Seite 12
Ruckwarts:
(1:8)du
dx=
u(x� h)� u(x)
h+O(h)
Daraus ergeben sich durch Einsetzen der Gleichungen (1.7) und (1.8) die sog.zentralen Differenzen:
u00(x)V orw�arts� d
dx
u(x+ h)� u(x)
h
!
R�uckw�arts� u(x+ h) + u(x� h)� 2u(x)
h2= u00(x) +O(h2)
Vor- und Nachteile des Differenzenverfahrens:
� Das Verfahren ist einfach handhabbar Komplexe Geometrien lassen sich nur schwer umsetzen Das Verfahren ist nicht konservativ (flußerhaltend)
1.3.3 Finite Elemente
Ziel: Erzeugung eines linearen Gleichungssystems mit moglichst dunn besetzterMatrix A.
Die sogenannte”starke Formulierung“ ist zunachst punktweise gegeben mit
u 2 C2() \ C(�) und D > 0:
�divD � ru(x) = f(x) in ; u @ = 0
Wir definieren nun ein Skalarprodukt fur Funktionen g; h aus dem Raum der qua-dratintegrablen Funktionen H 0() := L2() = ff :
Rf2 <1g:
(g; h)0 =Z
g � h
und ersetzen die punktweise Gleichheit durch die Norm kgk0 = (g; g)120 .
Damit folgt: Z
j � divDru� f j2 = 0
Im Raum H01 () der einmal schwach differenzierbaren Funktionen mit kompak-
tem Trager in gilt allgemeiner 8v 2 H 01 ():Z
(�divDru� f) � v = 0
Wissenschaftliches Rechnen Seite 13
Nach C.S.U. bzw. der Normvorraussetzung kxyk = 0 , x = 0 oder y = 0erhalten wir:
�Z
divDru � v =Z
f � v partielleInt:()Z
Dru � rv +Z@
Dru � v =Z
f � v
Nach Voraussetzung ist v = 0 auf @. Damit folgt die Variations-oder schwacheFormulierung:
(1:9)Z
Dru � rv =Z
f � v; 8v 2 H01 ()
Diskretisierung: (Ritz-Verfahren)
Trianguliere das Gebiet in ein Gitter h mit den Gitterknotenpunkten x1; : : : xnund wahle als Basis die sog. Ansatzfunktionen(bzw. Formfunktionen,shape functions):
bi(xj) =�1 falls i = j;
0 sonst
x i
h
b > 0i
Ersetze außerdem H01 () =: V durch Vn := spanfb1; : : : ; bng � V mit
dimVn = n <1 und1Sk=1
Vn = V .
Zur Losung von (1.9) suchen wir nun eine Funktion u 2 Vn mitZ
Dru � rv =Z
f � v; 8v 2 Vn
mit dem Ansatz:
(1:10) u =nXi=1
�ibi
Wissenschaftliches Rechnen Seite 14
Durch Einsetzen erhalten wir:Z
DnXi=1
�irbi =Z
fbj ; j = 1; : : : ; n
()nXi=1
�i
Z
Drbirbj| {z }
=:aij
=Z
fbj
| {z }=:fj
; j = 1; : : : ; n
Wir definieren:
A = (aij)ni;j=1 ; x = (�1; : : : ; �n)
T ; g = (f1; : : : ; fn)T
Dies liefert uns schließlich das gewunschte LGS:
Ax = g
mit dessen Losung x nach Formel (1.10) die gesuchte Funktion u berechnet wer-den kann.
Eigenschaften von A :
� aij = aji , d.h. A = AT (A symmetrisch),
� Ist D(x; y) � � > 0 , dann ist �divDru elliptisch,
� aij =RDrbirbj ) xTAx =
RD|{z}>0
ruru| {z }>0
> 0
) A ist positiv definit,
� A ist dunn besetzt, d.h. ]faij 6= 0g = O(n), da die lokalen Tragerbj(xi) = Æij linear auf jedem Dreieck, also bj(x; y) = q1x+ q2y + q3) Variable in Knotenpunkten lokalisiert.
Vor- und Nachteile der Finiten Elemente:
� sehr flexibel� strenger mathematischer Zugang nicht konservativ
Wissenschaftliches Rechnen Seite 15
1.3.4 Finite-Volumen-Verfahren
Die Finite-Volumen-Methode kommt aus der Diskretisierung von Erhaltungssatzender Form
d
dt
ZG
�(w)dx+Z@G
~H(w) � ~n ds = 0
Dabei ist: G � � Rn,w : ! RN ; � : �RN ! RM ,~H : �RM �Rn (Flußvektor)
Indem wir nun G infinitesimal klein machen und den Gaußschen Integralsatz an-wenden, erhalten wir die punktweise Formulierung:
@
@t�(x;w(x; t)) + div ~H(x;w(x; t)) = 0
Beispiel:
� Skalare Erhaltungsgleichung:
w(x; t) : Volumenanteil einer Spezies
�(x;w) = %(x)w(x) : Masse dieser Spezies, wobei % die Dichtedes Gases/der Flussigkeit ist
H(x;w) = %(x)~v w(x) : Flußbei vorgegebenem Stromungsfeld ~S
oder:H(x;w) = %(x) k(x)rw(x) : diffusiver Fluß
� System von Erhaltungsgleichungen (Euler-Gleichungen, kompressibel,ohne Viskositat):
w =
0BBB@
%
u1u2#
1CCCA Dichte x-Komponente der Geschwindigkeit y-Komponente der Geschwindigkeit Temperatur
�(w) =
0BBB@
%
% u1% u2% �
1CCCA Massendichte Impuls in x-Richtung Impuls in y-Richtung innere Energie � = e(#) + 1
2 juj2 + P (%;#)%
H(w) =
0BBB@
% u1 % u2% u21 + P (%; #) % u1 u2
% u1 u2 % u22 + P (%; #)% � u1 % � u2
1CCCA mit @
@t%+div(% u) = 0.
Wissenschaftliches Rechnen Seite 16
Start der Finite-Volumen-Methode:
(1:11)d
dt
ZG
�(w)dx+Z@G
~H(w) � ~n ds = 0
� Unterteile das Gebiet in Kontrollvolumina Bi und fordere (1.11) 8Bi,
� Wahle w aus einem Ansatzraum (z.B. stuckweise polynomial auf dem Kon-trollvolumen) sowie passende Freiheitsgrade,
� Approximation vonR�ij
~H(w) � ~nds , mit �ij = Interface zw. Bi u.Bj
Beispiel: ut + ux = 0 mit Anfangs- und Randbedingungen
� -hxi xi+1
ui�1 ui+1 ui
Zellenzentrierte Finite Volumen:
ui =1
h
xi+1Zxi
u(x)dx =1
h
ZBi
u(x)dx ; xi+1 � xi = h
) d
dt
ZBi
u(x)| {z }=h�ui
dx =ZBi
uxdx =Z@Vi
u(x)ds = u(xi+1)� u(xi)
Mit H(u) x=xi =ui+ui�1
2ergibt sich die Gleichung fur ui:
@
@tui =
ui+1 � ui
2h� ui�1 + ui
2h=
ui+1 � ui�1
2h
D.h. der Konvektionsterm wird mit zentralen Differenzen diskretisiert�!KonsistenzordnungO(h2) , aber es konnen Oszillationen bei nicht glat-ten Losungen auftreten.
Alternative:”Upwind-Diskretisierung“ H(w) x=xi = wi�1
) @
@tui =
ui � ui�1
h(linksseitiger Differenzenquotient)
Diese Methode ist stabil und oszillationsfrei, hat aber eine großere numeri-sche Dissipation:
ui � ui�1
h=ui+1 � ui�1
2h� h
2� ui+1 + 2ui � ui�1
h2
�! Konsistenzordnung O(h).
Wissenschaftliches Rechnen Seite 17
Knotenzentrierte Finite Volumen - Box-Methoden
Ausgangspunkt: �div Dru = f
bzw. integriert:R�divDru =
Rf
Ziel: diskreter Erhaltungssatz, d.h. in jedem Volumen flußerhaltend
� Trianguliere �! Gitter h
Bemerkung: Eine Triangulierung heißt zulassig, falls:
–nSi=1
Ti = (Polygongebiet),
– je zwei Dreiecke sich entweder in einer gemeinsamen Kante oder ingemeinsamen Eckpunkten schneiden,
– der kleinste Winkel �min die Bedingung �min � �0 > 0 erfullt.
� Erzeuge zu h ein”duales“ Gitter (
”Vornoi-Diagramm“):
Verbinde fur jedes Element in h den Kantenmittelpunkt mit dem Element-schwerpunkt! Integrationsboxen Bi
�!m[i=1
Ti = ;m[i=1
Bi = (m = ] der Boxen)
B
Ti
x i
i
� Lokalisieren auf Bi:
�Z
divDru = �mXi=1
ZBi
div Dru =mXi=1
ZBi
f
) �ZBi
div Dru =ZBi
f; 8i = 1; : : : ;m (lokaler Erhaltungssatz)
Wissenschaftliches Rechnen Seite 18
Sei ~n Normalenvektor, dann gilt: 8i = 1; : : : ;m :
(1:12) �ZBi
divDru � 1 =ZBi
f
partielleInt:=) �
Z@Bi
div Dru � ~n| {z }= @u@n
+ZBi
div Dru � 0| {z }
=0
=ZBi
f
� Diskretisierung: Wahle einen Ansatzraum fur u, etwa lineare Finite Ele-mente auf Dreiecken.
x i+1xi-10
1
x
b
i
i
(1:13) u :=NXi=1
�i bi (N = ] der Ansatzfunktionen, bi = Basis von Vn)
Durch Einsetzen der Ansatzfunktionen fur u ergibt sich:
�NXi=1
�i
Z@Bi
D@bi
@~n=ZBi
f
Fur D = 1 liefert der Gaußsche IntegralsatzRBi
divru =R
@Bi
ru � ~n sowie der
ErhaltungssatzRBi
dudt
=R@Bi
du
d~n|{z}ru�~n
eine alternative Formulierung von (1.12):
ZBi
du
dt=Z@Bi
ru � ~n (vgl. auch (1.1) in Kap.1.1)
und mit den Ansatzfunktionen (1.13) fur u:
(1:14)NXj=1
�j
ZBi
dbj(x; y; t)
dt=
NXj=1
�j
Z@Bi
rbj � ~n; 8i = 1; : : : ;m
Wissenschaftliches Rechnen Seite 19
D.h.: m Gleichungen fur N Unbekannte ) m = N .(1.14) ist ein diskreter Erhaltungssatzfur u 2 Vn) Das Finite-Volumen-Verfahren ist flußerhaltend.
Ziel der Box-Methode ist es, fur elliptische Gleichungen die Konvergenzresultatefur Finite Elemente auf die Finite-Volumen-Diskretisierung zu uber- tragen.
Modellproblem:
�div(A(x)ru(x)) = f(x) auf ; u(x) = 0 auf @
Der AnsatzraumSh bestehe aus auf @ verschwindenden Funktionen, die stuckweiselinear auf einer Triangulierung Th sind ( Polygongebiet). Die Boxen seiengemaß dem Vornoi-Diagramm konstruiert. 'j seien Basisfunktionen zum Kno-ten j. Die Diskretisierung fuhrt 8 Boxen Bi auf:Z
Bi
div (Aruh)| {z }=Pj
uj'jh
=ZBi
f(x)A(x)
| {z }=fi
Damit ist mit dem Normalenvektor auf dem Rand ~n :Xj
AFVij uj = fi , AFV
ij =Z@Bi
~nA(x)r'j
Wir vergleichen dies mit dem Finite-Elemente Verfahren:
AFEij =
Ztr('i)\tr('j)
r'iA(x)r'j
und erhalten folgendes Ergebnis:
Satz 1.1.: (Vergleich der Matrizen des FE- und des FV-Verfahrens)
Tn bestehe aus Dreieckselementen, A sei konstant auf jedem Element. Danngilt:
AFE = AFV
Wissenschaftliches Rechnen Seite 20
Beweis: Es genugt, die Behauptung elementweise zu zeigen.
e sei ein Dreieckselement, J das Interface zwischen Bi und Bj .
i
e
h
x
i
J
B
AFV = AFE ,Z
@Bi\e
~n �Ar'j| {z }=:~c
=Ze
r'iAr'j| {z }=:~c
,Z
@Bi\e
~n � ~c = r'iZe
~c ,ZJ
~nJ = area(e) � r'i
, jJ j|{z}= 1
2c
~n =1
2c h � 1
h~n
pq:e:d:
Bemerkung: Es ist AFE = AFV , aber uh;FE 6= uh;FV , denn:
Vergleich der rechten Seiten liefert:
FE: fh =Z
f � bi FV: fh =ZBi
f � 1
=) fh;FV = fh;FE +O(h2):
Wissenschaftliches Rechnen Seite 21
2 Das HautproblemDiffusion von Arzneimitteln durch die menschli-che Haut
Vorteile der transdermalen Applikation:
� Erhohung der Akzeptanz durch den Patienten (Compliance),
� Umgehung des first-pass Metabolismus,
� Vermeidung vorzeitiger Arzneistoffinaktivierung im Gastrointetinaltrakt,! Reduzierung der Arzneistoffdosis,! Verminderung der Nebenwirkungen,
� Stabiler Blutspiegel uber langes Zeitintervall durch gleichmaßige Arznei-stoffinvasion.
Nachteil:Es ist nur eine geringe Arzneistoffaufnahme durch die Haut moglich, da die Horn-haut (stratum corneum) eine etwa 10 - 20 �m dicke Barriere, bestehend aus ca.10 - 15 toten abgeflachten Zellschichten (Korneozyten), darstellt.
2.1 Modellierung
?
6
6?
� -
?
1�m
30�mKorneozyt
LipidkanalDicke 0; 1�m
stratum corneum: etwa 10 Korneozytschichten
Wissenschaftliches Rechnen Seite 22
Fur u = u(x; y; t) ; (x; y) 2 und den Normalenvektor ~n gilt die Diffusionsglei-chung:
�div(Dru) + @c
@t= 0 ,
Z
du
dt=Z
Dru � ~n
mit der DiffusionsmatrixD = D(x; y) =�Dlip , falls (x; y) 2 LipidDcor , falls (x; y) 2 Korneozyt
und der Anfangsbedingung u(x; y; 0) = 0.
Wir fuhren nun noch die relative Permeabilit¨at des Korneozyten und der Lipid-phase
" =Dcor
Dlip
mit " � 1;
den VerteilungskoeffizientenKcor=lip, fur den gilt:
Kcor=lip � ulip(x; y; t)jh� = ucor(x; y; t)jh+sowie die Lag-Zeit
Tlag =d2
6D
ein, wobei d die Dicke der Stratum-Corneum-Membran und D der Diffusionsko-effizient der Dgl. ist.
An den Phasenrandern ist der Flußkontinuierlich, d.h.:
Dliprulip(x; y; t) � ~n = Dcorrucor(x; y; t) � ~n
Beispiel: Fur Dcor = Dlip = 10�8 cm2
s, d.h. " = 1, und Membrandicke
d = 11�m ergibt sich eine Lag-Zeit von etwa Tlag = 20s.
2.2 Diskretisierung
Die Zeitdiskretisierung der AusgangsgleichungZ
du
dt=Z
Dru � ~n
erfolgt mit dem impliziten Euler-Verfahren:
du
dt=
u(t+ � )� u(t)
�+O(� ) ) u(t+ � )� � � divDru(t+ � ) = u(t)
) (I � � � divDr)u(t+ � ) = u(t)
Wissenschaftliches Rechnen Seite 23
Im Ort diskretisieren wir uber das Finite-Volumen-Verfahren:Zbi
duh
dt=Z@Bi
Druh � ~n =mXk=1
Z@Bi;k
Druh � ~n;
mit dem DiffusionskoeffizientenD =�1 in Lipidschicht" in Korneozyten
.
2.3 Wahl des Gitters
Problem: Berechnung vonR@Bi
Drbj � ~n fur eine nicht konstante
DiffusionsmatrixD.
Losung durch Aufspaltung in Integrale uber die Randteilstucke von B i:
Z@Bi
Drbj � ~n =mXk=1
Z@Bi;k
Drbj � ~n
laßt sich dann leicht auswerten, wenn D auf @Bi;k konstant ist.
) Wahle ein moglichst grobes Gitter (”Beschreibungs-“ /
”Ausgangsgitter“ )
so, daßD konstant auf @Bi;k ist:
- �
6
?
� -
�
1�m
ca. 15�m0; 1�m
anisotrope ElementeSeitenverhaltnis 150
p
p p
p
pp
p
p
6
@Bi;k
) (Khij) =
ZBi
bj � � �Z@Bi
Drbj � ~n ; i; j = 1; : : : ; n
liefert das Losungssystem:Kh � uh = fh
Wissenschaftliches Rechnen Seite 24
Zu Problemen fuhrt die starke Anisotropie einiger Gitterelemente mit Seiten-verhaltnis 150.
Die”rote“ Gitterverfeinerung:
Dreiecke: Verbinde die Mittelpunkte der Dreieckskanten
! vier ahnliche Teildreiecke! keine extremen Winkel !
Vierecke: Verbinde Seitenmittelpunkte
s
ss
s c
c
c
c c
! erhalt die Eigenschaften des Grundgitters (bei unserem Problem).
Approximationsfehler bei anisotropen Gitterzellen:
h � "
h
��u � h�2h�1 2 �1
i| {z }
=c�h2
+h�2"�2 �264 �12�1
375
| {z }=c�h2"2
D.h. die Approximation in y-Richtung ist um "2 besser als in x-Richtung,oder: die Approximation in x-Richtung ist um "�2 schlechter !
Das Seitenverhaltnis geht also quadratisch ein, denn:
Wissenschaftliches Rechnen Seite 25
Transformiere: � = x; � = 1"y
� =@2
@x2+
@2
@y2=
@2
@�2+
@2
@�2� "2
Der Ausweg liegt in der anisotropen, sog.”blauen Verfeinerung“ :
� -
� -
hi+1
hi
Fur die Anzahl der Gitterpunkte Ni im i-ten Verfeinerungsschritt gilt allerdings:Ni+1 = 2 � Ni, was zu großen linearen Gleichungssystemen Kx = b fuhrt unddamit Loser mit einem Aufwand der Großenordnug O(N) erfordert(!Mehrgitterverfahren).
2.4 Wahl eines geeigneten L¨osungsverfahrens
Qualitativ am Beispiel”Laplace“ : ��h = h�2
0B@ �1�1 4 �1
�1
1CA
) K�h = h�2
0BBBBBBBBBBBBBBBB@
4 �1 �1�1 4 �1 �1
�1 4 �1�1 4 �1 �1
�1 �1 4 �1 �1�1 �1 4 �1
�1 4 �1�1 �1 4 �1
�1 �1 4
1CCCCCCCCCCCCCCCCA
Problem: fill-in! direkter Loser geht nicht O(N 3), auch nicht fur Binder O(N 2),! iteratives Verfahren notig, da das LGS schon grob fehlerbehaftet.
x(0)! x(1)! : : :! x(i)! : : :! x ; Fehlertoleranzvorgabe
Die Iterationsfunktion � : x(i) ! x(i+1) sollte also eine Schrittzahl von O(1)haben, wobei jeder Schritt mit O(N) geht.
Wissenschaftliches Rechnen Seite 26
Im LGS Kx = b ist K eine schwach besetzte (n� n)-Matrix mit O(n) Eintragen6= 0 ( Gitterpunkte N = n2). Wir setzen fur das Iterationsverfahren
(2:1) K = M �R;
wobei M die angenaherte Inverseund R die Restmatrixist. M ist leicht inver-tierbar, d.h. M � y = z ist 8y; z 2 Rn mit einem Aufwand O(N) losbar und wirerhalten Konvergenz durch die Fixpunktiteration
x(i+1) = x(i) �M�1 � (Kx(i) � b)
� Jacobi-Verfahren (1849)!M ist Diagonalmatrix
M = D = diag(Kii); mit (2.1): K = D �R
x(i+1) = x(i) �D�1 � (Kxi � b)
Jacobi fur ��hx = f :
Gegeben: x(i) ; definiere Defekt: d(i) = Kx(i) � b
h�2(4x(i)j � x(i)j+1 � x
(i)j�1 � x
(i)j+n � x
(i)j�n) = d
(i)j + bj
Iterationsschritt:
x(i+1)j = x
(i)j �
1
4(4x
(i)j � x
(i)j+1 � x
(i)j�1 � x
(i)j+n � x
(i)j�n � h2bj)
=1
4(x
(i)j+1 + x
(i)j�1 + x
(i)j+n + x
(i)j�n| {z }
Mittelwert der Nachbarn
+h2bj)
Wir erstzen also xj durch den Mittelwert der Nachbarn (”glatten“).
Aber: Jeder Schritt hat Aufwand O(N), Anzahl der Schritte: O(N)O(N2) .
� Gauß-Seidel-Verfahren (1823/70)!M = L+D,wobei Lij = Kij 8i > j; Lij = 0 sonst, D = diag(Kij),sowie Rij = Kij 8i < j; Rij = 0 sonst.
� M = (L+D)D�1(R+D), wobei L = untere Dreiecksmatrix,D = Diagonalmatrix, R = obere Dreiecksmatrix, wie bei Gauß-Seidel.
� Young (1952): Ersetze D durch D! = ! �D mit Dampfungsfaktor !.
Wissenschaftliches Rechnen Seite 27
� Incomplete Lower-Upper-Decomposition (ILU):
�14
+!
+!
0BBBBBBBBBBBBB@
4 �1 �1�1 4 �1 �1
�1 4 �1 �1�1 4 �1
�1 4 �1�1 �1 4 �1
�1 �1 4 �1�1 �1 4
1CCCCCCCCCCCCCA!
!
0BBBBBBBBBBBBB@
4 �1 �10 4�1
4�1 �1
4�1
�1 4 �1 �1�1 4 �1
0 �1
44 � 1
4�1
�1 �1 4 �1�1 �1 4 �1
�1 �1 4
1CCCCCCCCCCCCCA! usw.
Idee: Wir lassen die zusatzlichen Eintrage weg und erhalten die“ Falsche L-U-Zerlegung“ :
LU = K +R(2:1)=) M = LU
Wir konnen schließlich eine Reduktion der hochfrequenten Fehleranteile erlan-gen, indem wir durch Intervallhalbierung das Gitter verfeinern und die Funktions-werte auf den neuen Gitterpunkten eintragen, die durch Mittelung der alten Werteentstehen.
Wissenschaftliches Rechnen Seite 28
2.5 Algorithmus des Mehrgitterverfahrens
mgm(l; u; f )f if (l == 0) lose Kl � u = f
elsef u = glatte (�1; u); /* z.B. mit Gauß-Seidel */
d = defekt(l; u; f );d = r � d; /* Restriktion des Defekts */v = 0; /* v ist Losung der Defektgleichung */for i = 1 to do mgm(l� 1; v; d);u = u� p � v;u = glatte (�2; u); /*) O(N) */
ggend;
Wir haben nun das LGS Kx = b mit dem Glatter S �(x) = �x, dem Defektd = K � �x� b und mit der Korrektur v : x = �x� v, denn:
Kx = K(�x� v) = K�x� kv|{z}=d
= K�x� d = K�x�K�x+ b = b
Programmpaket zur Computersimulation: uG = unstrukturierte Gitter! adaptive Mehrgitter auf Parallelrechner,! verschiedene Gitter, verschiedene ",! Vergleich von intra- und interzellularen Wegen.
2.6 Charakterisierung der Membran durch Lag-Zeit
Flußdurch den oberen Rand �0:
F�0(t) =ZT0
Ddu
dtd� = f1 +
Xj
�je��jt
Masse:
M�0(T ) =
TZ0
F�0(t)dt = f1T +m0 +Xj
�je��jt
Lag-Time aus Masse-Zeit-Diagramm (Nullstellen der Asymptote):
Tlag;1 = �m0
f1
Wissenschaftliches Rechnen Seite 29
m
tLag-Time
Bemerkung:Im Falle einer homogenen Membran ist Tlag;1 = L2
6D, mit effektiver Weglange L
(richtig skaliert).
Fur ein Gleichungssystem mit etwa 5800 Unbekannten errechnen wir folgendeWerte:
" Tlag;10 1h
10�5 10:61h10�2 0:44h1 20s
! Sprung bei " = 0
>−− 0ε Limit für
10 h
14 h
1 h
10-210-510-8 1
ε
Lagtime steady state
d.h. Korneozyten verlieren”Schwamm“ -Effekt, da
lim"!0
limt!1
Tlag;t 6= limt!1
Tlag;tj"=0 (nicht gleichmaßig stetig).
Wir brauchen also bei " = 0 ein anderes Simulationsgebiet mit inneren R¨andern,an denen wir unterschiedliche Konzentrationen zulassen.
Wissenschaftliches Rechnen Seite 30
Lag-Time fur eine Versuchsdauer von t = 50h:
10-5 10-2
1 h
1
ε
11 h
10-7
Lagtime nach t=50 h
Es sei nun A" die Steifigkeitsmatrix im Falle " > 0, d.h. bei gegebener Randbe-dingung fur uj@ und
du
dt= div Dru ; mit D =
�1; Lipidschicht"; Korneozyt
;
und B sei Steifigkeitsmatrix im Falle " > 0, also fur
du
dt= �u in 0 = n fKorneozyteng mit
@u
@~n= 0 auf @fKorneozyteng
Es gilt aber:A"
"!0�! A"=0 6= B
Unser Modell ist allerdings nur sinnvoll, wenn A" = B fur " ! 0. Gesucht istdaher eine Modifikation vonA" so, daßfur "! 0 automatisch die Randbedingungfur die inneren Rander entsteht.
Wir benotigen also im Modell eine zusatzliche Bedingung an den inneren Randernzwischen Lipidschicht und Korneozyt, die Transmissionsbedingung:
ulip = C � ucor auf @fKorneozyteng
Der Verteilungskoeffizient C liefert einen fest vorgegebenen Sprung uber denRand der Korneozyten und es gilt fur "! 0:
1
C! 0 bzw. C !1
Wissenschaftliches Rechnen Seite 31
� -
!
Bemerkung:
Tlag(";C; !) hangt bei Versetzungen ! 2 [ 14;12) der Korneozyten nur schwach von
! ab (< 20%).
Anderer Ansatz: Reduktion des Problems auf die Rander (Interfaces)
Es sei P ein innerer Punkt im Korneozyten/Lipidschicht und I ein Punkt desInterface:
K =
KPP KPI
KIP KII
!! KPP KPI
0 KII �KIPK�1PPKPI
!
mit dem Poincare-Steklov-OperatorKIPK�1PPKPI .
Wissenschaftliches Rechnen Seite 32
3 Numerik der Homogenisierung
Ein alternativer Zugang zum Hautproblem, gelesen von Dr. Nicolas Neuß.
3.1 Problemstellung
Lochermodell:
G
bzw. Navier-Stokes-GleichungenStrömung, beschrieben durch Stokes-,
Stokes-Gleichungen im Gebiet G":
��~u" +rp" = ~f (Impulserhaltung)
div ~u" = 0
Randbedingungen: z.B.: ~u" = 0 auf Lochrandern, ansonsten vorgegeben.
Offensichtliche Schwierigkeiten bei der Simulation des Problems:
� Gittergenerierung;
� Trotz vieler Unbekannter wegen feinem Gitter ist die Approximation ziem-lich schlecht, wenn man nicht noch weiter verfeinert, da in die Fehlerschatzungder Faktor 1
"eingeht;
� Die Losung des Problems kann Schwierigkeiten machen, wie z.B. die großenAnisotropien beim Hautproblem.
Zur Beseitigung der Schwierigkeiten stellen wir uns folgende Fragen:
� Muß ich exakt dieses Problem tatsachlich losen?
� Was weißder Anwender und was will der Anwender von mir wissen?
Wissenschaftliches Rechnen Seite 33
In unserem Falle lautet die Antwort meist:
� Der Anwender (Geophysiker, Mediziner, ...) kennt das Medium/Gebietnicht genau,
� Der Anwender kennt die Randbedingungen nicht genau,
� Der Anwender ist eigentlich nicht an der exakten Losung des Problems in-teressiert, sondern stattdessen etwa:
– Fließt genug Grundwasser nach, damit die Menge x abgepumpt wer-den kann?
– Wieviel Prozent eines ausgelaufenen Schadstoffs sind in der flussigenPhase, wieviel Prozent sind schon in der Bodenmatrix absorbiert?
– Wie ist der zeitliche Verlauf des Medikamentenzuflusses in den Blut-kreislauf?
Wir mussen daher Ersatzmodelle finden, die weniger Parameter benotigen undweniger detallierte Informationen liefern. Im Grunde ist diese Vereinfachungalltaglich, denn alle physikalischen Gesetze sind Approximationen von genaue-ren Beschreibungen. Man muß lediglich ein der jeweiligen Situation angepasstesModell verwenden.
Ersatzmodelle lassen sich beispielsweise finden durch:
� Experimente (z.B. das Darcy-Gesetz wurde experimentell entdeckt),
� physikalische Intuition, Erfahrung (”Es gibt nichts, was Physiker noch nicht
wissen.“ )
� mehr oder weniger rigorose Ableitung aus anderen Prinzipien. Dies ge-schieht oft im Nachhinein (z.B.: Kepler-Bahnen durch Gravitationstheorem,Darcy-Gesetz durch (Navier-)Stokes im periodisch porosen Medium).
Bezeichnung: Die Herleitung von physikalischen Gesetzen auf einer groberenSkala aus Gesetzen auf einer feineren Skala (Makro- bzw. Mikroskala) nennt manHomogenisierung. Dieser Teil der Mathematik beschaftigt sich mit Beziehungenzwischen Mikro- und Makrobeschreibungen physikalischer Prozesse.
Wissenschaftliches Rechnen Seite 34
3.2 Homogenisierung eindimensionaler Diffusionsgleichungen
Wir betrachten das Gebiet = ]0; L[ und die Schar a"; u", wobei der"-Parameter die Heterogenitaten beschreibt. Auf gilt:
�(a"(x)u"0(x))0 = f(x)
mit den Randbedingungen:
u"(0) = u0 ; u"(L) = uL
Losung durch Integration:
a"(x)u"0(x) = c1 �xZ0
f(�)d�
Setze G(x) := �xR0f(�)d�.
) u"0(x) =c1
a"(x)+
G(x)
a"(x)
, u"(x) = c2 +
xZ0
c1
a"(�)d� +
xZ0
G(�)
a"(�)d�
Die Konstanten erhalten wir durch die Randbedingungen:
c2 = u0
uL = u0 + c1
LZ0
d�
a"(�)+
LZ0
G(�)
a"(�)d�
) c1 =uL � u0 �
LR0
G(�)a"(�)d�
LR0
d�a"(�)
:
Damit ergibt sich die Losung:
u"(x) = u0 +uL � u0 �
LR0
G(�)a"(�)d�
LR0
d�a"(�)
�xZ0
d�
a"(�)+
xZ0
G(�)
a"(�)d�
Wissenschaftliches Rechnen Seite 35
Im Prinzip ist dies ausrechenbar durch Quadratur.
Wir betrachten nun den Spezialfall, daß a"(x) = a( �"), wobei a : R ! R+ eine
1-periodische Funktion ist mit a(x+ k) = a(x) 8k 2 Z; x 2 R.
Fall 1: f = O ) G = 0
Wir verifizieren außerdem auf der Periodizitatszelle Y = [0; 1]:
LZ0
d�
a( �")� L �
ZY
a�1(y)dy
| {z }=:<a�1>Y
bis auf Fehler der Ordnung "
In die Losung der Dgl. eingesetzt erhalten wir:
u"(x) = u0 +uL � u0
L < a�1 >Y�
xZ0
d�
a"(�)
Beispielsweise hat a folgende Gestalt:
-
6
x"
a
0,1
1
0,5 1 1,5 2 2,5
In diesem Falle ist < a�1 >Y=12 � 10 + 1
2 = 112 , und mit L = k" und den Rand-
bedingungen u0 = 0 und uL = 1 ergibt sich fur u" folgende Naherungslosung:
0
1
uε
xε 2ε 3ε 4ε L
Wissenschaftliches Rechnen Seite 36
Fur die Steigung von u" gilt:
u"0 =uL � u0
L� < a�1 >Y� 1
a(x")=
8><>:
2011L
, falls a(x") = 0:1
211L
, falls a(x") = 1
Fall 2: u0 = uL = 0
(3:1) u"(x) =1
L< a�1 >�1
Y �LZ0
G(�)
a( �")d� �
xZ0
d�
a( �")+
xZ0
G(�)
a"(�)d�
Lemma: (Oszillationslemma)
Sei ' 2 C10(Y = [0; L]) und f 2 L1(R) = Raum der integrierbaren Funktio-
nen mitR jf j <1 und f sei 1-periodisch. Dann gilt:
jLZ0
f(x
") � '(x)dx� < f >Y
LZ0
'(x)dxj � c"j'jC10
Beweisidee: Unterteile Y in Intervalle der Lange ".Spalte ' auf solch einem Intervall auf in '(x) +O(").
Aus dem Oszillationslemma erhalten wir
LZ0
G(�)
a( �")d� � < a�1 >Y �
LZ0
G(�)d�
| {z }=L�<G>[0;L]
und verwenden dies in Gleichung (3.1):
) u"(x) � < G >[0;L] �xZ0
d�
a( �")+
xZ0
G(�)
a"(�)d� =
xZ0
G(�)� < G >[0;L]
a( �")
d�
) u"(x) � < a�1 >Y �xZ0
(G(x)� < G >[0;L]) d�
Wissenschaftliches Rechnen Seite 37
Der Vergleich mit der Losung uhom zum homogenen Problem
ahom � u00hom = �f(x) ; ahom = const:
namlich:
uhom = a�1hom �xZ0
(G(�)� < G >[0;L]) d�
zeigt, daßu" � uhom , wenn wir ahom =< a�1 >�1Y wahlen.
Satz 3.1.:
u" approximiert fur "! 0 die Losung uhom des Problems
ahom � u00hom = �f(x) in ]L; 0[ ; ahom = const:
mit den Randbedingungen uhom(0) = u0 und uhom(L) = uL. Mit dem Oszil-lationslemma gilt:
ju" � uhom(x)j � c" 8x 2 [0; L] mit c = c0jf jC0[0;L]
Bemerkung: Aber es gilt nicht, daßu"0 ! u0hom !Um so eine Approximation zu erhalten, mußman einenKorrekturterm hinzufugen, der die Oszillation beinhaltet.
Eine Anwendung ist die Bestimmung der effektiven Leitfahigkeit einer Sandwich-Struktur in senkrechter Richtung zu den Schichten.
3.3 Der mehrdimensionale Fall
Im Gebiet � Rn haben wir die Gleichung:nX
i;j=1
@i(aij(x
")@j u
"(x) = f(x); wobei @i :=@
@xi; f 2 L2()
aij sei 1-periodisch mit �j�j2 � nPi;j=1
aij�i�j < �j�j2 , und �;� > 0.
Außerdem fordern wir periodische oder Dirichlet-Randbedingungen@ = (@)per [ (@)Dir, sowie auf (@)Dir:
u" = gj(@)Dir; g 2 H2() = fu integrierbar mit u; ru und @ 2u 2 L2g
Wissenschaftliches Rechnen Seite 38
Ansatz fur das homogenisierte Problem:
�@i(a0ij @ju0(x)) = f(x); in
Randbedingung: u0(x) = 0; fur x 2 @.
(aij) ist also konstante Diffusionsmatrix ohne Oszillation.
2.)
u = 0 u = 1
periodisch
periodisch1.)
u = 0 u = 1
periodisch
periodisch
schlechtgut
leitendleitend
Im Fall 1.) ist (ahom11 )�1 =< a�111 >Y , wenn a(Y )ij = Æij;im Fall 2.) gilt dagegen ahom11 =< a11 >Y .
Vorgehen: Lose zunachst die sog. Zellprobleme.
Fur k = 1; : : : ; n :Finde 1-periodische Korrekturfunktionwk : Y = [0; 1]n ! R mit
�@i(aij(y)@j(wk(y) + yk)) = 0
Um eine eindeutige Losung zu erhalten, normiert man etwa durchRYw(y)dy = 0 und setzt dann:
(3:2) a0kl =ZY
aij(y)@j(wl + yl)@i(wk + yk) =ZY
aij(y)@j(wl + yl)Æik
Bemerkung:
� Im Eindimensionalen stimmt die Definition (3.2) mit dem harmonischenMittel uberein.
Beweis: Suche fur Zellproblem w periodisch auf [0; 1) mit:(hier sei �0 := @
@y):
(a(w+ y)0)0 = 0 , a(w + y)0 = const: , aw0 + a = c
Wissenschaftliches Rechnen Seite 39
, w0 =c� a
a
Durch Integration uber [0; 1) folgt mit periodischem w:
0 =
1Z0
dy
a(y)�
1Z0
dy , c = (
1Z0
dy
a(y))�1
Damit heißt (3.2) im eindimensionalen Fall:
a0 =
1Z0
a(w + y)0dy =
1Z0
c dy = c
� Die Losung der Zellprobleme kann mit dem Mehrgitteralgorithmus berech-net werden. Es ist allerdings wegen den periodischen Randbedingungenempfehlenswert, auf groben Gittern die Losbarkeitsbedingung
RYd(y)dy;
(d = Defekt) durch Projektion zu erzwingen.
Satz 3.2.:
Sei u0 2 H2() \H10 () Losung zu
�@i(a0ij@ju0) = f; in
und u" 2 H10 () sei Losung zu
�@i(aij(x")@ju
"(x)) = f(x); in :
Die wk und a0ij seien wie oben definiert, ~w sei der Vektor der wk. Dann giltdie Fehlerabschatzung
kru" � (ru0 + "ru0 � ~w)kL2() � C"12ku0kH2()
Beweis: siehe Jikov-Kozlov-Oleinik.
Bemerkung:
� Wenn keine Randschichten auftreten oder wenn man Randkorrekturtermeeinbezieht kann man die Abschatzung : : : � C"ku0kH2() zeigen.
Wissenschaftliches Rechnen Seite 40
� Auch die Flusse aij(x")@ju"(x) konvergieren im Mittel gegen aij@ju
0(x).Ebenso gilt:
ku" � u0kL2() � C"ku0kh2()
Anwendung: z.B. im Hautproblem
Fur periodische Heterogenitaten ist im Prinzip ein Modell mit homogenisiertenKoeffizienten zulassig, d.h. anstatt
Periodizitätszelle
10 Schichten
reicht die Berechnung von zwei Zellproblemen einer Periodizitatszelle aus:
periodisch
periodisch
Symmetrien: Am zweidimensionalen Beispiel gilt fur a ij(y) = a(y):
a heißt
8><>:
1-symmetrisch, falls a(y1; y2) = a(1� y1; y2)2-symmetrisch, falls a(y1; y2) = a(y1; 1� y2)symmetrisch, falls beides
.
Definition 3.1.: Symmetrie
aij : Y ! R heißt k-symmetrischgenau dann, wenn gilt:
aij(y) = (�1)Æikaij (y1; : : : ; yk�1; 1� yk; yk+1; : : : ; yn)| {z }=:�y
(�1)Æik
Wissenschaftliches Rechnen Seite 41
Satz 3.3.:
Falls die aij k-symmetrisch sind, dann gilt 8k 6= l:
1:) wk(�y) = �wk(y) 2:) wl(�y) = wl(y)
Damit ist a0kj = a0ik = 0 fur i; j 6= k. Falls (aij) symmetrisch ist, wenden wirdies auf alle k an und erhalten dadurch eine Diagonalmatrix(a0ij).
Anwendung: stationares Hautproblem
a0ij@ju0 = 0; auf
Die Losung u0(x1; x2) =x2H
ist unabhangig von aij .Der Fluß laßt sich berechnen durch:Z
x2=0
a22 � 1
H|{z}=@22u0
� 1|{z}=n2
dx1 =a22
H
Bei exakter Periodizitat der Zellen ware a022 berechenbar.
Bemerkung: In der Praxis wird man a22 bzw. den asymptotischen Flußdurch eine Hautschicht mittels Messungen bestimmen:! Tabellen (je nach Patient verschieden).
Leider zeigen Beobachtungen aber auch, daßdas zeitabhangige Problem betrach-tet werden muß. Die Zeit bis zum Erreichen des stationaren Zustands liegt namlichje nach aij(
x") im Bereich von Stunden bis Tagen und kann daher nicht ver-
nachlassigt werden.
3.4 Homogenisierung von zeitabh¨angigen Problemen
Problem: Suche u" mit:
@tu"(x; t)� @i(aij(
x
")@ju
"(x)) = f(x; t); x 2 ; t > 0
Randbedingung: u"(x; t) = g(x; t); x 2 @; t > 0Anfangsbedingung: u"(x; 0) = �u(x); x 2
Wissenschaftliches Rechnen Seite 42
homogenisiertes Problem: Suche u0 mit:
@tu0(x; t)� @i(a
0ij@ju
0(x; t)) = f(x; t); x 2 ; t > 0
Randbedingung: u0(x; t) = g(x; t); x 2 @; t > 0Anfangsbedingung: u0(x; 0) = �u(x); x 2
Wieder gilt: u" approximiert u0 im Limes "! 0.Fehlerabschatzungen sind unter geeigneten Regularitatsvoraussetzungen und Vor-aussetzungen an aij ahnlich wie im elliptischen Fall (O(") in L2, mit KorrektorO("
12 ) fur den Gradienten, O("
12 ) fur die Flusse im Mittel).
Bemerkung: Die Vorraussetzungen an aij sind:
�j�j2 � aij�i�j � �j�j2; 8� 2 Rn
d.h. in die Fehlerabschatzung geht der Faktor ��
ein.
Trotz dieser Bemerkung versuchen wir die Anwendung auf das Hautproblem:
Aufgrund der Randbedingungen:
periodischperiodisch
konstant
konstant
Y
hangt die Losung nur von x2 ab. Wir erhalten ein eindimensionales, homogeni-siertes Problem mit einem relevanten Koeffizienten, namlich a22:
@tu0(x2; t) = a022@x2x2u
0(x2; t)
Wir schreiben im folgenden fur x2 = x, a22 = D und fur u0 = u. Dann lautet dasProblem:
@tu(x; t) = D � @xxu(x; t); x 2]0;H[; t > 0
Randbedingung: u(0; t) = 0; u(1; t) = 1; 8t > 0Anfangsbedingung: u(x; 0) = 0; x 2]0;H[
Wissenschaftliches Rechnen Seite 43
Bemerkung:
� Die Randbedingung ist inkompatibel mit der Anfangsbedingung, d.h. dietheoretische Behandlung in Funktionenraumen mit von t abhangiger Regu-laritat ist besonders schwierig.
� Bei der numerischen Behandlung interessiert der Ubergangsbereich (t klein)oftmals nicht. Mit einem impliziten Verfahren kann daher trotzdem mitgroßen Zeitschritten gerechnet werden.
Analytische Losung:
Schritt 1: Ziehe die stationare Losung xH
von u ab: v := u� xH
) @tv = D � @xxv
mit Randbedingungen: v(0; t) = v(H; t) = 0und Anfangsbedingung: v(x; 0) = � x
H
Schritt 2: Die Eigenfunktionen von D � @xx in [0;H] mit Null-Randbedingungensind die Funktionen
sin(k�
Hx); k 2 Z
Dies liefert uns einen Ansatz, der die Gleichung lost, namlich:
v(x; t) =1Xk=1
ck � sin(k�Hx) � expf�Dk2�2
H2tg
Die ck ergeben sich nun aus den Anfangsbedingungen:
� x
H=
1Xk=1
ck � sin(k�Hx)
Da 8k 6= l mit dem Skalarprodukt< x; y >L2(I):=RIx � y gilt:
< sin(k�
Hx); sin(
l�
Hx) >L2([0;H])= 0
erhalten wir damit fur die ck:
ck =< � x
H; sin(k�
Hx) >L2([0;H])
< sin(k�Hx); sin(k�
Hx) >L2([0;H])
Wissenschaftliches Rechnen Seite 44
Durch Rechnung finden wir
ck =2
k�(�1)k
und damit die Losung:
u(x; t) =x
H+
1Xk=1
(�1)k 2
k�� sin(k�
Hx) � expf�Dk2�2
H2tg
Es gilt also fur den Flußdurch x = 0:
F (t) = D � ddx
u(x; t)jx=0 = D � ( 1H
+1Xk=1
(�1)k 2
H� expf�Dk2�2
H2tg)
Definiere F1 := limt!1
F (t) = DH
(es bleibt zu zeigen: F (t)! 0 fur t! 0).
ooF
t
Fluß
τchar
Nach einiger Zeit ist der Term fur k = 1 dominant und die charakteristische Zeit�char lasst sich wie folgt aus dessen Nullstelle berechnen:
F (t) � D
H� (1 � 2expf�Dk2�2
H2tg)
) 1� 2expf�Dk2�2
H2�charg = 0
, �char =ln 2
�2� H
D|{z}=F�11
�H =ln 2
�2� HF1
Mit Zahlen (aus dem Paper von Heisig et al):
F1 = 4:3 � 10�9 cms
und H = 10�3cm ! �char � 4:5h
Wissenschaftliches Rechnen Seite 45
Berechnung der Lag-Zeit �1:
Sei M(t) die Menge der ins Blut geflossenen Substanz:
M(t) =
tZ0
F (t)dt
oo
ooF = Steigung derAsymptote
t
M(t)
τ t
Aus dem Steigungsdreieck ist ersichtlich:
M(t) = (t� �1)F1 , �1 =F1t�M(t)
F1
Approximiere nun F (t) durch
(0, falls t � �charF1 � (1� 2expf�D k2�2
H2 tg), falls t � �char.
) �1 =1
F1(F1
tZ0
1dt
| {z }=
�charR0
dt+tR
�char
dt
�F1tZ
�char
(1� 2expf�Dk2�2
H2tg)dt =
= �char + 2
tZ�char
expf�Dk2�2
H2tg = �char � 2
H2
D�2� expf�Dk2�2
H2tg| {z }
= 12
fur t=�char
jt großt=�char =
= �char| {z }ln 2 H
F1�2
+H2
D�2= (1 +
1
ln 2) � �char
Wissenschaftliches Rechnen Seite 46
Dies liefert uns nun eine Naherung der Lag-Zeit mit dem zuvor berechneten�char = 4:5h:
�1 = (1 +1
ln 2) � 4:5h � 11h
Unter Umstanden reicht dieser einfache Zusammenhang schon fur praktische An-wendungen aus. Leider gibt es aber auch Falle, in denen das Modell erweitertwerden muß, wie z.B. beim Hautproblem der Fall " fast 0.
Wissenschaftliches Rechnen Seite 47
4 Eigenschwingungen von Seen
Die Oberflachenschwingungen von Seen heißen”Seiches“ (frz.), speziell am
Genfer See. Wir wollen hier die Seiches des Bodensees betrachten.
Wissenschaftliches Rechnen Seite 48
4.1 Modell
Wir betrachten den See als ein abgeschlossenes Wasserbecken, d.h. uber den Rand� existiert kein Zu- oder Abfluß:
@u
@~n= 0; auf �
Die Einflusse der Temperatur sowie der Stromung des Wassers konnen durch dieNavier-Stokes-Gleichungenmodelliert werden, was allerdings zu einem extremkomplizierten System fuhren wurde. Naherungsweise ist der See in der Tie-fe durch die Wassertemperatur zweigeschichtet in das
”Hyperlimnion“ (obere,
warmere Schicht) und das”Hypolimnion“ (untere, kaltere Schicht).
H [m]
Temp. [°C]124.5
10
20
40
Interface
Die internen Seiches am Temperatur-Interface heißen beim Bodensee”Rinnen“
und”Anlaufen“ , ubernommen aus der traditionellen Fischersprache, in der diese
Begriffe die auftretenden Oberflachenschwingungen beschrieben.
Die Stromung verursacht bei großen Geschwindigkeiten Turbulenzen, wahrenddas Gewasser bei kleinen Stromungsgeschwindigkeiten laminar (durchsichtig)bleibt. Die Turbulenzen sind Verwirbelungen im �m-Maßstab, besonders hoch-gradig in einer Grenzschicht am Seeboden mit einer Dicke im m-Bereich, wel-che durch den entstehenden Energieentzug zu einer Selbstbremsung der Stromungfuhrt.
Es entstunde dadurch bei der Flachenausdehnung des Bodensees uber etwa 120�30 km ein linearer Mehrskaleneffekt, weshalb wir zur Vereinfachung des Modellssowohl die Stromung vernachlassigen, als auch das Problem als eingeschichtet,also ohne den Einflußdes Hypolimnions, betrachten.
Wissenschaftliches Rechnen Seite 49
4.2 Herleitung der Flachwassergleichungen
Wir haben ein Kraftegleichgewicht zwischen den beteiligten Kraften:
� Tragheit: FT
� Corioliskraft: FC
� Gravitation: FG
� Druckkrafte: FD
und bekommen dadurch eine Gleichung fur die Impulserhaltung:
(4:1) FT + FC + FG + FD = 0
Konstitutive Beziehung:
Mit der Gravitationskonstante g und der Wasserdichte % gilt fur dieGravitationskraft:
(4:2) FG = �g �ZV
% dV �0B@ 0
01
1CA
Tragheit:
(4:3) FT =ZV
%@~u
@tdV
Corioliskraft mit Coriolisparameter f :
(4:4) FC = f �ZV
% � ~u�0B@ dx
dy
dz
1CA
Druckkrafte aus Taylorentwicklung des Drucks p = p0+dx �rp+ Terme hohererOrdnung:
(4:5) FD =ZV
rp dV
Durch Einsetzen der Gleichungen (4.2) bis (4.5) in (4.1) erhalten wir mit~u = (u; v; w)T :
ZV
%@~u
@tdV + f �
ZV
% � (u; v; 0)T �0B@ dx
dy
dz
1CA� g � Z
V
% dV �0B@ 0
01
1CA+
ZV
rp dV = 0
Wissenschaftliches Rechnen Seite 50
Daraus ergibt sich ein System von Bewegungsgleichungen (4.6) in den Einzelko-ordinaten, bestehend aus drei Gleichungen mit vier Unbekannten:
%@u
@t+ f%v +
@p
@x= 0 (i)
%@v
@t+ f%u+
@p
@y= 0 (ii)
%@w
@t+ g%+
@p
@z= 0 (iii)
Aus Gleichung (4.6)(iii) folgt durch Integration nach z mit @w@t
= 0, der Wasser-tiefe H = H(x; y) und der Oberflachenauslenkung � = �(x; y; t):
(4:7) p(z) = p0 + g% � (� +H � z); 0 � z � H
(x,y,t)ξ
x,y
z
H(x,y)
Mit der Annahme u(x; y; z; t) = u(x; y; t) folgt fur die Gleichungen (4.6.i) und(4.6.ii) durch Integration in z-Richtung von 0 bis H:
(4:8:i)
HZ0
@u
@tdz +
HZ0
fv dz +
HZ0
1
%
@p
@xdz = 0
(4:8:ii)
HZ0
@v
@tdz +
HZ0
fu dz +
HZ0
1
%
@p
@ydz = 0
Definiere U :=HR0u dz und V :=
HR0v dz und setze dies mit Gleichung (4.7) in
(4.8.i) ein:
0 =@U
@t+ f � V +
HZ0
1
%� (p0 + g% � (� +H � z)
@x)dz =
Wissenschaftliches Rechnen Seite 51
=@U
@t+ f � V +
1
%� @(p0 + g% � ((� +H)z � z2
2))jH0
@x=
=@U
@t+ f � V + gH � @�
@xmit
@%
@x= 0
Damit erhalten wir, analog fur Gleichung (4.8.ii), das System fur die Bewegungs-gleichungen (4.9) mit zwei Gleichungen fur drei Unbekannte:
@U
@t+ f � V + gH � @�
@x= 0
@V
@t+ f � U + gH � @�
@y= 0
Massenerhaltung:
ZV
@%
@tdV =
Z@V
~j � ~n d� =ZV
div~j dV
Der Ansatz: ~j = c � ~u (konstitutive Beziehung) liefert:ZV
@%
@tdV = c
ZV
div ~u dV
Aus der Inkompressibiltat folgt, daß @%@t
= 0:
, div ~u = 0
Integriere nun div ~u = @u(x;y;z;t)@x
+ @v(x;y;z;t)@y
+ @w(x;y;z;t)@z
= 0 von 0 bis H mit
~u = (u; v; w) und den Randbedingungen w(0) = @�@t
, w(H) = 0:
0 =
HZ0
div ~u =@U
@x+@V
@y+
HZ0
@
@zw dz =
@U
@x+@V
@y+ wjH � w0
) (4:10)@U
@x+@V
@y� @�
@t= 0
Aus den Bewegungsgleichungen (4.9) und aus (4.10) erhalten wir nun das dreidi-mensionale Differentialgleichungssystem:
(4:11)
0BBB@
@@t
f gH @@x
�f @@t
gH @@y
@@x
@@y
�
@@t
1CCCA �
0BBB@U
V
�
1CCCA = 0
Wissenschaftliches Rechnen Seite 52
Zeitabhangigkeit: periodische Losung in der Zeit
U(x; y; t) = U 0(x; y) � ei!t
V (x; y; t) = V 0(x; y) � ei!t
�(x; y; t) = �0(x; y) � ei!t
in (4.11) eingesetzt:
(4:12)
0BBB@
i! f gH @@x
�f i! gH @@y
@@x
@@y
�i!
1CCCA �
0BBB@U 0
V 0
�0
1CCCA = 0
Ab sofort lassen wir die Striche weg, woduch die Zeit eliminiert und durch diePeriodizitat ersetzt wird.
Fur die Coriolisparameter gilt am Bodensee: f � 10�5, weshalb wir die Coriolis-krafte ebenfalls vernachlassigen konnen und damit aus (4.12) folgendes Systemerhalten:
(4:13)
0BBB@i! 0 gH @
@x
0 i! gH @@y
@@x
@@y
�i!
1CCCA �
0BBB@U
V
�
1CCCA = 0
Durch Blockelimination eliminieren wir nun die ersten beiden Eintrage der letztenZeile und fuhren so die Systemmatrix uber in:
0B@ i! f gH @
@x
�f i! gH @@y
0 0 �i! + i!div gHr�
1CA
Aus der dritten Zeile des neuen Systems erhalten wir schließlich ein Eigenwert-problem fur das eingeschichtete Becken:
(4:14) (div gHr� !2)� = 0; in
div (gHr) ist elliptisch, falls H > 0, d.h. hat Eigenwerte > 0.Die Randbedingung lautet:
@�
@~u= 0 auf @:
Wissenschaftliches Rechnen Seite 53
4.3 Dreiecksgitter
Aufgrund der komplizierten Geometrie (Rand, Inseln) verwenden wir im Zu-ge des Mehrgitterverfahrens ein unstrukturiertes Gitter aus Dreiecken mit einemmoglichst groben Ausgangsgitter 0 (Bodensee: 88 Punkte), welches anschlie-ßend verfeinert wird.
Ausgangsgitter 0 und erster Verfeinerungsschritt:
Das folgend abgebildete Gitter (4337 Punkte) ist das feinste einer Folge von funfGittern, die durch fortgesetzte Verfeinerung aus dem groben Anfangsgitter 0
entstehen.
Wissenschaftliches Rechnen Seite 54
Wir benotigen hierzu zwei verschiedene Verfeinerungsmethoden, die sog.
”rote“ Verfeinerung und den
”grunen Abschluß“:
grün
rot
Bei der regularen, roten Verfeinerung erhalten wir kongruente Teildreiecke und esbleiben die Winkel erhalten, wahrend der grune Abschlußdie Winkel verschlech-tert, fur die Konsistenz des Endgitters aber vonnoten ist, d.h. es durfen keineisolierten Eckpunkte auf den Dreieckskanten entstehen.
Beispiel:
0 ! 1 ! 2 ! 3
Wegen der in einer Ecke konzentrierten Verfeinerung hat nun allerdings das Grund-gitter keinen gleichverteilten Fehler mehr, weshalb wir erneut
”rot“ verfeinern
mussen. Dazu entfernen wir den vorher gesetzten grunen Abschluss, um keineschlechteren Winkel zu erhalten.
Das Gitter 1 in obigem Beispiel sieht dann beispielsweise folgendermaßen aus:
Wissenschaftliches Rechnen Seite 55
Projektion auf den Rand:
Zur Verbesserung der Randauflosung werden im Zuge der Verfeinerung solchePunkte, die auf an den Rand grenzende Dreieckskanten fallen, auf die Randkurveselbst projiziert. Dabei wird der neue Randgitterpunkt so gewahlt, daßdie Winkelder neu entstehenden Dreiecke optimal sind.
Dies ist moglich, falls die Randkurve konvexist.
Bei nicht konvexen Randern dagegen konnen Probleme auftreten, indem die Drei-ecke
”umklappen“ , was das Gitter unbrauchbar macht.
Um dies zu verhindern, muß bereits das Grundgitter auf problematische Rand-punkte angepasst werden.
4.4 Das Algebraische Mehrgitterverfahren (AMG)
Auf dem gewonnenen Gitter diskretisieren wir nun die Differentialgleichung (4.14)mit Hilfe der Finiten Elemente:
�h 2 Vh := fvjT linear, v stetig auf g
Wissenschaftliches Rechnen Seite 56
Suche � 2 H1, fur das 8 Testfunktionen v 2 H1 gilt:Z
(div(gHr�)v � !2�) � v = 0
)Z
div(gHr�)v � !2Z
�v = 0
Durch partielle Integration folgt mit dem Gaußschen Integralsatz:Z
gHr�rv �Z@
gHr� � ~n| {z }= @�
@~n=0
�v � !2Z
�v = 0
)Z
gHr�rv � !2Z
�v = 0
Mit Basisdarstellung � =nPi=1
�ibi folgt 8j = 1; : : : ; n :
nXi=1
(Z
gHrbi � rbj| {z }
=:aij
) � �i � !2nXi=1
�i
Z
bi � bj| {z }=:mij
= 0
Dies liefert:Ax� !2Mx = 0
mit der SteifigkeitsmatrixA := (aij)ni;j=1 =
RgHrbi � rbj
!ni;j=1
und der Massenmatrix M := (mij)ni;j=1 =
Rbi � bj
!ni;j=1
.
Lumping-Verfahren: Wir betrachten nun die Massenmatrix M
1-dimensional:
0i-2 i-1 i
(bi; bj)0 =
8>>>>><>>>>>:
Rb2i , falls i = j
Rbibj 6= 0 , falls ji� jj = 1
0 , falls ji� jj � 2
Wissenschaftliches Rechnen Seite 57
D.h. M ist Tridiagonalmatrix, symmetrisch und positiv definit, denn mitmij =
Rbibj = (bi; bj)0 gilt:
xTMx =nXi=1
(�ibi; �jbj)0 =~�T ~� > 0; 8� 6= 0
Daraus ergibt sich fur die Eigenwerte �i von M :
0 < �0 < �min(M)| {z }!1(h!0)
� �max(M)| {z }!1(h!0)
< �1 <1
Die Massenmatrix M ist also”spektral aquivalent“ zur Identitat und kann daher
durch eine Diagonalmatrix ersetzt werden. Wir wahlen als Diagonalmatrix:
D = diag(dii) mit dii =nXj=1
jZ
bibjj
und erhalten schließlich das zu losende Eigenwertproblem:
(4:15) Ax� !2Dx = 0
Methoden zur Losung von Eigenwertproblemen:
QR-Verfahren:
Definiere mit regularer Matrix Q und oberer Dreiecksmatrix R:
A =: QR; sowie A(1) := RQ; und A(1) =: Q(1)R(1)
Daraus erhalten wir die Iterationsvorschrift:
A(i) := R(i�1)Q(i�1); mit A(i) =: Q(i)R(i)
Dieses direkte Verfahren konvergiert kubisch und liefert alle Eigenwerte. Es istaber aufgrund des großen Aufwands nur sinnvoll fur
”kleine“ Matrizen
(n � 1000).
Das QR-Verfahren ist also allenfalls als Grobgitterloser anwendbar, womit wireinmalig die Eigenwerte des Ausgangsgitters berechnen konnen.
Vektoriteration (v.Mises):
Die SteifigkeitsmatrixA ist symmetrisch und positiv definit; �i seien die n Eigen-
vektoren zu den Eigenwerten �i. Mit x =nPi=0
�i�i gilt:
Ax = A �nXi=1
�i�i =nXi=1
�iA�i =nXi=1
�i�i�i
Wissenschaftliches Rechnen Seite 58
Der Rayleigh-Quotient
�i =�iA�i
�Ti �i
liefert uns also den großten Eigenwert, da die restlichen Koeffizienten verschwin-den. Da wir aber zunachst den kleinsten Eigenwert suchen, betrachten wir dieInverse von A.
Inverse Iteration (Wielandt):
A�1x(i) = x(i+1) liefert kleinsten Eigenwert �min(A)
Mit einer Schatzung ~� fur den zweitkleinsten Eigenwert aus dem Rayleigh-Quotientenbetrachten wir die geshiftete Matrix
(A� ~�I)�1x(i) = x(i+1)
.........................................................................� -
.........................................................................
~�
�3�2�10
�1 � ~� �2 � ~� �3 � ~�
Die Iteration strebt gegen �i, falls gilt:
j~�i � �ij = minjj~�i � �j j
Sukzessive Inverse Iteration liefert zunachst ! �1; �1 ,dann mit AjRnnf�1g ! �2; �2 usw.
Realisiert wird der Ubergang x 2 Rn ! x0 2 Rn n f�g durch eine Projektion:
x0 = x� xT�
�T�� �
Angewandt auf Gleichung (4.15) konvergiert also !2i ! !2 und xi ! x, wobei
xi immer mehr im Kern liegt. Bei der Inversen Iteration mit Shift enstehen daherimmer
”singularere“ , d.h. schlechter konditionierte Gleichungssysteme Ax = b,
welche nur losbar sind, wenn b 2 Bild(A) = Kern(AT )?. Die rechten Seitenerfullen diese Kompatibilitatsbedingung immer weniger. Die Anwendung vonMehrgitterverfahren als Loser der Inversen Iteration hat daher nur Sinn, wenn einfester Shift � verwendet wird und � 6! !2
i .
Wissenschaftliches Rechnen Seite 59
4.5 Mehrgitterverfahren f ur Eigenwertprobleme
Wir betrachten das zu Gleichung (4.15) ahnliche Modellproblem:
(K � �I)� = f; K = KT
Falls gilt, daß f ? Kern(K) , so sind das Jacobi- bzw. das Gauß-Seidel-Verfahren anwendbar:
Jacobi: (mit Dampfungsfaktor �)
�neu = �alt � �D�1 � ((K � �I)�alt � f)
Gauß-Seidel: (mit K � �I = D + L+ U )
�neu = �alt � (D + L)�1 � ((K � �I)�alt � f)
Falls �alt ? Kern(K��I), so glatten das Gauß-Seidel bzw. das Jacobi-Verfahren.Falls �alt 2 Kern(K � �I) (d.h. �alt Eigenvektor von K zum Eigenwert �), soliefert �alt keinen Beitrag zum jeweiligen Verfahren.�alt ? Kern(K��I) mußalso nicht notwendig gelten und wir haben damit einenGlatter.
Definiere nun den Defekt:
(K � �I)�neu � f =: d = d1 � d2
wobei d1 := (K � �I)vl mit Korrekturvektor vl und d2 2 Kern(K � �I).
Wir stellen durch eine Projektion Qe sicher, daß nur der Anteil d1 weitergegebenwird. e sei Eigenvektor von K zum Eigenwert �:
Qe = � � (�T � e)e
Zweigitterverfahren fur EWe und EVen von K:
~�i;l :=�i;lKl�i;l�Ti;l�i;l
/* Rayleigh-Quotient, l=Gitterlevel */��i;l := S�
��i;l(�i;l) /* Glatter */
dl�1 := r(Kl � ~�i;lI)��i;l /* Defekt mit Restriktion r */d?l�1 := Qel�1dl�1 /* Projektion */Lose (Kl�1 � ~�i;lI)vl�1 = d?l�1�i+1;l = ��i;l � pQel�1 vl�1 /* mit Prolongation p */
Wissenschaftliches Rechnen Seite 60
Hier gilt nun, daßd?l�1 ? Kern(K � ~�I), also ist das Problem gutartig.
Die Eigenvektoren el�1 erhalten wir durch geschachtelte Iteration(nestediterations). Dazu berechnen wir zuerst die EVen auf dem grobsten Gitter mittelsQR-Zerlegung o.a. und speichern diese ab. Beim Ubergang auf das nachstfeinereGitter wird auf diese dann zuruckgegriffen (! e l�1).
Algorithmus des Mehrgitterverfahrens fur Eigenwertprobleme:
emgm(l; �; f )integer l;gitterfunktion �;beginf if (l == 0) lose Kl � � = f direkt
elsef � = �TKl�
�T �; /* Rayleigh-Quotient */
� = S��(�); /* glatten */
dl�1 = Qel�1r(Kl � �I)�; /* Berechnung des Defekts */v = 0; /* v ist Korrektur */for j = 1 to do smgm(l� 1; v; d);� = � � pQl�1 v;
ggend;
MGM fur singulares Problem:
smgm(l; v; f )beginf if (l == 0) lose Kl � v = f direkt
elsef v = S�
l (v); /* glatten */dl�1 = Qel�1r(Kl � �I)v; /* Berechnung des Defekts */w = 0; /* w ist Korrektur */for j = 1 to do smgm(l� 1; w; dl�1);v = v � pQl�1 w;
ggend;
Wissenschaftliches Rechnen Seite 61
Ergebnis: Die Eigenschwingungen des Bodensees
Benutzter Loser:
� Mehrgitterverfahren fur Eigenwertprobleme (emgm)
� ILU - Glatter (vgl. Kap.2.4)
Die Eigenwerte des Bodensees (in Minuten):
Nr. beobachtet EMGM1 55,64 56,662 37,72 40,053 28,28 27,834 19,13 20,745 15,196 15,007 15,00 13,718 14,60 13,209 12,42
10 11,60 11,6511 11,30 11,30
Wir vergleichen schließlich die Eigenwerte des Bodensee-Problems bzw. die zu-gehorigen Eigenschwingungsfrequenzen mit den Schwingungen einerInstrumentensaite:
Eigenwert Frequenz Faktor Ton exakt56,66 440,00 1,00 a 440,0037,72 660,93 1,50 e 659,2628,28 881,56 1,33 a’ 880,0019,13 1303,21 1,48 e’ 1318,5115,00 1662,03 1,28 gis 1661,2214,60 1707,56 1,03 gis+11,60 2149,17 1,26 c+ 2093,0011,30 2206,23 1,03 cis