Methoden: Diskretisierung, Sampling und Konvergenz

Dieses Kapitel legt im Stil eines Methodenteils die Rechenverfahren dar, mit denen eine Zahl aus der Schubart Master-Formel gewonnen wird. Es schreitet von der kontinuierlichen Grundbilanz zu ihrer praktischen Auswertung in fünf nummerierten Schritten fort: das exakte Integral, seine Reduktion auf eine repräsentative Scheibe, seine Diskretisierung als Quadratursumme, seine Schätzung durch Monte Carlo und die Fortpflanzung der Eingangsunsicherheit. Jede Schlüsselbeziehung ist nummeriert, sodass der umgebende Text und spätere Kapitel unmissverständlich als „Gl. (n)” darauf verweisen können. Die vollständige analytische Herleitung, die diese Schritte motiviert, findet sich in der vollständigen Herleitung; hier liegt der Schwerpunkt darauf, wie das Integral tatsächlich berechnet wird, und auf den Fehlern, die jede Methode trägt.

Gleichung (1) — die erweiterte Bilanz

Der Ausgangspunkt ist die erweiterte physikalische Darstellung der Schubart Master-Formel — die Bilanzidentität für die momentane gewandelte Leistung in ihrer spektral-transport-Form,

$$P(t) = \iint \Phi(E,\theta,t)\,\sigma(E)\,\Delta p\,v_{\text{ph}}\,\eta_{\text{conv}}\;dE\,d\theta \tag{1}$$

wobei die Integration über die einfallende Energie $E$ und die Winkelkoordinate $\theta$ läuft und die fünf Faktoren die Flussdichte, der Streuquerschnitt, der effektive Impulsübertrag, die Phononengeschwindigkeit und der Wandlungswirkungsgrad sind. Gl. (1) ist ein Doppelintegral über einen zweidimensionalen Bereich; jede nachfolgende Methode in diesem Kapitel ist eine Art, es auszuwerten. Die dimensionelle Buchführung, die zeigt, dass sich die rechte Seite auf Watt reduziert, wird in der Dimensionsanalyse durchgeführt und hier vorausgesetzt.

Gleichung (2) — repräsentative Scheibe

Wenn der Integrand über ein gut gewähltes Fenster $\Delta E \times \Delta\theta$ langsam veränderlich ist, oder wenn eine erste Schätzung gewünscht wird, reduziert sich Gl. (1) auf eine einzelne repräsentative Scheibe, indem das Integral durch den an repräsentativen Werten ausgewerteten Integranden mal der Fenstergröße ersetzt wird,

$$P \approx \Phi\,\sigma\,\Delta p\,v_{\text{ph}}\,\eta_{\text{conv}}\;\Delta E\,\Delta\theta . \tag{2}$$

Gl. (2) ist die Mittelwertform von Gl. (1): sie ist exakt, wenn die Faktoren über das Fenster konstant sind, und andernfalls eine kontrollierte Näherung. Ihre Gültigkeit ist genau die Bedingung der repräsentativen Mittelung aus Annahme C in den Annahmen und Grenzen — das Fenster muss breit genug sein, um repräsentativ zu sein, doch schmal genug, dass der Integrand nicht über es hinweg schwingt. Gl. (2) ist unschätzbar für Größenordnungsschätzungen und zur Prüfung, dass eine vollständigere Rechnung nicht in die Irre gegangen ist, darf aber nicht dort vertraut werden, wo der Integrand scharf gipfelt, etwa über eine Wirkungsquerschnitt-Resonanz.

Gleichung (3) — Quadratur-Diskretisierung

Um Gl. (1) ohne die Annahme einer einzelnen repräsentativen Scheibe auszuwerten, wird der Bereich in ein Gitter aus Zellen zerlegt, indiziert durch $j$ (in der Energie) und $k$ (im Winkel), und das Integral durch eine Riemann- oder Quadratursumme über diese Zellen angenähert,

$$P \approx \sum_{j,k} \Phi_{jk}\,\sigma_j\,\Delta p\,v_{\text{ph}}\,\eta_{\text{conv}}\;\Delta E_j\,\Delta\theta_k , \tag{3}$$

wobei $\Phi_{jk}$ und $\sigma_j$ die an den Zellmitten abgetasteten Faktorwerte und $\Delta E_j$, $\Delta\theta_k$ die Zellbreiten sind. Gl. (3) ist schlicht Gl. (2), summiert über viele Scheiben; wenn die Zellen schrumpfen, konvergiert sie gegen Gl. (1). Eine gleichförmige Mittelpunktsregel wie geschrieben hat einen Diskretisierungsfehler pro Zelle der Ordnung $(\Delta E)^3$ und $(\Delta\theta)^3$ in der zweiten Ableitung des Integranden, sodass der globale Fehler der Mittelpunkts- oder Trapezregel wie $O(h^2)$ in der Zellgröße $h$ skaliert; die Simpson-Regel verbessert dies auf glatten Integranden zu $O(h^4)$. Nichtgleichförmige Gitter, die Zellen dort dichter setzen, wo der Integrand am schnellsten variiert — etwa um einen Spektralgipfel —, gewinnen die Genauigkeit weit günstiger zurück als eine gleichförmige Verfeinerung. Diese gitterbasierte Auswertung ist das Rückgrat der deterministischen numerischen Simulation.

Gleichung (4) — der Monte-Carlo-Schätzer

Die Gitterquadratur wird teuer, wenn der Integrand hochdimensional, scharf gegipfelt oder unhandlich geformt ist, weil die Zahl der Zellen wie eine Potenz der Dimension wächst. Die Monte-Carlo-Alternative schreibt das Integral als Erwartungswert um und schätzt es durch zufälliges Sampling. Zieht man $N$ Punkte $x_n = (E_n,\theta_n)$ aus dem Bereich und schreibt $f(x)$ für den vollen Integranden von Gl. (1) mal dem Bereichsvolumen, so ist der Schätzer das Stichprobenmittel

$$\hat P = \frac{1}{N}\sum_{n=1}^{N} f(x_n), \qquad \mathrm{SE}\big[\hat P\big] \propto \frac{1}{\sqrt{N}} . \tag{4}$$

Gl. (4) ist erwartungstreu — ihr Erwartungswert gleicht dem wahren Integral — und ihr Standardfehler fällt wie $1/\sqrt{N}$, unabhängig von der Dimension des Bereichs. Diese Dimensionsunabhängigkeit ist der Grund, warum Monte Carlo die Tensorgitter-Quadratur bei komplizierten Integranden überholt, und sie ist auch der Grund, warum sich die Methode natürlich auf die Unsicherheitsfortpflanzung erweitert: indem man die Eingangsfaktoren selbst in einer äußeren Schleife aus ihren Unsicherheitsverteilungen zieht, liefert jeder Monte-Carlo-Durchlauf eine Realisierung von $P$, und das Ensemble reproduziert die Ausgangsverteilung. Das vollständige geschachtelte Verfahren ist in der Monte-Carlo-Validierung dargelegt.

Gleichung (5) — Fehlerfortpflanzung

Die letzte Methode betrifft nicht den Wert des Integrals, sondern seine Unsicherheit. Da Gl. (1) ein Produkt von Faktoren ist, addieren sich die relativen Varianzen der Faktoren. Für $P = \prod_i x_i$ mit statistisch unabhängigen Faktoren liefert eine Fortpflanzung erster Ordnung

$$\left(\frac{\sigma_P}{P}\right)^{2} = \sum_i \left(\frac{\sigma_{x_i}}{x_i}\right)^{2} . \tag{5}$$

Gl. (5) folgt aus dem Log-Differential $\mathrm{d}\ln P = \sum_i \mathrm{d}\ln x_i$ und wird vollständig, samt der Kovarianzterme für korrelierte Eingänge, in der Unsicherheitsanalyse hergeleitet. Ihre Rolle unter den Methoden besteht darin, uns vor jedem teuren Sampling zu sagen, welcher Eingang den Ausgangsfehler dominiert: weil sich die Terme quadratisch addieren, beherrscht die größte relative Unsicherheit das Ergebnis, und der Verfeinerungsaufwand sollte dorthin gelenkt werden. Im Grenzfall kleiner Unsicherheiten reproduziert die Monte-Carlo-Streuung aus Gl. (4) die Gl. (5), sodass die beiden als wechselseitige Konsistenztests dienen.

Konvergenz

Die drei Auswertungsmethoden — Gl. (2), (3) und (4) — tragen qualitativ verschiedenes Konvergenzverhalten, und die Wahl zwischen ihnen ist eine Frage der Anpassung der Methode an den Integranden.

Die repräsentative Scheibe, Gl. (2), konvergiert in keinem Verfeinerungssinn; sie ist eine feste Näherung, deren Fehler vollständig davon bestimmt wird, wie repräsentativ das gewählte Fenster ist. Sie ist eine erste Schätzung, kein konvergentes Schema.

Die Quadratur, Gl. (3), konvergiert wie eine Potenz der Zellgröße: $O(h^2)$ für die Mittelpunkts- und Trapezregel, $O(h^4)$ für die Simpson-Regel, auf glatten Integranden. Eine Halbierung der Zellgröße viertelt daher den Trapezfehler. Diese rasche Konvergenz hoher Ordnung macht die Quadratur zur Methode der Wahl für glatte, niedrigdimensionale Integranden — genau das Regime der zweidimensionalen Gl. (1), wenn die Faktoren gutartig sind. Ihre Schwäche ist, dass die Zellzahl wie eine Potenz der Dimension wächst und dass scharfe Gipfel sehr feine Gitter erzwingen.

Monte Carlo, Gl. (4), konvergiert wie $1/\sqrt{N}$: den Fehler zehnfach zu verringern verlangt eine hundertfache Erhöhung der Stichproben. Das ist in niedrigen Dimensionen langsam verglichen mit der Quadratur, doch der Exponent $-\tfrac12$ ist dimensionsunabhängig, sodass Monte Carlo für hochdimensionale oder scharf strukturierte Integranden gewinnt, wo das Quadraturgitter explodiert. Seine langsame Rate wird routinemäßig durch Varianzreduktion verbessert: Importance Sampling, das Ziehungen dort konzentriert, wo der Integrand am größten ist, und mit $f/g$ umgewichtet; stratifiziertes Sampling, das die Abdeckung wichtiger Teilbereiche garantiert; und Kontrollvariate, die eine korrelierte Größe mit bekanntem Integral abziehen. Jede dieser Techniken senkt die effektive Varianz $\mathrm{Var}[f]$ im Zähler des Standardfehlers, ohne die $1/\sqrt{N}$-Skalierung zu ändern, sodass ein moderates $N$ erreicht, was sonst um Größenordnungen mehr gleichförmige Stichproben erfordern würde.

In der Praxis nutzt eine robuste Auswertung von Gl. (1) alle drei: Gl. (2) für eine Plausibilitäts-Größenordnung, Gl. (3) für den deterministischen Wert, wenn der Integrand glatt ist, und Gl. (4) mit Varianzreduktion für die Unsicherheitsfortpflanzung und für jedes Regime, in dem das Gitter unpraktikabel wird. Das Fortpflanzungsgesetz Gl. (5) liefert dann das analytische Fehlerbudget, das das Sampling reproduzieren soll.

Weiterführend

• W. H. Press, S. A. Teukolsky, W. T. Vetterling und B. P. Flannery, Numerical Recipes: The Art of Scientific Computing, 3. Aufl. (Cambridge University Press, 2007) — zu Quadraturregeln, ihren Fehlerordnungen und Monte-Carlo-Integration mit Varianzreduktion.
• N. Metropolis und S. Ulam, „The Monte Carlo Method,” Journal of the American Statistical Association 44 (1949), 335–341 — die Gründungsarbeit der Methode, die Gl. (4) zugrunde liegt.
• C. P. Robert und G. Casella, Monte Carlo Statistical Methods, 2. Aufl. (Springer, 2004) — eine gründliche moderne Behandlung von Schätzern, Konvergenz und Varianzreduktion.
• J. R. Taylor, An Introduction to Error Analysis, 2. Aufl. (University Science Books, 1997) — zum Fortpflanzungsgesetz von Gl. (5).