The following is the text of a message from Baback Moghaddam to
Tim Wegner concerning the addition of Lyaponav fractal code to Fractint.
Following this is a note on the Lyapunov fractal and some optimizations
posted by Baback (or Bob as he calls himself). Bob gave me explicit
permission to upload these messages to Compuserve. Note: this code has
arrays larger than 64K and won't compile on a PC without changes.
				                 Tim Wegner
------------------------------------------------------------------------
  I think it's best if you get one of your guys to integrate it (I
really don't have a lot of free time to stumble around in you guys'
monumental code and probably would end up making a mess of it all!). An
experienced Fractint programmer who knows the entire thing and the
details of integration, modularity, menus, etc., would be best. Although
I wouldn't mind having the official source for 16.11!  ;-)

 Naturally, I'll be happy to explain/help out with my end of the code. 
I'm appending the C-code at then end of this mail. It should compile on
a DEC3100 running Ultrix, but then again it's real plain C-code so you
shouldn't have any trouble compiling it on different platforms. The
input options are pretty self-explanatory. You have the choice of B&W
output (PGM format), which is always best since you colormap it
later on, or else you can output PPM and you are asked to provide a
colormap (256 lines of ASCII triplets [0-255], like Fractint's but
without any text!).

 The magnitude/contour option means if you want a regular continuous
gray map or if you want to do a Mod-N pseudo-contouring (this gives
pretty nice B&W images of the "binary-decomposition" type but is not
suited for Markus fractals). As you might guess, this code is a rough
and dirty hack for the "backbone" of most of my iterated fractal
programs.

 If you compile it ok, test it out with my images, they used a
range of 3-4 on both axes and the chroma.map from Fractint. And oh, for
personal reasons I do not like the original technique of collapsing all
positive exponents into one color (my way you can see gradations of
chaos, which was of considerable interest to me when I was studying this
and other systems). I think this also has better visual appeal for Fractint
enthusiasts. A sample input file is shown following the code.    

 Good luck, and let me know if you have any questions.

--
Bob

--------------------------- CUT HERE ---------------------------------

/* ================================================================= */
/*  PROGRAM: lyap.c                Lyapunov Space of the             */
/*  Version: 1.1   09/02/91           Logistic Parabola              */
/*  By:      Baback Moghaddam      with Periodic Forcing             */
/* ================================================================= */

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

#define PI 3.1415926535
#define MAX 512
#define ABS(a) ((a)>0 ? (a):-(a))
#define SGN(a) ((a)>0 ?  1 : -1 )
#define FMIN 1e20
#define FMAX 1e-20

float v[MAX][MAX];                   /* Pixel Map             */
int colors[256][3];                /* Color Palette         */

/* ================================================================== */
/* ================================================================== */

main()
{

 FILE *fopen(),*fin,*fout;            /* Input & Output File Pointers */
 char in_file[50],out_file[50];       /* Input & Output File Names    */
 int Nx,Ny,No,Nc,L,P[30],cmode;
 register int i,j,k;
 int out,type,mod=2,cnt,c1,c2,c3,count;
 float x,y,dx,dy,lambda,f,prod,total;
 float x1,x2,y1,y2,vmin,vmax,z,lim,R[3];
 char s[30];
 
 printf("\n");
 printf("[*] Boundaries (Xmin,Xmax,Ymin,Ymax)  : ");
 scanf("%f %f %f %f",&x1,&x2,&y1,&y2);
 printf("[*] Spatial Resolutions      (Nx,Ny)  : ");
 scanf("%d %d",&Nx,&Ny);
 printf("[*] Convergence & Computaion Its      : ");
 scanf("%d %d",&No,&Nc);
 printf("[*] Period of Parameter Cycling       : ");
 scanf("%d",&L);
 printf("[*] Binary Parameter Sequence  (0,1)  : ");
 for(i=0;i<L;i++)
  scanf("%d",&P[i]);
 printf("\n");
 printf("[*] Magnitude or Contour Plot  (0,1)  : ");
 scanf("%d",&type);
 if(type==1) {
  cmode=0;
  printf("\n    > Contour Modulus   (0-255) : ");
  scanf("%d",&mod);
 }
 printf("\n");
 if(type==0) {
  printf("[*] B&W or Color Output        (0,1)  : ");
  scanf("%d",&cmode);
  if(cmode==1) {
   printf("[*] Color Map Filename                : ");
   scanf("%s",in_file);
  }
 }
 if(cmode==0)
  printf("[*] PGM Output File                   : ");
 if(cmode==1)
  printf("[*] PPM Output File                   : ");
 scanf("%s",out_file);
 printf("\n");
 
 fin=fopen(in_file,"r");             /* read color map from file */
 for(i=0;i<256;i++) {
  fscanf(fin,"%d %d %d",&c1,&c2,&c3);
  colors[i][0]=c1;
  colors[i][1]=c2;
  colors[i][2]=c3;
  }
 fclose(fin);

 dx=(x2-x1)/(Nx-1);
 dy=(y2-y1)/(Ny-1);
 vmin=1e30;
 vmax=0-1e30;
 for(i=0;i<Nx;i++)
 for(j=0;j<Ny;j++) {
   R[0]=x1+i*dx;
   R[1]=y1+j*dy;
   x=0.5;
   prod=1;
   count=0;
   total=0;
   for(k=0;k<No;k++)
    x=R[P[k%L]]*x*(1-x);
   for(k=No;k<No+Nc;k++) {
    x=R[P[k%L]]*x*(1-x);          /* Mapping Function  Xn+1=f(Xn)    */
    f=ABS(R[P[k%L]]*(1.0-2.0*x)); /* Derivative of f(Xn), ie. f'(Xn) */
    if (f>0) {
     prod*=f;
     if(prod>FMAX || prod<FMIN) {
      total+=log((double) prod);
      prod=1.0;
     }
     count++;
    }
   }
   total+=log((double) prod);
   if (count>0)
    lambda=total/(log((double) 2.0)*count);
   else
    lambda=vmin;
   f=lambda; 
   v[i][j]=f;
   if (f<vmin) vmin=f;
   if (f>vmax) vmax=f;
 }

 printf("(*) Min = %f \n(*) Max = %f\n",vmin,vmax);
 printf("    .... Rescaling to (0-255).\n\n");
 
 /* write intensity map to output file */
  
 fout=fopen(out_file,"w");
 if(cmode==0)
  fprintf(fout,"P5\n%d %d\n%d\n",Ny,Nx,255);
 else
  fprintf(fout,"P6\n%d %d\n%d\n",Ny,Nx,255);
 z=255.0/(vmax-vmin);
 for(j=Ny-1;j>=0;j--) 
 for(i=0;i<Nx;i++) {
  out=(int) z*(v[i][j]-vmin)+0.5;
  k=0;
  if ((out % mod)==0)  k=255;
  if (type==0)
   if(cmode==0)
    putc(out,fout);
   else {
    putc(colors[out][0],fout);
    putc(colors[out][1],fout);
    putc(colors[out][2],fout); }
  else
    putc(k,fout);
 }

 fclose(fout);
 return 0;
}

------------------------------- CUT HERE -----------------------------

[*] Sample input (for image lyap6.gif at the ftp site)

----- cut ----
512 512
600 6000
8
0 1 1 0 1 0 0 1
0
0
lyap6.pgm
----- cut ----

you can pipe this directly into the program. The periodic sequence is
the famous Morse-Thue sequence which is intimately related to the
bifurcation dynamics of the logistic parabola and believe it or not, the
Mandelbrot Set! 

(Below is the alt.fractals posting by Baback Moghaddam)
---------------------------------------------------------------------
 Being an occasional reader of SA (and even less frequently of Dewdney's
column) I was intrigued by all the articles referring to his September
column. So I finally went out and bought a copy to see what all the
"hoopla" was about. Naturally, the appeal of the images was enough to
make me spend the entire Labor Day weekend exploring this and other 
nonlinear maps with periodic forcing.

 I'd like to point out few typographical errors in the article which
were probably immediately apparent to those versed in the logicstic
parabola. Starting with the last paragraph on page 178 and onto the
second paragraph on page 179, the author erroneously associates the
parameter values 2.0 and 2.45 with the boundaries of Period-2 limit
cycle. The values are 3.0 and 3.45 (the latter is of course only
approximate). Therefore simply add 1.0 to these values (the accumulation
point of the period-doubling cascade Rinf=3.5699.. is however correct).

 A similar error appears on the next to last paragraph on page 179,
where he encourages the readers to try computing the Lyapunov exponent
for values or r=2.0 and r=3.0, claiming that they correspond to negative
(stability) and positive (possibly chaotic) exponents. The fact of the
matter is both parameters correspond to stable behaviour (with r=3.0
being slightly less stable - ie. less negative). In fact r=3.0 is,
loosely speaking "marginally stable" since it corresponds to the
"pitchfork" bifurcation (where the previously stable fixed point goes
unstable along with the simultaneous apperance of two marginally stable
points, marking the "birth" of the Period-2 cycle.

 I think what the author intended (and I realize how dangerous it is to
presume to know what someone else is thinking) was to suggest values 
like 3 and 4 (or even better, 2.5 and 4) which should give you values
like -1 and +1.

 A further point of clarification is the origin of the formula |r-2rx|.
Like me, you may have wondered how this "rate of growth" was derived. In
fact, it is simply the derivative of the mapping function Xn+1=f(Xn),
where

    f(Xn) = r(1-x)x 

which gives
     
    f'(Xn) = r - 2rx

This derivative characterizes how small perturbations to Xn will grow with
each iteration. For a discrete system, they scale geometrically
according to: error = slope^n for n iterations. This characterizes the
divergence/convergence of two trajectories with slightly different
initial conditions and hence quantifies the "sensitivity to initial
conditions". Slopes greater than 1.0 (which have a *positive* log) will
diverge while slopes with magnitude less than 1.0 (which have a
*negative* log) will die out. The Lyapunov exponent is therefore an
orbital "average" obtained by 

              N 
             ___
   L  = 1/N  \    log |f'(Xn)|
             /__ 
             n=1        

where positive exponents L signify an overall average tendency for
separation of nearby orbits (a potential indicator of chaotic dynamics).


 In regards to implementation of the above formula, there are several
important modifications. In his posting Dave Platt, correctly pointed
out that since the sum of logarithms is equal to the logarithm of
products, Dewdney's formula can be altered to give:

 Total=1
 for i=1:N
  x=r*x*(1-x)
  Total=Total*|r-2*r*x|
 end
 Lyap=log(Total)/log(2);
 Lyap=Lyap/N

Dave Platt goes on to point out that the *product* Total in this new
implementation can easily grow out of bounds, so he restricts it with:

 Prod=1
 Total=0
 for i=1 to N
  x=r*x*(1-x)
  Prod=Prod*|r-2*r*x|
  if Prod>1e20
   Total=Total+log(Prod)/log(2)
   Prod=1
  end
 end
 Total=Total+log(Prod)/log(2)


which gives a much more robust and efficient algorithm (efficient since
the computation of the logarithm is reduced significantly and depending
on your machine can make a worthwile difference).a

 In the time-honored tradition of "adding your two-cents worth", I'd
like to suggest the following modification to the above algorithm. Just
as the value of Prod is likely to grow out of bounds, it can also
converge on zero (especially for low values of r). Therefore one should
also limit it to a lower bound. Another important modification is to
trap cases in which a single zero value for |r-2rx| will nullify all
further computing (since it *multiplies* Prod). My modification to Dave
Platt's method is the following (slightly more robust) version:

 Prod=1
 Total=0
 Count=0
 for i=1 to N
  x=r*x*(1-x);
  slope=|r-2*r*x|
  if slope>0
   Prod=Prod*slope
   if Prod>1e20 OR Prod<1e-20
    Total=Total+log(Prod)/log(2)
    Prod=1
   end
  Count=Count+1
 end
 Total=Total+log(Prod)/log(2)
 if Count>0
  Lyap=Total/Count
 else
  Lyap= - Infinity
 end


The slope>0 loop insures that occasional zero slope values do not
nullify the computation. At the end if Count=0, which means all slopes
were zero (for example for a so-called "Superstable Orbit"), we know
that the exponent is in fact - infinity, otherwise we simply average by
the number of non-zero values (ie. Count).
   
 Armed with this new algorithm, I strongly suggest to the avid
experimenter to compute the Lyapunov plot for the simple Logistic
function (ie. without the periodic forcing). Vary r from 2.5 to 4 and
compute an L for each r, then plot the function L(r). It is highly
instructive to line up this plot against the bifurcation diagram in
order to illustrate the important "meaning" of the exponent. When this
is done, one can clearly see that stable regions (Limit Cycles)
correspond to negative values. In particular, the occasional "dips" in
L(r) correspond to the periodic windows in the chaotic region.

 I also suggest that you compute the two-dimensional "structural
stability" map of the forced systems with the above modification. In
particular take a look at the system with the forcing "aabab" for r values
in the range 1 to 4 (on both axes). Note that as Dave Platt points out
Dewdney's description of the colormap used is incomplete.

 Anyway, I hope that this rather long exposition will be of help to fellow
experimenters out there. Have fun!

------------------------------------------------------------------------
Baback Moghaddam                     |  Logistic <- logistique (French),
Image Sensing Lab / C3I Center       |  from the "lodgement" of troops!
Dept. of Electrical & Computer Engr. -----------------------------------
George Mason University              |  E-mail: bmoghad@sitevax.gmu.edu
Fairfax Va. 22030                    |     Tel: (703) 993-1697
------------------------------------------------------------------------