Sistemas dinámicos unidimensionales

Simulaciones dinámicas de la transición al Caos

Aprieta aquí para descargar la hoja de trabajo de Maple.

1. Introducción
2. Intervalos invariantes
3. Órbitas en el espacio de fases
4. Análisis gráfico
5. Comportamiento global
6. Comportamiento global de las gráficas
7. Conjuntos de Cantor

Introducción

El estudio de los sistemas dinámicos unidimensionales se ha desarrollado en años recientes debido principalmente a la posibilidad de hacer experimentos con ordenador de iteración de funciones. El germen de los experimentos que proponemos aquí no es por tanto nuevo aunque la manera en que Maple puede manipular las gráficas nos permite ir un poco más lejos en el tratamiento de los resultados. Esta característica visual hace que estos experimentos sean útiles para un curso introductorio a sistemas dinámicos discretos.

El proposito de los diferentes procedimientos que desarrollamos en este trabajo es experimentar la transición al caos en sistemas dinámicos unidimensionales. Mostramos primero como encontrar intervalos invariantes para polinomios, dejando abierta la cuestión de como extender este procedimiento a familias más generales de funciones. Damos a continuación la versión de Maple de un procedimiento para dibujar las órbitas en el espacio de fases, y otro para hacer el análisis gráfico de las órbitas. El análisis gráfico de las órbitas es ampliamente utilizado en estudio de sistemas dinámicos unidimensionales particulares. Los procedimientos que damos aquí (excepto el mencionado en primer lugar que solo es válido para polinomios) son aplicables a cualquier función real y son por tanto útiles para experimentar el comportamiento caótico de otros sistemas dinámicos unidimensionales, aparte de la conocida familia cuadrática.

A continuación, usando el procedimiento para el análisis gráfico de órbitas - un desarrollo similar se podría haber hecho con el procedimiento para calcular las órbitas en el espacio de fases - y aprovechando otra vez las facilidades de "display" de Maple, mostramos procedimientos para experimentar el comportamiento global de la familia, primero variando el punto pero dejando fijo el parámetro, luego fijando el punto y variando el parámetro. Después, mostramos como experimentar el cambio simultaneo de ambos valores. Finalmente vemos como experimentar con la familia de aplicaciones, incluso en el caso de no existir intervalos invariantes. En muchos casos, lo que se obtiene es un conjunto de Cantor invariante. Para hacer el trabajo homogeneo y no demasiado largo, no hemos incluido procedimientos para calcular valores numéricos de las órbitas.

Puesto que todos los procedimientos se basan en "displays insequence", permiten un tratamiento dinámico e interactivo de los resultados. Todos ellos tienen como variables la función generando la familia unidimensional y son por tanto aplicables a una amplia clase de funciones. Hemos usado como función de prueba la conocida familia cuadrática

> f:=x^2 +param:

y la función logística

> g:=param*x*(1-x):

Proponemos también aplicar los mismos procedimientos al estudio de familias más complicadas

> h:=x*(x-1)*(x-2)*(x-3)+param:

Intervalos invariantes

La primera información importante requerida para estudiar el sistema dinámico es la localización de intervalos invariantes para la familia. El siguiente procedimiento nos permite hacer una primera acotación de un intervalo donde estudiar la función, aunque no ha de ser necesariamente un intervalo invariante. Puesto que es dependiente de la resolución de una ecuación, sólo funciona bien para polinomios. Se basa en localizar los puntos fijos de la función y las preimagenes de éstos. Fuera del intervalo resultante, en el caso de polinomios y siempre que los límites en más infinito y menos infinito no sean respectivamente menos infinito y más infinito, el comportamiento de las iteraciones de la función está completamente determinado. La entrada del procedimiento es:

• f (familia de funciones dependiendo de un parámetro "param"),
• c (valor escogido del parámetro "param").

> interv:=proc(f,c)
> local g,fp,pfp,i,lll;
> g:=subs(param=c,f);
> fp:=[fsolve(g=x,x)];
> if nops(fp)=0 then lll:=[];
> else pfp:={};
> for i from 1 to nops(fp) do;
> pfp:=pfp union {fsolve(g=subs(x=fp[i],g))};
> od;
> pfp:=sort([op(pfp)]);
> lll:=[pfp[1],pfp[nops(pfp)]];
> fi;
> lll;
> end:

El siguiente procedimiento localiza intervalos invariantes maximales de f para un valor dado c del parámetro "param" pero, puesto que sigue siendo dependiente de la resolución de una ecuación, solo funciona bien para polinomios. Además hacemos uso en el procedimiento de la función "sturm" que calcula el número de raices de un polinomio en un intervalo. La entrada del procedimiento es la misma que la del procedimiento anterior. El resultado será una lista de tal forma que cada elemento de lugar impar y el elemento de lugar par siguiente son los extremos de un intervalo invariante.

> invinterv:=proc(f,c)
> local g,cp,fp,ll,pfp,i,j,pcp,s,r,l,ma,k,mi,lll;

Localizamos primero los puntos críticos de f para el valor actual del parámetro "param".

> g:=subs(param=c,f);
> cp:=[solve(diff(g,x)=0,x)];

localizamos después los puntos fijos de f,

> fp:=[fsolve(g=x,x)];

Si no existen puntos fijos es fácil ver que toda órbita diverge a más o menos infinito, según que el grafo de f esté por encima o por debajo de la diagonal,

> if nops(fp)=0 then ll:=[];

Si existe uno o más puntos fijos consideramos como candidatos a intervalos invariantes aquellos que tengan como extremos puntos fijos (al menos uno) o preimágenes de puntos fijos (obsérvese que todo intervalo invariante debe contener un punto fijo).

> else ll:=[]; pfp:={};
> for i from 1 to nops(fp) do;
> pfp:=pfp union {fsolve(g=subs(x=fp[i],g))};
> od;
> pfp:=sort([op(pfp)]);

Para cada uno de los puntos arriba mencionados, buscamos todos los posibles intervalos invariantes que tengan este punto como extremo inferior,

> for i from 1 to nops(pfp)-1 do;
> l:=[];
> for j from i+1 to nops(pfp) do;
> pcp:=[pfp[i],pfp[j]];
> readlib(sturm):
> s:=sturmseq(diff(g,x),x);
> r:=sturm(s,x,pfp[i],pfp[j]);
> if 0<r then
> pcp:=[op(pcp),fsolve(diff(g,x)=0,x,pfp[i]..pfp[j])];
> fi;
> ma:=max(seq(subs(x=pcp[k],g),k=1..nops(pcp)));
> mi:=min(seq(subs(x=pcp[k],g),k=1..nops(pcp)));
> if pfp[i]<=mi and ma<=pfp[j]
> then l:=[pfp[i],pfp[j]];
> fi;
> od;
> ll:=[op(ll),op(l)];
> od;
> fi;

En este punto hemos obtenido una lista cuyos pares de elementos son los intervalos invariantes maximales para cada posible extremo inferior en la lista de preimágenes de puntos fijos. En lo que sigue suprimimos solapamientos para obtener una familia de intervalos invariantes disjuntos,

> if ll=[] then lll:=[];
> else lll:=[ll[1],ll[2]];
> for i from 1 to nops(ll)/2-1 do;
> if ll[2*i]<>ll[2*(i+1)] then
> lll:=[op(lll),ll[2*(i)+1],ll[2*(i+1)]];
> fi;
> od;
> fi;
> lll;
> end:

Órbitas en el espacio de fases

El siguiente procedimiento representa de una forma dinámica los n primeros elementos de la órbita de un punto p. Unimos puntos consecutivos con un semiarco, superior si el movimiento es hacia la derecha e inferior si es hacia la izquierda. La entrada del procedimiento es:

• f (familia de funciones dependiendo de un parámetro "param"),
• c (valor escogido del parámetro "param"),
• p (punto cuya órbita se quiere calcular),
• n (número de elementos de la órbita),
• [d,e] (intervalo donde f está restringida).

> orbphase:=proc(f,c,p,n,d,e)
> local t,a,l,q,b,k;
> a:=evalf(p);
> l:=[plot([[a,0]],x=d..e,style=point,symbol=BOX)]:
> with(plots);
> q:=[plots[display](l)]:
> for k from 1 to n do;
> b:=evalf(subs(param=c,x=a,f));
> l:=[op(l),plot(signum(b-a)*sqrt(((a-b)/2)^2-(x-((a+b)/2))^2),x=min(a,b)..max(a,b))]:
> q:=[op(q),plots[display](l)]:
> a:=b;
> od;
> plots[display](q,insequence=true,scaling=constrained);
> end:

> orbphase(f,-1.5,0.1,64,-2,2);

> orbphase(g,3.5,0.1,64,0,1);

Análisis Gráfico

Damos ahora un procedimiento para hacer el análisis gráfico de las órbitas. La entrada del procedimiento es la misma que la del anterior.

> orbn:=proc(f,c,p,n,d,e)
> local p0,g,u,v,i,p1,l,a,b,k,p2;
> g:=subs(param=c,f);
> u:=evalf(d);
> v:=evalf(e);
> p0:=plot([[p,0]],x=u..v,u..v,style=point,symbol=BOX);
> p1:=plot({g,x}, x=u..v, y=u..v,color=blue):
> a:=evalf(p);
> l:=[[a,0]];
> for k from 1 to n while u<a and a<v do;
> b:=evalf(subs(x=a,g));
> l:=[op(l),[a,b],[b,b]];
> a:=b;
> od;
> p2:=plot(l,u..v,u..v,style=line,color=red):
> with(plots);
> plots[display]({p0,p1,p2},scaling=constrained);
> end:

> orbitera(f,-1.5,0.1,96,-2,2);

> orbitera(g,3.5,0.1,96,0,1);

Comportamiento global

Variando los diferentes parámetros de la función "orbn" obtenemos varios procedimientos que nos permiten estudiar el comportamiento global de la familia de funciones.

> orbitall(f,-1.5,64,-2,2,40);

> orbitall(g,3.5,64,0,1,40);

> orbifurc(f,0.25,-2,40,0,64,-2,2);

> orbifurc(g,0,4,40,0.5,64,0,1);

Comportamiento global de las gráficas

Finalmente podemos utilizar el siguiente procedimiento para experimentar el comportamiento asintótico global de toda la familia. Utilizaremos la siguiente tabla de colores:

> col := table([1=black, 2=blue, 3=navy, 4=coral, 5=cyan, 6=brown, 7=gold, 8=green, 9=gray, 10=grey, 11=khaki, 12=magenta, 13=maroon, 14=orange, 15=pink, 16=plum, 17=red, 18=sienna, 19=tan, 20=turquoise, 21=violet, 22=wheat, 23=white, 24=yellow, 0=aquamarine]):

En el primer procedimiento se representan "insequence" las iteraciones de la grafica de f. La entrada es:

• f (familia de funciones dependiendo de un parámetro "param"),
• c ( valor escogido del parámetro "param"),
• n (número de iteraciones de la gráfica),
• [d,e] intervalo donde f está restringida.

> grafanimate:=proc(f,c,n,d,e)
> local t,g,q,l,h,qq,ll,i;
> g:=subs(param=c,f):
> q:=[plot(x,x=d..e,y=d..e,color=col[1]),plot(g,x=d..e,y=d..e,color=col[2])];
> l:=[plots[display](q)];
> h:=g;
> for i from 1 to n-1 do:
> h:=subs(x=h,g):
> qq:=[plot(x,x=d..e,y=d..e,color=col[1]),plot(h,x=d..e,y=d..e,color=col[2],numpoints=400)]:
> l:=[op(l),plots[display](qq)]:
> od:
> plots[display](l,insequence=true,view=[d..e,d..e],scaling=constrained);
> end:

Podemos representar también las iteraciones simultaneamente. El siguiente procedimiento hace eso para las iteraciones entre la m y la n. El procedimiento, cuya entrada es:

• f (familia de funciones dependiendo de un parámetro "param"),
• c ( valor escogido del parámetro "param"),
• [m,n] (intervalo de iteraciones considerado),
• [d,e] (intervalo donde f está restringida),

es el siguiente.

> graf:=proc(f,c,m,n,d,e)
> local t,g,h,l,i,j;
> g:=subs(param=c,f):
> h:=g;
> for i from 1 to m-1 do: h:=subs(x=h,g): od:
> l:=[plot(x,x=d..e,y=d..e,color=col[1])]:
> l:=[op(l),plot(h,x=d..e,y=d..e,color=col[2])]:
> for j from m to n-1 do:
> h:=subs(x=h,g):
> l:=[op(l),plot(h,x=d..e,y=d..e,color=col[modp(j-m+3,24)],numpoints=200)]:
> od:
> plots[display](l,scaling=constrained);
> end:

> orbigraf(f,0.25,-2,17,1,10,-2,2);

> orbigraf(g,0,4,16,1,12,0,1);

Conjuntos de Cantor

Finalmente estudiamos el comportamiento de la familia cuando no hay intervalos invariantes. Para la familia cuadrática, esto ocurre para P mayor o igual que 4. Dado un intervalo apropiado ([0,1] para la familia cuadrática), existen puntos del intervalo cuya imagen bajo una iteración se sale del mismo. Otros puntos necesitan 2 iteraciones para salirse. Este comportamiento se puede visualizar con el siguiente procedimiento dinámico en el que se representan simultaneamente "todas" las órbitas y se observa como se escapan del intervalo. Las órbitas que permanecen completamente en el intervalo se puede ver que forman un conjunto de Cantor. La entrada del procedimiento es

• f (familia de funciones dependiendo de un parámetro "param"),
• c (valor escogido del parámetro "param"),
• r (número de puntos del dominio),
• n (número de iteraciones),
• [d,e] (intervalo donde f está restringida).

> orbcantor:=proc(f,c,r,n,d,e)
> local t,g,q,j,p1,p2,i,a,l,ll,lll,b,bb,k;
> g:=subs(param=c,f);
> q:=[];
> l:=[seq(d+(e-d)*i/r,i=0..r)];
> for j from 0 to n do;
> p1:=[plot({g,x},x=d..e,d..e,scaling=constrained,color=blue)]:
> p2:=[plot({g,x},x=d..e,d..e,scaling=constrained,color=blue)]:
> ll:=[];
> for i to nops(l) do;
> a:=evalf(l[i]);
> b:=evalf(subs(x=a,g));
> lll:=[[a,0],[a,b]];
> for k from 1 to j while d<b and b<e do;
> bb:=evalf(subs(x=b,g));
> lll:=[op(lll),[b,b],[b,bb]];
> b:=bb;
> od;
> p1:=[op(p1),plot(lll,x=d..e,style=line)]:
> if d<b and b<e then
> p2:=[op(p2),plot(lll,x=d..e,style=line)]:
> ll:=[op(ll),l[i]]:
> fi:
> od:
> q:=[op(q),plots[display](p1),plots[display](p2)]:
> l:=ll;
> od:
> plots[display](q,insequence=true,scaling=constrained);
> end:

> orbcantor(f,-2.2,100,9,-2.065247584,2.065247584);

> orbcantor(g,4.4,100,9,0,1);