Conference rusure::math

1037.0. "Henon graphs anyone??" by MYVAX::ROBERUC () Tue Mar 14 1989 18:57

    { 
    Hi.

    Below is program that prints a scaled plot of a Henon equation.
    This is an equation that produces Henon Attractor (or the squashed 
    banana that is in many Chaos books. Below are the basic equations.)
    
                        /  \2
    X     = 1 + Y  - A ( X  )
     n+1         n      \ n/

    Y     = B Y
     n+1       n
    
    
    Some very interesting pictures can be drawn from some very simple
    A, B, x, y, and number of Iterates to draw values.

     A     B     x    y   Iterates
    1.4   .3     0    0   5000 +            : Henon attractor (banana)
    .005  -1     0    0   5000              : a circle (I think)
    .01   -.999  75   5   5000              : a butterfly (??)
    .005  -.999  75   5   5000              : 4 circles
    -.4   -.999  1.5  0   7000              : a paw print
    -.6   -.999  0    0   5000              : a head to an egg (??)
    .5    -.99   1.8  0   2000              : a spiral ??
    .9    -.999  0    0   2000              : a triangular spiral
    
    Its all very pretty and nice. I especially like the one butterfly
    shape. But could someone tell me why the second one draws a very
    near to a perfect circle??
    
    Rich
    }
    
    [INHERIT ('SYS$LIBRARY:STARLET')] PROGRAM Loop (Input, Output);

VAR
  I : UNSIGNED;
  I_1, X, Y : INTEGER;
  S_1, S_2,
  Factor, X_Hi, X_Lo, Y_Hi, Y_Lo, A_1, B_1, X_1, X_2, Y_1, Y_2 : QUADRUPLE;
  Bel, St, Dcs, Esc, C : CHAR;
  Str : VARYING [80] OF CHAR;
  F : TEXT;                      
BEGIN
  Bel := CHR (7);
  Esc := CHR (27);
  Dcs := CHR (144);
  St  := CHR (156);
         
  WRITE ('A: ');  READLN (A_1); 
  WRITE ('B: ');  READLN (B_1);

  I_1 := 500;
  X_1 := 0.0; 
  Y_1 := 0.0; 
  WRITE ('X: ');  READLN (X_1); S_1 := X_1;
  WRITE ('Y: ');  READLN (Y_1); S_2 := Y_1;
  WRITE ('I: ');  READLN (I_1);

  WRITE ('Graph, File, or Points? [G, F, or P]: ');
  READLN (C);
  IF (C = 'f') THEN C := 'F';
  IF (C = 'g') THEN C := 'G';
  IF (C = 'p') THEN C := 'P';

  IF (C = 'G') OR 
     (C = 'F') THEN BEGIN
    X_Hi := X_1; X_Lo := X_1;
    Y_Hi := Y_1; Y_Lo := Y_1;
    FOR I := 1 TO I_1 DO BEGIN
      X_2 := 1 + Y_1 - (A_1 * (X_1 ** 2));
      Y_2 := X_1 * B_1;
      X_Hi := MAX (X_Hi, X_2);
      X_Lo := MIN (X_Lo, X_2);
      Y_Hi := MAX (Y_Hi, Y_2);
      Y_Lo := MIN (Y_Lo, Y_2);
      X_1 := X_2;
      Y_1 := Y_2;
    END;
    Factor := 478 / MAX (X_Hi - X_Lo, Y_Hi - Y_Lo);
  END;

  X_1 := S_1;
  Y_1 := S_2;

  IF (C = 'G') THEN WRITELN (Esc + 'Pp;S(E); S[0,0]');

  IF (C = 'F') THEN BEGIN
    WRITE ('Filename: ');
    READLN (Str);
    OPEN (F, Str, HISTORY:=NEW);
    REWRITE (F);
    WRITELN (F, Esc + 'Pp;S(E); S[0,0]');
  END;

  FOR I := 1 TO I_1 DO BEGIN
    IF (C = 'G') THEN BEGIN
      X := ROUND ((X_1 - X_Lo) * Factor);
      Y := ROUND ((Y_1 - Y_Lo) * Factor);
      WRITELN ('p[', X:3, ',', Y:3, ']');
      WRITELN ('v[', X:3, ',', Y:3, ']');
    END;
    IF (C = 'F') THEN BEGIN
      X := ROUND ((X_1 - X_Lo) * Factor);
      Y := ROUND ((Y_1 - Y_Lo) * Factor);
      WRITELN (F, 'p[', X:3, ',', Y:3, ']');
      WRITELN (F, 'v[', X:3, ',', Y:3, ']');
    END;
    IF (C = 'P') THEN BEGIN
      WRITELN (I, X_1, Y_1); 
    END;

    X_2 := 1 + Y_1 - (A_1 * (X_1 ** 2));
    Y_2 := X_1 * B_1;

    X_1 := X_2;
    Y_1 := Y_2;
  END;

  IF (C = 'F') THEN WRITELN (F, St);
  IF (C = 'G') THEN WRITELN (St);
END.

    Two asides:
    
    1: That program needs to be compiled in Pascal...
    2: It draws using ReGIS escape sequences
    
    Rich

�    But could someone tell me why the second one draws a very
�    near to a perfect circle??
    
    The Henon equations are

	x[k+1] = 1+y[k]-a*x[k]^2
	y[k+1] = b*x[k]

    If b = -1, and a is small then there will be a fixed point near (.5,-.5),
    given by a root of the equation

	x = 1-x-a*x^2

    and the recurrance will trace out an approximate circle about the fixed
    point.

    If you translate to the origin and scale the recurrance can be
    rewritten in the form

	x[k+1] =  y[k]-a*(1+x)^2
	y[k+1] = -x[k]

    This is a 90 degree clockwise rotation with a small nonlinear
    perturbation.  For the recurrance to stay "on" some closed circular
    shaped curve, the implicit equation for the curve, F(x,y) = 0, should be
    an identity upon the subtitution F(y-a*(1+x)^2, -x) = 0.

    For example a circle of radius r is F(x,y) = x^2+y^2-r^2 = 0.
    If we try substiting y-a*(1+x)^2 for x and -x for y we get

	F[k](x,y) = x^2 + y^2 - r^2 = 0

	F[k+1](x,y) = (y-a*(1+x)^2)^2 + x^2 - r^2

		    = y^2 - a*(1+x)^2*(y - a*(1+x)^2) + x^2 - r^2

	"error"	    = -a*(1+x)*(y-a*(1+x)^2)

    So there is a small error after one iteration.

    But suppose the real curve was expressed by adding a power series
    correction, E(x,y) to the equation for a circle.

	E(x,y) = E0 + Ex*x + Ey*y + Exx*x^2 + Exy*x*y + Eyy*y^2 + ...

    with as-yet undetermined coefficients.

    If we substitue for x and y in the series from the recurrance, expand
    out the powers of x and y, add in the "error" above and collect
    powers of x, y, and their products, there will result a set of
    equations relating the unknown coefficients for the power series, and
    it should be possible to solve for them.  Moreover, since high powers
    of "a" will occur, the coefficients will decrease and you should be
    able to show convergence of the series.

    It looks like the details are messy but I think this will work.

    Another idea along the same lines would be to assume the
    unknown function has a Fourier series expansion, plug in from the
    recurrance (involving trigonometric polynomials), match coefficients,
    and again get a solution.

    - Jim    

 Hi,

 I have written a DECwindows program that plots the difference equations
 that I believe are attributable to Michel H�non and Carl Heils.  This
 does not appear to be the problem disussed in relation to H�non and
 Heiles on page 149 of "Chaos" by James Gleick, but this is one I worked
 with in the old days:

                                      2
       X     = X  * COS(A) - ( Y  -  X  ) * SIN(A)
        N+1     N               N     N

                                      2
       Y     = X  * SIN(A) + ( Y  -  X  ) * COS(A)
        N+1     N               N     N

 My program allows setting various initial conditions, colors, the origin,
 and a scaling factor.  If you would like to try it look for
       
                      SOADC::SYS$PUBLIC:HENON.BACKUP_SAVESET

  You might have to link it if you don't use VMS 5.1.  There are example
  plots in the files henon*.six and henon*.ps.

 Regards,
 David

 [This reply was edited to included updated information.]