breakplaining zum #piTag: numerische Nullstellensuche //2698

Ja, schon wieder Jahrestag. Heute ist pi-Tag, den ich in den letzten Jahren auf vielfältige Weise zelebriert habe.
Beispielsweise hatte ich 2019 demonstriert, wie man pi numerisch über ein Monte-Carlo-Verfahren berechnen kann.
Heute möchte ich zur Feier des Tages zwei andere numerisches Verfahren vorstellen, mit deren Hilfe man pi berechnen kann.
Beide beruhen darauf, dass man Nullstellen einer gegebenen Funktion sucht. Wir wissen, dass sin(x) bei (Vielfachen von) pi eine Nullstelle hat. Also werden wir im Folgenden diese spezielle Nullstelle suchen. Dabei gehen wir einfach mal davon aus, dass wir Winkelfunktionen beliebig genau berechnen können, ohne den Wert von pi zu kennen.

Nur ganz kurz möchte ich das Intervallhalbierungsverfahren beschreiben. Wenn wir davon ausgehen, dass die gesuchte Nullstelle zwischen 2 und 4 liegt, halbieren wir dieses Intervall, und erhalten 3 als erste Approximation.
Der Trick ist es, jetzt herauszufinden, ob die Nullstelle zwischen 2 und 3, oder zwischen 3 und 4 liegt. Die Mathematik sagt uns, dass (unter gewissen Voraussetzungen, die hier zutreffen) in der Umgebung von (ungeraden) Nullstellen ein Vorzeichenwechsel stattfindet. Rechnerisch erhält man das, indem man die Funktionswerte der Intervallgrenzen multipliziert. Ist das Ergebnis negativ, muss die Funktion in diesem Intervall das Vorzeichen gewechselt haben, und (mindestens) eine Nullstelle haben. sin(3)*sin(4) < 0. Das Intervall hat sich jetzt halbiert. Unsere nächste Approximation befindet sich nun bei 3.5, die übernächste bei 3.25, dann 3.125, und so weiter und so fort. Der Wert nähert sich langsam, aber sicher immer mehr pi an.
Dieses Verfahren kann man noch verbessern, indem man nicht einfach stur den Mittelwert der Intervallgrenzen nimmt, sondern geschickt interpoliert (Regula falsi).

Etwas ausführlicher möchte ich das Newton-Verfahren darstellen, das ebenfalls der Nullstellensuche einer Funktion dient. Es hat seine Nachteile und Schwächen, auf die ich aber gar nicht eingehen will, ist dafür aber sehr einfach zu verstehen und zu realisieren.
Wir stellen uns eine Funktion f(x) vor, die in einem bestimmten Bereich eine Nullstelle hat. In diesem Bereich ist sie monoton und glatt, also differenzierbar, hat keine Extrema, Polstellen, Definitionslücken, Unstetigkeitsstellen, noch sonstige Gemeinheiten oder pathologisches Verhalten.
Wir raten jetzt einen Startwert x[0], von dem wir annehmen, dass er in der Nähe der Nullstelle ist. An den Startpunkt (x[0], f(x[0]) legen wir eine Tangente y(x) = m*x + t an. Dabei ist m die Steigung von f am Punkt x[0], nämlich f'(x) = df(x)/dx, und t der Ordinatenabschnitt der Tangente, den wir aber nicht zu kennen brauchen, da er bei der weiteren Berechnung wegfällt. Diese Tangente schneidet die x-Achse an einem Punkt, der i.A. näher an der Nullstelle liegt. Wir nehmen ihn als neuen Startwert x[i], legen wieder die Tangente an, suchen den Schnittpunkt mit der x-Achse, und wiederholen dies so oft, bis nach der i-ten Iteration f(x[i]) mit genügender Genauigkeit gleich 0 ist.
x[i+1] lässt sich nach geometrischer Betrachtung des Steigungsdreiecks (gegeben durch die Punkte (x[i],f(x[i])), (x[i],0) und (x[i+1],0)) rekursiv aus dem vorigen Punkt berechnen zu: x[i+1] = x[i] – f(x[i])/f'(x[i]). Das ist eine sehr einfache Formel, für die wir nur die Funktion selbst kennen müssen, sowie ihre Ableitung. Sinnvoll ist sie für Polynome, wobei zu berücksichtigen ist, dass Polynome i.A. mehrere Nullstellen haben, die durch Extrema voneinander getrennt sind.

Da sin(x) bei pi eine Nullstelle hat, können wir das Newton-Verfahren benutzen, um pi zu berechnen.
Die Funktion ist also f(x) = sin (x). Ihre Ableitung f'(x) = cos(x). Als Startwert nehmen wir x[0] = 2. Dann ist x[i+1] = x[i] – sin(x[i])/cos(x[i]) = x[i] – tan(x[i]).
Ich habe mal ein Progrämmchen geschrieben, das genau das umsetzt. Als Genauigkeit habe ich 1E-12 vorgesehen. Unten im Zusatzbereich seht ihr den Durchlauf tabellarisch dargestellt (wobei der Aufwand für eine einigermaßen ansprechende Formatierung um einiges größer war als für den mathematischen Part).
Wie man sieht, ist bereits nach wenigen Iterationen ein guter Näherungswert für pi erreicht.

So, darauf jetzt in aller Pietät ein Stück kugelige Schopilade!

---

i	xi	sin(xi)	xi – π
0	2.00000000000000	9.0929743E-01	-1.1415927E+00
1	4.18503986326152	-8.6414415E-01	1.0434472E+00
2	2.46789367451467	6.2388107E-01	-6.7369898E-01
3	3.26618627756911	-1.2427152E-01	1.2459362E-01
4	3.14094391231764	6.4874123E-04	-6.4874127E-04
5	3.1415926536808	-9.1010899E-11	9.1010899E-11
6	3.14159265358979	-5.4210109E-20	0.0000000E+00

Werbung