Pisałem o tym, jak metodą Archimedesa obliczano liczbę . Można tę starożytną metodę sformułować jako algorytm niezależny od geometrii. Bierzemy dwie dodatnie liczby rzeczywiste
(niech
). Potem rekurencyjnie obliczamy kolejne wartości ciągów
Kolejne wyrazy są tu średnimi harmonicznymi i średnimi geometrycznymi poprzednich wyrazów. Średnia harmoniczna dwóch liczb to np. średnia prędkość na pewnej drodze, gdy połowę drogi jedziemy z pierwszą prędością, a drugą połowę z drugą prędkością. Można pokazać bez trudu (*), że ciąg jest malejący, a ciąg
rosnący. Ponieważ oba są ograniczone, muszą być zbieżne, i to do wspólnej granicy równej
Wynik ten znał nauczyciel Carla Friedricha Gaussa Johann Friedrich Pfaff, uważany za najwybitniejszego niemieckiego matematyka epoki przed Gaussem.
Biorąc oraz
, otrzymujemy w granicy liczbę
.
Gdy , należy zamienić kolejność pod pierwiastkiem oraz arcus cosinus zamienić na arcosh. W roku 1880 odkrył ponownie ten algorytm Carl Wilhelm Borchardt, znany najbardziej jako redaktor „Journal für die reine und angewandte Mathematik”, zwanego też Żurnalem Crelle’a od nazwiska pierwszego redaktora. Dlatego w literaturze nazywa się go algorytmem Borchardta albo Archimedesa-Borchardta.
Młody Gauss od dziecka zapowiadał się na wyjątkowy talent matematyczny i w tym przypadku cudowne dziecko wyrosło na czołowego matematyka Europy. Podobnie jak Euler należał on do matematyków, którzy lubią i potrafią sprawnie wykonywać rozmaite rachunki numeryczne.
Zainteresował się on następującym algorytmem:
gdzie . Zauważmy, że jest to „kuzyn” algorytmu Archimedesa: średnie harmoniczne zostały zastąpione tu średnimi arytmetycznymi. Łatwo wykazać, że oba ciągi dążą do wspólnej granicy, którą oznaczymy
i będziemy nazywać średnią arytmetyczno-geometryczną obu wyjściowych liczb. Oto przykładowe rachunki Gaussa,
,
. Widać, że zbieżność jest niezwykle szybka (dokładność danych w tabeli odpowiada rachunkom Gaussa, który oczywiście musiał przeprowadzać je ręcznie z całym mozołem, ale też chyba i radością.
Mamy więc znakomity algorytm, ale nie wiemy, co jest jego granicą. Gauss potrafił udowodnić, że granica jest w tym przypadku związana z długością lemniskaty Bernoulliego, eleganckiej krzywej, przez którą bracia Jakob i Johann Bernoulli skłócili się śmiertelnie w 1694 roku (poszło o kwestie pierwszeństwa – nie tylko w tamtej epoce traktowane niezwykle ambicjonalnie, choć dziś trochę więcej wiemy o nieuniknionej równoległości pewnych odkryć i rozumowań). Oto lemniskata.
Jest to miejsce geometryczne punków, których iloczyn odległości od dwóch ognisk jest stały (przy warunku, że środek odcinka łączącego ogniska leży na krzywej, w przeciwnym razie otrzymamy owal Cassiniego). Równanie biegunowe lemniskaty ma postać (to lemniskata jednostkowa, wszelkie inne są do niej geometrycznie podobne). Możemy za jego pomocą wyrazić element łuku krzywej jako
Zatem długość całkowita lemniskaty jest równa
Gauss najpierw zauważył, porównując liczby, a następnie udowodnił, że
Wielkość ta przypomina liczbę : też jest stosunkiem długości krzywej do jej „promienia”. Jak wskazuje to postać całki dającej długość łuku lemniskaty, mamy tu do czynienia z pierwiastkiem z wielomianu czwartego stopnia. Wiadomo, że całki pierwiastków z wielomianu drugiego stopnia dają się wyrazić przez funkcje elementarne. W przypadku wielomianów stopnia trzeciego i czwartego otrzymujemy tzw. całki eliptyczne: jest to nazwa wspólna, wywiedziona z zagadnienia obliczania długości łuku elipsy. Tak się składa, że całki dające długość łuku elipsy są całkami eliptycznymi drugiego rodzaju. Całkę pierwszego rodzaju spotkaliśmy w zagadnieniu wahadła matematycznego. Możemy także za Gaussem dowieść, że całkę eliptyczną zupełną pierwszego rodzaju
można wyrazić przez średnią arytmetyczno-geometryczną.
Mamy równość
Można ją wyrazić:
gdzie min, max oznaczają najmniejszą i największą wartość funkcji podcałkowej na przedziale , a nawiasy kątowe oznaczają uśrednienie funkcji po przedziale całkowania. Twierdzenie to zapewnia szybki algorytm do obliczania całek eliptycznych, w istocie tak szybki, że można go stosować także do innych celów: obliczania liczby
albo funkcji w rodzaju
, jeśli powiąże się je odpowiednio z całkami eliptycznymi.
Jeden z takich algorytmów, podany przez Eugene’a Salamina, korzysta z trzech ciągów zdefiniowane jak wyżej oraz
przy
. Otrzymuje się wówczas nierówność, którą spełnia
:
Daje to w kolejnych iteracjach:
0 : 2.914213562373095048801689 < π < 4.000000000000000000000000
1 : 3.140579250522168248311331 < π < 3.187672642712108627201930
2 : 3.141592646213542282149344 < π < 3.141680293297653293918070
3 : 3.141592653589793238279513 < π < 3.141592653895446496002915
4 : 3.141592653589793238462643 < π < 3.141592653589793238466361
Dla porównania oryginalny algorytm Archimedesa:
0 : 3.0000000 < π < 3.4641017
1 : 3.1058285 < π < 3.2153904
2 : 3.1326286 < π < 3.1596600
3 : 3.1393502 < π < 3.1460863
4 : 3.1410319 < π < 3.1427146
Widzimy, jak bardzo średnie arytmetyczno-geometryczne przyspieszają zbieżność, przy czym algorytm Salamina pochodzi z roku 1976 i od tamtej pory przedstawiono znacznie szybsze.
Na koniec pokażemy, czemu całka eliptyczna pierwszego rodzaju daje się obliczać w sposób odkryty przez Gaussa (dla ścisłości historycznej należy dodać, że nie tylko Gauss odkrył takie podejście, trochę wcześniej był Joseph Lagrange, choć zdaje się tylko Gauss zrozumiał dobrze od razu aspekt numeryczny sprawy).
Oznaczmy (
) następującą całkę:
Pokażemy, że nie zmienia się, kiedy przechodzimy do kolejnych kroków rekurencyjnych:
A skoro tak jest (szczegóły niżej), to kolejne wartości są takie same i przechodząc do granicy otrzymamy
Jako przykład pokażemy, jak procedura tego rodzaju pozwala obliczać okres wahadła matematycznego „niemal stającego na głowie” dla amplitud bliskich . Np. dla amplitudy
otrzymamy
Obliczamy średnią
:
0,008726535498374 1,00000000000000
0,504363267749187 0,093415927434105
0,298889597591646 0,217061195105173
0,257975396348409 0,254710292798989
0,256342844573699 0,256337645964924
0,256340245269312 0,256340245256133
Okres wahadła wydłuża się przy tak dużym wychyleniu razy. Euler w pracy E503, na którą powoływaliśmy się w poprzednim poście, także pokazuje rachunki dla takiego wychylenia, jednak jego wynik jest błędny.
(*) Nasze liczby wyjściowe można zapisać jako
dla pewnych wartości i
. Pierwszy związek rekurencyjny daje nam
Drugi związek rekurencyjny przyjmuje postać
Widać, że wzór ogólny będzie
Oba ciągi zbieżne są do granicy Związek z geometrią wielokątów wpisanych i opisanych na okręgu przedstawia rysunek. Zaczynając od sześciokąta, otrzymamy wartości początkowe przytoczone w tekście.
Wykażemy tożsamość (**). Topornym sposobem jej udowodnienia jest odpowiednia zamiana zmiennych pod całką. Wadą tego podejścia jest to, że podstawienie pojawia się jako deus ex machina i nam zostaje tylko sprawdzenie rachunków. Można też tożsamość przekształcić do postaci bardziej przydatnej w naszym przypadku
Funkcja podcałkowa po obu stronach jest odwrotnością pierwiastka, będziemy więc korzystać z rozwinięcia dwumianowego
Dla przyjmujemy z definicji
.
Rozwinięcie ma postać
Całka, która się tu pojawia, może być zapisana w postaci
Mamy więc dla rozwinięcie
Funkcję podpierwiastkową po lewej stronie tożsamości możemy zapisać przy użyciu tożsamości jako
Otrzymujemy wówczas
Dwójce przed odpowiada podwojenie przedziału całkowania:
. Rozwijamy oba pierwiastki w szereg:
Tylko wyrazy diagonalne przeżywają całkowanie, zostaje nam
Całkę z sinusa w potędze łatwo znaleźć zauważając, że
a także rozszerzając przedział całkowania do . W rozwinięciu dwumianowym
całkowanie przeżywają tylko wyrazy, w których wykładniki są przeciwne, stąd przytoczony wyżej wzór. Korzystałem ze sformułowania Petera Durena w świetnej książce Invitation to Classical Analysis.