Funkcja dzeta - problem

[Rysiu]

Mamy taką oto funkcję:



Dla argumentu 0.5+14i dość dokładna aczkolwiek także przybliżona wartość tej funkcji to:

0.0222411426099935892462132 - 0.1032581232664500579023631i

Spróbujmy sami to policzyć dla kilku pierwszych wyrazów według wzoru.

Dla czterech wyrazów (razem z jedynką) mamy więc:

1 + 1/(2^(0.5+14i)) + 1/(3^(0.5+14i)) + 1/(4^(0.5+14i)) =
= (1+0i) + (1+0i)/((2+0i)^(0.5+14i)) + (1+0i)/((3+0i)^(0.5+14i)) + (1+0i)/((4+0i)^(0.5+14i)) =
= (1+0i) + (1+0i)/(-1.3594-0.3898i) + (1+0i)/(-1.6400+0.5569i) + (1+0i)/(1.6960+1.0599i) =
= (1+0i) - 0.6797+0.1949i - 0.5467-0.1856i + 0.4240-0.2649i =
= 0.3203-0.1949i - 0.5467-0.1856i + 0.4240-0.2649i =
= -0.2264-0.0093i + 0.4240-0.2649i =
= 0.1976-0.2742i

Dla pięciu wyrazów (suma czterech wcześniejszych plus piąty):

0.1976-0.2742i + 1/(5^(0.5+14i)) =
= 0.1976-0.2742i + (1+0i)/(-1.9167-1.1515i) =
= 0.1976-0.2742i - 0.3833+0.2303i =
= - 0.1857-0.5045i

Dla sześciu wyrazów:

- 0.1857-0.5045i + 1/(6^(0.5+14i)) =
= - 0.1857-0.5045i + (1+0i)/(6^(0.5+14i)) =
= - 0.1857-0.5045i + (1+0i)/(2.4466-0.1177i) =
= - 0.1857-0.5045i + 0.4077+0.0196i =
= 0.222-0.4849i

Dla siedmiu wyrazów:

0.222-0.4849i + 1/(7^(0.5+14i)) =
= 0.222-0.4849i + (1+0i)/(7^(0.5+14i)) =
= 0.222-0.4849i + (1+0i)/(-1.3584+2.2703i) =
= 0.222-0.4849i - 0.1940-0.3243i =
= 0.028-0.1606i

Dla ośmiu wyrazów:

0.028-0.1606i + 1/(8^(0.5+14i)) =
= 0.028-0.1606i + (1+0i)/(8^(0.5+14i)) =
= 0.028-0.1606i + (1+0i)/(-1.8923-2.1021i)
= 0.028-0.1606i - 0.2365+0.2627i =
= -0.2085-0.4233i

Dla dziewięciu wyrazów:

- 0.2085-0.4233i + 1/(9^(0.5+14i)) =
= - 0.2085-0.4233i + (1+0i)/(9^(0.5+14i)) =
= - 0.2085-0.4233i + (1+0i)/(2.3795-1.8269i) =
= - 0.2085-0.4233i + 0.2644+0.2029i =
= 0.0559-0.2204i

Dla dziesięciu wyrazów:

0.0559-0.2204i + 1/(10^(0.5+14i)) =
= 0.0559-0.2204i + (1+0i)/(10^(0.5+14i)) =
= 0.0559-0.2204i + (1+0i)/(2.1567+2.3126i) =
= 0.0559-0.2204i + 0.2156-0.2312i =
= 0.2715-0.4516i

Dla jedenastu wyrazów:

0.2715-0.4516i + 1/(11^(0.5+14i)) =
= 0.2715-0.4516i + (1+0i)/(11^(0.5+14i)) =
= 0.2715-0.4516i + (1+0i)/(-1.8281+2.7672i) =
= 0.2715-0.4516i - 0.1662-0.2515i =
= 0.1053-0.2001i

Suma dla pierwszych dwudzistu wyrazów wynosi: -0.2044+0.3320i, a dla dwustu już -0.9251+0.2422i. Jak widać wzrost ilości iteracji objawia się zwiększeniem błędu zamiast jego zmniejszeniem. Pytanie dlaczego tak jest? Gdzie mam błąd? Muszę jeszcze dodać, że akurat w tym przypadku błąd się zwiększa. Czasem wyniki są od razu złe, a czasem dobre (przy uwzględnieniu, że są one przybliżone). Jest to jednak chyba kwestią przypadku.

Plik z kodem źródłowym można podejrzeć tutaj: http://code.google.com/p/olib/source/browse/trunk/samples/mathematics/riemann_zeta.cpp. Jak łatwo zauważyć operacje na liczbach zespolonych sam implementowałem ale dobrze je przetestowałem i wygląda na to, że działają ok.

Co w tym wszystkim nie gra?

[stiven]

Może dlatego, że to nie funkcja a szereg i to tak nie działa, że dla liczb zespolonych, które tam wstawisz policzysz to w ten sposób. W jednym zdaniu nie da się tego opisać. Po prostu nie tędy droga zupełnie.

I zasadniczo z tego co ja wiem to aby ten szereg był zbieżny to część rzeczywista musi być większa od... 1. 

[Karlik]

Nawet dla rzeczywistych wartości suma szeregu nie jest oczywista i nie tak łatwo ją policzyć.
Przykładowo: szereg f(x) = 1-1+1-1+1-1+1....
Jest zbieżny? Ma sumę?
Suma tego szeregu to bodajże 1/2, ale licząc to Twoją metodą nigdy takiego wyniku nie osiągniesz.

[stiven]

Zdecydowanie musi być większa od 1. Zrobiłem symulacje w Originie i upewniłem się, że dla zespolonych jest tak samo jak dla rzeczywistych.

Tak więc Rysiu błąd w założeniu.

Zerknij sobie na wykresy część_rzeczywista_wyniku=f(liczba_iteracji) dla:

Co więcej widać wyraźnie, że im dalej zabrniesz przy błędnym założeniu tym głupszy wynik dostaniesz.

[Rysiu]

I zasadniczo z tego co ja wiem to aby ten szereg był zbieżny to część rzeczywista musi być większa od... 1.

Zgadza się. Właśnie sprawdziłem programem i dla części rzeczywistej większej od 1 wynik wychodzi z wielką precyzją dokładny. Potem jeszcze dokładniej to sprawdzę.

Tylko co jeżeli ja nie mogę tak się ograniczać?

[stiven]

Dorzuciłem obrazki zarówno dla części rzeczywistej jak i urojonej.

Co zrobić? Zmienić metodę. Szukam właśnie w originie czy ma funkcję aby od razu liczyć dzetę ale nie potrafię znaleźć. Nigdy mi to nie było potrzebne, szczerze to nawet liczb zespolonych rzadko używam i musiałem się chwilę dokształcić. Myślę, że są narzędzia do liczenia tego. Sprawdzałeś:

Kod [Zaznacz] Rozwiń

http://mathworld.wolfram.com/?

Skąd wziąłeś wynik, który podałeś? Mam dziwne wrażenie  ;) ;) ;), że skoro szereg nie jest zbieżny to raczej trudno oczekiwać czegoś na kształt "wyniku".

[Rysiu]

Wykorzystałem bibliotekę mpmath dla Pythona:

http://mpmath.googlecode.com/svn/trunk/doc/build/functions/zeta.html

Można ją łatwo zainstalować na Linuxie (jeżeli jest dostępna w repo):

Kod [Zaznacz] Rozwiń

apt-get install python-mpmath

Potem już według wcześniej podanej strony wystarczy tylko np.:

Kod [Zaznacz] Rozwiń

>>> from mpmath import *
>>> mp.dps = 25; mp.pretty = True
>>> zeta(-3+4j)
(-0.03373057338827757067584698 + 0.2774499251557093745297677j)

Wartości dla części rzeczywistej <1 muszą dać się policzyć. Funkcja dzeta ma miejsca zerowe. Pierwszym z nich jest 0.5+14.134725141734693790457251i i dla niego mpmath zwraca (-1.048365080558823465032955e-16 + 6.585259277605158331268578e-16j) - czyli prawdłowo.

[Karlik]

Rysiu: skoro masz gotową implementację (do tego z otwartym kodem źródłowym) to nie prościej po prostu zobaczyć jak to zostało zaimplementowane?

[stiven]

Rysiu: przecież to bardzo proste. Gamma Eulera/Borwein/Hurwitz. Przychylam się do zdania Karlika, że skoro masz gotowy program to zobacz jak on to liczy... albo trochę doczytaj o zagadnieniu bo nawet wikipedia zdaje się odpowiadać na Twoje pytania.

[stiven]

Pozwolę sobie odgrzebać temat.
Problemem zajął się zespół SZTAKI Desktop Grid i ruszył nowy podprojekt Zeta-Search.
Szczegóły tu:

Kod [Zaznacz] Rozwiń

http://riemann-siegel.com/
Aby dostawać próbki trzeba mieć zaptaszone w ustawieniach: "Riemann Zeta Research Project"

Piszą, że karmić będzie w tym tygodniu.