Danksagungen

An allererster Stelle m ochte ich Dr. Martin Hanke v om Institut f ur Praktische Mathematik der Universit at Karlsruhe sehr herzlich d a n k en, der durch s e i n e Hilfe sehr viel zur L osung des gestellten Problems beigetragen hat.

Mein Dank geht ferner an Patricia K. Lamm von der Michigan State University, die Gr underin und Organisatorin des IPNet, durch deren Hilfe der Kontakt uberhaupt erst zustande kam. zu Dr. Hanke

Weiterhin danke i c h Dipl.-Ing. Martin Haas vom Institut f ur Theorie der Elektrotechnik der Universit at Stuttgart, der mir zu Beginn der Arbeit einige wertvolle Hinweise geben konnte.

Schlielich danke i c h meinem Betreuer, Dipl.-Ing. Georg F assler, allen sonstigen Mitarbeitern und Studenten am Institut und insbesondere der UNO-Runde f ur die sehr angenehme und entspannte Atmosph are, in der die Arbeit durchgef uhrt werden konnte.

1

Inhaltsverzeichnis

  Danksagungen  1

  Inhaltsverzeichnis  2

  Legende  5

1  Einleitung  7

2  Herleitung der Matrixgleichungen  9

2.1  Ansatz und Einf uhrung der Bezeichnungen : : : : : : : : : : : :  9

2.2  Maxwellsche Gleichungen und Potentiale : : : : : : : : : : : : : :  10

2.3  Berechnung des Vektorpotentials  

A  : : : : : : : : : : : : : : : : :  12

2.4  Berechnung der magnetischen Feldst arke  

H  : : : : : : : : : : : :  13

3  Berechnung des Matrix-Vektor-Produkts  15

3.1  Der eindimensionale Fall : : : : : : : : : : : : : : : : : : : : : : :  15

3.1.1  Die Struktur der Matrix : : : : : : : : : : : : : : : : : : :  15

3.1.2  Die zyklische diskrete Faltung : : : : : : : : : : : : : : : :  16

3.1.3  Berechnung des Matrix-Vektor-Produkts : : : : : : : : : :  17

3.2  Der zweidimensionale Fall : : : : : : : : : : : : : : : : : : : : : :  18

3.2.1  Die Struktur der Matrix : : : : : : : : : : : : : : : : : : :  18

3.2.2  Die zweidimensionale zyklische diskrete Faltung : : : : : :  19

3.2.3  Berechnung des Matrix-Vektor-Produkts : : : : : : : : : :  20

3.3  Matrix-Vektor-Multiplikation in Bl ocken : : : : : : : : : : : : : :  21

4  Fredholmsche Integralgleichungen 1. Art  23

4.1  Einf uhrung : : : : : : : : : : : : : : : : : : : : : : : : : : : : : :  23

4.2  Allgemeine Schwierigkeiten schlecht gestellter Probleme : : : : :  25

4.2.1  Die Singul arwertentwicklung : : : : : : : : : : : : : : : :  25

2

INHALTSVERZEICHNIS

   3

4.2.2  Die gl attenden Eigenschaften von K : : : : : : : : : : : :  26

4.2.3  Die Picard-Bedingung : : : : : : : : : : : : : : : : : : : :  27

4.3  Analyse der Koezientenmatrix : : : : : : : : : : : : : : : : : :  28

4.3.1  Die Singul arwertzerlegung : : : : : : : : : : : : : : : : : :  28

4.3.2  Beziehungen zwischen der SVD und der SVE : : : : : : :  28

4.3.3  SVD-Analyse : : : : : : : : : : : : : : : : : : : : : : : : :  30

4.4  Regularisierung und Filterung : : : : : : : : : : : : : : : : : : : :  31

4.4.1  Regularisierung : : : : : : : : : : : : : : : : : : : : : : : :  31

4.4.2  Das Tikhonov-Verfahren : : : : : : : : : : : : : : : : : : :  32

4.4.3  Abgeschnittene Singul arwertzerlegung : : : : : : : : : : :  33

4.4.4  Iterative V erfahren : : : : : : : : : : : : : : : : : : : : : :  33

4.4.5  Geeignete Diskretisierung : : : : : : : : : : : : : : : : : :  34

4.5  Analyse regularisierter Probleme : : : : : : : : : : : : : : : : : :  35

4.5.1  Die verallgemeinerte Singul arwertzerlegung : : : : : : : :  35

4.5.2  Wichtige Beziehungen der GSVD : : : : : : : : : : : : : :  36

4.5.3  Untersuchungen mit der GSVD : : : : : : : : : : : : : : :  37

4.6  Die L-Kurve : : : : : : : : : : : : : : : : : : : : : : : : : : : : : :  38

4.6.1  Die Denition der L-Kurve : : : : : : : : : : : : : : : : :  39

4.6.2  Der " Knick in der L-Kurve : : : : : : : : : : : : : : : : :  39

4.6.3  Die Wahl des Regularisierungsparameters : : : : : : : : :  41

5  Das Programm CGFFT  44

5.1  Programmablauf : : : : : : : : : : : : : : : : : : : : : : : : : : :  44

5.2  Externe Unterprogramme : : : : : : : : : : : : : : : : : : : : : :  45

5.2.1  Basic Linear Algebra Subprograms : : : : : : : : : : : : :  45

5.2.2  LAPACK : : : : : : : : : : : : : : : : : : : : : : : : : : :  46

5.2.3  Transactions on Mathematical Software : : : : : : : : : :  47

5.2.4  Golden Oldies : : : : : : : : : : : : : : : : : : : : : : : : :  47

5.3  Compilierung des Programms : : : : : : : : : : : : : : : : : : : :  47

5.4  Das Hilfsprogramm DAT2INP : : : : : : : : : : : : : : : : : : : :  47

6  Ergebnisse  49

7  Versuche mit anderen Ansatzfunktionen  54

7.1  Das erste Verfahren : : : : : : : : : : : : : : : : : : : : : : : : :  54

7  2 Das zweite Verfahren : : : : : : : : : : : : : : : : : : : : : : : : :  57

INHALTSVERZEICHNIS

   4

8  Zusammenfassung und Ausblick  59

A  Informationen und Programme aus den Datennetzen  61

  A.1 Informationen : : : : : : : : : : : : : : : : : : : : : : : : : : : : :  61

  A.1.1 IPNet : : : : : : : : : : : : : : : : : : : : : : : : : : : : :  61

  A.1.2 NA-NET : : : : : : : : : : : : : : : : : : : : : : : : : : :  62

  A.1.3 Usenet Network News : : : : : : : : : : : : : : : : : : : :  62

  A.1.4 World Wide Web und Mosaic : : : : : : : : : : : : : : : :  63

  A.1.5 Electronic Transactions on Numerical Analysis : : : : : :  64

  A.1.6 Verzeichnisse mit Forschungsberichten : : : : : : : : : : :  64

  A.2 Programme : : : : : : : : : : : : : : : : : : : : : : : : : : : : : :  65

  A.2.1 netlib : : : : : : : : : : : : : : : : : : : : : : : : : : : : :  65

  A.2.2 Guide to Available Mathematical Software : : : : : : : :  65

  A.2.3 Archie : : : : : : : : : : : : : : : : : : : : : : : : : : : : :  66

B  N utzliche Programmierhilfen  67

  B.1 Emacs FORTRAN-Mode : : : : : : : : : : : : : : : : : : : : : :  67

  B.2 ftnchek : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : :  68

  B.3 s2d und d2s : : : : : : : : : : : : : : : : : : : : : : : : : : : : : :  69

  B.4 xfig : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : :  69

  B.5 Unigraph : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : :  70

C  Erzeugung der ben otigten Libraries  74

  C.1 Erzeugung des Programms s2d : : : : : : : : : : : : : : : : : : :  74

  C.2 Erzeugung der Library fft : : : : : : : : : : : : : : : : : : : : :  75

  C.3 Erzeugung der Library 691 : : : : : : : : : : : : : : : : : : : : :  77

D  Programm-Quelltexte  80

  D.1 Programm cgfft.f : : : : : : : : : : : : : : : : : : : : : : : : : :  80

  D.2 Programm dat2inp.f : : : : : : : : : : : : : : : : : : : : : : : :  97

  D.3 Matlab-Routine cgne : : : : : : : : : : : : : : : : : : : : : : : : :  99

  Literaturverzeichnis  101

Legende

Darstellungskonventionen

reelle skalare Gr oe ^F komplexe skalare Gr oe F Ortsvektor

^{~ r} komplexer Vektor eines linearen Gleichungssystems ^{~ F} x

A reelle Matrix A komplexe Matrix

^ Element einer periodischen Folge ^x n

^ ^X n diskrete Fouriertransformierte einer periodischen Folge arg É Argument einer komplexen Gr oe, arg ^F = = f ln F g

fÉg

=fÉg

j É j

jj É jj

É É

ij

5

LEGENDE 6

Liste der verwendeten Formelzeichen

a Abstand der Mepunkte

! Kreisfrequenz, ! = 2 f

Kapitel 1

Einleitung

Das Ziel dieser Diplomarbeit ist die Entwicklung eines numerischen Algorithmus', der es erlaubt, aus den oberhalb einer Leiterplatte gemessenen magnetischen Feldern die auf der Leiterplatte ieenden Str ome zu berechnen. Dieser Algorithmus soll anschlieend durch e i n F ORTRAN-Programm implementiert werden.

Um eine eindeutige L osung f ur das Problem zu erm oglichen, mu die Annahme getroen werden, da Str ome ausschlielich in der Platinenebene ieen. uberbestimmte lineare Gleichungssysteme, die so zu l osen sind, da die Summe der Fehlerquadrate minimiert wird. Auerdem mu ein Regularisierungsverfahren angewandt werden, um die numerische Instabilit at auszugleichen, die bei der L osung einer Integralgleichung erster Art regelm aig auftritt 14.

Die Gleichungssysteme werden mit der von R. F. Harrington 18 e n twickelten Momentenmethode aufgestellt. Dabei stellt man die Str ome als eine Summe von geeignet gew ahlten Ansatzfunktionen dar. Man fordert dann von den durch die Str ome verursachten elektromagnetischen Feldern die Erf ullung bestimmter Bedingungen, n amlich da ein Skalarprodukt mit einer sogenannten Testfunktion bestimmte Werte annimmt. Das sich so ergebende lineare Gleichungssystem stellt eine diskretisierte Version der zu l osenden Integralgleichung dar.

Es zeigt sich, da die das Gleichungssystem repr asentierende Matrix eine besondere Struktur hat, die man als Block-Toeplitz-Toeplitz-Block-Matrix bezeichnet. Diese Struktur ergibt sich h aug bei der Diskretisierung zweidimensionaler Probleme, denen eine Integralgleichung mit einem Verschiebungskern zugrunde liegt. Eine kennzeichnende Eigenschaft dieser Matrixstruktur ist es, da sich Produkte der Matrix mit einem Vektor sehr zeit- und speicherplatzsparend mit Hilfe einer zweidimensionalen diskreten Fouriertransformation berechnen lassen.

Damit bieten sich iterative V erfahren, wie etwa das Verfahren der konjugierten Gradienten, angewandt auf die Normalengleichungen CGNR, zur L osung des Gleichungssystems an. Dieses Iterationsverfahren hat, wie viele andere auch, auerdem den Vorteil, da es gleichzeitig regularisierend wirkt, es gen ugt dazu,

7

KAPITEL 1. EINLEITUNG 8

die Iterationen an einer geeigneten Stelle abzubrechen. Es ergibt sich, da mit einem solchen Verfahren auch relativ groe Gleichungssysteme in einer ange- messenen Zeit mit einem m aigen Speicherplatzbedarf l osbar sind.

Kapitel 2

Herleitung der

Matrixgleichungen

2.1 Ansatz und Einf uhrung der Bezeichnungen

F ur die Berechnung der auf der Leiterplatte ieenden Str ome wird die Leiterunnen Dr ahten angen ahert.

Die in diesem Gitter ieenden Str ome werden in den einzelnen Segmenten jeweils als konstant angenommen als Segment wird hier ein zwischen zwei Kreuzungspunkten des Gitters liegendes Drahtst uck bezeichnet. Es wird immer nur eine Frequenz betrachtet.

Da alle vorkommenden Gr oen zeitharmonisch sind, empehlt sich die Verwendung von komplexen Gr oen f ur die Rechnung. Eine zeitharmonische Gr oe mit der Kreisfrequenz ! = 2 f l at sich mit Hilfe der Eulerschen Formel

e ^j' = c o s ' + j sin ' 2.1

folgendermaen darstellen: n o = jF ~ r j É cos !t+ arg F ~ r :

F ~ r É e ^j!t 2.2

Wie man sieht, ist bei einer festgelegten Frequenz f die zeitharmonische Gr oe F ~ rrt durch d i e k omplexe Gr oe F ~ r eindeutig bestimmt. Bei der Verwendung von komplexen Gr oen ist @

= j! 2.3

@t

was die L osung partieller Dierentialgleichungen wie etwa d e r M a x w ellschen Gleichungen wesentlich erleichtert.

Da an den Kreuzungspunkten die Kontinuit atsbedingung 20

div ~ J = j! 2.4

erf ullt sein mu, m ussen neben den Str omen auerdem noch zeitharmonische Punktladungen an den Kreuzungspunkten angenommen werden. Diese sind f ur

9

KAPITEL 2. HERLEITUNG DER MATRIXGLEICHUNGEN 10

Abbildung 2.1: Bezeichnungen

die Berechnung der Str ome aus den magnetischen Feldern jedoch ohne Bedeutung, man ben otigt sie erst, wenn auch die elektrischen Felder ber ucksichtigt werden.

Das Gitter, mit dem die Leiterplatte modelliert wird, hat k Kreuzungspunkte in x-Richtung Index 1 : : : k und l Kreuzungspunkte in y-Richtung Index 1 : : : l . Die Kreuzungspunkte haben jeweils den Abstand a. B e i d e r B e z e i c hnung eines Kreuzungspunkts steht der Index f ur die x-Richtung vor dem f ur die y-Richtung. Die Str ome haben immer die Indizes des Knotens, von dem sie in Koordinatenrichtung ausgehen. Wie man in Abbildung 2.1 erkennt, sind k 1 Él Str ome in x-Richtung und k Él 1 Str ome in y-Richtung vorhanden.

Die magnetischen Feldst arken in x- u n d ^{y-Richtung H} x und ^H y werden uber den Kreuzungspunkten des Gitters in einer konstanten H ohe ^z 0 jeweils uber dem gemessen. Die Felder haben immer die Indizes des Kreuzungspunkts, sie gemessen wurden, zum Kreuzungspunkt g e h oren also die Werte ^{H x} und ^H y .

2.2 Maxwellsche Gleichungen und Potentiale

Die Maxwellschen Gleichungen in dierentieller Form f ur zeitharmonische Gr oen im Vakuum lauten 19, 200:

rot ^~ = ~ J + ^{j ! ~} 2.5

H E

KAPITEL 2. HERLEITUNG DER MATRIXGLEICHUNGEN 11

rot ~ E = j! ~ B 2.6

Daneben gelten noch die Materialgleichungenn das Vakuum ist linear, isotrop und homogen:

~ D = " 0 ~ E 2.9 ~ B = 0 ~ H: 2.10

Weil f ur beliebiges ~ A

= 0

rot ~ A 2.11 div

gilt, ist es zweckm aig, ~ B = r o t ~ A 2.12

zu setzen, da dadurch s i c hergestellt wird, da die Gleichung 2.7 stets erf ullt ist. Man nennt ~ A Vektorpotential, weil aus ihm durch eine Dierentiationsoperation eine Feldgr oe berechnet werden kann. Das Vektorpotential ist durch Gleichung 2.12 jedoch n i c ht eindeutig festgelegt, da wegen

= 0

rot 2.13

f ur eine beliebige skalare Funktion neben ~ A auch j e d e s ~ A + grad die Gleichung 2.12 erf ullt.

Mit Gleichung 2.12 erh alt man aus Gleichung 2.6

= 0 :

rot ~ E + j!rot ~ ~ E + j! ~ A = r o t A 2.14

Man setzt daher wegen Gleichung 2.13

~ E + j! ~ A = grad' 2.15

da dadurch die Gleichung 2.14 stets erf ullt ist. Man nennt ' das skalare Potential, es ist durch ~ E und die Wahl von ~ A eindeutig festgelegt. Kennt man die Potentiale ~ A und ', k ann man aus Gleichung 2.12 und

~ E = grad ' j! ~ A 2.16

die magnetischen und elektrischen Felder ~ B und ~ E berechnen.

Die verbleibende Freiheit bei der Wahl von ~ A benutzt man, um die Erf ullung der Gleichung div ~ A + 0 " 0 j!'= 0 2.17

zu fordern. Dies ist die sogenannte Lorentzeichung.

KAPITEL 2. HERLEITUNG DER MATRIXGLEICHUNGEN 12

Aus den bisher noch nicht v erwendeten Gleichungen 2.5 und 2.8 folgen mit diesen Festlegungen nach einigen Umformungen die Wellengleichungen f ur die Potentiale: div ~ rot

rot ~ ² ~ A = 0 ~ A A J grad 2.18

² : ^2.19 " 0

^{Dabei ist} 0 " 0 É ! = ! = 2

= 2.20 p

c 0

die Phasenkonstante.

Diese Wellengleichungen sind inhomogen, falls Str ome ieen oder Ladungen vorhanden sind. Ihre allgemeine L osung unter der Randbedingung, da die Potentiale im Unendlichen verschwinden, ist nach 20:

ZZ Z ~

A~ r = ⁰ É e dV

2.21 ^{0 0} j~ r ~ r 4 ZZ Z ~ r ⁰ j

V 0

É e '~ r = 1 j j~ r~ r ^{0 j} dV : 2.22

^{0 0} j~ r ~ r 4" 0

⁰ j

V 0

2.3 Berechnung des Vektorpotentials ~ A

Um im folgenden die Handhabung der Gleichungen zu erleichtern, werden alle auftretenden L angenmae auf den Gitterabstand a normiert:

x = x y = y z = z

~ = É a: ~ ~ ~ 2.23

x und ~ y identisch mit den Indizes und der Knoten-Dabei sind die Gr oen ~ punkte des Gitters nach Abbildung 2.1.

Wendet man die Gleichung 2.21 auf die Anordnung nach Abbildung 2.1 an, so erh alt man folgendes Ergebnis siehe auch h 2 :

Z ^k1 X X

p ^p ~ xs ² +~ y ² +~ z 2 ¹ 0 e j ^~ l

A z ~ xx ~ yy~ z = 0: 2.26

Diese Gleichungen lassen sich mit

Z p ^p ~ xs ² +~ y ² +~ z 2 ¹ yy~ z = 0 e

0

KAPITEL 2. HERLEITUNG DER MATRIXGLEICHUNGEN 13

auch s c hreiben als

^k1 X X l

A z ~ xx ~ yy~ z = 0: 2.30

2.4 Berechnung der magnetischen Feldst arke ~

H

Aus den Gleichungen 2.10 und 2.12 folgt: 0 1

^@A y ^@A z

: 2.31

^{@z @x} 0 0

^@A y ^{@A y} @x @y

F ur die Komponenten H x und H y erh alt man daraus:

X X ^{k l1} @f~ y ~ x ~ z 1

É I y 2.32

0 @z

^{=1 =1 k1} X X ^l @f~ x ~ y ~ z 1

H y ~ xx ~ yy~ z = É I x : 2.33

0 @z

=1 =1

Mit

@f~ xx ~ yy~ z @f~ xx ~ yy~ z yy~ z = 1 = 1

ergibt sich s c hlielich:

X X k l1

=1 =1

Bei der Anordnung nach Abbildung 2.1 werden diese beiden Gleichungen aus- z 0 = ^{z 0} x = 1 2 3 : : : k ~ y = 1 2 3 : : : l ~ z = ~ gewertet f ur ~ a .

KAPITEL 2. HERLEITUNG DER MATRIXGLEICHUNGEN 14

Die Funktion ^g ~ xx ~ yy ~ z besitzt die Symmetrien

^g 1 ~ xx ~ yy ~ = g ~ xx ~ yy ~ 2.37

die sich bei der Berechnung der notwendigen Funktionswerte benutzen lassen, um Rechenzeit zu sparen.

F ur die numerische Integration ben otigt man abschlieend noch den Real-und Imagin arteil von ^g ~ xx ~ yy ~

^z: Z n o = ~ 1

p 2.41

~ x s ² + ~ ² + ~ ^ds: 2

y z

Die linearen Gleichungssysteme 2.35 und 2.36 lassen sich auch mit Hilfe von Matrizen darstellen: ^H x = M 1 É ~

2.42

^{~ I y H} y = M 2 É ~

2.43

^{~ I} x :

^H y Vektoren der Dimension k Él mit den Werten der magne-Dabei sind ^{~ H} x und ~

tischen Feldst arke an den Kreuzungspunkten des Gitters nach Abbildung 2.1, ^I x ist ein Vektor der Dimension k 1 É l mit den Werten der Str ome in x-

^I y ein Vektor der Dimension k Él1 mit den Werten der Str ome Richtung und ~

in y-Richtung. M 1 und M 2 sind Matrizen, die die Abbildung der Str ome auf

die magnetischen Feldst arken beschreiben.

Mit den Gleichungen 2.42 und 2.43 kann man aus den Str omen die magnetischen Feldst arken berechnen. Bei dem hier gegebenen Problem sind jedoch die magnetischen Feldst arken bekannt und die Str ome gesucht. Man mu also die linearen Gleichungssysteme f ur die Str ome l osen. Wie dies geschieht, wird in den folgenden Kapiteln beschrieben.

Kapitel 3

Berechnung des

Matrix-Vektor-Produkts

In diesem Kapitel soll gezeigt werden, da Matrix-Vektor-Produkte, denen eine Struktur entsprechend den Gleichungen 2.35 und 2.36 zugrunde liegt, sehr zeit- und speicherplatzsparend mittels einer Schnellen Fouriertransformation FFT berechnet werden k onnen. Dies wird sich b e i d e r V erwendung von iterativen Verfahren zur L osung der Gleichungssysteme als sehr vorteilhaft erweisen. Die Idee f ur die hier vorgestellte Methode stammt a u s 1 , Abschnitt 6.22.

3.1 Der eindimensionale Fall

3.1.1 Die Struktur der Matrix

Der besseren Verst andlichkeit und der einfacheren Darstellung wegen wird zun achst nicht die dem zweidimensionalen Fall der Gleichungen 2.35 und 2.36 entsprechende Beziehung X X n N

y kkK = t kllKL É x llL k = 1 2 3 : : : m K = 1 2 3 : : : M 3.1 l=1 L=1

Aquivalent

X ⁿ y k = t kl É x l k = 1 2 3 : : : m : 3.2 l=1

Da in diesem Kapitel auer den Indizes s amliche Gr oen komplex sind, wird zugunsten einer einfacheren Darstellung hier auf den kennzeichnenden Unterstrich verzichtet.

Gleichung 3.2 ist mit der eindimensionalen diskreten Faltung ⁺¹ X

15

KAPITEL 3. BERECHNUNG DES MATRIX-VEKTOR-PRODUKTS 16

verwandt, die durch Diskretisierung der eindimensionalen kontinuierlichen Faltung

⁺¹ Z

1

entsteht. x n , y n und h n sind dabei Folgen, die durch die Diskretisierung aus den kontinuierlichen Funktionen xt, yt u n d ht h e r v orgehen.

ur die numerische Auswertung stets notwendige Beschr ankung auf endliche Teilfolgen x l l = 1 2 3 : : : nund y k k = 1 2 3 : : : mgeht

in t n u m benannt. Eine Gleichungsstruktur wie in 3.2 entsteht grunds atzlich bei der regelm aigen Diskretisierung von eindimensionalen Faltungen das heit Integralgleichungen mit Verschiebungskern nach G l e i c hung 3.4.

Schreibt man die Folgen x l und y k a l s V ektoren

~ x = x 1 x 2 x 3 : : : x n ^T 3.5

3.6

so l at sich die Gleichung 3.2 auch als Matrix-Vektor-Produkt

~ y = T É ~ x 3.7

darstellen mit 0 1 t 0 t 1 t 2 t 3 É É É t 1n

t m1 t m2 t m3 t m4 É É É t mn

Man erkennt leicht, da die Elemente der Matrix T auf einer Diagonalen jeweils gleich sind. Eine Matrix mit dieser Struktur heit Toeplitz-Matrix.

3.1.2 Die zyklische diskrete Faltung

Die zyklische diskrete Faltung wird auf periodische Folgen ^ f n a n g e w andt. Eine periodische Folge der Periodenl ange p ist durch

^ f n+iÉp = ^ f n i = : : : 2 1 0 1 2 : : : 3.9

deniert. Bei der zyklischen diskreten Faltung wird im Gegensatz zu Gleiuber eine Periodenl ange summiert:

^p1 X

x n ^Ë ^ h n = x m É ^ ^ y n = ^ h nm : ^ 3.10 m=0

KAPITEL 3. BERECHNUNG DES MATRIX-VEKTOR-PRODUKTS 17

Die Diskrete Fouriertransformation DFT periodischer Folgen ist deniert durch

^p1 X

^k=0 Neben ^ ^f n ist auch ^ ^F k periodisch:

^ ^F k+iÉp = ^ i = : : : 2 1 0 1 2 : : : : 3.13

^F k

Es l at sich zeigen 23, da man damit das Ergebnis der zyklischen diskreten Faltung nach G l e i c hung 3.10 auch folgendermaen erhalten kann: ^p1 X ^x n Ë ^ ^h n = 1

^ ^X k ^ j2 ^nk ^ ^y n = ^ 3.14

^{p : H} k É e

^p k=0

Die zyklische diskrete Faltung ergibt sich also durch die elementweise Multiplikation der transformierten Folgen und die anschlieende R ucktransformation des Ergebnisses. Wird f ur die Durchf uhrung der DFT ein geeigneter Algorithmus der Schnellen Fouriertransformation FFT verwendet, so l at sich d u r c h dieses Verfahren eine erhebliche Anzahl an Rechenoperationen einsparen.

3.1.3 Berechnung des Matrix-Vektor-Produkts

Die Berechnung des Matrix-Vektor-Produkts 3.7 mit Hilfe der zyklischen diskreten Faltung geschieht d a d u r c h, da man zun achst die periodischen Folgen ^ ^x k u n d ^ ^t k bildet, von denen eine Periode jeweils wie folgt lautet:

pn

D a b e i s i n d d i e x i die Komponenten des Vektors ~ x nach Gleichung 3.5 und die t i die Elemente der Matrix T nach Gleichung 3.8. Anschlieend berechnet x k Ë ^ man nach Gleichung 3.14 ^ y k = ^ t k . Wie man durch Einsetzen leicht

nachpr ufen kann, sind dann die ersten m Elemente der Periode von ^ y k die Elemente des Matrix-Vektor-Produkts ~ y nach G l e i c hung 3.6:

pm

^zz|| : : : 3.17

Beim Einsetzen stellt man ferner fest, da statt der Nullen in der Folge ^ t k a u c h beliebige andere Werte stehen k onnen, da sie die ersten m Ele- mente von ^ y k nicht beeinussen. Weiterhin ist auch die Periodenl ange p ohne

KAPITEL 3. BERECHNUNG DES MATRIX-VEKTOR-PRODUKTS 18

Bedeutung f ur die ersten m Elemente von ^ y k , solange p m + n 1 ist. Es uberhaupt keine weiteren Werte in der Mitte von ^ t k z u

stehen, t m1 kann direkt an t 1n anschlieen. Man sollte p so w ahlen, da sich ein f ur den FFT-Algorithmus g unstiger Wert ergibt. Selbstverst andlich m u p f ur ^ x k u n d ^ t k gleich s e i n .

3.2 Der zweidimensionale Fall

3.2.1 Die Struktur der Matrix

Die Behandlung des zweidimensionalen Falls nach Gleichung 3.1 erfolgt in vollkommener Analogie zum eindimensionalen Fall, wie er im vorhergehenden Kapitel dargestellt wurde.

Gleichung 3.1 ist mit der zweidimensionalen diskreten Faltung ⁺¹ X ⁺¹ X

y nnN = x nnN Ë h nnN = x mmM É h nmmNM 3.18 m=1 M=1

verwandt, die durch Diskretisierung der zweidimensionalen kontinuierlichen Faltung

⁺¹ Z ⁺¹ Z

x 1 2 Éht 1 1 t 2 2 dd 1 dd 2 3.19 yt 1 t 2 = xt 1 t 2 Ëht 1 t 2 = 1 1

entsteht. x nnN , y nnN u n d h nnN sind dabei Folgen, die durch die Diskre-

ur die numerische Auswertung stets notwendige Beschr ankung auf endliche Teilfolgen x llL l = 1 2 3 : : : n L = 1 2 3 : : : Nund y kkK k = 1 2 3 : : : m K = 1 2 3 : : : Mgeht G l e i c hung 3.18 schlielich uber dabei wurde die Folge h nnN i n t nnN u m benannt. Eine Gleichungsstruktur wie in 3.1 entsteht grunds atzlich bei der regelm aigen Diskretisierung von zweidimensionalen Faltungen das heit Integralgleichungen mit Verschiebungskern nach Gleichung 3.19.

Schreibt man die Folgen x llL und y kkK a l s V ektoren

~ x = x 11 x 21 : : : x nn1 x 12 x 22 : : : x nn2 : : : x 1N x 2N : : : x nnN ^T 3.20

~ y = y 11 y 21 : : : x mm1 y 12 y 22 : : : y mm2 : : : y 1M y 2M : : : y mmM ^T 3.21

so l at sich die Gleichung 3.1 auch als Matrix-Vektor-Produkt

~ y = T É ~ x 3.22

darstellen mit

KAPITEL 3. BERECHNUNG DES MATRIX-VEKTOR-PRODUKTS 19

0 1 ^t 10 É É É ^t 1nn0 ^t 11 É É É ^t 1nn1 É É É ^t 01N ^{t 11N É É Ét} 1nn1N ^t 00 ^t 01 B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B @ T

= 3.23

Man erkennt, da die einzelnen Bl ocke der Matrix T jeweils eine Toeplitz-Struktur haben, und da auerdem auch die Bl ocke selbst in einer Toeplitz-Struktur angeordnet sind. Eine Matrix mit dieser Struktur heit deshalb Block-Toeplitz-Toeplitz-Block-Matrix.

3.2.2 Die zweidimensionale zyklische diskrete Faltung

Die zweidimensionale zyklische diskrete Faltung wird auf zweidimensionale periodische Folgen ^ ^f nnN angewandt. Eine zweidimensionale periodische Folge mit den Periodenl angen p und P ist durch

^ ^f n+iÉppN+IÉP = ^ ii I = : : : 2 1 0 1 2 : : : 3.24

^f nnN

deniert. Bei der zweidimensionalen zyklischen diskreten Faltung wird im Geuber jeweils eine Periodenl ange summiert:

^P1 X ^p1 X

^x nnN Ë ^ ^h nnN = ^x mmM É ^ ^h nmmNM : 3.25

^ ^y nnN = ^ ^ m=0 M=0

Die zweidimensionale Diskrete Fouriertransformation 2D-DFT zweidimensionaler periodischer Folgen ist deniert durch

^P1 X ^p1 X

^f nnN É e ^p É e j 2 ^nk j 2 NK P ^ ^ ^F kkK = 3.26 ^{n=0 N=0 P1} X ^p1 X

^f nnN = 1 ^F kkK Ée j2 ^{nk p} É e j2 ^{NK P} : ^ ^ 3.27 p ÉP K=0 k=0

KAPITEL 3. BERECHNUNG DES MATRIX-VEKTOR-PRODUKTS 20

Neben ^ ^f nnN ist auch ^ ^F kkK periodisch:

^ ^F k+iÉppK+I ÉP = ^ ii I = : : : 2 1 0 1 2 : : : : 3.28

^F kkK

Analog zum eindimensionalen Fall kann man das Ergebnis der zweidimensionalen zyklischen diskreten Faltung auch mit Hilfe der zweidimensionalen diskreten Fouriertransformierten erhalten:

X X ^x nnN Ë ^ ^h nnN = 1 j 2 ^{nk p} É e p1 j 2 ^{NK P} : 3.29 P ¹ ^ ^ ^ ^y nnN = ^

^X kkK ^H kkK É e

K =0 p É P

k=0

Die zyklische diskrete Faltung ergibt sich also auch h i e r d u r c h die elementweise Multiplikation der transformierten Folgen und die anschlieende R ucktransformation des Ergebnisses. Wird f ur die Durchf uhrung der DFT ein geeigneter Algorithmus der zweidimensionalen Schnellen Fouriertransformation 2D-FFT verwendet, so l at sich durch dieses Verfahren eine noch g r oere Anzahl an Rechenoperationen als beim entsprechenden eindimensionalen Fall einsparen.

3.2.3 Berechnung des Matrix-Vektor-Produkts

Die Berechnung des Matrix-Vektor-Produkts 3.22 mit Hilfe der zweidimensionalen zyklischen diskreten Faltung geschieht dadurch, da man zun achst die zweidimensionalen periodischen Folgen ^ ^x k u n d ^ ^t k bildet, von denen eine Periode jeweils wie folgt lautet:

P ^N z |

^x 11 ^x 12 ^x 13 ^{É É É x} 1N 0 É É É 0 ^x 21 ^x 22 ^x 23 ^{É É É x} 2N 0 É É É 0 ^x 31 ^x 32 ^x 33 ^{É É É x} 3N 0 É É É 0

. . . . . . . . . . . . . . . . . . . . . . . .

^ ^x kkK = 3.30

^x nn1 ^x nn2 ^x nn3 ^{É É É x} nnN 0 É É É 0

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

tm10 tm11 tm12 É É É ^tm1M 1 0 É É É 0 tm11N tm12N É É É tm12 tm11

9

0 0 0 É É É 0 0 É É É 0 0 0 0 0 = É É É

^ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ^t kkK = pmn+1

.

0 0 0 É É É 0 0 É É É 0 0 0 0 0 É É É

t1nn0 t1nn1 t1nn2 É É É ^t1nnM 1 0 É É É 0 t1nn1N t1nn2N É É É t1nn2 t1nn1 t2nn0 t2nn1 t2nn2 É É É ^t2nnM 1 ⁰ É É É 0 t2nn1N t2nn2N É É É t2nn2 t2nn1

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

t20 t21 t22 ^{É É Ét2M} 1 0 É É É 0 ^{t21N t22N} É É Ét22 t21 t10 t11 t12 ^{É É Ét1M} 1 0 É É É 0 t11N t12N É É Ét12 t11

^{| z} 3.31

PMN+1

KAPITEL 3. BERECHNUNG DES MATRIX-VEKTOR-PRODUKTS 21

Dabei sind die x iiI die Komponenten des Vektors ~ x nach G l e i c hung 3.20 und die t iiI die Elemente der Matrix T nach G l e i c hung 3.23. Anschlieend berech- x kkK Ë^ net man nach Gleichung 3.29 ^ y kkK = ^ t kkK . Wie man durch Einsetzen leicht nachpr ufen kann, sind dann die ersten m É M Elemente der Periode von ^ y kkK die Elemente des Matrix-Vektor-Produkts nach G l e i c hung 3.21:

^PM zz||

Beim Einsetzen stellt man ferner fest, da statt der Nullen in der Folge ^ auch beliebige andere Werte stehen k onnen, da sie die ersten m É M Elemente von ^ y kkK n i c ht beeinussen. Weiterhin sind auch d i e P eriodenl angen p und P ohne Bedeutung f ur die ersten mÉM Elemente von ^ y kkK , solange p m+n1 uberhaupt keine weiteren Werte in der Mitte von ^ t kkK zu stehen, die t m1K k onnen direkt an die t 1nnK und die t kkM1 direkt an die t kk1N anschlieen. Man sollte p und P so w ahlen, da ur den 2D-FFT-Algorithmus g unstiger Wert ergibt. Selbstverst andlich m ussen p und P f ur ^ x kkK und ^ t kkK gleich sein.

3.3 Matrix-Vektor-Multiplikation in Bl ocken

Verwendet man man zur Berechnung der auf der Leiterplatte ieenden Str ome neben magnetischen auch elektrische Feldst arken, so erh alt man eine Matrix, die in ihrer Gesamtheit keine Block-Toeplitz-Toeplitz-Block-Struktur mehr aufweist 22. Man stellt jedoch fest, da die Matrix aus mehreren Bl ocken aufgebaut ist, die jeder f ur sich eine Block-Toeplitz-Toeplitz-Block-Struktur besitzen. Auch das Produkt einer solchen Matrix mit einem Vektor l at sich sehr zeit-und speicherplatzsparend mit Hilfe der beschriebenen Methode berechnen, wenn man die Multiplikation in Bl ocken durchf uhrt 17, A b s c hnitt 222. Die Matrix T sei aus M ÉN Bl ocken T kl mit jeweils m k Zeilen und n l Spalten

P M k=1 m k Zeilen und P N l=1 n l Spalten. aufgebaut, die gesamte Matrix besitze also

Dann unterteilt man den Vektor ~ x in N Teilvektoren ~ x l mit jeweils n l Elemen- P N n l Elementen. Anschlieend beten, der gesamte Vektor besteht also aus

rechnet man analog der bekannten Regel f ur die Matrix-Vektor-Multiplikation 0 1 1 0 1 0 P N l=1 T 1l É ~ x l

T 11 T 12 T 13 É É É T 1N ~ x 1

B B B B B B @ ~ y = T É ~ x =

KAPITEL 3. BERECHNUNG DES MATRIX-VEKTOR-PRODUKTS 22

F ur jedes einzelne Teilprodukt kann dabei der im vorhergehenden Abschnitt beschriebene Algorithmus verwendet werden, da es sich hierbei jeweils um ein Produkt einer Block-Toeplitz-Toeplitz-Block-Matrix mit einem Vektor handelt.

Kapitel 4

Fredholmsche

Integralgleichungen

erster Art

Die Gleichungen 2.21 und 2.22, die dem hier zu l osenden Problem zugrunde liegen, bezeichnet man auch als Fredholmsche Integralgleichungen erster Art. Gleichungen dieses Typs besitzen eine Reihe von Eigenschaften, die eine numerische Auswertung gelegentlich s c hwierig gestalten. Eine recht umfassende und leicht v erst andliche Einf uhrung in diese Problematik ist in 14 z u n d e n .

Die in diesem Kapitel gegebene Darstellung folgt weitestgehend 6. Dort sind auch zahlreiche weitere Literaturhinweise enthalten, die hier nicht n o c hmals angef uhrt werden sollen.

4.1 Einf uhrung

Fredholmsche Integralgleichungen erster Art spielen in vielen Problemen der Natur- und Ingenieurwissenschaften eine wichtige Rolle. Einige Beispiele sind der Entwurf von Antennen, Astrometrie, Computer-Tomographie, Bildr uckgewinnung, inverse Geo- und Helioseismologie und mathematische Biologie. Viele weitere Beispiele existieren. Eine Fredholmsche Integralgleichung erster Art hat im eindimensionalen Fall die Form

^b Z

^K ss tf t dt = gs c s dd 4.1

a

wobei die Funktionen ^K der Integralkern und g die rechte Seite zumindest prinzipiell bekannt sind, w ahrend ^f die unbekannte, gesuchte Funktion ist. In vielen, aber nicht allen praktischen Anwendungsf allen von 4.1 ist der Integralkern ^K durch das zugrunde liegende mathematische Modell genau gegeben, w ahrend die rechte Seite ^g typischerweise aus Mewerten besteht, das heit g ist

23

KAPITEL 4. FREDHOLMSCHE INTEGRALGLEICHUNGEN 1. ART 24

nur mit einer begrenzten Genauigkeit und nur an einer endlichen Anzahl von Punkten ^s 1 ^s 2 ^{: : : s} m bekannt.

Fredholmsche Integralgleichungen erster Art sind grunds atzlich s c hlecht g estellte Probleme, das heit die L osung ist auerst empndlich gegen uber beliebig kleinen St orungen des Systems. Das heit, da alle klassischen numerischen Verfahren wie LU- und Cholesky-Faktorisierung keine aussagekr aftige L osung berechnen k onnen, nachdem 4.1 diskretisiert wurde. Andererseits bedeutet das nicht, da eine L osung nicht berechnet werden k onnte. Die Entwicklung stabiler und verl alicher Verfahren, die f ur die L osung von 4.1 besonders geeignet sind, ist daher stets eine Herausforderung gewesen. Fr uher stellte die ziemlich langsame Rechengeschwindigkeit vieler Computer eine starke Begrenzung f ur die rechnerische Komplexit at und Verfeinerung der praktischen numerischen Verfahren dar, die f ur die L osung von 4.1 vorhanden waren. Moderne Computer und Workstations haben jedoch eine viel schnellere Rechenleistung und erlauben es dadurch dem heutigen Benutzer, wesentlich fortgeschrittenere numerische Verfahren zu verwenden, um die Integralgleichungen zu untersuchen u n d z u l osen.

uber solche fortgeschrittenen

numerischen Verfahren gegeben werden. Insbesondere soll gezeigt werden, wie die Singul arwertzerlegung und die verallgemeinerte Singul arwertzerlegung verwendet werden, um wichtige Einsichten in das schlecht gestellte Problem zu erhalten und dadurch helfen, eine aussagekr aftige L osung f ur die Integralgleichung zu berechnen. Ein anderes wichtiges Hilfsmittel, das vorgestellt werden soll, ist die sogenannte L-Kurve, die den Graph der Norm oder Seminorm der uber der Norm des Residuums darstellt. Die L-Kruve ist ein einfaches Mittel, um den Einu der Regularisierung zu zeigen und hilft dem Benutzer, einen guten Regularisierungsparameter zu w ahlen.

Es soll nochmals darauf hingewiesen werden, da in diesem Kapitel die numerischen Verfahren f ur die Behandlung schlecht gestellter Probleme vorgestellt werden sollen, wobei angenommen wird, da das Problem bereits diskretisiert wurde, so da man eine Matrixgleichung A~ x = ~ b vor sich h a t o d e r d e n T erm kA~ x ~ b k minimieren mu. Es soll nicht auf die Vielzahl der Diskretisierungsverfahren eingegangen werden, die in der Literatur ausf uhrlich behandelt werden. Oft ist zumindest eine der Dimensionen der Matrix A durch das gegebene Pro- festgelegt, zum Beispiel, wenn die rechte Seite ^{~ b} aus einer festen Anzahl von Messungen aus einem Versuch besteht, der nicht l e i c ht wiederholt werden kann. Der Hauptzweck d e s n umerischen Verfahrens ist es dann, aus den Daten uber das gegebene Problem zu erhalten. so viel Informationen wie m oglich

In Abschnitt 4.2 werden zun achst einige theoretische Ergebnisse zu Fredholmschen Integralgleichungen erster Art vorgestellt, weil diese die beste M oglichkeit bieten, die Schwierigkeiten schlecht gestellter Probleme zu verstehen. Ab Abschnitt 4.3 werden dann die numersichen Methoden zur L osung solcher Probleme vorgestellt. In Abschnitt 4.3 wird gezeigt, wie die Singul arwertzerlegung verwendet wird, um die Matrix zu analysieren, die durch die Diskretisierung von Gleichung 4.1 entsteht. Abschnitt 4.4 fat die verschiedenen Re- gularisierungsverfahren zusammen, die verwendet werden, um Gleichung 4.1

KAPITEL 4. FREDHOLMSCHE INTEGRALGLEICHUNGEN 1. ART 25

numerisch z u l osen, und in Abschnitt 4.5 wird gezeigt, wie die verallgemeinerte Singul arwertzerlegung ein Mittel daf ur bietet, diese regularisierten Probleme zu ur regularisierte Prouber der Norm des Residuums, uber das Problem

liefert und verwendet werden kann, um den Regularisierungsparameter geeignet zu w ahlen.

4.2 Allgemeine Schwierigkeiten schlecht gestellter

Probleme

Wie bereits erw ahnt wurde, kann die Integralgleichung 4.1 auerst schwierig zu l osen sein, weil sie ein schlecht gestelltes Problem darstellt. Solche Probleme sind dadurch g e k ennzeichnet, da beliebig kleine St orungen der rechten Seite ^g zu beliebig groen St orungen der L osung f f uhren k onnen. Die L osung ist also auerst empndlich gegen uber St orungen.

Diese Schwierigke i t e n s i n d u n trennbar verbunden mit der Kompaktheit des Operators, der mit dem Integralkern ^K in Verbindung steht. Physikalisch gedeuubt die Integration mit ^K in Gleichung 4.1 einen gl attenden Eekt auf ^f aus in dem Sinne, da hochfrequente Anteile wie etwa Spitzen und Spr unge durch d i e I n tegration herausgegl attet werden. Es ist daher zu erwarten, da der umgekehrte Vorgang, das heit die Berechnung von ^f aus g, dazu neigt, hochfrequente Anteile in ^g zu verst arken. Wie weiter unten dargestellt wird, ist dies tats achlich d e r F all.

4.2.1 Die Singul arwertentwicklung

Das am besten geeignete Hilfsmittel f ur die Analyse Fredholmscher Integralgleichungen erster Art ist die Singul arwertentwicklung singular value expansion, SVE von ^K. Mit Hilfe der SVE kann jeder quadratisch i n tegrierbare Integralkern ^K geschrieben werden als die unendliche Summe

X

1 ^K ss t = ^v i t 4.2

i=1

bei einem degenerierten Integralkern ist 1 durch den Rang des Kerns zu er- Die Funktionen ^u i und ^v i heien linke und rechte singul are Funktionen von ^K . Sie sind orthonormal in Bezug auf das ubliche Skalarprodukt, das heit

^u i ^u j = ^v i ^v j = ij , w obei É É deniert ist als

^Z 4.3

= t dt:

Die Gr oen i sind die Singul arwerte von ^K , sie sind stets gr oer oder gleich Null und k onnen stets in monoton fallender Folge angeordnet werden, so da

1 2 3 : : : 0:

KAPITEL 4. FREDHOLMSCHE INTEGRALGLEICHUNGEN 1. ART 26

Die Tripel ^f i ^u i ^v i g stehen mit den beiden mit K verkn upften Eigenwertproblemen in folgender Verbindung: f ² i ^u i g sind die Eigenl osungen Eigenwerte ^{R b} und Eigenfunktionen des symmetrischen Kerns a Kss xKtt x dx, w ahrend ^R d c Kxx sKxx t dx darstellen. Dies macht d e u t - i ^v i g die Eigenl osungen von f 2

lich, da die Tripel ^f i ^u i ^v i g charakteristisch und im wesentlichen einzigartig f ur einen gegebenen Integralkern ^K sind.

Vielleicht d i e w i c htigste Beziehung zwischen den Singul arwerten und den singul aren Funktionen ist die Beziehung

b

a

die zeigt, da jede singul are Funktion ^v i auf die entsprechende Funktion ^u i abgebildet wird, wobei der Singul arwert i der Vergr oerungsfaktor dieser speziellen Abbildung ist.

Setzt man diese Beziehung zusammen mit der SVE 4.2 in die Integralgleichung 4.1 ein, so erh alt man die Gleichung

1 1 X X

i ^{f v} i ^u i s = ^{g u} i ^u i s 4.5

i=1 i=1

ur die L osung von 4.1 f uhrt:

1 X

^ft = ^v i t: 4.6 i

i=1

Es wird betont, da ^f nur existiert, falls die rechte Seite von 4.6 tats achlich konvergiert. Dies ist auf jeden Fall dann gegeben, wenn die Gleichung 4.1 f ur die gegebene Funktion ^g tats achlich e i n e L osung f besitzt. Auch i n a nderen F allen kann 4.6 konvergieren, sofern die Picard-Bedingung nach A bschnitt 4.2.3 erf ullt ist. Es l at sich zeigen, da die durch 4.6 gegebene Funktion ^f in diesem Fall eine L osung im quadratischen Mittel f ur die Gleichung 4.1 darstellt. Man erkennt, da in Gleichung 4.6 die Funktion ^f mittels der sin- aren Funktionen v i und der entsprechenden Entwicklungskoezienten

i

dargestellt wird. Die L osung ^f l at sich daher vollst andig durch eine Analyse

der Koezienten i und der Funktionen v i charakterisieren.

4.2.2 Die gl attenden Eigenschaften von K

Das Verhalten der Singul arwerte i und der singul aren Funktionen u i und v i ist keineswegs " beliebigihrVerhalten ist stark verbunden mit den Eigenschaften des Integralkerns ^K. Es gilt folgendes:

Je " glatter der Kern ^K ist, desto schneller fallen die Singul arwerte ⁱ gegen Null wobei die " Glattheit d u r c h die Anzahl der stetigen partiellen Ableitungen gemessen wird.

KAPITEL 4. FREDHOLMSCHE INTEGRALGLEICHUNGEN 1. ART 27

Je kleiner die i sind, desto mehr Oszillationen oder Nulldurchg ange treten in den singul aren Funktionen ^u i und ^v i auf. Dies ist eine Folge der Tatsache, da die Fourierentwicklung von ^K solche Oszillationseigenschaften aufweisen mu.

Die praktische Bedeutung dieser Eigenschaften der Tripel f i ^u i ^v i g ist, da

ur ^f als eine Spektralentwicklung betrachtet werden kann,

in der die Koezienten i die spektralen Eigenschaften der L osung f beschreiben. Man sieht aus Gleichung 4.5, da die Integration mit ^K tats achlich eine gl attende Wirkung aus ubt: je h oher die spektralen Anteile in ^f sind, desto mehr werden sie in ^g durch die Multiplikation mit i ged ampft. Dar uber hinaus zeigt Gleichung 4.6, da das inverse Problem, n amlich das der Berechnung von ^f aus g, t a t s achlich die " umgekehrte Wirkung auf die Oszillationen in ^g aus ubt, n amlich e i n e V erst arkung der spektralen Anteile ^{g u} i mit dem Faktor 1 i . D i e s v erst arkt selbstverst andlich die hochfrequenten Anteile.

4.2.3 Die Picard-Bedingung

Wenn man das gerade angesprochene Verhalten bedenkt, ist es oensichtlich, da wegen der Verst arkungsfaktoren 1 i nicht j e d e r e c hte Seite ^g zu einer

" glatten, quadratisch i n tegrierbaren L osung ^f f uhren wird. Tats achlich m u die rechte Seite ^g etwas " glatter als die gew unschte L osung f sein, damit die rechte Seite in Gleichung 4.6 tats achlich g e g e n ^f konvergiert. Damit eine quadratisch i n tegrierbare L osung ^f f ur die Integralgleichung 4.1 existiert, mu die rechte Seite die sogenannte Picard-Bedingung

X

1 ² 1 4.7 i

i=1

erf ullen. Die Picard-Bedingung besagt, da von einem gewissen Punkt der Summation in 4.6 an die Betr age der Koezienten ^{g u} i schneller fallen m ussen als die entsprechenden Singul arwerte i , damit eine quadratisch i n tegrierbare L osung ^f existiert. Da diese Bedingung in Verbindung mit Fredholmschen Integralgleichungen erster Art so wesentlich ist, sollte sie, wo i m e r m oglich, zun achst gepr uft werden, bevor man versucht, die Integralgleichung zu l osen. Die numerischen Aspekte einer solchen Untersuchung werden im n achsten Abschnitt beschrieben.

Der Abfall der Singul arwerte i ist so grundlegend f ur das Verhalten von schlecht gestellten Problemen, da es sinnvoll ist, diesen Abfall zu benutzen, um den Grad der Schlechtgestelltheit des Problems zu kennzeichnen. Dies wurde zuerst von Wabha erw ahnt. Hofmann hat die folgende Denition vorgeschlagen: falls eine positive reelle Zahl existiert, so da die Singul arwerte die Bedingung ullen, dann heit der Grad der Schlechtgestelltheit, und das

Problem wird als leicht b z w . m aig schlecht gestellt bezeichnet, falls 1 ur alle 0 d. h. = 1, dann

wird das Problem stark schlecht gestellt genannt.

KAPITEL 4. FREDHOLMSCHE INTEGRALGLEICHUNGEN 1. ART 28

Man beachte, da f ur Probleme mit endlichem Rang einschlielich jeder Diskretisierung von 4.1 von einem rein mathematischen Gesichtspunkt aus die Picard-Bedingung stets erf ullt ist, die L osung ist stabil, und es wird keine Regularisierung ben otigt. Diskrete Probleme leiden jedoch stets unter einer Kombination von Mefehlern, Diskretisierungsfehlern und Rundungsfehlern, und die uber diesen

Fehlern. Daher wird von einem praktischen Gesichtspunkt aus Regularisierung weiterhin ben otigt, um den Einu dieser Fehler herauszultern, damit eine n utzliche L osung f ur das diskrete System berechnet werden kann. Diese Aspekte werden in Abschnitt 4.5 n aher behandelt.

4.3 Analyse der Koezientenmatrix

4.3.1 Die Singul arwertzerlegung

Das beste Hilfsmittel f ur die Analyse der Koezientenmatrix A, die aus der

Diskretisierung des Integralkerns K hervorgegangen ist, ist die Singul arwertzerlegung singular value decomposition, SVD, die das diskrete Analogon zur SVE darstellt. Die SVD einer Matrix A mit m Zeilen und n Spalten hat die

^minmn X Form

A = i ~ u i ~ v 4.8 ^{Ë i=1} ubliche

Skalarprokukt sind, das heit ~ u i ~ u j = ~ u ^T u j =~ v i ~ v j = ~ v ^T v j = ij und die Singul arwerte i von A die Bedingung

1 2 : : : minmmn 0

erf ullen. Analog zur SVE sind die Paare f ² i ~ u i g und f ² i ~ v i g die Eigenl osungen Eigenwerte und Eigenvektoren der positiv semideniten Matrizen AA ^Ë bzw. A ^Ë A. Die grundlegende Denition der SVD, das heit das Analogon zu

Gleichung 4.4 f ur die SVE, hat die Form

A~ v i = i ~ u i i = 1 2 : : : minmm n 4.9

was zeigt, da jeder Vektor ~ v i auf den entsprechenden Vektor ~ u i abgebildet wird mit i als Vergr oerungsfaktor.

4.3.2 Beziehungen zwischen der SVD und der SVE

Die grundlegenden numerischen Schwierigkeiten, die mit der L osung eines linearen Gleichungssystems A~ x = ~ b oder der Minimierung von kA ~ x ~ b k also der

L osung des Gleichungssystems im quadratischen Mittel verbunden sind, wenn dieses Gleichungssystems durch Diskretisierung aus Gleichung 4.1 entstanden ist, r uhren daher, da die berechnete L osung ~ x sehr empndlich gegen uber St orungen von A oder ~ b ist. Dies spiegelt sich i n d e r T atsache wieder, da die

KAPITEL 4. FREDHOLMSCHE INTEGRALGLEICHUNGEN 1. ART 29

Konditionszahl von A, die durch d a s V erh altnis ¹ n gegeben ist, sehr gro ist.

Dar uber hinaus w achst die Konditionszahl von A sowohl mit der Ordnung n

als auch mit der Anzahl der Datenwerte m. Diese numerischen Schwierigkeiten sind direkt mit der Schlechtgestelltheit des Problems 4.1 verbunden: je besser die Diskretisierung ist, das heit je besser das lineare Gleichungssystem die Integralgleichung wiedergibt, desto mehr ahnelt die Schlechtgestelltheit des Gleichungssystems der Schlechtgestelltheit der Gleichung 4.1.

Der Grund f ur die Schlechtgestelltheit des linearen Gleichungssystems liegt in der Tatsache, da die SVD der Matrix A sehr eng mit der SVE des Integralkerns K verkn upft ist. Tats achlich sind die Singul arwerte i von A in

vielen F allen N aherungen der Singul arwerte i des Integralkerns K, w ahrend uber die singul aren Funktionen von K geben. Insbesondere kann f ur ein Galerkin-Entwicklungsverfahren mit orthonormalen Basisfunktionen i und i gezeigt werden, da

0 i i k K ~ 4.10

wobei die Norm von K ~ K n gegeben ist durch

v

und wobei ~ K n ein degenerierter Kern vom Rang n ist, der durch

n n

~ ^{X X} K n ss t = i t 4.11

i=1 j=1

gegeben ist, wobei die a ij = i t ds dt die Elemente der Ma- ^R b ^R d

^a trix A sind. Entsprechend gilt, wenn u ij und v ij die Elemente der singul aren

Vektoren ~ u j und ~ v j bezeichnen, und wenn ~ u j und ~ v j N aherungen der singul aren Funktionen u i und v i sind, die durch

n n

^{X X} u j s = u ij i s ~ v ij i s v j s = ~ 4.12

^{i=1 i=1} u j und ~ v j durch gegeben sind, da die N aherungsfehler in ~

v t 2kK ~

beschr ankt sind. Als Folge aus dieser engen Beziehung zwischen der SVD und der SVE ist es nun oensichtlich, da f ur jede Diskretisierung eines schlecht gestellten Problems gilt:

1. Die Matrix A ist schlecht k onditioniert, das heit die Konditionszahl ¹ n

ist gro.

KAPITEL 4. FREDHOLMSCHE INTEGRALGLEICHUNGEN 1. ART 30

2. Die Konditionszahl von A steigt mit n.

Dar uber hinaus haben die singul aren Vektoren ^{~ u} j und ^{~ v} j eine steigende Anzahl von Vorzeichen anderungen, wenn j w achst, und es kann gezeigt werden, da ur groe j zutrit, wo die Schranke in Gleichung 4.13 weniger eng dies auch f wird.

4.3.3 SVD-Analyse

Die Schlufolgerung aus der Betrachtung im letzten Abschnitt ist, da klassische Methoden f ur die L osung von A~ x = ~ b oder die Minimierung von kA ~ x ~ b k,

wie zum Beispiel Cholesky-, LU- oder QR-Faktorisierung, nicht v erwendet werden k onnen, um die numerische L osung einer Fredholmschen Integralgleichung erster Art nach Gleichung 4.1 zu berechnen. Einerseits ist die Konditionszahl von A stets so gro, da Rundungsfehler die Berechnung einer exakten nume- L osung ^{~ x} verhindern. Andererseits w urde man wegen der Oszillationen der singul aren Vektoren selbst dann keine " glatte L osung ^{~ x} erhalten, wenn man in der Lage w are, das lineare Gleichungssystem ohne Rundungsfehler zu l osen, weil das diskretisierte Problem stets durch N aherungs- und Diskretisierungsfehler gest ort ist, die Anteile entlang aller singul aren Vektoren haben. Es gibt mehrere Gr unde, die SVD einer Matrix A zu berechnen. Wab- schlug dies im Jahr 1980 in ihrem einureichen Bericht v or, in dem die Verf ugbarkeit numerischer Programme betont w i r d . W enn die SVD einmal berechnet ist, gestattet sie dem Anwender, mittels der Singul arwerte und der singul aren Vektoren die " spektralen Eigenschaften des Operators ^K in den Begrien der SVE zu studieren, wie sie in Abschnitt 4.2 vorgestellt wurde.

uberpr ufen, ob die Picard-Man kann die SVD auch v erwenden, um zu

Bedingung f ur die Integralgleichung nach Gleichung 4.7 f ur das zugrunde liegende, ungest orte Problem erf ullt zu sein scheint, indem man die Schau-

i betrachtet.

die Fehler in der rechten Seite ^{~ b} in der Regel Anteile entlang allen singul aren Vektoren. Vorausgesetzt, da die Picard-Bedingung erf ullt ist, werden die Fehler daher nur auf die Koezienten ^{~ b ~ u} i einen starken Einu haben, die kleinen Singul arwerten entsprechen. Falls daher f ur kleine Indizes i die Koezienten ^{~ b ~ u} i im Mittel st arker abfallen als die Singul arwerte i , ist es in der Tat sinnvoll anzunehmen, da die Picard-Bedingung erf ullt ist. Falls dar uber hinaus die Koezienten ^{~ b ~ u} i f ur steigendes i schlielich ein Plateau erreichen, dann ist dieses Plateau eine Absch atzung f ur die Gr oe des Fehlers der rechten Seite. In diesem Zusammenhang ist es wichtig, sich klarzumachen, da die N aherungsfehler in den berechneten singul aren Funktionen f ur kleine i am geringsten sind, vergleiche Gleichung 4.13.

KAPITEL 4. FREDHOLMSCHE INTEGRALGLEICHUNGEN 1. ART 31

4.4 Regularisierung und Filterung

4.4.1 Regularisierung

Ein vern unftiger Weg, um eine aussagekr aftige, " glatte L osung der Integralgleichung 4.1 zu berechnen, das heit eine L osung, die einige n utzliche Eigenschaften gemeinsam hat mit der exakten L osung des zugrundeliegenden und unbekannten ungest orten Problems, ist es, auf irgendeine Weise die hochfrequenten Anteile herauszultern, die mit den kleinen Singul arwerten in Verbindung stehen. Dieser Ansatz ist selbstverst andlich n ur dann hilfreich, wenn die verlangte L osung tats achlich eine gewisse Glattheit besitzt. Obwohl das nicht immer der Fall ist, zum Beispiel bei der Bildr uckgewinnung, wo m a n c hmal Kanten und Unstetigkeiten in der L osung gewollt sind, gibt es eine groe Klasse von Problemen, bei denen es vern unftig ist, eine " glatte,n aherungsweise L osung zu suchen, bei der die hochfrequenten Anteile herausgeltert sind.

Die klassische Vorgehensweise, um die hochfrequenten Anteile, die mit den kleinen Singul arwerten in Verbindung stehen, herauszultern, ist es, das Problem einer Regularisierung zu unterziehen. Der Begri Regularisierung war urspr unglich v erkn upft mit einer bestimmten von Tikhonov v orgeschlagenen Techublichen Terminologie wird jedoch jedes Verfahren, das nik, nach der heutzutage

zur Berechnung einer " glatten L osung dient, als Regularisierungsverfahren bezeichnet. Bei dem urspr unglichen Verfahren wurde die Regularisierung direkt auf die Integralgleichung 4.1 angewandt. Eine Regularisierung kann jedoch g enauso gut angewandt werden auf das lineare Gleichungssystem, das durch d i e Diskretisierung aus 4.1 entstanden ist. Dies ist in der Praxis oft wesentlich einfacher durchzuf uhren und wegen der engen Verbindung zwischen dem urspr unglichen Problem und dem diskretisierten Problem wie im vorhergehenden Abschnitt erw ahnt ist die Auswirkung der Regularisierung auf die berechnete " glatte L osung im wesentlichen dieselbe. In dieser Darstellung beschr ankt sich die Diskussion daher auf die Regularisierung der linearen Probleme A~ x = ^{~ b} und min kA~ ^{x ~} b k.

Die grundlegende Schwierigkeit bei diesen Systemen ist die Instabilit at, in dem Sinne, da die L osung ^{~ x} durch Beitr age, die den kleinen Singul arwerten entsprechen, beherrscht w i r d , w obei diese Beitr age weitgehend aus Fehlern bestehen. Die geringe Menge an Information in der rechten Seite, die mit den kleinen Singul arwerten in Verbindung steht, geht d u r c h d i e F ehler verloren. Daher kann man sagen, da die Systeme im Grunde unterbestimmt sind, weil man nur in der Lage ist, die Information wiederzugewinnen, die mit den groen Singul arwerten von A verkn upft ist. Um eine eindeutige L osung berechnen zu

k onnen, mu man die Stabilit at daher erzwingen, indem man zus atzliche Information bereitstellt, die

1. genau eine L osung aussondert und

2. eine L osung auszusondern versucht, die in einem gewissen Sinne nahe bei der gew unschten, aber unbekannten exakten L osung liegt.

KAPITEL 4. FREDHOLMSCHE INTEGRALGLEICHUNGEN 1. ART 32

Kennt man von vorneherein einen Sch atzwert ^{~ x} 0 der L osung, so ist es nat urlich, wenn man die Seminorm der Dierenz zwischen der berechneten L osung ^{~ x} und dem Sch atzwert ^{~ x} 0 zu minimieren versucht:

min kL~ x ~ ^x 0 k 4.14

wobei die Matrix L eine geeignet gew ahlte Matrix ist. Typischerweise ist L entweder die Identit atsmatrix I n I n bezeichnet die Identit atsmatrix mit n Zeilen

und n Spalten oder eine diskrete N aherung f ur einen Dierentiationsoperauber die gew unschte L osung vorhanden tor. Falls keine besonderen Kenntnisse

sind, ist es vern unftig, ^{~ x} 0 = 0 zu setzen, w ahrend es schwieriger ist, eine spezielle Matrix L zu empfehlen. Versuche mit mehreren verschiedenen L sind

m oglicherweise notwendig.

4.4.2 Das Tikhonov-Verfahren

Es gibt viele M oglichkeiten, die Randbedingung 4.14 mit dem urspr unglichen linearen Problem zu verkn upfen. Vielleicht der bekannteste Ansatz ist derjenige, der unabh angig von Tikhonov und Phillips vorgeschlagen wurde, n amlich, die regularisierte L osung ^{~ x} als L osung des Problems n kA~ ^{x ~} b k o

kL~ x ~ ^x 0 k min + 4.15 2 2 2

zu denieren. Hier ist die Gr oe der Regularisierungsparameter, der das relative Gewicht z w i s c hen der Minimierung der Randbedingung und der Minimierung der Norm des Residuums bestimmt. Wenn der Schnitt der Nullr aume von A und L trivial ist, ist die L osung ~ x eindeutig und formal durch

A + L ¹ A L~ x 0 ^{~ x} = A L L ^Ë ~ b +

Ë Ë Ë 2 2

gegeben. Diese Formel sollte jedoch n i c ht v erwendet werden, um ~ x tats achlich zu berechnen. Einerseits ist mit der Bildung der Matrix A ^Ë A wegen der endli- Rechengenauigkeit ein unvermeidlicher Verlust an Information verbunden. Wichtiger noch ist jedoch, da der Aufwand, der zur Berechnung der GSVD des Matrixpaars A,L n o t wendig ist, nicht w esentlich g r oer ist als der, der mit der Bildung der Matrix A ^Ë A + ² L L und der L osung des entsprechenden

Ë

Gleichungssystems verbunden ist. Daher stellt bei ann ahernd gleichem Rechenuber das zu l osende Problem aufwand die GSVD wesentlich mehr Informationen zur Verf ugung.

Man beachte, da die Tikhonov-Regularisierung in Gleichung 4.15 sowohl auf quadratische m = n als auch a u f uberbestimmte m m n Gleichungs- systeme angewandt werden kann. Man beachte weiterhin, da der grunds atzliche Gedanke bei der Tikhonov-Regularisierung ist, ein Verh altnis zwischen der Gr oe der Norm des Residuums und der Randbedingung einzuf uhren. Wenn man ein groes Gewicht auf die Randbedingung legt, bedeutet dies, da man ein gr oeres Residuum hinnehmen mu und umgekehrt. Indem man einen passenden Regularisierungsparameter w ahlt, kann man eine befriedigende L osung aussondern, f ur die die beiden Bedingungen ausgeglichen sind. Die geeignete Auswahl von wird in Abschnitt 4.6 weiter behandelt.

KAPITEL 4. FREDHOLMSCHE INTEGRALGLEICHUNGEN 1. ART 33

4.4.3 Abgeschnittene Singul arwertzerlegung

Ein anderes wohlbekanntes Regularisierungsverfahren ist es, einfach die auf der SVE basierende Entwicklung nach G l e i c hung 4.6 abzuschneiden, bevor die kleinen Singul arwerte zu dominieren beginnen. Bei den Gleichungssystemen wird die entsprechende Technik abgeschnittene SVD genannt, und die mit ihr verbundene regularisierte L osung ^{~ x} k ist deniert durch X k

^{~ x} k = 4.16 ^{~ v} i :

i i=1

Hier spielt die nat urliche Zahl k die Rolle des Regularisierungsparameters. Es kann gezeigt werden, da die abgeschnittene SVD im wesentlichen mit der Tikhonov-Regularisierung mit L = I n gleichwertig ist in dem Sinne, da f ur

jedes k ein existiert, so da k~ ^{x ~ x} k k = ^O k ist, was klein ist, falls

die Picard-Bedingung erf ullt ist. In Abschnitt 4.5 wird auf den allgemeinen Fall eingegangen werden, wo L 6 = I n ist.

4.4.4 Iterative V erfahren

Es gibt auch eine ganze Klasse von Regularisierungsverfahren f ur die L osung der linearen Gleichungssysteme A~ x = ~ b und min kA~ ^{x ~} b k, die mit iterativen

Verfahren in Verbindung stehen. Die vielversprechendsten hiervon sind semiiterative V erfahren und das Verfahren der konjugierten Gradienten 15, 166. Bei der hier vorgestellten Version des Verfahrens der konjugierten Gradienten, die auch als CGNR bezeichnet wird, wird der Algorithmus auf die Normalengleichung A ^Ë A~ x = A uberbestimmten Glei- ^Ë ~ angewandt, deren L osung bei

chungssystemen die L osung von A~ x = ~ b im quadratischen Mittel ist. Dabei wird jedoch das Matrixprodukt A ^Ë A nicht t a t s achlich berechnet, sondern die Terme A ^{Ë ~ b} A ^Ë A~ x werden in A ^{Ë ~ b} A~ x umgeformt, wobei die Subtraktion

vor der abschlieenden Matrix-Vektor-Multiplikation durchgef uhrt wird. Diese Variante wird allgemein f ur die numerisch stabilste gehalten. Der Algorithmus f ur L = I n f ur dieses Verfahren lautet nach h 1 , Algorithmus 6.33 und 16, Algo-

Man w ahle zun achst einen Sch atzwert ^{~ x} 0 f ur die L osung ~ ^x 0 = 0 , f a l l s k eine uber die L osung vorhanden sind und berechne

^d 0 = ~ b

Danach iteriere man f ur k = 1 2 3 : : : , bis die Abbruchbedingung erf ullt ist:

k = k~ ^r k1 k 2

KAPITEL 4. FREDHOLMSCHE INTEGRALGLEICHUNGEN 1. ART 34

^d k1 k A~ ^s k ^d k = ~ ~

Bei diesem Algorithmus ist ^{~ d} k der Defekt nach k Iterationen, das heit ^{~ d} k =

^{~ b} A~ ^x k . Der Vektor ^{~ s} k heit die k. Suchrichtung. Wie man sieht, ist dieses

Verfahren besonders dann g unstig, falls sich Matrix-Vektor-Produkte mit den Matrizen A und A ^Ë mit wenig Aufwand berechnen lassen, und falls bis zum

Abbruch des Verfahrens nicht allzu viele Schritte notwendig sind.

Es l at sich zeigen, da ^{~ x} k eine regularisierte L osung ist. Existiert zum Beispiel eine L ucke zwischen j ^{~ b ~ u} k j und j ^{~ b ~ u} k+1 j, s o l i e g t ^{~ x} kCGNR nahe bei der L osung der abgeschnittenen SVD ^{~ x} kTSVD wegen k~ ^x kCGNR ^{~ x} kTSVD k =

^{j ~ b~ u} k+1 ^j . Die abgeschnittene SVD ist andererseits gleichwertig mit der

^{O j ~ b~ u} k j

Tikhonov-Regularisierung vgl. Abschnitt 4.4.3, und der Abbruch des Verfahrens der konjugierten Gradienten ist daher aquivalent zu einer Tikhonov-Regularisierung mit L = I n . E s g i b t a u c h M oglichkeiten, das Verfahren der konjugierten Gradienten auf Probleme anzuwenden, bei denen L 6 = I n ist.

Um eine schnellere Konvergenz des Verfahrens der konjugierten Gradienten zu erreichen, was zu einem fr uheren Abbruch d e s V erfahrens und damit zu eiublich, eine Pr akonditionierung

der Matrix durchzuf uhren. Wendet man das Verfahren der konjugierten Gradienten zur Regularisierung an, so mu man dabei allerdings beachten, da die Unterr aume der Singal- und der St oranteile weiterhin getrennt bleiben, da sonst die Ergebnisse der Iteration von Anfang an St oranteile enthalten. Das hierzu notwendige Verfahren, bei dem auch die Vorteile der ezienten Berechnung von Matrix-Vektor-Produkten mit A und A ^Ë erhalten bleiben, ist in 2 a u s f uhrlich

beschrieben. Was geschieht, wenn man eine herk ommliche Pr akonditionierung durchf uhrt, ohne auf die Trennung der beiden Unterr aume zu achten, ist in 11, Abbildung 7 sehr deutlich z u e r k ennen.

4.4.5 Geeignete Diskretisierung

Neben den bisher dargestellten Verfahren l at sich eine Regularisierung auch dadurch erreichen, da man das kontinuierliche Problem auf eine geeignete Weise diskretisiert. Die in diesem Abschnitt gemachten Ausf uhrungen folgen dabei 1, Abschnitt 3.33.

Um eine numerische N aherung f ur ^f nach G l e i c hung 4.1 zu erhalten, mu man das urspr ungliche Problem auf irgendeine Weise diskretisieren. Es gibt mehrere M oglichkeiten, eine Integralgleichung auf eine endlich-dimensionale Matrixgleichung zu reduzieren, die bekanntesten sind dabei Quadraturverfahren, die Momentenmethode 18 und das Galerkin-Verfahren. Es sollen hier nicht die Vor- und Nachteile der einzelnen Verfahren dargestellt werden, wegen der

KAPITEL 4. FREDHOLMSCHE INTEGRALGLEICHUNGEN 1. ART 35

Schlechtgestelltheit des kontinuierlichen Problems scheint e s e h e r s o , a l s o b j edes Problem, m oglicherweise sogar jede rechte Seite eine andere Wahl f ur die bestm ogliche Diskretisierung verlangt. Stattdessen sollen einige grunds atzliche Gesichtspunkte von Diskretisierungen von Gleichung 4.1 betrachtet werden.

Es ist ein wichtiger Gesichtspunkt von Projektionsverfahren, wie dem Galerkin-Verfahren, da sie von sich aus regularisierende Eigenschaften aufweisen, vorausgesetzt, die Unterr aume werden passend gew ahlt. Dies wurde zuerst von Natterer beobachtet und wurde anschlieend von zahlreichen Autoren n aher analysiert. Die Schlufolgerung daraus ist folgende: Falls die Diskretisierung zu grob ist, ist die endlich-dimensionale Matrixgleichung ziemlich gut konditioniert, aber ihre L osung enth alt einen ziemlich groen Diskretisierungsfehler, dies

verwendet. Falls andererseits die Disretisierung zu fein ist, dann ist der Einu der kleinen Singul arwerte von K zu gro, und die L osung des diskretisierten Problems ahnelt einer Tikhonov-Regularisierung, bei der zu klein gew ahlt wurde. Irgendwo zwischen diesen beiden Extremen bendet sich ein optimaler Grad der Diskretisierung, bei dem die sich ergebende N aherung vergleichbar ist mit der L osung, die man aus einer optimalen Tikhonov-Regularisierung des kontinuierlichen Problems 4.1 erh alt.

Leider ist diese Theorie f ur praktische Zwecke s c hwierig anzuwenden. Zum einen ist die optimale Diskretisierung fast nie im voraus bekannt. Dar uber hinaus w are es selbst dann, wenn man den optimalen Grad der Diskretisierung zuf allig treen w urde, eine gute Idee f ur die numerische Stabilit at, das endlichdimensionale Problem zu regularisieren, so da die Konditionszahl verbessert wird. In der Praxis w ahlt man deshalb besser eine zu feine Diskretisierung und regularisiert dann die Matrixgleichung, als ob sie ein schlecht gestelltes Problem darstellen w urde. Es soll an dieser Stelle nochmals betont w erden, da das diskrete Problem von einem streng mathematischen Gesichtspunkt aus gesehen gut gestellt ist, da jede nichtsingul are Matrix automatisch eine stetige Inverse hat. Wegen der groen Konditionszahl der Matrix beobachtet man jedoch ahnliche Erscheinungen wie bei dem kontinuierlichen schlecht gestellten Problem, wenn man das endlich-dimensionale Gleichungssystem l osen m ochte.

4.5 Analyse regularisierter Probleme

4.5.1 Die verallgemeinerte Singul arwertzerlegung

Analog zur SVE und zur SVD, die die n utzlichsten Hilfsmittel zur Untersuchung der Integralgleichung und des diskretisierten Systems sind, gibt es auch ein sehr gut geeignetes Hilfsmittel f ur die Analyse von diskreten Regularisierungsproblemen mit allgemeinen Matrizen L. Dieses Hilfsmittel ist die verallgemeinerte

Singul arwertzerlegung generalized singular value decomposition, GSVD des Matrixpaars A,L, die urspr unglich v on Van Loan eingef uhrt wurde. F ur unseren Zweck ist es ausreichend, den Fall zu betrachten, in dem L vollen Rang hat und A m Zeilen und n Spalten sowie L p Zeilen und n Spalten hat, wobei

KAPITEL 4. FREDHOLMSCHE INTEGRALGLEICHUNGEN 1. ART 36

m n p gilt. Eine allgemeinere Formulierung f ur alle m, n und p macht n ur die Notation um einiges schwieriger, ohne wirklich mehr Licht auf das Problem zu werfen. Man deniert drei quadratische Matrizen

U ~ ~ u 1 ~ ~ u 2 : : : ~ ~ ~

W ~ w 1 ~ w 2 : : : ~

U = I m U ~ wobei die Vektoren ~ ~ u i und ~ ~ v i orthonormal sind das heit ~ ^Ë V ~ = I p , w ahrend die Vektoren ^~ V und ~ ^{Ë w} i linear unabh angig sind das heit W ist nichtsingul ar. Dann hat die GSVD von A,L d i e F orm

0 1 diag i 0

B @ C A W ¹ L = ~ A = ~ I np U

Hier sind diag i und diag i Diagonalmatrizen mit den nicht negativen Ele- i =1 menten 1 2 : : : p bzw. 1 2 : : : p , die jeweils die Gleichung ² i + ² f ur i = 1 2 : : : p erf ullen. Die verallgemeinerten Singul arwerte werden als die Verh altnisse i = ⁱ i deniert, und sie werden in monoton fallender Folge angeordnet, das heit 1 2 : : : p 0:

Falls insbesondere L = I n ist, dann ist ^~ ~ u i = ~ u i , ^~ ~ v i = ~ v i und i = i f ur

i = 1 2 : : : n , und Gleichung 4.17 wird zur normalen SVD von A. E s g i b t

andere verwandte Formulierungen der GSVDD die hier verwendete ist f ur den hier verfolgten Zweck geeignet, n amlich die oben erw ahnten diskreten Regularisierungsverfahren zu untersuchen.

4.5.2 Wichtige Beziehungen der GSVD

Ausgestattet mit der GSVD des Matrixpaars A,L ist es folgerichtig, die regularisierten L osungen auf dieselbe Weise zu untersuchen, wie man es mit Hilfe der SVD f ur die nicht regularisierten L osungen getan hat. Um dies f ur allgemeine Sch atzwerte ^{~ x} 0 also insbesondere ^{~ x} 0 6 = 0 zu zeigen, ist es hilfreich, den Vektor mit n Elementen

^{~ !} ! 1 ^! 2 : : : ! n ^T = W ^{1 ~ x} 0

zu denieren. Dann kann gezeigt werden, da f ur jedes lineare Regularisierungsverfahren eine Menge von Filterfaktoren f i i = 1 2 : : : p existiert, so da die regularisierte L osung ^{~ x} reg geschrieben werden kann als !

X X p n

+ 1 f i ^{! i ~ x} reg = w i + ^{~ b ~} ~ 4.18 f i ~ u i ~ w i :

i

i=1 i=p+1

Verwendet man diese Gleichung zusammen mit 4.17, so kann man leicht z e i - gen, da die Seminorm von L~ ^x reg ~ ^x 0 , das heit der Randbedingung, gegeben

KAPITEL 4. FREDHOLMSCHE INTEGRALGLEICHUNGEN 1. ART 37

ist durch

^{v 2} : ~ b ~ ~ p u i u

^{u X} kL ~ ^x reg ~ ^x 0 k = f ² i ! i 4.19 ^t i i

i=1

In gleicher Weise ist die Norm des Residuums gegeben durch

^v p m ² 1 f i ² j ~ b ~ ~ kA~ x reg ~ b k = ~ b ~ ~ ^{X X} u i i ! i j ² + u i : 4.20

^{u t} i=1 i=n+1

Die Filterfaktoren f i sind Funktionen des Regularisierungsparameters, der mit dem jeweiligen Regularisierungsverfahren in Verbindung steht. F ur viele ur die Filterfak-Regularisierungsverfahren gibt es ziemlich einfache Ausdr ucke f toren. Zum Beispiel ist bei der Tikhonov-Regularisierung

f i = ²

ⁱ 4.21

²

w ahrend f ur das Verfahren der konjugierten Gradienten, das nach q Iterationsschritten gestoppt wird, die Filterfaktoren

i P q ² f i = ² 4.22

sind, wobei P q die sogenannten Ritz-Polynome q. Ordnung bezeichnet, die mit dem Verfahren der konjugierten Gradienten in Verbindung stehen. Betrachtet man die abgeschnittene SVD, so ist es oensichtlich, da die entsprechenden Filterfaktoren einfach 0 u n d 1 s i n d . D i e V erallgemeinerung der abgeschnittenen SVD auf allgemeine Regularisierungsmatrizen L 6 = I n ist die abgeschnittene

ur andere Regularisie-GSVD, wiederum mit den Filterfaktoren 0 und 1. Auch f rungsverfahren k onnen die Filterfaktoren in dieser Weise angegeben werden.

Man beachte, da es, falls p p n ist, n p Komponenten der regularisierten L osung ~ x reg gibt, die nicht durch die Filterfaktoren f i beeinut werden und daher vom Regularisierungsparameter unabh angig sind. Diese Komponenten s i n d i n R i c htung der Vektoren ~ w p+1 ~ w p+2 : : : ~ w n gerichtet, die Nullvektoren der Matrix L sind, das heit

L ~ w i = 0 i = p + 1 p + 2 : : : n :

F ur jede sinnvolle Wahl der Regularisierungsmatrix L sind diese Nullvektoren

" glatt u n d b e n otigen daher keine Regularisierung.

4.5.3 Untersuchungen mit der GSVD

Man sieht n un, da man analog zu der fr uher erw ahnten SVD die GSVD von A,L v erwenden kann, um Einsicht in die Regularisierungseigenschaften der jeweils verwendeten Regularisierungsverfahren zu erhalten. Beispielsweise w i der Matrix W und der entsprechenden wird die Untersuchung der Spalten ~ Filterfaktoren f i Auskunft uber die spektralen Eigenschaften des Regularisie-

rungsverfahrens geben. Typischerweise stellt man fest, da groe Filterfaktoren

KAPITEL 4. FREDHOLMSCHE INTEGRALGLEICHUNGEN 1. ART 38

mit " glatten V ektoren ^{~ w} i verkn upft sind, w ahrend kleine Filterfaktoren, die eine D ampfung in die regularisierte L osung bringen, mit Vektoren ^{~ w} i in Verbindung stehen, die einen hohen Anteil an hochfrequenten Komponenten besitzen. Wenn dies der Fall ist, berechnet man tats achlich eine " glatte regularisierte L osung.

Das Konzept der diskreten Picard-Bedingung wurde von Varah vorgeschlagen und in 9 w eiter ausgearbeitet. Diese Bedingung besagt, da die Koezienten j ^{~ b ~} ~ ^u i j im Durchschnitt schneller fallen m ussen als die verallgemeinerten Singul arwerte i , damit sichergestellt ist,

1. da man eine " glatte regularisierte L osung berechnen kann,

2. da die regularisierte L osung eine vern unftige N aherung der exakten L osung des ungest orten Problems ist.

Die Untersuchung des Verhaltens der verallgemeinerten Singul arwerte i und der Koezienten j ^{~ b ~} ~ ^u i j ergibt daher eine tiefere Einsicht in das gestellte Prouberpr ufen, ob f ur kleine Indizes i, w o d i e F ehler blem. Insbesondere kann man in ^{~ b} das Skalarprodukt ^{~ b ~} ~ ^u i nicht beherrschen, die Koezienten j ^{~ b ~} ~ ^u i j

schneller fallen als die verallgemeinerten Singul arwerte i . F alls dies tats achlich der Fall ist, sagt man, da die diskrete Picard-Bedingung erf ullt ist.

Die Untersuchung der Koezienten der GSVD gibt auch eine wertvolle Hilfe bei der Auswahl eines guten Regularisierungsparameters. Es wird daran erinnert, da, weil die Filterfaktoren nur die Beitr age zur regularisierten L osung d ampfen sollen, die mit Koezienten in Verbindung stehen, die von Fehlern beherrscht w erden, man den Regularisierungsparameter so w ahlen mu, da die Filterfaktoren nur die gew unschten Koezienten d ampfen. Die Untersuchung ^u i und ^{~ b ~} ~

der durch die GSVD erhaltenen Gr oen i , ^{~ b ~} ~ i in Verbindung

mit den Filterfaktoren ^f i gibt daher einen guten Hinweis zur optimalen Wahl des Regularisierungsparameters. FORTRAN-Unterprogramme zur Berechnung der SVD und der GSVD sind in LAPACK 30, 311 siehe Abschnitt 5.2.2 vorhanden.

4.6 Die L-Kurve

Eines der wichtigsten Probleme in Verbindung mit der numerischen Behandlung linearer Gleichungssysteme, die von schlecht gestellten Problemen herr uhren, ist die Wahl des Regularisierungsparameters. Wie in der Einleitung erw ahnt, ist es das Ziel, so viel Information wie m oglich aus der gegebenen rechten Seite ^{~ b} herauszuholen. Eine zu starke Regularisierung unterdr uckt Information, die tats achlich i n ^{~ b} vorhanden w are, w ahrend eine zu geringe Regularisierung zu einer L osung f uhrt, die von Fehlern beherrscht wird. Daher sollte man idealerweise den Wert des Regularisierungsparameters nden, der ein ausgeglichenes Verh altnis herstellt zwischen dem Regularisierungsfehler das heit dem Fehler,

KAPITEL 4. FREDHOLMSCHE INTEGRALGLEICHUNGEN 1. ART 39

der durch d i e G l attung der Daten hervorgerufen wird und dem St orungsfehler, der von den Fehlern in ^{~ b} stammt. Man nennt diesen Wert des Regularisierungsparameters optimal wobei zu beachten ist, da diese Denition des optimalen Regularisierungsparameters abweicht v on der anderer Autoren, bei denen er f ur eine optimale Konvergenzrate sorgt, wenn der Fehler in ^{~ b} gegen Null geht. Der optimale Regularisierungsparameter im hier verwendeten Sinne P n

ist eng verkn upft mit dem " eektiven Rang,derals i=1 ^f i deniert ist und

angibt, wieviel verl aliche Information dem Gleichungssystem entnommen werden kannn beispielsweise ist f ur die abgeschnittene SVD der " eektive Rang einfach die Anzahl der Koezienten in der Entwicklung nach G l e i c hung 4.16, die ungleich Null sind. Ein n utzliches Hilfsmittel f ur die Untersuchung dieser Gesichtspunkte ist die sogenannte L-Kurve 77.

4.6.1 Die Denition der L-Kurve

Die L-Kurve ist der Graph der Randbedingung kL ~ ^x reg ~ ^x 0 k uber der Norm des

Residuums kA~ ^x reg ~ b k f ur ein bestimmtes Regularisierungsverfahren. Daher ist die L-Kurve i n d e r T at eine parametrisierte Kurve, deren Parameter der Regularisierungsparameter ist, beispielsweise im Falle der Tikhonov-Regularisierung oder k im Falle der abgeschnittenen SVD und des Verfahrens der konjugierten Gradienten. Der Name " L-Kurve stammt v on der Tatsache, da f ur viele Probleme die Kurve kA~ ^x reg ~ b k kL ~ ^x reg ^{~ x} 0 k einen L-f ormigen " Knick hat. Die L-Kurve ist zum einen Teil deshalb ein so n utzliches Hilfsmittel bei schlecht gestellten Problemen, weil sie eine gute M oglichkeit darstellt, die gegenseitige Abh angigkeit zwischen der Seminorm kL ~ ^x reg ~ ^x 0 k und der Norm des Residuums kA~ ^x reg ~ b k darzustellen, und zum anderen Teil deshalb, weil der L-f ormige " Knick oft einem Wert des Regularisierungsparameters entspricht, der ann ahernd optimal im obigen Sinne ist.

uck auf Miller sowie Lawson und Han-Die Verwendung der L-Kurve geht z u r

sonn sie war jedoch als ein Hilfsmittel zur Untersuchung von Regularisierungsproblemen nie weit verbreitet. Dennoch ist die L-Kurve ein praktisches Mituber das Regularisierungsproblem in einer kompakten

Form darzustellen. Zun achst illustriert die L-Kurve unmittelbar den notwendigen Kompromi zwischen der Minimierung der Randbedingung kL ~ ^x reg ^{~ x} 0 k

uber hinaus

zeigt das Auftreten eines " Knicks in der Kurve, da tats achlich ein spezieller Wert des Regularisierungsparameters existiert, der in einem gewissen Sinne optimal ist, weil er die beiden Minimierungsziele gegeneinander ausgleicht. Diese Information l at sich mit der L-Kurve l e i c ht darstellen, sie ist mit anderen Schaubildern schwieriger darzustellen.

4.6.2 Der " Knick in der L-Kurve

uberhaupt existiert und Um zu verstehen, warum ein " Knick in der L-Kurve warum er mit einem ann ahernd optimalen Wert des Regularisierungsparameters

KAPITEL 4. FREDHOLMSCHE INTEGRALGLEICHUNGEN 1. ART 40

. .

Abbildung 4.1: Das grunds atzliche Aussehen der L-Kurve f ur die Tikhonov-Regularisierung in einem doppelt logarithmischen Mastab. Die gestrichelte ur ein ungest ortes Problem, w ahrend die gepunkte-Linie zeigt die L-Kurve f

ur eine rechte Seite zeigt, die ausschlielich a u s F ehlern te Linie die L-Kurve f besteht.

verkn upft ist, ist es vern unftig, das Erscheinungsbild der L-Kurve detaillierter zu untersuchen. Es wird hier die Analyse aus 7 zusammengefat, weitere Einzelheiten sind in dem Artikel selbst zu nden. Geeigneterweise betrachtet man das folgende Modell, bei dem die Matrizen A und L exakt gegeben sind und sich

^~ ^ b plus einer zuf alligen St orung ~ e die rechte Seite ^{~ b} aus einem ungest orten Vektor bestehend aus Mefehlern, N aherungsfehlern usw. zusammensetzt. Falls man annimmt, da

1. die diskrete Picard-Bedingung erf ullt ist und

^~ ^ bk erf ullt, 2. die Norm von ^{~ e} die Bedingung k~ e k k

dann wird die L-Kurve einen L-f ormigen " Knick wie in Abbildung 4.1 aufweisen. Es gibt zwei unterschiedliche Teile der L-Kurve, einen achen Teil rechts vom " Knick und einen steilen Teil oberhalb des " Knicks.

Der ache Teil der L-Kurve rechts vom " Knick e n tspricht den regularisierten L osungen, bei denen der Regularisierungsfehler dominiert, das heit, wo die D ampfung so gro ist, da sie erfolgreich den meisten Einu der Fehler ^{~ e} d ampft und bis zu einem gewissen Ausma auch Information aus der zugrunde liegenden ungest orten rechten Seite ^{~ b} herausltert. Je weiter rechts man sich auf der Kurve bendet, desto mehr D ampfung ist vorhanden, und die Kurve biegt sich deshalb schielich n a c h u n ten, bis die Randbedingung vollst andig erf ullt ist: kL ~ ^x reg ~ ^x 0 k = 0. Die erste Annahme oben stellt sicher, da eine angemessen " glatte L osung des ungest orten Problems existiert, das heit, da die Seminorm kL ~ ^x reg ~ ^x 0 k klein bleibt, so da der Teil der L-Kurve, der mit dem ungest orten Problem in Verbindung steht die gestrichelte Kurve i n Abbildung 4.1 f ur ! 1 einen achen Teil besitzt.

KAPITEL 4. FREDHOLMSCHE INTEGRALGLEICHUNGEN 1. ART 41

Der steile Teil der L-Kurve oberhalb des " Knicks e n tspricht den regularisierten L osungen, bei denen die Beitr age der Fehler ^{~ e} dominieren, weil zu wenig D ampfung vorhanden ist. Aus Gleichung 4.19 erkennt man, da, wenn wenig Filterung vorhanden ist, die Seminorm kL ~ ^x reg ~ ^x 0 k wegen der Division durch kleine verallgemeinerte Singul arwerte i schnell ansteigt. Wendet man noch w eniger Filterung an, erreicht die Kurve s c hlielich ein Plateau an der Stelle, an der im wesentlichen der ganze Einu der Fehler extrahiert wurde dieses Plateau w urde bei unendlich-dimensionalen Problemen nicht auftreten. Die zweite Annahme oben stellt sicher, da der Teil der L-Kurve, der ausschlielich m i t den Fehlern in Verbindung steht die gepunktete Kurve in Abbildung 4.1 sich links von der Stelle nach u n ten biegt, an der der mit dem ungest orten Problem verkn upfte Teil dies tut, und die Abszisse bei k~ ^e k schneidet.

^~ ^ b + ~ e ist, entsteht ur das reale Problem, bei welchem ^{~ b} = ^~ ^ ^b und ~ e in selbstverst andlich aus der Kombination der beiden Teile, die mit Verbindung stehen. Da f ur die meisten Werte des Regularisierungsparameters einer der beiden Beitr age zu ^{~ x} reg vorherrschen wird, wird die L-Kurve als die durchgezogene Linie in Abbildung 4.1 erscheinen, mit einem steilen Teil dort, wo d i e S t orungsfehler in ^{~ x} reg vorherrschen, und einem achen Teil dort, wo d i e Regularisierungsfehler in ^{~ x} reg dominieren. Dar uber hinaus gibt es einen kleinen

sind. Es l at sich zeigen, da der " Knick in einem doppelt logarithmischen Mastab besonders ausgepr agt ist.

4.6.3 Die Wahl des Regularisierungsparameters

Aus dieser kurzen Untersuchung der L-Kurve ist es oensichtlich, da die optimale Wahl f ur den Wert des Regularisierungsparameters einem Punkt nahe beim " Knick der L-Kurve e n tspricht, weil dieser Punkt f ur eine L osung mit einem g unstigen Verh altnis zwischen den beiden Fehlerarten steht. Ist der Regularisierungsparameter kleiner als der optimale Wert, so wird die L osung zu sehr von Beitr agen der Fehler ^{~ e} beeinut, w ahrend andererseits ein Informationsverlust auftritt, falls man den Regularisierungsparameter wesentlich g r oer als den optimalen Wert w ahlt die Begrie " gr oer u n d " kleiner beziehen sich dabei auf den Parameter der Tikhonov-Regularisierung, bei der abgeschnittenen SVD und dem Verfahren der konjugierten Gradienten entspricht

und umgekehrt.

F ur die richtige Auswahl des Regularisierungsparameters ist es dabei wesentlich, da man einen doppelt logarithmischen Mastab w ahlt, um die L-Kurve darzustellen 1, Abschnitt 5.44. Der Grund daf ur ist, da bei dieser Darstellungs-Anderung von kL ~ ^x reg ~ ^x 0 k geteilt Anderung von kA~ ^x reg ~ b k darstellt, und dieses Verh altnis h angt grob gesagt davon ab, welche spektralen Anteile von ^{~ b} = ^~ ^ b + ~ e in ^{~ x} reg

vorherrschen. Im linearen Mastab dagegen stellt die Steigung das Verh altnis Anderungen in kL ~ ^x reg ^{~ x} 0 k und kA~ ^x reg ~ b k dar,

KAPITEL 4. FREDHOLMSCHE INTEGRALGLEICHUNGEN 1. ART 42

das praktisch unabh angig von ~ b ist, so da der lineare Mastab nicht v erwendet werden kann, um den Nutzanteil ~ ^ b vom Fehler ~ e zu trennen.

Es ist interessant festzustellen, da andere Verfahren f ur die Auswahl eines optimales Wertes des Regularisierungsparameters, die auf vollst andig anderen Kriterien als der L-Kurve basieren, oft zu einem Regularisierungsparameter f uhren, der nahe bei dem liegt, der aufgrund der L-Kurve g e w ahlt wurde. Dies ist beispielsweise der Fall f ur das Diskrepanzprinzip, das Quasi-Optimalit ats-Kriterium und die verallgemeinerte Kreuz-Validierung.

ur

die Norm des Fehlers angibt, e k ~ e k, und dann den Regularisierungsparameter so w ahlt, da die Norm des Residuums die Bedingung

kA~ x reg ~ b k = e

erf ullt. In Begrien der L-Kurve bedeutet dies einfach, den Schnittpunkt einer vertikalen Linie bei e mit der L-Kurve z u n d e n . F alls die obere Schranke ^e eine gute Sch atzung der Norm des Fehlers ist, wird dieser Schnittpunkt nahe beim " Knick der L-Kurve liegen, jedoch ein wenig rechts von ihm mit der sich daraus ergebenden Gefahr, da die L osung zu sehr gegl attet wird indem man die Koezienten zu sehr d ampft.

Eine gemeinsame Eigenschaft des Quasi-Optimalit ats-Kriteriums und der verallgemeinerten Kreuz-Validierung ist es, da sie keine zus atliche Information ben otigen, weil sie implizit die Norm k~ e k des Fehlers aus den Daten zu sch atzen versuchen. Das bei der Kreuz-Validierung zugrunde liegende Prinzip ist, da, falls eine beliebige Beobachtung weggelassen und dann aufgrund der verbleibenden m 1 B e o b a c htungen vorhergesagt wird, der optimale Regularisierungsparameter die Summe der Quadrate dieser Vorhersagefehler minimiert. Die verallgemeinerte Kreuz-Validierung basiert auf demselben Prinzip und stellt dar uber hinaus sicher, da der gefundene Regularisierungsparameter einige w unschenswerte Invarianzeigenschaften hat, wie beispielsweise Invarianz gegen uber einer orthogonalen Transformation der Daten was Permutationen einschliet. Bei der Tikhonov-Regularisierung f uhrt dies dazu, da der Wert des Regularisierungsparameters gew ahlt wird, der die Funktion

kA~ x ~ bk ² 4.23 ^Ë A + ² L ^Ë L ¹ A jSpurI m AA ^Ë j 2

minimiert. Die Funktion G hat idealerweise dort ein Minimum, wo i n d e r L osung ~ x die Vorherrschaft der Regularisierungsfehler in die Vorherrschaft ubergeht. Grunds atzlich das gleiche trit f ur die Quasi-

Optimalit ats-Funktion d

d die Ableitung von ~ x nach ist. Aus den vorhergehenden Ausf uhrunwobei

uber die L-Kurve und insbesondere ihren " Knick ist es oensichtlich, da sowohl die verallgemeinerte Kreuz-Validierung als auch das Quasi-Optimalit ats- kriterium zu L osungen f uhren, die nahe bei diesem " Knick liegen.

KAPITEL 4. FREDHOLMSCHE INTEGRALGLEICHUNGEN 1. ART 43

Die L-Kurve hat in letzter Zeit mehr Aufmerksamkeit gewonnen als ein eigenst andiges Verfahren zur Wahl des Regularisierungsparameters. Der Gedanke ist, den Wert des Regularisierungsparameters zu w ahlen, der genau dem " Knick in der L-Kurve e n tspricht. Hier deniert man den " Knick als den Punkt, an dem die Kurve die maximale Kr ummung aufweist. Falls die L-Kurve diskret ist beispielsweise f ur die abgeschnittene SVD und iterative V erfahren, n ahert man die L-Kurve durch eine Spline-Funktion an und berechnet die Kr ummung der Spline-Funktion. Der " Knick k ann dann mit Hilfe eines eindimensionalen Optimierungsverfahrens gefunden werden. Dieser Algorithmus ist mindestens uberlegen,

falls die rechte Seite stark miteinander korrelierte Fehler enth alt. Obwohl die L-Kurve ein sehr n utzliches und in der Praxis stabiles Verfahren f ur die Wahl des Regularisierungsparameters darstellt, fehlt gegenw artig allerdings noch e i n e strenge mathematische Konvergenzuntersuchung, wie sie beispielsweise f ur die Kreuz-Validierung durchgef uhrt wurde.

Kapitel 5

Das Programm CGFFT

Das Programm CGFFT dient der Berechnung der auf einer Leiterplatte ieenden St ome aufgrund der oberhalb der Leiterplatte gemessenen magnetischen Feldst arken. Dabei wird das Verfahren der konjugierten Gradienten verwendet, wie es in Abschnitt 4.4.4 beschrieben wurde, wobei die dabei auftretenden Matrix-Vektor-Produkte nach dem in Abschnitt 3.2 erl auterten Verfahren mit Hilfe der zweidimensionalen FFT berechnet werden. Der Quelltext des Programms bendet sich in Abschnitt D.1.

5.1 Programmablauf

Der Ablauf des Programms vollzieht s i c h in folgenden Schritten:

Deklaration der Konstanten, Unterprogramme und Variablen.

Eingabe des Dateinamens und Laden der Daten.

Berechnung der Matrixelemente nach G l e i c hung 2.34.

Jeweils f ur die Berechnung von I x aus H y und von I y aus H x :

Belegung der Felder f ur die Matrix und die adjungierte Matrix ent-

KAPITEL 5. DAS PROGRAMM CGFFT 45

des richtigen Regularisierungsparameters ein ausgefeilteres Verfah-

nauigkeit benutzt werden. Abspeichern der berechneten Koezienten der Str ome.

Abspeichern einer Diagnosedatei, falls dies gew unscht ist. Die Diagnosedatei kann von dem Programm Unigraph siehe Abschnitt B.5 als Report eingelesen werden. Die in Abschnitt 4.6 beschriebene L-Kurve e n tsteht, wenn man in doppelt logarithmischem Mastab die Norm der L osung uber der Norm des Defekts auftr agt. Daneben sind in der Diagnosedatei die mit einer einfachen N aherung durch Geradenst ucke berechnete Kr ummung der L-Kurve v orhanden sowie der Wert ' r , der in 3, Abschnitt 4 beschrieben ist und der ebenfalls f ur die Wahl des Regularisierungsparameters verwendet werden kann.

5.2 Externe Unterprogramme

In dem Programm werden mehrere externe Unterprogramme verwendet, die alle in der netlib siehe Abschnitt A.2.1 zu nden sind. Der Gr unde hierf ur sind zum einen die Zeitersparnis bei der Programmierung, zum anderen aber auch die h ohere Zuverl assigkeit und oft auch g r oere Ausf uhrungsgeschwindigkeit dieser Programme gegen uber selbst geschriebenen.

Die Ursache f ur die hohe Qualit at solcher Programme liegt darin, da sie in der Regel von Experten auf dem jeweiligen Gebiet geschrieben wurden, die sowohl mit den numerischen Algorithmen gut vertraut sind als auch Programmiererfahrung in der jeweiligen Programmiersprache meist FORTRAN besitzen. Ein weiterer nicht z u u n tersch atzender Punkt ist, da diese Programme dadurch, da sie frei erh altlich sind, eine sehr weite Verbreitung gefunden haben, wodurch sie sehr gr undlich i n v i e l e n v erschiedenen Anwendungen auf ihre Zuverl assigkeit getestet wurden. Im Programm CGFFT wurde daher so oft wie m oglich auf solche externen Unterprogramme zur uckgegrien.

In den folgenden Abschnitten werden die verschiedenen Unterprogramme ein wenig n aher beschrieben werden.

5.2.1 Basic Linear Algebra Subprograms

Die Basic Linear Algebra Subprograms BLAS 25, 26, 27, 28, 2 9 sind eine Sammlung von Unterprogrammen in FORTRAN, die elementare Vektor- und Matrixoperationen, wie zum Beispiel die Addition zweier Vektoren oder die

KAPITEL 5. DAS PROGRAMM CGFFT 46

Multiplikation einer Matrix mit einem Vektor durchf uhren. Die BLAS wurden eingef uhrt zur Vereinfachung der Programmierung sowie zur Verbesserung der Ubertragbarkeit von Programmen auf andere Rechner. Sie bilden inzwischen einen de-facto-Standard und werden von vielen Rechnerherstellern in optimierter Form mitgeliefert. Die Routinen von LAPACK siehe n achster Abschnitt bauen auf den BLAS auf.

Wegen der gr oeren Schnelligkeit der auf den Rechner abgestimmten BLAS sollte man die BLAS aus der netlib nur verwenden, falls der Rechnerhersteller keine optimierten BLAS zur Verf ugung stellt. Bei den am IHF installierten Rechnern vom Typ IBM RS6000 sollten die mitgelieferten BLAS aus der Datei libblibblas.a benutzt werden. Als Dokumentation sind die BLAS in der netlib jedoch sehr gut geeignet, wenn andere Quellen momentan nicht z u r Verf ugung stehen. F ur einige Unterprogramme der BLAS sind Manual-Pages vorhanden, die mit dem Befehl

man Name des Unterprogramms

abgerufen werden k onnen. Bei den am IHF installierten Rechnern vom Typ IBM RS6000 sind die BLAS aus der netlib auerdem im Verzeichnis usrrlocallLAPACKKBLASSSRC vorhanden.

Im Programm CGFFT werden von den BLAS die Unterprogramme DZNRM2, ZCOPY, ZAXPY und ZDSCAL verwendet.

5.2.2 LAPACK

LAPACK 30, 3 1 ist eine Sammlung von Unterprogrammen in FORTRAN, die zur L osung der h augsten Probleme der linearen Algebra dient: L osung von liuberbestimmter Systeme im quadratischen

Mittel, L osung von Eigenwert- und Singul arwertproblemen. LAPACK ist der Nachfolger der Programmsammlungen LINPACK und EISPACK und erweitert deren Funktionalit at. Da LINPACK und EISPACK nicht mehr weiterentwickelt werden, sollten Unterprogramme daraus nicht mehr verwendet werden. Auf den am IHF installierten Rechnern vom Typ IBM RS6000 sind auerdem die Quelltexte der LAPACK-Unterprogramme im Verzeichnis ^{usrrlocallLAPACKKSRC} verf ugbar. F ur alle Unterprogramme aus LAPACK und f ur LAPACK selbst sind Manual-Pages vorhanden, die mit dem Befehl

man Name des Unterprogramms

abgerufen werden k onnen.

Eine gute Zusammenfassung der zur Verf ugung stehenden Unterprogramme ist 31, Anhang A und BB. Auf den am IHF installierten Rechnern vom Typ IBM RS6000 ist 31 als L a T E X-Datei unter dem Namen usrrlocallLAPACKKINSTALLLinstall.tex verf ugbar.

Im Programm CGFFT werden aus LAPACK die Unterprogramme SECOND, ZLASCL, ZLAZRO und ZDRSCL verwendet.

KAPITEL 5. DAS PROGRAMM CGFFT 47

5.2.3 Transactions on Mathematical Software

In den Transactions on Mathematical Software TOMS der Association for Computing Machinery ACM werden regelm aig FORTRAN-Programme f ur die L osung mathematischer Probleme ver oentlicht. Diese Programme sind auch in der netlib im Verzeichnis toms archiviert.

Im Programm CGFFT wird aus dem TOMS-Algorithmus Nr. 691 32 d a s Unterprogramm DQXGS verwendet, das zur numerischen Integration einer Funktion dient. Die Erzeugung der entsprechenden Library ist im Abschnitt C.3 beschrieben.

5.2.4 Golden Oldies

Die Golden Oldies in der netlib sind eine Zusammenstellung von FORTRAN-Unterprogrammen, die weit verbreitet und allgemein anerkannt sind, aber nicht Teil einer Programmsammlung sind.

Im Programm CGFFT wird das Unterprogramm DFFT verwendet, das eine leicht abgewandelte Version des Unterprogramms FFT der Golden Oldies darstellt. Die Erzeugung der entsprechenden Library ist in den Abschnitten C.1 und C.2 beschrieben.

5.3 Compilierung des Programms

Unter der Annahme, da sich die Libraries libfft.a und lib691.a im Verzeichnis ihffuserrschmidttlib benden und das Programm selbst in einer Datei namens cgfft.f steht, mu der Befehl f ur die Compilierung des Programms auf den am IHF installierten Rechnern vom Typ IBM RS6000

xlf -O cgfft.f -ocgfft -llapack -lblas -LLihffuserrschmidttlib -l691 -lfft

lauten.

5.4 Das Hilfsprogramm DAT2INP

Das Meprogramm aus 21 liefert f ur die gemessenen magnetischen Feldst arken zwei getrennte Dateien, in denen jeweils Betrag und Phase der x- b z w . d e r ^y- Komponentengespeichert sind. Das Programm CGFFT erwartet jedoch Real-und Imagin arteil der Werte in einer einzigen Datei. Das Programm DAT2INP dient dazu, die beiden Datenformate entsprechend ineinander umzuwandeln. Der Quelltext zu DAT2INP bendet sich i m A b s c hnitt D.2.

Das Programm erwartet, da sich die Dateien mit den Mewerten im Verzeichnis ihffuserrwiedmannndata benden, wobei die Namen der zusammengeh orenden Dateien den gleichen Anfang und die Endungen hx.dat und hy.dat haben m ussen diese Vorgaben lassen sich selbstverst andlich problemlos im Quelltext andern. Die umgewandelte Datei wird im aktuellen Verzeichnis

KAPITEL 5. DAS PROGRAMM CGFFT 48

erzeugt und tr agt die Endung .inp, wie sie von CGFFT erwartet wird. Das Programm wird auf den am IHF installierten Rechnern vom Typ IBM RS6000 mit dem Befehl

xlf -O dat2inp.f -odat2inp

compiliert.

Kapitel 6

Ergebnisse

Um den entwickelten Algorithmus in der Praxis zu testen, wurden Berechnunuber

einer Leiterplatte gemessen worden waren. Die Leiterbahnen darauf hatten eine spiralf ormige Struktur, weil mit dieser die r aumliche Aufl osung des Verfahrens besonders gut dargestellt werden kann. Die Messung erfolgte bei der Frequenz f = 10 MHz, die Mepunkte hatten einen Abstand von a = 2 mm. Der Mebereich umfate 125 Mepunkte in x-Richtung und 75 Mepunkte in ^y-Richtung, entsprech e n d e i n e m B e r e i c h v on 25 cm Ê 15 cm.

Die gemessenen magnetischen Feldst arken in x- u n d y-Richtung zeigen die Abblidungen 6.1 und 6.2. Da gegenw artig die Sensoren f ur die magnetische uber die relatiuber ihren absoluten Wert

m oglich. Bei Messung mu beachtet werden, da die Werte f ur die x- u n d ^y- ubereinanderliegenund nicht r aumlich gegeneinander verschoben sind, da dies ung unstige Auswirkungen auf die berechneten Str ome hat, insbesondere dann, wenn auch die elektrischen Feldst arken in die Berechnung einbezogen werden.

Die Abbildung 6.3 zeigt die Daten aus der Diagnosedatei, wenn das Programm CGFFT mit n = 10000 Iterationsschritten durchgef uhrt wird. Selbstverst andlich wird man in der Praxis nie so viele Schritte durchf uhren, aber die Darstellung der L-Kurve ist bei einer gr oeren Anzahl von Schritten deutlicher ausgepr agt. Man beachte, da mit einer wachsenden Anzahl von Iterationsschritten die Norm des Defekts ^{~ d} k siehe Abschnitt 4.4.4 abnimmt, w ahrend

die Norm der L osung ^{~ x} k w achst. Mit einer zunehmenden Anzahl an Iterationsschritten setzt sich die L-Kurve a l s o n a c h links oben und nicht e t wa n a c h rechts fort. Zu Beginn der Iteration liegen die den einzelnen Iterationsschritten entsprechenden Kurvenpunkte noch recht w eit auseinander, im Verlauf der Iteration n ahern sie sich j e d o c h immer mehr einander an. Der Knick i n d e r L -Kurve wird hier nach e t wa 60 Iterationsschritten erreicht. Da an dieser Stelle die Kurvenpunkte bereits recht d i c ht beieinanderliegen, haben einige Schritte mehr oder weniger keine gr oeren Auswirkungen.

49

KAPITEL 6. ERGEBNISSE 50

Gemessene magnetische Feldstärken

Abbildung 6.1: Gemessene magnetische Feldst arken in x-Richtung

KAPITEL 6. ERGEBNISSE 51

Gemessene magnetische Feldstärken

Abbildung 6.2: Gemessene magnetische Feldst arken in y-Richtung

KAPITEL 6. ERGEBNISSE 52

Ströme in x-Richung Ströme in y-Richtung

L-Kurve L-Kurve

^-1 10

^-2 10 ^{Norm der Lösung -3} 10

^-4 10

^-5 10 0.05 Krümmung

0.0

^{-0.05 -1} 10

^-2 10

-3 ^r 10

^-4 10

^-5 10

Abbildung 6.3: Die L-Kurve und der Sch atzer ' r bei 10000 Iterationsschritten

KAPITEL 6. ERGEBNISSE 53

Berechnete Ströme, f=10 MHz

Grundlage: gemessene magnetische Feldstärken oberhalb der Leiterplatte

Die Schaubilder in der Mitte zeigen die Kr ummung der L-Kurve, die aufgrund einer einfachen N aherung mittels Geradenst ucken errechnet wurde. Man sieht, da dieselbe rechts vom Knick der Kurve stark zu oszillieren neigt, w ahrend sie oberhalb des Knicks praktisch auf Null zur uckgeht. Um einen Kr ummungsverlauf zu erhalten, der besser dem optischen Gesamteindruck der L-Kurve e n tspricht, d urfte es sich empfehlen, diese durch Splines anzun ahern.

Auf den unteren beiden Schaubildern ist der in 3, A b s c hnitt 4 beschriebene Sch atzer ' r dargestellt. Anders als aufgrund 3 z u e r w arten, enth alt die Kurve in diesem Fall jedoch k ein globales Minimum, so da eine Auswahl aufgrund dieses Kriteriums hier nicht m oglich ist. Allerdings knickt der Graph von ' r in der N ahe der optimalen L osung ahnlich wie die L-Kurve s c harf nach o b e n a b , so da diese Tatsache gegebenenfalls f ur die Auswahl des Regularisierungsparameters genutzt werden k onnte.

Abbildung 6.4 zeigt schlielich die aufgrund der gemessenen magnetischen Feldst arken berechneten Str ome, wobei hier das Verfahren der konjugierten Gradienten nach 60 Iterationsschritten abgebrochen wurde. Der Verlauf des Stromusses ist deutlich z u e r k ennen. Die Rechenzeit des Programms betr agt bei 60 Iterationsschritten auf dem am IHF installierten Rechnern vom Typ IBM RS6000 etwa 8 M i n uten.

Kapitel 7

Versuche mit anderen

Ansatzfunktionen

Zu Beginn der Arbeit wurden Versuche unternommen, statt der auf den Segmenten konstanten Str ome auf den Segmenten lineare Stromverl aufe als Ansatzfunktionen zu verwenden. Bei diesen Versuchen wurde nicht wie sp ater das Verfahren der konjugierten Gradienten, sondern die abgeschnittene Singul arwertzerlegung als Regularisierungsverfahren verwendet. Die Ergebnisse d urften aber ubertragbar sein.

Da in diesem fr uhen Stadium der Arbeit noch k eine brauchbaren Meergebnisse zur Verf ugung standen, wurden diese Modelle mit Werten getestet, die mit Hilfe des Programms FEKO O 2 4 berechnet worden waren. Mefehler wurden dadurch simuliert, da zu den berechneten Feldst arken zuf allige komplexe Werte addiert wurden, deren Betrag standardverteilt war mit einer Standardabweichung von einem Prozent der maximalen berechneten Feldst arke, und deren Phase gleichverteilt war.

Da das Verfahren der abgeschnittenen Singul arwertzerlegung wesentlich mehr Speicherplatz ben otigt als das sp ater verwendete Verfahren der konjugierten Gradienten unter Benutzung der FFT, wurde hier ein groberes Gitter mit einem Abstand von a = 1 cm zwischen den Kreuzungspunkten gew ahlt.

7.1 Das erste Verfahren

Die Ansatzfunktionen f ur die Str ome beim ersten Verfahren zeigt Abbildung 7.1. Jedem Kreuzungspunkt des Gitters wurden drei dreieckf ormige Funktionen zugeordnet, eine in x-Richtung, eine in y-Richtung und eine, die von der ^x- Richtungin die y-Richtung wechselt. Mit diesen Funktionen, deren Koezienten zu bestimmen sind, lassen sich alle auf dem Gitter st uckweise linearen Stromverl aufe darstellen, bei denen die Summe der zuieenden Str ome jedes Kreuzungspunkts zu jedem Zeitpunkt gleich der der abieenden ist, wo a n d e n Kreuzungspunkten also keine Punktladungen vorhanden sind.

54

KAPITEL 7. VERSUCHE MIT ANDEREN ANSATZFUNKTIONEN 55

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

berechnete Ströme, f=10kHz

Abbildung 7.2: Berechnete Str ome f ur das erste Verfahren, exakte Werte

KAPITEL 7. VERSUCHE MIT ANDEREN ANSATZFUNKTIONEN 56

berechnete Ströme, f=10kHz

1. Verfahren, 1% Datenfehler, mit Regularisierung

Abbildung 7.3: Berechnete Str ome f ur das erste Verfahren, Daten mit einem Prozent F ehleranteil, mit Regularisierung

Verwendet man mit diesen Ansatzfunktionen die genauen berechneten Werte f ur die Feldst arken, so erh alt man ein ausgezeichnetes Ergebnis wie in Abbildung 7.2. Das Bild andert sich j e d o c h sofort, sobald zu den Werten die erw ahnten Fehler von einem Prozent des Maximalwerts hinzugef ugt werden. Ohne Regularisierung erh alt man hier ein vollst andig nutzloses Ergebnis, die Koefzienten f ur die Str ome oszillieren sehr stark und haben keine Aussagekraft mehr. Wendet man nun jedoch ein Regularisierungsverfahren an, so erh alt man stets einen charakteristischen, s agezahnf ormigen Verlauf wie in Abbildung 7.3 f ur die berechneten Str ome.

Dies l at sich dadurch erkl aren, da durch die oben getroene Auswahl der Ansatzfunktionen eine gewisse Asymmetrie in das Problem gebracht wird. Das Regularisierungsverfahren d ampft die Koezienten aller drei Ansatzfunktionen in gleichem Mae. Da jedoch die eine H alfte der Segmente diejenigen, die in der die Richtung wechselnden Ansatzfunktion enthalten sind in zwei Ansatzfunktionen vorkommt, w ahrend die andere H alfte nur in einer auftaucht, wird diese erste H alfte entsprechend st arker ged ampft als die zweite. Daraus folgt, da die obige Auswahl der Ansatzfunktionen in dieser Form f ur ein Regulari- sierungsverfahren nicht geeignet ist.

KAPITEL 7. VERSUCHE MIT ANDEREN ANSATZFUNKTIONEN 57

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Abbildung 7.4: Ansatzfunktionen f ur das zweite Verfahren

7.2 Das zweite Verfahren

Die Ansatzfunktionen f ur die Str ome beim zweiten Verfahren zeigt Abbildung 7.4. Jedem Kreuzungspunkt werden vier dreieckf ormige Funktionen zugeordnet, wobei es sich diesmal um rechtwinklige und nicht wie beim ersten Verfahren um gleichschenklige Dreiecke handelt. Durch diese Wahl der Ansatzfunktionen ist die Symmetrie wiederhergestellt, allerdings ist es nun nicht mehr automatisch g e w ahrleistet, da f ur jeden Kreuzungspunkt die Summe der zuieenden und abieenden Str ome zu jedem Zeitpunkt gleich gro ist. Es m ussen daher an den Kreuzungspunkten zeitharmonische Punktladungen angenommen werden, um die Kontinuit atsbedingung zu erf ullen, wie das auch b e i dem sonst in dieser Arbeit verwendeten Modell mit den auf einem Segment jeweils konstanten Str omen der Fall ist.

Abbildung 7.5 zeigt die L osung f ur den Fall aus Abbildung 7.3, die man bei dieser Wahl der Ansatzfunktionen erh alt. Die durch die Unsymmetrie verursachte S agezahnstruktur ist nun verschwunden, allerdings zeigen sich dort, wo der Stromverlauf seine Richtung andert, deutliche Absenkungen, die durch d i e Regularisierung verursacht w erden. Solche Absenkungen treten zwar auch b e i der Verwendung konstanter Funktionen auf, sie sind dort aber nicht s o a u s g epr agt.

Insgesamt b r a c hte die Verwendung linearer Funktionen statt konstanter Funktionen keine so bedeutenden Verbesserungen, da diese den durch d i e erh ohte Anzahl der Koezienten verursachten gr oeren numerischen Aufwand rechtfertigen w urden. Dieser Aufwand sollte besser in eine Erh ohung der r aum- lichen Au osung investiert werden.

  KAPITEL 7. VERSUCHE MIT ANDEREN ANSATZFUNKTIONEN  58

  berechnete Ströme, f 10kHz  

2.  Verfahren, 1 Datenfehler, mit Regularisierung  

I  in A  

  Über 2.5 10 -6  

2.  4 10 -6 - 2.5 10 -6  

2.  3 10 -6 - 2.4 10 -6  

2.  2 10 -6 - 2.3 10 -6  

2.  1 10 -6 - 2.2 10 -6  

2.  0 10 -6 - 2.1 10 -6  

1.  9 10 -6 - 2.0 10 -6  

1.  8 10 -6 - 1.9 10 -6  

1.  7 10 -6 - 1.8 10 -6  

1.  6 10 -6 - 1.7 10 -6  

1.  5 10 -6 - 1.6 10 -6  

1.  4 10 -6 - 1.5 10 -6  

1.  3 10 -6 - 1.4 10 -6  

1.  2 10 -6 - 1.3 10 -6  

1.  1 10 -6 - 1.2 10 -6  

1.  0 10 -6 - 1.1 10 -6  

9.  0 10 -7 - 1.0 10 -6  

8.  0 10 -7 - 9.0 10 -7  

7.  0 10 -7 - 8.0 10 -7  

6.  0 10 -7 - 7.0 10 -7  

5.  0 10 -7 - 6.0 10 -7  

4.  0 10 -7 - 5.0 10 -7  

3.  0 10 -7 - 4.0 10 -7  

2.  0 10 -7 - 3.0 10 -7  

1.  0 10 -7 - 2.0 10 -7  

  Unter 1.0 10 -7  

  Abbildung 7.5: Berechnete Str ome f ur das zweite Verfahren, Daten mit einem  

  Prozent F ehleranteil, mit Regularisierung  

Kapitel 8

Zusammenfassung und

Ausblick

Das Verfahren der konjugierten Gradienten erlaubt es zusammen mit dem auf der zweidimensionalen FFT basierenden Algorithmus zur Matrix-Vektor-Multiplikation von Block-Toeplitz-Toeplitz-Block-Matrizen, die Str ome auf der Leiterplatte aus den gemessenen Feldst arken auf sehr eziente Weise zu berechnen. Gleichzeitig wird bei diesem Verfahren eine Regularisierung durchgef uhrt, wie sie bei einem solchen schlecht gestellten inversen Problem unabdingbar ist. Die Auswahl des Regularisierungsparameters, in diesem Fall also die Anzahl der Iterationsschritte, nach denen das Verfahren der konjugierten Gradienten abgebrochen wird, geschieht mit Hilfe des heuristischen Kriteriums der L-Kurve e 7 , f ur das zwar noch k eine streng mathematischen Konvergenzuntersuchungen vorliegen, das sich aber in der Praxis bestens bew ahrt hat.

Um die f ur die Berechnung der Str ome notwendige Rechenzeit noch w eiter zu verringern, bietet es sich an, eine Pr akonditionierung durchzuf uhren, so wie s i e i n n 2 b e s c hrieben ist. Ferner d urfte es angebracht sein, den gegenw artig recht einfach gehaltenen Algorithmus zur Bestimmung des Punktes der st arksten Kr ummung der L-Kurve w eiter zu verfeinern, etwa d u r c h die Verwendung von Splines. Weiterhin k onnen zuk unftig auch w eitere Mewerte, wie beispielsweise elektrische Feldst arken, in die Berechnung miteinbezogen werden, um dadurch m oglicherweise noch genauere Ergebnisse zu erhalten.

Dar uber hinaus sollte man in Erw agung ziehen, die Sensoren f ur die magnetische und die elektrische Feldst arke z u e i c hen, um nicht n ur qualitative, sonuber die auf der Leiterplatte ieenden Str ome dern auch quantitative Aussagen

machen zu k onnen. Falls man die Sensoren noch genauer ausmit, kann man auerdem noch die Empfangscharakteristik der Sensoren in die Rechnung einbeziehen, was aber vermutlich n ur geringf ugige Verbesserungen ergeben d urfte. uberlegen, eine andere Anordnung der Mepunkte oder andere Ansatzfunktionen f ur die Str ome zu verwenden, beispielsweise Fl achenstr ome statt den hier verwendeten Linienst omen. Allerdings werden sich der Charakter des Problems und seine L osung durch s o l c he Variationen voraussichtlich n ur wenig andern. Der gr ote Zuwachs an Genauigkeit d urfte sich w ohl

59

KAPITEL 8. ZUSAMMENFASSUNG UND AUSBLICK 60

durch e i n e V erfeinerung der r aumlichen Au osung sowohl bei der Messung als auch bei der rechnerischen Auswertung erreichen lassen.

Anhang A

Informationen und

Programme aus den

Datennetzen

Mit der st andigen Ausdehnung der Datennetze ist es wesentlich l e i c hter geworden, notwendige Informationen zu erhalten und mit anderen Wissenschaftlern aus dem gleichen Forschungsbereich i n V erbindung zu treten. Dar uber hinaus sind viele Programme zu numerischen Problemen frei erh altlich, die nicht zuletzt wegen ihrer weltweiten h augen Verwendung oft von hervorragender Qualit at und sehr zuverl assig sind. Einige Zugrism oglichkeiten auf diese Informationen und Programme werden im folgenden beschrieben.

A.1 Informationen

A.1.1 IPNet

IPNet ist ein Netz f ur Forscher, die an inversen oder schlecht gestellten Problemen arbeiten. Es hat zum Ziel, den Informationsaustausch z w i s c hen den Wissenschaftlern, die in diesen Gebieten arbeiten, zu f ordern, einen Informaur Notizen und wissenschaftliche Fragen von all- " IPNet Digest f

gemeinem Interessse herauszugeben, und eine zentrale Anlaufstelle f ur aktuelle E-Mail-Adressen der Mitglieder bereitzustellen. Die Organisatorin und gegenw artige Herausgeberin des IPNet Digest ist Patricia Lamm von der Michigan State University.

Der IPNet Digest wird in der Regel einmal monatlich per E-Mail verschickt und enth alt alle im vorherigen Monat eingegangenen Beitr age der Mitglieder. Mit dem Befehl

echo help|mail ipnet-request@math.msu.edu

erh alt man eine E-Mail, in der alle weiteren Informationen zum IPNet enthalten sind.

61

ANHANG A. INFORM. UND PROGR. AUS DEN DATENNETZEN 62

A.1.2 NA-NET

NA-NET erf ullt zwei Funktionen: zum einen dient es zum einfachen Informationsaustausch der Mitglieder durch eine entsprechende Weiterleitung von E-Mails, zum anderen wird etwa w ochentlich per E-Mail der NA-Net News Digest versandt, der Mitteilungen und Anfragen von allgemeinem Interesse aus dem Gebiet der numerischen Mathematik und der mathematischen Software enth alt, ubermittelt wurden. die von den Mitgliedern

NA-NET wurde an der Stanford University v on Gene Golub gegr undet. Der gegenw artige Leiter des NA-NET ist Jack Dongarra, der Herausgeber des NA-NET News Digest ist Cleve Moler. Mit dem Befehl

echo|mail na.help@na-net.ornl.gov

erh alt man eine E-Mail, in der alle weiteren Informationen zum NA-NET enthalten sind.

A.1.3 Usenet Network News

Usenet Network News ist ein weltweites Diskussionsforum, in den die Beitr age uber die Rechnernetze verbreitet werden. Das System ist nach der Benutzer

sogenannten Newsgroups organisiert, in denen jeweils bestimmte Themen behandelt werden. Der Zugang erfolgt mit Hilfe der Programme tin oder ^trn. F ur das Thema dieser Arbeit waren vor allem die Newsgroups

comp.lang.fortran sci.math.num-analysis sci.math

interessant. F ur Neulinge sind auerdem die Newsgroups

news.announce.newusers news.newusers.questions news.answers sci.answers comp.answers

empfehlenswert.

Bestimmte Fragen, die in einer Newsgroup immer wieder gestellt werden, ur Frequently Asked Que- in der Regel in einer Sammlung namens FAQ f stions einschlielich der jeweiligen Antworten zusammengestellt. Die FAQ w erden normalerweise einmal monatlich in die betreende Newsgroup eingespielt, auerdem enthalten die *.answers Newsgroups Sammlungen der FAQs, wobei die Newsgroup news.answers alle FAQ e n th alt und die anderen die FAQ d e r ur wissenschaftliche Themen jeweiligen Hierarchie, also sci.answers die FAQ f ur Themen, die Computer betreen. Regelm ai-und comp.answers die FAQ f

ge Beitr age wie die FAQ sind auerdem auf dem Rechner rtfm.mit.edu im

ANHANG A. INFORM. UND PROGR. AUS DEN DATENNETZEN 63

Verzeichnis pubbusenet archiviert und k onnen von dort mit anonymem ^ftp abgerufen werden beim Einloggen den Namen anonymous oder ftp angeben.

Bevor man eine Frage in einer Newsgroup stellt, sollte man zun achst die FAQ und eine Anzahl anderer Beitr age lesen, um mit den Sitten und Gebr auchen in der betreenden Newsgroup vertraut zu sein.

A.1.4 World Wide Web und Mosaic

Das World Wide Web WWW ist ein weltweites Netz von miteinander verbundenen Text-, Bild- und Tondokumenten, die untereinander durch Q u e r v erweise, sogenannte Hyperlinks, verkn upft sind. F ur den Benutzer ist es dabei gleichg ultig, wo auf der Welt ein solches Dokument gespeichert ist, es gen ugt, den entsprechenden Querverweis anzuw ahlen. Die Idee f ur das WWW entstand am Conseil Europ een pour la Recherche Nucl eaire CERN in Genf.

Mosaic ist ein Programm, das den Zugang zum WWW erm oglicht. Es wurde am National Center for Supercomputing Applications NCSA in Champaign, Illinois USA entwickelt. F ur die Benutzung ist ein Bildschirm erforderlich, auf dem die Grakober ache X l auftt das Programm wird dort mit dem Befehl

Mosaic &

aufgerufen.

Die Adresse eines Dokuments im WWW wird mit dem sogenannten Uniform Resource Locator URL angegeben, der das Ubertragungsprotokoll, den

Rechner, das Verzeichnis und den Namen des Dokuments angibt. Ein Dokument, dessen URL man kennt, kann man anw ahlen, indem man in Mosaic das Feld Open... anklickt und dann den entsprechenden URL eingibt. Die URLs einiger besonders n utzlicher Seiten werden im folgenden angegeben werden.

Eine gute Einf uhrung in WWW und Mosaic ndet man ausgehend von der Seite, die beim Aufruf von Mosaic erscheint. Weitergehende Informationen zu WWW und Mosaic sind auf den Seiten

http:::info.cern.chhhypertexttWWWTheProject.html http:::www.ncsa.uiuc.eduuSDGGSoftwareeMosaiccDocsshelp-about.html

zu nden.

Gute Ausgangspunkte f ur mathematische Themen sind die Seiten

http:::euclid.math.fsu.eduScienceemath.html http:::www.csc.fiimath_topicssGeneral.html

utzliche Seite ist Eine auerordentlich n

http:::cui_www.unige.chhmeta-index.html

Hier sind viele verschiedene Suchfunktionen zusammengefat, die es einem erlauben, innerhalb kurzer Zeit Informationen zu einem bestimmten Thema zu nden.

ANHANG A. INFORM. UND PROGR. AUS DEN DATENNETZEN 64

A.1.5 Electronic Transactions on Numerical Analysis

Die Electronic Transactions on Numerical Analysis ETNA sind eine elektronische Zeitschrift, die sich mit Themen der Numerischen Mathematik befat. Sie wird viertelj ahrlich herausgegeben von L. Reichel, R. S. Varga und A. Ruttan von der Kent State University, K e n t, Ohio USA. Die einzelnen Artikel werden im PostScript-Format zur Verf ugung gestellt und k onnen entweder auf dem Bildschirm gelesen oder ausgedruckt werden.

uber die Seite Der Zugang zu ETNA erfolgt entweder mit Mosaic

http:::etna.mcs.kent.eduu

uber anonymes ftp an die Adresse etna.mcs.kent.edu. A u c h der Zugang uber einen Mail-Server ist m oglich. N ahere Informationen hier uber sind mit dem Befehl

echo send index|mail mailer@etna.mcs.kent.edu

zu erhalten. Artikel in ETNA, die m oglicherweise f ur die zuk unftige Arbeit interessant s e i n k onnten, sind 4 und 12. Es d urfte sich lohnen, auch zuk unftig die dort erscheinenden Artikel im Auge zu behalten.

A.1.6 Verzeichnisse mit Forschungsberichten

Viele Forschungseinrichtungen stellen ihre Forschungsberichte als PostScript-Dateien der Allgemeinheit zur Verf ugung. Der Zugri auf die Dateien uber anonymes ftp m oglich. Beispielhaft soll hier der Rechner deacon.mthcsc.wfu.edu genannt w erden, wo i m V erzeichnis ^pubbplemmons mehrere Artikel zum Thema des Verfahrens der konjugierten Gradienten auf Grundlage der FFT zu nden sind.

uber das WWW m oglich,

uber entsprechende Querverweise oder mit Hilfe der Suchfunktionen gefunden werden k onnen. Der Zugri auf die Artikel ist mit Mosaic auch sehr komfortabel m oglich, so gelangt man beispielsweise mit dem URL

ftp:::deacon.mthcsc.wfu.eduupubbplemmons

ubrigens als Datei

preg.ps.Z und der Artikel 11 als Datei decon.ps.Z zu nden.

Weitere Quellen, die ebenfalls als PostScript-Datei erh altlich sind, sind 15 auf dem Rechner reports.adm.cs.cmu.edu im Verzeichnis 1994 als Datei CMU-CS-94-125.ps sowie 16 auf dem Rechner ftp.desy.de im Verzeichnis pubbmgghmg als Datei cg.ps, die beide leicht v erst andliche und anschauliche Einf uhrungen in das Verfahren der konjugierten Gradienten darstellen.

ANHANG A. INFORM. UND PROGR. AUS DEN DATENNETZEN 65

A.2 Programme

A.2.1 netlib

Die netlib ist eine Sammlung von mathematischen Programmen, Daten, Dokumenten, Adressenlisten und anderen n utzlichen Informationen. Sie existiert auf einer ganzen Anzahl von Rechnern, von denen der wohl bekannteste von den AT&T Bell Labs in Murray Hill, New Jersey USA unterhalten wird. Die im Rahmen dieser Diplomarbeit verwendeten externen Unterprogramme sind alle in der netlib zu nden.

Mit dem Befehl

echo send index|mail netlib@research.att.com

erh alt man eine E-Mail, die weitere Informationen zur netlib und ein Inhaltsverzeichnis enth alt. Uber anonymes ftp ist die netlib unter der Adresse netlib.att.com zu erreichen. Der komfortabelste Zugang ist mittels Mosaic unter dem URL

ftp:::netlib.att.commnetlibbmasterrreadme.html.Z

m oglich.

A.2.2 Guide to Available Mathematical Software

Der Guide to Available Mathematical Software GAMS erleichtert den Zugri auf die in der netlib vorhandenen Programme wesentlich. Er wurde am National Institute of Standards and Technology NIST in Gaithersburg, Maryland USA entwickelt. Neben den Programmen der netlib enth alt der GAMS noch Dokumentationen zu verschiedenen am NIST installierten Programmen, die Programme selbst sind in der Regel jedoch n i c ht z u g anglich, da es sich hier um kommerzielle Software handelt.

Die Programme sind im GAMS nach mathematischen Teilgebieten geordnet, so da sich e i n e L osung f ur ein bestimmtes Problem sehr leicht nden l at. Dabei mu man allerdings beachten, da gegenw artig zwar alle wichtigeren Programme der netlib in den GAMS aufgenommen wurden, einzelne Teile jedoch noch fehlen.

Der Zugang zum GAMS erfolgt mit dem Befehl

telnet gams.nist.gov

wobei beim Einloggen der Name gams angegeben werden mu. Die komfortabelste M oglichkeit ist auch hier wieder die Verwendung von Mosaic, als URL ist dabei

http:::gams.nist.govv

anzugeben.

ANHANG A. INFORM. UND PROGR. AUS DEN DATENNETZEN 66

A.2.3 Archie

Falls man den Namen eines frei erh altlichen Programms kennt, aber nicht w ei, woher man dieses erhalten kann, ist die Suche mittels Archie angebracht. Dort sind die Verzeichnisse vieler Rechner gespeichert, von denen man mittels anonymem ftp Programme und sonstige Dateien laden kann. Diese Verzeichnisse werden nach dem gew unschten Begri durchsucht und das Ergebnis der Suche anschlieend angezeigt. Mit dieser Information kann man dann das Programm von einem der angezeigten Rechner laden.

Der Zugang zu Archie erfolgt mit dem Befehl

telnet archie.th-darmstadt.de

wobei beim Einloggen der Name archie angegeben werden mu. Auch h i e r i s t jedoch d i e k omfortabelste M oglichkeit die Verwendung von Mosaic, als URL ist dabei

http:::cui_www.unige.chhOSGGarchieplexform.html

anzugeben. Hier kann man anschlieend die angezeigten Programme direkt durch Anklicken auf den eigenen Rechner laden.

Anhang B

utzliche Programmierhilfen

B.1 Emacs FORTRAN-Mode

Emacs ist ein sehr leistungsf ahiger und komfortabler Editor, der auerdem frei erh altlich ist. Er wurde geschrieben von Richard Stallman und der Free Software Foundation. Die jeweils aktuelle Version bendet sich auf dem Rechner prep.ai.mit.edu im Verzeichnis pubbgnu, v on wo man sie mit anonymem ftp laden kann. Auf den am IHF installierten Workstations ist dieser Editor bereits vorhanden. Zu Emacs gibt es eine Manual-Page, die mit dem Befehl man emacs abgerufen werden kann. Auerdem sind im Programm selbst sehr umfangreiche Hilfsfunktionen vorhanden. Bei Problemen, die man anderweitig nicht l osen kann, kann man sich a u c h an die Usenet Network News siehe Abschnitt A.1.3 Newsgroup gnu.emacs.help wenden.

In Emacs ist unter anderem ein FORTRAN-Mode enthalten, der die Erstellung von FORTRAN-Programmen sehr vereinfacht, da er die in FORTRAN so ubernimmt

und auerdem das Einr ucken von Schleifen und bedingt auszuf uhrenden Proubernimmt.

Falls man das FORTRAN-Programm dem Standard entsprechend ublichen Compiler verstehen auch Kleinschreibung, machen aber in der Regel keinen Unterschied zwischen Gro- und Kleinbuchstaben, kann man sich die Datei capslock.el vom Rechner archive.cis.ohio-state.edu aus dem Verzeichnis pubbgnuuemacsselisp-archiveemisc mit anonymem ftp besorgen.

Unter der Annahme, da sich die Datei capslock.el im Verzeichnis ihffuserrschmidttemacs bendet, k onnten n utzliche Einstellungen in der Datei ^.emacs die sich d a n n i m V erzeichnis ihffuserrschmidt benden mu f ur die FORTRAN-Programmierung und das Schreiben von Texten mit L a T E X etwa folgendermaen aussehen:

autoload 'caps-lock-mode ""ihffuserrschmidttemacsscapslock.el" "Toggle caps-lock mode." t

67

UTZLICHE PROGRAMMIERHILFEN 68

setq transient-mark-mode 1

setq-default indent-tabs-mode nil

setq fortran-tab-mode-default nil

setq fortran-continuation-string "&"

setq fortran-comment-line-extra-indent -4

add-hook 'fortran-mode-hook 'lambda setq abbrev-mode 1

add-hook 'fortran-mode-hook 'lambda fortran-auto-fill-mode 1

add-hook 'fortran-mode-hook 'lambda caps-lock-mode 1

make-variable-buffer-local 'line-number-mode

add-hook 'fortran-mode-hook 'lambda line-number-mode 1

setq tex-dvi-view-command "xdvi"

setq tex-dvi-print-command "dvips -otttmp.ps **ghostview -a4 tttmp.psssusrrbinrm tttmp.ps"

setq tex-alt-dvi-print-command "dvips -f *|lp -dps"

setq tex-show-queue-command "lpstat -vps"

setq tex-default-mode 'latex-mode

add-hook 'tex-mode-hook 'lambda local-unset-key """"

ftnchek

B.2

Das Programm ftnchek dient dazu, bestimmte Fehler in FORTRAN-Programmen aufzusp uren, die ein FORTRAN-Compiler normalerweise nicht entdeckt. Solche Fehler sind beispielsweise Variablen, die zwar deklariert, aber nie benutzt werden, Variablen, die vor ihrer ersten Verwendung nicht initialisiert wurden, oder Variablen, die vor ihrer Verwendung nicht deklariert wurden.

Das Programm ftnchek wurde von Dr. Robert Moniot geschrieben und ist frei erh altlich, man ndet es in der netlib siehe Abschnitt A.2.1 im Verzeichnis fortran. Zusammen mit dem Programm erh alt man auch eine ausf uhrliche Anleitung.

ANHANG B. N UTZLICHE PROGRAMMIERHILFEN 69

s2d und d2s

B.3

Die Programme s2d und d2s dienen dazu, FORTRAN-Programme mit Variablen einfacher Genauigke i t i n s o l c he mit Variablen doppelter Genauigkeit umzuwandeln und umgekehrt. Die Programme wurden unter anderem von Jim Meyering entwickelt und sind von der netlib siehe Abschnitt A.2.1 frei erh altlich. Die genaue Erzeugung der Programme ist im Abschnitt C.1 ausf uhrlich beschrieben.

Falls ein mit s2d umgewandeltes Programm oder Unterprogramm implizit deklarierte Variablen einfacher Genauigkeit enthielt, so ist es notwendig, am Beginn des jeweiligen Programms oder Unterprogramms die Zeile

IMPLICIT DOUBLE PRECISION A-H, O-Z

einzuf ugen.

xfig

B.4

Das Programm xfig ist ein recht k omfortables Zeichenprogramm, wobei es den besonderen Vorteil besitzt, da sich die damit erzeugten Bilder sehr gut in L a T E X-Dokumente einbinden lassen. Das Programm wurde urspr unglich v on Supoj Sutanthavibul geschrieben und seither von vielen anderen Autoren entscheidend erweitert. Es ist frei erh altlich, man ndet es auf dem Rechner ^ftp.x.org im Verzeichnis contribbR5fixessxfig-patches und kann es von dort mit anonymem ftp laden. Auf den am IHF installierten Rechnern ist xfig bereits vorhanden. Zu xfig ist eine ausf uhrliche Manual-Page vorhanden, die sich m i t dem Befehl man xfig abrufen l at.

Man kann die Zeichnung unter anderem als L a T E X-picture-Umgebung ausurfte jedoch die Ausgabe als P I CT E X-Macros geeignet sein. geben, besser noch d

Hierzu mu man sich z u n achst die Dateien prepictex.tex, pictex.tex und postpictex.tex mit anonymem ftp besorgen, die man beispielsweise auf dem Rechner ftp.uni-stuttgart.de im Verzeichnis ^{pubbtexxgraphicsspictex} ndet. In der L a T E X-Datei m ussen dann nach d e r A n weisung ^{beginndocumentt} die Anweisungen

input prepictex input pictex input postpictex

stehen. Nach diesen Vorbereitungen kann man die von xfig erzeugten P I CT E X-Dateien mit einem input-Befehl einlesen. Einzelheiten zu P I CT E X ndet man in 34, Kapitel 44.

Der Aufruf von xfig mit geeineten Voreinstellungen f ur L a T E X lautet

xfig -mo -me -P -pw 21 -ph 29.7 -lat -sp -e pictex Dateiname &

ANHANG B. N UTZLICHE PROGRAMMIERHILFEN 70

B.5 Unigraph

Unigraph ist ein Programm, das es erlaubt, auf einfache Weise Schaubilder anhand von Datenwerten zu erstellen. Das Programm ist am IHF auf den Rechnern vom Typ HP 9000 installiert und ist nicht f r e i e r h altlich. Der Aufruf des Programms erfolt mit den Befehlen

. usrrlpppunirass6v3aabaseeuni.profile unigraph -m -d mx11 &

Eine Dokumentation zu Unigraph ist am IHF vorhanden, auerdem sind auch im Programm selbst Hilfsfunktionen eingebaut.

Eine Besonderheit von Unigraph sind die Command-Dateien, die es erlauben, eine Anzahl von Funktionen automatisch auszuf uhren, was besonders dann von Nutzen ist, wenn verschiedene Datens atze immer wieder auf die gleiche Art dargestellt werden sollen. Command-Dateien lassen sich automatisch erstellen, indem man im Men u OPTIONSSDIALOG den Men upunkt ^{Command logging} auf "ON" setzt. Dadurch w erden alle eingegebenen Kommandos in der Datei gespeichert. Hat man die gew unschten Operationen durchgef uhrt, setzt man Command logging wieder auf "OFF". Die so erzeugte Datei kann anschlieend editiert werden und so jeweils den aktuellen Erfordernissen angepat werden.

Im folgenden wird beispielhaft eine Command-Datei dargestellt, die dazu dient, die berechnete Stromverteilung anzuzeigen.

10000 DATA 11000 IMPORT 11155 ASCII 11 File name "..beispiel.sol" 22 First line in file 1 33 Last line in file 9300 66 Conversion filter " ,, " 77 Dataset names "IXRE IXIM" 99 Dataset rows 124 100 Dataset columns 75 122 File layout "VERTICAL" -11 Done 11155 ASCII 22 First line in file 9301 33 Last line in file 18550 77 Dataset names "IYRE IYIM" 99 Dataset rows 125 100 Dataset columns 74 -11 Done 10000 DATA 17000 OPERATIONS 17100 LET "NULLX" "IXRE * 0"

  UTZLICHE PROGRAMMIERHILFEN  

  71  

  17100 LET "NULLY" "IYRE 0"  

  30000 GRAPH  

  30500 CREATE  

11  Expressionndataset "IXRE"  

22  Group name "X"  

33  Graph type "2D VECTOR"  

  199 Y component "NULLX"  

77  X "range1.555124.5,124"  

88  Y "range11175,75"  

  -11 Done  

  38255 VECTOR  

33  Justification "MIDDLE"  

88  Lineeframe color  1

  -11 Done  

  30500 CREATE  

11  Expressionndataset "NULLY"  

22  Group name "Y"  

  199 Y component "IYRE"  

77  X "range111125,125"  

88  Y "range1.55574.5,74"  

  -11 Done  

  38255 VECTOR  

33  Justification "MIDDLE"  

88  Lineeframe color  1

  -11 Done  

  40000 GROUP  

  42500 LEGEND  

  42522 ENABLE "YES"  

  42544 TYPE "VECTOR"  

  42566 POSITION 1030000E 03, 1000000E  03

  42588 LAYOUT  

11  Origin alignment "UPPER LEFT"  

77  Number of decimals  2

88  Floating format "X.X 10"  

  -11 Done  

  40000 GROUP  

  41000 AXES  

  41100 SWITCHES  

44  1ST 2D X-AXIS "OFF"  

88  1ST 2D Y-AXIS "OFF"  

  122 2ND 2D X-AXIS "OFF"  

  166 2ND 2D Y-AXIS "OFF"  

  -11 Done  

  30000 GRAPH  

  32000 SELECT  2

  40000 GROUP  

  41000 AXES  

ANHANG B.  N  

  UTZLICHE PROGRAMMIERHILFEN  

  72  

  41100 SWITCHES  

44  1ST 2D X-AXIS "ON"  

88  1ST 2D Y-AXIS "ON"  

  122 2ND 2D X-AXIS "ON"  

  166 2ND 2D Y-AXIS "ON"  

  -11 Done  

  41200 ATTRIBUTES  

  41244 AXLE  

66  Cross value 0000000E  00

  -11 Done  

  41222 SELECT "1ST 2D Y-AXIS"  

  41244 AXLE  

66  Cross value 0000000E  00

  -11 Done  

  41222 SELECT "2ND 2D X-AXIS"  

  41244 AXLE  

66  Cross value 7600000E  02

  -11 Done  

  41222 SELECT "2ND 2D Y-AXIS"  

  41244 AXLE  

66  Cross value 1.2600000E  02

  -11 Done  

  50000 LAYOUT  

  52000 VIEWPORT  

  52155 SELECT  0

  52500 TITLES  

  52700 MODIFY  

22  Title text "Stromverteilung"  

  -11 Done  

  40000 GROUP  

  40500 LIMITS  

44  X minimum 0000000E  00

55  X maximum 1.2600000E  02

66  Y minimum 0000000E  00

77  Y maximum 7600000E  02

  -11 Done  

  41000 AXES  

  41266 LABELS  

11  Labels "OFF"  

  466 Level 1 labels "OFF"  

  477 Level 2 labels "OFF"  

  488 Level 3 labels "OFF"  

  499 Level 4 labels "OFF"  

  -11 Done  

  41322 TICKMARKS  

11  Major tickmarks "OFF"  

  -11 Done  

  UTZLICHE PROGRAMMIERHILFEN  

  73  

  41344 TICKLINES  

11  Major ticklines "OFF"  

  -11 Done  

  30000 GRAPH  

  38255 VECTOR  

88  Lineeframe color  1

55  Size in "VALUEEPERCENT"  

77  Unit valuee 2E-06  

  -11 Done  

  122 REPAINT  

  Mit dieser Command-Datei werden die Daten der Datei beispiel.sol im  

  aktuellen Verzeichnis geladen und die entsprechende Stromverteilung darge-  

  stellt. Falls diese Datei den Namen strom.com tr agt und sich im beim Program-  

  maufruf aktuellen Verzeichnis bendet, lautet die einzugebende Anweisung nach  

  dem Start von Unigraph einfach strom. Der Wert in der drittletzten Zeile  

77  Unit valuee 2E-06  

  bestimmt den Mastab, mit dem die Daten dargestellt werden und mu jeweils  

  an die aktuellen Werte angepat werden  

Anhang C

Erzeugung der ben otigten

Libraries

F ur das Programm CGFFT werden zwei Unterprogramm-Bibliotheken Libra-

das bei den Libraries lapack und blas der Fall ist. Die Erzeugung dieser Libraries wird im folgenden beschrieben. Alle Angaben beziehen sich dabei auf die Rechner vom Typ IBM RS6000, wie sie am Institut f ur Hochfrequenztechnik installiert sind.

C.1 Erzeugung des Programms s2d

Zu Beginn soll die Erzeugung des Programms s2d beschrieben werden, das

sp ater f ur die Erstellung der Library fft ben otigt wird.

Es wird zun achst mit den Befehlen

mkdir s2d cd s2d

ein neues Verzeichnis erzeugt und dieses zum aktuellen Verzeichnis gemacht. Danach wird mit dem Befehl

ftp netlib.att.com

eine Verbindung zu dem Rechner aufgebaut, auf dem die Programme vorhanden sind. Als Benutzername beim Einloggen mu anonymous oder ftp angegeben werden. Anschlieend folgen die Befehle

cd netlibbfortran bin get f-s2d1.Z get f-s2d2.Z

74

OTIGTEN LIBRARIES 75

get f-s2d3.Z quit uncompress *.Z sh f-s2d1 sh f-s2d2 sh f-s2d3 cd f-s2d-1.1

Damit werden die Dateien auf den eigenen Rechner geholt und entpackt, danach wird in das Verzeichnis mit den entpackten Programmen gewechselt.

Als n achstes mu man mit einem Editor in der Datei Makefile das erste Zeichen # in den Zeilen

#CC = cc #CFLAGS = -g $defs

entfernen sowie ebenfalls mit einem Editor den Inhalt der Dateien ^system.h und proto.h l oschen, so da diese beiden Dateien zwar noch existieren, aber leer sind.

Schlielich w erden mit den Befehlen

make mv s2d ~~bin mv d2s ~~bin

die Programme erzeugt und in das Verzeichnis f ur die ausf uhrbaren Dateien verlegt.

Erzeugung der Library fft

C.2

Es wird zun achst mit den Befehlen

mkdir fft cd fft

ein neues Verzeichnis erzeugt und dieses zum aktuellen Verzeichnis gemacht. Danach wird mit dem Befehl

ftp netlib.att.com

eine Verbindung zu dem Rechner aufgebaut, auf dem die Programme vorhanden sind. Als Benutzername beim Einloggen mu anonymous oder ftp angegeben werden. Anschlieend folgen die Befehle

ANHANG C. ERZEUGUNG DER BEN OTIGTEN LIBRARIES 76

cd netlibbgo bin get fft.f.Z quit uncompress fft.f.Z

mit denen die Datei auf den eigenen Rechner geholt und entpackt wird.

Als n achstes wird mit einem Editor in die Datei fft.f nach der Zeile

maxp=209

die Zeile

maxf=23

eingef ugt und die Zeile

dimensiona1,b1

in

dimensiona*,b*

ur das Leerzeichen. ge andert das Symbol steht d a b e i f

Danach wird mit dem Befehl

s2d fft.f dfft.f

die Version des Unterprogramms f ur doppelte Genauigkeit erzeugt. Anschlieend mu mit einem Editor in der Datei dfft.f die Zeile

subroutineffta,b,ntot,n,nspan,isn

in

subroutinedffta,b,ntot,n,nspan,isn

ge andert und nach dieser Zeile die Zeile

implicitdoubleprecisiona-h,o-z

eingef ugt werden. Auerdem werden die Zeilen

c72=0.30901699437494742d0 s72=0.95105651629515357d0 s120=0.86602540378443865d0 rad=6.2831853071796d0

OTIGTEN LIBRARIES 77

in

c72=0.309016994374947424102293417183d0 s72=0.951056516295153572116439333378d0 s120=0.866025403784438646763723170755d0 rad=6.28318530717958647692528676656d0

ge andert, wodurch die vom Programm ben otigten Konstanten mit einer f ur doppelte Genauigkeit ausreichenden Stellenzahl angegeben werden.

Schlielich wird mit den Befehlen

xlf -cO *.f ar r libfft.a *.o mv libfft.a ~lib

die Library erzeugt und in das f ur die selbst erzeugten Libraries vorgesehene Verzeichnis verlegt. Die Option O beim Compilieren ist hierbei wichtig, da sich dadurch die Zeit f ur die Programmausf uhrung gegen uber einer Compilierung ohne diese Option ann ahernd halbiert.

Erzeugung der Library 691

C.3

Es wird zun achst mit den Befehlen

mkdir 691 cd 691

ein neues Verzeichnis erzeugt und dieses zum aktuellen Verzeichnis gemacht. Danach wird mit dem Befehl

ftp netlib.att.com

eine Verbindung zu dem Rechner aufgebaut, auf dem die Programme vorhanden sind. Als Benutzername beim Einloggen mu anonymous oder ftp angegeben werden. Anschlieend folgen die Befehle

cd netlibbtoms bin get 691.Z quit uncompress 691.Z fsplit 691

Dadurch wird die Datei auf den eigenen Rechner geholt, entpackt und in einzelne Unterprogramme zerlegt.

Als n achstes werden mit den Befehlen

OTIGTEN LIBRARIES 78

cat main001.f f.f dtest.f cat main002.f zzz000.f stest.f rm main001.f main002.f f.f zzz000.f

einige der Dateien neu geordnet. Die Datei main000.f mu mit einem Editor von Hand in eine Datei readme mit dem Anleitungstext und eine Datei dqpsrt.f mit dem entsprechenden Unterprogramm aufgeteilt werden.

In der Datei r1mach.f m ussen anschlieend noch die Zeilen

CDATARMACH1Z00100000 CDATARMACH2Z7FFFFFFF CDATARMACH3Z3B100000 CDATARMACH4Z3C100000 CDATARMACH5Z41134413

durch die Zeilen

DATARMACH1Z00100000 DATARMACH2Z7FFFFFFF DATARMACH3Z3B100000 DATARMACH4Z3C100000 DATARMACH5Z41134413

u n d i n d e r D a t e i d1mach.f die Zeilen

CDATASMALL1,SMALL2Z00100000,Z00000000 CDATALARGE1,LARGE2Z7FFFFFFF,ZFFFFFFFF CDATARIGHT1,RIGHT2Z33100000,Z00000000 CDATADIVER1,DIVER2Z34100000,Z00000000 CDATALOG101,LOG102Z41134413,Z509F79FF

durch die Zeilen

DATASMALL1,SMALL2Z00100000,Z00000000 DATALARGE1,LARGE2Z7FFFFFFF,ZFFFFFFFF DATARIGHT1,RIGHT2Z33100000,Z00000000 DATADIVER1,DIVER2Z34100000,Z00000000 DATALOG101,LOG102Z41134413,Z509F79FF

ersetzt werden.

Danach wird mit den Befehlen

xlf -c *.f rm ?test.o ar r lib691.a *.o

OTIGTEN LIBRARIES 79

die Library erzeugt. Die Option O darf hierbei bei der Compilierung nicht v erwendet werden, da damit das Unterprogramm DQXGS nicht die volle Genauigkeit erreicht, mit man mit Hilfe der Testprogramme feststellen kann.

Schlielich w erden mit den Befehlen

xlf -O stest.f -ostest -L . -l691 xlf -O dtest.f -odtest -L . -l691

die Testprogramme erzeugt und deren Ergebnis mit

..stest|more ..dtest|more

betrachtet. Wenn die Ergebnisse in Ordnung sind, wird zuletzt die erzeugte Library mit dem Befehl

mv lib691.a ~lib

in das f ur die selbst erzeugten Libraries vorgesehene Verzeichnis verlegt.

Anhang D

Programm-Quelltexte

D.1 Programm cgfft.f

PROGRAM CGFFT

************************************************************************ * *

* Programm: CGFFT * Autor: * Betreuer: Georg Faessler * Datum: * * Das Programm berechnet aus den ueber einer Platine gemessenen * magnetischen Feldern die auf der Platine fliessenden Stroeme. *

* Aufgrund der besonderen Struktur der Matrix, die bei einem solchen * regelmaessig diskretisierten zweidimensionalen Problem entsteht *

* Toeplitz-Block Block-Toeplitz, auch doppelt Toeplitz genannt, *

* lassen sich Produkte der Matrix mit einem Vektor sehr zeit- und *

* speicherplatzsparend mit Hilfe der zweidimensionalen FFT berechnen. * * *

* Damit bietet sich ein iteratives Verfahren zur Loesung der Matrix-

* gleichung an. Hier wurde das Verfahren der konjugierten Gradienten, * * angewandt auf die Normalengleichungen, gewaehlt. * * *

* Wegen der bei der Messung der Felder unweigerlich auftretenden Mess- * * fehler ist es guenstig, ein regularisierendes Verfahren zu verwenden.* * Dies wird durch den Abbruch des Verfahrens an einer geeigneten *

* Stelle erreicht. Das Kriterium hierfuer ist die Stelle der staerk-

* sten Kruemmung der L-Kurve, wobei die Tatsache ausgenutzt wird, *

* dass die Kruemmung an dieser Stelle stark zu oszillieren neigt. * * *

************************************************************************

C In diesem Programm werden folgende Befehle verwendet, die nicht dem C FORTRAN77-Standard entsprechen, jedoch von allen ueblichen Compilern C verstanden werden: C

80

ANHANG D. PROGRAMM-QUELLTEXTE

81

C IMPLICIT NONE fuer die Abschaltung der impliziten Variablendeklaration C DOUBLE COMPLEX fuer komplexe Variablen doppelter Genauigkeit C DCONJG als intrinsische Funktion C DO ohne Zeilennummer ... ENDDO

C FORMAT..., $ fuer die Unterdrueckung des Zeilenvorschubs

C Massnahme gegen Tippfehler bei Variablen C alle Variablen muessen explizit deklariert werden IMPLICIT NONE

C Angabe, ob eine Diagnose-Datei ausgegeben werden soll LOGICAL DIAG PARAMETER DIAG=.TRUE.

C Angabe, ob die Loesung mit Hilfe der L-Kurve ausgewaehlt werden soll C anderenfalls wird das Ergebnis nach N Iterationsschritten ausgegeben LOGICAL LCURVE PARAMETER LCURVE=.FALSE.

C Anzahl der Iterationsschritte INTEGER N PARAMETER N=60 INTEGER NM1 PARAMETER NM1=N-1

C Frequenz in Hz, fuer die die Matrix berechnet wird DOUBLE PRECISION F PARAMETER F=1D7

C Abstand der Messpunkte in m DOUBLE PRECISION A PARAMETER A=2D-3

C Hoehe in m, in der die Felder gemessen werden DOUBLE PRECISION Z PARAMETER Z=1D-2

C Bereich, fuer den die Matrix berechnet wird: C Feld mit XMAX*YMAX Messwerten XMAX YMAX INTEGER XMAX, YMAX PARAMETER XMAX=125, YMAX=75 INTEGER XMAXP1, YMAXP1 PARAMETER XMAXP1=XMAX+1, YMAXP1=YMAX+1 INTEGER XMAXM1, YMAXM1 PARAMETER XMAXM1=XMAX-1, YMAXM1=YMAX-1 INTEGER XMAXM2 PARAMETER XMAXM2=XMAX-2 INTEGER XMAX2 PARAMETER XMAX2=2*XMAX INTEGER XMA2M1 PARAMETER XMA2M1=XMAX2-1 INTEGER XMA2M2 PARAMETER XMA2M2=XMAX2-2

ANHANG D. PROGRAMM-QUELLTEXTE

82

INTEGER YMAX2 PARAMETER YMAX2=2*YMAX INTEGER YMA2M1 PARAMETER YMA2M1=YMAX2-1 INTEGER MN PARAMETER MN=XMAX*YMAX INTEGER INX, INY

PARAMETER INX=XMAXM1*YMAX, INY=XMAX*YMAXM1

C Konstanten DOUBLE COMPLEX ZONE PARAMETER ZONE=1D0, 0D0 DOUBLE COMPLEX J PARAMETER J=0D0, 1D0 DOUBLE COMPLEX ZNULL PARAMETER ZNULL=0D0, 0D0 DOUBLE PRECISION PI PARAMETER PI=3.1415926535897932384626433D0 DOUBLE PRECISION C0 PARAMETER C0=299792458D0 DOUBLE PRECISION BETA0 PARAMETER BETA0=2D0*PI*F*AAC0

C Groesse der Arbeitsbereiche fuer DQXGS INTEGER LIMIT, LENIW, LENW PARAMETER LIMIT=100, LENIW=300, LENW=4600

C nicht standardgemaesse interne Funktion DOUBLE COMPLEX DCONJG INTRINSIC DCONJG

C Funktionen und Unterprogramme aus BLAS, LAPACK und TOMS 691 REAL SECOND EXTERNAL SECOND DOUBLE PRECISION DZNRM2 EXTERNAL DZNRM2

EXTERNAL DQXGS EXTERNAL ZLASCL EXTERNAL ZCOPY EXTERNAL ZLAZRO EXTERNAL ZDRSCL EXTERNAL ZAXPY EXTERNAL ZDSCAL

C eigene Unterprogramme

C diese verwenden das Unterprogramm FFT aus netlibbgo EXTERNAL ZFFT2D EXTERNAL ZFFTMV

C Funktionen zur Berechnung der Integrale DOUBLE PRECISION FRE, FIM EXTERNAL FRE, FIM

ANHANG D. PROGRAMM-QUELLTEXTE

83

C Dimensionierung der Arbeitsvektoren fuer DQXGS INTEGER IWORKLENIW DOUBLE PRECISION QWORKLENW

C Dimensionierung der Matrizen

DOUBLE COMPLEX MFT0:XMA2M1, 0:YMA2M1, MHFT0:XMA2M1, 0:YMA2M1 DOUBLE PRECISION MFT22, 0:XMA2M1, 0:YMA2M1 DOUBLE PRECISION MHFT22, 0:XMA2M1, 0:YMA2M1 EQUIVALENCE MFT, MFT2, MHFT, MHFT2

DOUBLE COMPLEX HXMN, HYMN DOUBLE COMPLEX IXINX, IYINY DOUBLE COMPLEX WORKXMAX2, YMAX2 DOUBLE COMPLEX K-XMAXM2:XMAXM1, 0:XMAXM1

C Deklaration der Variablen fuer das Verfahren der konjugierten Gradienten DOUBLE COMPLEX XINX DOUBLE COMPLEX RINX DOUBLE COMPLEX SINX DOUBLE COMPLEX DMN DOUBLE COMPLEX ASMN DOUBLE PRECISION ALPHA, ALALT, BETA DOUBLE COMPLEX ZALPHA DOUBLE PRECISION RHO, RHODIF DOUBLE PRECISION RNORM, RNNEU DOUBLE PRECISION XNORMXN, XNORMYN DOUBLE PRECISION DNORMXN, DNORMYN DOUBLE PRECISION PHIXN, PHIYN DOUBLE PRECISION XNMLOGN, DNMLOGN DOUBLE COMPLEX XALTINX, 2 INTEGER L

C Deklaration der Variablen fuer die Kruemmung DOUBLE PRECISION CURVX2:NM1, CURVY2:NM1 DOUBLE PRECISION ANGLE1:NM1 DOUBLE PRECISION CURMAX INTEGER ICUMAX LOGICAL CURNEG

C Deklaration der sonstigen Variablen CHARACTER*20 DATEI DOUBLE PRECISION DX, DY, DZ, BETAI DOUBLE PRECISION IRE, IIM DOUBLE PRECISION ERREST INTEGER INFO, LAST INTEGER L1, L2 INTEGER IDUMMY REAL STIME

C COMMON-Block zur Uebergabe der Parameter an FRE und FIM COMMON FPARAMM DX, DY, DZ, BETAI

ANHANG D. PROGRAMM-QUELLTEXTE

84

C Formatangaben fuer die Ausgabe der CPU-Zeit, der Iterationsschritte C und die Eingabe des Dateinamens 1 FORMAT T2, A, F10.2 2 FORMAT T2, A, I3 3 FORMAT T2, A, $

C Messung der CPU-Zeit starten STIME=SECOND

C Eingabe des Dateinamens PRINT *

PRINT *, 'Programm CGFFT - Konjugierte Gradienten mit FFT' PRINT *

PRINT 3, 'Dateiname ohne Erweiterung: ' READ *, 'A' DATEI PRINT *

PRINT *, 'Laden der Daten'

C Laden der Daten

OPEN 1, FILE=DATEI1:INDEXDATEI,' '-1'.inp', & FORM='UNFORMATTED' READ 1 HX, HY CLOSE 1

PRINT 1, 'CPU-Zeit in Sekunden: ', SECOND-STIME PRINT *

PRINT *, 'Berechnung der Matrixelemente'

C Berechnung der Matrixelemente

C DZ gibt die relative Hoehe der Messpunkte, C bezogen auf den Gitterabstand, an DZ=ZZA BETAI=BETA0

C Berechnung der Integrale DO L1=1, XMAXM1 DX=DBLEL1 & &

ANHANG D. PROGRAMM-QUELLTEXTE

85

ENDDO ENDDO

C Kopieren der Werte fuer die negativen Abstaende C der mittlere Wert tritt doppelt auf, deshalb waere die C Verwendung der Betragsfunktion recht umstaendlich DO L1=0, XMAXM1

CALL ZCOPYXMAXM1, K1, L1, 1, K-XMAXM2, L1, -1 ENDDO

C Multiplikation der Integrale mit -DZZ4*PI*A CALL ZLASCL'G', IDUMMY, IDUMMY, 4D0*PI*A, -DZ, & XMA2M2, XMAX, K, XMA2M2, INFO IF INFO .NE. 0 THEN PRINT *, 'FEHLER IN ZLASCL, INFO=', INFO ENDIF

PRINT 1, 'CPU-Zeit in Sekunden: ', SECOND-STIME PRINT *

PRINT *, 'Belegung der Matrix fuer die FFT'

C Belegung der Matrix fuer die FFT DO L1=0, XMA2M1 DO L2=0, YMA2M1

ENDDO ENDDO PRINT 1, 'CPU-Zeit in Sekunden: ', SECOND-STIME PRINT *

PRINT *, 'Durchfuehrung der FFT'

C Durchfuehrung der FFT

ANHANG D. PROGRAMM-QUELLTEXTE

86

CALL ZFFT2DXMAX2, YMAX2, MFT2, .TRUE.

PRINT 1, 'CPU-Zeit in Sekunden: ', SECOND-STIME PRINT *

PRINT *, 'Belegung der adjungierten Matrix fuer die FFT'

C Belegung der adjungierten Matrix fuer die FFT DO L1=0, XMA2M1 DO L2=0, YMA2M1

ENDDO ENDDO PRINT 1, 'CPU-Zeit in Sekunden: ', SECOND-STIME PRINT *

PRINT *, 'Durchfuehrung der FFT'

C Durchfuehrung der FFT CALL ZFFT2DXMAX2, YMAX2, MHFT2, .TRUE.

PRINT 1, 'CPU-Zeit in Sekunden: ', SECOND-STIME PRINT *

PRINT *, 'Verfahren der konjugierten Gradienten fuer IX'

C Verfahren der konjugierten Gradienten fuer IX C Dieser Algorithmus folgt der Matlab-Funktion cgne Stand 29.12.93 aus: C Martin Hanke, Habilitationsschrift, Karlsruhe 1994

C Initialisierung CURNEG=.TRUE.

CALL ZLAZROINX, 1, ZNULL, ZNULL, X, INX

ANHANG D. PROGRAMM-QUELLTEXTE

87

CALL ZCOPYMN, HY, 1, D, 1 CALL ZFFTMVXMAX2, YMAX2, MHFT, & XMAX, YMAX, D, XMAXM1, YMAX, R, WORK CALL ZCOPYINX, R, 1, S, 1 RNORM=DZNRM2INX, R, 1**2

C Iteration DO L=1, N

& & PHIXL=SQRTRHO*DNORMXL C Ende des Algorithmus nach Martin Hanke

C Berechnung der Kruemmung der L-Kurve IF LCURVE .OR. DIAG THEN

& ENDIF C Suchen der Stelle mit der staerksten Kruemmung IF LCURVE THEN

ANHANG D. PROGRAMM-QUELLTEXTE

88

IF CURVX2 .GT. 0D0 CURNEG=.FALSE. ELSEIF L .GE. 5 THEN C Aktualisierung des Ergebnisses,

C falls die Kruemmung groesser ist als die bisherige maximale Kruemmung C und mindestens einer der benachbarten Kruemmungswerte negativ ist. C Die negativen Kruemmungswerte zu Beginn werden uebersprungen. IF CURNEG THEN

& & ENDIF C Zwischenspeichern der letzten beiden Iterationsschritte

ENDIF ENDDO

C Speichern des Ergebnisses des letzten Iterationsschrittes, C falls die Auswahl nicht mit Hilfe der L-Kurve erfolgt IF .NOT. LCURVE THEN CALL ZCOPYINX, X, 1, IX, 1 ENDIF

PRINT 2, 'Anzahl der durchgefuehrten Iterationsschritte: ' , N

IF LCURVE THEN

PRINT 2, 'maximale Kruemmung der L-Kurve bei Schritt ', ICUMAX ENDIF

PRINT 1, 'CPU-Zeit in Sekunden: ', SECOND-STIME PRINT *

PRINT *, 'Belegung der Matrix fuer die FFT'

C Belegung der Matrix fuer die FFT DO L1=0, XMA2M1 DO L2=0, YMA2M1

ANHANG D. PROGRAMM-QUELLTEXTE

89

ELSE

ENDDO ENDDO PRINT 1, 'CPU-Zeit in Sekunden: ', SECOND-STIME PRINT *

PRINT *, 'Durchfuehrung der FFT'

C Durchfuehrung der FFT CALL ZFFT2DXMAX2, YMAX2, MFT2, .TRUE.

PRINT 1, 'CPU-Zeit in Sekunden: ', SECOND-STIME PRINT *

PRINT *, 'Belegung der adjungierten Matrix fuer die FFT'

C Belegung der adjungierten Matrix fuer die FFT DO L1=0, XMA2M1 DO L2=0, YMA2M1 ENDDO

ANHANG D. PROGRAMM-QUELLTEXTE

90

ENDDO

PRINT 1, 'CPU-Zeit in Sekunden: ', SECOND-STIME PRINT *

PRINT *, 'Durchfuehrung der FFT'

C Durchfuehrung der FFT CALL ZFFT2DXMAX2, YMAX2, MHFT2, .TRUE.

PRINT 1, 'CPU-Zeit in Sekunden: ', SECOND-STIME PRINT *

PRINT *, 'Verfahren der konjugierten Gradienten fuer IY'

C Verfahren der konjugierten Gradienten fuer IY C Dieser Algorithmus folgt der Matlab-Funktion cgne Stand 29.12.93 aus: C Martin Hanke, Habilitationsschrift, Karlsruhe 1994

C Initialisierung CURNEG=.TRUE.

CALL ZLAZROINY, 1, ZNULL, ZNULL, X, INY CALL ZCOPYMN, HX, 1, D, 1 CALL ZFFTMVXMAX2, YMAX2, MHFT, & XMAX, YMAX, D, XMAX, YMAXM1, R, WORK CALL ZCOPYINY, R, 1, S, 1 RNORM=DZNRM2INY, R, 1**2

C Iteration DO L=1, N

& &

ANHANG D. PROGRAMM-QUELLTEXTE

91

PHIYL=SQRTRHO*DNORMYL C Ende des Algorithmus nach Martin Hanke

C Berechnung der Kruemmung der L-Kurve IF LCURVE .OR. DIAG THEN

& ENDIF C Suchen der Stelle mit der staerksten Kruemmung IF LCURVE THEN

IF CURVY2 .GT. 0D0 CURNEG=.FALSE. ELSEIF L .GE. 5 THEN C Aktualisierung des Ergebnisses,

C falls die Kruemmung groesser ist als die bisherige maximale Kruemmung C und mindestens einer der benachbarten Kruemmungswerte negativ ist. C Die negativen Kruemmungswerte zu Beginn werden uebersprungen. IF CURNEG THEN

& & ENDIF C Zwischenspeichern der letzten beiden Iterationsschritte

ENDIF ENDDO

C Speichern des Ergebnisses des letzten Iterationsschrittes, C falls die Auswahl nicht mit Hilfe der L-Kurve erfolgt IF .NOT. LCURVE THEN CALL ZCOPYINY, X, 1, IY, 1

ANHANG D. PROGRAMM-QUELLTEXTE

92

ENDIF

PRINT 2, 'Anzahl der durchgefuehrten Iterationsschritte: ' , N

IF LCURVE THEN

PRINT 2, 'maximale Kruemmung der L-Kurve bei Schritt ', ICUMAX ENDIF

PRINT 1, 'CPU-Zeit in Sekunden: ', SECOND-STIME PRINT *

PRINT *, 'Abspeichern des Ergebnisses'

C Abspeichern des Ergebnisses

C Abspeichern der Loesung fuer die Stroeme OPENUNIT=1, FILE=DATEI1:INDEXDATEI,' '-1'.sol', & FORM='FORMATTED' DO L=1, INX WRITE 1, * IXL ENDDO DO L=1, INY WRITE 1, * IYL ENDDO CLOSE 1

C Abspeichern der Diagnosedatei, falls gewuenscht IF DIAG THEN

& & & & & & & & & & & & & CLOSE 1 ENDIF

PRINT 1, 'CPU-Zeit in Sekunden: ', SECOND-STIME PRINT *

ANHANG D. PROGRAMM-QUELLTEXTE

93

PRINT *, 'Programmende'

END

DOUBLE PRECISION FUNCTION FRES C FUNCTION: Realteil des Integranden

C Massnahme gegen Tippfehler bei Variablen C alle Variablen muessen explizit deklariert werden IMPLICIT NONE

C Variable, ueber die integriert wird DOUBLE PRECISION S

C weitere Funktions-Parameter die Funktion darf nur ein Argument haben DOUBLE PRECISION DX, DY, DZ, BETA COMMON FPARAMM DX, DY, DZ, BETA

C lokale Hilfsvariablen DOUBLE PRECISION DIST, ARG

C Berechnung des Funktionswerts DIST=DX-S**2D0 + DY**2D0 + DZ**2D0 ARG=BETA*SQRTDIST FRE=COSARGSQRTDIST + SINARG*BETADIST

END

DOUBLE PRECISION FUNCTION FIMS C FUNCTION: Imaginaerteil des Integranden

C Massnahme gegen Tippfehler bei Variablen C alle Variablen muessen explizit deklariert werden IMPLICIT NONE

C Variable, ueber die integriert wird DOUBLE PRECISION S

C weitere Funktions-Parameter die Funktion darf nur ein Argument haben DOUBLE PRECISION DX, DY, DZ, BETA COMMON FPARAMM DX, DY, DZ, BETA

C lokale Hilfsvariablen DOUBLE PRECISION DIST, ARG

C Berechnung des Funktionswerts DIST=DX-S**2D0 + DY**2D0 + DZ**2D0 ARG=BETA*SQRTDIST FIM=COSARG*BETA - SINARGSQRTDISTDIST

ANHANG D. PROGRAMM-QUELLTEXTE

94

END

SUBROUTINE ZFFT2DM, N, A, FORWD C Zweidimensionale FFT unskaliert C Bei einer Hin- und Ruecktransformation C wird die Matrix mit M*N multipliziert. C C Variablen: C M Zahl der Zeilen der Matrix A

C N C A C C F O R W D Hintransformation falls .TRUE., C Ruecktransformation falls .FALSE. C C Anmerkung:

C Anstelle einer Variablen DOUBLE PRECISION A2, M, N kann auch eine Variable C DOUBLE COMPLEX ZAM, N uebergeben werden. Die richtige Zuordnung der C Werte geschieht dabei automatisch. Das ist zwar nicht voellig konform C mit dem FORTRAN Standard, funktioniert aber bei allen ueblichen Compilern. C

C Mit dem Standard konform ist folgende Vorgehensweise: C Man definiert sowohl A als auch ZA und setzt sie mittels der Anweisung C EQUIVALENCE A, ZA miteinander gleich. An das Unterprogramm wird dann C statt ZA die Variable A uebergeben. Diese Vorgehensweise ist leider nicht C moeglich, falls man ZFFT2D aus einer SUBROUTINE heraus aufruft, da in C dieser keine neuen Felder mit variablen Grenzen definiert werden koennen, C weil FORTRAN keine dynamische Variablendeklaration zulaesst.

C Massnahme gegen Tippfehler bei Variablen C alle Variablen muessen explizit deklariert werden IMPLICIT NONE

C Deklaration der Parameter INTEGER M, N LOGICAL FORWD DOUBLE PRECISION A2, M, N

C verwendetes Unterprogramm C DOUBLE PRECISION Version von FFT aus netlibbgo, C umgewandelt mit s2d aus netlibbfortran EXTERNAL DFFT

C lokale Hilfsvariable INTEGER MN

MN=M*N

C Durchfuehrung der zweidimensionalen FFT IF FORWD THEN

CALL DFFTA1, 1, 1, A2, 1, 1, MN, M, M, -2 CALL DFFTA1, 1, 1, A2, 1, 1, MN, N, MN, -2

ANHANG D. PROGRAMM-QUELLTEXTE

95

ELSE

CALL DFFTA1, 1, 1, A2, 1, 1, MN, M, M, 2 CALL DFFTA1, 1, 1, A2, 1, 1, MN, N, MN, 2 ENDIF

END

SUBROUTINE ZFFTMVMA, NA, AFFT, MX, NX, X, MY, NY, Y, WORK C Multiplikation einer doppelten Toeplitz-Matrix mit einem Vektor, C von der Matrix wird die zweidimentionale Fouriertransformierte C der zugehoerigen Verschiebungsmatrix benoetigt. C Die Vektoren X und Y werden als zweidimensionale Felder uebergeben, C die der raeumlichen Anordnung der Elemente entsprechen. C Die Routine berechnet Y=A*X. C C Variablen: C M A

C N A C A F F T C M X C N X C X C M Y C N Y C Y C W O R K C C C Anmerkung: C Die Felder X und Y koennen dem Unterprogramm auch als eindimensionale C Felder uebergeben werden. Dies ist konform mit dem FORTRAN Standard. C Die Behandlung ist die gleiche, als wenn das eindimensionale Feld C zuvor einem zweidimensionalen mit den entsprechenden Dimensionen C mittels EQUIVALENCE gleichgesetzt wuerde und dann das zweidimensionale C Feld uebergeben wuerde.

C Massnahme gegen Tippfehler bei Variablen C alle Variablen muessen explizit deklariert werden IMPLICIT NONE

C Deklaration der Parameter INTEGER MA, NA, MX, NX, MY, NY DOUBLE COMPLEX AFFTMA, NA DOUBLE COMPLEX XMX, NX DOUBLE COMPLEX YMY, NY DOUBLE COMPLEX WORKMA, NA

C verwendete Unterprogramme ZFFT2D im Programmtext, andere aus LAPACK EXTERNAL ZLAZRO EXTERNAL ZLACPY EXTERNAL ZFFT2D EXTERNAL ZDRSCL

ANHANG D. PROGRAMM-QUELLTEXTE

96

C Konstante DOUBLE COMPLEX ZNULL PARAMETER ZNULL=0D0, 0D0

C lokale Hilfsvariablen INTEGER L1, L2 DOUBLE PRECISION SCAL

SCAL=MA*NA

C zweidimensionale Fouriertransformation des Vektors X CALL ZLAZROMA, NA, ZNULL, ZNULL, WORK, MA CALL ZLACPY'A', MX, NX, X, MX, WORK, MA CALL ZFFT2DMA, NA, WORK, .TRUE.

C Multiplikation der Komponenten der Fouriertransformierten DO L1=1, MA DO L2=1, NA WORKL1, L2=AFFTL1, L2*WORKL1, L2 ENDDO ENDDO

C Ruecktransformation und Skalierung des Vektors Y C Skalierung notwendig, da ZFFT2D die unskalierte FFT berechnet CALL ZFFT2DMA, NA, WORK, .FALSE. CALL ZLACPY'A', MY, NY, WORK, MA, Y, MY CALL ZDRSCLMY*NY, SCAL, Y, 1

END

ANHANG D. PROGRAMM-QUELLTEXTE

97

D.2 Programm dat2inp.f

PROGRAM DAT2INP

************************************************************************ * *

* Programm: DAT2INP * Autor: * Betreuer: Georg Faessler * Datum: * * * Programm zur Umwandlung der Messdaten *

* in den Vektor, der von dem Programm CGFFT benoetigt wird * * *

************************************************************************

C In diesem Programm werden folgende Befehle verwendet, die nicht dem C FORTRAN77-Standard entsprechen, jedoch von allen ueblichen Compilern C verstanden werden: C

C IMPLICIT NONE fuer die Abschaltung der impliziten Variablendeklaration C DOUBLE COMPLEX fuer komplexe Variablen doppelter Genauigkeit C DO ohne Zeilennummer ... ENDDO

C FORMAT..., $ fuer die Unterdrueckung des Zeilenvorschubs

C Massnahme gegen Tippfehler bei Variablen C alle Variablen muessen explizit deklariert werden IMPLICIT NONE

C Bereich, aus dem Messwerte vorhanden sind INTEGER XMAX, YMAX PARAMETER XMAX=125, YMAX=75

C Konstanten DOUBLE COMPLEX J PARAMETER J=0D0, 1D0 DOUBLE PRECISION PI PARAMETER PI=3.1415926535897932384626433D0 DOUBLE PRECISION GRAD PARAMETER GRAD=PII180D0

C Variablen

DOUBLE COMPLEX HXXMAX, YMAX, HYXMAX, YMAX INTEGER L1, L2 CHARACTER*20 DATEI DOUBLE PRECISION BETRAG, PHASE

C Eingabe des Dateinamens 1 FORMAT T2, A, $ PRINT *

PRINT *, 'Programm DAT2INP - Umwandlung der Daten fuer CGFFT' PRINT *

PRINT 1, 'Dateiname ohne Erweiterung: ' READ *, 'A' DATEI

ANHANG D. PROGRAMM-QUELLTEXTE

98

C Einlesen der Daten und Umwandlung in komplexe Zahlen OPEN 9, FILE=''ihffuserrwiedmannndataa' &

DO L2=1, YMAX ENDDO ENDDO CLOSE 9 OPEN 9, FILE=''ihffuserrwiedmannndataa' &

DO L2=1, YMAX ENDDO ENDDO CLOSE 9 C Abspeichern des Ergebnisses OPEN 10, FILE=DATEI1:INDEXDATEI,' '-1'.inp', & FORM='UNFORMATTED' WRITE 10 HX, HY CLOSE 10

END

ANHANG D. PROGRAMM-QUELLTEXTE

99

Matlab-Routine cgne

D.3

function xk,ek,dk,hk,kstar = cgneA,b,x,k,tol xk,ek,dk,hk,kstarr = cgneA,b,x,k,tol

Performs k steps of cg applied to the normal equations system A'Ax=A'b. A' always denotes the conjugate transpose of A

The routine returns the errors, the residuals and the heuristic comparison sequence in vectors ek, dk, and hk. xk contains the kth iterate.

If tol is specified then kstar is the iteration index where the dk-sequence drops below tol

Reference: M. Hanke, Habilitationsschrift, 1994.

Martin Hanke, Karlsruhe, 29912293.

m,nn = sizeA ek = zerosk,1 dk = ekk hk = dkk if nargin 5 t o l = 0 else kstar = 00 endd

Prepare for CG iteration. xk = zerosn,1 r = b Ar = A'*rr d = A r normAr = normAr^22 beta = 0

Iterate. for j=1:k Ad = A*dd alpha = normArrnormAd^2 if j==1 d_abl = alphaa abl = d_abll alpha_alt = alphaa else

d_abl = alpha + alpha*beta*d_ablalpha_alt abl = abl + d_abll alpha_alt = alphaa endd xk = xk + alpha*dd r = r - a l p h a * A d Ar = A'*rr normAr_new = normAr^22 beta = normAr_newwnormArr

ANHANG D. PROGRAMM-QUELLTEXTE

100

normAr = normAr_neww d = Ar + beta*dd ekj = normxk-x dkj = normr hkj = sqrtabl * dkj if dkj = tol kstar = jj tol = 00 endd end

Literaturverzeichnis

11 M. Hanke, P. C. Hansen, Regularization methods for large-scale problems, Surv. Math. Ind. 3 1993, S. 253-315

22 M. Hanke, J. G. Nagy und R. J. Plemmons, Preconditioned Iterative Regularization for Ill-Posed P r oblems, in Numerical Linear Algebra and Scientic Computing, L. Reichel, A. Ruttan und R. S. Varga Hrsg., de Gruyter, Berlin, 1993

33 M. Hanke, T. Raus, A G e n e r al Heuristic for Choosing the Regularization Parameter in Ill-Posed P r oblems, eingereicht bei SIAM J. Sci. Comput.

44 M. Hanke, M. Hochbruck, A Chebyshev-Like Semiiteration for Inconsistent Linear Systems, Elec. Trans. Numer. Anal. 1 1993, S. 89-103

55 M. Hanke, Habilitationsschrift, Institut f ur Praktische Mathematik, Universit at Karlsruhe, 1994

66 P. C . H a n s e n , Numerical Tools for Analysis and Solution of Fredholm Integral Equations of the First Kind, Inverse Problems 8 1992, S. 849-872

77 P. C . H a n s e n ,Analysis of Discrete Ill-Posed P r oblems by Means of the L-Curve, SIAM Rev. 34 1992, S. 561-580

88 P. C. Hansen, Truncated Singular Value Decomposition Solutions to Discrete Ill-Posed P r oblems with Ill-Determined Numerical Rank, SIAM J. Sci. Stat. Comput. 11 1990, S. 503-518

99 P. C. Hansen, The Discrete Picard Condition for Discrete Ill-Posed P r oblems, BIT 30 1990, S. 658-672

100 P. C. Hansen, The Truncated SVD as a Method for Regularization, BIT 27 1987, S. 534-553

111 J. G. Nagy, R. J. Plemmons, T. C. Torgersen, A N e w P r econditioner for Iterative Deconvolution, Vorabdruck, 5. April 1994

122 R. H. Chan, J. G. Nagy, R. J. Plemmons, Displacement Preconditioner for Toeplitz Least Squares Iterations, Elec. Trans. Numer. Anal 2 1994, S. 44-56

101

LITERATURVERZEICHNIS 102

133 A. K. Louis, Inverse und schlecht gestellte Probleme, Teubner, Stuttgart,

1989

144 G. M. Wing, A Primer on Integral Equations of the First Kind: the Problem of Deconvolution and Unfolding, SIAM, Philadelphia, 1991

155 J. R. Shewchuk, An Introduction to the Conjugate Gradient Method Without the Agonizing Pain, Report CMU-CS-94-125, School of Computer Science, Carnegie Mellon University, Pittsburgh, 7. M arz 1994

166 J. A. Scales, Iterative Methods for Large, Sparse, Inverse Calculations, Version 1.0, April 1993

177 R. Zurm uhl, S. Falk, Matrizen und ihre A nwendungen: f ur Ingenieure, Physiker und Angewandte Mathematiker, 6. Auage, Springer, Berlin, Heidelberg, 1992

188 R. F. Harrington, Field Computation by Moment Methods, Macmillan Company, N e w Y ork, 1968

199 F. Landstorfer, Hochfrequenztechnik I-III und Mikrowellentechnik, Vorlesungen, Universit at Stuttgart

200 G. Lehner, Elektromagnetische Feldtheorie, Springer, Berlin, Heidelberg,

1990

211 A. Sautter, Feldverteilung auf Platinen, Erste Semesterarbeit, Institut f ur Hochfrequenztechnik, Universit at Stuttgart, Mai 1994

222 P. S c haich, Momentenmethode zur numerischen Berechnung von Leiterplattenstr omen, Erste Semesterarbeit, Institut f ur Hochfrequenztechnik, Universit at Stuttgart, Juli 1993

233 M. Wintermantel, Faltung und Korrelation, Versuch Nr. 4, Fachpraktikum der digitalen Signalverarbeitung, Institut f ur Netzwerk- und Systemtheorie, Universit at Stuttgart

244 U. Jakobus, Bedienungsanleitung zum Programm FEKO, Institut f ur Hochfrequenztechnik, Universit at Stuttgart

255 C. L. Lawson, R. J. Hanson, D. R. Kincaid und F. T. Krogh, Basic Linear Algebra Subprograms for Fortran Usage, ACM Trans. Math. Soft. 5 1979, S. 308-323

266 J. Dongarra, J. Du Croz, S. Hammarling und R. Hanson, An Extended Set of Fortran Basic Linear Algebra Subprograms, ACM Trans. Math. Soft. 14 1988, S. 1-17

277 J. Dongarra, J. Du Croz, S. Hammarling und R. Hanson, An Extended S e t of Fortran Basic Linear Algebra Subprograms: Model Implementation and Test Programs, ACM Trans. Math. Soft. 14 1988, S. 18-32

LITERATURVERZEICHNIS 103

288 J. Dongarra, J. Du Croz, I. Du und S. Hammarling, A S e t o f L evel 3 Basic Linear Algebra Subprograms, ACM Trans. Math. Soft. 16 1990, S. 1-17

299 J. Dongarra, J. Du Croz, I. Du und S. Hammarling, A S e t o f L evel 3 Basic Linear Algebra Subprograms: Model Implementation and Test Programs, ACM Trans. Math. Soft. 16 1990, S. 18-28

300 E. Anderson, Z. Bai, C. Bischof, J. Demmel, J. Dongarra,

J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, S. Ostrouchov und D. Sorensen, LAPACK Users' Guide, SIAM, Philadelphia, 1992

311 E. Anderson, J. Dongarra, S. Ostrouchov, LAPACK Working Note 41, Installation Guide for LAPACK, Department of Computer Science, University o f T ennessee, Knoxville, Version 1.1, 31. M arz 1993

322 P. F ava t i e t a l . , A CM Trans. Math. Soft. 17 1991, S. 218-232

333 H. Kopka, L a T E X, eine Einf uhrung, 4. Auage, Addison-Wesley, Bonn, M unchen, Paris, 1992

344 H. Kopka, L a T E X-Erweiterungsm oglichkeiten mit einer Einf uhrung in Me- tafont, 2. Auage, Addison-Wesley, Bonn, M unchen, Paris, 1991