4
$\begingroup$

This question is similar to this one: Properties of King's Dream fractal

The fractal is described here: Softology - Visions Of Chaos 2D Strange Attractor Tutorial

I accidentally modefied the formula to the following:

$$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$$ picture 1

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

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

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

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

Now as you can see, sometimes the fractal is really spread out, sometimes it is squished together, but sometimes the attractor collapses to really small loops or even just a few of points. If you change to values of one attractor continuously to the values of another attractor, the changes are high that in between for a brief moment the attractor collapses and then changes back.

To check at what points that would happen a crude method was applied. The number of individual pixels that were hit was counted. To speed things up the grid was scaled down.

When this value is plotted on a map where $a$ and $b$ are varied and $c$ and $d$ are fixed the following maps are obtained, where $a$ is varied from $-4$ to $4$ on the $x$-axis and $b$ is varied from $-4$ to $4$ on the $y$-axis. Ignore the values of $a$ amd $b$, they are not exact due to rounding errors accumulating. White means high number of pixels hit, black means low number of pixels hit:

map 1

map 2

map 3

Is there anything known about this map? It looks a bit like a fractal. Due to this shape it seems like there is no simple formula to determine the state of the attractor beforhand. Also it seems to be symmetric.

Obviously this method is very crude in determining the state of the fractal. Also it is slow, the last picture took over 4 hours to render. Is there a better method?

For anyone trying to replicate the pictures of the map, here is the processing code for that:

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); 
}

And for anyone trying to test out various values of $a$, $b$, $c$ and $d$ in the actual attractor itself, here is some basic processing code for a black/white image:

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);
}

Thanks a bunch!


These animations show one variable continuously being varied:

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

Here is a crude processing port for rendering the attractor: https://pastebin.com/nVCJipkA


Some efforts were made to draw a bifurcation diagram, which is obviusly hard because it is a 2D-attractor.

A crude solution is simply dropping the $x$-coordinate and instead plotting $a$ on the $x$-axis.

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

Zoom in the first upper little bifurcation tree:

bifurcation 2

Here is a processing script for drawing bifurcations: https://pastebin.com/9ndL2dH3

And here is an animation of the bifurcation plot: https://youtu.be/tRgM-HhtXVk

And some more: https://youtu.be/6WBgtXTPgbA https://youtu.be/ptnNFpxixow

Upon viewing a logarithmic plot it's visible that not every point in the dense regions are hit equally often:

logplot 1

Zoom:

logplot 2

Another plot:

logplot 3

Code for logarithmic colorscale: https://pastebin.com/Xgjkg8tM


Thanks to Claude I could make some Lyapunov plots.

IMPORTANT NOTE! All these maps show the value of the LARGEST Lyapunov exponent. As it turns out that is not always the important thing in determining if this attractor is chaotic or not! What seems to be the case is that THE SUM of the Lyapunov exponents tells us if the system is chaotic or not!

Lyapunov 1

But then I noticed the noise on the left side. At first it seemed just a precision problem, but it was very persistent.

But then, I realised something.

With a random starting point, there was noise. With the set starting point, there was no noise anymore, but there also were hard bounderies.

Random starting position:

Lyapunov 2

Starting position always at $x_0 = -2$ and $y_0 = -2$:

Lyapunov 3

Further investigation shows strange behaviour: There seem to be multiple cycles, and which cycle is reached during the iteration is determined by the starting position.

Animation: https://youtu.be/AME2A2DhSxg

Pictures of cycles (previous point is connected to current point, 50 starting positions overlayed. Three individual cycles can be seen):

Cycles 1

Cycles 2

This poses a problem, because if the starting position determines if the iteration is producing chaos or not, it becomes a difficult task to determine the Lyapunov exponent. Usually given enough iterations the starting position is irrelevant, but not here.

Here is a map that shows where each starting point ends up after the iteration. Be aware that there is a 2-cycle and two 4-cycles, therefore adding a single iteration shifts colors around in the 2- and 4-cycle. As shown before though, being in one of the cycles means being stuck there.

Color-coding of the points: Map 1

Map for starting points: Map 2

(rendered with 100 iterations)

$\endgroup$
  • $\begingroup$ It would be interesting to make a Lypunov map for an a-b-plot or similar; however, I don't know how that would be accomplished, since I only ever made a Lyapunov fractal for the logistic map, which looks like this: en.wikipedia.org/wiki/Lyapunov_fractal Any help on calculating the Lyapunov exponent for a 2D-map like this is appreciated! Thanks! $\endgroup$ – G. Ünther Aug 5 '18 at 18:42
1
$\begingroup$

Here is a Lyapunov fractal plotted for $d = 1, c = -4$ with $a,b \in [-4,4]$. Yellow is stable (largest Lyapunov exponent $\lambda_1 < 0$) and blue is chaotic ($\lambda_1 > 0$). There is some code that computes a fractal dimension estimate from the Lyaponov exponents $\lambda_k$, but the preconditions ($\lambda_1>0 \wedge \lambda_1+\lambda_2<0$) are not always met which means it gives 0 too often to be useful.

Lyapunov fractal

Source code (which is this BASIC code for Henon map converted to C with minor changes):

// 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;
}
$\endgroup$
  • $\begingroup$ Thank you very much for the code! I made a crude port to processing (not finished yet though), that would not have been possible without you. I gotta say though, that code contains some really strange bits. One example being m, declared and used for counting up, but never displayed anywhere else. Another weird one is fminf(fmaxf(l * 0, 0), 255); that would be zero no matter what. But those are just some small things I noticed. $\endgroup$ – G. Ünther Aug 6 '18 at 21:35
  • $\begingroup$ Also, there is a mistake in the formula I just noticed. xnew[1]=sinf(a * x[0]) + d * sinf(b * x[1]); should really be xnew[1]=sinf(a * x[0]) + d * sinf(a * x[1]); Seems like a typo to me. $\endgroup$ – G. Ünther Aug 7 '18 at 20:03
  • $\begingroup$ oh that a vs b typo is bad, I'll fix it and re-render the image shortly. porting to processing (Java based iirc), from C, ported from BASIC, ported from FORTRAN - this code has legs... $\endgroup$ – Claude Aug 7 '18 at 20:24
  • $\begingroup$ It seems like the chosen starting point has a lot of influence over the behaviour of the system. One starting value might produce chaos, another one produces a trajectory that ends in a 2-cycle or even a fixed point (for the same set of parameters $a,\ b,\ c,\ d$). Choosing a random starting point therefore results in noise. How can this be fixed? $\endgroup$ – G. Ünther Aug 12 '18 at 19:37
  • $\begingroup$ @G.Ünther A common method is to continue attractor each time when you change parameter values. For example, if you know any point on attractor (or very close to it) for $p_0$, choose this point as starting for computing attractor at $p_0 + \Delta p$. You have to let some transient to die out and then in principle you get an attractor of perturbed system close to the previous one. $\endgroup$ – Evgeny Aug 31 '18 at 11:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.