Comportamiento de sueño atractor fractal

Question

Esta pregunta es similar a este: Propiedades de Sueño del Rey fractal

El fractal se describe a continuación: Softology - Visiones Del Caos 2D Atractor Extraño Tutorial

Yo accidentalmente modefied la fórmula a la siguiente:

$$x_{n+1}=\sin(by_{n})-c\sin(bx_{n})$$ $$y_{n+1}=\sin(ax_{n})+d\sin(ay_{n})$$

So just the sign of $c$ is flipped, so the fractals are really the same.

I'll give a short description of how it is generated.

Given a random starting position of $x_0$ and $y_0$, $x_0 = -2$ and $y_0 = -2$ work quite well, the number of times every pixel is hit. This number, if plotted on a logarithmic scale, gives rise to the strange attractor.

Here are some pictures rendered at 8000000 iterations each, $x$ and $y$ from $-3$ to $3$ each:

$$a = 1.5,\ b= -0.7,\ c=-1.4,\ d=1.2$$

$$a = -2.5,\ b= -1.5,\ c=1.9,\ d=1.9$$

$$a = 2,\ b= -1.5,\ c=1,\ d=1$$

$$a = 2,\ b= 0.3,\ c=3.4,\ d=-0.9$$

$$a = -1.5,\ b= 0.7,\ c=1.7,\ d=-2$$

Ahora, como usted puede ver, a veces, el fractal se extiendan, a veces es apretados, pero a veces el atractor se derrumba realmente pequeños bucles, o incluso sólo un par de puntos. Si cambia los valores de un atractor continuamente a los valores de otro atractor, los cambios son altas de que entre por un breve momento el atractor se derrumba y, a continuación, cambia.

Para comprobar en qué puntos que iba a suceder de un crudo método fue aplicado. El número de píxeles individuales que fueron éxito fue contado. Para acelerar las cosas en la parrilla, a una escala menor.

Cuando este valor se representa en un mapa donde $a$ $b$ son variadas y $c$ $d$ son fijos los siguientes mapas se obtienen, donde $a$ es variada de $-4$ $4$$x$- eje y $b$ es variada de $-4$ $4$$y$-eje. Ignorar los valores de $a$ amd $b$, que no son exactos debido a la acumulación de errores de redondeo. Blanco significa alto número de píxeles de golpe, negro por el bajo número de píxeles de golpe:

¿Se sabe algo de este mapa? Se parece un poco a un fractal. Debido a esta forma parece que no hay ninguna fórmula simple para determinar el estado de la atractor beforhand. También parece ser simétrica.

Obviamente, este método es muy crudo en la determinación del estado de la fractal. También es lento, la última foto que tomó más de 4 horas de render. Hay un método mejor?

Para cualquier persona tratando de reproducir las imágenes del mapa, aquí está el código de procesamiento para que:

float a, b, c, d, x, y, xp, yp;
int hit[][] = new int[200][200];
PFont f;

void setup() {
  size(800, 800);
  f = createFont("Arial", 16, true);  
  a = -4.2;
  b = -0.5;
  c = -1;
  d = -2.2;  
  for (int i = 0; i < 200; i++) {
    for (int o = 0; o < 200; o++) {
      hit[i][o] = 0;
    }
  }
}

void draw() {
  background(0); 
  a = -4;
  while (a < 4) {
    b = -4;
    while (b < 4) { 
      for (int i = 0; i < 200; i++) {
        for (int o = 0; o < 200; o++) {
          hit[i][o] = 0;
        }
      }
      xp = -2;
      yp = -2;
      for (int i = 0; i < 100000; i++) {
        x = sin(yp * b) - c * sin(xp * b);
        y = sin(xp * a) + d * sin(yp * a);
        hit[int(constrain(map(x, -5, 5, 0, 100), 0, 100))][int(constrain(map(y, -5, 5, 0, 100), 0, 100))] = 1;
        xp = x;
        yp = y;
      } 
      int sum = 0;
      for (int i = 0; i < 200; i++) {
        for (int o = 0; o < 200; o++) {
          sum += hit[i][o];
        }
      }
      fill(constrain(10*map(float(sum)/40000, 0, 1, 0, 255), 0, 255));
      noStroke();
      rect(int(map(a, -4, 4, 0, 800)), int(map(b, -4, 4, 0, 800)), 100, 100); 
      b += 0.01;
    }
    a += 0.01;
  } 
  fill(255, 0, 0);
  noStroke();
  textFont(f, 20);
  text("a = " + str(a) + "\nb = " + str(b) + "\nc = " + str(c) + "\nd = " + str(d), 10, 30); 
}

Y para cualquier persona que está tratando de probar varios valores de $a$, $b$, $c$ y $d$ en el atractor sí, aquí hay algunos básicos de procesamiento de código para un negro/blanco de la imagen:

float a, b, c, d, x, y, xp, yp;
PFont f;

void setup() {
  size(800, 800);
  f = createFont("Arial", 16, true); 
  a = -4.2;
  b = -0.5;
  c = -1;
  d = -2.2;
  background(0);
}

void draw() {
  background(0);  
  xp = random(-2, 2);
  yp = random(-2, 2);  
  for (int i = 0; i < 10000; i++) {
    x = sin(yp * b) - c * sin(xp * b);
    y = sin(xp * a) + d * sin(yp * a);  
    set(int(map(x, -5, 5, 0, 800)), int(map(y, -5, 5, 0, 800)), color(255));  
    xp = x;
    yp = y;
  }  
  fill(255, 0, 0);
  noStroke();
  textFont(f, 20);
  text("a = " + str(a) + "\nb = " + str(b) + "\nc = " + str(c) + "\nd = " + str(d), 10, 30);
}

Muchísimas gracias!

---

Estas animaciones muestran una variable continua siendo muy variadas:

https://youtu.be/3iYnrgRpMX8 https://youtu.be/ZveuYBnm6Us https://youtu.be/-PgzRCYboeU

Aquí es un crudo de procesamiento de puerto para la prestación del atractor: https://pastebin.com/nVCJipkA

---

Se hicieron algunos esfuerzos para dibujar un diagrama de bifurcación, que es obviamente difícil, porque es un 2D-atractor.

Un crudo solución es simplemente bajando el $x$-coordinar y, en su lugar de trazado $a$ $x$- eje.

$$a \in [0.59; 0.75],\ b = 3,\ c = -1.4,\ d = 1.2$$

Zoom en el primer superior de la pequeña bifurcación del árbol:

Aquí es un script de procesamiento para la elaboración de bifurcaciones: https://pastebin.com/9ndL2dH3

Y esta es una animación de la bifurcación de la trama: https://youtu.be/tRgM-HhtXVk

Y algunos más: https://youtu.be/6WBgtXTPgbA https://youtu.be/ptnNFpxixow

Tras la visualización de un logarítmica de la parcela es visible que no todos los puntos en las regiones densas son golpeados con la misma frecuencia:

Zoom:

Otra parcela:

Código para logarítmica colorscale: https://pastebin.com/Xgjkg8tM

Answer 1

Aquí está una de Lyapunov fractal trazan para $d = 1, c = -4$$a,b \in [-4,4]$. El amarillo es estable (el más grande exponente de Lyapunov $\lambda_1 < 0$) y el azul es caótico ($\lambda_1 > 0$). Hay algún código que calcula una estimación de la dimensión fractal de la Lyaponov exponentes $\lambda_k$, pero las condiciones previas ($\lambda_1>0 \wedge \lambda_1+\lambda_2<0$) no siempre se cumple lo que significa que se le da a 0 con demasiada frecuencia como para ser útil.

Código fuente (que es este código BÁSICO para el mapa de Henon convierte en C, con cambios menores):

// usage: gcc -O3 l.c -lm -fopenmp -march=native && ./a.out > l.ppm
// converted from henon bas at http://sprott.physics.wisc.edu/chaos/lespec.htm
#include <math.h>
#include <stdio.h>
typedef float R; // if you change this, change sinf() to sin() or sinl() etc
void dream(R x[], R xnew[], R a, R b, R c, R d)
{
  xnew[0]=sinf(b * x[1]) + c * sinf(b * x[0]);
  xnew[1]=sinf(a * x[0]) + d * sinf(a * x[1]);
  R J[2][2];
  J[0][0]=b * c * cosf(b * x[0]);
  J[0][1]=b     * cosf(b * x[1]);
  J[1][0]=a     * cosf(a * x[0]);
  J[1][1]=a * d * cosf(a * x[1]);
  R D[2][2] = { { x[2], x[3] }, { x[4], x[5] } };
  xnew[2] = J[0][0] * D[0][0] + J[0][1] * D[1][0];
  xnew[3] = J[1][0] * D[0][0] + J[1][1] * D[1][0];
  xnew[4] = J[0][0] * D[0][1] + J[0][1] * D[1][1];
  xnew[5] = J[1][0] * D[0][1] + J[1][1] * D[1][1];
}
#define n 2
#define nn (n * (n + 1))
#define H 1024
#define W 1024
unsigned char out[H][W][3];
int main(int argc, char **argv)
{
  R d = 1;
  R c = 4;
  fprintf(stdout, "P6\n%d %d\n255\n", W, H);
  for (int py = 0; py < H; ++py)
  {
    int progress = 0;
    R b = -4 + (py + 0.5) / H * 8;
    #pragma omp parallel for
    for (int px = 0; px < W; ++px)
    {
      R a = -4 + (px + 0.5) / W * 8;
      R dky = 0;
      R x[nn], xnew[nn], v[nn], ltot[n], znorm[n], gsc[n];
      int irate = 10;
      int io = 10000;
      for (int i = 0; i < n; ++i) v[i] = -200;
      int m = 0;
      int t = 0;
      for (int i = n; i < nn; ++i) v[i] = 0;
      for (int i = 0; i < n; ++i) { v[(i + 1) * (n + 1) - 1] = 1; ltot[i] = 0; }
      while (1)
      {
        for (int j = 0; j < irate; ++j)
        {
          for (int i = 0; i < nn; ++i) x[i] = v[i];
          dream(x, xnew, a, b, c, d);
          for (int i = 0; i < nn; ++i) v[i] = xnew[i];
          t += 1;
        }
      znorm[0] = 0;
      for (int j = 0; j < n; ++j)
        znorm[0] += v[n * (j + 1)] * v[n * (j + 1)];
      znorm[0] = sqrtf(znorm[0]);
      for (int j = 0; j < n; ++j)
        v[n * (j + 1)] /= znorm[0];
      for (int j = 1; j < n; ++j)
      {
        for (int k = 0; k < j; ++k)
        {
          gsc[k] = 0;
          for (int l = 0; l < n; ++l)
            gsc[k] += v[n * (l + 1) + j] * v[n * (l + 1) + k];
        }
        for (int k = 0; k < n; ++k)
          for (int l = 0; l < j; ++l)
            v[n * (k + 1) + j] -= gsc[l] * v[n * (k + 1) + l];
        znorm[j] = 0;
        for (int k = 0; k < n; ++k)
          znorm[j] += v[n * (k + 1) + j] * v[n * (k + 1) + j];
        znorm[j] = sqrtf(znorm[j]);
        for (int k = 0; k < n; ++k)
          v[n * (k + 1) + j] /= znorm[j];
      }
        for (int k = 0; k < n; ++k)
          if (znorm[k] > 0) ltot[k] += logf(znorm[k]);
        ++m;
        if (t >= 10100) {
          R lsum = 0;
          R lsum0 = 0;
          int kmax = 0;
          for (int k = 0; k < n; ++k)
          {
            R le = ltot[k] / t;
            lsum += le;
            if (lsum > 0) { lsum0 = lsum; kmax = k + 1; }
          }
          dky = ltot[0] > 0 && kmax < n ? kmax + t * lsum0 / fabsf(ltot[kmax]) : 0;
          break;
        }
      }
      R l = ltot[0] / t;
      if (l > 0)
      {
        out[py][px][0] = fminf(fmaxf(l *   0, 0), 255);
        out[py][px][1] = fminf(fmaxf(l * 128, 0), 255);
        out[py][px][2] = fminf(fmaxf(l * 255, 0), 255);
      }
      else
      {
        out[py][px][0] = fminf(fmaxf(-l * 255, 0), 255);
        out[py][px][1] = fminf(fmaxf(-l * 128, 0), 255);
        out[py][px][2] = fminf(fmaxf(-l *   0, 0), 255);
      }

      #pragma omp critical
      fprintf(stderr, "%8d%8d\r", py, progress++);
    } 
    fwrite(&out[py][0][0], W * 3, 1, stdout);
    fflush(stdout);
  }
  return 0;
}