R. Schaback

Die Numerische Mathematik befaßt sich mit der Entwicklung und der mathematischen Analyse von Verfahren, die zahlenmäßige Lösungen mathematischer Probleme berechnen.

Letztere stammen aus Anwendungsbereichen der Mathematik, z. B. in den Ingenieurwissenschaften, der Physik, der Ökonomie, usw.. Das Wissenschaftliche Rechnen ist mit der Numerischen Mathematik sehr eng verwandt, und die Grenzen sind nicht klar zu ziehen. Etwas pointiert ausgedrückt, ist man im Wissenschaftlichen Rechnen eher an der konkreten technischen Produktion der Lösung als an der mathematischen Analyse des Lösungsverfahrens interessiert, und deshalb kann man das Wissenschaftliche Rechnen auch als ein Anwendungsgebiet der Informatik sehen, in dem spezielle Rechnerstrukturen, z. B. massiv parallele Systeme, oft eine starke Rolle spielen.

Wegen der Vielfalt der technischen und wissenschaftlichen Aufgabenstellungen ist die Numerische Mathematik ein sehr breites Gebiet, das unter anderem die numerische Lösung von

• gewöhnlichen und partiellen Differentialgleichungen,
• Integralgleichungen,
• Optimierungsaufgaben,
• Approximationsproblemen für Funktionen, Kurven und Flächen,
• Eigenwert– und Verzweigungsproblemen

sowie alle numerischen Simulationen im weitesten Sinne umfaßt.

Weil sich die genannten Problemkreise in bezug auf ihren mathematischen Hintergrund und damit auch in bezug auf die sachgerechten Lösungsmethoden stark unterscheiden, erfordert die Lehre und Forschung in Numerischer Mathematik ein breites mathematisches Grundwissen. Es ist deshalb nicht möglich, im Rahmen eines kurzen Artikels alle Richtungen der Numerischen Mathematik darzustellen. Es können nur einige allgemeine, für alle numerischen Verfahren gültige Gesichtspunkte herausgearbeitet und typische Beispiele angegeben werden.

Die Verfahren der Numerischen Mathematik verwenden auf digitalen Rechenanlagen Gleitkommazahlen mit fester relativer Genauigkeit. Weil diese Zahlenmengen endlich sind, müssen numerische Verfahren zwangsläufig Fehler produzieren, und es gehört zu den zentralen Problemen der Numerischen Mathematik, diese Fehler abzuschätzen und mit dem Rechenaufwand zu vergleichen. Oft kann man durch erhöhten Rechenaufwand die Fehler reduzieren oder abschätzen, und das gilt besonders für spezielle Formen der Rechnerarithmetik, die es erlauben, gesicherte Fehlerabschätzungen numerisch zu berechnen.

Ferner sind numerische Rechnungen nicht immer stabil, d. h. kleine Fehler in den Eingabedaten eines Problems können eventuell große Fehler in der näherungsweisen Lösung bewirken. Deshalb ist die Stabilitätsanalyse ein weiteres wichtiges Aufgabenfeld der Numerischen Mathematik. Bei der Lösung linearer Gleichungssysteme \begin{eqnarray}\begin{array}{cc}\begin{array}{cc}Ax=b, & x,b\in {{\mathbb{R}}}^{n}\backslash \{0\}\end{array},A\in {{\mathbb{R}}}^{n\times n} & (1)\end{array}\end{eqnarray}

ist beispielsweise der relative Fehler der Lösung x im wesentlichen beschränkt durch den relativen Fehler der Eingabedaten A und b multipliziert mit der Kondition\begin{eqnarray}\begin{array}{cc}\kappa (A):=\Vert A\Vert \cdot \Vert {A}^{-1}\Vert \ge 1 & (2)\end{array}\end{eqnarray}

der Matrix, gemessen in einer Matrixnorm. Vom Gesichtspunkt der Stabilität her erweist es sich dann als sinnvoll, die Lösung des Systems (1) statt mit dem wohlbekannten Gaußschen Eliminationsverfahren durch eine QR–Zerlegung A = Q · R mit einer Orthogonalmatrix Q und einer oberen Dreiecksmatrix R über Rx = Q^Tb zu berechnen. Dieses Verfahren wird uns weiter unten noch einmal in ganz anderem Zusammenhang begegnen.

Die Stabilität eines numerischen Verfahrens zur Lösung eines Anwendungsproblems kann natürlich nicht besser sein als die des Problems selbst. Leider sind aber manche wichtigen Probleme schlecht gestellt, d. h. auch die exakte Lösung hängt nicht stetig von den Daten ab, von einer numerischen ganz zu schweigen. Dies gilt insbesondere für inverse Probleme, bei denen man typischerweise aus bestimmten Beobachtungen eines Systems auf dessen Eigenschaften schließen möchte (z. B. bei der Computertomographie oder der Bestimmung der Form von Körpern aus der Streuung von Schall- oder elektromagnetischen Wellen). In solchen Fällen wird durch Regularisierung eine künstliche Verbesserung der Stabilität vorgenommen, die dann auch eine stabile numerische Behandlung möglich macht, was aber auf Kosten eines zusätzlichen Fehlers oder einer Veränderung des Problems geschieht.

Viele Probleme aus den Anwendungen lassen sich nicht durch endlich viele Zahlen beschreiben, etwa dann, wenn die Eingangsdaten oder die Ausgangsdaten aus reellen Funktionen einer oder mehrerer reeller Variablen bestehen. In solchen Fällen wird durch Diskretisierung zu einem neuen Problem übergegangen, das nur auf endlich vielen Daten arbeitet. Statt einer Funktion f benutzt man endlich viele Funktionswerte f(x₁), …, f(x_n), oder man ersetzt f näherungsweise durch eine Linearkombination \begin{eqnarray}\begin{array}{cc}f\approx \displaystyle \sum _{j=1}^{n}{\alpha }_{j}{\phi }_{j} & (3)\end{array}\end{eqnarray}

von Basisfunktionen ϕ_j, die beispielsweise als algebraische oder trigonometrische Polynome gewählt werden können. In beiden Fällen kann man mit Vektoren der Länge n weiterarbeiten, aber man hat natürlich zu untersuchen, wie sich der durch diese Vereinfachung bedingte Diskretisierungsfehler auswirkt. Eine typischer Fall ist die näherungsweise Berechnung eines bestimmten Integrals durch eine gewichtete Linearkombination der Funktionswerte, etwa durch die Simpson-Regel \begin{eqnarray}\displaystyle \underset{a}{\overset{b}{\int }}f(x)dx\approx \frac{b-a}{6}\left(f(a)+4f\left(\frac{a+b}{2}\right)+f(b)\right),\end{eqnarray}

wobei die Eingabefunktion durch drei Funktionswerte diskretisiert wurde.

Die näherungsweise Darstellung von Funktionen durch Linearkombinationen (3) ist aus der Theorie der Potenz– oder Fourierreihen wohlbekannt. Mit der schnellen Fouriertransformation lassen sich letztere sehr effizient und stabil auswerten, und Varianten dieser Technik (mit diskreter schneller Cosinustransformation) bilden die Grundlage von Kompressionsverfahren wie JPEG und MPEG für Bild– und Tonsignale.

Glatte Ansatzfunktionen wie algebraische oder trigonometrische Polynome führen aber nur bei sehr glatten Funktionen f zu brauchbaren Abschätzungen des Diskretisierungsfehlers, wie man beispielsweise am Restglied der Taylorformel ablesen kann. Deshalb studiert man in der Numerischen Mathematik bevorzugt Ansatzfunktionen mit begrenzter Glätte. Deren wichtigste Vertreter sind

Splines bzw. finite Elemente als uni- bzw. multivariate stückweise polynomiale Funktionen. Das einfachste Beispiel sind stetige Splines ersten Grades, also univariate Polygonzüge. Wie alle anderen Funktionen dieser Art erfordern sie eine Diskretisierung oder Triangulation ihres Definitionsbereichs, was bei großen Raumdimensionen problematisch werden kann. Abhilfe schaffen dann Ansätze \begin{eqnarray}f(x)\approx \displaystyle \sum _{j=1}^{n}{\alpha }_{j}\phi ({\Vert x-{x}_{j}\Vert }_{2}),x,{x}_{j}\in {{\mathbb{R}}}^{d}\end{eqnarray}

mit radialen Basisfunktionen φ : [0, ∞) → ℝ, deren Theorie zwar vielversprechend, aber noch nicht weit genug entwickelt ist. Dabei tritt in gewissen Fällen das erstaunliche Phänomen auf, daß die Ergebnisse mit steigender Raumdimension d immer besser werden.

Bei einer Diskretisierung (3) ist es numerisch sinnvoll, nach der besten Darstellung zu fragen, die nur k von Null verschiedene Koeffizienten hat, damit man die Funktion mit kleinstmöglichem Aufwand speichern oder weiterverarbeiten kann. Ist f als zeitabhängiges Signal aufzufassen, so bekommt man dadurch eine optimale Datenkompression für Speicher– oder Wiedergabezwecke. Bilden die Ansatzfunktionen ϕ_j ein vollständiges Orthonormalsystem bezüglich eines Skalarprodukts (.,.) in einem Hilbertraum H, so ist für jedes f ∈ H wegen \begin{eqnarray}||f-\displaystyle \sum _{j}{\alpha }_{j}{\varphi }_{j}|{|}^{2}=\displaystyle \sum _{j}{((f,{\varphi }_{j})-{\alpha }_{j})}^{2}\end{eqnarray}

die Wahl α_j = (f, ϕ_j) für die k betragsgrößten Werte (f, ϕ_j) optimal. Diese nichtlineare Strategie wird erfolgreich im Zusammenhang mit wavelets verwendet, die eine hochinteressante moderne Erweiterung der Fourieranalyse in Richtung auf verbesserte Signaldarstellungen (3) bilden.

In vielen Anwendungen hat man Operatorgleichungen mit Integraloperatoren oder mit gewöhnlichen oder partiellen Differentialoperatoren zu lösen. Nach Diskretisierung bleibt dann normalerweise ein System von endlich vielen Gleichungen mit endlich vielen Unbekannten übrig, aber bei physikalischtechnischen Problemen kann die Anzahl der Unbekannten gigantisch sein. Davon weiter unten mehr.

Ist das resultierende System nicht linear, so wird ein weiterer Standardtrick der Numerischen Mathematik angewendet: die Linearisierung. Dabei ersetzt man ein nichtlineares Problem durch eine Folge von linearen Problemen, die sich durch lokale Annäherung des ursprünglichen Problems durch lineare Ersatzprobleme ergeben. Die Lösungen der linearen Teilprobleme sind dann oft als Iteration einer festen Abbildung zu schreiben, und man hat die Konvergenzgeschwindigkeit einer solchen Iteration zu untersuchen.

Der typische Fall einer Linearisierung und Iteration tritt schon bei einer einzigen nichtlinearen Gleichung der Form f(x) = 0 mit einer differenzierbaren reellen Funktion f auf. Das Newtonverfahren linearisiert die Funktion f an einer Stelle x₀ durch die Tangente \begin{eqnarray}{T}_{{x}_{0}}(x):=f({x}_{0})+{f}^{\prime}({x}_{0})(x-{x}_{0})\end{eqnarray}

und löst dann das lineare Ersatzproblem \({T}_{{x}_{0}}(x)=0\). Iterativ angewandt, ergibt sich damit das Verfahren \begin{eqnarray}{x}_{j+1}:={x}_{j}-f({x}_{j})/{f}^{\prime}({x}_{j}),j=0,1,2,\ldots,\end{eqnarray}

und man hat die Konvergenz sowie die Konvergenzgeschwindigkeit der Folge {x_j}_j zu untersuchen. Verallgemeinert man die Situation auf ein nichtlineares System f(x) = 0 mit f: ℝⁿ → ℝⁿ, so bekommt man eine Iteration, bei der man in jedem Schritt ein lineares Gleichungssystem \begin{eqnarray}f({x}_{j})=(\text{grad}f({x}_{j}))({x}_{j}-{x}_{j+1}),j=0,1,2,\ldots \end{eqnarray}

zu lösen hat.

Die Lösung von linearen Gleichungssystemen der allgemeinen Form (1), die entweder auf direktem Wege oder nach einer Linearisierung auftreten, bildet eine der zentralen Aufgaben der Numerischen Mathematik. Weil solche Probleme auch innerhalb von Iterationen und mit großen Raumdimensionen n auftreten können, sind Effizienzgesichtspunkte besonders wichtig. Der Rechenaufwand des klassischen Gaußschen Eliminationsverfahrens steigt mit n wie n³, und es ist ein offenes Problem, das aufwandsoptimale Rechenverfahren zu finden. Nach einer bahnbrechenden Arbeit von Volker Strassen (Numer. Math. 13, 1969) gibt es ein Verfahren mit Aufwand \({n}^{{\mathrm{log}}_{2}7}\approx {n}^{2,807}\), aber das Rennen nach immer kleineren Exponenten ist noch im Gange; siehe hierzu auch Algebra und Algorithmik. Weil das Problem ja insgesamt n² + n Eingabedaten hat, kann man erwarten, daß der Aufwand eines Optimalverfahrens mindestens proportional zu n² sein muß.

Die obige Situation betrifft nur die direkten Verfahren, die beim Rechnen mit reellen Zahlen die Lösung fehlerlos nach endlich vielen arithmetischen Einzeloperationen berechnen würden. Wenn man aber nur eine vorgegebene relative Genauigkeit ϵ einer durch Iteration einer linearen Abbildung berechneten Näherungslösung anstrebt, kann man wieder auf einen Aufwand der Größenordnung n² hoffen. Weil jede Matrix-Vektor-Multiplikation den Aufwand n² hat, darf man in solchen Verfahren nur eine begrenzte Anzahl von Matrix-Vektor-Multiplikationen einsetzen, und diese Anzahl darf nicht von n, sondern nur von ϵ abhängen. Der Einfluß der Genauigkeit auf die Iterationszahl hat die Konsequenz, daß die Stabilität der Lösung eines solchen Systems eine wesentliche Rolle spielt, und deshalb darf die Kondition der Matrix nicht allzu groß sein. Man hat hier ein typisches Beispiel für die Rückwirkung der Stabilität eines Verfahrens auf dessen Effizienz.

Ein Verfahren, das für symmetrische und positiv definite Koeffizientenmatrizen mit fester Kondition den obigen Bedingungen genügt, ist das Verfahren der konjugierten Gradienten. Es reduziert bei jedem Schritt den relativen Fehler etwa um den Faktor \begin{eqnarray}\frac{\sqrt{\kappa (A)}-1}{\sqrt{\kappa (A)}+1}\end{eqnarray}

in Abhängigkeit von der Kondition (2) der Koeffizientenmatrix. Ferner benötigt es bei jedem Schritt nur eine Matrix-Vektor-Multiplikation und ein paar Skalarprodukte. Die Analyse dieses Verfahrens ist übrigens ein schönes Beispiel für das Zusammenwirken von Argumenten aus der Optimierung, der Approximationstheorie und der linearen Algebra, aber die Details würden hier zu weit führen.

Die großen Systeme, die durch Diskretisierung partieller Differentialoperatoren entstehen, haben nur dünn besetzte Koeffizientenmatrizen, d. h. von den n² möglichen Matrixelementen sind nur höchstens k n mit einem von n unabhängigen k ≪ n von Null verschieden, wobei allerdings n in die Millionen gehen kann. Die Matrix-Vektor-Multiplikation hat dann nur einen Aufwand k n, und bei fester Kondition und vorgegebener Genauigkeit würde das Verfahren konjugierter Gradienten einen nur zu n proportionalen Aufwand haben, was alle denkbaren Effizienzwünsche erfüllen würde.

Die bei der Diskretisierung von Differentialoperatoren auftretenden Matrizen haben aber leider im Normalfall eine mit einer positiven Potenz von n anwachsende Kondition, was die obige Argumentation zunichte macht. Deshalb verwendet man raffinierte Zusatzstrategien (Mehrgittermethoden und Präkonditionierung), die unter gewissen Voraussetzungen zu einer von n unabhängig beschränkten Kondition führen. Der neue Ansatz algebraischer Mehrgitterverfahren versucht, die erzielten Effizienzgewinne auch auf allgemeinere Situationen zu erweitern.

Eine weitere interessante Problemstellung tritt z. B. nach Diskretisierung von technischen Schwingungsproblemen auf: die Berechnung der Eigenwerte symmetrischer Matrizen A ∈ ℝ^n×n. Führt man, beginnend mit A₀ := A, eine Iteration \begin{eqnarray}{A}_{j}={Q}_{j}\cdot {R}_{j},{A}_{j+1}={R}_{j}\cdot {Q}_{j}\end{eqnarray}

aus, die eine QR–Zerlegung berechnet und dann einfach die Faktoren „verkehrt“ wieder zusammenmultipliziert, so konvergieren unter gewissen Zusatzvoraussetzungen die Matrizen A_j gegen eine obere Dreiecksmatrix, die alle n Eigenwerte von A der Größe nach geordnet (!) in der Diagonale enthält. Dieses erstaunliche Verhalten eines so simplen Algorithmus ist intensiv weiter studiert worden. Mit diversen Zusatzstrategien (Reduktion auf Tridiagonalform, Spektralverschiebung, Abdividieren) kann man sowohl die Effizienz als auch die Konvergenzgeschwindigkeit ganz erheblich steigern. Unter praktisch durchaus realistischen Voraussetzungen ist der Aufwand, alle n Eigenwerte mit einer fest vorgegebenen relativen Genauigkeit auszurechnen, proportional zu n³, also vergleichbar mit dem Aufwand des Gaußschen Eliminationsverfahrens zur Lösung eines linearen Gleichungssystems. Erstaunlich ist, daß man für die Lösung eines schwierigen nichtlinearen Problems größenordnungsmäßig denselben Aufwand benötigt wie für ein direkt lösbares lineares Problem mit vergleichbaren Eingabedaten.

Weiteres kann man den vielen Lehrbüchern zur Numerischen Mathematik entnehmen, die jeweils unterschiedliche Schwerpunkte innerhalb des oben umrissenen Spektrums setzen, sowie den diversen Stichworteinträgen zu den angesprochenen Themenbereichen im vorliegenden Lexikon.

Das Fernziel der Weiterentwicklung der Numerischen Mathematik ist natürlich, zu allen Proble men ein stabiles Lösungsverfahren zu finden, dessen Rechenaufwand nach Möglichkeit nur proportional zur Anzahl der Eingabedaten ist, wenn eine feste relative Genauigkeit des Ergebnisses verlangt wird.