Procedura iteracje rysuje uciąglony wykres funkcji, która liczbie naturalnej j przyporządkowuje j-tą itercję funkcji f na punkcie startu x0. Zakres zmienności j - 1..n

> iteracje:=proc(f,x0,n)
local i:
x(0):=x0:
for i from 1 to n do
x(i):=f(x(i-1))
od:
PLOT(CURVES([seq([j,x(j)],j=0..n)]));
end:

> iteracje(x->2*x*(1-x),.1,50);
iteracje(x->3*x*(1-x),.1,50);

[Maple Plot]

[Maple Plot]

> iteracje(x->3.5*x*(1-x),.1,50);
iteracje(x->4*x*(1-x),.1,50);

[Maple Plot]

[Maple Plot]

Teraz to samo, ale dla dwóch punktów startu x0 i y0

> iteracje2:=proc(f,x0,y0,n)
local i:
x(0):=x0:y(0):=y0:
for i from 1 to n do
x(i):=f(x(i-1)):
y(i):=f(y(i-1))
od:
PLOT(CURVES([seq([j,x(j)],j=0..n)]),POINTS(seq([j,y(j)],j=0..n)));
end:

> iteracje2(x->4*x*(1-x),.1,.10001,30);

[Maple Plot]

Procedura gęstość zaznacza na ukadzie wspólrzędnych (\mu,x) trajektorie ukadu dyskretnego zadanego przez funkcję f zależną od parametru \mu. Dla kolejnych wartości parametru \mu od r0 do rk (m - liczba pośrednich wartości tego parametru) obliczane są trajektorie o punktach startu k*siec, przy czym k zmienia się od k=1 do k=dlugosc/siec (dlugosc to dlugość przedzialu, na którym dziala uklad dynamiczny). Trajektorie mają dlugość n, a 10 początkowych punktów trajektorii jest opuszczanych. W rezultacie, jesli trajektoria zbliża się do jakiegoś punktu (czy trajektorii okresowej), to na wykresie zobaczymy tylko ten punkt (cz też trajektorię okresową) za to ,,pogrubiony". Jeśli dla jakiejś wartości parametru trajektorie równomiernie rozkadają się na calym odcinku, to cay przedzial nad tą wartościa \mu jest zaczerniony.

> gestosc:=proc(f,dlugosc,r0,rk,m,siec,n)
local i,j,k,step:
r(0):=r0:step:=(rk-r0)/m:
for j from 1 to m do
r(j):=r(0)+j*step:
for k from 1 to dlugosc/siec do
x(j,0,k):=k*siec:
for i from 1 to n do
x(j,i,k):=f(r(j),x(j,i-1,k))
od:
od:
od:
PLOT(POINTS(seq(seq(seq([r(j),x(j,i,k)],i=20..n),k=2..dlugosc/siec-1),j=1..m),SYMBOL(POINT))):
end:

> gestosc((mu,x)->mu*x*(1-x),1.,3.,4.,100,.1,100);

[Maple Plot]

Procedura rozklad rysuje rozklad funkcji gęstości prawdopodobieństwa rozkladu trajektorii dla dyskretnego ukladu dynamicznego danego przez funkcję g:[0,1]->[0,1] przy punkcie startu x0, dlugości trajektorii - n i liczbie m kawalków, na które dzielimy przedzial [0,1].

> rozklad:=proc(g,x0,n,m)
local x,i,j,licz:
x(0):=x0:
for j from 1 to m do
licz(j):=0
od:
for i from 1 to n do
x(i):=g(x(i-1)):
for j from 1 to m do
if x(i)>=(j-1)/m and x(i)<j/m then licz(j):=licz(j)+1 fi:
od:
od:
PLOT(CURVES([seq([evalf(j/m),licz(j)/n], j=1..m)]));
end:

> rozklad(x->4*x*(1-x),.2,10000,40);

[Maple Plot]

Procedura lapunov rysuje wykres wykladniak Lapunowa dla ukladu dyskretnego danego przez funkcję f w zależności od parametru r zmienającego się od r0 do rk (liczba punktów pośrednich - m). Punkt startu - x0. Liczba n - parametr procedury pochodzi przybliżenia wartości granicy w definicji wykladnika Lapunowa przez n-ty wyraz ciągu. Dodatkowy parametr epsilon jest dodany, aby wyeliminować przypadek, gdy trajektoria za bardzo zbliży się do punktu, gdzie pochodna f znika (bliżej niż 10^(-3)). Skladnik związany z takim punktem trajektorii, albo nie móglby być obliczony (logarytm) albo bylby zbyt bliski -nieskończoności, co zniszczyloby dokladność poprzednich obliczeń.

la:=lim sup {n->infinity)1/n sum(i=1,i=n) ln |f ' (x(i)|

> lapunov:=proc(f,r0,rk,m,x0,n,epsilon)
local i,j,step,r,la:
step:=(rk-r0)/m:r(0):=r0:x(0):=x0:
for j from 1 to m do
r(j):=r(0)+j*step:
for i from 1 to n do
x(i):=f(r(j),x(i-1)):
if abs(evalf(subs(x=x(i),diff(f(r(j),x),x))))<10^(-3) then x(i):=x(i)+epsilon fi:
od:
la(j):=1/n*sum(ln(abs(subs(x=x(k),diff(f(r(j),x),x)))),k=1..n):
od:
PLOT(CURVES([seq([r0+k*step,la(k)],k=1..m)])):
end:

> lapunov((r,x)->r*x*(1-x),2.5,4,100,.2,100,.001);

[Maple Plot]

Jeśli utworzymy średnie trajektorii

x^ := lim(n->infinity) 1/n sum (i=0,i=n-1) f^i(x0)= lim(n->infinity) 1/n sum (i=0,i=n-1) f(x(i)),

to możemy zdefiniować funkcję autokorelacji wzorem

C(m):= lim sup 1/n sum(i=0,i=n-1) ((x(i+m)-x^)(x(i)-x^), m=1,2,3,...

Dla trajektorii okresowej o okresie p C(m)=0 dla m=/ k*p (k=1,2,3,...), natomiast C(k*p) jest liczbą dodatnią Dla trajektorii bliskich od pewnego miejsca okresowym jest podobnie

(jeden ,,pik" funkcji C). Dla ukladów chaotycznych funkcja autokorelacji ma wiele ,,pików".

Procedura autokorel rysuje wykres funkcji m->C(m). Jej argumentami jest funkcja g definiująca uklad dynamiczny, punkt startu z0 i liczba iteracj n.

> autokorel:=proc(g,z0,n)
local i,m,xsr,C:
x(0):=z0:
for i from 1 to 2*n do
x(i):=g(x(i-1))
od:
xsr:=sum(x(j),j=0..2*n)/(2*n+1):
for m from 1 to n do
C(m):=sum((x(j+m)-xsr)*(x(j)-xsr),j=0..n)/(n+1)
od:
PLOT(CURVES([seq([j,C(j)],j=1..n)]));
end:

> autokorel(x->3.81*x*(1-x),.2,100);

[Maple Plot]

>