Egzotyczna hiperkula - czyli o pomiarach w przestrzeni wielowymiarowej

"Jak oczami wyobraźni zobaczyć 4 wymiary? - zapytano matematyka.
To proste - odpowiedział - wystarczy wyobrazić sobie n-wymiarów i podstawić n=4"
🙂

Hipersześcian i Hiperkula - rzut

Dzisiejszy wpis poświęcę pomiarom odległości, powierzchni i pojemności w przestrzeniach wielowymiarowych. N-wymiarowa przestrzeń euklidesowa dostarcza dosyć oczywistą metrykę - a przez to wydawałoby się - bardzo intuicyjną. To wrażanie jest jednak mylne, co łatwo pokazać analizując wpływ zwiększania liczby wymiarów na dokonywane pomiary. Jak w zależności od liczby wymiarów zmienia się powierzchnia i objętość kuli? Analogicznie - jak zmienia się maksymalna odległość pomiędzy wierzchołkami kostki? Obiecuję - odpowiedzi będą zaskakujące 🙂

Możesz mieć wrażenie, że to wyłącznie abstrakcyjne rozważania. Czy na pewno? Ja w zasadzie na co dzień analizuję Klientów opisanych szeregiem miar. Poszukiwanie podobieństw, skupień, segmentów czy "najbliższych sąsiadów" niemal w całości opiera się na wielowymiarowej metryce euklidesowej. Zapraszam do pogłębienia wiedzy w tym obszarze:-) Zapewniam - warto!

Przestrzeń kartezjańska \mathbb{R}^n

Przestrzeń euklidesowa R^n

Przestrzeń kartezjańska \mathbb{R}^n to przestrzeń współrzędnych rzeczywistych

(x_1, x_2, \ldots, x_n)\in\mathbb{R}^n

z określoną metryką

d(\mathbb{x},\mathbb{y})=\sqrt{\sum_{i=1}^n(x_i-y_i)^2}

gdzie

\mathbb{x}=(x_1, x_2, \ldots, x_n)\in\mathbb{R}^n

\mathbb{y}=(y_1, y_2, \ldots, y_n)\in\mathbb{R}^n

Można powiedzieć, że metryka euklidesowa to pewna forma twierdzenia Pitagorasa.

Hipersześcian

Dla uproszczenia zdefiniujemy hipersześcian o długości boku a, którego jeden z wierzchołków to punkt 0, a inne wierzchołki składają się ze współrzędnych nieujemnych.

Hipersześcian o boku długości a to zbiór punktów

\mathbb{x}=(x_1, x_2, \ldots, x_n)\in\mathbb{R}^n

Spełniających zależność

\begin{cases}0\leq x_1\leq a\\0\leq x_2\leq a\\ \cdots\\0\leq x_n\leq a\end{cases}

Hipersześcian - wymiar od 1 do 4

Odcinek powstaje z "przesunięcia" punktu w dodatkowym nowym wymiarze. Kwadrat powstaje z odcinka poprzez "przesunięcie" odcinka w dodatkowym nowym wymiarze. Kostka powstaje z kwadratu poprzez przesunięcie kwadratu w dodatkowym nowym wymiarze. W tym postępowaniu uwidacznia się zależność rekurencyjna.

Hiperkula

Hiperkula to uogólnienie pojęcia kuli na n-wymiarów przy rozważaniu punktów w przestrzeni euklidesowej \mathbb{R}^n. Ponownie, dla uproszczenia, założymy, że środkiem kuli jest punkt 0.

Hiperkula o promieniu r to zbiór punktów

\mathbb{x}=(x_1, x_2, \ldots, x_n)\in\mathbb{R}^n

Spełniających zależność

d(0,\mathbb{x})=\sqrt{\sum_{i=1}^n x_i^2}\leq r

Hiperkula - wymiar od 1 do 4

Odcinek powstaje z "przesunięcia" punktu w dodatkowym nowym wymiarze. Koło powstaje z odcinka poprzez "przesunięcie" odcinka (jednocześnie odcinek zmienia długość) w dodatkowym nowym wymiarze. Kula powstaje z koła poprzez przesunięcie koła (jednocześnie koło zmienia promień) w dodatkowym nowym wymiarze. Ponownie w tym postępowaniu uwidacznia się zależność rekurencyjna, którą niebawem wykorzystamy.

Przydatna całka \displaystyle\int_{-r}^r\Big(r^2-x^2\Big)^y dx

Wszelkie niezbędne wzory będę wyznaczał korzystając z elementarnych całek. W okręgach, kołach, kulach często pojawiają się równania postaci (r^2-x^2)^y. Poniżej wyznaczę ogólny i pomocny wzór na całkę oznaczoną:

\displaystyle\int_{-r}^r\Big(r^2-x^2\Big)^y dx=\displaystyle\int_{-r}^r\Bigg[r^2\bigg(1-\Big(\frac{x}{r}\Big)^2\bigg)\Bigg]^y dx=

=r^{2y}\displaystyle\int_{-r}^r\bigg[1-\Big(\frac{x}{r}\Big)^2\bigg]^y dx=

={\scriptsize\begin{bmatrix}\frac{x}{r}=\sin t && t = \arcsin\frac{x}{r}\\x=r\sin t && t_{-r}=\arcsin\frac{-r}{r}=\arcsin(-1)=-\frac{\pi}{2}\\dx=r\cos t dt && t_r=\arcsin\frac{r}{r}=\arcsin(1)=\frac{\pi}{2}\end{bmatrix}}=

=r^{2y}\displaystyle\int_{-\frac{\pi}{2}}^\frac{\pi}{2}\Big(1-\sin^2t\Big)^y r\cos t dt=

=r^{2y+1}\displaystyle\int_{-\frac{\pi}{2}}^\frac{\pi}{2}\cos^{2y}t\cos t dt=

=r^{2y+1}\displaystyle\int_{-\frac{\pi}{2}}^\frac{\pi}{2}\cos^{2y+1}t dt

Teraz skorzystamy z poniższej formuły redukcyjnej

\displaystyle\int\cos^nxdx=\frac{1}{n}\cos^{n-1}x\sin x+\frac{n-1}{n}\displaystyle\int\cos^{n-2}xdx

Z wyprowadzeniem powyższego wzoru można zapoznać się np. w filmie "Calculus - Reduction Formula for Powers of Cosine".

Podstawiając otrzymujemy

\displaystyle\int_{-r}^r\Big(r^2-x^2\Big)^y dx=

=\frac{r^{2y+1}}{2y+1}\Bigg(\Big[\cos^{2y}t\sin t\Big]_{-\frac{\pi}{2}}^\frac{\pi}{2}+2y\displaystyle\int_{-\frac{\pi}{2}}^\frac{\pi}{2}\cos^{2y-1}t dt\Bigg)=

={\scriptsize\frac{r^{2y+1}}{2y+1}\Bigg(\cos^{2y}\frac{\pi}{2}\sin\frac{\pi}{2}-\cos^{2y}\big(-\frac{\pi}{2}\big)\sin\big(-\frac{\pi}{2}\big)+2y\displaystyle\int_{-\frac{\pi}{2}}^\frac{\pi}{2}\cos^{2y-1}t dt\Bigg)}=

=\frac{r^{2y+1}}{2y+1}\Bigg(0\cdot 1-0\cdot\big(-1\big)+2y\displaystyle\int_{-\frac{\pi}{2}}^\frac{\pi}{2}\cos^{2y-1}t dt\Bigg)=

=\frac{2y}{2y+1}r^{2y+1}\displaystyle\int_{-\frac{\pi}{2}}^\frac{\pi}{2}\cos^{2y-1}t dt

Ostatecznie

{\small\displaystyle\int_{-r}^r\Big(r^2-x^2\Big)^y dx=\frac{2y}{2y+1}r^{2y+1}\displaystyle\int_{-\frac{\pi}{2}}^\frac{\pi}{2}\cos^{2y-1}t dt}

Jest to wzór warty zapamiętania 🙂

Uogólnienie wzoru na objętość kuli

W uogólnieniu wzoru na objętość kuli n-wymiarowej pomaga zrozumienie zależności rekurencyjnej (pisałem o tym definiując hipersześcian i hiperkulę), dosyć oczywistej gdy się rozważa postać i objętość kuli krok po kroku, tzn. zaczynając od wymiaru 1, przechodząc przez wymiary 2 i 3.

Objętość kuli o wymiarze n = 1

Można powiedzieć, że kula o wymiarze 1 (czyli odcinek) powstaje z przesunięcia punktu (który ma wymiar 0), wzdłuż nowego wymiaru X w zakresie x\in[-R,R]

Kula - wymiar 1

{\Large V_R^{\dim 1}=2R}

Objętość kuli o wymiarze n = 2

Mając odcinek (kulę o wymiarze 1), przesuwając go wzdłuż nowego wymiaru X w zakresie x\in[-R,R], dostosowując jednocześnie jego długość (czyli promień kuli o wymiarze 1), otrzymujemy koło - tzn. kulę o wymiarze 2.

Kula - wymiar 2

V_R^{\dim 2}=\displaystyle\int_{-R}^R V_{r(x)}^{\dim 1}dx=

=\displaystyle\int_{-R}^R 2r(x)dx=2\displaystyle\int_{-R}^R \Big(R^2-x^2\Big)^\frac{1}{2}dx=...

Stosujemy wcześniej wyprowadzony wzór dla y=\frac{1}{2}

...=2\frac{2\cdot\frac{1}{2}}{2\cdot\frac{1}{2}+1}R^{2\cdot\frac{1}{2}+1}\displaystyle\int_{-\frac{\pi}{2}}^\frac{\pi}{2}\cos^{2\cdot\frac{1}{2}-1}t dt=

=2\frac{1}{2}R^2\displaystyle\int_{-\frac{\pi}{2}}^\frac{\pi}{2}\cos^0t dt=

=R^2\displaystyle\int_{-\frac{\pi}{2}}^\frac{\pi}{2}dt=R^2\Big[t\Big]_{-\frac{\pi}{2}}^\frac{\pi}{2}=

=R^2\Big(\frac{\pi}{2}+\frac{\pi}{2}\Big)=\pi R^2

{\Large V_R^{\dim 2}=\pi R^2}

Objętość kuli o wymiarze n = 3

Procedura jest analogiczna. Przesuwając koło (czyli kulę o wymiarze 2) wzdłuż nowego wymiaru X w zakresie x\in[-R,R], jednocześnie zmieniając odpowiednio promień koła, otrzymujemy kulę o wymiarze 3.

Kula - wymiar 3

V_R^{\dim 3}=\displaystyle\int_{-R}^R V_{r(x)}^{\dim 2}dx=

=\displaystyle\int_{-R}^R \pi r^2(x)dx=\pi\displaystyle\int_{-R}^R \Big(\sqrt{R^2-x^2}\Big)^2dx=

=\pi\displaystyle\int_{-R}^R \big(R^2-x^2\big)dx=...

Jak widać jest to całka zupełnie elementarna, nie ma konieczności stosowania wcześniej wyprowadzonego wzoru. Ja to jednak uczynię aby potwierdzić, że poprzednio otrzymana formuła jest prawdziwa w każdym przypadku 🙂

Stosując elementarną całkę

...=\pi\Bigg(\Big[R^2x\Big]_{-R}^R-\Big[\frac{x^3}{3}\Big]_{-R}^R\Bigg)

=\pi\Bigg(\Big[R^2R-R^2(-R)\Big]-\frac{R^3-(-R^3)}{3}\Bigg)=

\pi\bigg(2R^3-\frac{2R^3}{3}\bigg)=\pi\bigg(\frac{6R^3}{3}-\frac{2R^3}{3}\bigg)=\pi\frac{4R^3}{3}

{\Large V_R^{\dim 3}=\frac{4}{3}\pi R^3}

Stosując wcześniej wyprowadzony wzór dla y=1

...=\pi\displaystyle\int_{-R}^R \big(R^2-x^2\big)^1dx=

=\pi \frac{2\cdot 1}{2\cdot 1+1}R^{2\cdot 1+1}\displaystyle\int_{-\frac{\pi}{2}}^\frac{\pi}{2}\cos^{2\cdot 1-1}t dt=

=\pi \frac{2}{3}R^3\displaystyle\int_{-\frac{\pi}{2}}^\frac{\pi}{2}\cos t dt=\pi \frac{2}{3}R^3\Big[\sin t\Big]_{-\frac{\pi}{2}}^\frac{\pi}{2}=

=\pi \frac{2}{3}R^3\Big(\sin\frac{\pi}{2}-\sin(-\frac{\pi}{2})\Big)=\pi \frac{2}{3}R^3\Big(1-(-1)\Big)=

=\pi \frac{2}{3}R^3\big(2\big)=\frac{4}{3}\pi R^3

{\Large V_R^{\dim 3}=\frac{4}{3}\pi R^3}

Objętość kuli o wymiarze n - uogólnienie

Na bazie przypadków 1, 2, 3 wymiarowych ujawniła się procedura rekurencyjna pozwalająca tworzyć hiperkulę wymiaru n z hiperkuli wymiary n-1. Zwyczajnie "przesuwamy" hiperkulę wymiaru n-1 w nowym wymiarze X w zakresie x\in[-R,R], dostosowując jej promień r(x)=\sqrt{R^2-x^2}, gdzie R jest niezmiennym promieniem hiperkuli.

{\small V_R^{\dim n}=\begin{cases}2R&\text{dla}\quad n=1\\ \displaystyle\int_{-R}^R V_{r(x)=\sqrt{R^2-x^2}}^{\dim n-1}dx&\text{dla}\quad n>1\end{cases}}

Objętość kuli o wymiarze n = 4

V_R^{\dim 4}=\displaystyle\int_{-R}^R V_{r(x)}^{\dim 3}dx=

=\displaystyle\int_{-R}^R\frac{4}{3}\pi r^3(x)dx=\frac{4}{3}\pi\displaystyle\int_{-R}^R \Big(R^2-x^2\Big)^\frac{3}{2}dx=...

Stosujemy wzór dla y=\frac{3}{2}

...=\frac{4}{3}\pi\frac{2\cdot\frac{3}{2}}{2\cdot\frac{3}{2}+1}R^{2\cdot\frac{3}{2}+1}\displaystyle\int_{-\frac{\pi}{2}}^\frac{\pi}{2}\cos^{2\cdot\frac{3}{2}-1}t dt=

=\frac{4}{3}\pi\frac{3}{4}R^4\displaystyle\int_{-\frac{\pi}{2}}^\frac{\pi}{2}\cos^2t dt=

=\pi R^4\displaystyle\int_{-\frac{\pi}{2}}^\frac{\pi}{2}\cos^2t dt=...

Całkę \displaystyle\int_{-\frac{\pi}{2}}^\frac{\pi}{2}\cos^2t dt możemy oczywiście liczyć, ale łatwiej i szybciej zauważyć, że

\displaystyle\int_{-\frac{\pi}{2}}^\frac{\pi}{2}\cos^2t+\sin^2tdt=\displaystyle\int_{-\frac{\pi}{2}}^\frac{\pi}{2}dt=\pi

i że połowa tego "prostokąta" przypada dla \displaystyle\int_{-\frac{\pi}{2}}^\frac{\pi}{2}\cos^2t dt, zatem

\displaystyle\int_{-\frac{\pi}{2}}^\frac{\pi}{2}\cos^2t dt=\frac{\pi}{2}

Podstawiając do

...=\pi R^4\displaystyle\int_{-\frac{\pi}{2}}^\frac{\pi}{2}\cos^2t dt=\pi R^4\frac{\pi}{2}=\frac{\pi^2}{2}R^4

{\Large V_R^{\dim 4}=\frac{\pi^2}{2}R^4}

Hiperkula wpisana w hipersześcian

Wyznaczyliśmy, że:

  • dla 1 wymiaru V_R^{\dim 1}=2R
  • dla 2 wymiarów V_R^{\dim 2}=\pi R^2
  • dla 3 wymiarów V_R^{\dim 3}=\frac{4}{3}\pi R^3
  • dla 4 wymiarów V_R^{\dim 4}=\frac{\pi^2}{2}R^4

5-ty i większy wymiar można wyznaczyć analogicznie, poniżej podam gotowe wzory opublikowane na Wikipedii.

  • dla 5 wymiarów V_R^{\dim 5}=\frac{8}{15}\pi^2R^5
  • dla 6 wymiarów V_R^{\dim 6}=\frac{\pi^3}{6}R^6
  • dla 7 wymiarów V_R^{\dim 7}=\frac{16}{105}\pi^3R^7
  • dla 8 wymiarów V_R^{\dim 8}=\frac{\pi^4}{24}\pi^8R^5
  • dla 9 wymiarów V_R^{\dim 9}=\frac{32}{945}\pi^4R^9
  • dla 10 wymiarów V_R^{\dim 10}=\frac{\pi^5}{120}\pi^2R^{10}
  • ...
  • dla n wymiarów V_R^{\dim n}=\frac{\pi^\frac{n}{2}}{\Gamma(\frac{n}{2}+1)}R^n

Jeśli krawędź n-wymiarowego hipersześcianu ma długość a, to jego objętość wynosi a^n. Promień hiperkuli wpisanej w hipersześcian będzie stanowił połowę długości krawędzi hipersześcianu, tzn. R=\frac{a}{2}.

Zadziwiające, objętość kuli wpisanej w hipersześcian "znika" wraz ze wzrostem liczby wymiarów, i ten zanik postępuje na prawdę szybko!!

Aby wyjaśnić ten fenomen należy zauważyć, że wraz ze wzrostem liczby wymiarów rośnie odległość "po przekątnej" do wymiarów, przy jednoczesnym zachowaniu odległości wzdłuż wymiarów. Maksymalna długość "przekątnej" w hipersześcianie jednostkowym (o krawędzi długości 1) będzie równa \sqrt{n}, gdzie n oznacza liczbę wymiarów. Wynika to bezpośrednio z zaaplikowania metryki do wierzchołków \mathbb{0}=(0,0,\ldots,0) i \mathbb{1}=(1,1,\ldots,1).

W konsekwencji objętość hipersześcianu kumuluje się w okolicach jego wierzchołków, zaś kula wpisana leży w jego środku stykając się ze wszystkimi jego ścianami.

Objętość hiperkuli jednostkowej

Kolejna ciekawostka kryje się w objętości hiperkuli jednostkowej (o promieniu 1). Okazuje się, że ta objętość jest najwyższa dla 5 wymiarowej przestrzeni.

Objętość hiperkuli jednostkowej - zależność od liczby wymiarów

Do liczby wymiarów n=5 objętość hiperkuli jednostkowej rośnie, w 5-tym wymiarze osiąga maksimum, powyżej 5-tego wymiaru objętość systematycznie spada, w granicy osiągając wartość 0.

Ciekawostka: jeśli rozważyć powierzchnię hiperkuli n-wymiarowej (czyli hipersfery n-1 wymiarowej) to powierzchnia będzie najwyższa dla hipersfery 6-wymiarowej, czyli hiperkuli 7-wymiarowej. Powyżej 7 wymiaru powierzchnia zaczyna spadać z granicą w nieskończoności równą 0.  Powierzchnię hiperkuli można wyznaczyć poprzez pochodną objętości względem promienia.

Na koniec

Przestrzenie wielowymiarowe, nawet te z bardzo intuicyjną metryką, skrywają wiele tajemnic. Dziś odkryliśmy tylko kilka 🙂 Mam nadzieję, że Wam się podobało!

Pozdrowienia,

Mariusz Gromada

Spirala Ulama (spirala liczb pierwszych) - część 2 - czyli funkcja kwadratowa i zabawy z rekurencją (część 5)

W części 1 wpisu na temat Spirali Ulama zaznaczyłem, że efekt wizualnego ułożenia liczb pierwszych na diagonalach spirali kwadratowej jest konsekwencją głównie dwóch własności:

  1. Na przekątnych są albo same liczby parzyste, albo same liczby nieparzyste, zatem tylko diagonale z liczbami nieparzystymi będą agregować liczby pierwsze;
  2. Niektóre diagonale zagęszczają bardziej liczby pierwsze niż inne, co wynika z zależności pomiędzy przekątnymi i funkcją kwadratową oraz faktem, że niektóre funkcje kwadratowe generują więcej liczb pierwszych niż inne.

Dzisiejszy tekst poświęcę przybliżeniu własności nr 2.

Wielomiany i rekurencja

Funkcja kwadratowa i rekurencja

Dla wielomianów możemy zawsze podać ich postać rekurencyjną, Jest to własność mało znana, jednak dosyć prosta w uzasadnieniu. Pokażę to na przykładzie funkcji kwadratowej, jednocześnie wzbogacając cykl "Zabawy z rekurencją" 🙂

f(x) = ax^2+bx+c

Rozważmy następnie równanie

f(x+1)=f(x)+Bx+C

Podstawiając i upraszczając...

a(x+1)^2+b(x+1)+c=ax^2+bx+c+Bx+C

ax^2+2ax+a+bx+b+c=ax^2+bx+c+Bx+C

2ax+a+b=Bx+C

2a=B    oraz    a+b=C

otrzymujemy

a=\frac{B}{2}    oraz    b=C-a

Następnie analizując f(1) mamy

f(1)=a+b+c    zatem    c=f(1)-a-b

c=f(1)-a-(C-a)=f(1)-C

Wniosek: jeśli znana jest relacja rekurencyjna f(x+1)=f(x)+Bx+C oraz znamy wartość f(1) to jesteśmy w stanie jednoznacznie wskazać równanie kwadratowe ax^2+bx+c spełniające daną zależność rekurencyjna, gdzie

a=\frac{B}{2},    b=C-a,    c=f(1)-C

Uogólnienia dokonujemy na bazie dwumianu Newtona

(x+1)^n = {n\choose 0}x^n+{n\choose 1}x^{n-1}+{n\choose 2}x^{n-2}+\ldots+{n\choose {n-1}}x+{n\choose n}

Funkcja kwadratowa i linie proste / przekątne na spirali Ulama

Poniżej spróbuję pokazać w jaki sposób "nawiajnie prostej na kwadrat" sprawia, że parabole w efekcie otrzymują kształt linii prostych.

Przykład 1 - Pionowa prosta

Spirala Ulama - równanie prostej 1

Zapisujemy zależność rekurencyjną

f(1)=4

f(2)=f(1)+1+2+3+3+2

f(3)=f(2)+2+3+5+5+3

\ldots

f(n+1)=f(n)+n+2n+(2n+1)+(2n+1)+(n+1)

f(n+1)=f(n)+8n+3

Teraz, korzystając z wyprowadzonego wcześniej wzoru, wyznaczamy współczynniki a, b, c funkcji kwadratowej spełniającej f(n+1)=f(n)+8n+3.

a=\frac{8}{2}=4,    b=3-4=-1,    c=4-3=1

f(n)=4n^2-n+1

Dla testu czy wszystko jest ok sami podstawcie n=1,2,3\ldots

Przykład 2 - Przękątna (linia diagonalna)

Spirala Ulama - równanie prostej 2

Ponownie zapisujemy zależność rekurencyjną, tym razem nieco prostszą.

f(1)=3

f(2)=f(1)+2+2+3+3

f(3)=f(2)+4+4+5+5

\ldots

f(n+1)=f(n)+n+2n+2n+(2n+1)+(2n+1)

f(n+1)=f(n)+8n+2

Korzystając ze znanego wzoru wyznaczamy a, b, c dla f(n+1)=f(n)+8n+2.

a=\frac{8}{2}=4,    b=2-4=-2,    c=3-2=1

f(n)=4n^2-2n+1

W kolejnej części cyklu "o Spirali Ulama" podam wzory kilku przekątnych (czyli funkcji kwadratowych) generujących nienaturalnie dużo liczb pierwszych 🙂

Pozdrowienia,

Mariusz Gromada

Rekurencja pośrednia - czyli zabawy z rekurencją (część 4)

W pierwszych trzech częściach "Zabaw z rekurencją" skupialiśmy się na rekurencji bezpośredniej, tzn. na sytuacji, kiedy w ciele funkcji dochodzi do wywołania "siebie samej". Przebieg rekurencji bezpośredniej jest dość oczywisty, struktura wywołania, argumenty, jak też warunek stopu, są takie same dla wszystkich odwołań.

Rekurencja pośrednia

O rekurencji pośredniej mówimy w sytuacji"łańcucha wywołań". Przykładowo funkcja f(.) wywołuje funkcję g(.), następnie funkcja g(.) wywołuje f(.), zatem ponowne wywołanie funkcji f(.) realizowane jest bezpośrednio przez funkcję g(.), jednak pośrednio przez f(.), gdyż to f(.) wywołała g(.).

Typy rekurencji

Długość łańcucha nie musi być ograniczona, w rzeczywistości wywołania pośrednie mogą mieć nietrywialną strukturę, mogą "cofać się" do przednich elementów, "iść na skróty", "rozdzielać się", a w szczególności może dochodzić do wariantów mieszanych - tzn. wywołań bezpośrednich i pośrednich (różnego typu) w ramach jednej procedury. Dobrze to obrazuje poniższy schemat.

Rekurencja pośrednia

Aproksymacja funkcji sin(x) oraz cos(x) przy wykorzystaniu połączenia rekurencji bezpośredniej i rekurencji pośredniej

Przypomnijmy podstawowe tożsamości trygonometryczne dla wielokrotności kątów.

\sin(2x)=2\sin(x)\cos(x)

\cos(2x)=\cos^2(x)-\sin^2(x)

Równoważnie powyższe można zapisać jako

\sin(x)=2\sin\big(\frac{x}{2}\big)\cos\big(\frac{x}{2}\big)

\cos(x)=\cos^2\big(\frac{x}{2}\big)-\sin^2\big(\frac{x}{2}\big)

Zwróćmy uwagę, że znając rozwiązanie dla argumentu mniejszego \frac{x}{2} możemy podać rozwiązanie dla x – zatem tożsamości trygonometryczne są w istocie rekurencją z odwołaniami bezpośrednimi i pośrednimi! Funkcje \sin(x) w otoczeniu 0 można przybliżyć przez x, natomiast funkcję \cos(x) przez stałą wartość 1. Im mniejsze otoczenie 0 wybierzemy tym lepsza aproksymacja w zadanym przedziale, a w konsekwencji mniejszy błąd oszacowania w całości. Przyjęte wartości w otoczeniu 0 dają również pewny warunek stopu! Mamy więc wszystko co niezbędne do zastosowania strategii rekurencyjnej w aproksymacji.

Ustalmy stałą a>0 (reprezentującą otoczenie 0), następnie definiujemy dwie funkcje rekurencyjne

\text{s}(x)=\begin{cases}x&\text{dla}\quad |x|<a\\2\text{s}\big(\frac{x}{2}\big)\text{c}\big(\frac{x}{2}\big)&\text{dla}\quad |x|\geq a\end{cases}

\text{c}(x)=\begin{cases}1&\text{dla}\quad |x|<a\\\text{c}^2\big(\frac{x}{2}\big)-\text{s}^2\big(\frac{x}{2}\big)&\text{dla}\quad |x|\geq a\end{cases}

Podkreślmy ponownie, że funkcja \text{s}(x) wywołuje siebie bezpośrednio oraz wskazuje na funkcję \text{c}(x), która, oprócz bezpośredniego wywołania siebie samej, wskazuje ponownie na \text{s}(x). Jest to zatem ciekawa kombinacja rekurencji bezpośredniej z rekurencją pośrednią. Zapiszmy to w mXparser.

/* Definicja funkcji rekurencyjncyh */
Constant a = new Constant("a", 0.1);
Function s = new Function("s(x) =  if( abs(x) < a, x, 2*s(x/2)*c(x/2) )", a);
Function c = new Function("c(x) =  if( abs(x) < a, 1, c(x/2)^2-s(x/2)^2 )", a);

/* Wskazanie, ze 's' korzysta z 'c', a 'c' korzysta z 's' */
s.addDefinitions(c);
c.addDefinitions(s);

Oczekujemy, że im mniejszy parametr a>0 tym lepsza aproksymacja funkcji \sin(x) oraz \cos(x) przez odpowiednio \text{s}(x) oraz \text{c}(x). Poniżej wykresy dla a=0.5 oraz a=0.01.


/* Dane do wykresu s(x) vs sin(x) */
for (double x = -MathConstants.PI; x <= MathConstants.PI; x=x+0.02)
mXparser.consolePrintln("[ " + x +", " + MathFunctions.sin(x) + ", " + s.calculate(x) + " ],");

/* Dane do wykresu c(x) vs cos(x) */
for (double x = -MathConstants.PI; x <= MathConstants.PI; x=x+0.02)
mXparser.consolePrintln("[ " + x +", " + MathFunctions.cos(x) + ", " + c.calculate(x) + " ],");

Wniosek - proste zapisy rekurencyjne dają złożone wyniki! 🙂

Pozdrowienia,

Mariusz Gromada

Pamiętaj o:

import org.mariuszgromada.math.mxparser.*;
import org.mariuszgromada.math.mxparser.mathcollection.*;

Kod:

Zobacz również:

  1. Polowanie na czarownice - czyli zabawy z rekurencją (część 1)
  2. Prędkość ucieczki do nieskończoności - czyli zabawy z rekurencją (część 2)
  3. Naiwny test pierwszości - czyli zabawy z rekurencją (część 3)

Naiwny test pierwszości - czyli zabawy z rekurencją (część 3)

Liczby pierwsze

Jednym z najprostszych testów pierwszości jest weryfikacja czy dana liczba n posiada dzielnik z przedziału (2, \sqrt{n}) - takie podejście nazywane jest metodą naiwną - i niestety charakteryzuje się dużą złożonością obliczeniową. Nawet przy wykorzystaniu Sita Eratostenesa złożoność obliczeniowa sięga \frac{\sqrt{n}}{\log{n}}. Jednak w cyklu "Zabawy z rekurencją" nie bardzo zwracamy uwagę na złożoność 🙂 , bardziej chodzi o zobrazowanie jak całe algorytmy mogą być łatwo zapisane w postaci krótkich matematycznych funkcji rekurencyjnych - zatem do dzieła 🙂

Rekurencyjne poszukiwanie dzielników

Naszym zadaniem będzie zdefiniowanie funkcji zwracającej 1 jeśli podana liczba n jest liczbą pierwszą oraz 0 w przeciwnym wypadku. Zacznijmy jednak od podania funkcji weryfikującej czy liczba posiada dzielniki.

{\small\text{CzyDzielnik}(n, a, b)=}

{\small=\begin{cases}0&\text{dla}\quad a>b\\1&\text{dla}\quad n \mod a=0\\ \text{CzyDzielnik}(n, a+1, b)&\text{w inn. przyp.}\end{cases}}

Powyższa funkcja zwraca 1 jeśli liczba n posiada dzielnik z przedziału (a,b), oraz 0 w przeciwnym wypadku. Następnie definiujemy wyrażenie reprezentujące naiwny test pierwszości.

{\small\text{CzyPierwsza}(n)=}

{\small=\begin{cases}0&\text{dla}\quad n<2\\ \neg\text{CzyDzielnik}(n,2,\sqrt{n})&\text{dla}\quad n>=2\end{cases}}

Rolą funkcji "CzyPierwsza" jest jedynie "wprawienie algorytmu w ruch" oraz zwrócenie negacji wyniki funkcji "CzyDzielnik". Proste prawda? 🙂 Sprawdźmy więc w mXparser czy to faktycznie działa.

/* Definicja funkcji rekurencyjnych */
Function CzyDzielnik = new Function("CzyDzielnik(n, a, b) = if( a>b, 0, if( n%a = 0, 1, CzyDzielnik(n, a+1, b) ) )");
Function CzyPierwsza = new Function("CzyPierwsza(n) = if( n<2, 0, ~CzyDzielnik(n, 2, sqrt(n)) )", CzyDzielnik);

/* Obliczenie i wyświetlenie wartości */
mXparser.consolePrintln( "CzyPierwsza(1) = " + CzyPierwsza.calculate(1) + ", czas oblicz. = " + CzyPierwsza.getComputingTime() + " s." );
mXparser.consolePrintln( "CzyPierwsza(2) = " + CzyPierwsza.calculate(2) + ", czas oblicz. = " + CzyPierwsza.getComputingTime() + " s." );
mXparser.consolePrintln( "CzyPierwsza(3) = " + CzyPierwsza.calculate(3) + ", czas oblicz. = " + CzyPierwsza.getComputingTime() + " s." );
mXparser.consolePrintln( "CzyPierwsza(4) = " + CzyPierwsza.calculate(4) + ", czas oblicz. = " + CzyPierwsza.getComputingTime() + " s." );
mXparser.consolePrintln( "CzyPierwsza(5) = " + CzyPierwsza.calculate(5) + ", czas oblicz. = " + CzyPierwsza.getComputingTime() + " s." );
mXparser.consolePrintln( "CzyPierwsza(6) = " + CzyPierwsza.calculate(6) + ", czas oblicz. = " + CzyPierwsza.getComputingTime() + " s." );
mXparser.consolePrintln( "CzyPierwsza(7) = " + CzyPierwsza.calculate(7) + ", czas oblicz. = " + CzyPierwsza.getComputingTime() + " s." );
mXparser.consolePrintln( "CzyPierwsza(8) = " + CzyPierwsza.calculate(8) + ", czas oblicz. = " + CzyPierwsza.getComputingTime() + " s." );
mXparser.consolePrintln( "CzyPierwsza(9) = " + CzyPierwsza.calculate(9) + ", czas oblicz. = " + CzyPierwsza.getComputingTime() + " s." );
mXparser.consolePrintln( "CzyPierwsza(10) = " + CzyPierwsza.calculate(10) + ", czas oblicz. = " + CzyPierwsza.getComputingTime() + " s." );

+ wynik

CzyPierwsza(1) = 0.0, czas oblicz. = 0.08 s.
CzyPierwsza(2) = 1.0, czas oblicz. = 0.03 s.
CzyPierwsza(3) = 1.0, czas oblicz. = 0.026 s.
CzyPierwsza(4) = 0.0, czas oblicz. = 0.022 s.
CzyPierwsza(5) = 1.0, czas oblicz. = 0.038 s.
CzyPierwsza(6) = 0.0, czas oblicz. = 0.015 s.
CzyPierwsza(7) = 1.0, czas oblicz. = 0.028 s.
CzyPierwsza(8) = 0.0, czas oblicz. = 0.015 s.
CzyPierwsza(9) = 0.0, czas oblicz. = 0.053 s.
CzyPierwsza(10) = 0.0, czas oblicz. = 0.011 s.

Wygląda na to, że obliczenia są poprawne! Teraz możemy zweryfikować ile jest liczb pierwszych w podanym przedziale, definiując

\pi(n)=\sum_{i=1}^n \text{CzyPierwsza}(i)

/* Definicja wyrażenia sumującego wynik funkcji CzyPierwsza */
Expression pi100 = new Expression("sum(i, 1, 100, CzyPierwsza(i) )");
pi100.addFunctions(CzyPierwsza);

/* Obliczenie i wyświetlenie wyniku */
mXparser.consolePrintln( "Liczba liczb pierwszych w przedziale (1,100) = " + pi100.calculate());

+ wynik

Liczba liczb pierwszych w przedziale (1,100) = 25.0

Pozdrowienia,

Mariusz Gromada

Pamiętajcie, że uruchamiając kody mXparsera należy dodać w nagłówku:

import org.mariuszgromada.math.mxparser.*;

Kod:

 

Zobacz również:

  1. Polowanie na czarownice - czyli zabawy z rekurencją (część 1)
  2. Prędkość ucieczki do nieskończoności - czyli zabawy z rekurencją (część 2)
  3. Rekurencja pośrednia - czyli zabawy z rekurencją (część 4)

Prędkość ucieczki do nieskończoności - czyli zabawy z rekurencją (część 2)

Dziś ciekawostka w nawiązaniu do wpisu z dnia 20 października 2015 roku "Liczba PI ukryta w zbiorze Mandelbrota" - przeanalizujmy równanie rekurencyjne dla liczb rzeczywistych:

x_n=\begin{cases}x_{n-1}^2+\frac{1}{4}+\epsilon,&\text{dla}\quad n>0\\0,&\text{dla}\quad n=0\end{cases}

Zbliżanie się do "ostrza" zbioru Mandelbrota

Powyższe wyrażenie powstaje na bazie równania zbioru Mandelbrota z_n=z_{n-1}^2+c (w liczbach zespolonych), jeśli ograniczymy się do prostej rzeczywistej (dlatego użyłem zapisu x_n) i będzie nas interesowało zbliżanie się elementu x_1=\frac{1}{4}+\epsilon do "ostrza" (ang. "cusp") zbioru o współrzędnych (\frac{1}{4},0).

Mandelbrot - Ostrze

Szybkość ucieczki do nieskończoności

Ustalając odpowiednio małe \epsilon>0 decydujemy jak bardzo chcemy się zbliżyć do "ostrza", teraz zadanie polega na znalezieniu pierwszego n, dla którego x_n>=2. Takie minimalne n jest dobrą miarą prędkości ucieczki x_n do nieskończoności w zależności od wybranego \epsilon. Na marginesie dodam, że zbiór Juli dla równania Mandelbrota jest na powyższym obrazku oznaczony kolorem czarnym, który reprezentuje punkty "nieuciekające" do nieskończoności w trakcie nieskończonej iteracji (ta tematyka jest sama w sobie bardzo ciekawa i zapewne kiedyś coś napiszę o atraktorach).

Rekurencja na rekurencji

Zapiszmy nasze zadanie wykorzystując rekurencję w celu poszukiwania rozwiązania:

f(n,\epsilon)=\begin{cases}f(n+1,\epsilon),&\text{dla}\quad x_n<2\\n,&\text{dla}\quad x_n>=2\end{cases}

Niestety zdefiniowaliśmy rekurencję na rekurencji - to zły znak dla wydajności - cóż - sprawdźmy to używając mXparsera.

/* Definicja funkcji rekurencyjnej */
Function x = new Function("x(n,c) = if( n>0, x(n-1,c)^2+0.25+c, 0 )");
Function f = new Function("f(n,c) = if( x(n,c) >= 2, n, f(n+1, c) )", x);

/* Obliczenia i wyświetlenie wyniku */
mXparser.consolePrintln( "c = 0.01" + ", f(0,c) = " + f.calculate(0, 0.01) + ", czas = " + f.getComputingTime() + " s" );
mXparser.consolePrintln( "c = 0.0001" + ", f(0,c) = " + f.calculate(0, 0.0001) + ", czas = " + f.getComputingTime() + " s" );
mXparser.consolePrintln( "c = 0.000001" + ", f(0,c) = " + f.calculate(0,0.000001) + ", czas = " + f.getComputingTime() + " s" );
mXparser.consolePrintln( "c = 0.00000001" + ", f(0,c) = " + f.calculate(0,0.00000001) + ", czas = " + f.getComputingTime() + " s" );

+ wyczekiwany wynik


c = 0.01, f(0,c) = 30.0, czas = 0.224 s
c = 0.0001, f(0,c) = 312.0, czas = 1.532 s
c = 0.000001, f(0,c) = 3140.0, czas = 37.343 s
c = 0.00000001, f(0,c) = 31414.0, czas = 4068.338 s


Wzorzec prędkości ucieczki

Wow - jaki przedziwny wzorzec liczby iteracji wymaganych, aby przekroczyć 2!! Dostajemy coś, co przypomina \pi, jednak wymaga postawienia przecinka w odpowiednim miejscu! Można również zauważyć, że 100-krotne zmniejszenie \epsilon zwiększa niezbędną liczbę iteracji około 10-krotnie. Zmniejszając \epsilon otrzymamy liczbą coraz bardziej przypominającą \pi 🙂

Czasy wykonania

Tak jak oczekiwaliśmy - czasy wykonania są koszmarne - głównie z dwóch powodów:

  1. wykorzystanie generycznych funkcji rekurencyjnych (możliwe dowolne definicje, jednak kosztem wydajności - konieczność klonowania obiektów, brak możliwości zapamiętywania wartości z poprzednich przebiegów)
  2. definicji rekurencji na rekurencji,

mXparser dostarcza również nieco wydajniejszą wersję rekurencji - tzw. argumenty rekurencyjne z indexem (klasa RecursiveArgument). Wykorzystanie  RecursiveArgument jest ograniczone do jednego indexu, co umożliwia przechowywanie rezultatów z poprzednich iteracji. Zyskujemy na wydajności, tracąc na elastyczności i przejrzystości kodu. Działanie RecursiveArgument zaprezentuję w kolejnej części.

Pozdrowienia,

Mariusz Gromada

BTW: uruchamiając kody mXparsera nie zapomnijcie o dodaniu

import org.mariuszgromada.math.mxparser.*;

Zobacz również:

  1. Polowanie na czarownice - czyli zabawy z rekurencją (część 1)
  2. Naiwny test pierwszości - czyli zabawy z rekurencją (część 3)
  3. Rekurencja pośrednia - czyli zabawy z rekurencją (część 4)
  4. Liczba PI ukryta w zbiorze Mandelbrota

Polowanie na czarownice - czyli zabawy z rekurencją (część 1)

Rozpocznijmy branżowym humorem

Okres średniowiecza, kobieta winna uprawiania magii, kara straszna - spalenie na stosie! Nadszedł dzień, tłum gawiedzi, czarownica na stosie, płomienie, wiedźma krzyczy - więcej drewna! Więcej drewna! Tłum zdziwiony, mimo wszystko spełnia ostatnie życzenie opętanej. Wiedźma nie przerywa - jeszcze więcej drewna! Więcej drewna! Z oddali dobiega nagły i stanowczy sprzeciw - STOP! Czarownica chce przepełnić stos!

🙂

Czym jest rekurencja?

Zazwyczaj o rekurencji myślimy jako o procesie podziału zadania na mniejsze, następnie podziału na jeszcze mniejsze, i jeszcze mniejsze ... dochodząc do zadań, dla których rozwiązanie jest znane. Od tego momentu zaczyna się składanie "mniejszych" rozwiązań w "większe", następnie tych większych w jeszcze większe, ... i w jeszcze większe ... kończąc na rozwiązaniu zadania początkowego. Dla przykładu zapiszmy funkcję n! w postaci rekurencyjnej.

n!=\begin{cases}n\cdot(n-1)!&\text{dla}\quad n>0\\1,&\text{dla}\quad n=0\end{cases}

W celu zobrazowania reprezentacja powyższego podana w mXparser:

/* Definicja funkcji rekurencyjnej */
Function silnia = new Function("s(n) = if( n>0, n*s(n-1), 1 )");

/* Obliczenia i wyświetlenie wyniku */
System.out.println( "n = 0, s(n) = " + silnia.calculate(0) );
System.out.println( "n = 1, s(n) = " + silnia.calculate(1) );
System.out.println( "n = 2, s(n) = " + silnia.calculate(2) );
System.out.println( "n = 3, s(n) = " + silnia.calculate(3) );
System.out.println( "n = 4, s(n) = " + silnia.calculate(4) );
System.out.println( "n = 5, s(n) = " + silnia.calculate(5) );

+ wynik:

n = 0, s(n) = 1.0
n = 1, s(n) = 1.0
n = 2, s(n) = 2.0
n = 3, s(n) = 6.0
n = 4, s(n) = 24.0
n = 5, s(n) = 120.0

Wynik jest zgodny z oczekiwanym. Innym przykład rekurencji to iterowany operator sumowania, niech

A_n=a_1+a_2+\ldots+a_n=\sum_{i=1}^n a_i

Łatwo zauważyć, że

A_n=\begin{cases}a_n+A_{n-1},&\text{dla}\quad n>1\\a_1,&\text{dla}\quad n=1\end{cases}

Jak więc widzicie, rekurencja jest powszechna, często będąc nieco innym sposobem patrzenia na iteracje.

Formalna definicja rekurencji

O rekurencji mówimy jeśli metoda (funkcja, zachowanie, obiekt) może być opisana przez:

  1. elementy bazowe / rozwiązania bazowe;
  2. zestaw reguł, które redukuję (sprowadzają) każdy inny przypadek do (w kierunku) elementów bazowych.

Rekurencja
Powyższe określenie jest szerokie, ale takie być musi, bo i typów rekurencji jest wiele.

Rekurencja jako złożenie funkcji

Jednym (ale nie jedynym) sposobem zapisu ogólnych równań rekurencyjnych jest złożenie funkcji:

f_n=\begin{cases}F\big(f_{n-1},f_{n-2},\ldots,f_{n-k}\big)&\text{dla}\quad n>k\\f_1,f_2,\ldots,f_k&\text{el. baz. dla}\quad n<=k\end{cases}

Dobrą ilustracją powyższego jest ciąg Fibonacciego:

f_n=\begin{cases}0&\text{dla}\quad n=0\\1&\text{dla}\quad n=1\\f_{n-1}+f_{n-2}&\text{dla}\quad n>1\end{cases}

Zapiszmy ciąg Fibonacciego w mXparser:

/* Definicja funkcji rekurencyjnej */
Function fib = new Function("fib(n) = if( n>1, fib(n-1)+fib(n-2), if(n=1,1,0) )");

/* Obliczenia i wyświetlenie wyniku */
System.out.println( "fib(0) = " + fib.calculate(0) );
System.out.println( "fib(1) = " + fib.calculate(1) );
System.out.println( "fib(2) = " + fib.calculate(2) );
System.out.println( "fib(3) = " + fib.calculate(3) );
System.out.println( "fib(4) = " + fib.calculate(4) );
System.out.println( "fib(5) = " + fib.calculate(5) );

+ rezultat:

fib(0) = 0.0
fib(1) = 1.0
fib(2) = 1.0
fib(3) = 2.0
fib(4) = 3.0
fib(5) = 5.0

Rekurencja w roli pętli "For"

Podane wyżej przykłady zapisów rekurencyjnych (n!, suma n-pierwszych wyrazów ciągu, ciąg Fibonacciego) są tak naprawdę rekurencyjną realizacją pętli "for" - znamy przecież dokładnie liczbę niezbędnych operacji do wykonania, a i same operacje są raczej łatwe oraz czytelne - zatem zagnieżdżenie ich w pętli "for" nie powinno spowodować utraty przejrzystości kodu.

Poszukiwanie rozwiązania - czyli rekurencja w roli pętli "While/Until"

W metodach numerycznych często stosuje się strategie rekurencyjne - w takiej sytuacji, będąc w kroku n, weryfikujemy czy propozycja rozwiązania n spełnia kryterium stopu (np. jakość oszacowania), jeśli tak - kończymy z rozwiązaniem n, jeśli nie - przechodzimy do badania propozycji rozwiązania n+1. Procedurę rozpoczynamy od kroku 0 (zerowego).

Przykład: znając definicję silni chcemy znaleźć pierwsze n, dla którego n! >= 100 - takie zadanie formalnie możemy zapisać jako:

S_n=\begin{cases}S_{n+1},&\text{dla}\quad n!<100\\n,&\text{dla}\quad n!>=100\end{cases}

n_{100} = S(0)

Reprezentacja w mXparser:

/* Definicja funkcji rekurencyjnej */
Function S = new Function("S(n) = if( n! < 100, S(n+1), n )"); /* Obliczenia i wyświetlenie wyniku */ System.out.println( "Pierwsze n, że n! >= 100 to n = " + S.calculate(0) );

+ wynik:

Pierwsze n, że n! >= 100 to n = 5.0

Rekurencja bezpośrednia

Wszystkie omawiane wyżej typy rekurencji polegają na wywołaniu z ciała funkcji "siebie samej", dlatego należą do bardziej ogólnej klasy nazywanej rekurencją bezpośrednią.

Wydajność

Implementacje na bazie rekurencji są bardzo czytelne, o często minimalnym rozmiarze kodu. Jednak coś za coś - tracimy bardzo dużo na złożoności obliczeniowej (powtarzane operacje, dzielenie zadań, operacje na stosie) oraz wymogach pamięci (głównie struktura stosu) - to właśnie dlatego czarownica błagała o drewno - licząc na przerwanie procesu z tytułu przepełnienia stosu 🙂

Cdn 🙂

Pozdrowienia,

Mariusz Gromada

BTW: jeśli planujecie uruchamiać kody mXparsera nie zapomnijcie o dodaniu

import org.mariuszgromada.math.mxparser.*;

Zobacz również:

  1. Prędkość ucieczki do nieskończoności - czyli zabawy z rekurencją (część 2)
  2. Naiwny test pierwszości - czyli zabawy z rekurencją (część 3)
  3. Rekurencja pośrednia - czyli zabawy z rekurencją (część 4)