Od układu Lorenza do równania logistycznego

Edward Lorenz analizując zachowanie tego układu w 1963 roku, dokonał jeszcze jednego ważnego kroku - powiązał on dynamikę ciągłego trójwymiarowego układu dynamicznego z zachowaniem tzw. dyskretnego układu dynamicznego. Przeanalizujmy po kolei kroki, które wykonał Lorenz. Mając trajektorię układu, dla zmiennej \(z(t)\) obliczył on wartości \(z_i\) wszystkich lokalnych maksimów. Następnie na wykresie naniósł ich kolejne wartości tzn. narysował pary \((z_i,z_{i+1})\). Okazało się, że dla parametrów w których układ jest chaotyczny pary te układają się na pewnej krzywej \(F\) takiej, że \(z_{i+1}=F(z_i)\). Można by teraz zapomnieć skąd wzięły się wartości \(z_i\) bo dysponując jedynie krzywą \(F\), z jednego stanu możemy otrzymać kolejny. Układ po takiej operacji jest jednowymiarowy, gdyż stan określony jest przez jedną liczbę \(z_i\), ale ewolucja w czasie jest dokonywana w sposób skokowy, za pomocą przekształcenia \(z_{i+1}=F(z_i)\). Taki układ dynamiczny ze skokową ewolucją w czasie nazywa się właśnie dyskretnym układem dynamicznym. Można się też spotkać z określeniem „system funkcji iterowanych” ( z ang. iterated function system, IFS). Układy te stanowią są znaną już dziś z zaskakująco skomplikowanego zachowania, pomimo swojej prostoty. W dalszej części przejdziemy do analizy fascynujących własności tych układów. Zanim jednak to zrobimy, spróbujmy samodzielnie odtworzyć wynik Edwarda Lorenza.

Mając trajektorię układu Lorenza musimy się zastanowić jak z niej wyłowić lokalne maksima? Oczywiście ponieważ rozwiązanie układu Lorenza jest ciągłą funkcją czasu, powinniśmy zastosowywać metody badania przebiegu zmienności funkcji, czyli policzyć pierwszą pochodną, znaleźć jej wszystkie zera na zadanym odcinku i sprawdzić czy tak uzyskane ekstrema są maksimami. Niestety rozwiązanie układu Lorenza nie jest dane wzorem analitycznym. I tu jest pies pogrzebany, bo metodologia postępowania znana ze szkoły średniej wymaga algebraicznego obliczenia pochodnej. Dlatego zrobimy inaczej. Procedura desolve_odeint daje nam tabelę z wynikami. Zakładając ze odstępy pomiędzy kolejnymi punktami czasu w tej tabeli są odpowiednio małe, możemy policzyć lokalne maksima dla ciągu, zauważając, że punkt \(z_i\) jest lokalnym maksimum jeżeli jego otoczenie jest od niego mniejsze czyli zachodzi \(z_{i-1}<z_{i}\) i \(z_{i+1}<z_i\). Oczywiście nie będą to „prawdziwe” maksima funkcji \(z(t)\) a jedynie ich przybliżenie. Jedną z możliwości jest napisanie pętli (zachęcamy do zrobienia tego własnoręcznie), która dla każdego punktu z tabeli sprawdziła by czy zachodzą powyższe warunki i jeśli tak, to zapisałaby punkt na listę maksimów. Mając jednak do dyspozycji „oręż” w postaci biblioteki numpy możemy zrobić to w praktycznie jednej linii kodu. Oznaczając przez Z tablicę z wartościami trzeciej zmiennej układu Lorenza obliczamy najpierw tablicę różnic kolejnych elementów:

Zp = np.diff(Z)

następnie znajdujemy miejsca (np.nonzero) w których kolejne różnice mają przeciwny znak:

idx = np.nonzero(Zp[1:]*Zp[:-1]<0)[0]

i ostatecznie wyciągamy z tablicy Z te elementy:

Zm = Z[idx+1]

(pytanie do czytelnika - skąd to +1?)

Wypróbujmy czy taka procedura zadziała np. na funkcji \(\sin(x)\):

Zaskakujące jest to, że wszystkie punkty znajdują się na jednej linii. Oznacza to bowiem, że otrzymany obrazek jest wykresem funkcji przeprowadzającej jedno maximum w kolejne. Możemy więc zapomnieć o układzie równań różniczkowych a jedynie badać jak jedno maksimum przechodzi w drugie. Skoro wyjściowy układ Lorenza ma własności chaotyczne: nieregularność i czułość na warunki początkowe, to powinniśmy zaobserwować to również w zachowaniu się ciągu \(z_i\) generowanym przez iteracje funkcji, której wykres wcześniej narysowaliśmy:

(1)\[z_{i+1} = F(z_{i})\]

Można również odwrócić problem i zastanowić się dla jakich funkcji \(F\), układ dynamiczny ze skokowym, lub jak częściej się mówi, dyskretnym czasem, będzie posiadał własności chaotyczne. Przekonamy się, że takie rozważania doprowadzą nas do kolejnej rodziny zagadnień matematycznych. Zanim jednak do tego przejdziemy, musimy się zastanowić nad populacją karaluchów.