[MODELADO Y SIMULACIÓN] Sistemas caóticos

Diagrama de Bifurcación:

Un diagrama de bifurcación de un sistema dinámico es una estratificación de su espacio de parámetros inducida por la equivalencia topológica, junto con los retratos de fase representativos de cada estrato. Las soluciones estables suelen representarse mediante líneas continuas, mientras que las soluciones inestables se representan con líneas punteadas.

Código en python

import math
max = 480
a = open("datos.txt", "w")

xa = 2.9
xb = 4.0
maxit = 250

for i in range(max):
    r = xa + (xb - xa) * float(i) / (max - 1)
    x = 0.5
    for j in range(maxit):
        x = r * x * (1 - x)
     a.write(str(x)+"\n")
a.close

Código gnuplot

set term postscript eps color 30
set key outside Right
set size 2, 1.5
set output "simugeom.eps"
set pointsize 1
plot "datos.txt" with points lt 2

Billar dinámico:

(Este ejemplo lo pongo por que me pareció divertido)

Un billar dinámico es un sistema dinámico en el cual una partícula alterna entre movimiento rectilíneo y reflexiones especulares en un contorno o frontera.

Los billares capturan toda la complejidad de los sistemas hamiltonianos, desde integrabilidad a movimiento caótico, sin las dificultades de integrar las ecuaciones de movimiento para determinar su mapa de Poincaré.

Código en python

import math
import random
from PIL import Image, ImageDraw
imgx = 480
imgy = 320
image = Image.new("RGB", (imgx, imgy))
draw = ImageDraw.Draw(image)

# Only 1 ball is used!
maxSteps = 20000 # of steps of ball motion (in constant speed)

n = random.randint(1, 7) # of circular obstacles
crMax = int(min(imgx - 1, imgy - 1) / 4) # max circle radius
crMin = 10 # min circle radius

# create circular obstacle(s)
cxList = []
cyList = []
crList = []
for i in range(n):
    while(True): # circle(s) must not overlap
        cr = random.randint(crMin, crMax) # circle radius
        cx = random.randint(cr, imgx - 1 - cr) # circle center x
        cy = random.randint(cr, imgy - 1 - cr) # circle center y
        flag = True
        if i > 0:
            for j in range(i):
                if math.hypot(cx - cxList[j], cy - cyList[j]) < cr + crList[j]:
                    flag = False
                    break
        if flag == True:
            break
    draw.ellipse((cx - cr, cy - cr, cx + cr, cy + cr))
    cxList.append(cx)
    cyList.append(cy)
    crList.append(cr)

# initial location of the ball must be outside of the circle(s)
while(True):
    x = float(random.randint(0, imgx - 1))
    y = float(random.randint(0, imgy - 1))
    flag = False
    for i in range(n):
        if math.hypot(x - cxList[i], y - cyList[i]) <= crList[i]:
            flag = True
            break
    if flag == False:
        break

# initial direction of the ball
a = 2.0 * math.pi * random.random()
s = math.sin(a)
c = math.cos(a)

counter = 0
j=0

for i in range(maxSteps):
    counter = counter +1
    if counter == 15:
        j=j+1
        image.save(str(j)+".png", "PNG")
        counter = 0
    image.putpixel((int(x), int(y)), (255, 255, 255))
    xnew = x + c
    ynew = y + s
    
    # reflection from the walls
    if xnew < 0 or xnew > imgx - 1:
        c = -c
        xnew = x
    if ynew < 0 or ynew > imgy - 1:
        s = -s
        ynew = y
    
    # reflection from the circle(s)
    for i in range(n):
        if math.hypot(xnew - cxList[i], ynew - cyList[i]) <= crList[i]:
            # angle of the circle point
            ca = math.atan2(ynew - cyList[i], xnew - cxList[i])
            # reversed collision angle of the ball
            rca = math.atan2(-s, -c)
            # reflection angle of the ball
            rab = rca + (ca - rca) * 2
            s = math.sin(rab)
            c = math.cos(rab)
            xnew = x
            ynew = y
    
    x = xnew
    y = ynew

Cuando encontre el código originalmente generaba una imagen al finalizar el programa, pero lo modifique para que generara una imagen cada 15 iteraciones y de esa manera observar como se desarrolla.

Para generar imagenes se utiliza la libreria PIL y para generar un video de una secuencia de imagenes se utiliza el comando "ffmpeg -f image2 -i imagen%d.png video.mpeg" donde %d es la numeración de imagenes, ejemplo: imagen1.png, imagen2.png, etc.

Este es el resultado de 1333 png's generados: