Numerische Simulation

Das Doppelintegral der Schubart Master-Formel lässt sich in den meisten Fällen nicht geschlossen auswerten und muss numerisch behandelt werden. Dieses Kapitel beschreibt die Methodik einer solchen numerischen Auswertung: die Diskretisierung der Integrationsdomänen über Energie und Winkel, die Wahl geeigneter Quadraturverfahren und die Kriterien zur Beurteilung der Konvergenz. Die Darstellung bleibt generisch und methodisch; sie nennt keine erfundenen konkreten Ergebniszahlen, und wo Wertebereiche genannt werden, sind sie ausdrücklich als illustrativ gekennzeichnet.

Diskretisierung der Integrationsdomänen

Die numerische Auswertung beginnt mit der Zerlegung der kontinuierlichen Integration in eine endliche Summe. Die Energiedomäne wird in eine Reihe von Stützstellen $E_j$ unterteilt, die Winkeldomäne in Stützstellen $\theta_k$. Der Integrand wird an jedem Gitterpunkt ausgewertet und mit dem zugehörigen Quadraturgewicht multipliziert. In der einfachsten Form ergibt sich die diskrete Näherung

$$P(t) \approx \sum_{j}\sum_{k} \Phi(E_j,\theta_k,t)\,\sigma(E_j)\,\Delta p\,v_{\text{ph}}\,\eta_{\text{conv}}\;w_j\,w_k,$$

wobei $w_j$ und $w_k$ die Quadraturgewichte in Energie- und Winkelrichtung sind. Die Feinheit dieser Diskretisierung bestimmt die Genauigkeit der Näherung: Ein zu grobes Gitter kann scharfe Strukturen im Fluss oder im Streuquerschnitt übersehen, während ein zu feines Gitter den Rechenaufwand unnötig erhöht.

Die Struktur der Doppelsumme entspricht einem Tensorprodukt zweier eindimensionaler Quadraturregeln. Legt man in Energierichtung $N_E$ Stützstellen und in Winkelrichtung $N_\theta$ Stützstellen an, so umfasst das Gitter $N_E\cdot N_\theta$ Auswertungspunkte, an denen der volle Integrand berechnet werden muss. Dieser produktartige Aufwand ist der Grund, warum die Wahl der Auflösung in jeder Dimension mit Bedacht zu treffen ist: Eine Verdopplung der Stützstellen in beiden Richtungen vervierfacht die Zahl der Funktionsauswertungen. Ein sinnvoller erster Schritt jeder Untersuchung ist daher, die beiden Richtungen zunächst getrennt zu betrachten und die jeweils erforderliche Mindestauflösung separat zu bestimmen.

Zwischenschritte einer konkreten Auswertung

Eine praktische Auswertung durchläuft typischerweise vier Zwischenschritte. Zuerst werden die Integrationsgrenzen in Energie und Winkel festgelegt; die Energieobergrenze ist so zu wählen, dass der abgeschnittene Beitrag jenseits der Grenze vernachlässigbar bleibt. Zweitens werden Stützstellen und Gewichte des gewählten Quadraturverfahrens erzeugt. Drittens wird der Integrand an jedem Gitterpunkt ausgewertet, wobei tabellierte Eingangsgrößen wie $\sigma(E)$ gegebenenfalls interpoliert werden. Viertens wird die gewichtete Summe gebildet. Jeder dieser Schritte trägt seine eigene Fehlerquelle bei, die getrennt zu kontrollieren ist.

Wahl des Quadraturverfahrens

Die Wahl des Quadraturverfahrens richtet sich nach der Glattheit des Integranden. Für glatte Funktionen eignen sich Verfahren höherer Ordnung, etwa Gauß-Quadratur, die mit wenigen Stützstellen hohe Genauigkeit erreichen. Weist der Fluss oder der Streuquerschnitt scharfe Resonanzen oder Kanten auf, so sind adaptive Verfahren vorzuziehen, die die Stützstellendichte lokal an die Struktur des Integranden anpassen.

Besondere Aufmerksamkeit verdient die Energieabhängigkeit des Streuquerschnitts $\sigma(E)$. Besitzt dieser ausgeprägte Maxima, so muss die Energiediskretisierung fein genug sein, um diese Bereiche aufzulösen. Eine ungeeignete Diskretisierung in der Nähe solcher Strukturen ist eine häufige Fehlerquelle und würde zu einer systematischen Über- oder Unterschätzung der Leistung führen.

Die Konvergenzordnung eines Quadraturverfahrens beschreibt, wie schnell der Diskretisierungsfehler mit wachsender Stützstellenzahl abnimmt. Für ein Verfahren der Ordnung $p$ skaliert der Fehler bei glattem Integranden asymptotisch wie

$$\varepsilon(h) \sim C\,h^{p},$$

wobei $h$ die typische Gitterweite und $C$ eine von den höheren Ableitungen des Integranden abhängige Konstante ist. Die zusammengesetzte Trapezregel besitzt die Ordnung $p=2$, die Simpson-Regel die Ordnung $p=4$; Gauß-Quadratur mit $n$ Knoten integriert Polynome bis zum Grad $2n-1$ exakt und erreicht damit für glatte Funktionen eine sehr rasche Konvergenz. Entscheidend ist die Voraussetzung der Glattheit: An Kanten oder Sprüngen des Integranden bricht die theoretische Ordnung zusammen, und der beobachtete Fehler fällt nur noch langsam. In solchen Fällen ist die Integrationsdomäne an den Unstetigkeitsstellen zu zerlegen, sodass innerhalb jedes Teilintervalls wieder ein glatter Integrand vorliegt.

Konvergenzanalyse

Die zentrale Frage jeder numerischen Auswertung ist, ob das Ergebnis mit zunehmender Verfeinerung des Gitters gegen einen stabilen Wert konvergiert. Die übliche Vorgehensweise besteht darin, die Diskretisierung schrittweise zu verfeinern und die relative Änderung des Ergebnisses zu verfolgen. Sinkt diese Änderung unter eine vorgegebene Toleranz, so gilt das Ergebnis als konvergiert.

Als rein illustrative Größenordnung — nicht als gemessenes Ergebnis — sei genannt, dass in typischen glatten Fällen eine relative Änderung im Bereich von etwa einem Prozent oder darunter als hinreichende Konvergenz betrachtet werden kann; der konkrete Schwellenwert hängt jedoch vom Anwendungsfall ab und ist stets gesondert festzulegen. Die Konvergenzanalyse sollte getrennt für die Energie- und die Winkelrichtung durchgeführt werden, da beide Domänen unterschiedliche Auflösungsanforderungen stellen können.

Eine praktische Fehlerabschätzung nutzt die bekannte Konvergenzordnung des Verfahrens. Vergleicht man das Ergebnis auf einem Gitter der Weite $h$ mit dem auf einem halbierten Gitter $h/2$, so lässt sich für ein Verfahren der Ordnung $p$ der verbleibende Diskretisierungsfehler durch Richardson-Extrapolation schätzen:

$$\varepsilon(h/2) \approx \frac{P_{h/2} - P_{h}}{2^{p} - 1}.$$

Diese Beziehung liefert nicht nur eine Fehlerschranke, sondern zugleich einen genaueren extrapolierten Wert $P_{h/2} + \varepsilon(h/2)$. Beobachtet man beim Halbieren der Gitterweite nicht die theoretisch erwartete Fehlerreduktion um den Faktor $2^{p}$, so ist dies ein deutliches Warnzeichen dafür, dass der Integrand die vorausgesetzte Glattheit verletzt oder eine scharfe Struktur noch unzureichend aufgelöst wird.

Behandlung der Zeitabhängigkeit

Da $\Phi(E,\theta,t)$ zeitabhängig ist, muss die Integration für jeden betrachteten Zeitpunkt erneut ausgewertet werden. Bei der Berechnung eines zeitlichen Verlaufs von $P(t)$ wird die räumliche Integration in einer äußeren Zeitschleife wiederholt. Die Wahl der zeitlichen Auflösung richtet sich nach der Variabilität des Flusses, wie in Systemgrenzen im Zusammenhang mit der quasistationären Näherung diskutiert.

Fehlerquellen und Validierung

Neben dem Diskretisierungsfehler treten weitere Fehlerquellen auf: Rundungsfehler bei der Summation vieler kleiner Beiträge, Fehler bei der Interpolation tabellierter Eingangsgrößen und mögliche systematische Fehler durch eine unzureichende Abdeckung der Integrationsdomäne. Eine sorgfältige numerische Untersuchung dokumentiert diese Fehlerquellen und schätzt ihren Einfluss ab.

Ein wirksames Mittel gegen Rundungsfehler bei der Summation vieler kleiner Beiträge ist die kompensierte Summation, etwa nach dem Kahan-Verfahren, das den bei jeder Addition verlorenen niederwertigen Anteil mitführt und im nächsten Schritt zurückaddiert. Zur Validierung der Implementierung selbst empfiehlt sich die Methode der hergestellten Lösung: Man setzt einen künstlichen Integranden ein, dessen Integral geschlossen bekannt ist, und prüft, ob der Code diesen Wert innerhalb der erwarteten Ordnung reproduziert. Erst nach bestandenem Test dieser Art ist eine Auswertung mit den tatsächlichen Eingangsgrößen belastbar.

Die numerische Integration liefert einen deterministischen Schätzwert für ein gegebenes Set von Eingangsgrößen. Um die Robustheit dieses Schätzwerts gegenüber Unsicherheiten in den Eingangsgrößen zu beurteilen, ergänzt man die deterministische Integration durch statistische Verfahren. Diese werden in Monte-Carlo-Validierung beschrieben. Die mathematische Grundlage der hier diskretisierten Beziehung findet sich in Mathematische Architektur; interaktiv lässt sich das Verhalten des Integrals im Playground erkunden. Die Methodik dieses Kapitels ist damit der Ausgangspunkt jeder quantitativen Untersuchung des Rahmenwerks, ohne selbst empirische Aussagen über konkrete Zahlenwerte zu treffen.

Weiterführende Literatur

• Press, W. H., Teukolsky, S. A., Vetterling, W. T. & Flannery, B. P. (2007). Numerical Recipes: The Art of Scientific Computing (3. Aufl.). Cambridge University Press. — Standardwerk zu Quadratur, adaptiver Integration und Fehlerkontrolle.
• Davis, P. J. & Rabinowitz, P. (1984). Methods of Numerical Integration (2. Aufl.). Academic Press. — Umfassende Behandlung ein- und mehrdimensionaler Quadraturverfahren.
• Stoer, J. & Bulirsch, R. (2002). Introduction to Numerical Analysis (3. Aufl.). Springer. — Konvergenztheorie, Richardson-Extrapolation und Fehleranalyse.