breakplaining zum pi-Tag: Integration nach Romberg //2864

Den pi-Tag habe ich schon des öfteren zum Anlass genommen, um mathematische Verfahren vorzustellen. Letztes Jahr etwa hatte ich pi als Nullstelle mit Hilfe des Newton-Verfahrens berechnet. Etwas länger her ist es, dass ich pi über die Fläche eines Viertelkreises mit einer Monte-Carlo-Simulation abgeschätzt habe.
Auch heuer möchte ich für die Berechnung über die Fläche gehen, setze dafür aber unseren Integrationskurs fort und nutze eine klassische numerische Integration.

Numerische Integrationsverfahren beruhen darauf, dass der Integrationsbereich in senkrechte, sehr dünne Streifen aufgeteilt wird, deren (ungefähre) Flächen mit dem Computer berechnet, und dann aufsummiert werden. Gegebenenfalls kommen noch Gewichtungsfaktoren dazu, um die Berechnung zu beschleunigen.
Das einfachste Verfahren nutzt die Trapezregel. Dabei wird die Fläche eines Streifens berechnet als Produkt der Breite mit dem arithmetischen Mittel der Funktionswerte ganz links und ganz rechts. Das Verfahren ist allerdings ziemlich ungenau, bzw. braucht relativ lang, um einen befriedend genauen Wert zu ermitteln.
Etwas besser ist die Integration mit der Simpsonregel. Dabei werden für einen Streifen drei Funktionsauswertungen benötigt. Durch diese drei Punkte lässt sich eine quadratische Parabel legen, die die Kurve passgenauer beschreibt als die Sekante bei der Trapezregel. Die Fläche unter einer Parabel lässt sich analytisch berechnen, so dass bei der Aufsummierung der berechneten Funktionswerte Gewichtsfaktoren hinzukommen, die die Konvergenz verbessern.

Um pi zu berechnen nutzen wir wieder, dass die Fläche unter f(x) = sqrt(1-x^2) im 1. Quadranten gleich pi/4 ist. Die Funktion beschreibt nämlich den Einheitskreis x^2 + y^2 = 1 (bzw. dessen obere Hälfte). Wenn wir sie also von 0 bis 1 integrieren, und das Ergebnis mal 4 nehmen, erhalten wir (einen Näherungswert für) pi.
Da die Gesamtbreite gleich 1 ist, brauchen wir uns im Folgenden nicht mehr um diesen konstanten Faktor zu kümmern. Wir bestimmen einfach die mittlere Höhe. [Das erinnert mich gerade ein wenig an die Quadratur des Kreises, wenn das Runde ins Eckige transformiert wird, aber hier muss leider ein Rechteck reichen.]
Für die Berechnung nutze ich ein (von mir modifiziertes, vereinfachtes) Romberg-Verfahren.
Bei diesem wird die Fläche zunächst mit breiten Streifen berechnet, deren Breite dann schrittweise halbiert, und schließlich auf Breite 0 extrapoliert.

Zur Veranschaulichung beginnen wir mit unserer Funktion und einem Streifen, der die ganze Breite ausmacht: Da f(0) = 1 und f(1) = 0, wäre (laut Trapezregel) die Fläche 0.5, für pi käme also im ersten Schritt das Vierfache, nämlich 2 heraus.
Der nächste zu berechnende Funktionswert ist in der Mitte von 0 und 1, also f(0.5)=sqrt(3)/2. Nach ein paar einfachen, aber langweiligen Rechnungen erhalten wir dann für die Fläche (1+sqrt(3))/4 mit einem ungefähren Zahlenwert von 0.683. Das entspricht einer weiteren Näherung für pi von 2.732 (was zumindest schon näher dran ist als 2).
Wenn wir jetzt durch diese beiden Punkte (Breite=1, Naeherung=2) und (Breite=0.5, Naeherung=1+sqrt(3)) eine Gerade Naeherung(Breite) = m*Breite + t legen, so erhalten wir für Breite=0 den Ordinatenabschnitt t als neue und bessere Näherung für pi. Das wäre in diesem Fall (wieder nach einfacher, aber langweiliger Rechnung) 2*sqrt(3) = 3.464..
Diesmal sind wir über den tatsächlichen Wert von pi hinausgeschossen, aber wenn wir das Verfahren nur oft genug wiederholen, nähern wir uns pi wieder (im Rahmen der Rechengenauigkeit) beliebig an.

Ich denke, ich konnte es begreiflich machen, wie das Romberg-Verfahren funktioniert.
Bei meinem modifizierten Verfahren berechne ich nicht die beiden Funktionswerte am Rand jedes Intervalls, sondern nur den Funktionswert in der Mitte.
Bei unserem Beispiel wäre das f(1/2) = sqrt(3)/2. Als erste Näherung für pi wäre das 2*sqrt(3), also ungefähr 3.464. (Dieses gleiche Zwischenergebnis erhielten wir nach dem üblichen Verfahren, wie weiter oben dargestellt, erst nach längerer Rechnung.)
Beim nächsten Schritt ergibt sich für f(1/4) = sqrt(15)/4 und f(3/4) = sqrt(7)/4, woraus sich als vierfacher Mittelwert als Näherung für pi (sqrt(15)+sqrt(7))/2 = 3.259 ergibt, was schon gar nicht mehr so weit von pi entfernt ist.
Legen wir wieder eine Gerade durch die beiden Ergebnispunkte und extrapolieren für die Intervallbreite 0 so ergibt sich 3.055.

Ich hab‘ den Computer mal rechnen lassen, und erhalte für die ersten Schritte n (mit 2^(n-1) Punkten, also Intervallbreite 2^(-n+1)): 3.46410162, 3.25936733, 3.18392922, 3.15668693, 3.14695181, 3.14349140, 3.14226467, 3.14183037, 3.14167672, 3.14162238, 3.14160316, ..
Wie ihr seht, nähern wir uns pi immer mehr an.
Leider hatte ich diesmal jedoch nicht mehr Zeit, um detailliertere Rechenergebnisse aufzubereiten und anschaulich darzustellen.

Üblicherweise erfolgt die Extrapolation über ein Polynom, das mit dem Neville-Aitken-Algorithmus berechnet wird. Eine einfachere und hier zur Demonstration ausreichende Möglichkeit ist es, durch die zwei neuesten Punkte eine Gerade zu legen, und den Ordinatenabschnitt zu berechnen. Wer Lust hat, kann auch eine Regressionsgerade durchlegen, aber da insbesondere die ersten Punkte, deren Breiten-Koordinaten noch weit von Null entfernt liegen, große Abweichungen haben, ist dies nicht so geschickt. Zwar gibt es dafür auch Tricks, aber die führen hier zu weit. Ihr könnt ja den Integrationshelfer eures Vertrauens fragen.

Werbung