Chaos Theory - Lorenz Attractor on the Oscilloscope

Expanded on the X-Y oscilloscope control idea from my last project and have programmed the arduino to display a Lorenz strange attractor on an Oscilloscope.. The poor arduino does struggle with the calculations but the chaotic effect can be seen quite nicely - If I got a Maple would it help? I heard the processing power is much stronger than the arduino and it's shield compatible? Is there something I'm doing wrong? How could I get it to run faster?

/* 
Chaos Theory on your Scope! 
Displays the lorenz attractor on your Oscilloscope.
Author: Chris Gozzard (cgozzard@yahoo.com)
Created: 25-5-13
License: This code CC-BY-SA 3.0 and is unsupported.
(see creativecommons.org/licenses for info)

the following circuit is on both PWM ports (5 and 6)

	       R
PWM OUT ----/\/\/\-----+------------ OUTPUT
		       |
		      === C
                       |
                      GND

R = 10k
C = 0.1uF		

Many thanks to John M. De Cristofaro for the original idea (based on his Oscilloscope Christmas Tree circuit)  
*/
float x_pos;         
float y_pos; 
int x_out = 5;
int y_out = 6;
float x = 0.7;
float y = 0;
float z = 0;

float alpha  = 15.6;
int beta   = 28; 
float m0     = -1.143;
float m1     = -0.714;
float xdot;
float ydot;
float zdot;
float h;
float dt = 0.01;

void setup()  { 

  pinMode(x_out, OUTPUT);
  pinMode(y_out, OUTPUT);
  // this next section changes the PWM frequency - don't mess with it!
    TCCR0A = (	1<<COM0A1 | 0<<COM0A0 |		// clear OC0A on compare match (hi-lo PWM)
		1<<COM0B1 | 0<<COM0B0 |		// clear OC0B on compare match (hi-lo PWM)
		1<<WGM01  | 1<<WGM00);		// set PWM lines at 0xFF

  TCCR0B = (	0<<FOC0A | 0<<FOC0B |		// no force compare match
		0<<WGM02 |			// set PWM lines at 0xFF
		0<<CS02	 | 0<<CS01 |		// use system clock (no divider)
		1<<CS00 );

  TIMSK0 = (	0<<OCIE0B | 0<<TOIE0 |
		0<<OCIE0A ); 
} 

void loop()  {

xdot = alpha*(y-x-h);
ydot = x - y + z;
zdot  = (-beta)*y;
h = (m1*x)+(0.5*(m0-m1))*(abs(x+1)-abs(x-1));
x = x + (xdot*dt);
y = y + (ydot*dt);
z = z + (zdot*dt);

  
x_pos = 128+(50 * x);
y_pos = 128+(50 * y);

analogWrite(x_out, x_pos);
analogWrite(y_out, y_pos);
}

could floating point be an issue with something like this??

removing floating point math would surely help.

make all numbers long.

multiply values with 100, and divide them at appropriate places with correction values.

I'll be back

Replaced some float math with int math, but it was tricky not to "loose the attractor".
Still I got loop down from ~240 micros to ~180 micros ==> -60 so about 25% faster.

give it a try

/* 
 Chaos Theory on your Scope! 
 Displays the lorenz attractor on your Oscilloscope.
 Author: Chris Gozzard (cgozzard@yahoo.com)
 Created: 25-5-13
 License: This code CC-BY-SA 3.0 and is unsupported.
 (see creativecommons.org/licenses for info)
 
 the following circuit is on both PWM ports (5 and 6)
 
 	       R
 PWM OUT ----/\/\/\-----+------------ OUTPUT
 		       |
 		      === C
 |
 GND
 
 R = 10k
 C = 0.1uF		
 
 Many thanks to John M. De Cristofaro for the original idea (based on his Oscilloscope Christmas Tree circuit)  
 */

int x_pos;         // can be int
int y_pos;          // can be int
int x_out = 5;
int y_out = 6;

float x = 0.70;
float y = 0;
float z = 0;

float dt = 0.01;
float alpha  = 15.60 *dt;   // precalculate the *dt
float beta   = -28*dt;    // precalculate the * dt and the - sign
float m0     = -1.143;
float m1     = -0.714;

float m2 = 0.5 * (m0-m1);  // precalculated const

float xdot;
float ydot;
float zdot;

float h;


void setup()  
{ 
  Serial.begin(115200);
  pinMode(x_out, OUTPUT);
  pinMode(y_out, OUTPUT);
  
  // this next section changes the PWM frequency - don't mess with it!
   TCCR0A = (	1<<COM0A1 | 0<<COM0A0 |		// clear OC0A on compare match (hi-lo PWM)
   1<<COM0B1 | 0<<COM0B0 |		// clear OC0B on compare match (hi-lo PWM)
   1<<WGM01  | 1<<WGM00);		// set PWM lines at 0xFF
   
   TCCR0B = (	0<<FOC0A | 0<<FOC0B |		// no force compare match
   0<<WGM02 |			// set PWM lines at 0xFF
   0<<CS02	 | 0<<CS01 |		// use system clock (no divider)
   1<<CS00 );
   
   TIMSK0 = (	0<<OCIE0B | 0<<TOIE0 |
   0<<OCIE0A ); 
   
} 

void loop()  
{
    float a = y-x;      // precalculate for reuse
    xdot = alpha * (a-h); 
    ydot = (z-a) * dt; 
    zdot = beta * y;

    h = (m1 * x) + m2 * (abs(x + 1) - abs(x - 1)); 

    x = x + xdot;  
    y = y + ydot;
    z = z + zdot;

    x_pos = 128 + (int)(50*x); 
    y_pos = 128 + (int)(50*y);

    analogWrite(x_out, x_pos);
    analogWrite(y_out, y_pos);
}

"simplified" the abs formula after some tinkering..

loop takes now ~146micros() so almost 40% faster

/* 
 Chaos Theory on your Scope! 
 Displays the lorenz attractor on your Oscilloscope.
 Author: Chris Gozzard (cgozzard@yahoo.com)
 Created: 25-5-13
 License: This code CC-BY-SA 3.0 and is unsupported.
 (see creativecommons.org/licenses for info)
 
 the following circuit is on both PWM ports (5 and 6)
 
 	       R
 PWM OUT ----/\/\/\-----+------------ OUTPUT
 		       |
 		      === C
 |
 GND
 
 R = 10k
 C = 0.1uF		
 
 Many thanks to John M. De Cristofaro for the original idea (based on his Oscilloscope Christmas Tree circuit)  
 */

int x_pos;         
int y_pos; 
int x_out = 5;
int y_out = 6;

float x = 0.70;
float y = 0;
float z = 0;

float dt = 0.01;
float alpha  = 15.60*dt;
float beta   = -28*dt; 
float m0     = -1.143;
float m1     = -0.714;

float m2 = 0.5*(m0-m1);
float m3 = 2 * m2;

float xdot;
float ydot;
float zdot;

float h;


void setup()  
{ 
  Serial.begin(115200);
  pinMode(x_out, OUTPUT);
  pinMode(y_out, OUTPUT);

  // this next section changes the PWM frequency - don't mess with it!
   TCCR0A = (	1<<COM0A1 | 0<<COM0A0 |		// clear OC0A on compare match (hi-lo PWM)
   1<<COM0B1 | 0<<COM0B0 |		// clear OC0B on compare match (hi-lo PWM)
   1<<WGM01  | 1<<WGM00);		// set PWM lines at 0xFF
   
   TCCR0B = (	0<<FOC0A | 0<<FOC0B |		// no force compare match
   0<<WGM02 |			// set PWM lines at 0xFF
   0<<CS02	 | 0<<CS01 |		// use system clock (no divider)
   1<<CS00 );
   
   TIMSK0 = (	0<<OCIE0B | 0<<TOIE0 |
   0<<OCIE0A ); 
} 

void loop()  
{
  //  unsigned long start = micros();
  //
  //  for (int i=0; i < 10000; i++)
  //  {
  float a = y-x;
  xdot = alpha * (a-h); 
  ydot = (z-a) * dt; 
  zdot = (beta) * y;

  // the formula esp the abs() function part can be replaced by the code below.
  // draw a graph of y = (abs(x+1) - abs(x-1))  to get the transition
  // m3 is a precalculated value
  if (x <= -1) h = (m1*x) - m3; 
  else if (x >= 1) h = (m1*x) + m3; 
  else h = (m1 + m3) * x; 

  // original formula
  // h = (m1*x) + m2 * (abs(x + 1) - abs(x - 1));   

  //Serial.print(h);
  //Serial.print("\t");

  x = x + xdot;  
  y = y + ydot;
  z = z + zdot;

  x_pos = 128 + (int)(50*x);
  y_pos = 128 + (int)(50*y);

  //Serial.print(x_pos);
  //Serial.print("\t");
  //Serial.println(y_pos);
  analogWrite(x_out, x_pos);
  analogWrite(y_out, y_pos);
  //  }
  //  Serial.println((micros()-start)/10000.0);
}

Please give it a try

Sorry but it does not look like a Lorentz attractor to me, is is supposed to look like a butterfly. I have plotted 3D versions in the past.

Another thing you could optimize is the analogWrite().
As it has fixed pin numbers, the math that's involved with the mapping on the HW ports can be precalculated.
Drawback is that you might need to dig deep ...

The 2 analogWrites() take now 23 usec , that could be much less. Left as an exercise :wink:
// see number of similar project here - FastDigitalWrite | Mecharobotics's Blog - )

Thanks for your help and guidance with this project Rob! Looking forward to trying it out on Monday..

Grumpy_Mike:
Sorry but it does not look like a Lorentz attractor to me, is is supposed to look like a butterfly. I have plotted 3D versions in the past.

I could be wrong but I think on the arduino it's just very slow to process and therefor the butterfly effect is lost.. I have a MEGA here too. Isn't the clock speed of the MEGA faster? I might try it on that and see if there is a difference.

. Isn't the clock speed of the MEGA faster?

No

minor improvement (1-2 usec), one addition precalculated

  if (x <= -1) h = (m1*x) - m3; 
  else if (x >= 1) h = (m1*x) + m3; 
  else h = (m1 + m3) * x;

==>

float m4 = m1 + m3;
....

    if (x <= -1) h = (m1*x) - m3; 
    else if (x >= 1) h = (m1*x) + m3; 
    else h = m4*x;

Next step is to see if the multiplication by 50 can be "assimilated"

*50 removed to the precalculated values. (it is a tricky one)

Loop clocks now at about 130-131 uSec, ==> 46% faster.
Time to stop optimizing and test if it still works (debug output seems reasonable similar to original sketch - note: h is factor 50 larger.).

In summary, the optimizations reduced the number of float multiplications from 10 to 4 and the number of additions from 15 to 10.

Note: these kind of optimizations should only be done if the additional speed adds value for an application.

// http://forum.arduino.cc/index.php?topic=168375.0
/* 
 Chaos Theory on your Scope! 
 Displays the lorenz attractor on your Oscilloscope.
 Author: Chris Gozzard (cgozzard@yahoo.com)
 Created: 25-5-13
 License: This code CC-BY-SA 3.0 and is unsupported.
 (see creativecommons.org/licenses for info)
 
 the following circuit is on both PWM ports (5 and 6)
 
 	       R
 PWM OUT ----/\/\/\-----+------------ OUTPUT
 		       |
 		      === C
 |
 GND
 
 R = 10k
 C = 0.1uF		
 
 Many thanks to John M. De Cristofaro for the original idea (based on his Oscilloscope Christmas Tree circuit)  
 */

int x_pos;         
int y_pos; 
int x_out = 5;
int y_out = 6;

float x = 0.70;
float y = 0;
float z = 0;

float dt = 0.01;
float alpha  = 15.60*dt;
float beta   = -28*dt; 
float m0     = -1.143;
float m1     = -0.714;

float m2 = 0.5*(m0-m1);
float m3 = 2 * m2 *50;
float m4 = (m1*50+m3);

float xdot;
float ydot;
float zdot;

float h;

void setup()  
{ 
  Serial.begin(115200);
  pinMode(x_out, OUTPUT);
  pinMode(y_out, OUTPUT);

  // this next section changes the PWM frequency - don't mess with it!
   TCCR0A = (	1<<COM0A1 | 0<<COM0A0 |		// clear OC0A on compare match (hi-lo PWM)
   1<<COM0B1 | 0<<COM0B0 |		// clear OC0B on compare match (hi-lo PWM)
   1<<WGM01  | 1<<WGM00);		// set PWM lines at 0xFF
   
   TCCR0B = (	0<<FOC0A | 0<<FOC0B |		// no force compare match
   0<<WGM02 |			// set PWM lines at 0xFF
   0<<CS02	 | 0<<CS01 |		// use system clock (no divider)
   1<<CS00 );
   
   TIMSK0 = (	0<<OCIE0B | 0<<TOIE0 |
   0<<OCIE0A ); 
} 

void loop()  
{
//  unsigned long start = micros();
//
//  for (int i=0; i < 10000; i++)
//  {
    float a = y-x;
    xdot = alpha * (a-h); 
    ydot = (z-a) * dt; 
    zdot = beta * y;

    if (x <= -50) h = (m1*x) - m3; 
    else if (x >= 50) h = (m1*x) + m3; 
    else h = m4*x; 
    //    Serial.print(h);
    //    Serial.print("\t");

    x = x + xdot;  
    y = y + ydot;
    z = z + zdot;

    x_pos = 128 + (int)(x); 
    y_pos = 128 + (int)(y);

    //    Serial.print(x_pos);
    //    Serial.print("\t");
    //    Serial.println(y_pos);
    analogWrite(x_out, x_pos);
    analogWrite(y_out, y_pos);
//  }
//  Serial.println((micros()-start)/10000.0);
}

Hi Rob. This is looking much faster. Not tried on the scope yet but can see from the serial data it seems much faster. my only concern is the x_out seems to go over 256 (the limit of the PWM output) is there some kind of mapping you could recommend? at the moment my solution is simply to apply x/3 to cut it down but I think you might have a more elegant solution... If I leave it out I think the the scope image will appear like its clipping..? anyway, we shall see tomorrow. If I can borrow my bosses iPad I'll take some video.. thanks for all your help and advice.. This is a great hobby. I'm really enjoying these little projects!

simply to apply x/3 to cut it down

shifting is a much more efficient way of getting a division however that restricts you to dividing by powers of two.

Hi Rob. This is looking much faster. Not tried on the scope yet but can see from the serial data it seems much faster. my only concern is the x_out seems to go over 256 (the limit of the PWM output) is there some kind of mapping you could recommend?

I could not test intermediate versions except for serial output as I don't have a scope available at the moment.
That is asking for failure I know :wink: However you get an idea how to optimize the code as the versions are not too far apart.

I guess we will find out tomorrow!

hmm.. somethings gone wrong. Let me go through the changes when I have time and see where we have slipped up.. :astonished: