Image:Lyapunov-fractal.png

From Wikipedia, the free encyclopedia

• Image
• File history
• File links

Size of this preview: 800 × 533 pixels
Full resolution‎ (3,000 × 2,000 pixels, file size: 4.61 MB, MIME type: image/png)

[edit] Summary

Author: --Wickerprints 10:49, 10 July 2006 (UTC)

Summary: Lyapunov fractal for the string bbbbbbaaaaaa, for the growth parameter region . Generated with Mathematica 5

Code used to generate the image (~30 hours on a 1 GHz PPC G4):

In[1]:=
F[a_, b_, x_, n_, k_] := Log[2, Times @@ (Abs[#2(1 - 2#1)] & /@ Transpose[{Drop[FoldList[#2#1(1 - #1) &, 0.5, #], 1 + k*Length[x]], 
Drop[#, k*Length[x]]} &[Flatten[Part[{a, b}, #] & /@ Table[x, {n}]]]])]/((n - k)*Length[x])

In[2]:=
u[m_] := Table[F[a, b, IntegerDigits[222222111111], 200, 50], {b, 2.5, 3.4, (3.4 - 2.5)/(3m/2)}, {a, 3.4, 4, (4 - 3.4)/m}] 

In[3]:=
fu = u[1500]; fmin = Min[fu]; fmax = Max[fu];

In[4]:=
ListDensityPlot[Transpose[fu], Mesh -> False, AspectRatio -> Automatic, ColorFunctionScaling -> False, 
ColorFunction -> (If[# > 0, RGBColor[0, 0, #/fmax], RGBColor[1 - (#/fmin)^(2/3), 1 - (#/fmin)^(1/3), 0]] &)]

[edit] Recomputation with a C program

I (Bernard) have recomputed the image with a C program. The image has been downsampled from a 9000x6000 image.

Note: the original image was a png file with gamma correction included. In order to have the same colors on all displays, I had to set the same gamma value in my file, although I don't like this whole png gamma correction thing.

Following is a simplified version of the code used to generate the image. The complete, optimized code has become messy (but you can ask me if you want it). The main differences are the tests to avoid computing all the rounds when not necessary. In the complete code, the computation for a point is stopped if a convergent cycle is detected (yellow area), or if the mean derivatives become "high enough" (blue area).

Compiles under Linux with gcc -O3 -lm lyapunov. The execution produces a file lyapunov.ppm.

#include <stdio.h>
#include <math.h>

#define CLAMP(x) (((x) > 255) ? 255 : ((x) < 0) ? 0 : (x))

int seq_length = 12;
int seq[12] = {1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0};
double amin = 3.4;
double amax = 4.0;
double bmin = 2.5;
double bmax = 3.4;
int asize = 1000;
int bsize = 1500;
int nmax = 1000;                /* number of rounds */

/* for color generation; somewhat empirical, in order to match Wickerprint's original colors */
double lambda_min = -2.55;
double lambda_max = 0.3959;

int main()
{
  int ai, bi, n, m;
  double a, b, lambda, sum_log_deriv, prod_deriv, r, x, rgb_f[3];
  unsigned char rgb[3];
  FILE *fp;
  
  fp = fopen("lyapunov.ppm", "w");
  fprintf(fp, "P6\n%d %d\n255\n", bsize, asize); /* ppm header */
  
  for (ai = asize-1; ai >= 0; ai--) {
    fprintf(stderr, "\r%d/%d", asize-1-ai, asize);
    for (bi = 0; bi < bsize; bi++) {
      a = amin + (amax-amin)/asize*(ai+0.5);
      b = bmin + (bmax-bmin)/bsize*(bi+0.5);
      x = 0.5;
      /* one round without derivating, so that the value 0.5 is avoided */
      for (m = 0; m < seq_length; m++) {
        r = seq[m] ? b : a;
        x = r*x*(1-x);
      }
      sum_log_deriv = 0;
      for (n = 0; n < nmax; n++) {
        prod_deriv = 1;
        for (m = 0; m < seq_length; m++) {
          r = seq[m] ? b : a;
          /* avoid computing too many logarithms. One every round is acceptable. */
          prod_deriv *= r*(1-2*x); 
          x = r*x*(1-x);
        }
        sum_log_deriv += log(fabs(prod_deriv));
      }
      lambda = sum_log_deriv / (nmax*seq_length);
      
      if (lambda > 0) {
        rgb_f[2] = lambda/lambda_max;
        rgb_f[0] = rgb_f[1] = 0;
      } else {
        rgb_f[0] = 1 - pow(lambda/lambda_min, 2/3.0);
        rgb_f[1] = 1 - pow(lambda/lambda_min, 1/3.0);
        rgb_f[2] = 0;
      }
      rgb[0] = CLAMP(rgb_f[0]*255);
      rgb[1] = CLAMP(rgb_f[1]*255);
      rgb[2] = CLAMP(rgb_f[2]*255);
      fwrite(rgb, 3, 1, fp);
    }
  }
  fprintf(stderr, "\n");
  fclose(fp);
  return 0;
}

[edit] Licensing

This work is licensed under the Creative Commons Attribution-ShareAlike 2.5 License.

File history

Click on a date/time to view the file as it appeared at that time.

• Search for duplicate files
• Edit this file using an external application

  See the setup instructions for more information.

File links