1 Worum geht es?

Das wissenschaftliche Rechnen (auch: Scientific Computing) als Fachgebiet umfasst alle Aspekte des Einsatzes von Computern zur Klärung von Fragestellungen aus Naturwissenschaften, Technik und Wirtschaft. Das Spektrum reicht von der Modellierung realer Vorgänge über die Implementierung am Computer bis zur Visualisierung der erzielten Resultate.

Der erfolgreiche Einsatz des wissenschaftlichen Rechnens erfordern die Zusammenarbeit von Anwendern aus dem jeweiligen Wissenschaftsbereich, Mathematikern und Informatikern. Alle Beteiligten sollten neben Detailwissen in ihrem Gebiet auch einen Überblick über den Gesamtprozess haben. Nur so sind Missverständnisse in der Kommunikation und Fehlinterpretationen der jeweiligen Anforderungen zu vermeiden.

1.1Zu lösende Teilprobleme¶

Wesentliche Schritte beim wissenschaftlichen Rechnen sind:

• Modellierung (Anwender + Mathematiker)

  • Welche mathematischen Modelle für den betrachteten Vorgang gibt es?

  • Welches davon nutzen? Einfaches oder komplexes Modell?

  • Eigenes/neues Modell entwickeln?

  • Wie gut ist das Modell (Modellierungsfehler)?

• Diskretisierung (Mathematiker + Informatiker)

  • Wie kann das Modell in computertaugliche Form gebracht werden?

  • Ggf. ist die Auswahl von Algorithmen im nächsten Schritt schon hier einzubeziehen!

  • Welcher Fehler entsteht dabei (Diskretisierungsfehler)?

• Implementierung (Informatiker + Mathematiker)

  • Welche Algorithmen wählen?

  • Gibt es geeignete Bibliotheken?

  • Unnötige Ungenauigkeiten in Berechnungen verweiden (Approximationsfehler, Rundungsfehler, vgl. Veranstaltung zur numerischen Mathematik)!

  • In welcher Form sollen die Ergebnisse geliefert werden (Tabellen, Visualisierungen,...)?

• Einsatz der Software (Anwender)

  • Interpretation der erzielten Ergebnisse. Sind die plausibel?

Am Ende sollte ein Software-Produkt vorliegen, welches Ergebnisse mit bekannter Genauigkeit liefert. Ein Schwerpunkt des wissenschaftlichen Rechnens wird also die Fehleranalyse sein. Ohne solide Fehleranalyse sind die erzielten Ergebnisse nicht vertrauenswürdig und damit im Wesentlichen unbrauchbar! Manchmal kann man die Ergebnisse durch reale Messungen validieren, oft aber nicht oder nicht in hinreichendem Umfang.

Das wissenschaftliche Rechnen ist also eng verbunden mit der numerischen Mathematik. Eine genaue Abgrenzung ist nicht möglich. Üblicherweise zählt man das Entwickeln und Untersuchen von Algorithmen für klar abgegrenzbare, unabhängig vom konkreten Anwendungsfall immer wiederkehrende Teilaufgaben zur numerischen Mathematik. Komplexere Algorithmen und die Koordinierung der Einzelschritte zu einem Gesamtprozess werden hingegen dem wissenschaftlichen Rechnen zugeordnet.

1.2Beispiele¶

Die folgenden Beispiel dienen zur Verdeutlichung der im wissenschaftlichen Rechnen zu lösenden Probleme und Teilprobleme. Später werden wir einige der Beispiele detaillierter untersuchen.

1.2.1Brachistochrone¶

Anwendungsproblem: Eine Kugel soll sich aus “eigener” Kraft (Schwerkraft) entlang einer Bahn von einem höher gelegenen Startpunkt zu einem niedriger gelegenen Zielpunkt bewegen. Wie ist die Bahn zu Formen, damit die Kugel in möglichst kurzer Zeit den Zielpunkt erreicht? Diese optimale Bahn heißt Brachistochrone (Wikipedia).

Modell: Sind $(x_1,y_1)$ und $(x_2,y_2)$ Start- und Zielpunkt, so kann die Bahn als Funktion $h:[x_1,x_2]\to\bbR$ beschrieben werden, insbesondere gilt also $h(x_1)=y_1$ und $h(x_2)=y_2$. Man kann zeigen, dass die benötigte Zeit $T(h)$ für das Durchlaufen der Bahn $h$ mittels

berechnet werden kann. Dabei ist $g$ die wirkende Fallbeschleunigung. Die optimale Bahn ist somit als Lösung des Minimierungsproblems

gegeben, wobei $H$ die Menge aller stetig differenzierbaren Funktionen mit $h(x_1)=y_1$ und $h(x_2)=y_2$ sein soll.

An diesem recht einfachen Problem und seiner Modellierung treten schon verschiedene Aspekte des wissenschaftlichen Rechnens zutage:

• Obwohl das Anwendungsproblem im dreidimensionalen Raum formuliert ist, arbeitet das Modell im zweidimensionalen Raum. Dass diese Vereinfachung zulässig ist, also keinen (wie hier) oder nur einen vernachlässigbar kleinen Fehler verursacht, muss man bei der Modellierung prüfen!

• Die Formel für $T(h)$ nimmt keine Rücksicht auf die Tatsache, dass die Kugel auf der Bahn rollen und nicht gleiten wird. Man kann sich auch hier überlegen, dass diese Vereinfachung keinen Einfluss auf die Lösung hat.

• Die Vernachlässigung der Reibung zwischen Kugel und Bahn sowie zwischen Kugel und Luft wird hingegen einen Fehler verursachen. Dieser sollte (!) vernachlässigbar sein, wenn die Kugel schwer und sowohl Kugel als auch Bahn sehr glatt sind. Diese Aussage ist keine exakte Fehlerbetrachtung, sondern eine bei der Modellierung leider manchmal nötige Abschätzung nach “Augenmaß”.

• Wir haben als Wertebereich für $h$ alle reellen Zahlen zugelassen, obwohl man in Versuchung kommen könnte, nur das Intervalle $[y_2,y_1]$ zuzulassen. Es wird sich allerdings herausstellen, dass für die optimale Bahn auch $h(x)<y_2$ für manche $x$ gelten kann. Also: Vereinfachungen kritisch auf Zulässigkeit prüfen.

• Im Integrand wird durch $h(x)$ geteilt, der Integrand hat also eine Polstelle. Da das Integral minimiert wird, wird die Polstelle der Lösung aber so beschaffen sein, dass ein endlicher Wert für das Integral entsteht. Modelle auf mathematische Korrektheit prüfen!

• Das Modell wird als Optimierungsproblem über einer Menge von Funktionen formuliert. Diese Art mathematischer Probleme taucht in den Grundveranstaltungen zur Mathematik üblicherweise nicht auf. Alle praktisch relevanten und interessanten Probleme erfordern zusätzliche Mathekenntnisse!

Diskretisierung: Die Menge aller stetig differenzierbaren Funktionen ist nicht mit dem Computer darstellbar. Somit müssen wir uns auf eine hinreichend repräsentative, aber endlic hdimensionale Menge von Funktionen als Suchraum bei der Optimierung beschränken. Kandidaten sind Polynome oder stückweise Polynome, jedoch auch eine ganze Reihe anderer Ansätze. Man kann zeigen, dass die optimale Bahn $h(x_1)=-\infty$ erfüllt. Mit Polynomen wird man diese Bedingung nie erfüllen können. Verwendet man stückweise Polynome, so sollten die Teilintervalle in der Nähe von $x_1$ sehr klein sein, damit der starke Abfall von $h$ dort gut angenähert werden kann. Obwohl die gesuchte Funktion differenzierbar ist, können auch nicht differenzierbare Ansätze wie stückweise lineare Funktionen zum Ziel führen (keep it simple). Dabei ist die Frage zu klären, was eigentlich das Ziel ist. Geht es nur um die Visualisierung der optimalen Bahn? Dann reicht eine stückweise lineare Näherungslösung aus. Oder soll die optimale Bahn in Folgeprozessen verwendet werden, die auf Ableitungen beruhen (z.B. tiefsten Punkt der Bahn analytisch berechnen)? Dann muss sollte auch die Näherungslösung differenzierbar sein.

Implementierung: Welchen Algorithmus zur Lösung des diskretisierten Optimierungsproblems wählen? Gibt es anwendbare Standardalgorithmen aus Bibliotheken oder muss ein auf die konkrete Aufgabe maßgeschneiderter Algorithmus entwickelt werden? In welcher Form das das Ergebnis geliefert werden? Nur die optimale Bahn $h$ oder beispielsweise eine Animation der Bewegung der Kugel entlang der Bahn mit korrekter Geschwindigkeit bzw. Beschleunigung? Beachte: Die gewählte Diskretisierung hat Einfluss auf die Art des entstehenden endlichdimensionalen Optimierungsproblems und damit auf die Auswahl der Algorithmen!

Randbemerkung: Das Brachistochrone-Problem kann man analytisch lösen, sodass eine numerische Lösung eigentlich nicht nötig ist. Dennoch sieht man an diesem einfachen Beispiel gut, wie vielfältig die zu klärenden Fragen beim wissenschaftlichen Rechnen selbst für einfache Problem sind.

1.2.2Laplace-DLTS¶

Anwendungsproblem: Neben Silizium werden zunehmend alternative Materialien für Solarzellen untersucht und teils auch schon eingesetzt. Dazu gehört das recht preiswerte Material Perowskit. Durch noch nicht hinreichend präzise kontrollierbare Produktionsprozesse kann es zu Materialfehlern kommen, welche den Wirkungsgrad der Solarzellen verringern. Deshalb ist man an experimentellen Verfahren interessiert, die Einblick in die Art der Materialfehler geben, um diese anschließend die Produtionsprozesse gezielt zu verbessern. Ein solches Verfahren ist Laplace-DLTS (Laplace deep-level transient spectroscopy).

Die für Laplace-DLTS verwendeten experimentellen Messergebnisse sind ohne weitere algorithmische Auswertung nicht verwertbar. Vereinfacht dargestellt: Sogenannte Elektronenfallen im Material (Materialfehler) werden mit Ladungung gefüllt und anschließend beobachtet man das relativ langsame Abfließen der Ladungen aus den Fallen. Dieser Vorgang wird bei verschiedenen Temperaturen wiederholt. Es entsteht eine Folge von sogenannten Transienten.

Jede gemessene Transiente entsteht durch Überlagerung der Transienten mehrerer Elektronenfallen im Material. Da das Abklingverhalten einer Einzeltransiente bekannt ist (Exponentialfunktion), kann die Überlagerung rückgängig gemacht werden. Man erhält als Modell für die Entstehung der Transienten die Laplace-Transformation

Dies ist eine lineare Integralgleichung. Dabei ist $C$ die in Abhängigkeit der Zeit $t$ gemessene Gesamttransiente zu einer festen Temperatur und $g$ ist die Verteilung der Abklinggeschwindigkeiten in der Menge aller Einzeltransienten. Typischerweise wird $g$ aus wenigen sehr scharfen Peaks bestehen. Aus Position, Höhe und Breite der Peaks kann auf Art und Anzahl der Materialfehler geschlossen werden. Trägt man die Peak-Positionen bei verschiedenen Temperturen in ein Diagramm ein, so erhält man einen sogenannten Arrhenius-Plot. Dieser zeigt typischerweise abfallende Gerade, aus deren Lage und Neigung auf weitere Materialeigenschaften geschlossen werden kann.

Modellierung: Grundfrage ist hier, ob mehrere Einzelmodelle verwendet werden sollen, also

• Invertieren der Laplace-Transformation zum Erhalt der Abklingratenfunktion $g$,

• Ermitteln der Peak-Positionen in $g$ bei verschiedenen Temperaturen,

• Finden von Geraden im Arrhenius-Plot,

oder ob das Problem “Ende-zu-Ende” formuliert werden soll: Finde Lage und Anstieg der Geraden im Arrhenius-Plot aus den $C$-Messungen bei verschiedenen Temperaturen. Die erste Variante ist einfacher zu modellieren, da jedes Teilmodell eine klar abgegrenzte Aufgabe löst. Die zweite Variante vermeidet eventuelle Ungenauigkeiten in den Einzelschritten (Peak-Position bei breiten Peaks,...), ist aber in der mathematischen Formulierung sehr anspruchsvoll. Die zweite Variante ist ein typischer Kandidat für mittels maschinellem Lernen erstellte Modelle, sofern genügend Trainingsdaten beschafft werden können (siehe später). Verfolgen hier nur die erste Variante weiter.

Diskretisierung: Klar ist, dass sowohl $g$ als auch $C$ durch Vektoren zu ersetzen sind, z.B. als Werte von stückweise linearen Funktionen an den Stützstellen. Da $C$ direkt aus Messungen entsteht, sollte sich die Diskretisierung hier an den messtechnischen Gegebenheiten orientieren. Im konkreten Anwendungsfall besteht eine $C$-Messung aus ca. 300000 Messwerten an äquidistanten Stützstellen. Wird diese hohe Auflösung beibeihalten, entsteht ein extrem hoher Rechenaufwand! Gleichmäßiges Ausdünnen der Messpunkte führt allerdings zu Datenverlust. Ob dieser hinnehmbar ist, muss mit den Anwendern geklärt werden. Es stellt sich heraus, dass diese hohe Auflösung am Anfang des Messzeitraums notwendig ist (sehr schnell abfallende Einzeltransienten), im weiteren Zeitverlauf aber zunehmend ausgedünnt werden kann (langsam abfallende Einzeltransienten). Für $C$ ist also eine geeignete nichtäquidistante Diskretisierung zu wählen. Für die Abklingratenfunktion $g$ sieht es ähnlich aus. Werte Nahe der Null müss präzise darstellbar sein; größere Werte können mit geringerer Genauigkeit wiedergegeben werden. Alternative: Da aufgrund physikalischer Gegebenheit $g$ nur aus wenigen Peaks besteht, kann man $g$ auch als Summe von Einzelpeaks darstellen. Der Zusammenhang zwischen den entsprechenden Parametern von $g$ (Peak-Positionen, PeakHöhen) und der Funktion $C$ ist dann allerdings nichtlinear, was das algorithmische Invertieren der Laplace-Transformation deutlich erschwert.

Implementierung: Das numerische Invertieren der Laplace-Transformation ist sehr schwierig, obwohl es sich (nach Diskretisierung) um ein lineares Gleichungssystem handelt. Problem ist die extrem hohe Kondition dieses Gleichungssystems, welche um so höher ist, je feiner (!) diskretisiert wird. Einzige Rettung ist das Einbringen von Zusatzinformationen über die Lösung (z.B. dass diese nur wenige scharfe Peaks hat), was algorithmisch sehr anspruchsvoll werden kann. Siehe Improved evaluation of deep-level transient spectroscopy on perovskite solar cells reveals ionic defect distribution für einen konkreten Algorithmus.

Sehen auch hier typische Aspekte des wissenschaftlichen Rechnens:

• Bei der Modellierung gibt es viele Freiheiten, ober auch Unschärfen. Bereits bei der Modellierung sollte man die weiteren Schritte zur Lösung im Blick haben.

• Diskretisierung unterliegt vielen praktischen Zwängen wie Art und Anzahl verfügbarer Messungen oder der Verfügbaren Computerressourcen.

• Manchmal gibt es keine Standardalgorithmen, obwohl das gewählte Modell seit Jahrzehnten vielfach im Einsatz ist bei ähnlichen Fragestellungen.

1.2.3SD-SPIDER¶

Die kürzesten je von Menschen erzeugten Ereignisse sind mit Femtosekundenlasern erzeugte ultrakurze Laserpulse. Jeder Puls dauert nur wenige Femtosekunden (die Lichtgeschwindigkeit liegt bei 0.3 Mikrometer pro Femtosekunde). Zur Untersuchung solcher Laserpulse können keine einfachen Verfahren eingesetzt werden, da diese praktisch immer auf Abtastung des Signals basieren. Das Signal ist hier aber kürzer als der Abstand zwischen zwei Abtaststellen selbst bei optischen Verfahren.

Eine Möglichkeit zur Charakterisierung von ultrakurzen Laserpulsen ist SD-SPIDER (self-diffraction spectral phase interferometry for direct electric field reconstruction). Kurzfassung: Der zu untersuchende Laserpuls wird auf verschiedene Weisen transformiert und dann mittels eines nichtlinearen optischen Materials (Bariumfluorid) wieder zusammengeführt. Das Ergebnis kann dann durch einen Spektrograf aufgenommen werden. Die ursprüngliche Pulsform kann daraus jedoch nicht direkt abgelesen werden, sondern erfordert relativ komplexe Berechnungen.

Modellierung: Aus physikalischen Zusammenhängen kann man ableiten, dass der Zusammenhang zwischen gemessener Funktion $y:\bbR\to\bbC$ (aus Spektrograf) und gesuchter Funktion $x:[0,1]\to\bbC$ (auf Einheitsintervall transformierte Pulsform) durch die nichtlineare Integralgleichung

gegeben ist. Dabei ist $k$ eine bekannte Funktion, die sogenannte Kernfunktion. Diese enthält Informationen über das eingesetzte nichtlineare Material und weitere Parameter des Experiments. Gleichungen dieser Form heißen Selbstfaltungsgleichungen.

Diskretisierung: Für die Diskretisierung können Standardverfahren eingesetzt werden (z.B. stückweise konstante Funktionen). Allerdings ist zu beachten, dass die Funktionen $x$ und $y$ komplexwertig sind, Funktionswerte also als Paar von reellen Zahlen anzusehen sind mit entsprechenden Anpassungen bei den Rechenoperationen. Ergebnis ist ein nichtlineares Gleichungssystem mit sehr schlechter Kondition (je feiner diskretisiert wird, desto schlechter die Kondition).

Implementierung: Für das zu lösende nichtlineare Gleichungsystem ist kein Algorithmus bekannt, der garantiert die korrekte Lösung oder wenigstens eine hinreichend gute Näherungslösung liefert. Standardverfahren (z.B. Newton-Verfahren) sind im Prinzip anwendbar, aber keine der dafür bekannten Konvergenzaussagen ist im konkreten Fall zutreffend, da das Gleichungssystem die in diesen Aussagen formulierten Voraussetzungen nicht erfüllt. Auf der anderen Seite hat das Gleichungssystem viel Struktur, die sich ggf. zur Entwicklung von Algorithmen nutzen lässt. In der Tat stellt sich heraus, dass man die hier vorhandene “quadratische” Struktur nutzen kann für eine Zerlegung in zwei Teilprobleme. Eins davon ist ein (großes, aber leicht parallelisierbares) lineares Gleichungssystem. Das andere ist ein Eigenwertproblem, welches mit Standardverfahren gelöst werden kann. Für Details siehe Deautoconvolution: A new decomposition approach versus TIGRA and local regularization

Wir sehen:

• Schlecht konditionierte Probleme treten häufiger auf im wissenschaftlichen Rechnen als man erwarten würde. Insbesondere in Physik und Technik sind sie sehr präsent. Entsprechend vorsichtig muss bei der Auswahl von Lösungsverfahren und konkreten Algorithmen vorgehen.

• Für aktuelle Aufgabenstellungen aus der Forschung gibt es meist noch keine Standardverfahren. Bestenfalls kann man bekannt Verfahren neu kombinieren oder anpassen.

• Die Aufgabenverteilung zwischen Mathematik (nichtlineares Problem in einfachere Teilprobleme zerlegen) und Informatik (effiziente Lösbarkeit der Teilprobleme und damit praktische Durchführbarkeit) ist fließend.

1.2.4Lorenz-System¶

Schon in der Anfangszeit der Computer wollte man Wetterphänomene simulieren und zukünftige Entwicklungen in der Atmosphre vorhersagen. Aufgrund technischer Beschränkungen mussten die Fragestellungen mit sehr einfachen Modellen formuliert und beantwortet werden. Eins davon ist das Lorenz-System, welches Aussagen zum vertikalen Wärme- und Feuchtigkeitstransport in der Atmosphäre liefert.

Modellierung: Sind $x$, $y$, $z$ von der Zeit $t$ abhängige Atmosphärenparameter (vertikaler Feuchtetransport pro Zeiteinheit, Temperaturdifferenz zwischen auf- und absteigenden Luftmassen, Abweichung des vertikalen Temperaturprofils von einer linearen Funktion), so ist deren Zusammenhang durch ein System nichtlinearer gewöhnlicher Differentialgleichungen gegeben:

Dabei sind $a$ und $b$ Materialparameter der strömenden Luft und $c$ gibt gewissen geometrische Eigenschaften des betrachteten Luftvolumens an.

Diskretisierung: Für das Lösen von gewöhnlichen Differentialgleichungen existieren zahlreiche Verfahren. Die meisten nutzen zur Diskretisierung stückweise lineare Funktionen.

Implementierung: Es können Standardalgorithmen für gewöhnliche Differentialgleichungen genutzt werden. Allerdings ist bei der Auswahl zu beachten, diese mit der Tatsache umgehen können, dass das Differentialgleichungssystem schlecht konditioniert ist, also kleine Änderungen in $a$, $b$, $c$ zu extrem großen Änderungen in der Lösung führen können (schlechte Kondition). Art und Umfang dieser Abweichungen sehen sehr unsystematisch aus, weshalb man den Lösungen des Lorenz-Systems auch von chaotischem Verhalten spricht. In Kontext des Lorenz-Systems wurde auch erstmals der Begriff Schmetterlingseffekt genutzt.

Bemerkenswerte Aspekte:

• Auch einfache Modelle können sehr komplexes Lösungsverhalten zeigen. Man unterschätze nie die Schwierigkeiten, die auch scheinbar einfache Aufgabenstellungen mit sich bringen können!

1.2.5Wärmeleitung¶

Ein während des Produktionsprozesses durch und durch erhitztes Stahlteil muss vor dem nächsten Verarbeitungsschritt unter eine zulässig Höchsttemperatur abkühlen. Wie lange dauert es bis auch im Inneren des Bauteils in jedem Punkt diese Temperatur unterschritten wird?

Modellierung: Bezeichnen mit $B\subseteq\bbR^3$ das Bauteil als Teilmenge des Raumes und mit $\partial B$ dessen Oberfläche. Sei $u:B\times[0,\infty)\to\bbR$ die zeitabhängige Temperaturverteilung im Bauteil; $u(x,t)$ gibt also die Temperatur zur Zeit $t$ im Punkt $x\in B$ an. Der Abkühlprozess beginne bei $t=0$. Aus physikalischen Überlegungen und Gesetzmäßigkeiten erhält man den Zusammenhang

für alle $x\in B$ und alle $t\geq 0$, wobei $a>0$ ein Materialparameter ist, der im Wesentlichen die Wärmeleitfähigkeit des Materials wiedergibt. Sind $T_{\mathrm{U}}$ und $T_0$ die Umgebungstemperatur und die Temperatur des Bauteils zur Zeit $t=0$, so kommen noch die Anfangsbedingung

und die Randbedingung

zu dieser partiellen Differentialgleichung hinzu.

Diskretisierung: Die Diskretisierung kann mittels Standardverfahren erfolgen. Diese führen auf ein lineares Gleichungssystem. Zu beachten ist, dass die Geometrie des Körpers wesentlichen Einfluss auf die Diskretisierung hat. Das Diskretisieren von partiellen Differentialgleichungen ist wesentlicher Bestandteil dieser Veranstaltung.

Implementierung: Grundsätzlich stehen vielfältige Algorithmen in Software-Bibliotheken zur Verfügung. Je nach Geometrie des Bauteils können aber Anpassungen oder Erweiterungen nötig sein.

Relevante Aspekte:

• Der Großteil der Gesetze der (klassischen) Physik sind in Form von partiellen Differentialgleichungen formuliert. Entsprechend häufig tauchen partielle Differentialgleichungen als Modelle in der Anwendung auf.

• Für das Formulieren von partiellen Differentialgleichungen sind solide Fähigkeiten in der mathematischen Sprache sehr hilfreich. Das analytische Untersuchen (und Lösen, sofern möglich) von partiellen Differentialgleichung erfordert sehr tiefe Mathematikenntnisse.

1.2.6Elektrische Impedanz-Tomografie¶

Möchte man das Innere eines Testobjekts zerstörungsfrei untersuchen, kommen verschiedene Tomografieverfahren in Betracht. Eins davon, welches insbesondere in der medizinischen Bildgebung im Einsatz ist, aber mit leichten Modifikationen auch zur Untersuchung von Bodenverhältnissen und Rohstofflagerstätten dient, ist die elektrische Impedanz-Tomografie. Hier werden an der Oberfläche des Testobjekts Elektroden angebracht und in verschiedenen Kombinationen Spannungen angelegt. Die daraus an den anderen Elektroden resultierenden Ströme können gemessen werden und geben Aufschluss über die (ortsveränderliche) Leitfähigkeit und damit über Materialstrukturen im Inneren des Testobjekts.

Modellierung: Sei $B\subseteq\bbR^2$ das Testobjekt und seien die Elektroden gleichmäßig über dessen Rand $\partial B$ verteilt (ringförmige Anordnung der Elektroden an Testkörper, somit reicht 2D-Betrachtung). Die Leitfähigkeit im Testkörper sei durch eine unbekannte Funktion $a:B\to\bbR$ gegeben. Bezeichnen wir mit $u:B\to\bbR$ das elektrische Potential im Testobjekt, so folgt aus physikalischen Überlegungen und Gesetzmäßigkeiten, dass die partielle Differentialgleichung

erfüllt sein muss. Das Anlegen von Spannungen an verschiedenen Elektrodenkombinationen entspricht der Vorgabe verschiedener Randwerte $f$ für Lösungen $u$ dieser Gleichung:

Die daraus resultierenden Strommessungen entlang des Randes von $B$ liefern jeweils ein Funktion $g$ mit

wobei $\frac{\partial}{\partial n}u(x)$ die Richtungsableitung von $u$ am Randpunkt $x$ in Richtung des Normalenvektors des Randes an diesem Punkt bezeichnet. Eingangsdaten in das Modell der elektrischen Impedanz-Tomografie sind nun hinreichend viele Paare $(f,g)$ aus angelegten Spannungen und daraus resultierenden Strömen. Die Ausgabe des Modells ist die ortsabhängige Leitfähigkeit $a$. Aus der Kenntnis von Eigenschaften der Lösungen (Randwerte, Fluss über den Rand) soll also auf einen Parameter der zugrunde liegenden partiellen Differentialgleichung geschlossen werden. Diese Art von Aufgabenstellungen nennt man deshalb auch Parameteridentifikationsprobleme.

Diskretisierung: Die Diskretisierung kann mit Standardverfahren erfolgen. Allerdings muss unbedingt die Geometrie der Elektroden mit eingebracht werden um brauchbare Resultate erhalten zu können.

Implementierung: Algorithmen für das Lösen des Impedanz-Tomografie-Problems sind in gewissem Maße noch Gegenstand der Forschung. Letztlich kann man Standardverfahren kombinieren, erhält aber nur sehr schwache Aussagen über die Vertrauenswürdigkeit der Lösungen. Zu allen angelegten Spannungen kann man bei gegebener Leitfähigkeitsverteilung durch Lösen der Differentialgleichung die zu erwartenden Ströme berechnen. Diese vergleicht man mit den gemessenen Strömen und passt die Leitfähigkeitsverteilung ausgehend von einer initialen Verteilung so lange an, bis die Abweichung zwischen berechneten und gemessenen Strömen klein genug ist. Dieser algorithmische Ansatz löst also ein nichtlineares Minimierungsproblem iterativ, wobei jede einzelne Auswertung der Zielfunktion das Lösen mehrerer partieller Differentialgleichungen erfordert (eine Gleichung pro Spannungs-Strom-Paar).

1.3Fallstudien¶

Wissenschaftliches Rechnen ist ein sehr umfassendes und fachübergreifendes Thema. In dieser Veranstaltung werden wir keine systematische Abhandlung vornehmen, sondern uns anhand mehrerer umfangreicher Fallstudien mit exemplarischen Anwendungsszenarien und deren Umsetzung von der Modellierung bis zur Implementierung befassen.

Nur an einigen Punken, insbesondere mit Blick auf partielle Differentialgeleichungen, werden wir tiefer in die theoretische Fundierung und Motivation von Algorithmen einsteigen. Partielle Differentialgleichungen sind das Hauptwerkzeug für die Modellierung. Es wird sich allerdings zeigen, dass man große Klassen partieller Differentialgleichungen durch Transformation in Integralgleichungen lösen kann (und dies auch tut), sodass die Diskussion von Integralgleichungen nicht ausgleiben kann.

Hauptziele der Veranstaltung sind:

• einen Überblick über die verschiedenen Schritte von der Problemformulierung bis zur Lösung am Computer zu bekommen,

• Begriffe einordnen zu können, die bei der Beschreibung von Algorithmen aus dem Umfeld des wissenschaftlichen Rechnens auftauchen können,

• das Auge für Feinheiten zu schärfen, die durch ungeeignete Implementierung von Algorithmen zu (teils groben) Fehlern in der Lösung führen können.

1.4Sleipner A zur Mahnung¶

Als Negativbeispiel sei der Untergang der norwegischen Ölbohrplattform “Sleipner A” in der Nordsee am 23. August 1991 noch während des Baus erwähnt. Die tragende Betonkonstruktion wurde ausschließlich durch numerische Berechnungen am Computer geplant. Zum Einsatz kam die damals übliche und auchheute noch genutzte Software NASTRAN der NASA zum Lösen partieller Differentialgleichungen. Durch ungünstige Diskretisierung der zugrunde liegenden partiellen Differentialgleichung, die die in den Betonteilen wirkenden Kräfte beschreiben sollte, wurden einige Bereiche mit besonders hohen Lasten zu schwach dimensioniert. Die berechneten Lasten waren nur etwa halb so hoch wie die tatsächlichen. Die Plattform versank im 200 Meter tiefen Meer und löste durch Implosionen ein schwaches Erdbeben aus.

Die numerischen Berechnungen wurden nicht auf plausibilität geprüft und nicht mit klassischen Berechnungsverfahren für Betontragwerke verglichen. Die Berechnungssoftware trifft hier übrigens keine Schuld. Es handelt sich um einen Fehler der Bediener, die ungeeignete Parameter gewählt haben. Für mehr Informationen zu den Ursachen des Untergangs siehe The sinking of the Sleipner A offshore platform und Die Ursache für den Totalverlust der Betonplattform Sleipner A und Verweise darin.