Conference rusure::math

819.0. "Lorenz Attractor program" by CADM::ROTH (If you plant ice you'll harvest wind) Fri Jan 22 1988 06:58

    Here is a simple program that shows the Lorenz attractor on the
    Vaxstation, as an illustration of chaotic dynamics and non periodic
    motion.  It asks for a time step (.01 is a good number), and a starting
    point in 3 space, (.01,.01,.01) is OK, and plots the y and z projections
    of the resulting curve.  Answering y or n shows line segments or not
    connecting the successive solution points.  The resulting curve
    lies mostly in a plane, I didn't bother to allow rotating the curve
    in 3 dimensions to show this but it could easily be added.

    Note that the equations are really very simple; most of the code is
    concerned with the numerical integration, and I probably used a fancier
    method than was really necessary.

    - Jim

/*
 *  Show Lorenz attractor, using RK4(5) integration
 */

#include stdio
#include math
#include ssdef
#include descrip

#include "vsi$library:uisentry.h"
#include "vsi$library:uisusrdef.h"

#define DEFAULT_ATTR 0
#define OVER_ATTR 1
#define ERAS_ATTR 2

static float wc_xmin, wc_ymin, wc_xmax, wc_ymax;
static float vd_width, vd_height;
static int vd_id, wd_id;

static $DESCRIPTOR(ws_devname, "SYS$WORKSTATION");
static $DESCRIPTOR(vd_title, "Lorenz Attractor");

static setup_window(size, range)
double size, range;
{
  vd_width = size;
  vd_height = size;

  wc_xmin = -range;
  wc_ymin = -range;
  wc_xmax =  range;
  wc_ymax =  range;

  vd_id = uis$create_display(&wc_xmin, &wc_ymin, &wc_xmax, &wc_ymax,
			     &vd_width, &vd_height);

  uis$disable_display_list(&vd_id);

  wd_id = uis$create_window(&vd_id, &ws_devname, &vd_title,
			    &wc_xmin, &wc_ymin, &wc_xmax, &wc_ymax,
			    &vd_width, &vd_height);

  uis$set_writing_mode(&vd_id, &DEFAULT_ATTR, &OVER_ATTR, &UIS$C_MODE_OVER);
  uis$set_writing_mode(&vd_id, &DEFAULT_ATTR, &ERAS_ATTR, &UIS$C_MODE_ERAS);
}

static clnup_window()
{
  uis$delete_display(&vd_id);
}

int ask_yn(query)
char *query;
{
  static response[132];

  printf(query);
  gets(response);
  return response[0] == 'Y' || response[0] == 'y';
}

double ask_d(query)
char *query;
{
  static response[132];
  double answer;

  while(1) {
    printf(query);
    gets(response);
    if (sscanf(response, "%lf", &answer) > 0) return answer;
    }
}

main(ac,av)
int ac;
char *av[];
{
  double x0, y0, z0, x1, y1, z1;
  double x2, y2, z2, x3, y3, z3;
  double x4, y4, z4, x5, y5, z5;
  double dx0, dy0, dz0, dx1, dy1, dz1;
  double dx2, dy2, dz2, dx3, dy3, dz3;
  double dx4, dy4, dz4, dx5, dy5, dz5;
  double dt;
  double s, r, b;
  static double b10 =   2.0/9.0;
  static double b20 =   1.0/12.0,  b21 =    1.0/4.0;
  static double b30 =  69.0/128.0, b31 = -243.0/128.0, b32 = 135.0/64.0;
  static double b40 = -17.0/12.0,  b41 =   27.0/4.0;
  static double b42 = -27.0/5.0,   b43 =   16.0/15.0;
  static double b50 =  65.0/432.0, b51 =   -5.0/16.0,  b52 =  13.0/16.0;
  static double b53 =   4.0/27.0,  b54 =    5.0/144.0;
  static double c0 = 47.0/450.0, c1 = 0.0,      c2 = 12.0/25.0;
  static double c3 = 32.0/225.0, c4 = 1.0/30.0, c5 = 6.0/25.0;
  int i, j, k, n;
  int showlines = 0;
  double y0_old, z1_old;
 
  s = 10.0; b = 8.0/3.0; r = 28.0;

  dt = ask_d("enter dt: ");
  n = ask_d("enter number of iterations: ");

  x0 = ask_d("enter x0: ");
  y0 = ask_d("enter y0: ");
  z0 = ask_d("enter z0: ");

  showlines = ask_yn("show connecting lines? ");

  setup_window(30.0, 30.0);

  uis$plot(&vd_id, &OVER_ATTR, &0.0, &0.0);

  y0_old = y0; z1_old = z0-28.0;

  for (i = 0; i < n; i++) {
    dx0 = s*(y0-x0);
    dy0 = r*x0-y0-x0*z0;
    dz0 = x0*y0-b*z0;

    x1 = x0+dt*b10*dx0;
    y1 = y0+dt*b10*dy0;
    z1 = z0+dt*b10*dz0;

    dx1 = s*(y1-x1);
    dy1 = r*x1-y1-x1*z1;
    dz1 = x1*y1-b*z1;

    x2 = x0+dt*(b20*dx0+b21*dx1);
    y2 = y0+dt*(b20*dy0+b21*dy1);
    z2 = z0+dt*(b20*dz0+b21*dz1);

    dx2 = s*(y2-x2);
    dy2 = r*x2-y2-x2*z2;
    dz2 = x2*y2-b*z2;

    x3 = x0+dt*(b30*dx0+b31*dx1+b32*dx2);
    y3 = y0+dt*(b30*dy0+b31*dy1+b32*dy2);
    z3 = z0+dt*(b30*dz0+b31*dz1+b32*dz2);

    dx3 = s*(y3-x3);
    dy3 = r*x3-y3-x3*z3;
    dz3 = x3*y3-b*z3;

    x4 = x0+dt*(b40*dx0+b41*dx1+b42*dx2+b43*dx3);
    y4 = y0+dt*(b40*dy0+b41*dy1+b42*dy2+b43*dy3);
    z4 = z0+dt*(b40*dz0+b41*dz1+b42*dz2+b43*dz3);

    dx4= s*(y4-x4);
    dy4 = r*x4-y4-x4*z4;
    dz4 = x4*y4-b*z4;

    x5 = x0+dt*(b50*dx0+b51*dx1+b52*dx2+b53*dx3+b54*dx4);
    y5 = y0+dt*(b50*dy0+b51*dy1+b52*dy2+b53*dy3+b54*dy4);
    z5 = z0+dt*(b50*dz0+b51*dz1+b52*dz2+b53*dz3+b54*dz4);

    dx5= s*(y5-x5);
    dy5 = r*x5-y5-x5*z5;
    dz5 = x5*y5-b*z5;

    x0 += dt*(c0*dx0+c1*dx1+c2*dx2+c3*dx3+c4*dx4+c5*dx5);
    y0 += dt*(c0*dy0+c1*dy1+c2*dy2+c3*dy3+c4*dy4+c5*dy5);
    z0 += dt*(c0*dz0+c1*dz1+c2*dz2+c3*dz3+c4*dz4+c5*dz5);

    z1 = z0-28.0;
    if (showlines) {
      uis$plot(&vd_id, &OVER_ATTR, &y0_old, &z1_old, &y0, &z1);
      y0_old = y0; z1_old = z1;
      }
    else  uis$plot(&vd_id, &OVER_ATTR, &y0, &z1);
    }

  printf("Final [x, y, z] = [%16.12f %16.12f %16.12f]\nAt t = %f\n",
	 x0, y0, z0, n*dt);

  { char dummy[132]; gets(dummy); }

  clnup_window();
}