________________________________________________________________________________
Inhalt....: Mathematik entdecken - Iterative Berechnung von PI mit MuPAD
Kategorie.: Unterrichtsmaterial
Mathematik: Analysis, Numerik, Geometrie R^2
MuPAD.....: 3.1.0
Datum.....: 2002-01-17
Autoren...: Kai Gehrs <acrowley@mupad.de>
Funktionen: |, ->, plot, plot::Circle2d, plot::Polygon2d, int, sum,
Funktionen: student::riemann, student::trapezoid, student::plotRiemann,
Funktionen: student::plotSimpson, student::plotTrapezoid, DIGITS
________________________________________________________________________________
Mathematik entdecken - Iterative Berechnung der
Zahl p
Die Berechnung des Umfangs des Einheitskreises, das bedeutet die Bestimmung
von 2π, ist ein sehr altes Problem, das vor über 4000 Jahre erstmals in Erscheinung
getreten ist und schon im Alten Testament zu finden ist. Dort verwendete man π=3.
In China, 5. Jahrhundert, fand man eine erstaunliche Näherung von π durch die
rationale Zahl 355/113, deren ersten 10 Dezimalziffern 3,141592920 lauten und
von denen die ersten 6 Nachkommstellen korrekt sind (im Vergleich: π = 3,141592654...).
Lambert bewies 1761, daß π eine irrationale Zahl ist, und Lindemann zeigte 1882,
daß π transzendent über den rationalen Zahlen ist.
Wir werden in diesem Notebook einige iterative Verfahren zur Berechnung von π
kennlernen.
Die Methode von Archimedes
________________________________________________________________________________
Archimedes war der erste, der eine systematische Lösung des Problems
gefunden hat.
Er legte in den Einheitskreis ein gleichseitiges Dreieck und verdoppelte
dann fortlaufend die Zahl der Ecken. Der Umfang der Drei-, Sechs-,
Zwölf-, ... Ecke nähert den Kreisumfang immer besser an.
Dem Schüler kann folgende Aufgabe gegeben werden: Zeige, daß im
Einheitskreis:
1. die Seitenlänge des gleichseitigen Dreiecks
beträgt.
2. das regelmäßige 2n-Eck die Seitenlänge
hat, wenn sn die Seitenlänge des n-Ecks ist.
Mit diesen Voraussetzungen können wir das Problem der Approximation
von π über den Umfang des n-Ecks in MuPAD formulieren:
Die Seitenlänge des n-Ecks ergibt sich nach oben zu:
Seitenlaenge:= n -> (
if n > 3 then sqrt(2-sqrt(4-Seitenlaenge(n/2)^2))
else sqrt(3)
end
)
n -> (if 3 < n then
sqrt(2 - sqrt(4 - Seitenlaenge(n/2)^2))
else
sqrt(3)
end_if)
Archimedes verdoppelte die Eckenzahl mit jedem Näherungsschritt:
Eckenzahl:= j -> (
if j > 1 then 2*Eckenzahl(j-1)
else 3
end
)
j -> (if 1 < j then
2*Eckenzahl(j - 1)
else
3
end_if)
Dann ergibt sich der Umfang des n-Ecks zu:
nEckUmfang:= j->Eckenzahl(j)*Seitenlaenge(Eckenzahl(j))
Diesen nehmen wir nun als Annäherung von 2π. Wir führen die ersten 10
Schritte der Methode von Archimedes durch:
n:= 10:
for j from 1 to n do
print(Unquoted, "Umfang des ".Eckenzahl(j)."-Ecks: ".
expr2text(float(nEckUmfang(j))/2)
)
end:
Umfang des 3-Ecks: 2.598076211
Umfang des 6-Ecks: 3.0
Umfang des 12-Ecks: 3.105828541
Umfang des 24-Ecks: 3.132628613
Umfang des 48-Ecks: 3.139350203
Umfang des 96-Ecks: 3.141031951
Umfang des 192-Ecks: 3.141452472
Umfang des 384-Ecks: 3.141557608
Umfang des 768-Ecks: 3.141583892
Umfang des 1536-Ecks: 3.141590463
Erhöht man sukzessive die Eckenzahl, so läßt sich ein interessanter
Vorgang beobachten. Nehmen wir beispielsweise:
n:= 40:
for j from 25 to n do
print(Unquoted, "Umfang des ".Eckenzahl(j)."-Ecks: ".
expr2text(float(nEckUmfang(j))/2)
)
end:
Umfang des 50331648-Ecks: 3.141608696
Umfang des 100663296-Ecks: 3.141674265
Umfang des 201326592-Ecks: 3.142023942
Umfang des 402653184-Ecks: 3.14307274
Umfang des 805306368-Ecks: 3.148660429
Umfang des 1610612736-Ecks: 3.159806165
Umfang des 3221225472-Ecks: 3.181980515
Umfang des 6442450944-Ecks: 3.354101966
Umfang des 12884901888-Ecks: 4.242640687
Umfang des 25769803776-Ecks: 6.0
Umfang des 51539607552-Ecks: 12.0
Umfang des 103079215104-Ecks: 24.0
Umfang des 206158430208-Ecks: 48.0
Umfang des 412316860416-Ecks: 96.0
Umfang des 824633720832-Ecks: 192.0
Umfang des 1649267441664-Ecks: 384.0
Man erkennt, daß die Näherungsfolge zu divergieren scheint! Der Grund
liegt daran, daß obige Formel zur Berechnung desUmfangs des n-Ecks
numerisch instabil ist. Es tritt Auslöschung auf, was man sich leicht
überlegen kann. Erhöht man beispielsweise die Genauigkeit der numerischen
Rechnungen, beispielsweise auf 50 Dezimalstellen, so erhält man für den 50
Iterationsschritt noch ein "vernünftiges" Ergebnis:
DIGITS:= 50:
float( nEckUmfang(50) ) / 2
delete DIGITS:
Dieses "Phänomen" bietet interessanten Diskussionsstoff für den Unterricht
und motiviert die Suche nach einer numerisch besser geeigneten Formel für
den Umfang des n-Ecks, was uns auf die Formel
führen kann.
Abschließend sollen die visuellen Fähigkeiten von MuPAD genutzt werden,
um die Methode von Archimedes grafisch darzustellen:
Kreis:= plot::Circle2d(1, [0,0], Color = RGB::Black):
pi:= float(PI):
n:= 10:
nEck:= n -> plot::Polygon2d([[sin((2*PI*k)/n), cos((2*PI*k)/n)]
$ k=0..n]):
plot( Kreis,nEck(n), Scaling = Constrained )
Die Formel von John Machin
________________________________________________________________________________
Mitte des 17. Jahrhunderts gab es die ersten Ansätze der -Berechnung über
Kettenbrüche, unendliche Produkte und Reihen.
Ein Ansatz stammt von Gregory (1671), der von der geometrischen Reihe
ausging:
delete n:
G:= 1/(1+x) = freeze(sum)( (-1)^k*x^k,k=0..n)
Nun ersetzte er x durch 
und erhielt:
G2:= G | x=x^2
Integriert man nun beide Seiten, so erhält man unmittelbar die folgende
Gleichung:
G3:= int( lhs(G2),x ) = int( rhs(G2),x )
Bekanntlich ist arctan(1)=, d.h. die rechte Seite liefert uns eine beliebig
genaue Näherung von , wenn wir für x den Wert 1 einsetzen und n
hinreichend groß wählen.
Wir nehmen uns die rechte Seite dieser Gleichung also vor und fassen
sie als Funktion in x und n auf:
Pi:= fp::unapply( rhs(G3),x,n ):
Damit erhalten wir eine erste Näherung von über:
4.0*Pi(1,2)
Error: 2nd argument: expecting an equation such as x = g(y) [intlib::alt]
Die Konvergenz ist sehr langsam, wie folgende Näherungen zeigen:
[4.0*Pi(1,n) $ n=1..50]
Eine numerisch wesentlich günstigere Formel ist
= 16 arctan(1/5) - 4 arctan(1/239),
die wir aus unseren Berechnungen von oben mit n Schritten direkt wie folgt
ausdrücken können:
Pi2:= n -> 16*Pi(1/5,n) + (-4)*Pi(1/239,n)
Sie stammt von John Machin aus dem Jahr 1706.
Im Vergleich zu Pi(50) - der letzte Eintrag der Liste - sind bereits
10 Stellen erfüllt:
float( Pi2(50) )
Damit rechnete John Machin als erster auf 100 Dezimalen aus.
Wir können in wenigen Sekunden herausfinden, wie viele Schritte er
rechnen mußte, um diesen Wert zu erhalten:
DIGITS:= 100: pi:= float(PI): n:= 1:
eps:= 10^(-DIGITS):
while abs(pi-Pi2(n)) > eps do n:= n+1 end:
n
Pi2(n)
float(%)
Neumann rechnete mit der Eniac, dem ersten elektronischen Rechenautomaten,
auf 2035 Stellen aus (er wollte die Ziffernhäufigkeiten mit statistischen Mitteln
untersuchen).
Das dauerte 70 Stunden. Soviel Zeit würden Sie uns sicher nicht zubilligen.
Mit MuPAD ist das alles einfacher und schneller: Das Ergebnis lautet
DIGITS:= 2035: float( PI )
delete DIGITS:
Numerische Integration
________________________________________________________________________________
Ein weiteres Verfahren zur Approximation von ist über das Integral
mit Hilfe numerische Verfahren wie beispielsweise über Riemannsche
Summen, die Trapez- oder die Simpsonregel.
Auch hier treten symbolische und numerische Rechnungen sowie die
Visualisierungsmöglichkeiten von CAS in Erscheinung.
Wir geben obige Funktion in MuPAD ein:
f:= x -> sqrt(1-x^2):
und erhalten die Riemann-Summen mit n Iterationen über die von MuPAD
mitgelieferte Funktion riemann( f(x),x=a..b,n ):
export(student):
F:= riemann( f,x=0..1,4 )
unfreeze( F )
In Gleitkommadarstellung:
4 * float( F )
Entsprechend bestimmen wir eine Approximation von über die
Trapezregel:
F:= trapezoid( f,x=0..1,4 )
4 * float( F )
und abschließend im Vergleich die Simpsonregel:
F:= simpson( f,x=0..1,4 )
4 * float( F )
Natürlich können diese Formeln erst im Unterricht hergeleitet, um sie dann
in MuPAD in Form von kleinen Prozeduren zu implementieren. Aufgrund der
Kürze der Zeit habe ich die von MuPAD bereits implementierten Prozeduren
verwendet.
Solche numerischen Berechnungen geben Aufschluß über die Güte der
Integrationsverfahren:
for n in [6,10,20,50,100] do
print(Unquoted,
"n = ".n.":\n".
" Rechteck: ".expr2text(4*float(riemann(f,x=0..1,n)))."\n".
" Trapez : ".expr2text(4*float(trapezoid(f,x=0..1,n)))."\n".
" Simpson : ".expr2text(4*float(simpson(f,x=0..1,n)))."\n".
"--------------------------"
)
end:
Stoff, der zum Experimentieren und zu Diskussionen anregt.
Zum Abschluß verwenden wir die Visualisierungsmöglichkeiten von MuPAD,
um den Prozeß der numerischen Approximation anschaulich zu verdeutlichen:
plot( student::plotRiemann( f(x),x=0..1,8 ),
Scaling = Constrained )
plot( student::plotTrapezoid( f(x),x=0..1,8 ),
Scaling = Constrained )
plot( student::plotSimpson( f(x),x=0..1,8 ),
Scaling = Constrained )
Moderne Methoden der Computer Algebra
________________________________________________________________________________
Mit modernen Methoden der Mathematik aus dem Bereich der Computer-
algebra und schnellen Algorithmen für Ganzzahl- und Gleitkommaarithmetik
mit hoher Genauigkeit, basierend auf der "Fast Fourier Transformation" und
"Schneller Division" wird heute die Jagd nach dem Weltrekord für die Nach-
kommastellen von betrieben.
Derzeit liegt der Weltrekord bei mehreren Billionen Ziffern, der aber wohl bald
wieder übertroffen werden dürfte.
Solche Berechnungen dienen als gute Tests für Hardware von Supercomputern,
bevor sie zum Verkauf angeboten werden.
Mit Computer-Algebra-Systemen wie MuPAD lassen sich bereits tausende von
Nachkommstellen von in sekundenschnelle bestimmen.
In 200 Millisekunden berechnet MuPAD beispielsweise auf 10.000 Stellen genau,
wobei bei folgender Eingabe die meiste Zeit für die Ausgabe des Ergebnisses
benötigt wird:
DIGITS:= 10000:
time(float(PI));
delete DIGITS:
Bis heute ist unklar, ob die Dezimalstellen systematisch angeordnet sind oder
aber zufallsverteilt sind.
Wir wissen ebenso wenig, ob beispielsweise die Ziffer 1 unendlich oft in der
Dezimalentwicklung von vorkommt.
________________________________________________________________________________
Anmerkungen:
1. Weitere Anregungen finden Sie unter: http://schule.mupad.de bzw. http://studium.mupad.de
_______________________________________________________________________________